1. はじめに
量子コンピュータが話題にのぼることが多くなりました。機械学習の世界も例外ではなく量子コンピュータを用いた機械学習の実装も多く考えられているようです。
本稿は量子コンピュータの中でもゲート計算方式と呼ばれるものを用いた「量子パーセプトロン」の論文を紐解きその実装に挑戦します。
元の論文は「Y. Cao et al. (2017) Quantum neuron: an
elementary building block for machine learning on quantum computers
」https://arxiv.org/abs/1711.11240 です。
その他に以下の日本語の記事を参考にさせていただいております。
- https://blueqat.com/yuichiro_minato2/d56b78ac-0ab4-4310-a5a4-6b1dca6b18cb
- https://qiita.com/piyo7/items/2104fe7084c95ed4b97b
特に2番目の記事は丁寧に計算を紐解いており、大変参考になりました。
2. パーセプトロン
パーセプトロンとはニューロンと呼ばれる神経回路モデルをもとにRosenblatt(1962)により開発された機械学習モデルです。ニューラルネットワークから現在のDeepLearningに至るまで、このモデルを基礎としています。
2.1 パーセプトロン定義
パーセプトロンは線形識別モデルの一つです。
入力ベクトル$\mathbf{x}$、重みパラメータ$\mathbf{w}$、バイアスパラメータ$b$に対して
f(\mathbf{x}) = \mathbf{w}^t\mathbf{x} + b
と定義し2クラス識別ルールを、$f(\mathbf{x})\geq0のとき1$、$f(\mathbf{x})<0のとき-1$としたモデルをパーセプトロンと呼びます。
2.2 線形分離可能
2クラスD次元入力データのクラスを完全に分離するD-1次元超平面が存在するとき線形分離可能と呼びます。線形分離可能なデータが与えられた場合、ある手続きを繰り返すことでパーセプトロンが入力データを全て正しく識別する具体的な$\mathbf{w},b$が得られることが保証されています(パーセプトロンの収束定理)。
2.3 多層パーセプトロン
一般的に与えられるデータが線形分離可能であることは稀です。そのため複数のパーセプトロンの出力を再度入力としてパーセプトロンを構築する多層パーセプトロンというものが考えられました。これを用いると線形分離不可能な入力データを場合によっては完全に分類することが可能となります。
他の線形識別モデルであるサポートベクトルマシンではカーネル関数を用いて線形分離不可能データに対処しています。
3. 量子コンピューティング
量子力学の原理を用いた計算を量子コンピューティングと一般的に呼ぶようです。計算方式であるため量子コンピュータの実機ではなく古典コンピュータでのシミュレートする場合もこのように言われることが多いです。
3.1 量子コンピューティングの種類
計算方式には量子ゲート方式、量子アニーリング方式があります。量子ゲート方式は量子ゲートと呼ばれる演算子を幾層にも重ねて回路を作成し計算結果を得ようとするものです。量子アニーリング方式は設定したハミルトニアンの最小エネルギー(実際は最小に近い場所)を実現する固有状態を得るアルゴリズムです。横磁場をかけた状態から段々横磁場を小さく、目的ハミルトニアンの寄与を大きく変化させていくことでその基底状態を得ます。組み合わせ最適化問題などへの応用が研究されています。
その他にも量子アニーリングと類似した断熱計算方式やクラスター状態という初期状態を用意し観測により演算結果を得る測定型量子計算と言われるものもあります。量子アニーリングと断熱計算の違い(私は区別がついていなかったのですが)についてはこちらの記事が参考になります。https://qiita.com/Kashalpha/items/9337c4f9fe4fbbb636fe
量子ゲート方式は古典コンピュータ上でのシミュレータが多数開発されています。本稿でもその一つであるIBM Qiskitを使用します。
3.2 量子機械学習
量子コンピューティングを用いた機械学習を量子機械学習と一般的に呼びます。しかしあまり厳密な定義はないようです。パラメータの更新などの部分を古典コンピュータに任せる量子古典ハイブリッド方式もこの中に含めます。
4. 量子パーセプトロン
元の論文を紐解いていきます。なお元の論文ではneuronとなっていますが、パーセプトロンと置き換えて問題ないようですので、本稿ではパーセプトロンとします。
4.1 構造
元の論文では、Repeat until success (RUS) 回路というものを使ってRy(シグモイド)の作用を計算しています。入力の線形結合に応じて、Ryの回転角が大きくなりシグモイド関数のような働きをするところが肝です。
- 図の引用元 「Y. Cao et al. (2017) Quantum neuron: an
elementary building block for machine learning on quantum computers
」https://arxiv.org/abs/1711.11240 ]
4.2 Repeat until success (RUS) 回路
量子回路計算のひとつです。アンシラを観測し0が観測されればOK、1の場合それを破棄してやり直す手法です。成功するまで繰り返す量子ゲートを用いて、シグモイド関数のような形をした非線形関数を得ます。σを欲しい関数として、以下のような結果を得る回路を構成、|0>が観測されればOK。|1>であればやり直し、これを成功するまで繰り返します。
$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$
A\ket{0}\otimes R_{y}(\sigma(\sum_{i=1}^{n}x_{i}+b))\ket{0}+B\ket{1}\otimes\ket{other}
A^2+B^2=1
4.3 計算
4.3.1 step1 線形結合を計算
制御$R_{y}$ゲートを利用し計算する。入力が1であれば重み$w_{i}$だけ加算され、0であれば何も作用しないことで線形結合計算します。この場合入力が${0,1}$の2値であることが仮定されている。
4.3.2 step2 活性化関数を計算
入力の線形結合を引数とする活性化関数を計算します。4.2でいうところの$\sigma$ですここの計算はhttps://qiita.com/piyo7/items/2104fe7084c95ed4b97b をほとんど真似させていただきました。
図のInput Registerはstep1で計算されており、その線形結合$2\theta$のY回転からスタートします。
\ket{00} \xrightarrow{qubit1\quad R_{y}(2\theta)} (\cos\theta I-i\sin\theta Y)\ket{0}\otimes\ket{0} = (\cos\theta\ket{0}+\sin\theta\ket{1}\otimes \ket{0} \\
\xrightarrow{Controll:-iY}\cos\theta\ket{00}+\sin\theta \ket{1} \otimes(-iY)\ket{0} =\cos\theta\ket{00}-i\sin\theta\ket{1}\otimes Y\ket{0} \\
\xrightarrow{qubit1\quad R_{y}^{\dagger}(2\theta)} \cos\theta \begin{pmatrix}
\cos\theta & \sin\theta \\
-\sin\theta & \cos\theta
\end{pmatrix} (1,0)^{t}\otimes\ket{0}-i\sin\theta \begin{pmatrix}
\cos\theta & \sin\theta \\
-\sin\theta & \cos\theta
\end{pmatrix} (0,1)^{t}\otimes Y\ket{0} \\
=\cos\theta(\cos\theta\ket{0}-\sin\theta\ket{1})\otimes\ket{0}-i\sin\theta(\sin\theta\ket{0}+\cos\theta\ket{1})\otimes Y\ket{0} \\
= (\cos^2\theta\ket{0}-\sin\theta\cos\theta\ket{1})\otimes\ket{0}-(i\sin^2\theta\ket{0}+i\sin\theta\cos\theta\ket{1})\otimes Y\ket{0} \\
=\cos^2\theta\ket{00}-\sin\theta\cos\theta\ket{10}-i\sin^2\theta\ket{0}\otimes Y\ket{0} + i\sin\theta\cos\theta\ket{1}\otimes Y\ket{0} \\
= \ket{0}\otimes(\cos^2\theta\ket{0}-i\sin^2\theta Y \ket{0})-\sin\theta\cos\theta\ket{1}(\ket{0} + iY\ket{0})
ここで
\cos^2\theta = \sqrt{\frac{\cos^4 \theta+\sin^4 \theta}{1+\tan^4 \theta}} \\
\sin^2\theta = \sqrt{\frac{\cos^4 \theta+\sin^4 \theta}{(1+\tan^4 \theta)/\tan^4\theta}} \\
(I+iY) = \sqrt{2}R_{y}(\frac{\pi}{2})
を利用すると与式は
\sqrt{\cos^4 \theta+\sin^4 \theta}\ket{0}(\frac{1}{\sqrt{1+\tan^4\theta}}\ket{0}-i\frac{\tan^2\theta}{\sqrt{1+\tan^4\theta}}Y\ket{0})-\frac{\sqrt{2}}{2}\sin2\theta\ket{1}\otimes R_{y}(\frac{\pi}{2})\ket{0} \\
=\sqrt{\frac{\cos^2\theta+1}{2}}\ket{0} R_{y}(2\arctan(\tan^2\theta))\ket{0}-\frac{\sqrt{2}}{2}\sin2\theta\ket{1}\otimes R_{y}(\frac{\pi}{2})\ket{0}
最後に$\arctan$が出てくる仕組みは
\frac{1}{\sqrt{1+\tan^4\theta}}\ket{0}-i\frac{\tan^2\theta}{\sqrt{1+\tan^4\theta}}Y\ket{0} \\
= \begin{pmatrix}
\frac{1}{\sqrt{1+\tan^2\theta}} & \frac{-\tan^2\theta}{\sqrt{1+\tan^2\theta}} \\
\frac{\tan^2\theta}{\sqrt{1+\tan^2\theta}} & \frac{1}{\sqrt{1+\tan^2\theta}}
\end{pmatrix} \ket{0}
= \begin{pmatrix}
\cos\phi & -\sin\phi \\
\sin\phi & \cos \phi
\end{pmatrix} \ket{0} = R_{y}(2\phi)\ket{0}
$\phi$ について解くと
\tan\phi = \frac{\sin\phi}{\cos\phi}=\tan^2\theta \Leftrightarrow \phi = \arctan(\tan^2 \theta)
$R_{y}(2\arctan(\tan^2\theta))\ket{0}$の項は$\theta$の値に応じてY軸方向へ回転します。$R_{y}(\pi)$の時に$R_{y}(\pi)\ket{0}$=$\ket{1}$となり確実にラベル1が観測できます。(通常のシグモイド関数の右端に相当)。
4.3.3 k回繰り返す
アンシラビットで0が観測された場合、もう一つの(式内の二番目の量子ビット)ビットは$R_{y}(\arctan(\tan^2\theta))$だけ回転したもの入力として同じ処理を繰り返します。1が観測された場合はもう一度最初からやり直します。(元論文ではもとに戻す角度で回転させるとあるが、当方の理解が追い付いていないため、のちの実装では0が観測された場合のみで計算する)
k回繰り返した場合前項の$R_{y}(2\arctan(\tan^2\theta))\ket{0}$は$R_{y}(2\arctan(\tan^{2^k}\theta))\ket{0}$となり、より急勾配なシグモイド関数になります。
k=2、k=3、kが一般の場合の回路を以下に記します。
5. 実装
2022年現在、複数の量子スタートアップなどから量子シミュレーターがリリースされています。本稿ではIBMのQiskitを使用します。
Qiskitの使用方法はIBMのHPをご確認ください。
5.1 k=1の場合の回路を実装
必要なライブラリをimportします。
import numpy as np
from qiskit import ClassicalRegister, QuantumRegister, QuantumCircuit, execute, Aer, IBMQ, assemble
from qiskit.visualization import plot_histogram
import itertools
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import minimize
k=1の場合の量子回路を作成する関数を用意します。
# 入力ベクトルの初期化関数 0の場合そのまま1の場合Xゲートで反転
def encode(qc, x):
for i in range(len(x)):
if x[i]==1:
qc.x(i)
else:
pass
def neuron(w,b,x):
qr = QuantumRegister(n+2+k-1, name="q") # 量子ビットを準備
c0 = ClassicalRegister(1, name="c0") # 観測用の古典レジストを準備
c1 = ClassicalRegister(1, name="c1") # 観測用の古典レジストを準備
qc = QuantumCircuit(qr, c0, c1) # 量子回路を作成
#print (x)
encode(qc, x) # 入力ベクトル初期化
qc.barrier()
#qc.x(0)
#qc.x(1)
for i in range(n):
qc.cry(2*w[i],i,n) #Ry
qc.ry(2*b,n) #切片
qc.cy(n,n+1) #Cy
qc.rz(-np.pi/2,n) #-i
#転置共役
qc.ry(-2*b,n)
for i in reversed(range(n)):
qc.cry(-2*w[i],i,n)
qc.measure(n, 0)
qc.barrier()
qc.measure(n+1, 1)
return qc
量子回路を作成します。
# 回路作成
n = 4
k = 1
w = [0.1, 0.4, 0.7, 0.15] # ニューロンの重み(適当な値)
b = 0.1 # バイアス(適当な値)
x =[1,0,1,1] # 入力データ
qc = neuron(w,b,x)
qc.draw('mpl',fold = 40)
計算時には簡便のために$R_{z}(\frac{-\pi}{2})$下の回路を組み込んで-iYとしましたが、実装時にはそのままにしています。
この入力に対して出力を計算してみます。量子回路モデルの出力値は観測したい量子状態$\psi$に対してMが出現する確率$P(M) = \bra{\psi}\ket{M}\bra{M}\ket{\psi}$ (M=0 or 1)に基づいて分布します。実機の場合実際の量子状態を観測しますが、bit数が少ない場合シミュレータを用いて代数計算からこの分布を計算しサンプリングし結果を得ます。
実際に上記コードの入力で計算をシミュレートしてみましょう。
backend = Aer.get_backend('qasm_simulator') # シミュレータを選択
shots = 1024 # サンプリング回数を指定
results = execute(qc, backend=backend, shots=shots).result()
answer = results.get_counts() 結果の経験分布を計算
分布を確認します。
plot_histogram(answer)
repeat untill success法は0の場合を成功、それ以外はやり直しとするものでしたので、観測が0であるものに限って確率を計算します。つまり上の図から出力が1になる確率は$P(1|観測ビットが0) = P(10) /(P(10)+P(00)) = 0.583/(0.583+0.060) $ということになります。
ではこの回路に入力データを変えながら出力を計算してみます。
# 回路にループでデータを投入
temp = np.array(list(itertools.product(range(2), repeat=n)))
qcl = []
for x in temp:
qc = neuron(w,b,x)
qcl.append(qc)
l1 = [] # 格納リスト
for i in range(len(qcl)):
backend = Aer.get_backend('qasm_simulator')
shots = 1024
results = execute(qcl[i], backend=backend, shots=shots).result()
answer = results.get_counts()
counts = 0
z_counts = 0
for k,v in answer.items():
k = k.replace(' ',"")
if k[-1]=='0': #0を観測できたもの
z_counts += v
if k[0]=='1':
counts += v
else:
pass
print (np.dot(temp[i],w)+b,counts / shots, counts / z_counts)
l1.append([np.dot(temp[i],w)+b,counts / z_counts])
入力の線形結合を横軸に縦軸に$P(1)$としたグラフを書くと以下のようになります。
sigmoid関数っぽくなっているのがわかります。
5.2 理論上の出力と合致しているか？
セクション4の式より理論上の1の出現分布は以下の通りです。
p(1|w,b,x) = \bra{0}R_{y}(2\arctan^{2^k}\theta) ^\dagger\ket{1}\bra{1}R_{y}(2\arctan^{2^k}\theta) \ket{0}
この理論的な式から得られる分布と一致しているか確認してみます。
test = QuantumCircuit(1, 1)
theta =np.pi /4
test.ry(2*np.arctan(np.tan(theta)**2),0)
test.measure(0,0)
thetas = np.linspace(0,np.pi/2,100)
l = []
for theta in thetas:
test = QuantumCircuit(1, 1)
test.ry(2*np.arctan(np.tan(theta)**2),0)
test.measure(0,0)
backend = Aer.get_backend('qasm_simulator')
shots = 1024
results = execute(test, backend=backend, shots=shots).result()
answer = results.get_counts()
try:
l.append([theta,answer['1'] / shots])
except:
pass
plt.scatter(np.array(l1)[:,0],np.array(l1)[:,1],color='red',label='RUS method')
plt.scatter(np.array(l)[:,0],np.array(l)[:,1],label='Arctan',marker='.')
plt.legend()
大体一致していることがわかります。
5.3 学習
これまでは与えられた$(w,b)$に対して出力を計算してきました。学習データを与えて$(w,b)$を学習することを考えます。
5.3.1 対象データ
以下のようなデータに対してラベル(0 or 1)を識別するように学習します。余談ですが以下の入出力のルールはOR回路の論理式に対応しています。
| x1 | x2 | y |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 1 | 1 |
| 1 | 0 | 1 |
| 1 | 1 | 1 |
# 学習データ作成
train = np.array([[0, 0, 0],
[0, 1, 1],
[1, 0, 1],
[1, 1, 1]])
5.3.2 ロス関数を作成
ニューロン関数は5.1で作成したk=1の場合を使用します。入力は2次元なのでn=2とします。loss関数は$1-(2y-1)(2p(1|w,b,x)-1)$です。$2*y-1$は${0,1}→{-1,1}$と変換するためです。
def loss_f(w):
#w=[w1, w2, w3, w4]
b = w[-1] # バイアス
w = w[0:n] # 係数
qcl = []
y = []
# 学習データ1個ずつに対して回路を作成 (ここらへん重ね合わせでできそうな気がするが...)
for x in train: # 学習データ(外部から与える)
# print (w,b)
qc = neuron(w,b,x[0:n]) # 回路作成
qcl.append(qc)
#print (x[-1])
y.append(x[-1]*2-1) #{0,1}→{-1,1}に変換
loss = []
for i in range(len(qcl)): # 学習データに対応した回路ごとにループして処理
backend = Aer.get_backend('qasm_simulator')
shots = 1024
results = execute(qcl[i], backend=backend, shots=shots).result()
answer = results.get_counts()
counts = 0
z_counts = 0
for k1,v in answer.items():
k1 = k1.replace(' ',"")
# print (k1)
if np.sum(int(k1[1:]))==0: #0を観測できたもの
z_counts += v
if k1[0]=='1':
counts += v
else:
pass
#print (np.dot(temp[i],w)+b,counts / shots, counts / z_counts)
# loss:=1-y*P(1,w,b,x) 真の値に近ければlossは小さくなる
loss.append(1-y[i]*((counts / (z_counts + 0.01))*2-1)) # avoid zero divide
losses = np.sum(loss)
return losses
5.3.3. 学習
パラメータを初期化します。乱数で初期化したいところですが、一旦以下のようにします。
w = [0.1,0.9]
b = 0.1
w = w + [b]
print (w)
学習はloss関数が小さくなるようにパラメータ(w,b)を更新するように行います。量子コンピューティングといいつつここの最適化計算は古典コンピューティングソルバーで行います。量子コンピューティングの世界で有名なQAOAなども一部の計算は古典コンピュータに投げています。(ハイブリッド方式とか言うようです。)
Scipyのminimizeオプティマイザーを使用します。元の論文ではnelder mead法を使用していたががうまくいかなかったためPowell法というもので行いました。
Nfeval = 1
def callbackF(w):
global Nfeval
f = loss_f(w)
print ('{0:4d} {1: 3.6f} {2: 3.6f} {3: 3.6f} {3: 3.6f}'.format(Nfeval, w[0], w[1] ,w[2], f))
Nfeval += 1
plt.scatter(Nfeval, f, color='black')
res = minimize(loss_f, w, method="Powell", callback=callbackF) # 関数の戻り値が一つ
得らえた解で線形分離超平面を描画します。通常の超平面と異なり+$\pi/4$とするのは この点が通常のパーセプトロンの閾値に該当し、超平面は2次元だと$ax+by=\pi/4$と記述されるからです。
# 超平面描画関数
def line(w,x):
y=-x*w[0]/w[1]-w[2]+np.pi/4
return (y)
# 学習前パラメータ
w = [0.1,0.9]
b = 0.1
w = w + [b]
# 学習前パラメータの超平面
plt.figure(figsize=(7,4))
x = np.linspace(-0.5,1.5)
y = line(w,x)
plt.plot(x,y,label='init separating hyperplane')
plt.scatter(train[:,0][train[:,2]==1],train[:,1][train[:,2]==1],label='Train label Postive',color='red')
plt.scatter(train[:,0][train[:,2]==0],train[:,1][train[:,2]==0],label='Train label Negative')
# 学習後パラメータの超平面
x = np.linspace(-0.5,1.5)
y = line(res.x,x)
plt.plot(x,y,label='Trained separating hyperplane')
plt.legend()
上手く分類できています。
6. 今後の展開
量子ニューラルネットワークと称される論文はいくつか発表されていますが、それぞれが独自の方式をとっているものが多く古典のように統一的な理論が整備されているわけではないのが現状です。
今回ご紹介した論文のように量子パーセプトロンに着目した研究は、非線形な関係を量子ニューラルネットワークに導入するにはどうするかという観点から注目されることが多いようです。
個人的な興味としては、この関数を部品に多層化を目指し線形分離不可能なものに拡張するなどがあります。ほかにもパーセプトロンの識別容量の理論的な解析なども面白そうです。
ほかに読んでみたい論文は、https://arxiv.org/abs/1803.00745 です。これはパラメータシフトルールという画期的な微分システムを開発し深層学習ライクな回路を構成しています。
では！