$$
\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}}
$$
はじめに
自作の量子計算シミュレータqlazyで、量子アルゴリズムの基本を確認するよ。今回は「位相推定アルゴリズム」です。参考にさせていただいた記事を以下に挙げておきます。
位相推定アルゴリズム
問題設定 (やりたいこと)
ユニタリ演算子$U$とそれに対する一つの固有状態$\ket{\psi}$が与えられたとき、
U \ket{\psi} = e^{i 2\pi \phi} \ket{\psi}
が成り立ちます(後の議論の都合により、$2\pi$を単位とした位相を$\phi$と定義しておきます)。ここで、$e^{i 2\pi \phi}$はこの固有状態に対応した$U$の固有値です。「位相推定アルゴリズム」というのは、この$\phi$の値を推定するアルゴリズムのことです。
前々回の記事で「アダマールテスト」を確認したときに、測定を大量に繰り返すことで、同様の位相推定ができるという説明をしました。その際、収束精度について説明しませんでしたが、「位相推定アルゴリズム」は、それよりも効率の良いアルゴリズムとして提案されたもののようです。
解決方針 (どんな考え方で実行するか)
$\phi$を以下のように2進小数表現します。
\phi = 0.\phi_{1} \phi_{2} \cdots \phi_{n} \qquad (\phi_{i} = \{0,1\})
与えられた$U$と$\ket{\psi}$を使って量子計算をガチャガチャやりながら、$\ket{\phi_{1} \phi_{2} \cdots \phi_{n}}$に対応した確率振幅のみ大きくなる状態が生成できたら、位相推定ができたといえそうです。なぜなら、この状態を観測することで、{0,1}を要素とする系列($\phi_{1}, \phi_{2}, \cdots , \phi_{n}$)が高い確率でわかり、これがすなわち位相を2進小数表現した際の各桁の値に相当するので。
ということで、$e^{i 2\pi \phi}$の$\phi = 0.\phi_{1} \phi_{2} \cdots \phi_{n}$に含まれる数列{$\phi_{i}$}を状態$\ket{\phi_{1} \phi_{2} \cdots \phi_{n}}$として、どう埋め込むのかを考えます。具体的にどうしたら良いでしょうか。ヒントは「量子フーリエ変換」です。量子フーリエ変換の定義は、
\ket{j} \rightarrow \frac{1}{\sqrt{2^{n}}} \sum_{k=0}^{2^{n}-1} exp(i \frac{2 \pi kj}{2^{n}}) \ket{k}
でした。この$\ket{j}$を$\ket{\phi}=\ket{\phi_{1} \phi_{2} \cdots \phi_{n}}$と思うことにすると、
\frac{j}{2^{n}} = 0.j_{1} j_{2} \cdots j_{n} = \phi
となるので、上の量子フーリエ変換の定義式は、
\ket{\phi} \rightarrow \frac{1}{\sqrt{2^{n}}}\sum_{k=0}^{2^{n}-1} exp(i 2 \pi k\phi) \ket{k}
と書き直せます。位相$\phi$を持った固有状態と関係づけるため、$\ket{\psi}$とのテンソル積の形にします。
\ket{\phi} \otimes \ket{\psi} \rightarrow \frac{1}{\sqrt{2^{n}}} \sum_{k=0}^{2^{n}-1} exp(i 2 \pi k \phi) \ket{k} \otimes \ket{\psi} \tag{1}
この右辺の状態をなんとか作り出すことができれば、逆量子フーリエ変換をやることで、$\ket{\phi}$をつくることができそうです。これが、位相推定アルゴリズムの基本的な考え方です。
アルゴリズムの導出
さて、それでは、実際にガチャガチャ量子計算をやりながら、この(1)式の右辺をこうしたら作り出せますよ、ということを見ていきます。まず、n量子ビットの状態$\ket{0}^{\otimes n} = \ket{00\cdots0}$を補助量子ビットとして用意します(以後、簡単のため単に$\ket{0}$と書きます)。で、各々の量子ビットにアダマールゲートをかけて、各基底状態が等確率で重ね合わさった状態を作ります。
\ket{0} \rightarrow \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} \ket{k}
これに$\ket{\psi}$に対する量子ビットを加えます。つまり、テンソル積をとります。
\ket{0} \otimes \ket{\psi} \rightarrow \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} \ket{k} \otimes \ket{\psi} \tag{2}
さらに、これを以下のように展開しておきます。
\begin{align}
&\frac{1}{\sqrt{2^{n}}} \sum_{k=0}^{2^{n-1}} \ket{k} \otimes \ket{\psi} \\
&= \frac{1}{\sqrt{2^{n}}} \sum_{k_{0}=0,1}\sum_{k_{1}=0,1} \cdots \sum_{k_{n-1}=0,1} \ket{k_{n-1} \cdots k_{0}} \otimes \ket{\psi} \\
&= \frac{1}{\sqrt{2^{n}}} \sum_{k_{0}=0,1}\ket{k_{0}} \sum_{k_{1}=0,1}\ket{k_{1}} \cdots \sum_{k_{n-1}=0,1}\ket{k_{n-1}} \otimes \ket{\psi} \tag{3}
\end{align}
量子回路で書くと、
|0> ----H---
...
|0> ----H---
|0> ----H---
|0> ----H---
|psi> -/------
です(ここで、|psi>に対する線は1本しか書かれていませんが、"-/--"で複数本あることを表していると見てください)。
(3)式の最後の式から何とか$exp(i 2\pi k \phi)$をひねり出したいので、ユニタリ演算子$U$を使って、以下に示す(4)式のような演算子を構成して状態全体にかけるようにします。ここで、$\Lambda_{i}(U^{2^i})$は、$\ket{k_i}$に対する量子ビット(上の量子回路で言うと、n個の補助量子ビットの中の一番下を0番目、一番上をn-1番目と見たときのi番目)を制御ビットとしたユニタリ演算子を状態$\ket{\psi}$に適用する演算子と定義しています。
\Lambda_{0}(U^{2^0}) \cdot \Lambda_{1}(U^{2^1}) \cdots \Lambda_{n-1}(U^{2^{n-1}}) \tag{4}
これを量子状態全体に適用します。
\begin{align}
&\frac{1}{\sqrt{2^n}} \sum_{k_0=0,1}\ket{k_0} \sum_{k_1=0,1}\ket{k_1} \cdots \sum_{k_{n-1}=0,1}\ket{k_{n-1}} \otimes \ket{\psi} \\
&\rightarrow \frac{1}{\sqrt{2^n}} (\Lambda_{0}(U^{2^0}) \sum_{k_0=0,1}\ket{k_0}) (\Lambda_{1} (U^{2^1})\sum_{k_1=0,1}\ket{k_1}) \cdots (\Lambda_{n-1} (U^{2^{n-1}})\sum_{k_{n-1}=0,1}\ket{k_{n-1}}) \\
&= \frac{1}{\sqrt{2^n}} \bigr[ (\ket{0}\bra{0} \otimes I+\ket{1}\bra{1} \otimes U^{2^0}) \sum_{k_0=0,1}\ket{k_0} \cdot (\ket{0}\bra{0} \otimes I+\ket{1}\bra{1} \otimes U^{2^1}) \sum_{k_1=0,1}\ket{k_1} \cdots (\ket{0}\bra{0} \otimes I+\ket{1}\bra{1} \otimes U^{2^{n-1}}) \sum_{k_{n-1}=0,1}\ket{k_{n-1}} \bigl] \ket{\psi} \\
&= \frac{1}{\sqrt{2^n}} (\ket{0}\bra{0} \otimes I + \ket{1}\bra{1} \otimes U^{2^0})(\ket{0}\bra{0} \otimes I + \ket{1}\bra{1} \otimes U^{2^1}) \cdots (\ket{0}\bra{0} \otimes I + \ket{1}\bra{1} \otimes U^{2^{n-1}})\ket{\psi} \\
&= \frac{1}{\sqrt{2^n}} \sum_{k_0=0,1} \sum_{k_1=0,1} \cdots \sum_{k_{n-1}=0,1} \ket{k_{n-1} \cdots k_1 k_0} \otimes U^{2^{k_0}} U^{2^{k_1}} \cdots U^{2^{k_{n-1}}} \ket{\psi} \\
&= \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} \ket{k} \otimes exp(i 2\pi \sum_{l=0}^{n-1} 2^l k_l \phi) \ket{\psi} \\
&= \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} \ket{k} \otimes exp(i 2\pi k\phi) \ket{\psi} \\
&= \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} exp(i 2\pi k\phi) \ket{k} \otimes \ket{\psi}
\end{align}
というわけで、(1)式が導けました。全体を量子回路で書くと、
|0> ----H----------------------- ... ---*----------
... ...
|0> ----H------------------*---- ... --------------
|0> ----H-----------*----------- ... --------------
|0> ----H----*------------------ ... --------------
|psi> -/------U(1)---U(2)---U(4)-- ... --U(2^(n-1))--
となります(ここで、U(i)は$U^{i}$を表すものとします)。あとは先述した(1)式の逆、つまり逆量子フーリエ変換(IQFT)を実行します。
\frac{1}{\sqrt{2^{n}}} \sum_{k=0}^{2^{n}-1} exp(i 2 \pi k\phi) \ket{k} \otimes \ket{\psi} \rightarrow \ket{\phi} \otimes \ket{\psi}
これで、欲しかった$\ket{\phi}$を得ることができました。最終的な量子回路は、補助量子ビットに対するIQFTを追加して、以下のように書けます。
|0> ----H----------------------- ... ---*-----------|IQFT|---
... ... ... ...
|0> ----H------------------*---- ... ---------------|IQFT|--- |phi>
|0> ----H-----------*----------- ... ---------------|IQFT|---
|0> ----H----*------------------ ... ---------------|IQFT|---
|psi> -/------U(1)---U(2)---U(4)-- ... --U(2^(n-1))------------ |psi>
シミュレータで確認
それでは、qlazyで、位相推定アルゴリズムの動作を確認します。今回、対象としたユニタリ演算子ですが、2量子ビットに対する2つの位相ゲートのテンソル積$P_0(\pi/2) \otimes P_1(\pi/4)$としました。行列表現にすると、
\begin{pmatrix}
1 & 0 \\
0 & e^{i \pi/2}
\end{pmatrix}
\otimes
\begin{pmatrix}
1 & 0 \\
0 & e^{i \pi/4}
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & e^{i\frac{\pi}{4}} & 0 & 0 \\
0 & 0 & e^{i\frac{\pi}{2}} & 0 \\
0 & 0 & 0 & e^{i\frac{3\pi}{4}}
\end{pmatrix}
です(Quantum Native Dojoに例として挙げられているものと同じです)。すでに対角化されているので、固有値と固有状態は一目瞭然ですが、念のため書き出すと、
| 固有値 | 固有状態 |
|---|---|
| 1 | $\ket{00}$ |
| $e^{i\frac{\pi}{4}}$ | $\ket{01}$ |
| $e^{i\frac{\pi}{2}}$ | $\ket{10}$ |
| $e^{i\frac{3\pi}{4}}$ | $\ket{11}$ |
という4パターンになります。位相推定をやらなくてもすでにわかっちゃっているのですが、ここは一旦わかっていない振りをして、アルゴリズムが正しく動作するかを確認します。
さて、位相推定アルゴリズムのスタート地点は、「固有状態はわかっているが固有値(の位相)がわからない」ということでした。今回、どれかの固有状態を適当に決めてアルゴリズムを走らせてみて推定結果を確認、としても良かったのですが、ちょっとひとひねり加えてみました。固有状態がわかっていないとして、ユニタリ演算子から推定することから始めてみました。前々回の記事で、アダマールテストを繰り返すことで確率的にどれかの固有状態に収束するということが確認できたので、同じことを今回の位相推定の前段階でやってみたというわけです。
では、以下にPythonコードの全体を示します。
【2021.9.5追記】qlazy最新版でのソースコードはここに置いてあります。
from qlazypy import QState
def iqft(self,id=None):
dim = len(id)
iid = id[::-1]
for i in range(dim):
phase = -1.0/2**i
for j in range(i):
self.cp(iid[j],iid[i],phase=phase)
phase *= 2.0
self.h(iid[i])
return self
def c_unitary(self,id=None): # id[0] = controll qubit id
self.cp(id[0],id[1],phase=0.5).cp(id[0],id[2],phase=0.25)
return self
def eigen_state_ctr(psi_dim):
# initial |psi> and controll qubit
qs = QState(psi_dim+1)
for i in range(psi_dim):
qs.h(i)
# hadamard test iteratively
id = [psi_dim] + list(range(psi_dim))
iter = 100
for _ in range(iter):
qs.h(psi_dim).c_unitary(id).h(psi_dim).m(id=[2],shots=1)
print("== eigen state (convergence) ==")
id = list(range(psi_dim))
qs.show(id)
return qs
def main():
# add methods
QState.iqft = iqft
QState.c_unitary = c_unitary
# prepare initial state
sub_dim = 3 # order of precision
psi_dim = 2 # qubit number for |psi> (must be qubit number for the unitary)
qs_sub = QState(sub_dim)
qs_psi_ctr = eigen_state_ctr(psi_dim) # estimate the eigen state of the unitary
qs_total = qs_sub.tenspro(qs_psi_ctr)
# qubit id list (for the following steps)
psi_qid = list(range(sub_dim,sub_dim+psi_dim))
sub_qid = list(range(sub_dim))
# hadamard operation to sub-qubits
for i in range(sub_dim):
qs_total.h(i)
# c-unitary operation to the eigen state:|psi> iteratively
for i in range(sub_dim):
times = 2**i
for _ in range(times):
id = [i] + psi_qid # qubit id list for eigen state:|psi>
qs_total.c_unitary(id)
# inverse QFT
qs_total.iqft(sub_qid)
# show final quantum state
print("== final quantum state ==")
qs_total.show(sub_qid)
# mesurement
print("== result of measurement ==")
md = qs_total.m(sub_qid)
md.show()
# phase estimation
print("== phase estimation ==")
print(2*md.lst/md.state_num,"* PI")
# free memory
qs_sub.free()
qs_psi_ctr.free()
qs_total.free()
if __name__ == '__main__':
main()
何をやっているか、簡単に説明します。
def iqft(self,id=None):
dim = len(id)
iid = id[::-1]
for i in range(dim):
phase = -1.0/2**i
for j in range(i):
self.cp(iid[j],iid[i],phase=phase)
phase *= 2.0
self.h(iid[i])
return self
で、逆量子フーリエ変換の関数を定義しています。前回の記事で「量子フーリエ変換」をやってみたましたが、そのとき掲載していた量子フーリエ変換の逆バージョンです。
def c_unitary(self,id=None): # id[0] = controll qubit id
self.cp(id[0],id[1],phase=0.5).cp(id[0],id[2],phase=0.25)
return self
で、今回対象としているユニタリ演算を制御ユニタリ化した関数を定義しています。上に示した量子回路を見ればわかる通り、ユニタリではなくそれを制御ユニタリにしたものしか登場してこないので、この関数があれば十分です。中身は制御位相ゲートを2つつなげているだけです。こうすれば位相ゲートを2つつなげたものの制御ユニタリ化が実現できますね。関数の引数idは量子ビット番号のリストで、id[0]を制御量子ビット番号と想定しています。
def eigen_state_ctr(psi_dim):
# initial |psi> and controll qubit
qs = QState(psi_dim+1)
for i in range(psi_dim):
qs.h(i)
# hadamard test iteratively
id = [psi_dim] + list(range(psi_dim))
iter = 100
for _ in range(iter):
qs.h(psi_dim).c_unitary(id).h(psi_dim).m(id=[2],shots=1)
print("== eigen state (convergence) ==")
id = list(range(psi_dim))
qs.show(id)
return qs
で、固有状態を推定してアウトプットする関数を定義しています。引数psi_dimは$\ket{\psi}$に対する量子ビット数(今回は2です)を表しています。が、アダマールテストで推定するため、制御量子ビットのためのビットを1個追加しています。なので、アウトプットする固有状態の量子ビットは本来2で良いのですが、3ビットの量子状態をアウトプットしています(最後の量子ビットは不要なので、この後の処理では無視します)。アダマールテストの実体は、
qs = QState(psi_dim+1)
for i in range(psi_dim):
qs.h(i)
で、初期状態を用意し(全ビットにアダマールをかけて各基底が等確率に重ね合わさるようにしました)、
iter = 100
for _ in range(iter):
qs.h(psi_dim).c_unitary(id).h(psi_dim).m(id=[2],shots=1)
で繰り返します。繰り返し回数は100回に決め打ちしました。これで4つある固有状態のどれかに(確率的に)収束するようになります。
これらを定義した上で、mainの処理です。
# add methods
QState.iqft = iqft
QState.c_unitary = c_unitary
で、逆量子変換と制御ユニタリゲートをQStateクラスのメソッドに加えます(安易にやっちゃっていますが本当は要注意、既存メソッドと衝突しないように)。
# prepare initial state
sub_dim = 3 # order of precision
psi_dim = 2 # qubit number for |psi> (must be qubit number for the unitary)
qs_sub = QState(sub_dim)
qs_psi_ctr = eigen_state_ctr(psi_dim) # estimate the eigen state of the unitary
qs_total = qs_sub.tenspro(qs_psi_ctr)
で、全体の初期状態を用意します。なるべく汎用的に使えるように、精度を表す補助量子ビット数をsub_dim変数として可変設定できるようにしました。psi_dimは対象としているユニタリ演算子の量子ビット数(=固有状態$\ket{\psi}$の量子ビット数)です。c_unitary関数の中でユニタリ演算をどう設定したかに応じて、この数を正しく入れます(今回は2です)。あとは補助量子ビットの状態を$\ket{0}$にして、固有状態$\ket{\psi}$を推定して、それらをテンソル積で合成します(qs_total)。
# qubit id list (for the following steps)
psi_qid = list(range(sub_dim,sub_dim+psi_dim))
sub_qid = list(range(sub_dim))
で、この後の処理で必要になる$\ket{\psi}$と補助量子ビットに対する、量子ビット番号のリストをつくっています。
あとは、位相推定アルゴリズムに対応した量子回路の実行です。
# hadamard operation to sub-qubits
for i in range(sub_dim):
qs_total.h(i)
で、補助量子ビットに対するアダマール。
# c-unitary operation to the eigen state:|psi> iteratively
for i in range(sub_dim):
times = 2**i
for _ in range(times):
id = [i] + psi_qid # qubit id list for eigen state:|psi>
qs_total.c_unitary(id)
で、制御ユニタリゲートを連続実行。
# inverse QFT
qs_total.iqft(sub_qid)
で、最後の逆量子フーリエ変換。
# show final quantum state
print("== final quantum state ==")
qs_total.show(sub_qid)
で、一応、測定前の量子状態を表示しています。
位相推定結果は、これを測定することで得られますので、測定します。
# mesurement
print("== result of measurement ==")
md = qs_total.m(sub_qid)
md.show()
測定結果は変数mdに格納されます。
# phase estimation
print("== phase estimation ==")
print(2*md.lst/md.state_num,"* PI")
mdのプロパティlstに測定結果(3ビット2進数の10進整数)、state_numに量子状態数(今回は3ビットなので8)が入っているので、それを使って最後に推定された位相を表示しています。
以下に、実行結果を示します。
== eigen state (convergence) ==
c[00] = +0.0000+0.0000*i : 0.0000 |
c[01] = +0.0000+1.0000*i : 1.0000 |+++++++++++
c[10] = -0.0000+0.0000*i : 0.0000 |
c[11] = -0.0000-0.0000*i : 0.0000 |
== final quantum state ==
c[000] = +0.0000+0.0000*i : 0.0000 |
c[001] = -0.0000-1.0000*i : 1.0000 |+++++++++++
c[010] = +0.0000-0.0000*i : 0.0000 |
c[011] = -0.0000-0.0000*i : 0.0000 |
c[100] = +0.0000-0.0000*i : 0.0000 |
c[101] = +0.0000-0.0000*i : 0.0000 |
c[110] = +0.0000+0.0000*i : 0.0000 |
c[111] = +0.0000+0.0000*i : 0.0000 |
== result of measurement ==
direction of measurement: z-axis
frq[001] = 100
last state => 001
== phase estimation ==
0.25 * PI
ということで、この場合、固有状態$\ket{01}$に収束し、測定結果は"001"でした。これに対応した位相は2進小数で言うと$2\pi \times 0.001$です。10進小数で書くと$\pi/4=0.25\pi$となるので、正しく位相推定がなされたことになります。めでたしめでたし。
おわりに
今回、アダマールテストで固有状態を推定してから、位相推定をしてみましたが、単純に勉強したことをつなげてみたかったからです。おそらくアダマールテストを使った収束は、効率がとてもよろしくないです。誤差のないシミュレータで、しかも今回のようなとても小さい演算子で、10O回繰り返してようやくきれいに収束したレベルです。誤差のある実機では、ちょっと現実的ではない気がします。なので、NISQ前提の量子アルゴリズムでは、量子化学計算で位相推定を使う場合も、固有状態を仮設定して古典コンピュータとのハイブリッドで実現するのが現実的とされているのだと思います(このあたりまだ勉強中なので、間違ったことを言っていたらご指摘ください)。
以上