Gate Reordering and Approximate Comparison (#13)

* ✨ add reorder operations pass to pre-processing
* ✨ add functionality to compare the resulting DD to the identity with a given threshold
* ✨ introduce new optimization that tries to cancel identical operations in the alternating strategy whenever possible
* ⬆️ update qfr

Author: burgholzer

Diff for: README.md

@@ -70,11 +70,12 @@ The verification procedure can be configured with the following settings and opt
         - simulation
     - `tolerance`: Numerical tolerance used during computation (`1e-13` per default)
 - Settinggs for the ![G \rightarrow \mathbb{I} \leftarrow G'](https://render.githubusercontent.com/render/math?math=G%20%5Crightarrow%20%5Cmathbb%7BI%7D%20%5Cleftarrow%20G') method:
-    - `strategy`: strategy to use for the scheme
+    - `strategy`: Strategy to use for the scheme
         - naive
         - proportional (*default*)
         - lookahead
         - compilationflow
+    - `identity_threshold`: Numerical tolerance used for checking the similarity to the identity matrix (`1e-10` per default)
 - Settings for the simulation-based method:
     - `fidelity`: Fidelity limit for comparison (`0.999` per default)
     - `max_sims`: Maximum number of simulations to conduct (`16` per default)
@@ -89,6 +90,7 @@ The verification procedure can be configured with the following settings and opt
     - `fuse_single_qubit_gates`: Fuse consecutive single qubit gates (*on* per default)
     - `remove_diagonal_gates_before_measure`: Remove diagonal gates before measurements (*off* by default)
     - `tranform_dynamic_circuit`: Transform dynamic to static circuit (*off* by default)
+    - `reorder_operations`: Reorder operations in order to eliminate unnecessary dependencies (*on* by default)
 
 The `qcec.Results` class that is returned by the `verify` function provides `json()` and `csv()` methods to produce JSON or CSV formatted output.
 

Diff for: include/EquivalenceChecker.hpp

@@ -19,15 +19,17 @@ namespace ec {
                             RIGHT = false };
 
     struct Configuration {
-        ec::Method   method    = ec::Method::G_I_Gp;
-        ec::Strategy strategy  = ec::Strategy::Proportional;
-        dd::fp       tolerance = dd::ComplexTable<>::tolerance();
+        ec::Method   method            = ec::Method::G_I_Gp;
+        ec::Strategy strategy          = ec::Strategy::Proportional;
+        dd::fp       tolerance         = dd::ComplexTable<>::tolerance();
+        dd::fp       identityThreshold = 1e-10;
 
         // configuration options for optimizations
         bool fuseSingleQubitGates             = true;
         bool reconstructSWAPs                 = true;
         bool removeDiagonalGatesBeforeMeasure = false;
         bool transformDynamicCircuit          = false;
+        bool reorderOperations                = true;
 
         // configuration options for PowerOfSimulation equivalence checker
         double      fidelity_limit = 0.999;
@@ -40,7 +42,8 @@ namespace ec {
             nlohmann::json config{};
             config["method"] = ec::toString(method);
             if (method == ec::Method::G_I_Gp) {
-                config["strategy"] = ec::toString(strategy);
+                config["strategy"]           = ec::toString(strategy);
+                config["identity threshold"] = identityThreshold;
             }
             config["tolerance"]                                   = tolerance;
             config["optimizations"]                               = {};
@@ -49,6 +52,7 @@ namespace ec {
             optimizations["reconstruct swaps"]                    = reconstructSWAPs;
             optimizations["remove diagonal gates before measure"] = removeDiagonalGatesBeforeMeasure;
             optimizations["transform dynamic circuit"]            = transformDynamicCircuit;
+            optimizations["reorder operations"]                   = reorderOperations;
             if (method == ec::Method::Simulation) {
                 config["simulation config"]               = {};
                 auto& simulation                          = config["simulation config"];

Diff for: include/ImprovedDDEquivalenceChecker.hpp

@@ -33,14 +33,72 @@ namespace ec {
         /// \return initial matrix
         qc::MatrixDD createInitialMatrix();
 
-        /// Create the goal matrix used for the G->I<-G' scheme.
-        /// [1 0] if the qubit is no ancillary
-        /// [0 1]
-        ///
-        /// [1 0] for an ancillary that is present in either circuit
-        /// [0 0]
-        /// \return goal matrix
-        qc::MatrixDD createGoalMatrix();
+        bool gatesAreIdentical(const qc::Permutation& perm1, const qc::Permutation& perm2) {
+            // exit prematurely whenever one of the circuits is already finished
+            if (it1 == end1 || it2 == end2)
+                return false;
+
+            const auto& op1 = **it1;
+            const auto& op2 = **it2;
+
+            // check type
+            if (op1.getType() != op2.getType()) {
+                return false;
+            }
+
+            // check number of controls
+            const auto nc1 = op1.getNcontrols();
+            const auto nc2 = op2.getNcontrols();
+            if (nc1 != nc2) {
+                return false;
+            }
+
+            // SWAPs are not handled by this routine
+            if (op1.getType() == qc::SWAP && nc1 == 0) {
+                return false;
+            }
+
+            // check parameters
+            const auto param1 = op1.getParameter();
+            const auto param2 = op2.getParameter();
+            for (std::size_t p = 0; p < qc::MAX_PARAMETERS; ++p) {
+                // it might make sense to use fuzzy comparison here
+                if (param1[p] != param2[p]) {
+                    return false;
+                }
+            }
+
+            // check controls
+            if (nc1 != 0) {
+                dd::Controls controls1{};
+                for (const auto& control: op1.getControls()) {
+                    controls1.emplace(dd::Control{perm1.at(control.qubit), control.type});
+                }
+                dd::Controls controls2{};
+                for (const auto& control: op2.getControls()) {
+                    controls2.emplace(dd::Control{perm2.at(control.qubit), control.type});
+                }
+                if (controls1 != controls2) {
+                    return false;
+                }
+            }
+
+            // check targets
+            std::set<dd::Qubit> targets1{};
+            for (const auto& target: op1.getTargets()) {
+                targets1.emplace(perm1.at(target));
+            }
+            std::set<dd::Qubit> targets2{};
+            for (const auto& target: op2.getTargets()) {
+                targets2.emplace(perm2.at(target));
+            }
+            if (targets1 != targets2) {
+                return false;
+            }
+
+            // operations are identical
+            return true;
+        }
 
     public:
         ImprovedDDEquivalenceChecker(qc::QuantumComputation& qc1, qc::QuantumComputation& qc2):

Diff for: jkq/qcec/bindings.cpp

@@ -126,6 +126,10 @@ PYBIND11_MODULE(pyqcec, m) {
                            R"pbdoc(
 					Numerical tolerance used during computation
 				)pbdoc")
+            .def_readwrite("identity_threshold", &ec::Configuration::identityThreshold,
+                           R"pbdoc(
+					Numerical threshold used for checking the similarity to the identity matrix (default: 1e-10)
+				)pbdoc")
             .def_readwrite("reconstruct_swaps", &ec::Configuration::reconstructSWAPs,
                            R"pbdoc(
 					Optimization pass reconstructing SWAP operations
@@ -138,6 +142,10 @@ PYBIND11_MODULE(pyqcec, m) {
                            R"pbdoc(
 					Optimization pass that transforms dynamic circuits to static circuits for equivalence checking
 				)pbdoc")
+            .def_readwrite("reorder_operations", &ec::Configuration::reorderOperations,
+                           R"pbdoc(
+					Optimization pass that reorders operations to eliminate unnecessary dependencies
+				)pbdoc")
             .def_readwrite("remove_diagonal_gates_before_measure", &ec::Configuration::removeDiagonalGatesBeforeMeasure,
                            R"pbdoc(
 					Optimization pass removing diagonal gates before measurements

Diff for: src/CompilationFlowEquivalenceChecker.cpp

@@ -33,6 +33,12 @@ namespace ec {
             }
 
             if (it1 != end1 && it2 != end2) {
+                if (results.result.p->ident && gatesAreIdentical(perm1, perm2)) {
+                    ++it1;
+                    ++it2;
+                    continue;
+                }
+
                 auto cost1 = costFunction((*it1)->getType(), (*it1)->getControls().size());
                 auto cost2 = costFunction((*it2)->getType(), (*it2)->getControls().size());
 
@@ -78,8 +84,15 @@ namespace ec {
         results.result = dd->reduceAncillae(results.result, ancillary1, LEFT);
         results.result = dd->reduceAncillae(results.result, ancillary2, RIGHT);
 
-        auto goal_matrix    = createGoalMatrix();
-        results.equivalence = equals(results.result, goal_matrix);
+        if (dd->isCloseToIdentity(results.result, config.identityThreshold)) {
+            if (!results.result.w.approximatelyOne()) {
+                results.equivalence = Equivalence::EquivalentUpToGlobalPhase;
+            } else {
+                results.equivalence = Equivalence::Equivalent;
+            }
+        } else {
+            results.equivalence = Equivalence::NotEquivalent;
+        }
 
         auto                          endVerification   = std::chrono::steady_clock::now();
         std::chrono::duration<double> preprocessingTime = endPreprocessing - start;

Diff for: src/EquivalenceChecker.cpp

@@ -10,7 +10,7 @@ namespace ec {
 
     EquivalenceChecker::EquivalenceChecker(qc::QuantumComputation& qc1, qc::QuantumComputation& qc2):
         qc1(qc1), qc2(qc2) {
-        // currently this modifies the underlying quantum circuits
+        // currently, this modifies the underlying quantum circuits
         // in the future this might want to be avoided
         qc1.stripIdleQubits();
         qc2.stripIdleQubits();
@@ -266,6 +266,11 @@ namespace ec {
             qc::CircuitOptimizer::singleQubitGateFusion(qc2);
         }
 
+        if (config.reorderOperations) {
+            qc::CircuitOptimizer::reorderOperations(qc1);
+            qc::CircuitOptimizer::reorderOperations(qc2);
+        }
+
         it1  = qc1.begin();
         it2  = qc2.begin();
         end1 = qc1.cend();

Diff for: src/ImprovedDDEquivalenceChecker.cpp

@@ -42,16 +42,6 @@ namespace ec {
         return e;
     }
 
-    qc::MatrixDD ImprovedDDEquivalenceChecker::createGoalMatrix() {
-        auto goalMatrix = dd->makeIdent(nqubits);
-        dd->incRef(goalMatrix);
-        goalMatrix = dd->reduceAncillae(goalMatrix, ancillary2, RIGHT);
-        goalMatrix = dd->reduceGarbage(goalMatrix, garbage2, RIGHT);
-        goalMatrix = dd->reduceAncillae(goalMatrix, ancillary1, LEFT);
-        goalMatrix = dd->reduceGarbage(goalMatrix, garbage1, LEFT);
-        return goalMatrix;
-    }
-
     /// Use dedicated method to check the equivalence of both provided circuits
     EquivalenceCheckingResults ImprovedDDEquivalenceChecker::check(const Configuration& config) {
         EquivalenceCheckingResults results{};
@@ -99,8 +89,16 @@ namespace ec {
         results.result = dd->reduceAncillae(results.result, ancillary1, LEFT);
         results.result = dd->reduceAncillae(results.result, ancillary2, RIGHT);
 
-        results.equivalence = equals(results.result, createGoalMatrix());
-        results.maxActive   = std::max(results.maxActive, dd->mUniqueTable.getMaxActiveNodes());
+        if (dd->isCloseToIdentity(results.result, config.identityThreshold)) {
+            if (!results.result.w.approximatelyOne()) {
+                results.equivalence = Equivalence::EquivalentUpToGlobalPhase;
+            } else {
+                results.equivalence = Equivalence::Equivalent;
+            }
+        } else {
+            results.equivalence = Equivalence::NotEquivalent;
+        }
+        results.maxActive = std::max(results.maxActive, dd->mUniqueTable.getMaxActiveNodes());
 
         auto                          endVerification   = std::chrono::steady_clock::now();
         std::chrono::duration<double> preprocessingTime = endPreprocessing - start;
@@ -114,6 +112,11 @@ namespace ec {
     /// Alternate between LEFT and RIGHT applications
     void ImprovedDDEquivalenceChecker::checkNaive(qc::MatrixDD& result, qc::Permutation& perm1, qc::Permutation& perm2) {
         while (it1 != end1 && it2 != end2) {
+            if (result.p->ident && gatesAreIdentical(perm1, perm2)) {
+                ++it1;
+                ++it2;
+                continue;
+            }
             applyGate(qc1, it1, result, perm1, LEFT);
             ++it1;
             applyGate(qc2, it2, result, perm2, RIGHT);
@@ -123,13 +126,17 @@ namespace ec {
 
     /// Alternate according to the gate count ratio between LEFT and RIGHT applications
     void ImprovedDDEquivalenceChecker::checkProportional(qc::MatrixDD& result, qc::Permutation& perm1, qc::Permutation& perm2) {
-        auto ratio  = static_cast<unsigned int>(std::round(
-                static_cast<double>(std::max(qc1.getNops(), qc2.getNops())) /
-                static_cast<double>(std::min(qc1.getNops(), qc2.getNops()))));
+        auto ratio  = static_cast<unsigned int>(std::round(static_cast<double>(std::max(qc1.getNops(), qc2.getNops())) / static_cast<double>(std::min(qc1.getNops(), qc2.getNops()))));
         auto ratio1 = (qc1.getNops() > qc2.getNops()) ? ratio : 1;
         auto ratio2 = (qc1.getNops() > qc2.getNops()) ? 1 : ratio;
 
         while (it1 != end1 && it2 != end2) {
+            if (result.p->ident && gatesAreIdentical(perm1, perm2)) {
+                ++it1;
+                ++it2;
+                continue;
+            }
+
             for (unsigned int i = 0; i < ratio1 && it1 != end1; ++i) {
                 applyGate(qc1, it1, result, perm1, LEFT);
                 ++it1;