量子回路学習の実装

0. はじめに

本稿では、量子機械学習のフレームワークであるQuantum Circuit Learning(K. Mitarai et al. 2019)を紹介します。この論文は少々古いですが、量子機械学習のエッセンスが詰まっています。

元の論文に加えて、筆者らのHPを参考にしました。以下が参考文献です。

[1] 元論文 https://arxiv.org/pdf/1803.00745.pdf
[2] 筆者らの実装記事 https://dojo.qulacs.org/ja/latest/notebooks/5.2_Quantum_Circuit_Learning.html
[3] 量子ゲートの代数表現を参考にした https://utokyo-icepp.github.io/qc-workbook/dynamics_simulation.html

筆者らは独自のシミュレータであるQulacsを使用していますが、本稿ではIBMのQiskitを使用します。

1. 仕組み

まず、量子回路の構成を俯瞰し、その後学習アルゴリズムを説明します。

これ以降、nはqubitの数を表します。

1.1 回路の構成

3つのコンポーネントで構成されます。1つ目は入力データをエンコードし、2つ目はqubitを相互に干渉させ、3つ目はパラメータに応じた回転ゲートを作用させます。

1.1.1 入力データのエンコーディング

適当にデータを角度エンコーディングします。実装ではデータ$x$を$R_{y}(\arcsin(x)),R_{z}(\arccos(x^2))$を用いて以下のように角度エンコーディングさせています(例はn=3の場合です)。

1.1.2 ランダムハミルトニアンによる干渉

作用させるハミルトニアンを$H = \sum_{i<j}Z_iZ_j+\sum_{i} X_{i}$と定義します。$X,Z$はパウリオペレータです。添え字は、作用するqubitの位置を表します。私は勘違いしていましたが、この式に出てくるパウリ行列は通常の$2×2$の行列ではなく、$2^n×2^n$の行列です。たとえば、3つのqubitで1番目と2番目に作用する行列を正確に書くと$Z_{1}\otimes Z_{2}\otimes I$で、これは$8×8$のエルミート行列です。

シュレーディンガー方程式を解くと$\exp(-iHt)$が得られ、これはユニタリ行列です。（エルミート行列の指数はユニタリ行列であるため）tは時間を表しこの実装では適当な定数が与えられます。

参考ページでは対角化して実装していますが、この文書では、qiskitのハミルトニアンから直接$\exp(-iHt)$をゲートとして実装できる関数を使用しています。$\exp(-iHt)$をそのまま行列としゲートに埋め込んでいるため(おそらく3-qubit程度なのでできた？)、qiskitではその詳細は表示されないようです。

ここで中身がどのように相互作用しているかを考えます。元のハミルトニアンは行列の指数であるため、$\exp(a+b) = \exp(a)\exp(b)$は一般的には成立しません。そこで、鈴木・トロッター分解なるものを行うことで、大きな数$M$(トロッタ数)を用いて$H$を$\left( \prod_{i<j} \exp(-itZ_{i}Z_{j}/M)\prod_{i} \exp(-itX_{i}/M) \right)^M$と近似することができます。$\exp(-itZ_{i}Z_{j}/M)$は、一般的に$CNOT_{ij}R_{z}(2*t/M)j CNOT_{ij}$となることが知られています(参考:[3]や筆者らの元記事のhttps://dojo.qulacs.org/ja/latest/notebooks/4.2_trotter_decomposition.html の下の方)。つまり、上記の図中の2つ目のコンポーネントは、近似して中身を見てみると、2量子ビットに作用するユニタリ変換の積となっていると言えます。

1.1.3 回転角をパラメータとしたゲートで回転させる

最後のコンポーネントは、パラメータ付きの回転ゲートで構成されます。このパラメータを学習誤差が小さくなるように学習します。$i$番目のビットに対して、$R_{x,i}(θ_{i,1})R_{z,i}(θ_{i,2})R_{z,i}(θ_{i,3})$と定義します。$n=3$の場合、$\prod_{i=1}^{3} R_{x,i}(θ_{i,1})R_{z,i}(θ_{i,2})R_{z,i}(θ_{i,3})$となります。各回転ゲートも$2×2$ではなく、$2^n×2^n$の行列です。この時点で最適化すべきパラメータ数は$3×n$であることがわかります。

1.1.4 コンポーネントを繰り返す

元論文の筆者らの実装をまねして、エンコーディング以外の2つのコンポーネント群が3回繰り返されるようにします。そして1番上のビットのz方向の値を観測するために、その値を格納する古典レジスタを用意して、

とします。これで完成です。最適化するために動かすパラメータは$3×3×3=27$です。

1.2 学習

1.2.1 誤差関数

$Z$測定で一番上のbitの$Z$方向の物理量の期待値を$\hat{y}$とし、この値と教師データ$y$の2乗誤差を小さくするようにパラメータを学習します。$-1≤E[Z]≤1$であるので、一般的なデータの値に対応するために、$\hat{y}$に定数倍$a$を乗じ($-a≤E[Z]≤a$となる)、この$a$も学習します。

今回は学習データが$sin$関数なので$a$はパラメータではなく定数を与えます。$-1≤sinx≤1$なので$a=1$で良いと思うのですが、参考ページでは$a=2$を与えています。実際$a=1$だと学習が上手くいきません($y=-1,1$付近で誤差が大きいまま学習が終わってしまう)。学習が上手くいくようにヒューリスティックに定めた値なのでしょう。私は$2$でもうまくいかなかったので$2.5$で設定しました。

1.2.2 学習アルゴリズム

• 2乗誤差の最小化のために古典の勾配フリーソルバーを使用します。ただし元の論文ではパラメータθの微分を計算する方法が提示されています。ご興味ある方はそちらを参照願います。

2. 実装

python上で実行可能なIBMのqiskitで実装します。qubit=3、コンポーネントの繰り返しも3で行います。設定部分など[2]の記事のコードをそのまま使わせて頂いています。

python

# 必要ライブラリのインストール
import numpy as np
from qiskit import *
import matplotlib.pyplot as plt
from qiskit.quantum_info.operators import Operator
from qiskit import Aer
from qiskit.quantum_info.operators import Operator, Pauli
from qiskit.extensions import HamiltonianGate # Hからexp(-itH)を計算
# 最適化計算のため
from scipy.optimize import minimize

使用するqubit数や、コンポーネントの繰り返し数を設定します。また学習するサンプルデータを作成します。ここでは$\sin(x)$但し$(-1≤x≤1)$を学習します。定義域を$(-1≤x≤1)$としているのは、エンコーディング角度を入力の逆三角関数としているからです。より一般の定義域を持つ関数の場合、入力を適当な値でスケールすれば良さそうです。

python

# パラメータ設定
nqubit = 3 # qubitの数
c_depth = 3 # circuitの深さ
time_step = 0.77  # ハミルトニアンによる時間発展の経過時間

# [x_min, x_max]のうち, ランダムにnum_x_train個の点をとって教師データとする.
x_min = - 1.; x_max = 1.;
num_x_train = 50

# 学習したい1変数関数
func_to_learn = lambda x: np.sin(x*np.pi)

# 乱数シード設定
np.random.seed(0)
# 教師データ
x_train = x_min + (x_max - x_min) * np.random.rand(num_x_train)
y_train = func_to_learn(x_train)

# 現実のデータを用いる場合を想定し、きれいなsin関数にノイズを付加
mag_noise = 0.1 # 参考先よりやや大きいばらつきとする
y_train = y_train + mag_noise * np.random.randn(num_x_train)

plt.plot(x_train, y_train, "o")# ; plt.show()

エンコードするコンポーネントを作成します。

python

# xをエンコードするゲートを作成する関数
def U_in(x):
    cirq = QuantumCircuit(nqubit)

    angle_y = np.arcsin(x)
    angle_z = np.arccos(x**2)

    for i in range(nqubit):
        cirq.ry(angle_y, i)
        cirq.rz(angle_z, i)

    return cirq

ハミルトニアン部分を作成します。

python

# ハミルトニアンの部分を作成
np.random.seed(0)
ham = Operator(np.zeros((2**nqubit, 2**nqubit))) 

# l = []
for i in range(nqubit):
    init_x = list('I'*nqubit) # Xはnqubit通り 全部単位行列で初期化 I⊗I⊗I
    init_x[i] = 'X' # i番目だけXに変える I⊗X⊗I
    
    print ("".join(init_x))
    jx = -1. + 2.* np.random.rand() # ランダム磁場
    # print (jx)
    ham += jx * Operator.from_label("".join(init_x))
    
    for j in range(i+1, nqubit): # Zはnqubit_C_2通り
        init_z = list('I'*nqubit) # 全部単位行列で初期化 I⊗I⊗I
        # print (i, j)
        init_z[i] = 'Z' # i番目Z
        init_z[j] = 'Z' # j番目Z Z⊗Z⊗I
        print ("".join(init_z))
        J_ij = -1. + 2.* np.random.rand() # ランダム相互作用
        #print (J_ij)
        ham += J_ij * Operator.from_label("".join(init_z))

# H→exp(itH)
gate = HamiltonianGate(ham, time_step) 

回転ゲートも合わせて作成し、コンポーネント群を3回繰り返す関数を作成します。ハミルトニアン部分だけ関数ではなく、上記で作成したオブジェクトを使います(特に入力がないため)。

python

# パラメータを初期化
params_init = 2.0 * np.pi * np.random.rand(c_depth*nqubit*3) 
# 入力x_iからモデルの予測値y(x_i, theta)を返す関数
def qcl_pred(x, params):

    # 入力状態計算
    circ = U_in(x)

    # 出力状態計算
    for k in range(c_depth): # c_depth回繰り返す
        circ.append(gate, range(nqubit)) # あらかじめ定義しておいたハミルトニアンgate
        # 回転ゲート作成
        qc = QuantumCircuit(nqubit) # nqubitの回路作成
        pmt = params.reshape(c_depth, nqubit ,3) # c_depth×nqubit×3のパラメータ配列 
        for i in range(nqubit):
            qc.rx(pmt[k,i,0], i)
            qc.rz(pmt[k,i,1], i)
            qc.rx(pmt[k,i,2], i)
        circ.append(qc, range(nqubit))
    
    # 古典レジスタ
    cr = ClassicalRegister(1,'creg')
    circ.add_register(cr)
    circ.measure(0,0) # 0番目の状態をz方向に観測

    # モデルの出力
    emulator = Aer.get_backend('qasm_simulator')
    job = execute(circ, emulator, shots=1024 ) # 1024回観測
    hist = job.result().get_counts()
    # print (hist)
    # z方向の期待値計算
    if len(hist.keys())==2: # 0(or1)が100%の場合だけ処理かえる 要領悪いが...
        res = (1 * hist['0'] -1 * hist['1']) / (hist['0'] + hist['1'])
    else:
        res = ((int(list(hist.keys())[0])+1) % 2)*2-1
        # print (res)
    res = 2.5 * res # この2.5倍はヒューリスティックであるが重要            

    return res

optimizeに与える目的関数(2乗誤差関数)を定義します。

python

# cost function Lを計算
def cost_func(params):

    # num_x_train個のデータについて計算
    y_pred = [qcl_pred(x, params) for x in x_train]

    # quadratic loss
    L = ((y_pred - y_train)**2).mean()

    return L

scipyのoptimizeソルバーで解きます。

python

%%time
# 学習 2時間程度
opt_result = minimize(cost_func, params_init, method='Nelder-Mead')

学習したパラメータで$\sin$を再現できているか確認しましょう！

python

# パラメータthetaの初期値のもとでのグラフ
xlist = np.arange(x_min, x_max, 0.02)
y_init = [qcl_pred(x, params_init) for x in xlist]
plt.plot(xlist, y_init)

xlist = np.arange(x_min, x_max, 0.02)
y_init = [qcl_pred(x, opt_result.x) for x in xlist]
plt.plot(xlist, y_init)
plt.plot(x_train, y_train, "o")# ; plt.show()

# plt.plot(x,y)

青色の線が初期値で、オレンジの線が学習後のパラメータで描画したものです。$y=-1,1$に近いところでの誤差が依然として残っていますが、まあまあいい感じです。

3. 考察など

驚くべきは、入力50事例に対して27個ものパラメータを用いていながら、過剰な適合をしていないところです。しかしこれを非効率(多くのパラメータが必要)とみるか、汎化性能が高いとみるかは議論が残ります。

入力事例1つにつき量子回路を作成する必要があります。そのため量子重ね合わせの計算の利点が失わているとも言えます。この問題は他の量子回路機械学習にみられる問題です。

今度はクラス分類にも挑戦したいです。では！