289 add backprop subbranch

code-spell-checker-dictionary.txt

@@ -1,6 +1,7 @@
 arnoldi
 autoapi
 autoremove
+backprop
 BFLOAT
 braket
 ccache
@@ -38,8 +39,8 @@ GROUPNAME
 gtest
 habs
 Hadamard
-hessenberg
 hendrikmuhs
+hessenberg
 IBMQ
 IBMQ's
 icircuit

include/scaluq/circuit/circuit.hpp

@@ -7,6 +7,7 @@
 
 #include "../gate/gate.hpp"
 #include "../gate/param_gate.hpp"
+#include "../operator/operator.hpp"
 #include "../types.hpp"
 
 namespace scaluq {
@@ -101,6 +102,16 @@ class Circuit {
         const std::map<std::string, double>& parameters = {},
         std::uint64_t seed = 0) const;
 
+    template <ExecutionSpace Space>
+    std::unordered_map<std::string, double> backprop_inner_product(
+        StateVector<Prec, Space>& state,
+        StateVector<Prec, Space>& bistate,
+        const std::map<std::string, double>& parameters);
+
+    template <ExecutionSpace Space>
+    std::unordered_map<std::string, double> backprop(
+        const Operator<Prec, Space>& observable, const std::map<std::string, double>& parameters);
+
 private:
     std::uint64_t _n_qubits;
 
@@ -229,6 +240,17 @@ void bind_circuit_circuit_hpp(nb::module_& m) {
             "parameters"_a = std::map<std::string, double>{},
             "seed"_a = std::nullopt,
             "Simulate noise circuit. Return all the possible states and their counts.")
+        .def("backprop_inner_product",
+             &Circuit<Prec>::template backprop_inner_product<ExecutionSpace::Host>,
+             "state"_a,
+             "bistate"_a,
+             "parameters"_a,
+             "Compute gradients of inner product between state and bistate using back propagation.")
+        .def("backprop",
+             &Circuit<Prec>::template backprop<ExecutionSpace::Host>,
+             "observable"_a,
+             "parameters"_a,
+             "Compute gradients of expectation value of observable using back propagation.")
         .def(
             "update_quantum_state",
             [&](const Circuit<Prec>& circuit,
@@ -306,6 +328,17 @@ void bind_circuit_circuit_hpp(nb::module_& m) {
             "parameters"_a = std::map<std::string, double>{},
             "seed"_a = std::nullopt,
             "Simulate noise circuit. Return all the possible states and their counts.")
+        .def("backprop_inner_product",
+             &Circuit<Prec>::template backprop_inner_product<ExecutionSpace::HostSerial>,
+             "state"_a,
+             "bistate"_a,
+             "parameters"_a,
+             "Compute gradients of inner product between state and bistate using back propagation.")
+        .def("backprop",
+             &Circuit<Prec>::template backprop<ExecutionSpace::HostSerial>,
+             "observable"_a,
+             "parameters"_a,
+             "Compute gradients of expectation value of observable using back propagation.")
 #ifdef SCALUQ_USE_CUDA
         .def(
             "update_quantum_state",
@@ -383,6 +416,17 @@ void bind_circuit_circuit_hpp(nb::module_& m) {
             "parameters"_a = std::map<std::string, double>{},
             "seed"_a = std::nullopt,
             "Simulate noise circuit. Return all the possible states and their counts.")
+        .def("backprop_inner_product",
+             &Circuit<Prec>::template backprop_inner_product<ExecutionSpace::Default>,
+             "state"_a,
+             "bistate"_a,
+             "parameters"_a,
+             "Compute gradients of inner product between state and bistate using back propagation.")
+        .def("backprop",
+             &Circuit<Prec>::template backprop<ExecutionSpace::Default>,
+             "observable"_a,
+             "parameters"_a,
+             "Compute gradients of expectation value of observable using back propagation.")
 #endif  // SCALUQ_USE_CUDA
         .def("copy",
              &Circuit<Prec>::copy,

include/scaluq/operator/apply_pauli.hpp

@@ -220,6 +220,8 @@ void apply_pauli_rotation(std::uint64_t control_mask,
                           Complex<Prec> coef,
                           Float<Prec> angle,
                           StateVector<Prec, Space>& state_vector) {
+    const auto raw = state_vector._raw;
+    const std::uint64_t dim = state_vector.dim();
     std::uint64_t global_phase_90_rot_count = std::popcount(bit_flip_mask & phase_flip_mask);
     Complex<Prec> true_angle = angle * coef;
     Complex<Prec> half_angle = true_angle / Float<Prec>{2};
@@ -230,43 +232,40 @@ void apply_pauli_rotation(std::uint64_t control_mask,
         const Complex<Prec> cval_pls = cosval + Complex<Prec>(0, 1) * sinval;
         Kokkos::parallel_for(
             "apply_pauli_rotation",
-            Kokkos::RangePolicy<SpaceType<Space>>(
-                0, state_vector.dim() >> std::popcount(control_mask)),
+            Kokkos::RangePolicy<SpaceType<Space>>(0, dim >> std::popcount(control_mask)),
             KOKKOS_LAMBDA(std::uint64_t i) {
                 std::uint64_t state_idx =
                     insert_zero_at_mask_positions(i, control_mask) | control_value_mask;
                 if (Kokkos::popcount(state_idx & phase_flip_mask) & 1) {
-                    state_vector._raw[state_idx] *= cval_min;
+                    raw(state_idx) *= cval_min;
                 } else {
-                    state_vector._raw[state_idx] *= cval_pls;
+                    raw(state_idx) *= cval_pls;
                 }
             });
         return;
     }
     std::uint64_t pivot = Kokkos::bit_width(bit_flip_mask) - 1;
     Kokkos::parallel_for(
         "apply_pauli_rotation",
-        Kokkos::RangePolicy<SpaceType<Space>>(
-            0, state_vector.dim() >> (std::popcount(control_mask) + 1)),
+        Kokkos::RangePolicy<SpaceType<Space>>(0, dim >> (std::popcount(control_mask) + 1)),
         KOKKOS_LAMBDA(std::uint64_t i) {
             std::uint64_t basis_0 =
-                internal::insert_zero_at_mask_positions(i, control_mask | 1ULL << pivot) |
-                control_value_mask;
+                insert_zero_at_mask_positions(i, control_mask | 1ULL << pivot) | control_value_mask;
             std::uint64_t basis_1 = basis_0 ^ bit_flip_mask;
 
             int bit_parity_0 = Kokkos::popcount(basis_0 & phase_flip_mask) & 1;
             int bit_parity_1 = Kokkos::popcount(basis_1 & phase_flip_mask) & 1;
 
             // fetch values
-            Complex<Prec> cval_0 = state_vector._raw[basis_0];
-            Complex<Prec> cval_1 = state_vector._raw[basis_1];
+            Complex<Prec> cval_0 = raw(basis_0);
+            Complex<Prec> cval_1 = raw(basis_1);
 
             // set values
-            state_vector._raw[basis_0] =
+            raw(basis_0) =
                 cosval * cval_0 +
                 Complex<Prec>(0, 1) * sinval * cval_1 *
                     PHASE_M90ROT<Prec>()[(global_phase_90_rot_count + bit_parity_0 * 2) % 4];
-            state_vector._raw[basis_1] =
+            raw(basis_1) =
                 cosval * cval_1 +
                 Complex<Prec>(0, 1) * sinval * cval_0 *
                     PHASE_M90ROT<Prec>()[(global_phase_90_rot_count + bit_parity_1 * 2) % 4];
@@ -309,8 +308,7 @@ void apply_pauli_rotation(std::uint64_t control_mask,
             {0, 0}, {states.batch_size(), states.dim() >> (std::popcount(control_mask) + 1)}),
         KOKKOS_LAMBDA(const std::uint64_t batch_id, const std::uint64_t i) {
             std::uint64_t basis_0 =
-                internal::insert_zero_at_mask_positions(i, control_mask | 1ULL << pivot) |
-                control_value_mask;
+                insert_zero_at_mask_positions(i, control_mask | 1ULL << pivot) | control_value_mask;
             std::uint64_t basis_1 = basis_0 ^ bit_flip_mask;
 
             int bit_parity_0 = Kokkos::popcount(basis_0 & phase_flip_mask) & 1;
@@ -372,9 +370,9 @@ void apply_pauli_rotation(std::uint64_t control_mask,
                     Kokkos::TeamThreadRange(team,
                                             states.dim() >> (std::popcount(control_mask) + 1)),
                     [&](const std::uint64_t i) {
-                        std::uint64_t basis_0 = internal::insert_zero_at_mask_positions(
-                                                    i, control_mask | 1ULL << pivot) |
-                                                control_value_mask;
+                        std::uint64_t basis_0 =
+                            insert_zero_at_mask_positions(i, control_mask | 1ULL << pivot) |
+                            control_value_mask;
                         std::uint64_t basis_1 = basis_0 ^ bit_flip_mask;
                         int bit_parity_0 = Kokkos::popcount(basis_0 & phase_flip_mask) & 1;
                         int bit_parity_1 = Kokkos::popcount(basis_1 & phase_flip_mask) & 1;

include/scaluq/operator/operator.hpp

@@ -46,6 +46,8 @@ class Operator {
 
     [[nodiscard]] StdComplex get_expectation_value(
         const StateVector<Prec, Space>& state_vector) const;
+    [[nodiscard]] std::vector<StdComplex> get_expectation_values(
+        const StateVector<Prec, Space>& state_vector) const;
     [[nodiscard]] std::vector<StdComplex> get_expectation_value(
         const StateVectorBatched<Prec, Space>& states) const;
 
@@ -181,6 +183,13 @@ void bind_operator_operator_hpp(nb::module_& m) {
             },
             "state"_a,
             "Get the expectation value of the operator with respect to a state vector.")
+        .def(
+            "get_expectation_values",
+            [](const Operator<Prec, Space>& op, const StateVector<Prec, Space>& state) {
+                return op.get_expectation_values(state);
+            },
+            "state"_a,
+            "Get the expectation values of the pauli operators with respect to a state vector.")
         .def(
             "get_expectation_value",
             [](const Operator<Prec, Space>& op, const StateVectorBatched<Prec, Space>& states) {

include/scaluq/operator/pauli_operator.hpp

@@ -226,7 +226,7 @@ void bind_operator_pauli_operator_hpp(nb::module_& m) {
              "source"_a,
              "target"_a,
              "Get transition amplitude of measuring state vector. $\\bra{\\chi}P\\ket{\\psi}$.")
-#ifdef SCALUQ_ENABLE_CUDA
+#ifdef SCALUQ_USE_CUDA
         .def("get_expectation_value",
              nb::overload_cast<const StateVectorBatched<Prec, ExecutionSpace::Default>&>(
                  &PauliOperator<Prec>::template get_expectation_value<ExecutionSpace::Default>,

include/scaluq/util/utility.hpp

@@ -113,7 +113,7 @@ Complex<Prec> inner_product(const Kokkos::View<Complex<Prec>*, SpaceType<Space>>
     Kokkos::parallel_reduce(
         "inner_product",
         Kokkos::RangePolicy<SpaceType<Space>>(0, a.extent(0)),
-        KOKKOS_LAMBDA(std::uint64_t i, Complex<Prec> & lsum) {
+        KOKKOS_LAMBDA(std::uint64_t i, Complex<Prec>& lsum) {
             lsum += internal::conj(a(i)) * b(i);
         },
         result);
@@ -228,6 +228,20 @@ inline SparseComplexMatrix transform_sparse_matrix_by_order(
     return transform_dense_matrix_by_order(mat.toDense(), targets).sparseView();
 }
 
+template <Precision Prec>
+constexpr double get_epsilon() {
+    if constexpr (Prec == Precision::F16)
+        return 1.;
+    else if constexpr (Prec == Precision::F32)
+        return 1e-4;
+    else if constexpr (Prec == Precision::F64)
+        return 1e-12;
+    else if constexpr (Prec == Precision::BF16)
+        return 1.;
+    else
+        static_assert(internal::lazy_false_v<internal::Float<Prec>>, "unknown Precision");
+}
+
 }  // namespace internal
 
 }  // namespace scaluq

python/binding.cpp

@@ -11,6 +11,7 @@
 #include <nanobind/stl/string.h>
 #include <nanobind/stl/string_view.h>
 #include <nanobind/stl/tuple.h>
+#include <nanobind/stl/unordered_map.h>
 #include <nanobind/stl/variant.h>
 #include <nanobind/stl/vector.h>
 

src/circuit/circuit.cpp

@@ -2,6 +2,7 @@
 #include <scaluq/gate/gate_factory.hpp>
 #include <scaluq/gate/merge_gate.hpp>
 #include <scaluq/gate/param_gate_factory.hpp>
+#include <scaluq/util/utility.hpp>
 
 namespace scaluq {
 template <Precision Prec>
@@ -79,7 +80,7 @@ void Circuit<Prec>::update_quantum_state(StateVector<Prec, Space>& state,
         if (!parameters.contains(key)) {
             using namespace std::string_literals;
             throw std::runtime_error(
-                "Circuit::update_quantum_state(StateVector&, const std::map<std::string_view, double>&) const: parameter named "s +
+                "Circuit::update_quantum_state(StateVector&, const std::map<std::string, double>&) const: parameter named "s +
                 std::string(key) + "is not given.");
         }
     }
@@ -498,5 +499,120 @@ void Circuit<Prec>::check_gate_is_valid(const ParamGate<Prec>& gate) const {
     }
 }
 
+template <Precision Prec>
+template <ExecutionSpace Space>
+std::unordered_map<std::string, double> Circuit<Prec>::backprop_inner_product(
+    StateVector<Prec, Space>& state,
+    StateVector<Prec, Space>& bistate,
+    const std::map<std::string, double>& parameters) {
+    const auto eps = internal::get_epsilon<Prec>();
+    const std::uint64_t n_qubits = this->n_qubits();
+    const std::uint64_t n_gates = this->n_gates();
+
+    std::unordered_map<std::string, double> gradients;
+    const auto& key_set = this->key_set();
+    gradients.reserve(key_set.size());
+    for (const auto& key : key_set) {
+        gradients.emplace(key, 0.0);
+    }
+
+    StateVector<Prec, Space> Astate(n_qubits);
+
+    for (std::int64_t i = static_cast<std::int64_t>(n_gates) - 1; i >= 0; i--) {
+        const auto& cur_gate = this->get_gate_at(static_cast<std::uint64_t>(i));
+
+        if (cur_gate.index() == 1) {
+            const auto& [pgate, key] = std::get<1>(cur_gate);
+
+            Astate = state.copy();
+            double pauli_coef = 1.0;
+            if (pgate.param_gate_type() == ParamGateType::ParamPauliRotation) {
+                internal::ParamGatePtr<internal::ParamPauliRotationGateImpl<Prec>> ppr(pgate);
+
+                const auto c = ppr->pauli().coef();
+                const double cre = static_cast<double>(std::real(c));
+                const double cim = static_cast<double>(std::imag(c));
+                if (std::abs(cim) >= eps) {
+                    throw std::runtime_error("backprop_inner_product: pauli coef must be real");
+                }
+                if (std::abs(cre) < eps) {
+                    throw std::runtime_error("backprop_inner_product: pauli coef must be nonzero");
+                }
+                pauli_coef = cre;
+            }
+            const double scale = pgate->param_coef() * pauli_coef;
+            if (std::abs(scale) < eps) {
+                throw std::runtime_error(
+                    "backprop_inner_product: param_coef or pauli coef is zero");
+            }
+            pgate->update_quantum_state(Astate, -M_PI / scale);
+            const auto ip = internal::inner_product<Prec, Space>(bistate._raw, Astate._raw);
+            const double contrib = -scale * static_cast<double>(ip.real());
+
+            gradients[key] += contrib;
+        }
+
+        if (cur_gate.index() == 0) {
+            const auto& g = std::get<0>(cur_gate);
+
+            auto inv = g->get_inverse();
+            inv->update_quantum_state(bistate);
+            inv->update_quantum_state(state);
+        } else {
+            const auto& [pgate, key] = std::get<1>(cur_gate);
+
+            auto inv = pgate->get_inverse();
+            const auto it = parameters.find(key);
+            if (it == parameters.end()) {
+                throw std::runtime_error("backprop_inner_product: missing parameter for key=" +
+                                         key);
+            }
+            const auto param = it->second;
+
+            inv->update_quantum_state(bistate, param);
+            inv->update_quantum_state(state, param);
+        }
+    }
+    return gradients;
+}
+template std::unordered_map<std::string, double> Circuit<internal::Prec>::backprop_inner_product(
+    StateVector<internal::Prec, ExecutionSpace::Host>& state,
+    StateVector<internal::Prec, ExecutionSpace::Host>& bistate,
+    const std::map<std::string, double>& parameters);
+template std::unordered_map<std::string, double> Circuit<internal::Prec>::backprop_inner_product(
+    StateVector<internal::Prec, ExecutionSpace::HostSerial>& state,
+    StateVector<internal::Prec, ExecutionSpace::HostSerial>& bistate,
+    const std::map<std::string, double>& parameters);
+#ifdef SCALUQ_USE_CUDA
+template std::unordered_map<std::string, double> Circuit<internal::Prec>::backprop_inner_product(
+    StateVector<internal::Prec, ExecutionSpace::Default>& state,
+    StateVector<internal::Prec, ExecutionSpace::Default>& bistate,
+    const std::map<std::string, double>& parameters);
+#endif  // SCALUQ_USE_CUDA
+
+template <Precision Prec>
+template <ExecutionSpace Space>
+std::unordered_map<std::string, double> Circuit<Prec>::backprop(
+    const Operator<Prec, Space>& obs, const std::map<std::string, double>& parameters) {
+    const std::uint64_t n_qubits = this->n_qubits();
+    StateVector<Prec, Space> state(n_qubits);
+    this->update_quantum_state(state, parameters);
+
+    StateVector<Prec, Space> bistate = state.copy();
+    obs.apply_to_state(bistate);
+    return backprop_inner_product(state, bistate, parameters);
+}
+template std::unordered_map<std::string, double> Circuit<internal::Prec>::backprop<
+    ExecutionSpace::Host>(const Operator<internal::Prec, ExecutionSpace::Host>& obs,
+                          const std::map<std::string, double>& parameters);
+template std::unordered_map<std::string, double> Circuit<internal::Prec>::backprop<
+    ExecutionSpace::HostSerial>(const Operator<internal::Prec, ExecutionSpace::HostSerial>& obs,
+                                const std::map<std::string, double>& parameters);
+#ifdef SCALUQ_USE_CUDA
+template std::unordered_map<std::string, double> Circuit<internal::Prec>::backprop<
+    ExecutionSpace::Default>(const Operator<internal::Prec, ExecutionSpace::Default>& obs,
+                             const std::map<std::string, double>& parameters);
+#endif  // SCALUQ_USE_CUDA
+
 template class Circuit<internal::Prec>;
 }  // namespace scaluq