diff --git a/eval/eval_dd_package.cpp b/eval/eval_dd_package.cpp
index 15e9d3a00..a756577e3 100644
--- a/eval/eval_dd_package.cpp
+++ b/eval/eval_dd_package.cpp
@@ -71,25 +71,16 @@ MatrixDD buildFunctionality(const qc::Grover* qc, Package<Config>& dd) {
 
   auto e = iteration;
   dd.incRef(e);
-
   for (std::size_t i = 0U; i < qc->iterations - 1U; ++i) {
-    auto f = dd.multiply(iteration, e);
-    dd.incRef(f);
-    dd.decRef(e);
-    e = f;
-    dd.garbageCollect();
+    e = dd.applyOperation(iteration, e);
   }
+  dd.decRef(iteration);
 
   QuantumComputation setup(qc->getNqubits());
   qc->setup(setup);
-  auto g = buildFunctionality(&setup, dd);
-  auto f = dd.multiply(e, g);
-  dd.incRef(f);
-  dd.decRef(e);
+  const auto g = buildFunctionality(&setup, dd);
+  e = dd.applyOperation(e, g);
   dd.decRef(g);
-  e = f;
-
-  dd.decRef(iteration);
   return e;
 }
 
@@ -103,15 +94,13 @@ MatrixDD buildFunctionalityRecursive(const qc::Grover* qc,
   auto iter = buildFunctionalityRecursive(&groverIteration, dd);
   auto e = iter;
   std::bitset<128U> iterBits(qc->iterations);
-  auto msb = static_cast<std::size_t>(std::floor(std::log2(qc->iterations)));
+  const auto msb =
+      static_cast<std::size_t>(std::floor(std::log2(qc->iterations)));
   auto f = iter;
   dd.incRef(f);
   bool zero = !iterBits[0U];
   for (std::size_t j = 1U; j <= msb; ++j) {
-    auto tmp = dd.multiply(f, f);
-    dd.incRef(tmp);
-    dd.decRef(f);
-    f = tmp;
+    f = dd.applyOperation(f, f);
     if (iterBits[j]) {
       if (zero) {
         dd.incRef(f);
@@ -119,11 +108,7 @@ MatrixDD buildFunctionalityRecursive(const qc::Grover* qc,
         e = f;
         zero = false;
       } else {
-        auto g = dd.multiply(e, f);
-        dd.incRef(g);
-        dd.decRef(e);
-        e = g;
-        dd.garbageCollect();
+        e = dd.applyOperation(e, f);
       }
     }
   }
@@ -132,14 +117,8 @@ MatrixDD buildFunctionalityRecursive(const qc::Grover* qc,
   // apply state preparation setup
   qc::QuantumComputation statePrep(qc->getNqubits());
   qc->setup(statePrep);
-  auto s = buildFunctionality(&statePrep, dd);
-  auto tmp = dd.multiply(e, s);
-  dd.incRef(tmp);
-  dd.decRef(s);
-  dd.decRef(e);
-  e = tmp;
-
-  return e;
+  const auto s = buildFunctionality(&statePrep, dd);
+  return dd.applyOperation(e, s);
 }
 
 std::unique_ptr<SimulationExperiment>
@@ -302,9 +281,7 @@ class BenchmarkDDPackage {
       qc::QuantumComputation statePrep(qc->getNqubits());
       qc->setup(statePrep);
       auto s = buildFunctionality(&statePrep, *dd);
-      auto e = dd->multiply(s, dd->makeZeroState(qc->getNqubits()));
-      dd->incRef(e);
-      dd->decRef(s);
+      auto e = dd->applyOperation(s, dd->makeZeroState(qc->getNqubits()));
 
       qc::QuantumComputation groverIteration(qc->getNqubits());
       qc->oracle(groverIteration);
@@ -318,17 +295,10 @@ class BenchmarkDDPackage {
       dd->incRef(f);
       for (std::size_t j = 0U; j <= msb; ++j) {
         if (iterBits[j]) {
-          auto g = dd->multiply(f, e);
-          dd->incRef(g);
-          dd->decRef(e);
-          e = g;
-          dd->garbageCollect();
+          e = dd->applyOperation(f, e);
         }
         if (j < msb) {
-          auto tmp = dd->multiply(f, f);
-          dd->incRef(tmp);
-          dd->decRef(f);
-          f = tmp;
+          f = dd->applyOperation(f, f);
         }
       }
       dd->decRef(f);
diff --git a/include/mqt-core/dd/Operations.hpp b/include/mqt-core/dd/Operations.hpp
index f2a6b3aa9..38f05afcb 100644
--- a/include/mqt-core/dd/Operations.hpp
+++ b/include/mqt-core/dd/Operations.hpp
@@ -248,11 +248,7 @@ Edge<Node> applyUnitaryOperation(const qc::Operation* op, const Edge<Node>& in,
                                  const qc::Permutation& permutation = {}) {
   static_assert(std::is_same_v<Node, dd::vNode> ||
                 std::is_same_v<Node, dd::mNode>);
-  auto tmp = dd.multiply(getDD(op, dd, permutation), in);
-  dd.incRef(tmp);
-  dd.decRef(in);
-  dd.garbageCollect();
-  return tmp;
+  return dd.applyOperation(getDD(op, dd, permutation), in);
 }
 
 template <class Config>
diff --git a/include/mqt-core/dd/Package.hpp b/include/mqt-core/dd/Package.hpp
index a951b97ea..edbc009ec 100644
--- a/include/mqt-core/dd/Package.hpp
+++ b/include/mqt-core/dd/Package.hpp
@@ -1476,6 +1476,30 @@ template <class Config> class Package {
     }
   }
 
+  /**
+   * @brief Applies a matrix operation to a matrix or vector.
+   *
+   * @details The reference count of the input matrix or vector is decreased,
+   * while the reference count of the result is increased. After the operation,
+   * garbage collection is triggered.
+   *
+   * @tparam Node Node type
+   * @param operation Matrix operation to apply
+   * @param e Matrix or vector to apply the operation to
+   * @return The appropriately reference-counted result.
+   */
+  template <class Node>
+  Edge<Node> applyOperation(const mEdge& operation, const Edge<Node>& e) {
+    static_assert(std::disjunction_v<std::is_same<Node, vNode>,
+                                     std::is_same<Node, mNode>>,
+                  "Node must be a vector or matrix node.");
+    const auto tmp = multiply(operation, e);
+    incRef(tmp);
+    decRef(e);
+    garbageCollect();
+    return tmp;
+  }
+
   dEdge applyOperationToDensity(dEdge& e, const mEdge& operation) {
     const auto tmp0 = conjugateTranspose(operation);
     const auto tmp1 = multiply(e, densityFromMatrixEdge(tmp0), false);
@@ -1827,7 +1851,7 @@ template <class Config> class Package {
           "expectation value.");
     }
 
-    auto yPrime = multiply(x, y);
+    const auto yPrime = multiply(x, y);
     const ComplexValue expValue = innerProduct(y, yPrime);
 
     assert(RealNumber::approximatelyZero(expValue.i));
diff --git a/src/dd/Simulation.cpp b/src/dd/Simulation.cpp
index a1b5ef815..613fcb9be 100644
--- a/src/dd/Simulation.cpp
+++ b/src/dd/Simulation.cpp
@@ -285,10 +285,7 @@ void extractProbabilityVectorRecursive(const QuantumComputation* qc,
       if (RealNumber::approximatelyEquals(pone, 1.)) {
         const qc::MatrixDD xGate =
             dd.makeGateDD(X_MAT, static_cast<Qubit>(permutation.at(target)));
-        const qc::VectorDD resetState = dd.multiply(xGate, state);
-        dd.incRef(resetState);
-        dd.decRef(state);
-        state = resetState;
+        state = dd.applyOperation(xGate, state);
         continue;
       }