Circuit::optimize

scaluq/circuit/circuit.cpp

@@ -2,6 +2,9 @@
 
 #include <ranges>
 
+#include "../gate/gate_factory.hpp"
+#include "../gate/merge_gate.hpp"
+
 namespace scaluq {
 UINT Circuit::calculate_depth() const {
     std::vector<UINT> filled_step(_n_qubits, 0ULL);
@@ -143,4 +146,113 @@ void Circuit::check_gate_is_valid(const ParamGate& param_gate) const {
         throw std::runtime_error("Gate to be added to the circuit has invalid qubit range");
     }
 }
+
+void Circuit::optimize(UINT block_size) {
+    if (block_size >= 3) {
+        throw std::runtime_error(
+            "Currently block_size >= 3 is not supported because general matrix gate with qubits >= "
+            "3 is not implemented.");
+    }
+    std::vector<GateWithKey> new_gate_list;
+    double global_phase = 0.;
+    std::vector<std::pair<Gate, std::vector<UINT>>> gate_pool;
+    constexpr UINT NO_GATES = std::numeric_limits<UINT>::max();
+    std::vector<UINT> latest_gate_idx(_n_qubits, NO_GATES);
+    for (const GateWithKey& gate_with_key : _gate_list) {
+        if (gate_with_key.index() == 1) {
+            const auto& pgate = std::get<1>(gate_with_key).first;
+            for (UINT target : pgate->get_target_qubit_list()) {
+                latest_gate_idx[target] = NO_GATES;
+            }
+            for (UINT control : pgate->get_control_qubit_list()) {
+                latest_gate_idx[control] = NO_GATES;
+            }
+            new_gate_list.emplace_back(std::move(gate_with_key));
+            continue;
+        }
+        const auto& gate = std::get<0>(gate_with_key);
+        if (gate.gate_type() == GateType::I) {
+            continue;
+        }
+        if (gate.gate_type() == GateType::GlobalPhase) {
+            global_phase += GlobalPhaseGate(gate)->phase();
+            continue;
+        }
+        auto target_list = gate->get_target_qubit_list();
+        auto control_list = gate->get_control_qubit_list();
+        std::vector<UINT> targets;
+        targets.reserve(target_list.size() + control_list.size());
+        std::ranges::copy(target_list, std::back_inserter(targets));
+        std::ranges::copy(control_list, std::back_inserter(targets));
+        std::vector<UINT> previous_gate_indices;
+        std::vector<UINT> newly_applied_qubits;
+        for (UINT target : targets) {
+            if (latest_gate_idx[target] == NO_GATES) {
+                newly_applied_qubits.push_back(target);
+            } else {
+                previous_gate_indices.push_back(latest_gate_idx[target]);
+            }
+        }
+        previous_gate_indices.erase(std::ranges::unique(previous_gate_indices).begin(),
+                                    previous_gate_indices.end());
+        UINT merged_gate_size =
+            std::accumulate(previous_gate_indices.begin(),
+                            previous_gate_indices.end(),
+                            newly_applied_qubits.size(),
+                            [&](UINT sz, UINT idx) { return sz + gate_pool[idx].second.size(); });
+        auto is_pauli = [](const Gate& gate) {
+            GateType type = gate.gate_type();
+            return type == GateType::I || type == GateType::X || type == GateType::Y ||
+                   type == GateType::Z || type == GateType::Pauli;
+        };
+        bool all_pauli =
+            is_pauli(gate) && std::ranges::all_of(previous_gate_indices, [&](UINT idx) {
+                return is_pauli(gate_pool[idx].first);
+            });
+        if (!all_pauli && merged_gate_size > block_size) {
+            for (UINT idx : previous_gate_indices) {
+                for (UINT qubit : gate_pool[idx].second) {
+                    latest_gate_idx[qubit] = NO_GATES;
+                }
+                new_gate_list.emplace_back(std::move(gate_pool[idx].first));
+            }
+            UINT new_idx = gate_pool.size();
+            for (UINT qubit : targets) {
+                latest_gate_idx[qubit] = new_idx;
+            }
+            gate_pool.emplace_back(std::move(gate), std::move(targets));
+            continue;
+        }
+        Gate merged_gate = gate::I();
+        UINT new_idx = gate_pool.size();
+        std::vector<UINT> new_targets;
+        for (UINT idx : previous_gate_indices) {
+            double phase;
+            std::tie(merged_gate, phase) = merge_gate(merged_gate, gate_pool[idx].first);
+            global_phase += phase;
+            for (UINT qubit : gate_pool[idx].second) {
+                new_targets.push_back(qubit);
+                latest_gate_idx[qubit] = new_idx;
+            }
+        }
+        {
+            double phase;
+            std::tie(merged_gate, phase) = merge_gate(merged_gate, gate);
+            global_phase += phase;
+            for (UINT qubit : newly_applied_qubits) {
+                new_targets.push_back(qubit);
+                latest_gate_idx[qubit] = new_idx;
+            }
+        }
+        gate_pool.emplace_back(std::move(merged_gate), std::move(new_targets));
+    }
+    std::ranges::sort(latest_gate_idx);
+    latest_gate_idx.erase(std::ranges::unique(latest_gate_idx).begin(), latest_gate_idx.end());
+    for (UINT idx : latest_gate_idx) {
+        if (idx == NO_GATES) continue;
+        new_gate_list.emplace_back(std::move(gate_pool[idx].first));
+    }
+    if (std::abs(global_phase) < 1e-12) new_gate_list.push_back(gate::GlobalPhase(global_phase));
+    _gate_list.swap(new_gate_list);
+}
 }  // namespace scaluq

scaluq/circuit/circuit.hpp

@@ -53,6 +53,8 @@ class Circuit {
     Circuit copy() const;
     Circuit get_inverse() const;
 
+    void optimize(UINT block_size = 2);
+
 private:
     UINT _n_qubits;
 

scaluq/operator/operator.cpp

@@ -21,7 +21,7 @@ std::string Operator::to_string() const {
 
 void Operator::add_operator(const PauliOperator& mpt) { add_operator(PauliOperator{mpt}); }
 void Operator::add_operator(PauliOperator&& mpt) {
-    _is_hermitian &= mpt.get_coef().imag() == 0.;
+    _is_hermitian &= (mpt.get_coef().imag() == 0.);
     if (![&] {
             const auto& target_list = mpt.get_target_qubit_list();
             if (target_list.empty()) return true;

scaluq/util/utility.hpp

@@ -48,6 +48,19 @@ KOKKOS_INLINE_FUNCTION matrix_2_2 matrix_multiply(const matrix_2_2& matrix1,
             matrix1.val[1][0] * matrix2.val[0][1] + matrix1.val[1][1] * matrix2.val[1][1]};
 }
 
+template <class T, class... Args>
+inline std::vector<T> concat_vector(Args... args) {
+    std::initializer_list<const std::vector<T>&> vecs = {args...};
+    std::vector<T> res;
+    res.reserve(std::accumulate(vecs.begin(), vecs.end(), 0ULL, [](UINT acc, const auto& vec) {
+        return acc + vec.size();
+    }));
+    for (const auto& vec : vecs) {
+        std::ranges::copy(vecs.begin(), vecs.end(), std::back_inserter(res));
+    }
+    return res;
+}
+
 inline ComplexMatrix kronecker_product(const ComplexMatrix& lhs, const ComplexMatrix& rhs) {
     ComplexMatrix result(lhs.rows() * rhs.rows(), lhs.cols() * rhs.cols());
     for (int i = 0; i < lhs.rows(); i++) {

tests/CMakeLists.txt

@@ -3,6 +3,7 @@ cmake_minimum_required(VERSION 3.21)
 enable_testing()
 
 add_executable(scaluq_test EXCLUDE_FROM_ALL
+    circuit/circuit_optimize_test.cpp
     circuit/circuit_test.cpp
     circuit/param_circuit_test.cpp
     gate/gate_test.cpp

tests/circuit/circuit_optimize_test.cpp

@@ -0,0 +1,187 @@
+#include <gtest/gtest.h>
+
+#include <circuit/circuit.hpp>
+#include <gate/gate_factory.hpp>
+#include <gate/param_gate_factory.hpp>
+#include <state/state_vector.hpp>
+#include <types.hpp>
+#include <util/random.hpp>
+
+#include "../test_environment.hpp"
+#include "../util/util.hpp"
+
+using namespace scaluq;
+
+TEST(CircuitTest, Optimize1) {
+    Random random;
+    for ([[maybe_unused]] UINT _ : std::views::iota(0, 100)) {
+        UINT n = 5;
+        UINT param_id = 0;
+        Circuit circuit(n);
+        for ([[maybe_unused]] UINT _ : std::views::iota(0, 1000)) {
+            UINT kind = random.int32() % 33;
+            UINT q0 = random.int32() % n;
+            UINT q1 = random.int32() % (n - 1);
+            if (q1 == q0) q1 = n - 1;
+            double r0 = random.uniform() * PI() * 2;
+            double r1 = random.uniform() * PI() * 2;
+            double r2 = random.uniform() * PI() * 2;
+            auto random_pauli = [&] {
+                std::vector<UINT> target_list;
+                std::vector<UINT> pauli_id_list;
+                for (UINT q : std::views::iota(0ULL, n)) {
+                    if (random.int32() & 1) {
+                        target_list.push_back(q);
+                        pauli_id_list.push_back(random.int32() & 3);
+                    }
+                }
+                return PauliOperator(target_list, pauli_id_list, random.uniform());
+            };
+            auto gen_param = [&] { return std::to_string(param_id++); };
+            switch (kind) {
+                // Applying non-unitary gate causes precision issue, so P0 and P1 is skipped.
+                case 0:
+                    circuit.add_gate(gate::I());
+                    break;
+                case 1:
+                    circuit.add_gate(gate::GlobalPhase(random.uniform() * PI() * 2));
+                    break;
+                case 2:
+                    circuit.add_gate(gate::X(q0));
+                    break;
+                case 3:
+                    circuit.add_gate(gate::Y(q0));
+                    break;
+                case 4:
+                    circuit.add_gate(gate::Z(q0));
+                    break;
+                case 5:
+                    circuit.add_gate(gate::H(q0));
+                    break;
+                case 6:
+                    circuit.add_gate(gate::S(q0));
+                    break;
+                case 7:
+                    circuit.add_gate(gate::Sdag(q0));
+                    break;
+                case 8:
+                    circuit.add_gate(gate::T(q0));
+                    break;
+                case 9:
+                    circuit.add_gate(gate::Tdag(q0));
+                    break;
+                case 10:
+                    circuit.add_gate(gate::SqrtX(q0));
+                    break;
+                case 11:
+                    circuit.add_gate(gate::SqrtXdag(q0));
+                    break;
+                case 12:
+                    circuit.add_gate(gate::SqrtY(q0));
+                    break;
+                case 13:
+                    circuit.add_gate(gate::SqrtYdag(q0));
+                    break;
+                case 14:
+                    circuit.add_gate(gate::RX(q0, r0));
+                    break;
+                case 15:
+                    circuit.add_gate(gate::RY(q0, r0));
+                    break;
+                case 16:
+                    circuit.add_gate(gate::RZ(q0, r0));
+                    break;
+                case 17:
+                    circuit.add_gate(gate::U1(q0, r0));
+                    break;
+                case 18:
+                    circuit.add_gate(gate::U2(q0, r0, r1));
+                    break;
+                case 19:
+                    circuit.add_gate(gate::U3(q0, r0, r1, r2));
+                    break;
+                case 20:
+                    circuit.add_gate(gate::CX(q0, q1));
+                    break;
+                case 21:
+                    circuit.add_gate(gate::CZ(q0, q1));
+                    break;
+                case 22:
+                    circuit.add_gate(gate::Swap(q0, q1));
+                    break;
+                case 23: {
+                    UINT block_size = random.int32() % (n / 2) + 1;
+                    UINT q0 = random.int32() % (n - block_size * 2 + 1);
+                    UINT q1 = random.int32() % (n - block_size * 2 + 1);
+                    if (q0 > q1) std::swap(q0, q1);
+                    q1 += block_size;
+                    circuit.add_gate(gate::FusedSwap(q0, q1, block_size));
+                } break;
+                case 24: {
+                    ComplexMatrix mat = get_eigen_matrix_random_single_qubit_unitary();
+                    circuit.add_gate(
+                        gate::OneQubitMatrix(q0,
+                                             {std::array{Complex(mat(0, 0)), Complex(mat(0, 1))},
+                                              std::array{Complex(mat(1, 0)), Complex(mat(1, 1))}}));
+                } break;
+                case 25: {
+                    ComplexMatrix mat1 = get_eigen_matrix_random_single_qubit_unitary();
+                    ComplexMatrix mat2 = get_eigen_matrix_random_single_qubit_unitary();
+                    ComplexMatrix mat = internal::kronecker_product(
+                        mat1, mat2);  // This is not fully random two-qubit unitary matrix, but its
+                                      // generation is complex so currently kronecker product of two
+                                      // one-qubit matrix is adopted.
+                    circuit.add_gate(gate::TwoQubitMatrix(q0,
+                                                          q1,
+                                                          {std::array{Complex(mat(0, 0)),
+                                                                      Complex(mat(0, 1)),
+                                                                      Complex(mat(0, 2)),
+                                                                      Complex(mat(0, 3))},
+                                                           std::array{Complex(mat(1, 0)),
+                                                                      Complex(mat(1, 1)),
+                                                                      Complex(mat(1, 2)),
+                                                                      Complex(mat(1, 3))},
+                                                           std::array{Complex(mat(2, 0)),
+                                                                      Complex(mat(2, 1)),
+                                                                      Complex(mat(2, 2)),
+                                                                      Complex(mat(2, 3))},
+                                                           std::array{Complex(mat(3, 0)),
+                                                                      Complex(mat(3, 1)),
+                                                                      Complex(mat(3, 2)),
+                                                                      Complex(mat(3, 3))}}));
+                } break;
+                case 26:
+                    circuit.add_gate(gate::Pauli(random_pauli()));
+                    break;
+                case 27:
+                    circuit.add_gate(gate::PauliRotation(random_pauli(), r0));
+                    break;
+                case 28:
+                    circuit.add_param_gate(gate::PRX(q0, random.uniform()), gen_param());
+                    break;
+                case 29:
+                    circuit.add_param_gate(gate::PRY(q0, random.uniform()), gen_param());
+                    break;
+                case 30:
+                    circuit.add_param_gate(gate::PRZ(q0, random.uniform()), gen_param());
+                    break;
+                case 31:
+                    circuit.add_param_gate(gate::PPauliRotation(random_pauli(), random.uniform()),
+                                           gen_param());
+                    break;
+            }
+        }
+        std::map<std::string, double> params;
+        for (UINT pid : std::views::iota(0ULL, param_id)) {
+            params[std::to_string(pid)] = random.uniform() * PI() * 2;
+        }
+        auto state0 = StateVector::Haar_random_state(n);
+        auto state1 = state0.copy();
+        circuit.update_quantum_state(state0, params);
+        UINT ngates = circuit.gate_count();
+        circuit.optimize();
+        circuit.update_quantum_state(state1, params);
+        ASSERT_LT(circuit.gate_count(), ngates);
+        ASSERT_TRUE(same_state(state0, state1));
+    }
+}

tests/gate/gate_test.cpp

@@ -417,7 +417,7 @@ TEST(GateTest, ApplyProbablisticGate) {
     auto probgate = gate::Probablistic({.1, .9}, {gate::X(0), gate::I()});
     UINT x_cnt = 0, i_cnt = 0;
     StateVector state(1);
-    for (auto _ : std::views::iota(0, 100)) {
+    for ([[maybe_unused]] auto _ : std::views::iota(0, 100)) {
         UINT before = state.sampling(1)[0];
         probgate->update_quantum_state(state);
         UINT after = state.sampling(1)[0];

tests/gate/param_gate_test.cpp

@@ -101,7 +101,7 @@ TEST(ParamGateTest, ApplyPProbablisticGate) {
     auto probgate = gate::PProbablistic({.1, .9}, {gate::PRX(0), gate::I()});
     UINT x_cnt = 0, i_cnt = 0;
     StateVector state(1);
-    for (auto _ : std::views::iota(0, 100)) {
+    for ([[maybe_unused]] auto _ : std::views::iota(0, 100)) {
         UINT before = state.sampling(1)[0];
         probgate->update_quantum_state(state, scaluq::PI());
         UINT after = state.sampling(1)[0];