連載記事の２回目 １，０ではない重ね合わせでの初期状態
いきなりこの回を読むと意味不明かもしれない。その場合は、１回目 を読んで頂きたい。本連載記事の目的と概要を投稿している。今回から具体的なBlueqatでの量子コンピュータの記事となる。ただ本題となるセルオートマトンの量子アルゴリズムやその時間発展は4回目以降の予定。今回の目玉は１，０ではない確率初期状態のセッティング方法とBlueqatで確率だの状態ベクトルだのをどう料理するのかである。まず、その前提としてq-camの構造とモジュールを紹介する。
#2-1. q-camの構造とモジュール
Githubに公開したq-camsの中にあるプログラムとモジュールの表を示す。尚、以降、q-camsと表記した部分は、Githubのリポジトリ―にリンクしている。尚、御覧のとおりファイル名の拡張子はpy、つまり言語はpythonである。
| ファイル名 | 内容 | 説明のある回 |
|---|---|---|
| qcam.py | q-cam専用の、初期確率セット ベクトルと確率の抽出、出力、プロットなどのモジュール |
２回目 |
| qmcn.py | Blueqatで汎用的に使用できるモジュール マルチ制御NOTゲート CCCX～CCCCCCX |
３回目 |
| probshot.py | セットした確率を観測でみる | ２回目 |
| probvect.py | セットした確率を状態ベクトルでみる | ２回目 |
| q-cam30r.py | ルール30 セルオートマトン | ４回目 |
| q-cam60r.py | ルール60 セルオートマトン | ４回目 |
| q-cam90r.py | ルール90 セルオートマトン | ４回目 |
| q-cam102r.py | ルール102 セルオートマトン | ４回目 |
| q-cam110r.py | ルール110 セルオートマトン | ４回目 |
| q-cam184r.py | ルール184 セルオートマトン ルール | ５回目 |
| q-cam184i.py | ルール184 セルオートマトン 虚数版 | ５回目 |
| q-cambbr.py | ソリトン（箱玉系）セルオートマトン | ６回目 |
| q-camgr.py | グローバーアルゴリズム | ７回目 |
| q-frqca102.py | ルール102の量子位相推定 | ８回目 |
| q-frqga2.py | 2bitグローバーの量子位相推定 | ８回目 |
| q-frqga3.py | 3bitグローバーの量子位相推定 | ８回目 |
| q-frqga4.py | 4bitグローバーの量子位相推定 | ８回目 |
| q-frqga5.py | 5bitグローバーの量子位相推定 | ８回目 |
| q-frqga6.py | 6bitグローバーの量子位相推定 | ８回目 |
q-camの構造と言っているのは上の表のq-cam30r.py以下のプログラムにおける共通構造のことである。今回紹介する「１，０ではない確率初期状態のセッティング方法」は、qcam.pyのモジュール内の関数であり、またBlueqatで確率だの状態ベクトルだのをどう料理するのかは、probshot.pyとprobvect.pyと言うプログラムを用いての説明となる。モジュールqcamとqmcnはプログラムと同じフォルダにダウンロードして使用することが前提。各、プログラムはコンソールから動かす。さて、q-camの構造は単純で以下の通りである。
①モジュールの読み込み部
②定数セット部
③Blueqatで記述した量子プログラム部(Q-process)
④メインの制御部
⑤終了部分（プロットなど）
順に実際のコード例で説明しよう。
①モジュールの読み込み部
from blueqat import Circuit
import numpy as np
import qcam
import qmcn as q
一行目はBlueqatを使う場合、必要な読み込みで詳しくはこちらを参照。https://blueqat.readthedocs.io/ja/latest/
二行目は、numpyを使用するため。３行目のqcamは、q-camプログラムで必須のモジュール。
import qcam as qc などとし略称でモジュールを使用できるようにしても、勿論OK。ラストの行にあるqmcnは、マルチ制御NOTゲートを使わない限りは不要で、q-camsのプログラム群の中ではグローバーアルゴリズムとソリントセルオートマトン以外では使用しない。
②定数セット部
N=10 # number of cells
R=2 # number of registers
# initial probability distribution
initial_a=np.array([0.8,0,0.2,0.1,1,0.9,0.5,0.1,0.2,0.9],dtype='float')
Nはセルオートマトンでのセル数（下図参照）、Rはレジスタ数で、この積NRが全量子ビット数となる。セルに見立てたN個の量子ビット(0～N-1)もレジスタに含まれ、これをレジスタ0とする。セル数Nはユーザーが好きに書き変えられるが、レジスタ数はセルオートマトンのアルゴリズムに依存する。詳しくは４回目に説明する。Blueqatの場合、量子ビット数の上限は32bitなので、NR≦32だが、24ぐらいを超えると計算時間が長くなる。また、配列initial_a は、まさに今回の「０，１ではない確率の集合」を現しており、左から順に、セルに割り当てらる。従ってこの配列の要素数は必ず、上記のセル数と一致しなければならない。下図参照。これらの確率を実際の量子ビットにセットするのが、１，０ではない確率初期状態のセッティング方法であって2-2で説明する。
上のコード内では省略したが定数セット部はさらに続くが、主として各レジスタのスタート量子ビット番号とラスト量子ビット番号の割り当てである。
③Blueqatで記述した量子プログラム部(Q-process)
この部分は、まさにBlueqatで量子アルゴリズムを記述する部分であり、内容は４回目に説明する。最も簡単な例として、ルール102の場合の量子プログラム部を示す。
def qproc():
for i in range(N):
c.cx[1+i, i]
return
④メインの制御部
メインの制御部の代表的例
pstep=0
ret='y'
while ret=='y':
pstep+=1
c=Circuit(N)
pinitial_a=probability_a
if pstep==1:
pinitial_a=initial_a
qcam.propinit(N,c,pinitial_a)
qproc()
master_a=np.array(c.run())
num,vector_a,prob_a=qcam.qvextract(N,R,1,master_a)
stdev,probability_a=qcam.qcalcd(N, num,vector_a,prob_a)
qcam.qcresultout(N, pstep, num,pinitial_a,vector_a,prob_a,stdev,probability_a)
stepdist_a[pstep,0:N]=probability_a[0:N]
stdev_a[pstep]=stdev
ret=input(' NEXT(Y/N)?')
pstepは時間発展（１回目 参照）のカウンタで初期を０としている。配列pinitial_aはその時間ステップのセルオートマトンの遷移前の初期状態で、時刻1、つまりpstep=1のときは、initial_aそのものである。まずはじめにpinitial_a配列の確率を量子ビットにセットする。それが2-2で説明するqcamモジュール内のpropinit関数である。これでセル番号0～N-1に見立てたreg0の量子ビット番号0～N-1に確率（確率振幅）がセットされる。セット後、qproc()で遷移のためのBlueqatの量子プロセスとなる。実際の実行はBlueqatのc.run()である。このBlueqatの命令は2-3で説明する状態ベクトルを見る方法である。状態ベクトルといいつつ実はBlueqatではベクトルそのものがこれで分かるわけではなく、その状態ベクトルの確率振幅がベクトルの数だけ、複素数表記の一次元リストの形で戻ってくる。それをmaster_aという一次元の配列に取り込み、これから状態ベクトルだの、確率だのを抽出し、計算をする。
ある時刻ステップ（pstep）のセルの確率の結果は、次の時刻ステップ（pstep+1）の初期状態となるように代入している。stepdist_a[pstep,0:N]=probability_a[0:N] コンソール上で NEXT(Y/N)?と入力を促し、ｙと入力した場合は、次の時間ステップに進む。この繰り返しで、１ステップづつ時間発展をさせていく。ｙ以外を入力した場合は以下の⑤に続く。
⑤終了部分（プロットなど）
時間発展の結果をまとめて表示する関数qcam.qcfinalとセルオートマトンの時間発展を視覚的に表現するための関数qcam.qcamplotを実行し終了となる。尚、これらの出力内容に関しては４回目で詳述する。
#2-2. 1,0ではない確率的初期状態（ry 又は rx 又は u3 ゲートで作る）
さて、いよいよ、セルに見立てた量子ビットに「０，１ではない確率の集合」をセットする方法である。ある量子ビット、ここではi番目の量子ビットをqubit[i]と表現しよう。またqubit[i] の状態を示すのを qubit[i]⇒ と表現しよう。確率ゼロをセットするのは簡単で、BlueqatのCircuitオブジェクトを作成後、何もしなければゼロ状態である。つまり qubit[i]⇒0 である。状態１をセットするのも簡単でBlueqatのXゲート(反転あるいはフリップゲートとも称する）を使えばqubit[i]⇒1となる。コードで記せば、x[qubit[i]]である。では、次に０と１がちょうど半々で重ねあった状態、つまり確率0.5、qubit[i]⇒0.5は？ これも簡単でアダマールゲートを使えば良い。コードで記せば、h[qubit[i]]である。では、qubit[i]⇒0.33とするには？
今、y軸での回転ゲートryで0度から180度まで１度おきに回転させた時の一つのqubitの状態１と状態０のもっている割合、つまり重ね合わせの混合割合の平方根（＝確率振幅）と、それを二乗した値＝混合割合＝確率を見てみよう。（グラフ）
参考コード
from blueqat import Circuit
import numpy as np
for j in range(180):
c=Circuit(1)
theta=j*np.pi/180
c.ry(theta)[0]
result=c.run()
q0=round(result[0].real,2)
q1=round(result[1].real,2)
q2=round(result[0].real**2,2)
q3=round(result[1].real**2,2)
print(j,'{0:>6} {1:>6} {2:>6} {3:>6}'.format(q0,q1,q2,q3))
セルに見立てたひとつの量子ビットにセットするのは、車や玉がそのセルにどの程度存在しているのか、つまり「１（車あるいは玉などが有）」の確率である。（100%そこに車が存在していれば１、全くいなければ0%だが、33%そこに車がいるとは？）グラフを見てわかる通りryの確率振幅は単なるsin関数（実際には$\sinθ/2$）である。その確率（確率振幅の二乗）は $\ P_1=sin^2θ/2$ である。従って、ある確率Pが与えられた場合、この式から角度θが逆算できる。数式が出て申し訳ないが、$\ θ=2arcsin ( \sqrt{P_1}\ ) $、このθをryゲートの回転角として量子ビットに適用すれば、その量子ビットに所望の確率（確率振幅）をセットすることができる。確率0.5、つまりアダマールの場合、θ＝π/2、つまり90度となる。尚、「１」はπ（180度）、「０」はゼロ度である。
qcamモジュールのpropinit関数に引数として、セル数N、Blueqatのサーキットオブジェクトc、セルに見立てた量子ビットの各々の確率の配列を渡すと、ryゲートを用いて確率をセットできる。
def propinit(NN, c, ppinit_a):
for i in range(NN):
theta=2* np.arcsin(np.sqrt(ppinit_a[i])) # Correlation between Probability and Rotation Angle
c.ry(theta)[i] # You can change from ry to rx.
return
また、y軸での回転ゲートに代えてｘ軸での回転ゲートrxも全く同様に使用できる。（ryをrx書き換えるだけでOK) rxに置き換えた場合、「０」はryの場合と全く同じだが、「１」には実数での確率振幅はなくなり、虚数の確率振幅を持つようになる。だが、二乗すれば虚数（複素数）でも確率となり同じである。それは兎に角「１」の確率 $P_1$は、「０」の余事象であるので、$P_1＝1-P_0$。また、$P_1$ は $P_1=cos^2θ/2$ であるため、$P_1＝1-cos^2θ/2=sin^2θ/2$ となる。故に、これはryゲートの場合と全く同様となる。
u3ゲートで確率＋αを与えてみた
一つの量子ビットに対して、重ね合わせ的には「１」とその余事象である「０」しかなく、確率としての「１」、「０」を捉えているが、確率振幅であれば、もう一つ情報を追加できると考えた。「１」の実数と虚数、「０」の実数部である。「０」の虚数部は常にゼロである。確率以外の情報を増やすと何が生じるのか知りたくて考えたのだが... 当然、セルに見立てた量子ビットには単なる確率ではなくて、複素数をセットしなくてはならない。「１」の確率を $P_A$、 $P_B$の二つにわけ、前者を実数、後者を虚数に振り分けることとして、残りは余事象である、$P_0$とする。今、$P_0=0.3$とした場合、$P_A+P_B=0.7$。さらに$P_A=0.5$、$P_B=0.2$ とする場合、initial_aの配列要素としてnumpyのcomplex型を用いて、0.5+0.2j と記せばよい。これを量子ビットにセットするには、ユニバーサルゲートであるu3を使えばできる。u3のθ、φ、λをryやrxの時と同様に考えて以下の式を得る。 * $θ=2arccos\sqrt P_B$ * $φ＝arccos\sqrt{P_A/(1-P_0)}$ * $λ＝arcsin\sqrt{P_B/(1-P_0)}$ 以下はqcam.impropinitのコードで、q-cam184iのみに使用される。def rotangle(probi):
Areal=probi.real
Bimag=probi.imag
Zero=1-Areal-Bimag
if Zero==1:
phi=0
lam=0
else:
phi=np.arccos(np.sqrt(Areal/(1-Zero)))
lam=np.arcsin(np.sqrt(Bimag/(1-Zero)))
theta=2*np.arccos(np.sqrt(Zero))
return theta,phi,lam
def impropinit(NN, c, ppinit_a):
for i in range(NN):
ptheta,pphi,plam=rotangle(ppinit_a[i])
c.u3(ptheta,pphi,plam)[i]
return
#2-3. Blueqatの状態ベクトル表示からのベクトルと確率の抽出
Blueqatには、量子コンピューティングの結果をみる手段として、観測と状態ベクトルの「覗き見」が用意されている。観測に関しては後述。状態ベクトルののぞき見命令はc.run()である。状態ベクトルを見る方法ではあるが、ベクトルそのものがこれで分かるわけではなく、その状態ベクトルの確率振幅が状態ベクトルの数だけ、複素数表記の一次元リストの形で戻ってくる。状態ベクトルとその数って？ 例えば、全部で４量子ビットなら、状態ベクトルは、００００、１０００、０１００、１１００、００１０、１０１０、０１１０、１１１０、０００１、１００１、０１０１、１１０１、００１１、１０１１、０１１１、１１１１の１６個である。状態ベクトルの数は $2^4=16$ 、つまり、全部でn量子ビットならば $2^n$ 個となる。q-camの場合の全量子ビット数は大概20を超えるので、$2^{20}=1,048,576$ つまり百万個以上である。
Blueqatで確率振幅が状態ベクトルの数だけ複素数表記の一次元リストの形で戻ってくるとは、要は百万個以上の複素数の要素を持ったリストが答えとなっているということである。だが、幸いq-cmaの場合、その大多数の状態ベクトルは確率振幅ゼロ、つまり意味なしである。だから、確率振幅なり、その二乗である確率（その状態ベクトルが出現する確率）がゼロより大きいものだけを抽出すればよい。それは、しらみつぶしに c.run() の戻り値リストの先頭から順に値がゼロより大のものを抽出するということである。しかし、複素数の値そのものは確率振幅であって、それがどの状態ベクトルの値なのかが分からないのでは。・・・んなわけはないはず。
上述した４量子ビット時の状態ベクトルを左が小さい桁で表した二進数表記と考えると左から順に大きくなっていく。１０進数に直せば、０，１，２，３，４，５，６，７，８，９，１０，１１，１２，１３，１４，１５である。Blueqatの状態ベクトルの覗き見の結果は、まさにこの順で出力されている。例えば、全４量子ビットで、０番目と２番目の量子ビットにアダマールゲートを適用した時のc.run() の結果を見てみよう。
from blueqat import Circuit
c=Circuit(4)
c.h[0,2]
print(c.run())
上のコードをコンソールで実行した結果が以下である。
このc.run() の結果に、前から順番に１０進数を付けた後、それを二進数に変換し、さらに、その二進数を逆順にして、４桁に満たない場合は右に４桁になるように０をつける。（下表参照）これでできた二進数こそ、複素数のリストに対応した状態ベクトルである。
リストの先頭から１づつ増加するカウンタでループさせた場合、このカウンタこそが、確率振幅に対応した状態ベクトルである。前述したが、q-camでは状態ベクトルが百万以上になるため、確率振幅なり、その二乗である確率（その状態ベクトルが出現する確率）がゼロより大きいものだけを抽出し、その時のカウンタのみから状態ベクトルを求めれば良く、そのことを抽出と称している。その抽出関数がqcamモジュールのqcam.qvextractである。pythonで10進数のカウンタをbinで２進数化した時、先頭に二進数であることを明示する'0b'が付いてくる。この'0b'付き二進数をリスト化してdelで'0b'をカットする。さらに、qcamの場合、確率として必要なのは全量子ビット（NR個)の内、ゼロレジスタ（初めのN個）の部分のみであるため、全量子ビットのベクトルから一番末尾に来るN個分をスライスする。
def qvextract(NN,RN,RW,extract_a):
reg0_a=np.array([['0']*NN]*2**NN) # [result-number,cell-number]
tempo_list=['0']*NN
eprobv_a=np.array([0]*2**NN,dtype='complex')
ereprov_a=np.array([0]*2**NN,dtype='float')
nm=0
for i in range(len(extract_a)):
if abs(extract_a[i])>0.0001: # Decimal 'i' is the very vector!
tempo1_list=list(bin(i)) # vector extraction but reversed. (Decimal to Binary)
del tempo1_list[0:2] # Delete '0','b' in tempo list
tempo1_list[0:0]=['0']*(NN*RN-len(tempo1_list)) # Zero-fill to adjust digits
reg0_a[nm,0:NN]=tempo1_list[(1-RW)*NN-1:-RW*NN-1:-1] # slice and reverse
eprobv_a[nm]=extract_a[i] # Probability amplitude
nm+=1
for i in range(nm):
ereprov_a[i]=eprobv_a[i].real*eprobv_a[i].real+eprobv_a[i].imag*eprobv_a[i].imag
return nm,reg0_a,ereprov_a
#2-4. Blueqatのベクトルと観測
2-3.までで c.run()の結果から状態ベクトルとそれに対応したそのベクトルが出現する確率までが抽出できた。だが、知りたいのは量子ビットごとの「１（車あるいは玉などが有）」の確率である。2-3.での例では、０番目と２番目の量子ビットにアダマールゲートを適用した。ということは、2-2.で説明した「０と１がちょうど半々で重ねあった状態、つまり確率0.5」を０番目と２番目の量子ビットにセットしたことを意味している。つまり、qubit[0]⇒0.5 qubit[2]⇒0.5となっている筈である。
2-3.の結果の中から確率振幅>0のものだけを抽出すると、１０進数で言う０，１，４，５である。これらの状態ベクトルは下表に示した通りである。今、確率①（確率振幅の二乗）とベクトルの各要素②との積を③とする。各ベクトルごとのこの③の値を、量子ビット番号ごとに和をとると、表に示した通り、０番目と２番目の量子ビットが0.5となる。つまり、qubit[0]⇒0.5 qubit[2]⇒0.5であり、これでやっと c.run()の結果から知りたいことを知ることができた。
q-camで、この処理を行う関数が以下のqcam.qcalcdである。二次元配列 csum_a に、j番目のベクトルのi番目の量子ビットの値③を格納する。こうすることによりループを回さずにnumpyのsum関数np.sum(csum_a,axis=0)を用いて、量子ビットごとの確率を求め、それを一次元配列cproreal_aに格納している。こうすることにより、c.run()と言うBlueqatの「覗き見」命令から量子ビットごとの「１（車あるいは玉などが有）」の正確な確率を得ることができる。
def qcalcd(NN, cnum, cvect_a, cprobare_a):
cproreal_a=np.array([0]*NN,dtype = 'float')
csum_a=np.array([[0]*NN]*cnum,dtype = 'float')
for j in range(cnum):
for i in range(NN):
csum_a[j,i]=float(cvect_a[j,i])*cprobare_a[j]
cproreal_a=np.sum(csum_a,axis=0)
cstd=np.std(cproreal_a)
return cstd,cproreal_a
だが、リアルの量子コンピュータでは「覗き見」などできないはずである。代わりに「観測」を行う。この「観測」を行ってしまうと、残念なことに、ある量子ビットが１と０を同時にある割合で重ね合わされていたとしても、その割合を直接的に知ることはできず、「観測」をしたとたんに、結果は１か０のどちらかになってしまう。この量子力学上の残念な性質を波束の収縮と言う。でも心配することはない。同じ試行（この場合は量子コンピューティング）を何度も繰り返すと、１状態が出る回数と０状態が出る回数の割合を知ることができる。例えば、１００回繰り返して、１状態が６０回、０状態が４０回という結果であれば、大体、状態１の確率が0.6であると知ることができる。
Blueqatではこのような観測と試行の繰り替えしとして、 c.m[:].run(shots=number) （numberは試行回数）が準備されている。下のプログラムprobshot.pyではセル数を８とし、試行回数を1000回としている。0~N-1量子ビットの初期セット確率は [0.4,0,1,0.5,1,0,0,0.2] とした。これらは自由に書き換えて実行することができる。量子ビットにこの確率セットした後、何の量子コンピューティングもせずに観測のみ実行する。何の量子コンピューティングもせずに観測しているので、当然セットした確率と同じ確率が結果となるはずであるが・・・
尚、c.m[:].run(shots=number)の結果はdict（辞書）型で、キーは状態ベクトル、値はその状態ベクトルが出現した回数である。従って、この出現回数を試行回数で割れば確率となる。
コード
from blueqat import Circuit
import numpy as np
import qcam
# -------- Setting ------------------------------------------------------------------
N=8 # number of qubits :You can change.
nrun=1000 # number of runs :You can change.
initial_a=np.array([0.4,0,1,0.5,1,0,0,0.2],dtype='float') # initial probability distribution :You can change.
vector_a=np.array([[0]*N]*(2^N),dtype='int')
csum_a=np.array([[0]*N]*(2^N),dtype='float')
final_a=np.array([0]*N,dtype='float')
prob_a=np.array([0]*(2^N),dtype='float')
#---------- Main Body ----------------------------------------------------------------
ret='y'
while ret=='y':
c=Circuit(N)
qcam.propinit(N,c,initial_a)
ans=c.m[:].run(shots=nrun)
num=len(ans)
kk=list(ans.keys())
vv=list(ans.values())
for i in range(num):
vector_a[i,0:N]=list(kk[i]) # Result Vectors
prob_a[i]=vv[i]/nrun
stdev,final_a=qcam.qcalcd(N,num,vector_a, prob_a)
qcam.qcresultout(N, 0, num, initial_a, vector_a, prob_a,stdev, final_a)
ret=input(' NEXT(Y/N)?')
以下、コンソールでの実行結果である。NEXT(Y/N)?でｙを入力すると、q-camのような時間発展ではなく、ただ単に同じことを繰り返す。
この出力結果の見方を簡単に説明する。TimeとStepはq-camの用の時間発展を示しているのでここでは無視。sumは各量子ビットの確率の単純な和。またstdevはその標準偏差だが、やはりここでは関係ない。Initial Cell-Probability は初期にセットした確率。Result 1 からは確率が0より大きい結果となった状態ベクトルとそのベクトルの確率（＝出現回数／試行回数の％）をProbabilityとして出力している。最後の Final Cell-Probability は状態ベクトルとその確率から計算した、各量子ビットごとの確率である。
さて、当然、何も処理をしていないのだから Initial Cell-Probability と Final Cell-Probability は等しくなければならないが、どうであろう。１０００回の試行を三度繰り返しているが、小数点第２位では結構な誤差が出てしまっている上、毎回少しずつ異なる。観測とは、あるいは波束の収縮とは、このようなものだと理解するしかない。（この観測もリアル量子コンピュータではなく、それらしくシミュレートしたものであることは言うまでもない）
q-cmaでは正確な解が欲しい。となると、c.run()での覗き見を使うしかない。量子コンピュータと謳いながら・・・。シミュレーターさまさまである。念のために確認してみよう。
コード
from blueqat import Circuit
import numpy as np
import qcam
# -------- Setting ------------------------------------------------------------------
N=8 # number of cells :You can change.
R=1 # number of registries
initial_a=np.array([0.4,0,1,0.5,1,0,0,0.2],dtype='float') # initial probability distribution :You can change.
vector_a=np.array([[0]*N]*(2^N),dtype='int')
csum_a=np.array([[0]*N]*(2^N),dtype='float')
final_a=np.array([0]*N,dtype='float')
#---------- Main Body ----------------------------------------------------------------
ret='y'
while ret=='y':
c=Circuit(N)
qcam.propinit(N,c,initial_a)
master_a=np.array(c.run())
num,vector_a,prob_a=qcam.qvextract(N,R,1,master_a)
stdev,final_a=qcam.qcalcd(N, num,vector_a,prob_a)
qcam.qcresultout(N, 0, num, initial_a, vector_a, prob_a,stdev, final_a)
ret=input(' NEXT(Y/N)?')