Merge pull request qulacs#219 from Qulacs-Osaka/backprop_fix

backprop_fix

pysrc/qulacs/__init__.pyi

@@ -524,6 +524,10 @@ class ParametricQuantumCircuit(QuantumCircuit):
         """
         do backprop
         """
+    def backprop_inner_product(self, state: QuantumState) -> typing.List[float]: 
+        """
+        do backprop with innder product
+        """
     def copy(self) -> ParametricQuantumCircuit: 
         """
         Create copied instance

python/cppsim_wrapper.cpp

@@ -607,6 +607,7 @@ PYBIND11_MODULE(qulacs_core, m) {
         .def("add_parametric_multi_Pauli_rotation_gate", &ParametricQuantumCircuit::add_parametric_multi_Pauli_rotation_gate, "Add parametric multi-qubit Pauli rotation gate", py::arg("index_list"), py::arg("pauli_ids"), py::arg("angle"))
 
         .def("backprop",&ParametricQuantumCircuit::backprop,"do backprop",py::arg("obs"))
+        .def("backprop_inner_product",&ParametricQuantumCircuit::backprop_inner_product,"do backprop with innder product",py::arg("state"))
 
         .def("__repr__", [](const ParametricQuantumCircuit &p) {return p.to_string(); });
     ;

src/vqcsim/parametric_circuit.cpp

@@ -1,3 +1,4 @@
+#define _USE_MATH_DEFINES
 #include "parametric_circuit.hpp"
 
 #include <iostream>
@@ -168,66 +169,107 @@ void ParametricQuantumCircuit::add_parametric_multi_Pauli_rotation_gate(
         gate::ParametricPauliRotation(target, pauli_id, initial_angle));
 }
 
-// watle made
-using namespace std;
-std::vector<double> ParametricQuantumCircuit::backprop(
-    GeneralQuantumOperator* obs) {
+std::vector<double> ParametricQuantumCircuit::backprop_inner_product(
+    QuantumState* bistate) {
+    // circuitを実行した状態とbistateの、inner_productを取った結果を「値」として、それを逆誤差伝搬します
+    // bistateはノルムが1のやつでなくてもよい
     int n = this->qubit_count;
     QuantumState* state = new QuantumState(n);
+    //これは、ゲートを前から適用したときの状態を示す
     state->set_zero_state();
-    this->update_quantum_state(state);
-    // parametric bibunti tasu
-    std::vector<CPPCTYPE> bibun(1 << (this->qubit_count));
-    QuantumState* bistate = new QuantumState(n);
-    QuantumState* Astate = new QuantumState(n);  // for itizi work
-
-    obs->apply_to_state(Astate, *state, bistate);
-    bistate->multiply_coef(-1);
+    this->update_quantum_state(state);  //一度最後までする
 
-    double ansnorm = bistate->get_squared_norm();
-    if (ansnorm == 0) {
-        vector<double> ans(this->get_parameter_count());
-        return ans;
-    }
-    bistate->normalize(ansnorm);
-    ansnorm = sqrt(ansnorm);
-    int m = this->gate_list.size();
-    vector<int> gyapgp(m, -1);  // prametric gate position no gyaku
+    int num_gates = this->gate_list.size();
+    std::vector<int> inverse_parametric_gate_position(num_gates, -1);
     for (UINT i = 0; i < this->get_parameter_count(); i++) {
-        gyapgp[this->_parametric_gate_position[i]] = i;
+        inverse_parametric_gate_position[this->_parametric_gate_position[i]] =
+            i;
     }
-    vector<double> ans(this->get_parameter_count());
-    for (int i = m - 1; i >= 0; i--) {
-        auto gate = (this->gate_list[i])->copy();
-        if (gyapgp[i] != -1) {
-            Astate->load(bistate);
-            if (gate->get_name() != "ParametricRX" &&
-                gate->get_name() != "ParametricRY" &&
-                gate->get_name() != "ParametricRZ") {
-                std::cerr << "Error: " << gate->get_name()
-                          << " does not support backprop in parametric"
-                          << std::endl;
+    std::vector<double> ans(this->get_parameter_count());
+
+    /*
+    現在、2番のゲートを見ているとする
+    ゲート 0 1 2 3 4 5
+         state | bistate
+    前から2番までのゲートを適用した状態がstate
+    最後の微分値から逆算して3番まで　gateの逆行列を掛けたのがbistate
+
+    1番まで掛けて、YのΘ微分した行列を掛けたやつと、　bistateの内積の実数部分をとれば答えが出ることが知られている(知られてないかも)
+
+    ParametricR? の微分値を計算した行列は、Θに180°を足した行列/2 と等しい
+
+    だから、2番まで掛けて、 R?(π) を掛けたやつと、bistateの内積を取る
+
+    さらに、見るゲートは逆順である。
+    だから、最初にstateを最後までやって、　ゲートを進めるたびにstateに逆行列を掛けている
+    さらに、bistateが複素共役になっていることを忘れると、bistateに転置行列を掛ける必要がある。
+    しかしこのプログラムではbistateはずっと複素共役なので、　転置して共役な行列を掛ける必要がある。
+    ユニタリ性より、転置して共役な行列 = 逆行列
+    なので、両社にadjoint_gateを掛けている
+    */
+    QuantumState* Astate = new QuantumState(n);  //一時的なやつ
+    for (int i = num_gates - 1; i >= 0; i--) {
+        QuantumGateBase* gate_now = this->gate_list[i];  // sono gate
+        if (inverse_parametric_gate_position[i] != -1) {
+            Astate->load(state);
+            QuantumGateBase* RcPI;
+            if (gate_now->get_name() == "ParametricRX") {
+                RcPI = gate::RX(gate_now->get_target_index_list()[0], M_PI);
+            } else if (gate_now->get_name() == "ParametricRY") {
+                RcPI = gate::RY(gate_now->get_target_index_list()[0], M_PI);
+            } else if (gate_now->get_name() == "ParametricRZ") {
+                RcPI = gate::RZ(gate_now->get_target_index_list()[0],
+                    M_PI);  // 本当はここで2で割りたいけど、行列を割るのは実装が面倒
             } else {
-                double kaku = this->get_parameter(gyapgp[i]);
-                this->set_parameter(gyapgp[i], 3.14159265358979);
-                auto Dgate = (this->gate_list[i])->copy();
-                Dgate->update_quantum_state(Astate);
-                ans[gyapgp[i]] =
-                    (state::inner_product(state, Astate) * ansnorm).real();
-                this->set_parameter(gyapgp[i], kaku);
+                std::stringstream error_message_stream;
+                error_message_stream
+                    << "Error: " << gate_now->get_name()
+                    << " does not support backprop in parametric";
+                throw std::invalid_argument(error_message_stream.str());
             }
+            RcPI->update_quantum_state(Astate);
+            ans[inverse_parametric_gate_position[i]] =
+                state::inner_product(bistate, Astate).real() /
+                2.0;  //だからここで2で割る
+            delete RcPI;
         }
-
-        auto Agate = gate::get_adjoint_gate(gate);
+        auto Agate = gate::get_adjoint_gate(gate_now);
         Agate->update_quantum_state(bistate);
         Agate->update_quantum_state(state);
         delete Agate;
-        delete gate;
     }
     delete Astate;
     delete state;
-    delete bistate;
+    return ans;
+}  // CPP
+
+std::vector<double> ParametricQuantumCircuit::backprop(
+    GeneralQuantumOperator* obs) {
+    //オブザーバブルから、　最終段階での微分値を求めて、backprop_from_stateに流す関数
+    //上側から来た変動量 * 下側の対応する微分値 =
+    //最終的な変動量　になるようにする。
+
+    int n = this->qubit_count;
+    QuantumState* state = new QuantumState(n);
+    state->set_zero_state();
+    this->update_quantum_state(state);  //一度最後までする
+    QuantumState* bistate = new QuantumState(n);
+    QuantumState* Astate = new QuantumState(n);  //一時的なやつ　
 
+    obs->apply_to_state(Astate, *state, bistate);
+    bistate->multiply_coef(2);
+    /*一度stateを最後まで求めてから、さらにapply_to_state している。
+    なぜなら、　量子のオブザーバブルは普通の機械学習と違って、　二乗した値の絶対値が観測値になる。
+    二乗の絶対値を微分したやつと、　値の複素共役*2は等しい
+
+    オブザーバブルよくわからないけど、テストしたらできてた
+    */
+
+    //ニューラルネットワークのbackpropにおける、後ろからの微分値的な役目を果たす
+    auto ans = backprop_inner_product(bistate);
+    delete bistate;
+    delete state;
+    delete Astate;
     return ans;
-    // CPP
-}
+
+}  // CPP

src/vqcsim/parametric_circuit.hpp

@@ -3,6 +3,7 @@
 #include <cppsim/circuit.hpp>
 
 #include "cppsim/observable.hpp"
+#include "cppsim/state.hpp"
 class QuantumGate_SingleParameter;
 
 class DllExport ParametricQuantumCircuit : public QuantumCircuit {
@@ -49,4 +50,5 @@ class DllExport ParametricQuantumCircuit : public QuantumCircuit {
         std::vector<UINT> target, std::vector<UINT> pauli_id,
         double initial_angle);
     virtual std::vector<double> backprop(GeneralQuantumOperator* obs);
+    virtual std::vector<double> backprop_inner_product(QuantumState* bistate);
 };

test/vqcsim/test_backprop.cpp

@@ -29,25 +29,17 @@ TEST(Backprop, BackpropCircuit) {
     kairo.add_gate(gate::H(1));
     kairo.add_parametric_RX_gate(1, -1);
     kairo.add_parametric_RY_gate(2, 1);
-    //回路に適当にゲートを加える
-    // cout<<kairo<<endl;
-    QuantumState state(3);
-    state.set_zero_state();
-    kairo.update_quantum_state(&state);
-    // cout<<state<<endl;
+
     Observable observable(3);
-    observable.add_operator(1, "X 0");
+    observable.add_operator(1, "X 0 Z 2");
     observable.add_operator(1.2, "Y 1");
     observable.add_operator(1.5, "Z 2");
-    // cout<<observable.get_expectation_value(&state)<<endl<<endl;
+
     vector<double> kaku = {2.2, 0, 1.4, 1, -1, 1, 1, -1, 1};
-    // vector<double> kaku={2.2,1.4};
     GradCalculator wrakln;
     auto bibun = wrakln.calculate_grad(kairo, observable, kaku);
-    // for(auto it:bibun){cout<<it<<endl;}
-    // cout<<"de"<<endl;
     // culculate_gradした後は、パラメータがぐちゃぐちゃになるので、再セット
-    // a
+
     kairo.set_parameter(0, 2.2);
     kairo.set_parameter(1, 0);
     kairo.set_parameter(2, 1.4);
@@ -62,10 +54,53 @@ TEST(Backprop, BackpropCircuit) {
         cerr << bk[i] << " " << bibun[i].real() << endl;
         ASSERT_NEAR(bk[i], bibun[i].real(), 1e-10);
     }
-    // for(auto it:bk){cout<<it<<endl;}
-    // cout<<"de"<<endl;
-    // for(auto it:bibun){cout<<it<<endl;}
-    // cout<<"de"<<endl;
-    // cout<<observable.get_expectation_value(&state)<<endl<<endl;
-    // cerr<<(*gate::H(1))<<endl;
 }
+
+TEST(Backprop, BackpropCircuitInpro) {
+    ParametricQuantumCircuit kairo(3);
+    kairo.add_parametric_RX_gate(0, 2.2);
+    kairo.add_parametric_RY_gate(1, 0);
+    kairo.add_gate(gate::CNOT(0, 2));
+    kairo.add_parametric_RZ_gate(2, 1.4);
+    kairo.add_gate(gate::H(1));
+    kairo.add_parametric_RY_gate(0, 1);
+    kairo.add_gate(gate::CNOT(1, 0));
+    kairo.add_gate(gate::H(1));
+    kairo.add_parametric_RZ_gate(1, -1);
+    kairo.add_gate(gate::H(0));
+    kairo.add_gate(gate::CNOT(2, 0));
+    kairo.add_parametric_RX_gate(2, 1);
+    kairo.add_gate(gate::CNOT(1, 0));
+    kairo.add_parametric_RZ_gate(0, 1);
+    kairo.add_gate(gate::CNOT(0, 1));
+    kairo.add_gate(gate::H(1));
+    kairo.add_parametric_RX_gate(1, -1);
+    kairo.add_parametric_RY_gate(2, 1);
+    //回路に適当にゲートを加える
+
+    std::vector<CPPCTYPE> state_hai = {
+        1.0, 0.5, 3.0, -0.2, -2.0, 1.0, 0.7, 3.0};
+    QuantumState state_soku(3);
+    state_soku.load(state_hai);
+
+    QuantumState Astate(3);
+
+    auto bk = kairo.backprop_inner_product(&state_soku);
+    state_soku.load(state_hai);
+    vector<double> kaku = {2.2, 0, 1.4, 1, -1, 1, 1, -1, 1};
+
+    Astate.set_zero_state();
+    kairo.update_quantum_state(&Astate);
+    CPPCTYPE mto_sco = state::inner_product(&state_soku, &Astate);
+    for (int h = 0; h < 9; h++) {
+        kairo.set_parameter(h, kaku[h] + 0.0001);
+        Astate.set_zero_state();
+        kairo.update_quantum_state(&Astate);
+        CPPCTYPE gen_sco = state::inner_product(&state_soku, &Astate);
+
+        cerr << (gen_sco - mto_sco) * 10000.0 << " " << bk[h] << endl;
+        ASSERT_NEAR(((gen_sco - mto_sco) * 10000.0).real(), bk[h], 1e-2);
+
+        kairo.set_parameter(h, kaku[h]);
+    }
+}