[WIP] Partial equivalence checking

include/mqt-core/QuantumComputation.hpp

@@ -306,6 +306,14 @@ class QuantumComputation {
   [[nodiscard]] std::size_t getNqubitsWithoutAncillae() const {
     return nqubits;
   }
+  [[nodiscard]] std::size_t getNmeasuredQubits() const {
+    return static_cast<std::size_t>(
+        std::count(getGarbage().begin(), getGarbage().end(), false));
+  }
+  [[nodiscard]] std::size_t getNgarbageQubits() const {
+    return static_cast<std::size_t>(
+        std::count(getGarbage().begin(), getGarbage().end(), true));
+  }
   [[nodiscard]] std::size_t getNcbits() const { return nclassics; }
   [[nodiscard]] std::string getName() const { return name; }
   [[nodiscard]] const QuantumRegisterMap& getQregs() const { return qregs; }

include/mqt-core/dd/FunctionalityConstruction.hpp

@@ -8,11 +8,15 @@ namespace dd {
 using namespace qc;
 
 template <class Config>
-MatrixDD buildFunctionality(const QuantumComputation* qc, Package<Config>& dd);
+MatrixDD buildFunctionality(const QuantumComputation* qc, Package<Config>& dd,
+                            bool reduceAncillaryQubits = true,
+                            bool reduceGarbageQubits = true);
 
 template <class Config>
 MatrixDD buildFunctionalityRecursive(const QuantumComputation* qc,
-                                     Package<Config>& dd);
+                                     Package<Config>& dd,
+                                     bool reduceAncillaryQubits = true,
+                                     bool reduceGarbageQubits = true);
 
 template <class Config>
 bool buildFunctionalityRecursive(const QuantumComputation* qc,

include/mqt-core/dd/Package.hpp

@@ -2439,6 +2439,260 @@ template <class Config> class Package {
     nodes[index] = r.p;
     return r;
   }
+
+  // only keeps the first 2^d columns
+  mEdge setColumnsToZero(const mEdge& e, const Qubit d) {
+    if (e.isTerminal()) {
+      return e;
+    }
+    // the matrix of the current DD has dimensions 2^n x 2^n
+    const auto n = e.p->v + 1;
+
+    if (d >= n) {
+      return e;
+    }
+
+    std::array<mEdge, NEDGE> edges{};
+    edges[0] = setColumnsToZero(e.p->e[0], d);
+    edges[1] = mEdge::zero();
+    edges[2] = setColumnsToZero(e.p->e[2], d);
+    edges[3] = mEdge::zero();
+
+    auto f = makeDDNode(e.p->v, edges);
+    f.w = cn.lookup(e.w * f.w);
+    return f;
+  }
+
+  mEdge keepOnlyIthRow(const mEdge& e, const std::int64_t i) {
+    if (e.isZeroTerminal()) {
+      return e;
+    }
+    if (e.isTerminal()) {
+      if (i == 0) {
+        return e;
+      }
+      return mEdge::zero();
+    }
+    // the matrix of the current DD has dimensions 2^n x 2^n
+    const auto n = e.p->v + 1;
+    const auto twoPowNMinusOne = 1 << (n - 1);
+    std::array<mEdge, NEDGE> edges{};
+    if (i < twoPowNMinusOne) {
+      edges[0] = keepOnlyIthRow(e.p->e[0], i);
+      edges[1] = keepOnlyIthRow(e.p->e[1], i);
+      edges[2] = mEdge::zero();
+      edges[3] = mEdge::zero();
+    } else {
+      edges[0] = mEdge::zero();
+      edges[1] = mEdge::zero();
+      edges[2] = keepOnlyIthRow(e.p->e[2], i - twoPowNMinusOne);
+      edges[3] = keepOnlyIthRow(e.p->e[3], i - twoPowNMinusOne);
+    }
+
+    auto f = makeDDNode(e.p->v, edges);
+    f.w = cn.lookup(e.w * f.w);
+    return f;
+  }
+
+  UnaryComputeTable<mEdge, mEdge> shiftAllMatrixRows{};
+  Qubit shiftAllMatrixRowsM{0};
+  Qubit shiftAllMatrixRowsD{0};
+
+  /**
+    Equally divides the rows of the matrix represented by e (of size 2^n x 2^n)
+    into 2^g parts of size 2^m (where g = n - m).
+    For each part the (upperOffset + i)-th row is shifted by i*2^d columns.
+    **/
+  mEdge shiftAllRowsRecursive(const mEdge& e, const Qubit m, const Qubit d,
+                              std::int64_t upperOffset) {
+    if (e.isTerminal() && upperOffset == 0) {
+      return e;
+    }
+    if (e.isTerminal()) {
+      return mEdge::zero();
+    }
+    if (upperOffset == 0 && m == 0) {
+      return e;
+    }
+
+    // the matrix of the current DD has dimensions 2^n x 2^n
+    const auto n = e.p->v + 1;
+    if (upperOffset >= (1 << n) || upperOffset < 0) {
+      return mEdge::zero();
+    }
+
+    const mEdge eCopy{e.p, Complex::one()};
+    // check if it's in the compute table with an edge weight of one
+    // only store elements in the compute table if the offset is 0
+    if (upperOffset == 0) {
+      if (const auto* r = shiftAllMatrixRows.lookup(eCopy); r != nullptr) {
+        auto f = *r;
+        f.w = cn.lookup(e.w * f.w);
+        return f;
+      }
+    }
+
+    if (d >= n && m >= n) {
+      return keepOnlyIthRow(e, upperOffset);
+    }
+
+    std::array<mEdge, NEDGE> edges{};
+    const auto originalUpperOffset = upperOffset;
+
+    // if the current submatrix size is less than 2^m,
+    // then the offset of the lower half needs to be adapted
+    const bool adaptOffsetOfLowerMatrixHalf = n <= m;
+    edges[0] = shiftAllRowsRecursive(e.p->e[0], m, d, upperOffset);
+    if (n <= d) {
+      // CASE 1: this (sub)matrix has size < 2^d
+      // -> we need to keep all the columns of this submatrix
+      edges[1] = shiftAllRowsRecursive(e.p->e[1], m, d, upperOffset);
+      if (adaptOffsetOfLowerMatrixHalf) {
+        upperOffset -= (1 << (n - 1));
+      }
+      edges[2] = shiftAllRowsRecursive(e.p->e[2], m, d, upperOffset);
+      edges[3] = shiftAllRowsRecursive(e.p->e[3], m, d, upperOffset);
+    } else {
+      const std::int64_t addedOffset{1 << (n - 1 - d)};
+      if (upperOffset + addedOffset < (1 << m)) {
+        // CASE 2: this submatrix doesn't contain ancilla qubits
+        // -> we don't need to set any columns to zero
+        edges[1] =
+            shiftAllRowsRecursive(e.p->e[0], m, d, upperOffset + addedOffset);
+
+        if (adaptOffsetOfLowerMatrixHalf) {
+          upperOffset -= (1 << (n - 1));
+        }
+        edges[2] = shiftAllRowsRecursive(e.p->e[2], m, d, upperOffset);
+        edges[3] =
+            shiftAllRowsRecursive(e.p->e[2], m, d, upperOffset + addedOffset);
+      } else {
+        // CASE 3: the right half of the matrix represents the ancilla qubits,
+        // therefore it is set to zero
+        edges[1] = mEdge::zero();
+        if (adaptOffsetOfLowerMatrixHalf) {
+          upperOffset -= (1 << (n - 1));
+        }
+        edges[2] = shiftAllRowsRecursive(e.p->e[2], m, d, upperOffset);
+        edges[3] = mEdge::zero();
+      }
+    }
+
+    auto f = makeDDNode(e.p->v, edges);
+    // add to the compute table with a weight of 1
+    if (originalUpperOffset == 0) {
+      shiftAllMatrixRows.insert(eCopy, f);
+    }
+    f.w = cn.lookup(e.w * f.w);
+    return f;
+  }
+
+public:
+  /**
+      Equally divides the rows of the matrix represented by e into parts
+      of size 2^m.
+      For each part we keep only the first 2^d columns
+      and the i-th row of each part is shifted by i*2^d columns.
+  **/
+  mEdge shiftAllRows(const mEdge& e, const Qubit m, const Qubit d) {
+    if (shiftAllMatrixRowsM != m || shiftAllMatrixRowsD != d) {
+      shiftAllMatrixRows.clear();
+      shiftAllMatrixRowsM = m;
+      shiftAllMatrixRowsD = d;
+    }
+    return shiftAllRowsRecursive(e, m, d, 0);
+  }
+
+private:
+  mEdge partialEquivalenceCheckSubroutine(mEdge u, const Qubit m, const Qubit k,
+                                          const Qubit extra) {
+    // add extra ancillary qubits
+    if (extra > 0) {
+      if (u.p->v + 1U + extra > nqubits) {
+        resize(u.p->v + 1U + extra);
+      }
+      u = kronecker(makeIdent(extra), u);
+    }
+    if (u.isTerminal()) {
+      return u;
+    }
+    const auto n = static_cast<Qubit>(u.p->v + 1);
+    const Qubit d = n - k;
+    u = setColumnsToZero(u, d);
+    const auto u2 = shiftAllRows(u, m, d);
+    return multiply(conjugateTranspose(u), u2);
+  }
+
+public:
+  /**
+   Checks for partial equivalence between the two circuits u1 and u2,
+    where the first d qubits of the circuits are the data qubits and
+    the first m qubits are the measured qubits.
+   @param u1 DD representation of first circuit
+   @param u2 DD representation of second circuit
+   @param d Number of data qubits
+   @param m Number of measured qubits
+   @return true if the two circuits u1 and u2 are partially equivalent.
+   **/
+  bool partialEquivalenceCheck(mEdge u1, mEdge u2, const Qubit d,
+                               const Qubit m) {
+    if (m == 0) {
+      return true;
+    }
+    if (u1.isTerminal() && u2.isTerminal()) {
+      return u1 == u2;
+    }
+    if (u1.isZeroTerminal() || u2.isZeroTerminal()) {
+      return false;
+    }
+    // add qubits such that u1 and u2 have the same dimension
+    if (u1.isTerminal()) {
+      auto w = u1.w;
+      u1 = makeIdent(u2.p->v + 1);
+      u1.w = w;
+    } else if (u2.isTerminal()) {
+      auto w = u2.w;
+      u2 = makeIdent(u1.p->v + 1);
+      u2.w = w;
+    } else if (u1.p->v < u2.p->v) {
+      u1 = kronecker(makeIdent(u2.p->v - u1.p->v), u1);
+    } else if (u1.p->v > u2.p->v) {
+      u2 = kronecker(makeIdent(u1.p->v - u2.p->v), u2);
+    }
+
+    const Qubit n = u1.p->v + 1;
+    Qubit k = n - d;
+    Qubit extra{0};
+    if (m > k) {
+      extra = m - k;
+    }
+    k = k + extra;
+
+    const auto u1Prime = partialEquivalenceCheckSubroutine(u1, m, k, extra);
+    const auto u2Prime = partialEquivalenceCheckSubroutine(u2, m, k, extra);
+
+    return u1Prime == u2Prime;
+  }
+
+  /**
+      Checks for partial equivalence between the two circuits u1 and u2,
+       where all qubits of the circuits are the data qubits and
+       the first m qubits are the measured qubits.
+      @param u1 DD representation of first circuit
+      @param u2 DD representation of second circuit
+      @param m Number of measured qubits
+      @return true if the two circuits u1 and u2 are partially equivalent.
+      **/
+  bool zeroAncillaePartialEquivalenceCheck(const mEdge& u1, const mEdge& u2,
+                                           const Qubit m) {
+    auto u1u2 = multiply(u1, conjugateTranspose(u2));
+    const Qubit n = u1.p->v + 1;
+    std::vector<bool> garbage(n, false);
+    for (size_t i = m; i < n; i++) {
+      garbage[i] = true;
+    }
+    return isCloseToIdentity(u1u2, 1.0E-10, garbage, false);
+  }
 };
 
 } // namespace dd

include/mqt-core/dd/Verification.hpp

@@ -0,0 +1,94 @@
+#include "dd/FunctionalityConstruction.hpp"
+#include "dd/Package.hpp"
+
+namespace dd {
+/**
+    Checks for partial equivalence between the two circuits c1 and c2
+    that have no ancilla qubits.
+    Assumption: the input and output permutations are the same.
+
+    @param circuit1 First circuit
+    @param circuit2 Second circuit
+    @return true if the two circuits c1 and c2 are partially equivalent.
+    **/
+template <class Config>
+bool zeroAncillaePartialEquivalenceCheck(
+    qc::QuantumComputation c1, qc::QuantumComputation c2,
+    const std::unique_ptr<dd::Package<Config>>& dd) {
+  if (c1.getNqubits() != c2.getNqubits() ||
+      c1.getGarbage() != c2.getGarbage()) {
+    throw std::invalid_argument(
+        "The circuits need to have the same number of qubits and the same "
+        "permutation of input and output qubits.");
+  }
+  c2.invert();
+  for (auto& gate : c1) {
+    c2.emplace_back(gate);
+  }
+
+  const auto u = buildFunctionality(&c2, *dd, false, false);
+
+  return dd->isCloseToIdentity(u, 1.0E-10, c1.getGarbage(), false);
+}
+// get next garbage qubit after n
+inline Qubit getNextGarbage(Qubit n, const std::vector<bool>& garbage) {
+  while (n < static_cast<Qubit>(garbage.size()) && !garbage.at(n)) {
+    n++;
+  }
+  return n;
+}
+/**
+    Checks for partial equivalence between the two circuits c1 and c2.
+    Assumption: the data qubits are all at the beginning of the input qubits and
+    the input and output permutations are the same.
+
+    @param circuit1 First circuit
+    @param circuit2 Second circuit
+    @return true if the two circuits c1 and c2 are partially equivalent.
+    **/
+template <class Config>
+bool partialEquivalenceCheck(qc::QuantumComputation c1,
+                             qc::QuantumComputation c2,
+                             const std::unique_ptr<dd::Package<Config>>& dd) {
+
+  const auto d1 = c1.getNqubitsWithoutAncillae();
+  const auto d2 = c2.getNqubitsWithoutAncillae();
+  const auto m1 = c1.getNmeasuredQubits();
+  const auto m2 = c2.getNmeasuredQubits();
+  if (m1 != m2 || d1 != d2) {
+    return false;
+  }
+  const auto n1 = static_cast<Qubit>(c1.getNqubits());
+  const auto n2 = static_cast<Qubit>(c2.getNqubits());
+  if (d1 == n1 && d2 == n2) {
+    // no ancilla qubits
+    return zeroAncillaePartialEquivalenceCheck(c1, c2, dd);
+  }
+  // add swaps in order to put the measured (= not garbage) qubits in the end
+  const auto garbage1 = c1.getGarbage();
+
+  auto nextGarbage = getNextGarbage(0, garbage1);
+  // find the first garbage qubit at the end
+  for (std::int64_t i = std::min(n1, n2) - 1;
+       i >= static_cast<std::int64_t>(m1); i--) {
+    if (!garbage1.at(static_cast<Qubit>(i))) {
+      // swap it to the beginning
+      c1.swap(static_cast<Qubit>(i), nextGarbage);
+      c2.swap(static_cast<Qubit>(i), nextGarbage);
+      ++nextGarbage;
+      nextGarbage = getNextGarbage(nextGarbage, garbage1);
+    }
+  }
+
+  // partialEquivalenceCheck with dd
+
+  const auto u1 = buildFunctionality(&c1, *dd, false, false);
+  const auto u2 = buildFunctionality(&c2, *dd, false, false);
+
+  return dd->partialEquivalenceCheck(u1, u2, static_cast<Qubit>(d1),
+                                     static_cast<Qubit>(m1));
+}
+
+std::pair<qc::QuantumComputation, qc::QuantumComputation>
+generateRandomBenchmark(size_t n, Qubit d, Qubit m);
+} // namespace dd

src/dd/CMakeLists.txt

@@ -19,6 +19,7 @@ if(NOT TARGET ${MQT_CORE_TARGET_NAME}-dd)
     RealNumber.cpp
     RealNumberUniqueTable.cpp
     Simulation.cpp
+    Verification.cpp
     statistics/MemoryManagerStatistics.cpp
     statistics/Statistics.cpp
     statistics/TableStatistics.cpp