🧙
今更学ぶ魔法状態蒸留

2026/01/12に公開
 この記事の目的2024年9月に魔法状態栽培 [1] が提案されて以降、魔法状態蒸留は魔法状態を生成する手続きとしては非効率的なものとなり、研究も下火になってきたように感じます。ただ、魔法状態栽培の基礎や量子誤り訂正符号の応用としての側面があるため、2026年現在でも勉強することに意味はあると思います。本稿では、魔法状態蒸留について数式レベルで理解するための最短経路をたどります。プログラミングによる理解の補助も行います。具体的には、Bravyi-Kitaev (2005) [2] の 5-to-1 魔法状態蒸留を数式で追える形に整理します。

 背景FTQC では非 Clifford 操作の実装が高コストです。このため、非 Clifford 資源としての魔法状態を用意し、Clifford 操作と Pauli 測定のみで目的のゲートを実現する枠組みが研究されてきました。たとえば、 T ゲートは魔法状態 \ket{M} と Clifford 操作と Pauli 測定で実装できます [3]。
魔法状態蒸留は、精度の低い魔法状態を複数用意して高精度な状態へと濃縮する手続きです。必要に応じてこれを反復し、十分に高い精度の魔法状態を得ます。

 注意点文献ごとに記号や定義が統一されていない点には注意が必要です。特に Bravyi-Kitaev では T ゲートを現代の標準的定義と異なる形で定め、しかもそれが Clifford 操作です。また蒸留手続きには多様なバリエーションがあり、蒸留対象となる魔法状態、採用する符号、回路最適化の方針が文献ごとに異なるため、「魔法状態蒸留」の回路として見たことのない回路があらわれることはよくあります。

 記号と基本事項まず Pauli 演算子を X,Y,Z とします。アダマールゲートを H とし、 \frac{\pi}{2} 位相回転ゲートを S とします。長さ n の符号が + の Pauli 文字列全体の集合を S_+(n) と書きます。 S_+(n) の要素は I,X,Y,Z の n 個のテンソル積のため、要素数は 4^n です。Pauli 文字列 g \in S_+(n) の重み |g| は、g に含まれる非 I の個数とします。また、ビット列 \mathbb{x}=x_1 x_2 \cdots x_n のハミング重みを |\mathbb{x}| と書き、後で定義する量子状態 \ket{K_0},\ket{K_1} に対して \ket{K_\mathbb{x}} = \ket{K_{x_1}} \otimes \cdots \otimes \ket{K_{x_n}} と定義します。

 魔法状態の定義魔法状態 \ket{K_0} は Clifford 操作 K の固有状態として定義されます。

K = e^{i \frac{\pi}{4}} S H = \frac{e^{i \frac{\pi}{4}}}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ i & -i \\ \end{pmatrix}
とすると、共役作用は

\begin{align*}
K X K^{\dagger} &= Z \\
K Z K^{\dagger} &= Y \\
K Y K^{\dagger} &= X
\end{align*}
となります。K の固有値 e^{i \frac{\pi}{3}} の固有状態を \ket{K_0} と書き、固有値 e^{- i \frac{\pi}{3}} の固有状態を \ket{K_1} と書きます。密度行列は次の形になります。

\begin{align*}
\ket{K_0}\bra{K_0} &= \frac{1}{2}\left(I + \frac{1}{\sqrt{3}}(X+Y+Z)\right) \\
\ket{K_1}\bra{K_1} &= \frac{1}{2}\left(I - \frac{1}{\sqrt{3}}(X+Y+Z)\right)
\end{align*}
!Bravyi-Kitaev はこの K を T と定義しています。また、T に関連した魔法状態を T 型の魔法状態と呼んでいます。そのほか、アダマール操作 H に関連した魔法状態を H 型の魔法状態と呼んでいます。

 [[5,1,3]] perfect codeBravyi-Kitaev の魔法状態蒸留は [[5,1,3]] perfect code [4] を利用した手続きです。まずはこの量子誤り訂正符号についていくつかの性質を証明します。

 定義スタビライザー生成元を

\begin{align*}
S_1 &= X Z Z X I \\
S_2 &= I X Z Z X \\
S_3 &= X I X Z Z \\
S_4 &= Z X I X Z
\end{align*}
とします。S_1 S_2 S_3 S_4 = Z Z X I X からも対称性が見て取れます。スタビライザー群を G、スタビライザー空間を L、射影演算子を \Pi とします。
!補題①：射影演算子の数式表示

\Pi = \frac{1}{16} \prod_{j=1}^{4} (I + S_j) = \frac{1}{16} \sum_{h \in G} h
証明は次の通りです。\frac{1}{2}(I+S_j) は S_j の +1 固有空間への射影であり、L は S_1,\dots,S_4 の同時 +1 固有空間なので、積を取れば L への射影になります。さらに G の要素は S_1^{e_1} S_2^{e_2} S_3^{e_3} S_4^{e_4} で生成されるため、積の展開が \sum_{h \in G} h になります。
スタビライザー群の要素をすべて列挙してみましょう。
stabilizer.py
import itertools

from qiskit.quantum_info import Pauli

stabilizer_generators = ["XZZXI", "IXZZX", "XIXZZ", "ZXIXZ"]

stabilizer_group = []
for exps in itertools.product([0, 1], repeat=len(stabilizer_generators)):
    pauli = Pauli("IIIII")
    for gen, exp in zip(stabilizer_generators, exps):
        if exp == 1:
            pauli *= Pauli(gen)
    stabilizer_group.append(pauli.to_label())
stabilizer_group = sorted(list(set(stabilizer_group)))

print("stabilizer group size:", len(stabilizer_group))
print("stabilizer group:")
for h in stabilizer_group:
    print(h)
stabilizer group size: 16
stabilizer group:
IIIII
IXZZX
IYXXY
IZYYZ
XIXZZ
XXYIY
XYIYX
XZZXI
YIYXX
YXXYI
YYZIZ
YZIZY
ZIZYY
ZXIXZ
ZYYZI
ZZXIX
この出力から、非自明な 15 要素がすべて重み 4 であることが分かります。

 性質１（パウリの transversality）\hat{X}=X^{\otimes 5}, \hat{Y}=Y^{\otimes 5}, \hat{Z}=Z^{\otimes 5} は論理演算子です。実際、\hat{X},\hat{Y},\hat{Z} は S_1,\dots,S_4 と可換であり、しかもスタビライザー群の要素は重み 0 または 4 に限られるため、重み 5 のこれらは群に含まれません。

スタビライザー空間の基底を \hat{X},\hat{Y},\hat{Z} が論理 Pauli 演算子となるように定義します。

 性質２（K の transversality）\hat{K}=K^{\otimes 5} は論理演算子です。これを示すには \hat{K} がスタビライザー群 G を群として保つこと、すなわち \hat{K} G \hat{K}^\dagger = G であることを確認すれば良いです。\hat{K} の共役作用は前節の K と同様で、

\begin{align*}
\hat{K} \hat{X} \hat{K}^{\dagger} &= \hat{Z} \\
\hat{K} \hat{Z} \hat{K}^{\dagger} &= \hat{Y} \\
\hat{K} \hat{Y} \hat{K}^{\dagger} &= \hat{X}
\end{align*}
が成立します。これは transversal な定義から明らかです。このため、 \hat{K} は位相を除いてスタビライザー空間に論理 K として作用することがわかります。

 性質３（射影と固有状態）\Pi \ket{K_0}^{\otimes 5} は \hat{K} の固有状態になります。具体的に、

\ket{K_1^L} = \sqrt{6} \Pi \ket{K_{00000}}, \quad \ket{K_0^L} = \sqrt{6} \Pi \ket{K_{11111}}
が成り立ちます。ここで \ket{K_{00000}} = \ket{K_0}^{\otimes 5}、\ket{K_{11111}} = \ket{K_1}^{\otimes 5} と定めます。論理基底 \ket{K_0^L},\ket{K_1^L} はスタビライザー空間上の正規化された \hat{K} 固有状態です (\hat{K}\ket{K_0^L}=e^{i\frac{\pi}{3}}\ket{K_0^L}, \hat{K}\ket{K_1^L}=e^{-i\frac{\pi}{3}}\ket{K_1^L})。証明のためにいくつか補題を示します。
!補題②： \ket{K_\mathbb{x}} を \Pi で射影したら \hat{K} の固有状態もしくはゼロベクトルとなる。

\hat{K} \Pi \ket{K_{\mathbb{x}}} = e^{i \frac{\pi}{3} (5 - 2 |\mathbb{x}|)} \Pi \ket{K_{\mathbb{x}}}
証明は \hat{K} と \Pi の可換性から直ちに従います。

これにより、|\mathbb{x}|=0,3 の場合は (もしゼロベクトルに射影されていなければ) \ket{K_1^L}、|\mathbb{x}|=2,5 の場合は (もしゼロベクトルに射影されていなければ) \ket{K_0^L} に射影され、|\mathbb{x}|=1,4 の場合はゼロベクトルとなることが分かります。
!補題③

\ket{K_{00000}}\bra{K_{00000}} = \frac{1}{2^5} \sum_{g \in S_+(5)} 3^{-\frac{|g|}{2}} g
証明は \ket{K_0}\bra{K_0} = \frac{1}{2} (I + \frac{1}{\sqrt{3}} (X + Y + Z)) を用いて展開すれば得られます。
!補題④
任意の g,h \in S_+(5) について次が成立します。

\text{Tr} (gh) = 2^5 \delta_{g,h}
パウリ文字列のうちトレースが非ゼロなのは恒等演算子のみであり、S_+(5) の各要素は二乗すると恒等演算子になることから示されます。
これらを用いて命題を証明します。

\begin{align*}
\bra{K_{00000}} \Pi \ket{K_{00000}} 
&= \text{Tr}(\Pi \ket{K_{00000}}\bra{K_{00000}}) \\
&= \frac{1}{2^9} \sum_{h \in G} \sum_{g \in S_+(5)} 3^{-\frac{|g|}{2}} \text{Tr} (gh) \\
&= \frac{1}{2^4} \sum_{h \in G} \sum_{g \in S_+(5)} 3^{-\frac{|g|}{2}} \delta_{g,h} \\
&= \frac{1}{2^4} \sum_{h \in G} 3^{-\frac{|h|}{2}} \\
&= \frac{1}{2^4} (1 + 15 \times 3^{-2}) = \frac{1}{6}
\end{align*}
(無事にゼロベクトルに射影されないことが示されました。)よって \ket{K_1^L} = \sqrt{6} \Pi \ket{K_{00000}} が得られ、同様に \Pi \ket{K_{11111}} についても計算できます。補足として、\Pi \ket{K_{\mathbb{x}}}（|\mathbb{x}|=2,3）の係数は \mathbb{x} に依存しますが、次の恒等式から総和が決まります。\Pi \ket{K_{\mathbb{x}}} = a_{\mathbb{x}} \ket{K_0^L} (|\mathbb{x}| = 2), \Pi \ket{K_{\mathbb{x}}} = b_{\mathbb{x}} \ket{K_1^L} (|\mathbb{x}| = 3) とします。

\begin{align*}
\ket{K_0^L} \bra{K_0^L} + \ket{K_1^L} \bra{K_1^L}
&= \Pi \\
&= \Pi \Pi \\
&= \Pi (\sum_{\mathbb{x}} \ket{K_{\mathbb{x}}} \bra{K_{\mathbb{x}}}) \Pi \\
&= (\frac{1}{6} + \sum_{|\mathbb{x}| = 2} |a_{\mathbb{x}}|^2) \ket{K_0^L} \bra{K_0^L} + (\frac{1}{6} + \sum_{|\mathbb{x}| = 3} |b_{\mathbb{x}}|^2)\ket{K_1^L} \bra{K_1^L}  \\
\end{align*}
従って \sum_{|\mathbb{x}| = 2} |a_{\mathbb{x}}|^2 = \sum_{|\mathbb{x}| = 3} |b_{\mathbb{x}}|^2 = \frac{5}{6} が得られます。

 魔法状態蒸留の手順魔法状態蒸留は、精度の低い魔法状態を複数用意し、スタビライザー測定によってスタビライザー空間に射影することで精度を高める手続きです。

具体的には次のようにして蒸留を行います。

 エラー解析精度の低い魔法状態を \rho = (1 - \epsilon) \ket{K_0}\bra{K_0} + \epsilon \ket{K_1}\bra{K_1} とすると、入力状態 \rho^{\otimes 5} は

\begin{align*}
\rho^{\otimes 5}
&= (1 - \epsilon)^5 \ket{K_{00000}}\bra{K_{00000}} \\
&+ (1 - \epsilon)^4 \epsilon \sum_{|\mathbb{x}| = 1} \ket{K_{\mathbb{x}}}\bra{K_{\mathbb{x}}} \\
&+ (1 - \epsilon)^3 \epsilon^2 \sum_{|\mathbb{x}| = 2} \ket{K_{\mathbb{x}}}\bra{K_{\mathbb{x}}} \\
&+ (1 - \epsilon)^2 \epsilon^3 \sum_{|\mathbb{x}| = 3} \ket{K_{\mathbb{x}}}\bra{K_{\mathbb{x}}} \\
&+ (1 - \epsilon) \epsilon^4 \sum_{|\mathbb{x}| = 4} \ket{K_{\mathbb{x}}}\bra{K_{\mathbb{x}}} \\
&+ \epsilon^5 \ket{K_{11111}}\bra{K_{11111}}
\end{align*}
となります。スタビライザー測定でエラーが検出されなかった場合、測定後の状態は \Pi \rho^{\otimes 5} \Pi の定数倍（トレースが 1 になるよう正規化された状態）です。

\Pi \rho^{\otimes 5} \Pi を計算してみましょう。

\begin{align*}
\Pi \rho^{\otimes 5} \Pi
&= \Pi (1 - \epsilon)^5 \ket{K_{00000}}\bra{K_{00000}} \Pi \\
&+ \Pi (1 - \epsilon)^4 \epsilon \sum_{|\mathbb{x}| = 1} \ket{K_{\mathbb{x}}}\bra{K_{\mathbb{x}}} \Pi \\
&+ \Pi (1 - \epsilon)^3 \epsilon^2 \sum_{|\mathbb{x}| = 2} \ket{K_{\mathbb{x}}}\bra{K_{\mathbb{x}}} \Pi \\
&+ \Pi (1 - \epsilon)^2 \epsilon^3 \sum_{|\mathbb{x}| = 3} \ket{K_{\mathbb{x}}}\bra{K_{\mathbb{x}}} \Pi \\
&+ \Pi (1 - \epsilon) \epsilon^4 \sum_{|\mathbb{x}| = 4} \ket{K_{\mathbb{x}}}\bra{K_{\mathbb{x}}} \Pi \\
&+ \Pi \epsilon^5 \ket{K_{11111}}\bra{K_{11111}} \Pi \\
&= \frac{1}{6} (1 - \epsilon)^5 \ket{K_1^L}\bra{K_1^L}
+ \sum_{|\mathbb{x}| = 2} |a_{\mathbb{x}}|^2 (1 - \epsilon)^3 \epsilon^2 \ket{K_0^L}\bra{K_0^L} \\
&+ \sum_{|\mathbb{x}| = 3} |b_{\mathbb{x}}|^2 (1 - \epsilon)^2 \epsilon^3 \ket{K_1^L}\bra{K_1^L}
+ \frac{1}{6} \epsilon^5 \ket{K_0^L}\bra{K_0^L} \\
&= \frac{1}{6} (5 (1 - \epsilon)^3 \epsilon^2 + \epsilon^5) \ket{K_0^L}\bra{K_0^L} + \frac{1}{6} ((1 - \epsilon)^5 + 5 (1 - \epsilon)^2 \epsilon^3) \ket{K_1^L}\bra{K_1^L}
\end{align*}
このとき \ket{K_0^L}\bra{K_0^L} と \ket{K_1^L}\bra{K_1^L} の係数をそれぞれ A,B とおくと

\begin{align*}
A &= \frac{1}{6}(5 (1 - \epsilon)^3 \epsilon^2 + \epsilon^5) \\
B &= \frac{1}{6}((1 - \epsilon)^5 + 5 (1 - \epsilon)^2 \epsilon^3)
\end{align*}
となります。出力誤り率を \epsilon_{out} = A/(A+B) と定義し、t=\frac{\epsilon}{1-\epsilon} を導入すると

\epsilon_{out} = \frac{t^5+5t^2}{t^5+5t^3+5t^2+1}
が得られます。したがって \epsilon \ll 1 では \epsilon_{out} \approx 5 \epsilon^2 が成り立ち、「誤りが二次で抑制される」ことが見て取れます。
!この手続きで高い精度で得られる状態は \ket{K_1^L} ですが、デコードした \ket{K_1} に Clifford 操作 Y H を行えば \ket{K_0} が得られます。このため、魔法状態蒸留で高い精度で得られる状態が \ket{K_0^L} か \ket{K_1^L} かはどちらでも問題ありません。
成功確率は p_{success} = \text{Tr} (\Pi \rho^{\otimes 5} \Pi) = A+B であり、

p_{success} = \frac{1}{6} \left((1 - \epsilon)^5 + 5 (1 - \epsilon)^2 \epsilon^2 + \epsilon^5\right)
と書けます。
魔法状態蒸留の閾値は \epsilon_{out}=\epsilon を解いて \epsilon_{th} = \frac{1}{2} (1 - \sqrt{\frac{3}{7}}) \approx 0.173 で、これより小さい誤り率であれば \epsilon_{out} < \epsilon が成立します。
\epsilon_{out} と p_{success} をプロットすると次のようになります。

 具体的な回路魔法状態蒸留を回路として実装してみましょう。ここでは Qiskit の Clifford クラスを使って perfect code の encode/decode 回路を実装します。そのためにまず destabilizer を見つけます。実装は次のような流れで進めます。
destabilizer を見つける
encode/decode 回路を実装する
魔法状態蒸留を検証する

 destabilizer を見つけるここでは stabilizer を拡張し S_5=ZZZZZ とします。destabilizer D_1,D_2,D_3,D_4,D_5 は次の条件を満たすものです。
D_5=XXXXX

D_1,D_2,D_3,D_4,D_5 は可換

D_i は S_i と反可換で、 S_j (j \neq i) と可換
これを満たすパウリ文字列を全探索します。
destabilizer の探索コードdestabilizer.py
import itertools

from qiskit.quantum_info import Pauli

S = ["XZZXI", "IXZZX", "XIXZZ", "ZXIXZ"]
LOGICAL_Z = "ZZZZZ"
LOGICAL_X = "XXXXX"


def _commutes_all(labels):
    objs = [Pauli(s) for s in labels]
    for i in range(len(objs)):
        for j in range(i + 1, len(objs)):
            if not objs[i].commutes(objs[j]):
                return False
    return True


def _anticommutes_only_with(target, stabilizers, index):
    for j, stab in enumerate(stabilizers):
        if j == index:
            if not target.anticommutes(stab):
                return False
        else:
            if not target.commutes(stab):
                return False
    return True


def _find_good_destabilizer(destabilizers_candidates):
    rets = []
    for cand in destabilizers_candidates:
        # すべての destabilizer が XXZZI の順番を入れ替えたものであるものを探す
        x2z2i1 = True
        for pauli in cand:
            if not (
                pauli.count("X") == 2
                and pauli.count("Z") == 2
                and pauli.count("I") == 1
            ):
                x2z2i1 = False
                break
        if x2z2i1:
            rets.append(cand)
    return rets


def find_destabilizers():
    paulis = ["I", "X", "Y", "Z"]
    stabilizers = [Pauli(s) for s in S] + [Pauli(LOGICAL_Z)]
    logical_x = Pauli(LOGICAL_X)

    # destabilizer の候補を探す。
    # candidates_by_idx[i] は s[i] と非可換で、ほかのスタビライザーと
    # XXXXX とは可換なパウリ文字列のリストとなるよう計算する。
    candidates_by_idx = {i: [] for i in range(4)}
    for p in itertools.product(paulis, repeat=5):
        label = "".join(p)
        if label == "IIIII":
            continue
        pobj = Pauli(label)
        if not pobj.commutes(logical_x):
            continue
        for idx in range(4):
            if _anticommutes_only_with(pobj, stabilizers, idx):
                candidates_by_idx[idx].append(label)

    # destabilizer のすべての可能性を列挙し、 destabilizer が可換となる組み合わせを見つける
    destabilizers_candidates = []
    for d0 in candidates_by_idx[0]:
        for d1 in candidates_by_idx[1]:
            for d2 in candidates_by_idx[2]:
                for d3 in candidates_by_idx[3]:
                    cand = [d0, d1, d2, d3]
                    if _commutes_all(cand + [LOGICAL_X]):
                        destabilizers_candidates.append(cand)
                        break

    if len(destabilizers_candidates) == 0:
        raise RuntimeError("No commuting destabilizer set found.")

    return destabilizers_candidates


ds = find_destabilizers()
goods = _find_good_destabilizer(ds)

print("stabilizers:", S + [LOGICAL_Z])
print("Number of destabilizers candidates:", len(ds))
print("Number of good destabilizers:", len(goods))
print("destabilizers:", goods[0] + [LOGICAL_X])
実行結果は次のようになります。 destabilizer は 512 通りも選び方がありますが、今回はその中でも D_1=ZXIZX, D_2=ZIXXZ, D_3=IZZXX, D_4=XZXZI を採用します。
stabilizers: ['XZZXI', 'IXZZX', 'XIXZZ', 'ZXIXZ', 'ZZZZZ']
Number of destabilizers candidates: 512
Number of good destabilizers: 1
destabilizers: ['ZXIZX', 'ZIXXZ', 'IZZXX', 'XZXZI', 'XXXXX']

 encode/decode 回路を実装するQiskit の Clifford クラスを使って encode/decode 回路を実装します。decode 回路は encode の逆として実装できます。また、 encode 回路について以下を検証します。
encode 前の補助量子ビットへの Z 演算子が、 encode 後にスタビライザー S_1,S_2,S_3,S_4 に map される
encode 前の量子ビットへの X,Y,Z 演算子が、 encode 後に論理演算子 XXXXX,YYYYY,ZZZZZ に map される
encode/decode 回路の実装codec.py
from qiskit.quantum_info import Pauli, Clifford

S = ["XZZXI", "IXZZX", "XIXZZ", "ZXIXZ"]
D = ["ZXIZX", "ZIXXZ", "IZZXX", "XZXZI"]
LOGICAL_Z = "ZZZZZ"
LOGICAL_X = "XXXXX"


def build_encoder_decoder():
    stabilizers = [
        "+" + LOGICAL_Z,
        "+" + S[3],
        "+" + S[2],
        "+" + S[1],
        "+" + S[0],
    ]
    destabilizers = [
        "+" + LOGICAL_X,
        "+" + D[3],
        "+" + D[2],
        "+" + D[1],
        "+" + D[0],
    ]
    encoder = Clifford.from_dict(
        {"stabilizer": stabilizers, "destabilizer": destabilizers}
    ).to_circuit()
    decoder = encoder.inverse()
    return encoder, decoder


encoder, decoder = build_encoder_decoder()

print("encoder circuit:")
print(encoder)

print("decoder circuit:")
print(decoder)

##########################################
# encoder の検証
##########################################

mapping = {}
ancilla_z = ["ZIIII", "IZIII", "IIZII", "IIIZI"]
for i, lab in enumerate(ancilla_z):
    out = Pauli(lab).evolve(encoder, frame="s").to_label()
    mapping[f"Z on ancilla q{i}"] = out
mapping["logical Z"] = Pauli("IIIIZ").evolve(encoder, frame="s").to_label()
mapping["logical X"] = Pauli("IIIIX").evolve(encoder, frame="s").to_label()
mapping["logical Y"] = Pauli("IIIIY").evolve(encoder, frame="s").to_label()

for key, value in mapping.items():
    print(f"{key} -> {value}")
実行結果は次のようになります。 encode/decode 回路の取得と encode 回路の検証ができました。
encoder circuit:
          ┌───┐     ┌───┐     ┌───┐                                                                                           
q_0: ──■──┤ X ├─────┤ X ├──■──┤ Y ├───────────────────────────────────────────────────────────────────────────────────────────
     ┌─┴─┐└─┬─┘     └─┬─┘  │  └───┘                                      ┌───┐┌───┐     ┌───┐┌───┐┌───┐        ┌───┐┌───┐┌───┐
q_1: ┤ X ├──┼─────────┼────┼────■────────────────────────────────────────┤ X ├┤ H ├──■──┤ S ├┤ H ├┤ S ├─X───■──┤ H ├┤ S ├┤ Y ├
     └───┘  │         │  ┌─┴─┐  │  ┌───┐     ┌───┐┌───┐             ┌───┐└─┬─┘├───┤  │  └───┘└───┘└───┘ │ ┌─┴─┐├───┤└───┘└───┘
q_2: ───────┼────■────┼──┤ X ├──┼──┤ H ├──■──┤ H ├┤ S ├───────────X─┤ X ├──┼──┤ H ├──┼──────────────────X─┤ X ├┤ Y ├──────────
     ┌───┐  │  ┌─┴─┐  │  └───┘  │  └───┘  │  ├───┤├───┤     ┌───┐ │ └─┬─┘  │  └───┘┌─┴─┐┌───┐             └───┘└───┘          
q_3: ┤ H ├──┼──┤ X ├──■─────────┼─────────┼──┤ X ├┤ H ├──■──┤ S ├─X───■────■───────┤ X ├┤ Y ├─────────────────────────────────
     └───┘  │  ├───┤          ┌─┴─┐     ┌─┴─┐└─┬─┘└───┘┌─┴─┐├───┤                  └───┘└───┘                                 
q_4: ───────■──┤ H ├──────────┤ X ├─────┤ X ├──■───────┤ X ├┤ Y ├─────────────────────────────────────────────────────────────
               └───┘          └───┘     └───┘          └───┘└───┘                                                             
decoder circuit:
     ┌───┐                                                                                                        ┌───┐     ┌───┐     
q_0: ┤ Y ├─────────────────────────────────────────────────────────────────────────────────────────────────────■──┤ X ├─────┤ X ├──■──
     ├───┤┌─────┐┌───┐        ┌─────┐┌───┐┌─────┐     ┌───┐┌───┐                                               │  └─┬─┘     └─┬─┘┌─┴─┐
q_1: ┤ Y ├┤ Sdg ├┤ H ├──■───X─┤ Sdg ├┤ H ├┤ Sdg ├──■──┤ H ├┤ X ├─────────────────────────────────────■─────────┼────┼─────────┼──┤ X ├
     ├───┤└─────┘└───┘┌─┴─┐ │ └┬───┬┘└───┘└─────┘  │  └───┘└─┬─┘┌───┐   ┌─────┐┌───┐                 │  ┌───┐┌─┴─┐  │         │  └───┘
q_2: ┤ Y ├────────────┤ X ├─X──┤ H ├───────────────┼─────────┼──┤ X ├─X─┤ Sdg ├┤ H ├────────────■────┼──┤ H ├┤ X ├──┼────■────┼───────
     ├───┤            └───┘    └───┘             ┌─┴─┐       │  └─┬─┘ │ ├─────┤└───┘┌───┐┌───┐  │    │  └───┘└───┘  │  ┌─┴─┐  │  ┌───┐
q_3: ┤ Y ├───────────────────────────────────────┤ X ├───────■────■───X─┤ Sdg ├──■──┤ H ├┤ X ├──┼────┼──────────────■──┤ X ├──┼──┤ H ├
     ├───┤                                       └───┘                  └─────┘┌─┴─┐└───┘└─┬─┘┌─┴─┐┌─┴─┐┌───┐          └───┘  │  └───┘
q_4: ┤ Y ├─────────────────────────────────────────────────────────────────────┤ X ├───────■──┤ X ├┤ X ├┤ H ├─────────────────■───────
     └───┘                                                                     └───┘          └───┘└───┘└───┘                         
Z on ancilla q0 -> XZZXI
Z on ancilla q1 -> IXZZX
Z on ancilla q2 -> XIXZZ
Z on ancilla q3 -> ZXIXZ
logical X -> XXXXX
logical Y -> YYYYY
logical Z -> ZZZZZ

 魔法状態蒸留を検証する魔法状態蒸留の手続きを愚直に再現して \epsilon_{out}, p_{success} をシミュレートします。
魔法状態蒸留の実装まずは入力状態を作成するヘルパー関数 make_input を定義します。
import cmath
import math

import numpy as np


def k_projectors():
    z = 1 / math.sqrt(3)
    cos_t = math.sqrt((1 + z) / 2)
    sin_t = math.sqrt((1 - z) / 2)
    phase = cmath.exp(1j * math.pi / 4)

    k0 = np.array([cos_t, phase * sin_t], dtype=complex)
    k1 = np.array([sin_t, -phase * cos_t], dtype=complex)

    p0 = np.outer(k0, k0.conj())
    p1 = np.outer(k1, k1.conj())
    return p0, p1


def kron_n(mat, n):
    out = mat
    for _ in range(n - 1):
        out = np.kron(out, mat)
    return out


def make_input(eps):
    p0, p1 = k_projectors()
    rho = (1 - eps) * p0 + eps * p1
    return kron_n(rho, 5)
蒸留手続きを distill() 関数で実装します。 \epsilon と decoder 回路を受け取って、 \epsilon_{out}, p_{success} を計算する関数です。実装のポイントは次の 3 点です。
decoder 回路の適用 : Qiskit の Operator クラスを使って Clifford クラスに対応する行列を作成します。
射影演算子の適用: decoder 回路を適用した後にスタビライザー空間に射影します。 decoder 回路はスタビライザーを補助量子ビットの Z に map するため、+1 固有空間への射影演算子は非常に単純な形をしています。

\epsilon_{out} の計算: \ket{K_0}\bra{K_0}, \ket{K_1}\bra{K_1} の係数を求めるために、射影演算子を適用してからトレースを計算しています。
from qiskit.quantum_info import Operator, DensityMatrix, partial_trace


def distill(eps, decoder):
    # 入力状態 \rho^{\otimes 5} の作成
    rho5 = make_input(eps)

    # decode 回路の適用
    u = Operator(decoder).data
    rho_dec = u @ rho5 @ u.conj().T

    # 射影演算子の適用
    # decode 回路の適用後なので、論理演算子におけるスタビライザーは補助量子ビットの Z に map されている
    # このため、射影演算子は下の proj のように実装できる
    proj_anc = np.zeros((16, 16), dtype=complex)
    proj_anc[0, 0] = 1.0
    proj = np.kron(proj_anc, np.eye(2))
    rho_post = proj @ rho_dec @ proj

    # 魔法状態蒸留の成功確率の計算
    p_success = np.trace(rho_post).real

    # \epsilon_{out} の計算
    rho_post = rho_post / p_success
    dm = DensityMatrix(rho_post)
    rho_data = partial_trace(dm, [1, 2, 3, 4]).data

    p0, p1 = k_projectors()
    e0 = np.real(np.trace(p0 @ rho_data))
    e1 = np.real(np.trace(p1 @ rho_data))
    eps_out = e0 / (e0 + e1)
    return eps_out, p_success
これまでに実装した decoder を使って distill() 関数を様々な \epsilon で実行し、得られた \epsilon_{out},p_{success} をプロットしたら次のようになり、理論値と完全に一致することが確かめられます。

 まとめperfect code では \hat{K}=K^{\otimes 5} が論理演算子になるため、T 型魔法状態の蒸留に適しています。出力誤り率は \epsilon_{out}\approx 5\epsilon^2 と二次で抑制されます。一方で成功確率は p_s\sim 1/6 と低く、蒸留コストの主因になります。本稿の内容についてきちんと学ばれたい方は [2:1] や [5] を読むことをお勧めします。
魔法状態蒸留に関する素晴らしい理論である Bravyi-Haah [6] や、中性原子型量子コンピュータでの魔法状態蒸留の実証実験 [7] 等についても今後紹介できればと考えています。

脚注
Gidney, Craig, Noah Shutty, and Cody Jones. "Magic state cultivation: growing T states as cheap as CNOT gates." arXiv preprint arXiv:2409.17595 (2024). ↩︎
Bravyi, Sergey, and Alexei Kitaev. "Universal quantum computation with ideal Clifford gates and noisy ancillas." Physical Review A—Atomic, Molecular, and Optical Physics 71.2 (2005): 022316. ↩︎ ↩︎
https://pennylane.ai/qml/glossary/what-are-magic-states ↩︎
https://errorcorrectionzoo.org/c/stab_5_1_3 ↩︎
https://www.iqst.ca/events/csqic05/talks/nathan b.pdf ↩︎
Bravyi, Sergey, and Jeongwan Haah. "Magic-state distillation with low overhead." Physical Review A—Atomic, Molecular, and Optical Physics 86.5 (2012): 052329. ↩︎
Sales Rodriguez, Pedro, et al. "Experimental demonstration of logical magic state distillation." Nature 645.8081 (2025): 620-625. ↩︎
今更学ぶ魔法状態蒸留

この記事の目的

背景

注意点

記号と基本事項

魔法状態の定義

[[5,1,3]] perfect code

定義

性質１（パウリの transversality）

性質２（ $K$ の transversality）

性質３（射影と固有状態）

魔法状態蒸留の手順

エラー解析

具体的な回路

destabilizer を見つける

encode/decode 回路を実装する

魔法状態蒸留を検証する

まとめ

Discussion
