量子情報理論の基本：マジック状態蒸留

$$
\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}}
$$

はじめに

誤り訂正符号を使ってユニバーサル量子計算を行うためには、「(1)論理パウリ演算」「(2)論理CNOT演算」「(3)論理アダマール演算」「(4)論理位相シフト演算」の各々を誤り耐性がある形で構成する必要があります。このうち、(1)(2)(3)は1量子ビットに対する（単体の）クリフォード演算(X,Y,Z,H,CNOTゲート)だけで実現できます¹。が、(4)については単体のクリフォード演算(X,Y,Z,H,CNOTゲート)以外に、十分に精度の高い特別な状態「マジック状態(magic state)」をアンシラとして用意しないと実現できません²。今回の記事では、このマジック状態を作成する手法である「マジック状態蒸留(magic state distillation)」について勉強します。「蒸留」という名前の通り、精度の良くない状態を沢山用意しておいて、それを入力にしてフィルターで濾し取るような作業を何回も繰り返すことで精度の高い状態を取得する手法です。どのような理屈でこんなことが実現できるのか理解できたところで、量子計算シミュレータqlazyを使ってそのプロセスをシミュレーションしてみます。

参考にさせていただいたのは、以下の文献です。

1. Austin G. Fowler, Matteo Mariantoni, John M. Martinis, Andrew N. Cleland, "Surface codes: Towards practical large-scale quantum computation", arXiv:1208.0928 [quant-ph]
2. K.Fujii,"Quantum Computation with Topological Codes - from qubit to topological fault-tolerance",arXiv:1504.01444v1 [quant-ph] 7 Apr 2015
3. イエルン・ユステセン、トム・ホーホルト「誤り訂正符号入門（第2版）」森北出版（2019年4月）

理論説明

マジック状態の定義

まず復習です。作成したいマジック状態は、

\ket{Y} = \ket{0} + e^{i\pi/2} \ket{1}  \tag{1}

と

\ket{A} = \ket{0} + e^{i\pi/4} \ket{1}  \tag{2}

です。$\ket{Y}$が用意できればパウリゲートとアダマールゲートとCNOTゲートだけで$S$ゲートを実現することができます。また、$\ket{A}$が用意できれば同じくパウリゲートとアダマールゲートとCNOTゲートだけで$T$ゲートを実現することができます。論理$S$ゲートと論理$T$ゲートを実現したければ$\ket{Y}, \space \ket{A}$の論理状態版である、

\ket{\bar{Y}} = \ket{\bar{0}} + e^{i\pi/2} \ket{\bar{1}}  \tag{3}

\ket{\bar{A}} = \ket{\bar{0}} + e^{i\pi/4} \ket{\bar{1}}  \tag{4}

を用意します。ここで記号の上のバーは論理状態を意味していると思ってください（以下同様）。表面符号（Braiding）の上で式(3)や式(4)のような論理的なマジック状態を作りたい場合は、「表面符号によるユニバーサル量子計算（３）〜Braidingで論理位相シフト演算〜」で説明したように、どこかの量子ビットを$X$測定し最小の欠陥対を作っておいてから$X$測定した量子ビットに単体の$S$ゲート（または単体の$T$ゲート）を適用し、欲しい符号距離になるまで欠陥対を拡大します。ここで利用する$S$ゲート（または$T$ゲート）の精度が十分に高ければ出来上がる$\ket{\bar{Y}}$（または$\ket{\bar{A}}$）の精度は高くなります（もちろん欠陥拡大の途中で誤り訂正をきちんとやる前提です）。$\ket{\bar{Y}}$（または$\ket{\bar{A}}$）の精度が高くなれば、それを使って実現される論理$S$ゲート（または論理$T$ゲート）の精度は高くなります。すなわち、精度の良い論理シフト演算が実現できるようになります。ということで、以下では大本の$S$ゲート（または$T$ゲート）を精度良く実現するために必要となる単体の$\ket{Y}, \space \ket{A}$をいかにして蒸留して取得するかについて説明していきます³。

マジック状態蒸留の量子回路：|Y>

$\ket{Y}$を蒸留する回路を示します。

ここで、右下の$\ket{\psi}$が蒸留されて出てくる$\ket{Y}$を表しており、青枠で示した$7$個の$S$ゲートは$\ket{Y}$と$H$ゲートとCNOTゲートを使って以下のように構成したものです。

この回路からわかる通り、$7$個の$\ket{Y}$を入力して$1$個の$\ket{Y}$を出力する格好になっています。後の節で説明しますが、この回路によって精度の悪い入力$\ket{Y}$から少しだけ精度の高い$\ket{Y}$を得ることができます。これを繰り返し繰り返し適用することで十分に精度の高い$\ket{Y}$を得るという手法がマジック状態蒸留です。一回あたり$7$個の入力から1個の出力しか得られないので、最初に大量の$\ket{Y}$を用意しておく必要があります。非常に単純化して言うと$7^2=49$個の$\ket{Y}$を用意すると上の回路を2回繰り返すことができます。$7^3=343$個の$\ket{Y}$を用意すると3回繰り返すことができます。$7^N個$の$\ket{Y}$を用意すると$N$回繰り返すことができます⁴。ここで、十分な精度を得るために無限に繰り返さないといけないんじゃね？と心配される方がいるかもしれませんが、ご安心ください。繰り返し数$N$に対して指数的に精度が向上するという良い性質をもった手法になっています、というあたりも後で説明します。

完璧なマジック状態を入力した場合

さて、これでどうして精度が上がるのかという説明をしたいのですが、これ一見したところ何をやっているのか全く謎な回路ですよね。というわけで、入力の$\ket{Y}$が完璧な精度だった場合、右下から完璧な精度の$\ket{Y}$が出力される回路であるということを最初におさえておきたいと思います。まず、上の蒸留回路の一番下の$2$つの量子ビットの入力側を見てください。入力状態$\ket{0}\ket{+}$にCNOTがかかっているのでCNOT直後のこの$2$量子ビット状態はベル状態$\ket{0}\ket{0} + \ket{0}\ket{0}$になります⁵。その後、そのベル状態の上の方の量子ビット（上から数えて$7$番目の量子ビット）と$1$番目から$6$番目の量子ビットとがCNOTゲートで何やら複雑にエンタングルされています。実はこれ、$7$番目の量子ビットをSteane符号で符号化する回路になっています。つまり、もし$7$番目の量子ビットが$a \ket{0} + b \ket{1}$だった場合、複雑なCNOTの後の$7$量子状態は$a \ket{\bar{0}} + b \ket{\bar{1}}$になっています。ちなみに、このSteane符号に対応したスタビライザーは、

XIIXIXX \\
IXIIXXX \\
IIXXXXI \\
ZIIZIXX \\
IZIIZZZ \\
IIZZZZI \tag{5}

で、$3$つの$X$スタビライザーと$3$つの$Z$スタビライザーから成ります⁶。Steane符号に対するスタビライザーとして、

IIIXXXX \\
IXXIIXX \\
XIXIXIX \\
IIIZZZZ \\
IZZIIZZ \\
ZIZIZIZ \tag{6}

という定義をよく見かけますが（ニールセン・チャンとか）、量子ビットの順番を入れ替えれば式(5)の形になります。

これを検査行列で書くと、

H =
\begin{pmatrix}
H_x & 0 \\
0   & H_z
\end{pmatrix}  \tag{7}

H_x =
\begin{pmatrix}
1 & 0 & 0 & 1 & 0 & 1 & 1 \\
0 & 1 & 0 & 0 & 1 & 1 & 1 \\
0 & 0 & 1 & 1 & 1 & 1 & 0
\end{pmatrix}  \tag{8}

H_z =
\begin{pmatrix}
1 & 0 & 0 & 1 & 0 & 1 & 1 \\
0 & 1 & 0 & 0 & 1 & 1 & 1 \\
0 & 0 & 1 & 1 & 1 & 1 & 0
\end{pmatrix}  \tag{9}

です。後で使うので検査行列までやや詳しく述べましたが、結局、青枠の$S$ゲート直前まで来たとき全体の量子状態がどうなっているかというと、

\ket{\bar{0}}\ket{0} + \ket{\bar{1}}\ket{1}  \tag{10}

です。各項の前半のケットは$1$番目から$7$番目までの（論理）量子状態で後半のケットは$8$番目の（単体の）量子状態を表しています。この論理量子状態（前半）の方に$S$ゲートをトランスバーサルに適用します。Steane符号の論理状態に対する$S$ゲートの効果は、

\begin{align}
S^{\otimes 7} \ket{\bar{0}} &= \ket{\bar{0}} \\
S^{\otimes 7} \ket{\bar{1}} &= e^{-i\pi/2} \ket{\bar{1}} \tag{11}
\end{align}

となるので⁷、これは論理$S^{\dagger}$を適用したということに相当します。つまりトランスバーサルな$S$ゲートは論理$S^{\dagger}$ゲートになります。したがって、この直後の全体の量子ビット状態は、

\begin{align}
& \ket{\bar{0}}\ket{0} + e^{-i\pi/2} \ket{\bar{1}}\ket{1}  \\
&= (\ket{\bar{+}} + \ket{\bar{-}}) \ket{0} + e^{-i\pi/2} (\ket{\bar{+}} - \ket{\bar{-}}) \ket{1} \\
&= \ket{\bar{+}} (\ket{0} + e^{-i\pi/2} \ket{1}) + \ket{\bar{-}} (\ket{0} - e^{-i\pi/2} \ket{1}) \\
&= \ket{\bar{+}} (\ket{0} - e^{i\pi/2} \ket{1}) + \ket{\bar{-}} (\ket{0} + e^{i\pi/2} \ket{1}) \tag{12}
\end{align}

となります。この後、$1$番目から$7$番目の量子ビットに対して$X$測定をするということは論理$X$を測定するということなので、その値が$+1$（各測定値の積が$+1$）の場合$8$番目の量子状態、すなわち、$\ket{\psi}$は、

\ket{\psi} = \ket{0} - e^{i\pi/2} \ket{1}  \tag{13}

となり、$-1$の場合（各測定値の積が$-1$）、

\ket{\psi} = \ket{0} + e^{i\pi/2} \ket{1} = \ket{Y}  \tag{14}

となります。したがって、論理$X$測定値が$-1$の場合何もしないでも$\ket{Y}$が得られますが、$+1$の場合は$\ket{1}$に対する位相が反転しているのでパウリ$Z$を適用すれば$\ket{Y}$を得ることができます。というわけで、完璧な$\ket{Y}$を入力すれば完璧な$\ket{Y}$が出力されるということがわかりました。

精度の悪いマジック状態を入力した場合

それでは完璧ではない精度の悪い$\ket{Y}$を入力した場合にどうなるでしょうか。

これまで「精度の悪い」というやや曖昧な言い方をしてきたので、ここできちんと定義しておきます。精度の悪い$\ket{Y}$というのは、$\ket{Y}$を用意しようとして外界の影響等でノイズが加わりエラーが発生した状態ということです⁸。代表的なエラーとして「ビット反転エラー」「位相反転エラー」「ビット・位相反転エラー」があり、各々パウリ$X,Z,Y$ゲートが演算された状態に等しいです。実際の物理過程として現れるエラーは一般に、

\sum_{i=1}^{N} E_{i}^{\dagger} E_{i} = I \tag{15}

を満たすクラウス演算子${E_1,E_2, \dots ,E_n}$として表現することができます。ビット反転、位相反転、ビット・位相反転が各々$p_X, \space p_Z, \space p_Y$で起きる場合のエラー集合（クラウス演算子）は、

\begin{align}
E_0 &= \sqrt{p_{I}} I = \sqrt{1-p_X-p_Y-p_Z} I \\
E_1 &= \sqrt{p_{X}} X \\
E_2 &= \sqrt{p_{Y}} Y \\
E_3 &= \sqrt{p_{Z}} Z \tag{16}
\end{align}

です。ここで、$E_0$は何もエラーが起きないことを表している演算子です。この物理過程によって1量子ビットの密度演算子$\rho$は、

\rho \rightarrow \sum_{i=0}^{3} E_{i} \rho E_{i}^{\dagger}  \tag{17}

のように変化します。最初の$\rho$が純粋状態$\ket{\Psi}$だった場合は、確率$p_I$で$\ket{\Psi}$、確率$p_X$で$X \ket{\Psi}$、確率$p_Y$で$Y \ket{\Psi}$、確率$p_Z$で$Z \ket{\Psi}$に量子状態が変化します。というわけで、ここでは精度の悪い$\ket{Y}$というものを純粋状態$\ket{Y}$に対して確率${p_X, p_Y, p_Z}$で${X,Y,Z}$が演算された混合状態として定義することにします。つまり、

\ket{\Psi} \bra{\Psi} \rightarrow p_I \ket{\Psi} \bra{\Psi} + p_X X \ket{\Psi} \bra{\Psi} X + p_Y Y \ket{\Psi} \bra{\Psi} Y + p_Z Z \ket{\Psi} \bra{\Psi} Z  \tag{18}

のように変化した右辺のような混合状態として定義することにします。

さて、$\ket{Y}$を使って実現した$S$ゲートの図を上で示しましたが、$\ket{Y}$の精度が悪かった場合の効果について考えてみます。つまり、入力側の$\ket{Y}$が$X \ket{Y}, \space Y \ket{Y}, \space Z \ket{Y}$だった場合の出力状態がどうなるかについてです。少し計算すればわかりますが、各々$ZS \ket{Y}, \space S \ket{Y}, \space ZS \ket{Z}$となります。$X$または$Z$が$\ket{Y}$に加わった場合、青枠の$S$ゲートの後に$Z$が加わったことと同じ効果を及ぼしますが、$Y$が加わった場合は青枠の$S$ゲートの後に何の効果も及ぼしません。したがって、この蒸留回路で考慮しないといけないのは$X$エラーと$Z$エラーのみということになります。が、どっちのエラーが来たとしても青枠$S$ゲート直後には$Z$エラーにすり替わります。

精度の悪い$\ket{Y}$を入力状態にした回路の例を以下に示します。

上図はほんの一例であり、これ以外にも起こり得るエラーの場所のバリエーションが沢山あって、それらが確率的に発生しているとイメージしておいてください。これらのエラー確率がどのように出力側$\ket{\psi}$に伝わるのかをこれから解析していこうと思います。

エラー確率の解析

まず、上図の蒸留回路の$1$番目から$7$番目までの量子ビットはSteane符号で符号化されている状態になっていることを思い出しましょう。トランスバーサルな$S$ゲートは論理$S^{\dagger}$ゲートだったので、それを通ったとしてもSteane符号の符号空間の中の論理状態になっています。その論理状態に対して$Z$エラーがいろんな形で加わることになるわけです。$Z$エラーが加わることで起き得る状況は以下のように$3$つのパターンに分かれます。

• [パターン1] 論理状態は変化しない
• [パターン2] 論理状態は変化するが符号空間内の状態になっている（つまり何らかの論理演算子が掛かった状態）
• [パターン3] 符号空間からはみ出している

ここで、マジック状態蒸留の手法においては、[パターン3]の符号空間からはみ出してしまった場合というのは無視します、ということに注意しておきます。つまり、こういう状態が発生したら破棄してやり直すわけです。なんとも乱暴なと思われるかもしれませんが、ここでゴニョゴニョと誤り訂正の回路をもってくるよりも、潔くやり直したほうがよっぽど楽ちんという考え方なのだと思います。具体的に、$Z$エラーが加わった結果が符号空間からはみ出しているかどうかは、$X$スタビライザーを測定してみればわかります⁹。式(5)から$X$スタビライザーは$X_1 X_4 X_6 X_7, \space X_2 X_5 X_6 X_7, \space X_3 X_4 X_5 X_6 $でした。トランスバーサルな$X$測定の測定値を${m_1, m_2, m_3, m_4, m_5, m_6, m_7 } \space (m_i = -1,1)$としたとき、$m_1 m_4 m_6 m_7, \space m_2 m_5 m_6 m_7, \space m_3 m_4 m_5 m_6 $の値の中に一つでも$-1$があったら、その状態を破棄してやり直すということです。ということで、割り切ってしまったので考えるべきなのは[パターン1]と[パターン2]になります。

エラー確率を定量的に解析したいので、まず[パターン1]または[パターン2]が発生する確率（つまり[パターン3]が発生しない確率）について考えてみます。まず$7$個の量子ビットに加わるエラーを

E(\boldsymbol{c}) = \prod_{i=1}^{7} Z_{i}^{(\boldsymbol{c})_{i}}  \tag{19}

と定義しておきます。ここで$\boldsymbol{c}$は$\{0,1\}^{\otimes 7}$の要素で各量子ビットに$Z$エラーがあるかどうかを表すバイナリベクトルです。さらに、（古典）符号理論でよく使われる「重み分布母関数(weight enumerator)」というものを導入します。それは、$x,y$を何らかの実数値として、

W_{C}(x,y) = \sum_{\boldsymbol{c} \in C} x^{N-wt(\boldsymbol{c})} y^{wt(\boldsymbol{c})}  \tag{20}

で定義されるものです。ここで、$C$は符号長$N$の古典的な意味での何らかの符号空間を表します。つまり、その符号空間を定義する生成行列の列ベクトルを基底として構成されるベクトル空間です¹⁰。$wt$は重み関数で、バイナリベクトルに含まれる$1$の個数を返す関数です。でも、これ、ちょっと謎な表式ですよね。何が便利なのかさっぱり意味がわかりません。ということで、少し横道にそれて噛み砕いて説明します。

$C$に含まれる符号を全部もってきて、そこに含まれる$1$の個数が$i$である符号数を$A_i$と表すことにすると、式(20)は、

W_{C}(x,y) = \sum_{i=0}^{N} A_{i} x^{N-i} y^{i}  \tag{21}

のように表現を変えることができます（$N$は今の場合$7$です）。ここで、変数$x,y$を2種類の符号「$0$」と「$1$」の各々が生起する確率だと思うことにします。どうでしょうか。じんわりとわかってきますよね。$W_{C}(x,y)$というのは、$2$の$N$乗個ある$N$次元バイナリベクトルの要素（符号語）のうち、この符号空間$C$の要素になっている割合を表す量になっています。デタラメに$N$個の$0$,$1$からなるベクトルを作ったときにそれが$C$の要素（符号語）になっている確率と言い換えても良いです。「$1$」が生起する確率を$p$とすると、

W_{C}(1-p, p) = \sum_{i=0}^{N} A_{i} (1-p)^{n-i} p^{i}  \tag{22}

と書き換えることができますので、よりわかりやすい表式になりますね。この$A_i$は「重み分布(weight distribution)」と呼ばれています。

さて、横道は以上ということにして、何をやりたかったかというと、[パターン1]または[パターン2]が発生する確率でした。その確率$p_{pass}$は、

p_{pass} = W_{V_{H_x}^{\perp}}(1-p, p)  \tag{23}

と表すことができます。さて、いきなりですが大丈夫でしょうか？ここで、$V_{H_x}$は式(8)の$H_x$の行ベクトルを基底としてその線形結合として構成される符号空間であり、$V_{H_x}^{\perp}$はその双対符号です。つまり、$V_{H_x}^{\perp}$は$H_x$の行ベクトルと直交する基底で張られる符号空間を表しています。それらの（古典）符号を使って$E(\boldsymbol{c})$を作るとするなら、それは$X$スタビライザーと可換な$Z$エラーを表すことになります¹¹。$X$スタビライザーと可換な$Z$エラーということは、そのエラーが加わっても（量子的な意味での）符号空間からはみ出さないということです。したがって、上の式(23)は$X$スタビライザーと可換な$Z$エラーが発生する確率ということになり、（量子的な意味での）符号空間からはみ出さない確率を表しています。つまり、式(23)は$Z$エラーが加わった結果、その状態が[パターン1]または[パターン2]になる確率を表しています。どうでしょう。わかりましたでしょうか？多分一読ではなかなか腹落ちしないと思いますので、じっくりと反芻してみてください。

ここで、式(23)を$p$の関数として具体的に表すために、（古典）符号理論でよく知られた「MacWilliamsの恒等式(MacWilliams relation)」

W_{C^{\perp}}(x,y) = \frac{1}{|C|} W_{C}(x+y, x-y)  \tag{24}

を召喚します。立て続けに謎な表式登場でうんざりって感じになりそうなので、この証明は後回しにします。一旦こういうもんだということを認めてください。そうすると、$p_{pass}$は、

\begin{align}
p_{pass} &= W_{V_{H_x}^{\perp}}(1-p, p) = \frac{1}{|V_{H_x}|} W_{V_{H_x}}(1, 1-2p) \\
&= \frac{1}{2^3} (\sum_{i=0}^{7} A_{i} (1-2p)^{i})  \tag{25}
\end{align}。

のように変形できます。ここで、重み分布$A_{i}$がわかれば良いのですが、これは地道にカウントすると、

\begin{align}
A_0 &= 1 \\
A_4 &= 7 \\
A_1 &= A_2 = A_3 = A_5 = A_6 = A_7 = 0 \tag{26}
\end{align}

となります¹²。

そうすると、式(25)は、

p_{pass} = \frac{1 + 7 (1-2p)^{4}}{8}  \tag{27}

となります。これで、[パターン1]または[パターン2]が発生する確率がわかりました。

次に[パターン1]と[パターン2]の中身について考えます。両方ともSteane符号の符号空間内の状態なのですが[パターン2]の方は論理演算子がかかっていて元の論理状態からずれています。その結果、正しい$\ket{Y}$が得られません。かと言って、こうなったという判定もできないので、これは一旦出力するしかないのですが、この[パターン2]が発生する確率が$p$よりも小さくなってくれていれば、入力の精度の悪さが多少改善されたと考えるわけです。そしてこの出力を入力にして再度蒸留回路に通せば、さらに精度が向上していきます。これを繰り返すことで十分な精度をもったのものが得られます。という仕掛けがマジック状態蒸留の正体です。

では、実際に本当に入力のエラー確率$p$が出力側で改善されるかどうか確認してみます。具体的に[パターン2]が発生する確率$p_{fail}$を計算してみます。$E(\boldsymbol{c})$の$\boldsymbol{c}$が$H_z$で張られる（古典的な意味での）符号空間に含まれない確率を計算すれば良いので¹³、

p_{fail} = W_{V_{H_z}} (p, 1-p) = \frac{1}{|V_{H_z}^{\perp}|} (1, 2p-1) \tag{28}

の右辺を計算すれば良いです（$2$番めの等号でMacWilliamsの恒等式を使いました）。$V_{H_z}^{\perp}$は$H_z$と直交する空間なので、その基底は以下の行列の$4$つの行ベクトルで表せます。

H_{z}^{\perp} =
\begin{pmatrix}
1 & 0 & 0 & 1 & 0 & 1 & 1 \\
0 & 1 & 0 & 0 & 1 & 1 & 1 \\
0 & 0 & 1 & 1 & 1 & 1 & 0 \\
1 & 1 & 1 & 1 & 1 & 1 & 1 \\
\end{pmatrix}  \tag{29}

これの重み分布は、

\begin{align}
A_0 &= 1 \\
A_3 &= 7 \\
A_4 &= 7 \\
A_7 &= 1 \tag{30}
\end{align}

で、これ以外はすべて$0$なので、式(28)は、

p_{fail} = \frac{1 + 7(2p-1)^{3} + 7(2p-1)^{4} + (2p-1)^{7}}{16} \tag{31}

となります。これで、[パターン2]が発生する確率がわかりました。

[パターン1]または[パターン2]が起きた前提で[パターン2]が起きる確率は$p^{\prime} = p_{fail}/p_{pass}$で計算できます。これが入力のエラー確率$p$よりも小さくなってくれていれば、蒸留回路がうまく機能するということになります。というわけで、$p^{\prime}$を計算すると、

p^{\prime} = \frac{1 + 7(2p-1)^{3} + 7(2p-1)^{4} + (2p-1)^{7}}{2 (1 + 7 (1-2p)^{4})}  \tag{32}

となり、$p^{\prime} < p$という条件を式(32)に当てはめて数値計算すると、

p < 0.293...  \tag{33}

となります。つまり、入力のエラー確率が閾値$0.293...$よりも小さくなっていれば、この蒸留回路はうまく機能するということになります。ちなみに、この数値計算は、以下のPythonプログラムで簡単に計算できます¹⁴。

from scipy import optimize
def func(p):
    u, v = 2 * p - 1, 1 - 2 * p
    return (1 + 7 * u**3 + 7 * u**4 + u**7) / (2 * (1 + 7 * v**4)) - p
print(optimize.fsolve(func,0.2))

さらに、$p^{\prime}$が$p$に対してどのくらいのオーダーの量なのか見てみます。$p_{fail}$の$0$次、$1$次、$2$次の項を丹念に計算するとすべて$0$になって、$3$次の項の係数が$7$となるので、

p_{fail} = 7p^{3} + O(p^4)  \tag{34}

となり、それに、

\frac{1}{p_{pass}} = 1 + O(p)  \tag{35}

が掛かっているとということなので、結局、

p^{\prime} = 7p^{3} + O(p^4)  \tag{36}

ということになります。これを$2$回繰り返すと$7(7p^3)^3 = 7^{4} p^{9}$、$3$回繰り返すと$7(7^{4} p^{9})^3 = 7^{13} p^{27}$、$n$回繰り返すと$p^{3^n}$のオーダーでエラー確率は急激に減少します。例えば、最初のエラー確率が$1%$だったとして蒸留を$3$回繰り返すと、エラー確率は$10^{-15}$くらいになります。

改めてマジック状態蒸留の手順をまとめます。

• [step.1] 大量の$\ket{Y}$（エラー確率は$p$はある閾値$0.293...$以下になっている必要があります）を発生させて、7個ずつ順に蒸留回路に通します。
• [step.2] 出力側の7個の測定値${m_1,m_2,m_3,m_4,m_5,m_6,m_7}$を見て、$m_1 m_4 m_6 m_7, \space m_2 m_5 m_6 m_7, \space m_3 m_4 m_5 m_6$のどれかに$-1$が含まれている場合、結果を破棄してやり直します。破棄されなかったものに対して$m_1 m_2 m_3 m_4 m_5 m_6 m_7$の値が$-1$だった場合、$\ket{\psi}$をそのまま出力し、$+1$だった場合パウリ$Z$を演算したものを出力します。
• [step.3] この出力を沢山集めて、7個ずつ順に蒸留回路に通します。
• [step.4] [step.2]に戻ります。

という具合に繰り返します。

マジック状態蒸留の量子回路：|A>

次に$\ket{A}$の蒸留についてです。その回路を示します。

ここで、右下の$\ket{\psi}$が蒸留されて出てくる$\ket{A}$を表しており、青枠で示した$15$個の$T^{\dagger}$ゲートは$\ket{A}$、CNOTゲート、パウリ$X$ゲート、$S^{\dagger}$および$Z$測定とを使って以下のように構成したものです¹⁵。

この回路からわかる通り、$15$個の$\ket{A}$を入力して$1$個の$\ket{A}$を出力する格好になっています。$\ket{Y}$を蒸留したときと同様に、最初に大量の$\ket{A}$を用意しておいて繰り返しこの蒸留回路を通すことで十分に精度の高い$\ket{A}$を得ることができます。

完璧なマジック状態を入力した場合

さて、$\ket{Y}$の蒸留回路のときと同様、一見したところ何をやっている回路だかさっぱりわからないので再び噛み砕きます。入力の$\ket{A}$が完璧な精度だった場合、右下から完璧な精度の$\ket{A}$が出力される回路であるということを最初に確認しておきます。一番下の$2$つの量子ビットの入力状態$\ket{0}\ket{+}$にCNOTがかかっているので、まずベル状態$\ket{0}\ket{0} + \ket{0}\ket{0}$が生成されます。その後、そのベル状態の上の方の量子ビット（上から数えて$15$番目の量子ビット）と$1$番目から$14$番目の量子ビットとがCNOTゲートで複雑にエンタングルされています。実はこれ、$15$番目の量子ビットをReed-Muller符号で符号化する回路になっています。Reed-Muller符号に対するスタビライザーを検査行列で書くと、以下のようになります。

H =
\begin{pmatrix}
H_x & 0 \\
0   & H_z \\
\end{pmatrix}  \tag{37}

H_x =
\begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1 & 0 & 1 \\
0 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 1 \\
0 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 \\
0 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 1 & 1 & 1 \\
\end{pmatrix}  \tag{38}

H_z =
\begin{pmatrix}
0 & 0 & 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
1 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
1 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \\
1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 \\
1 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 \\
1 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
\end{pmatrix}  \tag{39}

Reed-Muller符号で符号化した$\ket{0}, \space \ket{1}$を各々$\ket{\bar{0}}, \space \ket{\bar{1}}$とおくと、トランスバーサルな$T^{\dagger}$直前の状態は、

\ket{\bar{0}} \ket{0} + \ket{\bar{1}} \ket{1}  \tag{40}

となります。これに$T^{\dagger}$をトランスバーサルに適用します。Reed-Muller符号の論理状態に対する$T$ゲートの効果は、

\begin{align}
(T^{\dagger})^{\otimes 15} \ket{\bar{0}} &= \ket{\bar{0}} \\
(T^{\dagger})^{\otimes 15} \ket{\bar{1}} &= e^{i\pi/4} \ket{\bar{1}} \tag{41}
\end{align}

となるので¹⁶、これは論理$T$を適用したことに相当します。つまりトランスバーサルな$T^{\dagger}$ゲートは論理$T$ゲートになります。したがって、この直後の全体の量子ビット状態は、

\begin{align}
& \ket{\bar{0}}\ket{0} + e^{i\pi/4} \ket{\bar{1}}\ket{1}  \\
&= (\ket{\bar{+}} + \ket{\bar{-}}) \ket{0} + e^{i\pi/4} (\ket{\bar{+}} - \ket{\bar{-}}) \ket{1} \\
&= \ket{\bar{+}} (\ket{0} + e^{i\pi/4} \ket{1}) + \ket{\bar{-}} (\ket{0} - e^{i\pi/4} \ket{1}) \tag{42}
\end{align}

となります。この後、$1$番目から$15$番目の量子ビットに対して$X$測定をするということは論理$X$を測定するということなので、その値が$+1$（各測定値の積が$+1$）の場合$16$番目の量子状態、すなわち、$\ket{\psi}$は、

\ket{\psi} = \ket{0} + e^{i\pi/4} \ket{1} = \ket{A} \tag{43}

となり、$-1$の場合（各測定値の積が$-1$）、

\ket{\psi} = \ket{0} - e^{i\pi/4} \ket{1}  \tag{44}

となります。したがって、論理$X$測定値が$+1$の場合何もしないでも$\ket{A}$が得られますが、論理$X$測定値が$-1$の場合は$\ket{1}$に対する位相が反転しているのでパウリ$Z$を適用すれば$\ket{A}$を得ることができます。というわけで、完璧な$\ket{A}$を入力すれば完璧な$\ket{A}$が出力されるということがわかりました。

精度の悪いマジック状態を入力した場合

それでは、完璧ではない精度の悪い$\ket{A}$を入力した場合にどうなるか考えてみます。先ほどと同様に$X$エラー、$Y$エラー、$Z$エラーが$\ket{A}$に確率的に加わった状態を考えたいところなのですが、ここでは簡単のため（というか$Y$エラー、$Z$エラーの解析は面倒そうなので）$Z$エラーについてのみ考えることにします¹⁷。

$Z$エラーが$\ket{A}$に加わった効果は上図の$T^{\dagger}$の回路から容易にわかる通り$\ket{\chi}$にパウリ$Z$が加わった状態となります。したがって、精度の悪い$\ket{A}$を入力状態にした場合、蒸留回路は例えば以下のようになります。

先ほどと同様、上図はほんの一例です。これ以外にも起こり得るエラーの場所が沢山あって、それらが確率的に発生しているとイメージしておいてください。これらのエラー確率がどのように出力側に伝わるのかをこれから解析していこうと思います。、

エラー確率の解析

先ほどと同様に$Z$エラーが加わったことで起き得る状況は、

• [パターン1] 論理状態は変化しない
• [パターン2] 論理状態は変化するが符号空間内の状態になっている（論理演算子が掛かった状態）
• [パターン3] 符号空間からはみ出している

です。このうち[パターン3]になったかどうかは、$X$スタビライザーを測定してみればわかります。式(38)から$X$スタビライザーは$X_1 X_6 X_7 X_{10} X_{11} X_{12} X_{13} X_{15}$ , $X_2 X_5 X_7 X_9 X_{11} X_{12} X_{14} X_{15}$, $X_3 X_4 X_7 X_9 X_{10} X_{13} X_{14} X_{15}$ , $X_3 X_5 X_6 X_8 X_{11} X_{13} X_{14} X_{15}$の$4$つです。トランスバーサルな$X$測定の測定値を${m_1, m_2, m_3, m_4, m_5, m_6, m_7, m_8, m_9, m_{10}, m_{11}, m_{12}, m_{13}, m_{14}, m_{15}} \space (m_i = -1,1)$としたとき、$m_1 m_6 m_7 m_{10} m_{11} m_{12} m_{13} m_{15}$, $m_2 m_5 m_7 m_9 m_{11} m_{12} m_{14} m_{15}$ , $ m_3 m_4 m_7 m_9 m_{10} m_{13} m_{14} m_{15}$, $m_3 m_5 m_6 m_8 m_{11} m_{13} m_{14} m_{15}$の値の中に一つでも$-1$があったら、その状態を破棄してやり直します。

これで[パターン1]または[パターン2]が得られるのですが、それが発生する確率$p_{pass}$は、

\begin{align}
p_{pass} &= W_{V_{H_x}^{\perp}}(1-p, p) = \frac{1}{|V_x|} W_{V_{H_x}}(1, 1-2p) \\
&= \frac{1}{2^4} (\sum_{i=0}^{15} A_{i} (1-2p)^{i})  \tag{45}
\end{align}

で計算できます。ここで、$A_{i}$がわかれば良いのですが、地道にカウントすると、

\begin{align}
A_0 &= 1 \\
A_8 &= 15 \tag{46}
\end{align}

となります（上記以外の$A_i$はすべて$0$）。そうすると、式(45)は、

p_{pass} = \frac{1 + 15 (1-2p)^{8}}{16}  \tag{47}

となります。これで、[パターン1]または[パターン2]が発生する確率$p_{pass}$がわかりました。

次に、[パターン2]が発生する確率$p_{fail}$ですが、$E(\boldsymbol{c})$の$\boldsymbol{c}$が$H_z$で張られる（古典的な意味での）符号空間に含まれない確率を計算すれば良いので、

p_{fail} = W_{V_{H_z}} (p, 1-p) = \frac{1}{|V_{H_z}^{\perp}|} (1, 2p-1) \tag{48}

です。$V_{H_z}^{\perp}$は$H_z$と直交する空間なので、その基底は以下の行列の$5$つの行ベクトルで表せます。

H_{z}^{\perp} =
\begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1 & 0 & 1 \\
0 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 1 \\
0 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 \\
0 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
\end{pmatrix}  \tag{49}

これの重み分布は、

\begin{align}
&A_0 = 1 \\
&A_7 = 15 \\
&A_8 = 15 \\
&A_{15} = 1 \tag{50}
\end{align}

で、これ以外はすべて$0$なので、式(48)は、

p_{fail} = \frac{1 + 15(2p-1)^{7} + 15(2p-1)^{8} + (2p-1)^{15}}{32} \tag{51}

となります。これで、[パターン2]が発生する確率$p_{fail}$がわかりました。

あとは、$p^{\prime} = p_{fail}/p_{pass}$が入力のエラー確率$p$よりも小さくなってくれていれば、蒸留回路がうまく機能するということなので、$p^{\prime}$を計算すると、

p^{\prime} = \frac{1 + 15(2p-1)^{7} + 15(2p-1)^{8} + (2p-1)^{15}}{2 (1 + 15 (1-2p)^{8})}  \tag{52}

となり、$p^{\prime} < p$という条件を上の式に当てはめて数値計算すると、

p < 0.141...  \tag{53}

が得られます。つまり、入力のエラー確率が閾値$0.141...$よりも小さくなっていれば、この蒸留回路はうまく機能するということになります。ちなみに、この数値計算は、以下のPythonプログラムで簡単に計算できます。

from scipy import optimize
def func(p):
    u, v = 2 * p - 1, 1 - 2 * p
    return (1 + 15 * u**8 + 15 * u**7 + u**15) / (2 * (1 + 15 * v**8)) - p
print(optimize.fsolve(func,0.2))

さらに、$p^{\prime}$が$p$に対してどのくらいのオーダーの量なのか見てみます。$p_{fail}$の$0$次、$1$次、$2$次の項を丹念に計算するとすべて$0$になって、$3$次の項の係数が$35$となるので、

p^{\prime} = 35p^{3} + O(p^4)  \tag{54}

したがって、この蒸留回路を$n$回繰り返すと$p^{3^n}$のオーダーでエラー確率は急激に減少することがわかります。

MacWilliamsの恒等式

それでは、先ほど後回しにすると言っていた「MacWilliamsの恒等式」を証明します¹⁸。

$C$を$[n,k]$符号とし$C^{\perp}$を$C$の双対符号¹⁹とすると、両者の重み分布母関数の間に以下のような関係が成り立ちます。これをMacWilliamsの恒等式と呼びます。

W_{C^{\perp}}(x,y) = \frac{1}{|C|} W_{C}(x+y, x-y)  \tag{55}

【証明】

$H$を符号$C$のパリティ検査行列とします。$H_{ext}$を$H$の行のあらゆる線形結合からなる行列、すなわち、双対符号$C^{\perp}$のすべての符号語を行としてもつ行列とします²⁰。この$H_{ext}$を使って

x \in \boldsymbol{F}_{2}^{n}

に対する拡張シンドローム$s_{ext} = (s_1, s_2, \dots s_{2^{n-k}})$を以下で定義します²¹。

s_{ext} = H_{ext} x  \tag{56}

このとき、$x \in C$とすると、

s_{ext} = (s_1, s_2, \dots s_{2^{n-k}}) = (0,0, \dots 0)  \tag{57}

のようにすべての拡張シンドロームの値は$0$になります。一方、$x \notin C$とすると、

s_{ext} = (s_1, s_2, \dots s_{2^{n-k}}) = (1,0,0,1, \dots 1)  \tag{58}

のように拡張シンドロームに$1$が混じることになるのですが、どんな場合でも$0$と$1$の個数は各々$2^{n-k-1}$になるという性質があります。ん？ほんと？という声が聞こえてきそうなので、なぜそうなるかを説明します。$x \notin C$である$x$を一つもってきて、

C^{\prime} = \{y \in C^{\perp} \space | \space y^{T} x = 0\}  \tag{59}

のような集合を定義します。そうするとこれは$C^{\perp}$の部分集合になっており、その次元は$C^{\perp}$の次元よりも一つ小さい$n-k-1$になっています。$N$次元ベクトル空間を考えてその中の一つのベクトルに直交する空間を考えると、それは次元が一つ小さくなった$N-1$次元の部分空間になっていますよね。それと同じことです。今の場合ベクトルを構成する要素は$0$または$1$の2種類しかないのでベクトル空間を構成する要素の数を数え上げることができます。すなわち、$C^{\perp}$の中で$x \notin C$と直交するもの（拡張シンドローム値が$0$）の個数は$2^{n-k-1}$個あり、直交しないもの（拡張シンドローム値が$1$）の個数は$2^{n-k-1}$個あるということになります。ということで、式(58)の右辺を構成する$0$と$1$の数は半々になることがわかりました。

さて、そうしたときに、$E_j$という集合を以下のように定義します。

E_j = \{x \in \boldsymbol{F}_{2}^{n} \space | \space s_j = 0 \}, \space j = 1,2, \dots 2^{n-k}  \tag{60}

つまり、$C^{\perp}$の$j$番目の基底（= $H_{ext}$の$j$行目の行ベクトル）と直交する符号語からなる集合を$E_j$と定義します。$C$は$C^{\perp}$に含まれるすべての符号語と直交するので、この$E_j$を使って、

C = E_1 \cap E_2 \cap \dots \cap E_{2^{n-k}}  \tag{61}

のように表現し直すことができます。$E_j$の補集合、

\bar{E_j} = \boldsymbol{F}_{2}^{n} \verb|\| E_j = \{ x \in \boldsymbol{F}_{2}^{n} \space | \space s_j = 1 \}  \tag{62}

を使うと式(61)は、

C = \boldsymbol{F}_{2}^{n} \verb|\| (\bar{E_1} \cup \bar{E_2} \cup \dots \cup \bar{E_{2^{n-k}}})  \tag{63}

と書くこともできます（集合論の基本法則より）。ここで、右辺の和集合を、

\bar{E} = \bar{E}_1 \cup \bar{E}_2 \cup \dots \cup \bar{E}_{2^{n-k}}  \tag{64}

と書くことにします。$\bar{E}$は$C$の補集合なのでその重み分布母関数は、$\boldsymbol{F}_{2}^{n}$に対する重み分布関数、

W_{\boldsymbol{F}_{2}^{n}} = (x + y)^{n}  \tag{65}

を使って、

W_{C}(x,y) = (x+y)^{n} - W_{\bar{E}}(x,y)  \tag{66}

のように表すことができます。

ここで、式(58)の$s_{ext}$に含まれる$0$と$1$の個数は各々$2^{n-k-1}$個だったことを思い出しましょう。そして、$\bar{E}$は式(64)で定義されている通り$\bar{E}_j$の和集合でした。ということは、$\bar{E}$に含まれる任意の一つの符号語を取り上げると、それは$2^{n-k-1}$個の$\bar{E}_j$に含まれるということが言えます²²。とすると$\bar{E}$と$\bar{E}_j$の重み分布母関数には以下のような関係式が成り立つことがわかります²³。

2^{n-k-1} W_{\bar{E}}(x,y) = \sum_{j=1}^{2^{n-k}} W_{\bar{E}_j}(x,y)  \tag{67}

ここで、右辺について考えます。$H_{ext}$の$j$番目の符号語のハミング重み（$1$の個数）を$w_j$とすると、$W_{\bar{E}_j}(x,y)$は、

W_{\bar{E}_j}(x,y) = (x+y)^{n-w_j} \sum_{k=1,odd}^{w_j}
\begin{pmatrix}
w_j \\
k
\end{pmatrix} x^{w_{j}-k} y^{k}  \tag{68} 

と書けます。これ、わかりますでしょうか？$H_{ext}$の$j$番目の符号語の中で$w_j$個の桁に$1$が並んでいます。$\bar{E_j}$であるためにはその桁と同じ桁に奇数個の$1$が並んでいれば良いです（内積が$1$になるためには）。その他$n-w_{j}$個の桁については$0$でも$1$でもどっちでも良いです。ということを重み母関数の式で表すと式(68)のようになります。この式は以下のように変形することができます。

\begin{align}
W_{\bar{E}_j}(x,y) &= (x+y)^{n-w_j} \frac{1}{2} \bigl[ (x+y)^{w_j} - (x-y)^{w_j} \bigr] \\
&= \frac{1}{2} (x+y)^{n} - \frac{1}{2} (x+y)^{n-w_j} (x-y)^{w_j}  \tag{69}
\end{align}

これを式(67)に代入すると、

\begin{align}
2^{n-k-1} W_{\bar{E}}(x,y) &= 2^{n-k-1} (x+y)^{n} - \frac{1}{2} \sum_{j=1}^{2^{n-k}} (x+y)^{n-w_j} (x-y)^{w_j}  \\
W_{\bar{E}}(x,y) &= (x+y)^{n} - \frac{1}{2^{n-k}} \sum_{j=1}^{2^{n-k}} (x+y)^{n-w_j} (x-y)^{w_j}  \\
&= (x+y)^{n} - \frac{1}{|C^{\perp}|} W_{C^{\perp}}(x+y, x-y)  \tag{70}
\end{align}

さらに、これを式(66)に代入すると、

W_{C}(x,y) = \frac{1}{|C^{\perp}|} W_{C^{\perp}}(x+y, x-y)  \tag{71}

となります。$C$と$C^{\perp}$の立場を入れ替えると、

W_{C^{\perp}}(x,y) = \frac{1}{|C|} W_{C}(x+y, x-y)  \tag{72}

となり、MacWilliamsの恒等式が導けました。（証明終）

シミュレーション

マジック状態蒸留のプロセスが理解できたので、量子計算シミュレータqlazyを使ってエラー確率が理論通りに減少するのかを確認してみます。

実装

$\ket{Y}$の蒸留シミュレーションを以下のようにPythonで実装しました²⁴。$p$の確率を適当に変化させ各々500回試行したときのエラー頻度から確率$p^{\prime}$を計算してmatplotlibでプロットしています。ソースコードはここにあります。

import random
import numpy as np
import matplotlib.pyplot as plt
import logging
from qlazy import QState
from qlazy.tools.Register import CreateRegister, InitRegister

EPS = 1.e-8
QS_Y = QState(qubit_num=1).h(0).s(0)

def init_qs_list(N, prob):

    return [QState(qubit_num=1).h(0).s(0).z(0) if random.uniform(0, 1) <= prob else QState(qubit_num=1).h(0).s(0) for i in range(N)]

def measure_X_stab(mval_list):

    p0 = (mval_list[3-1] + mval_list[4-1] + mval_list[5-1] + mval_list[6-1]) % 2
    p1 = (mval_list[2-1] + mval_list[5-1] + mval_list[6-1] + mval_list[7-1]) % 2
    p2 = (mval_list[1-1] + mval_list[4-1] + mval_list[6-1] + mval_list[7-1]) % 2
    if p0 == 0 and p1 == 0 and p2 == 0: return True
    else: return False

def measure_logical_X(mval_list):

    if sum(mval_list) % 2 == 0: return True
    else: return False

def distillation_Y(qs_list):

    data = CreateRegister(1)  # data qubit
    code = CreateRegister(7)  # logical qubit for Steane code
    anci = CreateRegister(1)  # ancilla qubit for S gate
    qubit_num = InitRegister(data, code, anci)
    qs = QState(qubit_num=qubit_num-len(anci))

    # bell state
    qs.h(data[0]).cx(data[0], code[6])

    # steane code
    qs.h(code[0]).h(code[1]).h(code[2])
    qs.cx(code[6], code[3]).cx(code[6], code[4])
    qs.cx(code[2], code[3]).cx(code[2], code[4]).cx(code[2], code[5])
    qs.cx(code[1], code[4]).cx(code[1], code[5]).cx(code[1], code[6])
    qs.cx(code[0], code[3]).cx(code[0], code[5]).cx(code[0], code[6])

    # S gates and measurements
    mval_list = []
    for i in range(7):
        qs_total = qs.tenspro(qs_list[i])
        qs_total.cx(code[i], anci[0]).h(anci[0]).cx(code[i], anci[0]).h(anci[0])
        mval = int(qs_total.mx(qid=[code[i]]).last)
        mval_list.append(mval)
        qs = qs_total.partial(qid=data+code)

    if measure_X_stab(mval_list) == False:
        return None

    qs_pat = qs_total.partial(qid=data)
    if measure_logical_X(mval_list) == True: qs_pat.z(0)
    
    return qs_pat

def prob_simulation(prob, trial):

    fail = 0
    for _ in range(trial):
        while True:
            qs_list = init_qs_list(7, prob)
            qs = distillation_Y(qs_list)
            if qs is not None: break
        if abs(qs.fidelity(QS_Y) - 1.0) > EPS: fail += 1
    return fail / trial

def prob_theoretical(p):

    u, v = 2 * p - 1, 1 - 2 * p
    return (1 + 7 * u**3 + 7 * u**4 + u**7) / (2 * (1 + 7 * v**4))

if __name__ == '__main__':

    # logging
    logger = logging.getLogger('ErrorSimulation')
    logger.setLevel(20)
    sh = logging.StreamHandler()
    fh = logging.FileHandler('error_simulation_Y.log', mode='w')
    logger.addHandler(sh)
    logger.addHandler(fh)
    
    # error probability (theoretical)
    prob_in = np.linspace(0.0, 0.4, 50)
    prob_out_linear = np.array([p for p in prob_in])
    prob_out_theoretical = np.array([prob_theoretical(p) for p in prob_in])

    # error probability (simulation)
    trial = 500
    prob_in_simulation = np.linspace(0.0, 0.4, 9)
    logger.info("{0:}, {1:}, {2:}".format("in", "out(simulation)", "out(theoretical)"))
    prob_out_simulation = []
    for p_in in prob_in_simulation:
        p_out = prob_simulation(p_in, trial)
        logger.info("{0:.6f}, {1:.6f} {2:.6f}".format(p_in, p_out, prob_theoretical(p_in)))
        prob_out_simulation.append(p_out)

    # plot
    fig, ax = plt.subplots()
    ax.set_xlabel("error probability (in)")
    ax.set_ylabel("error probability (out)")
    ax.set_title("magic state distillation: |Y>")
    ax.grid()
    ax.plot(prob_in, prob_out_theoretical, color='orange', label='theoretical')
    ax.plot(prob_in, prob_out_linear, color='gray', linestyle='dashed')
    ax.plot(prob_in_simulation, prob_out_simulation, color='green', label='simultion', marker='s')
    ax.legend(loc=0)
    fig.tight_layout()
    plt.savefig('error_simulation_Y.png')
    plt.show()

長くなるので細かくは説明しません。distillation_Y関数が蒸留を実行する関数で、$7$個の量子状態（精度の悪い$\ket{Y}$）をリストにして与えることで$1$個の蒸留結果（精度の良い$\ket{Y}$）を得るものです。それをprob_simulation関数から呼んでいます。prob_simulation関数は各$p$に関して500回の蒸留を試行して確率$p^{\prime}$を得る関数です。

$\ket{A}$の蒸留シミュレーションは、$\ket{Y}$と類似しているので掲載しません。ソースコードをここに置きましたので、ご興味ある方はご参照ください。

結果

$\ket{Y}$の蒸留シミュレーション結果は以下のようになりました。

ここで、横軸は入力のエラー確率$p$、縦軸は出力のエラー確率$p^{\prime}$を表しています。オレンジ色の曲線は理論値、

p^{\prime} = \frac{1 + 7(2p-1)^{3} + 7(2p-1)^{4} + (2p-1)^{7}}{2 (1 + 7 (1-2p)^{4})}  \tag{32}

のグラフで、緑色の折れ線は上記プログラムでシミュレーションしたエラー確率をプロットしたものです。一見してわかるように理論曲線にだいたい沿ったプロットが得られています。グレーの点線は$p^{\prime}=p$の直線ですがそれとの交点が閾値$0.293$付近に来ていることもわかります。

次に、$\ket{A}$の蒸留シミュレーション結果です。以下のようになりました。

オレンジ色の曲線は理論値、

p^{\prime} = \frac{1 + 15(2p-1)^{7} + 15(2p-1)^{8} + (2p-1)^{15}}{2 (1 + 15 (1-2p)^{8})}  \tag{52}

のグラフで、緑色の折れ線はシミュレーション結果です。こちらも理論から予想される曲線付近にだいたいプロットされていて、閾値も$0.141$付近にあることが見て取れます。というわけで、両方とも蒸留シミュレーション成功！という結果になりました。

おわりに

誤り耐性のあるユニバーサル量子計算を行うためになくてはならない（とされている）マジック状態蒸留の原理が理解できました。今回理論の理解も実装も結構苦労したので、理論曲線に寄り添うプロットが得られたときは、おおっ！て感じで気持ち良かったです（当然の結果ではあるのですが）。qlazyが正しく動作している確認にもなったと思います。

さて、ショアの符号からはじまり、およそ1年半にわたり量子誤り訂正符号について勉強してきて、大体一区切り付いたかなーということで次の話題をどうするか？目下考え中です。一区切り付いたとはいえ積み残しもあるので（例えばLattice Surgeryとか）、そういった誤り訂正符号まわりのあれこれを取り上げていくか、あるいは、方向性を大きく変えて量子暗号とか量子インターネットとかの通信系に舵を切るか、はたまた、テンソルネットワークとか圏論の話とかも面白そうなので興味の赴くまま、そっちの話題をちょっとずつつまんでみたり、ということもあるかもしれません。いずれにしろ量子情報理論のネタはまだまだ尽きることはありません。例によって予定は未定ですが、引き続き悪しからず。

以上

1. Braidingを使った構成法については、「表面符号によるユニバーサル量子計算(Braiding)」 「表面符号によるユニバーサル量子計算（１）〜Braidingで論理CNOT演算〜」 「表面符号によるユニバーサル量子計算（２）〜Braidingで論理アダマール演算〜」 という一連の記事で説明しました。 ↩

2. 「表面符号によるユニバーサル量子計算（３）〜Braidingで論理位相シフト演算〜」という記事で説明しました。 ↩

3. 最終的には精度の高い論理$S$ゲート（または論理$T$ゲート）を実現できれば良いので、精度の悪い$\ket{Y}, \space \ket{A}$で作られた精度の悪い$\ket{\bar{Y}}, \space \ket{\bar{A}}$を起点にしても良いです。その場合、論理量子ゲートで構成されたマジック状態蒸留回路を通して精度の高い$\ket{\bar{Y}}, \space \ket{\bar{A}}$を得るというやり方になります。本記事では簡単のため1量子ビット状態の$\ket{Y}, \space \ket{A}$を得る手法としてマジック状態蒸留を説明しますが、単体量子ゲートを論理量子ゲートに置き換えれば全く同じ意味合いの回路になりますので、適宜読み替えてください。 ↩

4. 実際には蒸留失敗する場合もあるので（後述）、それを見越してもっと沢山の$\ket{Y}$を用意する必要があります。 ↩

5. 面倒なので規格化定数は省略しています（以下同様）。 ↩

6. $7$番目の量子ビットが$\ket{0}$だったとすると複雑なCNOT直前の最初のスタビライザー状態は${XIIIIII,IXIIIII,IIXIIII,IIIZIII,IIIIZII,IIIIIZI,IIIIIIZ}$です。スタビライザー形式でこれらCNOT後の状態を計算すると${XIIXIXX,IXIIXXX,IIXXXXI,IZZZIIZ,ZIZIZIZ,ZZZIIZI,ZZIIIIZ}$となります。XスタビライザーとZスタビライザーの格好が微妙に違っていて気持ち悪いのでZスタビライザー同士を適当に掛け合わせると${XIIXIXX,IXIIXXX,IIXXXXI,ZIIZIXX,IZIIZZZ,IIZZZZI,ZZZZZZZ}$というきれいな格好にできます。最後の$ZZZZZZZ$は論理$Z$演算子を表しています。 ↩

7. Steane符号で$\ket{\bar{0}}$を物理量子ビットで表したとき重ね合わせ状態として各項に現れる単体$\ket{1}$の個数はすべて$4$個で、$\ket{\bar{1}}$の場合はすべて$3$個になることからわかると思います。 ↩

8. ということなので、$\ket{Y}$の中の位相$\pi/2$が$\pi/2 + \alpha$みたいにズレた純粋状態ではありません。 ↩

9. 誤り訂正をやりたい場合はここで間接測定するのですが、もう訂正しないと割り切っているので直接測定しちゃっています。 ↩

10. 古典符号理論に関しては「量子情報理論の基本：量子誤り訂正（古典線形符号）」に説明がありますので、ご参照ください。 ↩

11. スタビライザーを検査行列の形式（バイナリベクトル形式）で書くと$X$スタビライザーに対するバイナリベクトルと$Z$スタビライザーに対するバイナリベクトルの内積が$0$（ただし加算は2を法とする）のとき可換になります。 ↩

12. 式(8)の$H_x$の$3$つの行ベクトルを基底にした線型結合で構成したすべての符号$2^3=8$個をずらっと書き出してみて、各々$1$の個数がいくつになるか見て頻度を数えればわかります。紙と鉛筆を使ってひたすら目視で数えるというやり方でももちろん良いです（$8$個くらいだったら簡単ですし）。が、現代人だったらPythonを使いましょう。簡単なプログラムでさくっとカウントできます。任意のバイナリベクトルの基底を入力できるように汎用的に作っておけば後々便利に使えます（はい、これ宿題です、、とか言ってみる）。 ↩

13. $\boldsymbol{c}$で表現された$Z$エラーが今考えているスタビライザー群に含まれるかどうかは、その$Z$エラーがスタビライザー（$Z$エラーなので$Z$スタビライザーだけ考えれておけば良いです）の何らかの掛け算で表すことができるかどうかで判定できます。バイナリベクトル表現で言うと、$Z$エラーに対応したバイナリベクトルが$Z$スタビライザーに対するバイナリベクトル（つまり$H_z$の行ベクトル）の線形結合で表すことができるかどうかで判定できます。つまり、$V_{H_z}$に含まれるかどうかで判定できます。式(28)は$\boldsymbol{c}$が$V_{H_z}$に含まれない場合の重み分布母関数になっています。 ↩

14. optimizeの初期値は$0.0$とか$0.5$みたいな自明な解に落ち込まないように適当に0.2にしてみました。 ↩

15. 今回示した蒸留回路は参考文献1に掲載されていたものの引用なのですが、FIG.33(b)に示されている$\ket{A}$を使って実現した$T^{\dagger}$の回路だけは、どう計算しても出力結果が合いません。もしかすると論文がミスっているのかもと思ったので、自分で独自に修正しました。もしかすると自分の方が間違っているかもしれないので、その場合はご指摘いただけるとありがたいです。 ↩

16. Reed-Muller符号で$\ket{\bar{0}}$を物理量子ビットで表したとき重ね合わせ状態として各項に現れる単体$\ket{1}$の個数はすべて$8$個で、$\ket{\bar{1}}$の場合はすべて$7$個になります。$X$スタビライザーは$H_x$の行ベクトルに対応して$4$個あるので、その各々を$(S_x^{(1)}, S_x^{(2)}, S_x^{(3)}, S_x^{(4)})$として$\ket{\bar{0}} = (I+S_x^{(1)})(I+S_x^{(1)})(I+S_x^{(1)})(I+S_x^{(1)})\ket{0}^{\otimes 15}$を地道に計算してみればわかります。 ↩

17. 先ほどの$\ket{Y}$の蒸留の場合と違って、$X$エラー、$Y$エラーが加わったときにトランスバーサル$T^{\dagger}$直後の状態がどうなるか計算してみればわかるのですが、パウリ演算子だけでなく位相シフトもエラーの効果として加わってしまい解析がちょっと面倒そうなので、ここでは一旦追求はやめて、そっと見なかったことにしておきます、汗。 ↩

18. 参考文献3のp.13-15あたりにその証明があり参考にさせていただきました。が、重み分布母関数の定義が本記事での定義（=参考文献2の定義）と違って1変数関数として定義されているため、2変数関数に読み替えて証明を構成し直してみました。 ↩

19. 「量子情報理論の基本：量子誤り訂正（古典線形符号）」で説明したように、線形符号$C$の双対符号$C^{\perp}$というのは、$C$の生成行列とパリティ検査行列を各々$G,H$としたとき、$H^{T}$を生成行列、$G^{T}$をパリティ検査行列として定義される符号です。$H^{T}G = 0$という関係が成り立つことから$C$と$C^{\perp}$の要素同士は常に直交します。線形代数に馴染みのある人は、例えば$n$次元実ベクトル空間（または$n$次元複素ベクトル空間）の部分空間に対する直交補空間のようなものをイメージしておけば良いのですが、線形符号の場合、内積は$2$を法として計算するため、少し違うイメージを持っておいた方が良いです。線形符号の要素$x$に対して内積$x^{T} x = 0$になる場合があるので$C \cap C^{\perp} = \phi$とならない場合があります。つまり、$C$を構成する基底ベクトルの集合と$C^{\perp}$を構成する基底ベクトルの集合とが一部重なっていることがあります。それどころか$C \subseteq{C^{\perp}}$だったり$C = {C^{\perp}}$だったりということも十分あり得ます。前者を「弱い自己双対」後者を「厳密な自己双対」と言います。復習でした。 ↩

20. $C$は$[n,k]$符号なので$H$は$n-k$行$n$列、$H_{ext}$は$2^{n-k}$行$n$列の行列です。 ↩

21. $\boldsymbol{F}_{2}^{n}$はすべての$n$次元バイナリベクトルからなる集合です。ただし、加算は$2$を法として行われます。よく$\{ 0,1 \}^{\otimes n}$みたいな書き方をしますが、$2$を法とした加算をするという気持ち（=「$2$元体」という体をなす）が入っていないです。その気持ちを含ませるため$\boldsymbol{F}_{2}^{n}$という記号を使うと思っておけばとりあえず良いです。 ↩

22. ちょっとわかりにくいですが頑張ってついて来てください。$\bar{E}$に含まれる任意の一つの符号語を$x$とおくとすると$x \notin C$です。そういう$x$に対する拡張シンドロームを計算すると$0$と$1$が半分ずつ現れます。拡張シンドロームを計算するということは$H$の各行ベクトル、すなわち$C^{\perp}$と$\bar{E}$の各符号語との内積を計算するということです。$\bar{E}_j$は$C^{\perp}$の$j$番目と直交しない符号語の集合だったので、ある特定の$x \notin C$は$2^{n-k-1}$個の$\bar{E}_j$の中に含まれます。 ↩

23. ランダムに選んだある符号語が$\bar{E}$に含まれる確率が$W_{\bar{E}}(x,y)$で、$\bar{E_j}$に含まれる確率が$W_{\bar{E_j}}(x,y)$だとざっくりイメージすればわかると思います。 ↩

24. 本文とは関係ないですがqlazy最新版に関するアナウンスです。以前のqlazyでは量子状態のメモリをいちいち明示的にfreeメソッドで解放しないといけない面倒な仕様になっていましたが、最新版のv0.2.1から自動解放するように変更しました。freeはもう不要です（はじめからこうしておけ！というご指摘はごもっともですが、なんかメモリ関連でバグっぽいおかしな動作をする場合があったので、、でも修正されました）。 ↩