誤差逆伝播法と量子Fisher行列とデータセットについて

大阪大学のQIQBで量子コンピュータの研究開発をやっております森と申します。

今年の3月末に公開された国産初号機である理研の量子コンピュータの開発チームの一員でもあります。
（ソフトウェア技術者が足りないので入れてもらいました。）

去年も[量子コンピューターAdvent Calendar 2022 19日目](https://qiita.com/forest1040/items/5ad07b972650948a20f9)を書いたのですが、今年も量子機械学習関係で書いてみたいと思います。

去年も書いた誤差逆伝播法の説明をもう一度するのと、量子Fisher行列のことと、大阪大学で量子コンピュータ用のデータセットを作成しましたので、その事を書いてみようかと思います。
あとは、求人情報とかも載せるので興味がある方は、見てみてください。

# 誤差逆伝播法

今年も大阪大学の情報学科系の競プロ勢（プログラミングバイト）にお世話になりました。
プログラミングバイトの方に、誤差逆伝播法のロジックをPythonで書いてもらったので、それをちょっと整理して、再度説明しようかなと思います。
コードは[Qulacs](https://github.com/qulacs/qulacs)という大阪大学の競プロ勢を中心に開発保守している量子回路シミュレータを使っています。
今新しい量子回路シミュレータを競プロ勢と開発中なのですが、その話はまた今度。

```python
import math
from typing import List

from qulacs import Observable, ParametricQuantumCircuit, QuantumState, gate
from qulacs.state import inner_product


def python_backprop(circ: ParametricQuantumCircuit, obs: Observable) -> List[float]:
    n = circ.get_qubit_count()
    state = QuantumState(n)
    state.set_zero_state()
    circ.update_quantum_state(state)
    obs_state = QuantumState(n)
    work = QuantumState(n)
    obs.apply_to_state(work, state, obs_state)
    obs_state.multiply_coef(2)

    num_gates = circ.get_gate_count()
    inverse_parametric_gate_position = [-1] * num_gates
    for i in range(circ.get_parameter_count()):
        inverse_parametric_gate_position[circ.get_parametric_gate_position(i)] = i
    grad = [0.0] * circ.get_parameter_count()

    # パラメータゲート適用用
    temp_state = QuantumState(n)
    for i in range(num_gates - 1, -1, -1):
        gate_now = circ.get_gate(i)
        if inverse_parametric_gate_position[i] != -1:
            temp_state.load(state)
            # 回転ゲートの場合、Θの微分を求める
            # 微分は、Θに180°を足した行列/2 を適用すればよい
            # （行列の1/2を計算するのは面倒なので、最後に割る）
            if gate_now.get_name() == "ParametricRX":
                rcpi = gate.RX(gate_now.get_target_index_list()[0], math.pi)
            elif gate_now.get_name() == "ParametricRY":
                rcpi = gate.RY(gate_now.get_target_index_list()[0], math.pi)
            elif gate_now.get_name() == "ParametricRZ":
                rcpi = gate.RZ(gate_now.get_target_index_list()[0], math.pi)
            else:
                raise RuntimeError()
            rcpi.update_quantum_state(temp_state)
            # obs_stateと回転ゲートの微分を適用した状態の内積を取る
            grad[inverse_parametric_gate_position[i]] = inner_product(
                obs_state, temp_state
            ).real / 2.0 # ここで1/2をする
        # 逆行列をかけて１つ前のゲートに戻る
        inv_gate = gate_now.get_inverse()
        inv_gate.update_quantum_state(obs_state)
        inv_gate.update_quantum_state(state)
    return grad

```

`python_backprop()`が誤差逆伝播法で指定した量子回路の各パラメータの勾配を計算して、返してくれる関数です。
`ParametricQuantumCircuit`という、パラメータを埋め込むことができる形式の量子回路と、`Observable`を指定します。
`Observable`とは言ってしまえば行列なのですが、量子の世界の言葉でいうところの物理量です。
ハミルトニアンやコスト関数をPauli行列を指定する形で定義します。
ハミルトニアンとは？みたいな気持ちになりますがエネルギーの事と思ってもらえればと思います。物理屋は何でもエネルギーという抽象的な形で物事を扱って、美しい数式で世界を記述します。（多分）

量子回路を`Observable`を適用した量子状態と`Observable`を適用しない量子状態を用意して、最後のゲートから順に微分値を求めて、適用し内積を取ります。
回転ゲートの微分は、Θに180°を足した行列/2 を適用すればよいです。
あとは内積を取るだけです。

計算式は、[量子コンピューターAdvent Calendar 2022 19日目](https://qiita.com/forest1040/items/5ad07b972650948a20f9)も参照して見てください。

使用例は、以下のような感じになります。

```python
from qulacs import ParametricQuantumCircuit, GradCalculator, Observable
from backprop import python_backprop

n = 2
observable = Observable(n)
observable.add_operator(1.0, "X 0")
# observable.add_operator(1.0, "Z 0")
# observable.add_operator(1.0, "I 0 I 1")
circuit = ParametricQuantumCircuit(n)

theta = [2.2, 1.4, 0.8]
circuit.add_parametric_RX_gate(0, theta[0])
circuit.add_parametric_RY_gate(0, theta[1])
circuit.add_parametric_RZ_gate(0, theta[2])

# GradCalculatorの場合
gcalc = GradCalculator()
print(gcalc.calculate_grad(circuit, observable))
print(python_backprop(circuit, observable))

```

`Observable`を定義して、`ParametricQuantumCircuit`に回転ゲートを設定します。
`GradCalculator`はパラメータシフト法で勾配計算した値になります。
比較してみてください。

# 量子Fisher行列
量子Fisher行列とは、量子回路のモデルの形を表す行列になります。普通に学習させるよりも、速く最適解を見つけられるという自然勾配に使われたりします。
実際には、量子Fisher行列の計算にもコストがかかるので、モデルによっては、素朴に誤差逆伝播法で計算してしまった方が速い場合もある気がしています。
最適解を見つけられない場合にトライしてみてください。

大阪大学の藤井研に[https://github.com/quest-kit/QuEST](QuEST)という量子回路シミュレータを作っているTyson Jonesが来ていて、Qulacsにも量子Fisher行列を入れたらという話になって、ちょっとやってみました。

アルゴリズム的には、Tysonの[EFFICIENT CLASSICAL CALCULATION OF THE
QUANTUM NATURAL GRADIENT](https://arxiv.org/pdf/2011.02991.pdf)をそのまま実装しています。Qulacsは回転角の向きが一般の定義と逆なので注意が必要です。

```python
# https://arxiv.org/pdf/2011.02991.pdf 3:Algorithm
from typing import List
import math
import numpy as np
from qulacs import (
    QuantumCircuit,
    ParametricQuantumCircuit,
    QuantumState,
    gate,
)
from qulacs.state import inner_product

# 微分値の計算（素朴に計算している）
def get_differential_gate(g, theta):
    def _differential_gate(gate_matrix):
        return (
            -1 * math.sin(theta / 2) / 2 * np.array([[1, 0], [0, 1]])
            + -1 * -1.0j * math.cos(theta / 2) / 2 * gate_matrix
            # qulacsの回転角の方向が逆なので、-1をかけている
            # Because the direction of the quilacs rotation is reversed, multiplying it by -1.
            # + 1 * -1.0j * math.cos(theta / 2) / 2 * gate_matrix
        )

    if g.get_name() == "ParametricRX":
        matrix = _differential_gate(np.array([[0, 1], [1, 0]]))
        rcpi = gate.DenseMatrix(g.get_target_index_list()[0], matrix)
    elif g.get_name() == "ParametricRY":
        matrix = _differential_gate(np.array([[0, -1.0j], [1.0j, 0]]))
        rcpi = gate.DenseMatrix(g.get_target_index_list()[0], matrix)
    elif g.get_name() == "ParametricRZ":
        matrix = _differential_gate(np.array([[1, 0], [0, -1]]))
        rcpi = gate.DenseMatrix(g.get_target_index_list()[0], matrix)
    else:
        raise RuntimeError()
    return rcpi


def fisher(
    input_circuit: QuantumCircuit, ansatz: ParametricQuantumCircuit, theta: List[float]
):
    n = input_circuit.get_qubit_count()
    chi = QuantumState(n)
    input_circuit.update_quantum_state(chi)
    chi_string = "|in>"
    phi = chi.copy()
    phi_string = chi_string
    gate = ansatz.get_gate(0)
    gate.update_quantum_state(chi)
    chi_string = "U0" + chi_string
    psi = chi.copy()
    psi_string = chi_string
    rcpi = get_differential_gate(gate, theta[0])
    rcpi.update_quantum_state(phi)
    phi_string = "dU0" + phi_string

    num_param = ansatz.get_gate_count()
    # print(f"num_param: {num_param}")

    T = np.zeros(num_param, dtype=complex)
    T[0] = inner_product(chi, phi)
    L = np.zeros((num_param, num_param), dtype=complex)
    L[0][0] = inner_product(phi, phi)

    for j in range(1, num_param):
        print(f"j:{j}")
        lambda_state = psi.copy()
        lambda_string = psi_string
        print(f"lmd: {lambda_string}")
        phi = psi.copy()
        phi_string = psi_string
        # print(f"phi: {phi_string}")
        gate = ansatz.get_gate(j)
        rcpi = get_differential_gate(gate, theta[j])
        rcpi.update_quantum_state(phi)
        phi_string = f"dU{j}" + phi_string
        print(f"phi: {phi_string}")
        L[j][j] = inner_product(phi, phi)
        for i in range(j - 1, -1, -1):
            print(f"i:{i}")
            gate = ansatz.get_gate(i + 1).get_inverse()
            gate.update_quantum_state(phi)
            phi_string = f"U{i+1}^" + phi_string
            print(f"phi: {phi_string}")
            gate = ansatz.get_gate(i).get_inverse()
            gate.update_quantum_state(lambda_state)
            lambda_string = f"U{i}^" + lambda_string
            print(f"lmd: {lambda_string}")
            myu = lambda_state.copy()
            myu_string = lambda_string
            gate = ansatz.get_gate(i)
            rcpi = get_differential_gate(gate, theta[i])
            rcpi.update_quantum_state(myu)
            myu_string = f"dU{i}" + myu_string
            print(f"myu: {myu_string}")
            L[i][j] = inner_product(myu, phi)
            # print(f"<myu|phi>: {myu_string[::-1] + phi_string}")

        T[j] = inner_product(chi, phi)
        gate = ansatz.get_gate(j)
        gate.update_quantum_state(psi)
        psi_string = f"U{j}" + psi_string

    print(f"T: {T}")
    print(f"L: {L}")

    qfi = np.zeros((num_param, num_param))
    for i in range(num_param):
        for j in range(num_param):
            if i <= j:
                qfi[i][j] = L[i][j] - T[i].conj() * T[j]
            else:
                qfi[i][j] = L[j][i].conj() - T[i].conj() * T[j]
    return qfi


def main():
    n_qubit = 2
    x = 0.27392337

    theta = [
        4.002148315014479,
        1.6951199159934145,
        0.25744424357926954,
        0.10384619671527331,
        5.109927617709579,
        5.735012432197602,
    ]

    input_circuit = QuantumCircuit(n_qubit)

    for i in range(n_qubit):
        input_circuit.add_RY_gate(i, np.arcsin(x) * 2)
        input_circuit.add_RZ_gate(i, np.arccos(x * x) * 2)

    ansatz = ParametricQuantumCircuit(n_qubit)
    ansatz.add_parametric_RX_gate(0, theta[0])
    ansatz.add_parametric_RZ_gate(0, theta[1])
    ansatz.add_parametric_RX_gate(0, theta[2])
    ansatz.add_parametric_RX_gate(1, theta[3])
    ansatz.add_parametric_RZ_gate(1, theta[4])
    ansatz.add_parametric_RX_gate(1, theta[5])

    print("theta:", theta)
    # print("QFI")
    qfi = fisher(input_circuit, ansatz, theta)
    row_size, col_size = qfi.shape
    for i in range(row_size):
        tmp = ""
        for j in range(col_size):
            tmp += str("{:.08f}, ".format(qfi[i][j]))
        print(tmp)

    # print(qfi[1][3], qfi[1][4], qfi[1][5])


main()

```

2つのパラメータの組み合わせの行列を作る感じになります。
逆行列で消し合って、対角要素は0になります。
`get_differential_gate()`は、回転ゲートの微分を計算しています。こっちでは素朴に計算しています。
ロジック的には、TとLを求めて、最後にqfi行列を作ります。
量子Fisher行列の式の定義とロジックを見比べると、式が言っている事がわかってくるような気がしてくると思いますので追いかけてみてください。

双対空間等、数学的なところは北野先生に教えて頂きました。
（北野先生はたくさんの量子人材を排出している先生です！）

# MNISQ(MNISTの量子回路データセット)
藤井研のドクターの学生が中心となって構築したMNISTを量子状態に埋め込んだデータセットです。
バイト軍団と共に私も少し手伝ったので紹介しておきます。
[FujiiLabCollaboration/MNISQ-quantum-circuit-dataset: MNISQ circuit dataset for machine learning and quantum machine learning](https://github.com/FujiiLabCollaboration/MNISQ-quantum-circuit-dataset)

MNISTの画像データを量子状態に埋め込みそれを量子回路で再現するために、[AQCE](https://arxiv.org/abs/2112.14524)というアルゴリズムで作成したデータセットです。
量子回路は、OpenQASM形式で作成しました。
（ちなみにAQCEは大阪大学QIQBの上田先生達が考案したアルゴリズムです。）

この量子回路を実行すると、MNIST画像が量子状態に埋め込まれます。
それを元に量子機械学習等に使ってもらえればと思っています。
[MNISQ-quantum-circuit-dataset/doc/source/notebooks/0_quickstart.ipynb at main · FujiiLabCollaboration/MNISQ-quantum-circuit-dataset](https://github.com/FujiiLabCollaboration/MNISQ-quantum-circuit-dataset/blob/main/doc/source/notebooks/0_quickstart.ipynb)にチュートリアルがあり、これを見て頂ければ大体使い方がわかると思います。

データセットは、クラウドで公開されており実行時にダウンロードされます。
結構、大量のデータセットを作成したので、藤井研のクラスタマシン（CPUコアが92x42個）あるやつが大活躍しました。
（活躍しすぎて部屋の温度が熱くなりすぎたのはナイショ）

PennyLaneのデータセットにも追加されました。[MNISQ](https://pennylane.ai/datasets/other/mnisq)

他のドクターと別の謎のデータセットも作成していて、そのうち公開できると思います。


# 求人情報的な何か

ということで、今年自分が関わってやってきた量子機械学習に関係するもので、発表してもよさそうなものを書いてみました。

大阪大学には、量子コンピュータの実機もおいてあり、学園祭で一般公開したりもしています。
[国産量子コンピュータの中身を一般公開！＠まちかね祭 － 大阪大学](https://www.osaka-u.ac.jp/ja/event/2023/11/10578)
私は、最近少しずつ実機の制御系にも手を出しつつあるのですが、実機を触るのはすごくロマンがあります。
私は、量子コンピュータのスゴイところは、以下のようなところだと思っています。
「ミクロの世界に織込まれた広大な高次元空間（ヒルベルト空間）で計算して、マクロな世界に結果を取り出せる（測定）というアナログコンピュータのパワーとデジタルの良いところを併せ持つ事ができる（少なくとも理論上は！）」
普通のアナログコンピュータは、パワフルですがエラー訂正できない。
しかし量子コンピュータはアナログチックでありながら、エラー訂正の理論が存在し、使いこなすことができれば、ものスゴイパワーを秘めていると思っています。
今はまだコンピュータシステムというよりは、実験装置的なものだと思った方がしっくりくる量子コンピュータですが、アーキテクチャから研究してコンピュータシステムとして組み上げる事に凄く夢を感じます。
（自分を見出してくれた先生との飲み会で、現状の量子コンピュータについて率直な意見を言ったら、すごく悲観的ですね！と言われてしまったので、ロマンを語ってみました。。）

量子制御は、[NTT技術ジャーナルの2023年9月号](https://journal.ntt.co.jp/backnumber/2023/vol3509) のP17: 「超伝導量子コンピュータのシステムの設計と開発」がわかりやすいと思います。私もNTT鈴木さんのソフトウェアを使わせて頂いております。

大阪大学では、そんな量子コンピュータを実機からクラウド開発や量子アルゴリズム、量子化学計算等の量子アプリケーションまで、上から下までフルスタックで研究開発することができます。
しかしながら、ソフトウェア技術者が足りないのです。
クラウド開発や量子回路シミュレータの開発には、そんなに量子コンピュータ自体の知識は必要ありません。
必要な事は教えてくれる先生方や学生がいっぱいいます。
もし興味を持っている方がいらっしゃいましたら、ご連絡ください。
ぜひ、私とマイクロ波という熱い想いを実機にぶつけましょう！