8-1. はじめに
QIITA初投稿だった１回目は1年以上前の巣籠GW。それから時は流れ、再びの巣籠GWを超え、五輪もお盆もパラも越え。このシリーズは「新視点セルオートマトンとグローバー」と題する・・・故に最終回の今回はセルオートマトンとグローバーに量子位相推定を適用したら周期は求まるのか？をblueqatで実装し調べてみたという話。ただ周期を推定するにあたっての工夫や気づきは本投稿には書ききれないので別投稿予定？ セルオートマトンやグローバーって周期性あるの？ セルオートマトンではクラス2の場合パターンが繰り返す周期性があるし、グローバーでは振幅がステップの繰り返しに対し綺麗なサイン波となっている。が・・・
以降、量子位相推定をQPE（Quantum Phase Estimation）、セルオートマトンをCA（Cellular Automaton）、グローバーアルゴリズムをGA（Grover Algorithm）と略す。
8-2. QPEのCAやGAへの適用
概要
図１にQPEの概念的な量子回路図を示した。この回路（QPE)を適用すると図中の下部分のユニタリー演算Uの固有値($e^{2\pi iθ}$)の位相θがカウンタC（nビット）の二進数出力$C_{out}(=2^nθ)$として現れる。カウンタとは、ユニタリー演算における演算回数のカウンタであって、GAならアルゴリズムの適用回数（あるいはステップ数）であり、CAでは時間発展の単位数（あるいはステップ数）である。
図1
図ではカウンタの各線からカウンタ側が制御、標的側が制御ユニタリーゲート（U）に結ばれているためカウンタが'1'ならUが実行される。流れとしてカウンタが二進数で示している回数に応じた分のユニタリー演算を実施することを意味している。例えば、カウンタ下から2番目が1であれば$2^1=2$回分ユニタリー演算をすることになる。当然下位（一つ下の線）ビットも1だった場合、見ての通りユニタリー演算を$2^0+2^1=3$で3回繰り返すことになる。CAやGAではそれぞれ三回時間発展なり三回増幅（減少）させたこととなる。従ってユニタリー演算の結果がユニタリー側のビットには反映されていく。一方カウンタ側のビットは制御側のビットなので'1','0'は不変であり何も生じない。が、実はエンタングルメントによりブーリアン(1,0)とは別に位相キックバックと言って位相は変化している。そのビット毎に変わった位相を抽出し二進数にエンコードするのが量子逆フーリエ変換（$FT^{-1}$）部分である。カウンタの最初にアダマールで重ね合わせておけば、自動的にカウンタで表現できる全二進数の組み合わせ、つまり十進数で言うところの $0,1,2,3,$・・・ $2^n-1$の全てを計算することになる。当然のことではあるがもし量子逆フーリエ変換（$FT^{-1}$）をしないでカウンタを出力させると単にこのカウンタの $0,1,2,3,$・・・ $2^n-1$全てが$1/2^n$という等確率で出現するだけのこととなってしまう。要は量子逆フーリエ変換（$FT^{-1}$）で位相キックバックで生じた各qubitの位相方向情報を抽出しているわけだ。するとユニタリー演算の単なる結果とは別のこの演算に関する情報（固有値）が得られる。詳しいことは、ググれば結構な数の解説がでてくる。
位相から周期へ
$C_{out}$から位相は $θ=C_{out}/2^n$（nはカウンタのビット数）で求まる。ここから周期Tへはmを自然数とすれば$e^{2\pi iθ}=e^{2\pi i(m/T)}$より θ=m/T で目出度く求まることになる。が、しかし・・・出てくる位相は、つまり$C_{out}$は、一つではない。実行するたびに異なることが殆どであろう。何回も実行すると、ある$C_{out}$が出現する確率が分かる。θ=0や1を除く出現確率が大きいθほど答えに近いことが多い。よって最大確率の位相で目出度く求まることになる。が、しかし・・・ θ=m/T でmは1,2,3,4...・・・ショアのアルゴリズムの位数を導出するやり方としてよく解説されているが、位相の小数値から約分できない分数を求め、その分母を周期Tとして推定する方法である。例えば、θ＝0.666の確率が大きい場合、この分数表記は2/3、従ってθ=m/Tなのでm=2,T=3と求めるのである。幸いpythonのFractionモジュールのfractionsを用いれば任意の小数を約分できない分数に直してくれる。尚、答の精度はカウンタCのビット数nが大きければ大きいほど良いとされている。目出度し・・・とはいかない！
8-3. 実装
制御ユニタリーゲート
制御ユニタリーゲートを実装するにはどうすればよいのか？ GAやCAを制御ユニタリーゲート化しなければならない。図2に示した2ビットのGAを例で考えよう。
図2
制御ユニタリーゲートと言うとあたかも一つのゲートのような気がするがそうではない。何度ユニタリー変換してもユニタリー変換である故。
一つのカウンタのビットを制御ビットとして、図2のGAで使われている全てのゲート(図2ではZ、H、Xゲート）を制御ゲートに、つまり、Z→CZ、H→CH、X→CX とする。図3にその様子を示した。ただ全てのゲートと言っても制御ビットの状態に無関係、つまり制御ビットが'1'であろうが'0'であろうが結果に影響しないゲートは制御化する必要はない。
図3
これで実装可能に・・・が、しかし少なくともblueqatではCXやCCX（トフォリゲート）は準備されているが、CHやCZは準備されていない。ではどうする？
アダマールゲートとＺゲートの制御化
アダマールゲートもＺゲートもあるいは他の位相や回転ゲートでもユニバーサルゲートでも、（マルチ）制御NOTゲートを用いて（マルチ）制御化可能である。マルチ制御NOTゲートに関しては3回目の簡単に作成できるマルチ制御ゲート を参照されたし。Zゲートの場合は簡単である。パウリ演算変換Z=HXHを使えばよいので、（マルチ）制御NOTの前後の標的ビットにのみアダマールをすればよい。一方、アダマールの制御化の場合、まず（マルチ）制御NOT、次に標的ビットのみ $-\pi/4$回転させ、次に再び（マルチ）制御NOT、次に$\pi/4$回転、最後に再び（マルチ）制御NOTをすれば良い。
GAの制御ユニタリー
前述を基に図3の➀から➉までをblueqatのコードで示すと以下の通りとなる。ただし、量子ビットの割り当てとして、第1qubit,第2qubitをGA用としga0,ga1で表し、またカウンタの一つのビットを第0qubitとしてbfで表している。
def qproc(bf):
c.h[ga1].ccx[bf,ga0,ga1].h[ga1] #➀
c.cx[bf,ga0].ry(-np.pi/4)[ga0].cx[bf,ga0].ry(np.pi/4)[ga0].cx[bf,ga0] #➁
c.cx[bf,ga1].ry(-np.pi/4)[ga1].cx[bf,ga1].ry(np.pi/4)[ga1].cx[bf,ga1] #➂
c.cx[bf,ga0].cx[bf,ga1] #➃，➄
c.h[ga1].ccx[bf,ga0,ga1].h[ga1] #➅
c.cx[bf,ga0].cx[bf,ga1] #➆，➇
c.cx[bf,ga1].ry(-np.pi/4)[ga1].cx[bf,ga1].ry(np.pi/4)[ga1].cx[bf,ga1] #➈
c.cx[bf,ga0].ry(-np.pi/4)[ga0].cx[bf,ga0].ry(np.pi/4)[ga0].cx[bf,ga0] #➉
return
尚、➁、➂、➈、➉は制御アダマールゲートだが回転ゲートryが制御cryになっていない。制御化する必要がないためである。 4bit GA
3bitと4bitのGAの制御ユニタリーのコード
3bit GA
def qproc(bf):
c.h[ga2]
qmcn.cccx(c,bf,ga0,ga1,ga2)
c.h[ga2]
c.cx[bf,ga0].ry(-np.pi/4)[ga0].cx[bf,ga0].ry(np.pi/4)[ga0].cx[bf,ga0]
c.cx[bf,ga1].ry(-np.pi/4)[ga1].cx[bf,ga1].ry(np.pi/4)[ga1].cx[bf,ga1]
c.cx[bf,ga2].ry(-np.pi/4)[ga2].cx[bf,ga2].ry(np.pi/4)[ga2].cx[bf,ga2]
c.cx[bf,ga0].cx[bf,ga1].cx[bf,ga2]
c.h[ga2]
qmcn.cccx(c,bf,ga0,ga1,ga2)
c.h[ga2]
c.cx[bf,ga0].cx[bf,ga1].cx[bf,ga2]
c.cx[bf,ga2].ry(-np.pi/4)[ga2].cx[bf,ga2].ry(np.pi/4)[ga2].cx[bf,ga2]
c.cx[bf,ga1].ry(-np.pi/4)[ga1].cx[bf,ga1].ry(np.pi/4)[ga1].cx[bf,ga1]
c.cx[bf,ga0].ry(-np.pi/4)[ga0].cx[bf,ga0].ry(np.pi/4)[ga0].cx[bf,ga0]
return
def qproc(bf):
c.h[ga3]
q.ccccx(c,bf,ga0,ga1,ga2,ga3)
c.h[ga3]
c.cx[bf,ga0].ry(-np.pi/4)[ga0].cx[bf,ga0].ry(np.pi/4)[ga0].cx[bf,ga0]
c.cx[bf,ga1].ry(-np.pi/4)[ga1].cx[bf,ga1].ry(np.pi/4)[ga1].cx[bf,ga1]
c.cx[bf,ga2].ry(-np.pi/4)[ga2].cx[bf,ga2].ry(np.pi/4)[ga2].cx[bf,ga2]
c.cx[bf,ga3].ry(-np.pi/4)[ga3].cx[bf,ga3].ry(np.pi/4)[ga3].cx[bf,ga3]
c.cx[bf,ga0].cx[bf,ga1].cx[bf,ga2].cx[bf,ga3]
c.h[ga3]
q.ccccx(c,bf,ga0,ga1,ga2,ga3)
c.h[ga3]
c.cx[bf,ga0].cx[bf,ga1].cx[bf,ga2].cx[bf,ga3]
c.cx[bf,ga3].ry(-np.pi/4)[ga3].cx[bf,ga3].ry(np.pi/4)[ga3].cx[bf,ga3]
c.cx[bf,ga2].ry(-np.pi/4)[ga2].cx[bf,ga2].ry(np.pi/4)[ga2].cx[bf,ga2]
c.cx[bf,ga1].ry(-np.pi/4)[ga1].cx[bf,ga1].ry(np.pi/4)[ga1].cx[bf,ga1]
c.cx[bf,ga0].ry(-np.pi/4)[ga0].cx[bf,ga0].ry(np.pi/4)[ga0].cx[bf,ga0]
return
qmcnは3回目の簡単に作成できるマルチ制御ゲート を参照されたし。
CAの制御ユニタリー
CAの場合はクラス2でかつ補助レジスタを使用しないルールに限られる。とすると4回目セルオートマトンのアルゴリズムに記載したルール60かルール102あたりが最適である。結局、ここではルール102を採用しよう。ルール102はCXゲートのみの超簡単な演算故、制御ユニタリ化と言ってもcxに、制御をカウンタのビット(例えば上記同様bf)としてccx化すれば完了である。尚、n_counterはカウンタのビット数でCAのセルの先頭のqubit番号、またNは n_counter+セル数 で CAのセルの末端のqubit番号を意味している。
def qproc(bf):
for i in range(n_counter,N):
c.ccx[bf,1+i, i]
return
検証用の制御ユニタリー
QPEでどの程度周期が推定できるものか、周期が明確で単純な例で検証するためのユニタリーを考える。$2\pi$で1周期となる回転ゲートの回転角を変化させることが最も単純と思われるがこれは別の投稿に記す。ここでは、CAと類似したユニタリー演算として、m個のセルの端に車がいて、それが1ステップ毎に単に隣に進む演算を考える。車が端に達したら出だしに戻る。この場合、セル数mが周期Tとなる。図4(下の折りたたみ）参照。このアルゴリズムはステップ毎に隣とスワップすれば済む。以降、このモジュールをFRQXと称す。
図4 検証用ユニタリ時間発展の様子
def qproc(bf):
for i in range(n_counter,N-1):
c.ccx[bf,i, i+1].ccx[bf,i+1,i].ccx[bf,i,i+1]
return
return
量子逆フーリエ変換
QPEにとってもう一つの大事な要素が量子逆フーリエ変換$FT^{-1}$である。数理的な解説はググればたくさあるので省略。ただ、コードはシンプルで以下の通りである。尚、n_counterはカウンタCに用いるビット数nである。
def IQFT():
for i in range(n_counter):
for j in range(i):
c.cu1(-np.pi/(2**(i-j)))[i,j]
c.h[i]
return
プログラム
前述のqprocモジュールを含む実際に実行する2bitのGAにQPEを適用したプログラム全体を以下の折りたたみに示した。またGithubのリポジトリ名 q-cams の中に、上記プログラムq-frqga2.py及び、3bit~6bitまでのGAにQPEを適用したq-frqga3.py,q-frqga4.py,q-frqga5.py,q-frqga6.py、またルール102のCAにQPEを適用したq-frqca102.pyを公開している。qcamモジュールの使い方等は2回目の「2-1.q-camの構造とモジュール」に記載している。3つのプログラムともqprocと初期設定以外は同じで、Main Bodyでユニタリーのモジュール部分がカウンタに応じたべき乗回、回るようにしている。2bitのGA以外の場合、qproc(bf)に各ユニタリーのモジュールを2bitのGAの代わりにはめ込めば動作する。n_counterはカウンタのビット数、n_unitaryはユニタリー演算のビット数である。また同じ$C_{out}$であってもユニタリー側のベクトルが異なる場合、別個に出現確率を持つが、同じ$C_{out}$であると言うことは同じ位相となるので、xprobと言う配列を用いて$C_{out}$毎の確率にまとめている。
2bitのGAにQPEを適用したプログラムのコード
# q-frequency Applying QPE to Grover 2bit by Gigagulin 2021 Aug.
from blueqat import Circuit
import numpy as np
import qcam
from fractions import Fraction
# -------- Setting ------------------------------------------------------------------
n_counter,n_unitary=8,2
N=n_counter+n_unitary # number of qubit
initial_a=np.array([0.5]*n_counter+[0.5,0.5],dtype='float') # initial set
xprob=np.array([0]*2**n_counter,dtype='float') # Probability of counter
ga0,ga1=n_counter,n_counter+1
# -------- inversed QFT ------------------------------------------------------------------
def IQFT():
for i in range(n_counter):
for j in range(i):
c.cu1(-np.pi/(2**(i-j)))[i,j]
c.h[i]
return
# -------- q-process ------------------------------------------------------------------
def qproc(bf):
c.h[ga1].ccx[bf,ga0,ga1].h[ga1] #➀
c.cx[bf,ga0].ry(-np.pi/4)[ga0].cx[bf,ga0].ry(np.pi/4)[ga0].cx[bf,ga0] #➁
c.cx[bf,ga1].ry(-np.pi/4)[ga1].cx[bf,ga1].ry(np.pi/4)[ga1].cx[bf,ga1] #➂
c.cx[bf,ga0].cx[bf,ga1] #➃，➄
c.h[ga1].ccx[bf,ga0,ga1].h[ga1] #➅
c.cx[bf,ga0].cx[bf,ga1] #➆，➇
c.cx[bf,ga1].ry(-np.pi/4)[ga1].cx[bf,ga1].ry(np.pi/4)[ga1].cx[bf,ga1] #➈
c.cx[bf,ga0].ry(-np.pi/4)[ga0].cx[bf,ga0].ry(np.pi/4)[ga0].cx[bf,ga0] #⑩
return
# -------- Main Body --------------------------------------------------------------
c=Circuit(N)
qcam.propinit(N,c,initial_a)
for m in range(n_counter):
power=2**m
for k in range(power):
qproc(m)
IQFT()
master_a=np.array(c.run())
num,vector_a,prob_a=qcam.qvextract(N,1,1,master_a)
print(' counter phase probability guess')
for j in range(num):
xdeci=0
for i in range(n_counter):
xdeci=int(vector_a[j,i])*2**(N-n_unitary-1-i)+xdeci
xprob[xdeci]=prob_a[j]+xprob[xdeci]
for i in range(2**n_counter):
if xprob[i]>=0.00005:
phase=i/(2**n_counter)
guess=str(Fraction(phase).limit_denominator(100))
print('{:>5}'.format(i),' ','{:>5,.4f}'.format(phase),end=' ')
print(' ','{:>5,.5f}'.format(xprob[i]),end=' ')
spa=' '*(6-len(guess))
print(spa,':'+str(guess))
8-4. 検証
周期推定
検証用のユニタリーモジュールFRQXを前述のプログラムのq-processに組み込み、n_counter=10 (カウンタとして十進で0~1023)、n_unitary=7 (セル数7=周期T=7)とし、コマンドプロンプト上で実行すると図5（下の折りたたみ）の左に示した出力となる。元来、量子コンピュータであれば、1ショット実行した場合、ポッと一つだけ17桁の二進数が出力されるだけだろう（この17桁のうち上位10桁分が$C_{out}$に相当する）。そこで何度も何度もショットを繰り返す。すると同じ二進数が出現する確率が次第に分かってくる。これではあまりにも大変だ。幸いblueqatのシミュレーションとしてショットを繰り返さず初めからこの確率を求めることができる。しかけは、2回目の2-3. Blueqatの状態ベクトル表示からのベクトルと確率の抽出を参照されたし。前述のプログラムでそれを利用して図5にあるcounter（十進）とphase($=counter/2^{10}$)及びその位相の出現確率probablityを、さらにfractionモジュール(分母のリミットは100としている）でphaseの値の分数表記の推定をguessとして出力している。ご覧の通りcounterは0から順に全てを出力しておらず、これは画面が煩雑にならないように確率が0.00005未満のものを出力させていないからである。（そもそもの確率がゼロの場合も鼻から出力されてはいない）さて、周期を求めるには、まずこの出力の中から確率が最大のものあるいは確率が他より大きいものの位相を探すこととなるが、これは確率降順でソートすればよい。実は出力をそのままエクセルにコピペして確率降順で表としたものを図5の中央に示した。さらにまた図5の右には位相を横軸、確率を縦軸にした位相の確率分布を示した。ただし位相ゼロの確率が群を抜いて高いのでそこを除いてる。以降、位相をθ、確率をP(θ)、最大確率を$P_{max}$と略す。θの値からの推定分数をG分数、その分母の値をGとする。またQPEで予測した周期をTp、実際の周期をTaとしよう。
図5
図5の真ん中の表から$P_{max}$時のθはゼロとなっている。ゼロは周期推定にならない。次にP(θ)が高いものからは、同確率ものが二つづあり、確率の高いものから9番目までの片方のG分数の分母、つまりGは7となっている。原則的にはこの値が周期の推定値第一候補である。従ってこの例ではすんなりTp=7と推定できかつこれは実際の周期Taと一致している。尚、この検証用のユニタリーの場合、同確率がペアで出現するがそのペアは必ずθの区間 $2^{-(k+1)}$～$2^{-k}$ (k=0，1，2,…)内に現れ、区間の中央値のθに対して必ず対象の位置にきている。そのため片方が分かればもう一方も直ぐ手計算で分かる。詳細は別投稿予定。
実はセル数7は幸運な例かもしれない。表１にセル数2から11までの確率トップ3までのθとP(θ)とG分数を載せた。緑色で色付けしたG分数がTp=Taで正しく周期推定されている部分である。セル数つまり周期が大きくなるにつれて推定が難しくなっていくように見える。
表1
がしかし、P(0)に着目するとセル数が増大するにつれて小さくなっていることが分かる。セル数、つまり周期TaとP(0)をよく見ると
$$P(0)=1/Ta つまり Ta=Tp=1/P(0)$$となっている。従って$P_{max}$でもあるθ=0の確率P(0)が分かればこのユニタリーの系統では周期が求まる。
8-5.CAへのQPE適用の結果
CAルール102の周期
実はルール102でセルが無限の時は周期がない。周期があるのは周期的境界条件の時のみである。4回目セルオートマトンのアルゴリズム参照。図6（下の折りたたみ）に、セルの右端のみ'1'としてセル数4, 5, 7, 10 及びセル数は7だがセルの両端と中央を'1'とした場合（|u>違い）の時間発展の様子を示した。それぞれの周期を数えるとセル数4では Ta=4、セル数5と7と７の|u>違いではいずれも Ta=8、そしてセル数10では Ta=16である。
図6 ルール102の時間発展の様子
実は周期Taとセル数の関係は、セル数n_cellが $2^{i-1}$＜n_cell ≦$2^i$ (iはi>1自然数)の場合$T_a=2^i$となる。よってセル数7であればi=3の範囲内のため Ta=$2^3$=8。またセル数10であればi=4の範囲のため Ta=$2^4$=16 である。
CAの結果
CAのユニタリーモジュールをプログラムのq-processに組み込み、n_counter=10 (カウンタとして十進で0~1023)とし、図6に示した各パターンとなるようにn_unitaryにセル数をセット、さらにプログラムの#initial set の配列に|u>をセットした。例えば、セル数7で|u>違いはinitial_a=np.array([0.5]*n_counter+[1,0,0,1,0,0,1],dtype='float')のように。そして、セル数 4, 5, 7, 7の|u>違い 及び10でQPEを実行。
その結果・・・コントロールパネルの出力は煩雑なので、それぞれの位相の確率分布を図7（下の折りたたみ）の赤点プロットのグラフに示した。また、検証用のユニタリーでセル数4と8及び16での位相の確率分布を青点プロットのグラフで示している。
図7 ルール102 位相の確率分布
CAの場合、検証用のユニタリーでの位相の確率分布と周期が同じ場合、全く同じ確率分布となっていることが分かる。また各セル数でのCAのP(0)とそこからの推定周期は以下の通りである。従って CA(ルール102) にQPEを適用することで周期は推定できると言えよう。
| CA | P(0) | 推定周期Tp(=1/P(0)) | 実際の周期Ta |
|---|---|---|---|
| セル4 | 0.25 | 4 | 4 |
| セル5 | 0.125 | 8 | 8 |
| セル7 | 0.125 | 8 | 8 |
| セル7｜u>違い | 0.125 | 8 | 8 |
| セル10 | 0.0625 | 16 | 16 |
8-6.GAへのQPE適用の果て
GAの周期再考
GAはアルゴリズムの適用回数つまり何回繰り返すかで、マーキングしているベクトル、例えば3bitで |111>をマーキングした場合、繰り返し回数に対して|111>の出現確率が上下する。2bit, 3bit, 4bit, 5bit 及び 6bit のGAで初期状態を均一に重ね合わせた場合の繰り返し回数(以降t）に対する、|11>, |111>, |1111>, |11111> ,|111111> の出現確率の振動の様子を図8（下の折りたたみ）に示した。青点がtに対する出現確率であり、グレーの線はtを連続した変数とみなしたサイン関数である。具体的な関数に関しては7回目グローバーアルゴリズム一考を参照されたし。図8の左のグラフはtが20回まで、中央のグラフは単に200回まで、右のグラフはサイン関数を除いた1000回にまで拡張したグラフである。ちなみに世間のGA解説によると、nビットGAの周期は$T=2^{(n/2)-1}\pi$ で与えられるとしている。実際、図8のグレーの線のサイン関数の周期は概ねその通りとなっている。が、しかし・・・
図8
ところでWeblio辞書によると周期とは「定期的に同じことが繰り返される事象において、任意のある時点の状態に一度循環して戻るまでの期間（時間）または段数のことである」とされている。任意の時点の青点が一度循環して戻るのまでの段数は、果たしてグレーの線のサイン関数の周期と同じであろうか？図8の左のグラフの青点だけを見た時？2bitGAは確かに3回毎に青点は同じ確率に戻っているしグレーの線も同じ周期である。が、他のビットでは？ 20回までで循環して元に戻っているようには見えない。だとするとグレーの線のサイン関数での周期は本当に周期と言えるのであろうか？
図8の中央や右のグラフをやや遠くから眺めるとグレーの線とは別の循環が見えてくる。3bitで右のグラフの1000回までを見ると超長周期のサイン波が見て取れる。その半周期はほぼ696であろうか。一方、中央、つまり200回までのグラフでは周期があるのかよくわからない。4bitの場合はどうか？図8の中央、つまり200回までのグラフにおいて青点のピークとピークの周期が31と見え、また大波のサイン波の周期は174と読める。さらに5bitでは同様にピークトゥピークで53、大波で496と読める。6bitは少し様子が異なり分かりずらいが中央のグラフにおいて、青点のピークトゥピークは25である。一方、右のグラフでの大循環は確認できない。
一口に周期と言っても、いくつかの周期が考えられ、果たしてQPEを適用し得られる位相情報からの周期推定はどの周期なのだろうか？以下の表に各GAのbit数での周期をまとめてみた。尚、理屈上の周期を$T_{thy}$、グレーの線のサイン周期を$T_{short}$、青点のピーク間周期を$T_{mid}$、青点の大波の周期を$T_{long}$とする。
| GA bit数 | $T_{thy}$ | $T_{short}$ | $T_{mid}$ | $T_{long}$ |
|---|---|---|---|---|
| 2bit | 3.14 | 3.00 | ---- | ---- |
| 3bit | 4.44 | 4.35 | 113 | 1391 |
| 4bit | 6.28 | 6.22 | 31 | 174 |
| 5bit | 8.89 | 8.84 | 53 | 496 |
| 6bit | 12.57 | 12.50 | 25 | ---- |
GAの結果
さて、q-frqgr2.py, q-frqgr3.py,q-frqgr4.py,q-frqgr5.pyの実行結果は？カウンタビット数n_counterを14とし実行。コマンドプロンプト状態の結果ではわかりずらいので、図9（下の折りたたみ）にそれぞれの位相の確率分布のグラフを赤点プロットで示した。また、図9の2bitGAの横に、検証用のユニタリーでセル数3（Ta=3)の結果を示した。ご覧の通り2bitGAとセル数3の検証用のユニタリーのグラフはθ=0を除くと確率分布はほぼ同一となっている。つまり2bitGAの周期はTp=3と推定される。一方、他のbitの確率分布では検証用のユニタリーで得られている分布と類似のものは全くない。そこで、従前通りP(θ)の相対的に大きいもののθから作るG分数からの周期推定を考えることになる。
図9 GAの位相の確率分布
表2に確率トップ10(降順）のθとそのG分数を載せた。尚、図9から視覚的に分かるが、検証用ユニタリーの場合と同様、θの区間 $2^{-(k+1)}$～$2^{-k}$ (k=0，1，2,…)内に同確率θがペアで出現し、区間の中央値のθに対して必ず対象の位置にきている。トップ10と言ったがペアなので確率は5位まで、以降、上から$P_{max}$、$P_2$…$P_5$と記す。
表2の2行目に書かれている100だの2000だのの数値はG分数を計算するFractionモジュールで指定した分母のリミット数である。100を超える周期を算定するのに分母が100では不都合なためである・・・本来、周期を知らない筈であるから周期の都合で分母リミットを自在に変えることはできないのだが。
また表2でG分数を塗りつぶしている色に関しては、青系は$T_{short}$、紫系は$T_{mid}$、緑系は$T_{long}$の推定元であることを表し、色が濃いものはGそのものずばりがTpであること現している。逆にGそのものがTpではないとは以下の二点の推定方法の使用を意味する。
➀ $Tp=G±1$ 誤差ありとしての推定・・・QPEはカウンタビット数が少ないと誤差大故（条件次第では14bitといえども）
➁ $Tp=G/2^j$ （jは負を含む整数）として推定・・・詳細は別投稿予定・・・？？
以上を念頭に結果を見ていこう・・・
表2
3bitGAは残念な結果？ 先に示した3bitGAでの各種周期とずばり一致するGはない。まずは$P_{max}$と$P_3$の薄い青で示したG=73。これは推定方法➁でj=4と考えた場合、$T_p=73/2^4=4.56$ であり$T_{short}$ に近い（苦しい?）。そして薄い緑の1393。$T_{long}$ の1391に近い。が2異なる。もし事前に3bitGAの各種周期を一切知らず本当に推定したならば、$P_5$のものを$T_p$とするか疑義があろう。
4bitGAでは塗りつぶしが多い。$P_{max}$の薄い青のG=99は推定方法➁でj=4とすると$T_p=99/2^4=6.19$ であり$T_{short}$に近い。$P_2$と$P_3$にあるG=87も➁でk=-1とすると$T_p=174$であり$T_{long}$となる。また$P_4$の緑ではG=174で$T_{long}$、紫ではG=31で$T_{mid}$そのものとなっている。しかし同じθからの分数の作り方で両方の周期を現しているのは偶然の一致であろう。
5bitGAも塗りつぶしが多い。$P_{max}$から$P_3$にかけての紫でのG=53は、まさに$T_{mid}$であり信頼性が高い。一方G=31と62は➁推定法から$T_{long}$の496となるがjが大きい上にG=53のカウンターペアとして出ており怪しい。
6bitGAは3bitGAとともにあまり良い結果ではない。$P_{max}$において➀からG=50、さらに➁からj=1であれば$T_{mid}$、j=2であれば$T_{short}$と推定される。$T_{mid}=2T_{short}$故当然ではある。
結局、周期推定は可能であったと言えるのだろうか。完全とはほど遠いが部分的には可能と言ったところか。
不均一な重ね合わせのGAだったら
GAで初期状態|u>が均一重ね合わせでない場合どうなるであろうか。2bitGA、3bitGA、4bitGAの初期状態|u>をそれぞれ、[0.3,0.7]、[0.3,0.5,0.7]、[0.4,0.3,0.6,0.5]として、またn_counterを14とし実行。結果は図10の通りである。図9のグラフと較べれば一目瞭然だが、たった一点P(0)に値があることを除くと全く同一形状である。グラフの下に上位10の確率のθとそのG分数を示したが、これも表2とみ較べれば一目瞭然、θ=0を除くと値も順位も全く同一である。つまり、初期値の重ね合わせが均一であろうが不均一であろうが、周期は全く変わらないことを示している。実際、GAの周期はbit数のみで決まり、初期状態|u>には依存しない。ただし、振幅と位相は変わる。詳しくは7回目グローバーアルゴリズム一考を参照されたし。
逆に言えば、QPEのGAの適用の結果は周期にのみ依存していると言えよう。
8-7.おわりに
以上、結論は以下の二点である。
・CA(ルール102) にQPEを適用することで周期は推定できる
・QPEのGAの適用の結果は周期にのみ依存している
「新視点セルオートマトンとグローバー」と題するシリーズはこれで終わりとする。が別の量子コンピューティングのネタのはじまり。とりわけ投稿者自身が疑問に思っていることを投稿していきたい。で、最近、多くのショアアルゴリズムの解説や説明で核心部が不十分というか不満というか、そのことを疑問に思っており、そのあたりから忘れたころに投稿していこうと思う。その前にQPEでの周期推定に関しての補足を投稿するつもりではいる。忘れたころに。