古典ビットのシミュレートでは満足できない
以前、128古典ビットの加算器をBlueqatで実装しました。
けれど、まぁ、古典だしそりゃできるよね、つまんないよね。やっぱり重ね合わせたいよね。
そこで、今回はGottesman-Knillの定理を使ったバックエンドを作ってみました。
その場合、Toffoliゲートは残念ながら使えなくなってしまうのですが、X, Y, Z, H, S, S†, CNOTゲートが使えます。
重ね合わせ状態が作れるのは嬉しいですね！
なお、任意回転や$T$ゲートや$T^\dagger$ゲートがなく、特定の量子状態しか作れないので、このシミュレータではできない量子計算があります。Toffoliゲートも使えなくなった都合で、古典でできたものも書けなくなっている場合があります。
Gottesman-Knillを理解するための基礎知識
群論の知識がいくつか必要になります。
スタビライザー
$|\Psi\rangle$に対する群$G$の作用を考えます。$g \in G$について$g|\Psi\rangle = |\Psi\rangle$となるとき、$|\Psi\rangle$が$g$によってスタビライズされている、といいます。$|\Psi\rangle$をスタビライズするような群の元を集めてきて$Stab(G) = \{g \in G \ ;\ g|\Psi\rangle = |\Psi\rangle\}$を作ることができ、これも群をなします。この群をスタビライザーとよびます。
通常の量子計算と同じように、$|\Psi\rangle$はベクトル、作用する群の要素はユニタリ行列と考えましょう。$|\Psi\rangle$が0でない場合は、スタビライザーは可換群となり、また、スタビライザーに$-I$は含まれません。
さらに、今回は、考えるユニタリ行列は、単位行列、パウリ行列(とそのクロネッカー積)、パウリ行列(とそのクロネッカー積)のマイナス1倍だけに限定することにします。
P_n = \underbrace{\{I, \pm X, \pm Y, \pm Z \}\otimes \{I, \pm X, \pm Y, \pm Z \} \otimes \cdots}_{n\ times}
また、以下では
P_2 = \{II, IX, -IX, \cdots, XI, -XI, -XY, \cdots \}
のように、行列要素を続けて書くことでクロネッカー積を表す書き方や、IZを$Z_2$, XYを$X_1 Y_2$と表すような書き方をする場合があります。
生成元
群Gの元(ここでは行列)はいっぱいあっても、いくつかの元を持ってきて、掛け算をすれば、群全体を作ることができる場合があります。このような元の集まりを群の生成系Sと呼び、SはGを生成する、といいます。Sの元を生成元と呼びます。
例えば、$|00\rangle$のスタビライザーは、全部書くと$\{II, IZ, ZI, ZZ\}$ですが、IZとZIだけを用意すれば、掛け算で他のものも作ることができます。つまり、$\{II, IZ, ZI, ZZ\}$は$\{IZ, ZI\}$から生成されます。
ここで注意が必要なのは、生成元の取り方は一意ではない、ということです。
例えば、$\{II, IZ, ZI, ZZ\}$自体も$\{II, IZ, ZI, ZZ\}$の生成系です。
数が2個のものに限ったとしても$\{ZI, ZZ\}$や$\{IZ, ZZ\}$も生成系になります。
なお、今回考えるパウリ行列の群では、$n$量子ビットの状態のスタビライザーは、最小で$n$個の生成元から生成することができます。
量子計算では初期状態を$|0\cdots 0\rangle$と置くことが多いですが、そのスタビライザーの生成元を考えると、取り方の一例として$\{Z_1, Z_2, \cdots Z_n\}$があります。
Gottesman-Knillの定理の基本的な考え方
通常の量子コンピュータのシミュレーションでは、量子状態をメモリに保存し、ゲートを適用するたびに量子状態を表した配列を更新します。
Gottesman-Knillの定理を使ったシミュレータでは、量子状態ではなく、スタビライザーの生成系(つまり、生成元の配列)をメモリに保存し、ゲートを適用するたびに、スタビライザーの生成系を更新します。
例えば、$|\Psi\rangle$にゲートUをかけて、$U|\Psi\rangle$になった場合、$U|\Psi\rangle$のスタビライザーとその生成系は何でしょうか。$g$をスタビライザーの生成元のひとつとすると、
\begin{eqnarray}
U|\Psi\rangle &=& Ug|\Psi\rangle\\
&=& UgU^{\dagger}U|\Psi\rangle\\
&=& (UgU^{\dagger})U|\Psi\rangle
\end{eqnarray}
という式変形が成り立ちます。(1つ目: $g$はスタビライザーの元なので$|\Psi\rangle = g|\Psi\rangle$, 2つ目: $U^{\dagger}U = I$なので途中に入れてもよい, 3つ目: カッコを付けただけ)
この式をよく見ると、$U|\Psi\rangle$が$UgU^{\dagger}$によってスタビライズされていることが分かります。
この関係を利用して、ゲート$U$を適用するたびに、生成元の配列について、$g$を$UgU^{\dagger}$に更新します。
ゲートと、パウリ群$g$の関係については、予め計算したものがNielsen & Chuangなどに載っています。計算するだけなので、自分でしてもいいのですが、ネットで見つけた資料から引用します。
(出典: https://www.cs.umd.edu/~amchilds/teaching/w10/project-sample.pdf )
パウリなのに$Y$がないじゃないか！と思う方もいるかもしれませんが、$Y = iXZ$なので、$Y = iUXU^{\dagger}UZU^{\dagger}$が成り立ち、$X$のものと$Z$のものを掛けて$i$倍すればいいです。
どうやって測定するのか
スタビライザの生成元を追いかけるのはいいとして、どうやって測定するんだ、と疑問に思うかもしれません。
まず、測定もパウリ行列$M$で表します。通常の測定はZ基底での測定なので、$i$番目の測定は$M = Z_i$で表せます。
続いて、$M$がすべての生成元と可換であるかを確かめます。
可換でないものが含まれていた場合
可換でないものが複数ある場合は、可換でないものをひとつだけにします。
そのためにまず、可換でないものの中から、ひとつ適当に選びます。
他の可換でないものを、先ほど選んだものと掛け算した値に更新します。
こうすると、更新された値は測定のパウリ行列と可換になります。
可換でないものがひとつだけになったら、コインを投げて$|0\rangle$か$|1\rangle$かを決めます(確率1/2ずつ)
- $|0\rangle$に決まった場合は、可換でないもののスタビライザを、$M$に置き換えます
- $|1\rangle$に決まった場合は、可換でないもののスタビライザを、$-M$に置き換えます
これで、可換でないものの場合はおしまいです。
すべての生成元と可換だった場合
このとき、$M$または$-M$のどちらかは、スタビライザに含まれています。ですが、生成元の取り方はいろいろあるので、生成元の配列の中に含まれているかどうかは分かりません。
もし生成元の配列の中に$M$か$-M$が含まれていれば、$M$なら$|0\rangle$、$-M$なら$|1\rangle$を測定したことになります。
また、この場合は測定しても状態は変わりません。
含まれない場合、どうすればいいか具体的に書いてある資料が見つからなかったのですが、今回、生成元の配列要素同士を掛け算しながら、Mか-Mを作り出しました。
(生成元を選んで掛け算すれば、必ず、ちょうどぴったりあう組み合わせがあるはずです)
私は、掃き出し法のようなイメージで、ひとつずつ順にIに合わせていく形で作りました。
実際にやってみた
$ pip install blueqat_gk_backend
$ python -m blueqat_gk_backend
GHZ State
Circuit().h[0].cx[0, 1].cx[0, 2].m[:].run(shots=100)
numpy
Counter({'000': 54, '111': 46})
Gottesman-Knill
Counter({'000': 56, '100': 44})
2 qubit grover (Expected: |11>)
Circuit().h[:].h[1].cx[0, 1].h[1].h[:].x[:].h[1].cx[0, 1].h[1].x[:].h[:].m[:].run(shots=100)
numpy
Counter({'11': 100})
Gottesman-Knill
Counter({'11': 100})
Circuit (4 qubits)
Circuit().h[:].cx[0, 1].cx[1, 2].cx[2, 3].y[:].cx[1, 0].cx[2, 1].cx[3, 2].h[:].m[:]
numpy
Counter({'1010': 100})
Gottesman-Knill
Counter({'1010': 100})
Circuit (10 qubits)
Circuit().h[:].cx[0, 1].cx[1, 2].cx[2, 3].cx[3, 4].cx[4, 5].cx[5, 6].cx[6, 7].cx[7, 8].cx[8, 9].y[:].cx[1, 0].cx[2, 1].cx[3, 2].cx[4, 3].cx[5, 4].cx[6, 5].cx[7, 6].cx[8, 7].cx[9, 8].h[:].m[:]
numpy
Counter({'1010101010': 100})
Gottesman-Knill
Counter({'1010101010': 100})
Circuit (100 qubits)
Circuit().h[:].cx[0, 1].cx[1, 2].cx[2, 3].cx[3, 4].cx[4, 5].cx[5, 6].cx[6, 7].cx[7, 8].cx[8, 9].cx[9, 10].cx[10, 11].cx[11, 12].cx[12, 13].cx[13, 14].cx[14, 15].cx[15, 16].cx[16, 17].cx[17, 18].cx[18, 19].cx[19, 20].cx[20, 21].cx[21, 22].cx[22, 23].cx[23, 24].cx[24, 25].cx[25, 26].cx[26, 27].cx[27, 28].cx[28, 29].cx[29, 30].cx[30, 31].cx[31, 32].cx[32, 33].cx[33, 34].cx[34, 35].cx[35, 36].cx[36, 37].cx[37, 38].cx[38, 39].cx[39, 40].cx[40, 41].cx[41, 42].cx[42, 43].cx[43, 44].cx[44, 45].cx[45, 46].cx[46, 47].cx[47, 48].cx[48, 49].cx[49, 50].cx[50, 51].cx[51, 52].cx[52, 53].cx[53, 54].cx[54, 55].cx[55, 56].cx[56, 57].cx[57, 58].cx[58, 59].cx[59, 60].cx[60, 61].cx[61, 62].cx[62, 63].cx[63, 64].cx[64, 65].cx[65, 66].cx[66, 67].cx[67, 68].cx[68, 69].cx[69, 70].cx[70, 71].cx[71, 72].cx[72, 73].cx[73, 74].cx[74, 75].cx[75, 76].cx[76, 77].cx[77, 78].cx[78, 79].cx[79, 80].cx[80, 81].cx[81, 82].cx[82, 83].cx[83, 84].cx[84, 85].cx[85, 86].cx[86, 87].cx[87, 88].cx[88, 89].cx[89, 90].cx[90, 91].cx[91, 92].cx[92, 93].cx[93, 94].cx[94, 95].cx[95, 96].cx[96, 97].cx[97, 98].cx[98, 99].y[:].cx[1, 0].cx[2, 1].cx[3, 2].cx[4, 3].cx[5, 4].cx[6, 5].cx[7, 6].cx[8, 7].cx[9, 8].cx[10, 9].cx[11, 10].cx[12, 11].cx[13, 12].cx[14, 13].cx[15, 14].cx[16, 15].cx[17, 16].cx[18, 17].cx[19, 18].cx[20, 19].cx[21, 20].cx[22, 21].cx[23, 22].cx[24, 23].cx[25, 24].cx[26, 25].cx[27, 26].cx[28, 27].cx[29, 28].cx[30, 29].cx[31, 30].cx[32, 31].cx[33, 32].cx[34, 33].cx[35, 34].cx[36, 35].cx[37, 36].cx[38, 37].cx[39, 38].cx[40, 39].cx[41, 40].cx[42, 41].cx[43, 42].cx[44, 43].cx[45, 44].cx[46, 45].cx[47, 46].cx[48, 47].cx[49, 48].cx[50, 49].cx[51, 50].cx[52, 51].cx[53, 52].cx[54, 53].cx[55, 54].cx[56, 55].cx[57, 56].cx[58, 57].cx[59, 58].cx[60, 59].cx[61, 60].cx[62, 61].cx[63, 62].cx[64, 63].cx[65, 64].cx[66, 65].cx[67, 66].cx[68, 67].cx[69, 68].cx[70, 69].cx[71, 70].cx[72, 71].cx[73, 72].cx[74, 73].cx[75, 74].cx[76, 75].cx[77, 76].cx[78, 77].cx[79, 78].cx[80, 79].cx[81, 80].cx[82, 81].cx[83, 82].cx[84, 83].cx[85, 84].cx[86, 85].cx[87, 86].cx[88, 87].cx[89, 88].cx[90, 89].cx[91, 90].cx[92, 91].cx[93, 92].cx[94, 93].cx[95, 94].cx[96, 95].cx[97, 96].cx[98, 97].cx[99, 98].h[:].m[:]
Gottesman-Knill
Counter({'1010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010': 100})
効率を考えず、とりあえず動くことを優先したためか、100量子ビットだと少し待たされます。
ですが、100で動いているということは、指数関数的な計算時間からは開放されていることは間違いないですし、また、普通のシミュレータだと100量子ビットはそもそも動かないので、動くだけでも嬉しいかな、と思います。
コードは https://github.com/gyu-don/blueqat-gottesman-knill-backend に置いています。よければご覧ください。
まとめ
以前から気にはなっていたGottesman-Knillの定理を実装してみよう、という試みでした。
Nielsen & Chuang (日本語版だと3巻、まだ売ってるやつです)にもGottesman-Knillの定理は載っていて、今回記事を書くにあたり、大変参考になりましたが、最初にどういうスタビライザを用意するか、測定の具体的な手順といったものは書かれておらず、非常に苦労しました。
ところで、なんとなく、まだバグありそうなんだよなぁ……
もし何か見つけたら、バグの再現する回路をつけて、コメントかこちらにIssue登録お願いします。
通常のnumpyバックエンドと結果を比較するとバグかバグでないか分かりやすいです。