スワップテスト

アダマールゲートに続いてスワップテストです。 これはアダマールテストのユニタリーゲートについて、SWAPゲートを用いる場合です。 大雑把には制御量子ビット・入力量子ビット1(標的量子ビット1)・入力量子ビット2(標的量子ビット2)の3種類の役割を持つ量子ビットが存在します。 (入力量子ビットは1量子ビットとは限らず、複数量子ビットからなる場合もあります。) 入力量子ビット1, 2を2つのベクトルとした場合の2ベクトル間の内積を評価する量子回路です。

内積を評価する場面はさまざまあり、実装も簡単なため量子機械学習ほか様々な場面で使われるようです。

前提とする知識

• 線形代数がわかること
• 簡単な論理演算(AND, OR, NOT)が分かること
• 基本的な量子ゲート(X, CX, CCX)が分かること

バージョン情報

• Python 3.9.13
• Qiskit 0.36.2

目次

SWAPテスト回路

アダマールテストおさらい

アダマールテストの詳細については前回の解説を参照してください。

量子回路の形は下図の通りで、この回路を使うことで、ユニタリーU^$\hat{U}$によって第一量子ビットを測定した時の状態の確率分布を変化させられることがSWAPテストでは重要な意味があります。

図1. アダマールテストの量子回路

この時の第一量子ビットの|0⟩$|0⟩$及び|1⟩$|1⟩$状態の確率分布は以下です。

p0=1+Re⟨ψ|BU^|ψ⟩B2${\displaystyle p_0=\frac{1+\mathrm{R}\mathrm{e}⟨\psi {|}_B\hat{U}|\psi {⟩}_B}{2}}$

p1=1−Re⟨ψ|BU^|ψ⟩B2${\displaystyle p_1=\frac{1-\mathrm{R}\mathrm{e}⟨\psi {|}_B\hat{U}|\psi {⟩}_B}{2}}$

SWAPテスト回路の理論

SWAPテストによって二つのベクトルの内積を計算できます。 内積の値がアダマールテストの第一量子ビットの確率分布に反映されるのでその確率から内積を逆算できる、ということです。

入力する2つの量子状態を以下のように定義します。

|a⟩=∑i=0N−1ai|i⟩${\displaystyle |a⟩=\sum _{i=0}^{N-1}{a}_i|i⟩}$

|b⟩=∑i=0N−1bi|i⟩${\displaystyle |b⟩=\sum _{i=0}^{N-1}{b}_i|i⟩}$

この2つのベクトル|a⟩,|b⟩$|a⟩,|b⟩$を以下の量子回路に入力することで、その内積

|⟨a|b⟩|2=|∑i=0N−1aibi|2${\displaystyle |⟨a|b⟩{|}^2=|\sum _{i=0}^{N-1}{a}_i{b}_i{|}^2}$

を評価できます。

図2. SWAPテスト回路

第一量子ビットの状態を測定した時に、|0⟩$|0⟩$となる確率が以下となることから内積を評価できます。

p0=1+|⟨a|b⟩|22${\displaystyle p_0=\frac{1+|⟨a|b⟩{|}^2}{2}}$

SWAPテストの仕組み導出

まずアダマールゲートを実行すると

|0⟩|a⟩|b⟩⟶H112–√(|0⟩|a⟩|b⟩+|1⟩|a⟩|b⟩)${\displaystyle |0⟩|a⟩|b⟩\stackrel{H_1}{⟶}\frac{1}{\sqrt{2}}(|0⟩|a⟩|b⟩+|1⟩|a⟩|b⟩)}$

そして前回解説した制御SWAPゲートにより入力ビットの状態を反転させると以下になります。

12–√(|0⟩|a⟩|b⟩+|1⟩|b⟩|a⟩)${\displaystyle \frac{1}{\sqrt{2}}(|0⟩|a⟩|b⟩+|1⟩|b⟩|a⟩)}$

そして最後に第一量子ビットにアダマールゲートをかけると全体の量子状態の最終状態|ψ⟩$|\psi ⟩$は以下になります。

|ψ⟩=12{(|0⟩+|1⟩)|a⟩|b⟩+(|0⟩−|1⟩)|b⟩|a⟩}=12{|0⟩(|a⟩|b⟩+|b⟩|a⟩)+|1⟩(|a⟩|b⟩−|b⟩|a⟩)}$\begin{array}{rl}|\psi ⟩& =\frac{1}{2}\{(|0⟩+|1⟩)|a⟩|b⟩+(|0⟩-|1⟩)|b⟩|a⟩\}\\ & =\frac{1}{2}\{|0⟩(|a⟩|b⟩+|b⟩|a⟩)+|1⟩(|a⟩|b⟩-|b⟩|a⟩)\}\end{array}$

補助量子ビットの状態が|0⟩$|0⟩$となる確率はベクトル|ψ⟩$|\psi ⟩$の|0⟩1$|0{⟩}_1$成分の長さを求めれば良いので以下となります。

p0=|⟨0|ψ⟩|2=(⟨0|ψ⟩)†⟨0|ψ⟩=⟨ψ|0⟩⟨0|ψ⟩=14(⟨a|⟨b|+⟨b|⟨a|)(|a⟩|b⟩+|b⟩|a⟩)=14(⟨a|⟨b|a⟩|b⟩+⟨a|⟨b|b⟩|a⟩+⟨b|⟨a|a⟩|b⟩+⟨b|⟨a|b⟩|a⟩)$\begin{array}{rl}{p}_0& =|⟨0|\psi ⟩{|}^2\\ & =(⟨0|\psi ⟩{)}^{†}⟨0|\psi ⟩\\ & =⟨\psi |0⟩⟨0|\psi ⟩\\ & =\frac{1}{4}(⟨a|⟨b|+⟨b|⟨a|)(|a⟩|b⟩+|b⟩|a⟩)\\ & =\frac{1}{4}(⟨a|⟨b|a⟩|b⟩+⟨a|⟨b|b⟩|a⟩+⟨b|⟨a|a⟩|b⟩+⟨b|⟨a|b⟩|a⟩)\end{array}$

⟨a|b⟩$⟨a|b⟩$などは内積取ってしまったらただの数なので自由に交換できますので式を整理し

p0=14(⟨a|b⟩⟨b|a⟩+⟨a|a⟩⟨b|b⟩+⟨b|b⟩⟨a|a⟩+⟨b|a⟩⟨a|b⟩)=14(|⟨a|b⟩|2+1+1+|⟨a|b⟩|2)=14(2+2|⟨a|b⟩|2)=1+|⟨a|b⟩|22$\begin{array}{rl}{p}_0& =\frac{1}{4}(⟨a|b⟩⟨b|a⟩+⟨a|a⟩⟨b|b⟩+⟨b|b⟩⟨a|a⟩+⟨b|a⟩⟨a|b⟩)\\ & =\frac{1}{4}(|⟨a|b⟩{|}^2+1+1+|⟨a|b⟩{|}^2)\\ & =\frac{1}{4}(2+2|⟨a|b⟩{|}^2)\\ & =\frac{1+|⟨a|b⟩{|}^2}{2}\end{array}$

以上では愚直に計算しましたが、制御ユニタリで実行しているゲートU^$\hat{U}$に対応するものとしてSWAPゲートを示すU^SWAP${\hat{U}}_{\mathrm{S}\mathrm{W}\mathrm{A}\mathrm{P}}$にして実行したアダマールテストです。

アダマールテストだと思った時の入力量子ビットは|ψ⟩=|a⟩|b⟩$|\psi ⟩=|a⟩|b⟩$です。

SWAPゲートを示すユニタリーU^SWAP${\hat{U}}_{\mathrm{S}\mathrm{W}\mathrm{A}\mathrm{P}}$は以下です。

U^SWAP=⎛⎝⎜⎜⎜1000001001000001⎞⎠⎟⎟⎟${\displaystyle {\hat{U}}_{\mathrm{S}\mathrm{W}\mathrm{A}\mathrm{P}}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right)}$

この行列の固有ベクトルは

|x+⟩=⎛⎝⎜⎜⎜0110⎞⎠⎟⎟⎟=|1⟩|0⟩+|0⟩|1⟩$\begin{array}{rl}|x_{+}⟩& =\left(\begin{array}{c}0\\ 1\\ 1\\ 0\end{array}\right)\\ & =|1⟩|0⟩+|0⟩|1⟩\end{array}$

|x−⟩=⎛⎝⎜⎜⎜01−10⎞⎠⎟⎟⎟=|1⟩|0⟩−|0⟩|1⟩$\begin{array}{rl}|x_{-}⟩& =\left(\begin{array}{c}0\\ 1\\ -1\\ 0\end{array}\right)\\ & =|1⟩|0⟩-|0⟩|1⟩\end{array}$

のため、一般の入力状態

|ψ⟩=|a⟩|b⟩=(α|0⟩+β|1⟩)(u|0⟩+v|1⟩)=αu|0⟩|0⟩+αv|0⟩|1⟩+βu|1⟩|0⟩+βv|1⟩|1⟩=⎛⎝⎜⎜⎜αuαvβuβv⎞⎠⎟⎟⎟$\begin{array}{rl}|\psi ⟩& =|a⟩|b⟩\\ & =(\alpha |0⟩+\beta |1⟩)(u|0⟩+v|1⟩)\\ & =\alpha u|0⟩|0⟩+\alpha v|0⟩|1⟩+\beta u|1⟩|0⟩+\beta v|1⟩|1⟩\\ & =\left(\begin{array}{c}\alpha u\\ \alpha v\\ \beta u\\ \beta v\end{array}\right)\end{array}$

は固有状態とは限りません。

従って第一量子ビットの状態が|0⟩$|0⟩$である確率は以下となります。

p0=1+Re⟨ψ|U^SWAP|ψ⟩2=1+α2u2+2αβuv+β2v22=1+(αu+βv)22=1+∣∣∣(αβ)(uv)∣∣∣22=1+|⟨a|b⟩|22$\begin{array}{rl}{p}_0& =\frac{1+\mathrm{R}\mathrm{e}⟨\psi |{\hat{U}}_{\mathrm{S}\mathrm{W}\mathrm{A}\mathrm{P}}|\psi ⟩}{2}\\ & =\frac{1+{\alpha }^2{u}^2+2\alpha \beta uv+{\beta }^2{v}^2}{2}\\ & =\frac{1+(\alpha u+\beta v{)}^2}{2}\\ & =\frac{1+{\left|\left(\begin{array}{cc}\alpha & \beta \end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)\right|}^2}{2}\\ & =\frac{1+{|⟨a|b⟩|}^2}{2}\end{array}$

ということで、やはり愚直にやった計算と一致しました。
そのためユニタリーがわかっていれば、今後アダマールテストの計算する時は愚直に計算しなくとも、確率の式に当てはめれば良いです。

この確率の式から、入力ベクトルが直交している場合は確率0.5で平行であれば確率1であることがわかる。(ベクトルの正負の向きまではわからない。)

SWAPテスト回路の実装

実際に内積を計算してみます。

単純なベクトルの計算

今は|a⟩=α|0⟩+β|1⟩,|b⟩=u|0⟩+v|1⟩$|a⟩=\alpha |0⟩+\beta |1⟩,|b⟩=u|0⟩+v|1⟩$の二次元平面上のベクトルを考えています。
普段量子ビットを考えるときは[tex | 0 \rangle]と|1⟩$|1⟩$はBlock球面上の上下正反対を向いたベクトル(下図左)ですが、この二つのベクトルは正規直交基底(直交して、その平面・空間の軸に平行な単位ベクトル)なので直交座標上では直交しています(下図右)。

図3. | 0 >と| 1 >の位置関係 左: Bloch球面上 右: 直交座標上

わかりやすい例だと|a⟩=|0⟩,|b⟩=|0⟩$|a⟩=|0⟩,|b⟩=|0⟩$なら⟨a|b⟩=1$⟨a|b⟩=1$で |a⟩=|0⟩,|b⟩=|1⟩$|a⟩=|0⟩,|b⟩=|1⟩$なら⟨a|b⟩=0$⟨a|b⟩=0$です。

|a⟩$|a⟩$は|0⟩$|0⟩$に固定して、|b⟩$|b⟩$を|0⟩$|0⟩$から−|0⟩$-|0⟩$まで直交座標上でπ$\pi$ (Bloch球面上で2π$2\pi$)だけ連続的に変化させた時は

⟨a|b⟩=⟨0|(cosθ2|0⟩+sinθ2|1⟩)=cosθ2$\begin{array}{rl}⟨a|b⟩& =⟨0|(\cos \frac{\theta }{2}|0⟩+\sin \frac{\theta }{2}|1⟩)\\ & =\cos \frac{\theta }{2}\end{array}$

のため下図のようになります。(θ$\theta$はBloch球面上の|0⟩$|0⟩$となす角度)

図4. <a | b >の角度依存性

|0⟩$|0⟩$からπ$\pi$だけ角度が違うものが|1⟩$|1⟩$に一致して、確かに直交してcosπ/2=0$\cos \pi /2=0$になってますね。

スワップテスト回路で内積を計算

単純にベクトル計算した時と同じように、以下のように|a⟩$|a⟩$は固定して|b⟩$|b⟩$をY軸に周りに|0⟩$|0⟩$から|1⟩$|1⟩$を経由して一周して|0⟩$|0⟩$に戻るまで角度を変えていった場合に内積が計算できていることを確認します。

図5. スワップテスト回路の例

これで|b⟩$|b⟩$の角度を0から2π$2\pi$に変えた時に第一量子ビットで|0⟩$|0⟩$が観測される確率が下図です。

図6. 第一量子ビットで| 0 >が観測される確率

p0=1+|⟨a|b⟩|22${\displaystyle p_0=\frac{1+{|⟨a|b⟩|}^2}{2}}$

から|⟨a|b⟩|2${|⟨a|b⟩|}^2$を算出すると下図です。
ノイズの関係でマイナスになる関係で残念ながら自乗を外せません。

図7. |<a|b>|^2を測定と計算から求めた場合の比較

ベクトル計算から求めた場合と測定から求めた場合とでの誤差をRMSEで求めると0.02程度でした。
目盛の間隔が0.2ずつなのでその1/10です。

この誤差が大きいのか小さいのかは実際に解きたい問題に依存すると思うのでここでは評価しませんが、ほぼほぼベクトル計算だけから求めた場合と同等の値を測定から求められると思って良いと思います。

この時のソースです。

effs = np.linspace(0, 2, 51)
print(effs)
probs = []
for idx, eff in enumerate(effs):
    rad = np.pi * eff
    
    qr_control = QuantumRegister(1, 'control')
    qr_1 = QuantumRegister(1, 'a')
    qr_2 = QuantumRegister(1, 'b')
    cr = ClassicalRegister(1)
    qc = QuantumCircuit(qr_control, qr_1, qr_2, cr)

    qc.rx(rad, qr_2)

    qc.barrier()

    qc.h(qr_control)

    qc.cx(qr_1, qr_2)
    qc.ccx(qr_control, qr_2, qr_1)
    qc.cx(qr_1, qr_2)

    qc.h(qr_control)

    qc.barrier()

    qc.measure(qr_control, cr)
        
    backend = Aer.get_backend('qasm_simulator')
    job = execute(qc, backend, shot=1024)

    result = job.result()
    probs.append({
        'theta': rad,
        'prob': calc_p0(result)
    })

df_result = pd.DataFrame(probs)
display(df_result)

thetas = np.linspace(0, 2, 51)
thetas = np.asarray([theta * np.pi for theta in thetas])
cos = np.cos(thetas / 2)

fig = plt.figure()
plt.scatter(df_result['theta'] / np.pi, (df_result['prob']*2-1), marker='o', label='測定結果から求めた結果', color='blue')
plt.plot(thetas / np.pi, cos**2, label='代数的に求めた結果', color='red', linewidth=2)
rmse = sqrt(mean_squared_error(df_result['prob']*2-1, cos**2))
xlim = plt.xlim()
ylim = plt.ylim()
dy = ylim[1] - ylim[0]
plt.text(x=xlim[0], y=ylim[0]+0.1*dy, s=f'RMSE: {rmse:.3f}', fontsize=12)
plt.xlabel('θ / π', fontsize=15)
plt.ylabel('|<a|b>|^2', fontsize=15)
plt.legend(fontsize=12)
plt.show()
plt.close()

光の系でスワップテストを実装(ちゃんと読んでないですけど実験の提案だけかも)の論文がありました。

www.nature.com

まとめ

• スワップテストは2つのベクトル(2つの状態)の内積を評価できる
• スワップテストはアダマールテストのユニタリーがスワップゲートの場合に該当する

参考

• Swap test と quantum neuron について - 物理とか
• 量子計算のスワップテスト（swap test）を試してみる - Qiita

量子コンピューティング 基本アルゴリズムから量子機械学習まで [ 情報処理学会出版委員会 ]

SWAPゲート・制御SWAPゲート

先日アダマールゲートを解説したので続いてSWAPテストを解説したいのですが、その前にSWAPゲートと制御SWAPゲートを解説します。

対象とする読者

• 量子プログラミングに興味を持ち始めたばかりのエンジニア
• 量子計算に興味を持ち始めたばかりの高校生・大学生
• 数学・論理演算・電子回路に慣れていない方
  など

前提とする知識

• 簡単な論理演算(AND, OR, NOT)が分かること
• 基本的な量子ゲート(X, CX, CCX)が分かること

バージョン情報

• Python 3.9.13
• Qiskit 0.36.2

目次

SWAPゲート

SWAPゲート実装

英単語「swap」の意味は「交換する」という動詞です。

eow.alc.co.jp

2つの量子ビットの値を交換するというゲートです。

下図のようにCXゲートを3つ組み合わせます。

図1. SWAPゲート

この時の入出力の対応表は以下になっております。

入力		出力
量子ビット1(q1)	量子ビット2(q2)	量子ビット1(q1)	量子ビット2(q2)
0	0	0	0
0	1	1	0
1	0	0	1
1	1	1	1

それぞれの入力を用いて出力を測定した結果が下図です。

図2. SWAPゲートの実行結果(分かりにくいですが、グラフ中のビット列は右から第一量子ビットです。そのため、(1)の出力はそのまま|00>, (2)の出力は右から読み直して|10>, (3)の出力は右から読み直して|01>, (4)の出力はそのまま|11>です。)

図のキャプションに書いた通りグラフのx軸に書かれたビット列は右から第一量子ビットで読みます。 そうすると入力のビット列と出力のビット列が交換されていることが図の(2)(3)を読んでいただくとわかっていただけると思います。

この時のソースは以下です。

qrs = [
    [0, 0],
    [0, 1],
    [1, 0],
    [1, 1]
]

fig_result, axes = plt.subplots(2, 2, figsize=(10, 10))
for idx, (q1, q2) in enumerate(qrs):
    qr_1 = QuantumRegister(1, 'q1')
    qr_2 = QuantumRegister(1, 'q2')
    cr = ClassicalRegister(2)
    qc = QuantumCircuit(qr_1, qr_2, cr)

    if q1 == 1:
        qc.x(qr_1)
    if q2 == 1:
        qc.x(qr_2)
    qc.barrier()
    
    qc.cx(qr_1, qr_2)
    qc.cx(qr_2, qr_1)
    qc.cx(qr_1, qr_2)

    qc.barrier()
    qc.measure(qr_1, cr[0])
    qc.measure(qr_2, cr[1])

    fig, ax = plt.subplots(1, 1)
    qc.draw(output='mpl', ax=ax)
    plt.show()
    plt.close()
    
    backend = Aer.get_backend('qasm_simulator')
    job = execute(qc, backend, shot=1024)
    result = job.result()

    plt.figure(fig_result)
    plot_histogram(result.get_counts(), ax=axes[q1][q2])
    axes[q1][q2].set_title(f'({idx+1}) 入力が|{q1}{q2}>の場合', fontsize=12)
    
plt.figure(fig_result)
plt.tight_layout()
plt.show()
plt.close()

重ね合わせ状態のSWAP

重ね合わせ状態もその係数ごと交換されます。 つまり

(α|0⟩+β|1⟩)1⊗(u|0⟩+v|1⟩)2⟷(u|0⟩+v|1⟩)1⊗(α|0⟩+β|1⟩)2${\displaystyle (\alpha |0⟩+\beta |1⟩{)}_1\otimes (u|0⟩+v|1⟩{)}_2⟷(u|0⟩+v|1⟩{)}_1\otimes (\alpha |0⟩+\beta |1⟩{)}_2}$

ができます。

図3. 重ね合わせ状態のSWAP (左) 実行する量子回路 (右) 観測結果

入力の状態が交換されて、出力は以下の式に対応する確率になっていれば良いです。

(12|0⟩+3–√2|1⟩)1⊗(3–√2|0⟩+12|1⟩)2=3–√4|00⟩+14|01⟩+34|10⟩+3–√4|11⟩${\displaystyle {(\frac{1}{2}|0⟩+\frac{\sqrt{3}}{2}|1⟩)}_1\otimes {(\frac{\sqrt{3}}{2}|0⟩+\frac{1}{2}|1⟩)}_2=\frac{\sqrt{3}}{4}|00⟩+\frac{1}{4}|01⟩+\frac{3}{4}|10⟩+\frac{\sqrt{3}}{4}|11⟩}$

ビット列の読み順は式とグラフとで逆順なのでグラフ上での|01⟩$|01⟩$が式上の|10⟩$|10⟩$に対応していて確率振幅0.75の自乗で約0.56、グラフ上の|10⟩$|10⟩$が式上の|01⟩$|01⟩$に対応していて確率振幅0.25の自乗で約0.05、|00⟩$|00⟩$と|11⟩$|11⟩$はそのままで確率振幅0.43の自乗で0.18に近いので、計算通りでしょう。

ということで、重ね合わせ状態もSWAPできることを確かめました。

この時のソースです。

qr_1 = QuantumRegister(1, 'q1')
qr_2 = QuantumRegister(1, 'q2')
cr = ClassicalRegister(2)
qc = QuantumCircuit(qr_1, qr_2, cr)

qc.ry(np.pi / 3, qr_1)
qc.ry(2 * np.pi / 3, qr_2)

qc.cx(qr_1, qr_2)
qc.cx(qr_2, qr_1)
qc.cx(qr_1, qr_2)

qc.barrier()

qc.measure(qr_1, cr[0])
qc.measure(qr_2, cr[1])

qc.draw(output='mpl')

backend = Aer.get_backend('qasm_simulator')
job = execute(qc, backend, shot=1024)

result = job.result()
plot_histogram(result.get_counts())

複数量子ビットへのSWAPゲート実装

感覚的にはこういう交換をしたいですね。

|000⟩1⊗|011⟩2⟷|011⟩1⊗|000⟩2${\displaystyle |000{⟩}_1\otimes |011{⟩}_2⟷|011{⟩}_1\otimes |000{⟩}_2}$

これを行うには互いの量子ビットのiビット目を交換すれば良いです。 従って下図のような量子回路になります。

図4. 3量子ビットのSWAP

SWAPが起きるか確認した結果です。

図5. 図4の量子回路の実行結果

第一量子ビットの状態が|011⟩$|011⟩$で第二量子ビットの状態が|000⟩$|000⟩$になっていることを示しています。 よってちゃんと思った通りにSWAPできていますね。

この時のソースです。

qr_1 = QuantumRegister(3, 'q1')
qr_2 = QuantumRegister(3, 'q2')
cr = ClassicalRegister(6)
qc = QuantumCircuit(qr_1, qr_2, cr)

qc.x(qr_2[1])
qc.x(qr_2[2])

qc.barrier()

qc.cx(qr_1[0], qr_2[0])
qc.cx(qr_2[0], qr_1[0])
qc.cx(qr_1[0], qr_2[0])

qc.barrier()

qc.cx(qr_1[1], qr_2[1])
qc.cx(qr_2[1], qr_1[1])
qc.cx(qr_1[1], qr_2[1])

qc.barrier()

qc.cx(qr_1[2], qr_2[2])
qc.cx(qr_2[2], qr_1[2])
qc.cx(qr_1[2], qr_2[2])

qc.barrier()

qc.measure(qr_1, cr[:len(qr_1)])
qc.measure(qr_2, cr[len(qr_1):])

qc.draw(output='mpl')

backend = Aer.get_backend('qasm_simulator')
job = execute(qc, backend, shot=1024)
result = job.result()

plot_histogram(result.get_counts())

制御SWAPゲート

感覚的にはこんな感じですよね。

図6. 制御SWAPゲートのイメージ

SWAPゲートの中身はCXゲートなので、CXがかかるのに必要な制御量子ビットを一つ増やすことで実装可能です。

図7. 制御SWAPゲートの実装

SWAPが起きることを確認します。

まずはSWAPが起きない場合です。

図8. 制御SWAPが起きない場合 (左) 実行する量子回路 (右) 測定結果

入力状態が|010⟩$|010⟩$ (第一量子ビットが制御量子ビット)に対して、測定結果は|010⟩$|010⟩$のため、SWAPは起きていません。

この時のソースです。

qr_control = QuantumRegister(1, 'control')
qr_1 = QuantumRegister(1, 'q1')
qr_2 = QuantumRegister(1, 'q2')
cr = ClassicalRegister(3)
qc = QuantumCircuit(qr_control, qr_1, qr_2, cr)

qc.x(qr_1)

qc.barrier()

qc.ccx(qr_control, qr_1, qr_2)
qc.ccx(qr_control, qr_2, qr_1)
qc.ccx(qr_control, qr_1, qr_2)

qc.barrier()

qc.measure(qr_control, cr[0])
qc.measure(qr_1, cr[1])
qc.measure(qr_2, cr[2])

qc.draw(output='mpl')

backend = Aer.get_backend('qasm_simulator')
job = execute(qc, backend, shot=1024)

result = job.result()
plot_histogram(result.get_counts())

そしてSWAPが起きる場合です。

図9. 制御SWAPが起きる場合 (左) 実行する量子回路 (右) 測定結果

入力状態が|110⟩$|110⟩$ (第一量子ビットが制御量子ビット)に対して、測定結果は|101⟩$|101⟩$のため、SWAPが起きています。

この時のソースです。

qr_control = QuantumRegister(1, 'control')
qr_1 = QuantumRegister(1, 'q1')
qr_2 = QuantumRegister(1, 'q2')
cr = ClassicalRegister(3)
qc = QuantumCircuit(qr_control, qr_1, qr_2, cr)

qc.x(qr_control)
qc.x(qr_1)

qc.barrier()

qc.ccx(qr_control, qr_1, qr_2)
qc.ccx(qr_control, qr_2, qr_1)
qc.ccx(qr_control, qr_1, qr_2)

qc.barrier()

qc.measure(qr_control, cr[0])
qc.measure(qr_1, cr[1])
qc.measure(qr_2, cr[2])

qc.draw(output='mpl')

backend = Aer.get_backend('qasm_simulator')
job = execute(qc, backend, shot=1024)

result = job.result()
plot_histogram(result.get_counts())

以上より、制御SWAPゲートの動作を確認できました。

効率の良い制御SWAPゲート

よく考えたら制御SWAPゲートは以下で簡単に実装できました。

図10. 効率の良い制御SWAPゲート(SWAPしない場合) (左) 実行する量子回路 (右) 測定結果

図11. 効率の良い制御SWAPゲート(SWAPする場合) (左) 実行する量子回路 (右) 測定結果

ソースは割愛しますが、最初と最後のCCXをただのCXに直しました。

CCXや2つ以上の制御量子ビットを持つ制御ユニタリーゲートは実際には多数の制御ユニタリーゲートに分解して実行されます。 (これは本質的には物理では2体問題しか扱えないからですが、詳細は省きます。)

arxiv.org

この事情によって、ソフトウェア的にはどちらでも良い制御SWAPゲートの実装ですが、ハードウェアの事情を加味するとゲートの数が増えるほどゲートエラーが発生するリスクが高まります。

実際に比較してみました。

図12. 最初の制御SWAPゲートを実機で実行する場合 (左) 実行する量子回路 (右) 測定結果

図13. 高効率の制御SWAPゲートを実行する場合 (左) 実行する量子回路 (右) 測定結果

図14. 実行するibmq_quitoのトポロジー

not usedと記した量子ビットは実際には未使用ということではなく、実行するibmq_quitoの1の量子ビットを想定して、ancilaとして使われる想定です。 ただ、実行したjobを確認すると思ったようには動いてなさそうなので、transplieの仕組みを理解する必要があります。

ひとまず実行結果を確認すると期待していた|1001⟩$|1001⟩$状態の確率は後者の回路の方が高いことがわかりました。 やはり、なるべくゲート数を減らすことが重要であることがわかりました。

嶋田さんの本 を読むとこの回路はFredkinゲートというそうです。

• フレドキンゲート - Wikipedia
• https://www.science.org/doi/10.1126/sciadv.1501531

まとめ

• SWAPゲートは入力した量子ビットの状態を交換する
• 重ね合わせ状態の確率振幅も交換できる
• 複数量子ビットのSWAPゲートは対応する量子ビット同士でそれぞれSWAPゲートを実行する
• 制御SWAPゲートは制御量子ビットの状態に応じてSWAPゲートの実行を切り替える
• 制御SWAPゲートはFredkinゲートで実装できる

次回はこのSWAPゲート・制御SWAPゲートを使ってスワップテストを解説します。

アダマールテスト

前回までで測定・基底・期待値・密度演算子など基本的な理解が進んだので、アルゴリズムの理解をしたいと思います。

今回はアダマールテスト。
後々に量子位相推定やHHLアルゴリズムなど種々のアルゴリズムに発展する基本的なアルゴリズムのようなのでやってみます。

教科書ではサラッと書かれていますが、ちゃんと計算しようとすると密度演算子の知識が必要です。

対象とする読者

• 量子プログラミングに興味を持ち始めたばかりのエンジニア
• 量子計算に興味を持ち始めたばかりの高校生・大学生 など

前提とする知識

• 簡単な論理演算(AND, OR, NOT)が分かること
• 基本的な量子ゲート(X, CX, CCX)が分かること
• ブラケットでの計算がわかること
• 密度演算子がわかること

バージョン情報

• Python 3.9.13
• Qiskit 0.20.2

目次

アダマールテストとは

量子回路

図1. アダマールテストの量子回路

アダマールテストとは上図のように第一量子ビットは初期状態|0⟩$|0⟩$で、それ以外の任意の数の量子ビットは任意の状態|ψ⟩$|\psi ⟩$を初期状態としています。
そして第一量子ビットにアダマールゲートをかけた後、第一量子ビットを制御ビットとしそれ以外の量子ビットに対して制御ユニタリーゲート(controlled - U$U$)をかけます。
最後に第一量子ビットにもう一度アダマールゲートをかけて第一量子ビットを測定するというものです。

用途

この回路を使うとユニタリー行列U$U$の固有値がわかります。

ただ、後で詳細説明しますが、量子回路に使われるユニタリー行列の固有値は±1$\pm 1$なので、あえて量子コンピュータを使って求めるものではないです。 (私の私見のため、正しくないかもしれません。詳しい方がいたら補足お願いします。) アダマールテストを利用したスワップテストや量子位相推定アルゴリズムのことを考慮すると、むしろここで大事なのは「位相キックバック」と呼ばれる第一量子ビットの位相にユニタリーU$U$の固有値に相当する位相が乗る現象によって、第一量子ビットの確率分布がその分変化するということが重要と思います。

代数的に解いてみる

準備：ユニタリー行列とその固有値

ユニタリー行列とは

まず、ユニタリー行列は元々線形代数の言葉なので量子力学特有の名称・行列ではないです。

ユニタリ行列 - Wikipedia

によると、自身と自身のエルミート共役(転置して複素共役を取ったもの)の積が単位行列になるものです。

U^U^†=U^†U^=I^${\displaystyle \hat{U}{\hat{U}}^{†}={\hat{U}}^{†}\hat{U}=\hat{I}}$

数学的な定義だと僕も実はよくわからないのですが、物理学的には量子力学的に有効な操作(演算、つまり量子ビットをBloch球上で回転させること)であれば全てユニタリーです。 アダマールゲートやXゲートなどを行列でH^,X^$\hat{H},\hat{X}$と表すなら、それらの逆回転がエルミート共役H^†,X^†${\hat{H}}^{†},{\hat{X}}^{†}$です。

もう少し細かく話すと、量子状態の時間発展は以下のシュレディンガー方程式で記述されます。

iℏddt|ψ(t)⟩=H^|ψ(t)⟩${\displaystyle i\hslash \frac{d}{dt}|\psi (t)⟩=\hat{H}|\psi (t)⟩}$

この式のH^$\hat{H}$はハミルトニアンを示します。 この式の解は初期時刻t0$t_0$を用いて一般的に下記の形式で表されます。

|ψ(t)⟩=exp{1iℏH^(t−t0)}|ψ(t0)⟩${\displaystyle |\psi (t)⟩=\exp \{\frac{1}{i\hslash }\hat{H}(t-t_0)\}|\psi (t_0)⟩}$

この指数関数部分がユニタリー行列で表されます。 ハミルトニアンで記述される量子系はエネルギー保存している(だからハミルトニアンで書いているのですけど)ので、時間を巻き戻したら元に戻りますよねって話です。 物理学の言葉で「時間反転対称性」と言います。 巻き戻したら

exp{1iℏH^(t−t0)}⟶exp{−1iℏH^(t−t0)}${\displaystyle \exp \{\frac{1}{i\hslash }\hat{H}(t-t_0)\}⟶\exp \{\frac{-1}{i\hslash }\hat{H}(t-t_0)\}}$

となるだけで、これは

expiΔθ⟶exp(−iΔθ)${\displaystyle \exp i\mathrm{\Delta }\theta ⟶\exp (-i\mathrm{\Delta }\theta )}$

の変換なので、逆回転ということです。

時間反転について、ここでは複素共役しか言及してませんが、転置についてはハミルトニアンがエルミート行列(自身と自身のエルミート共役が同一)であることから言及しなくて良くなっています。

exp(iH^)=∑k=0∞1k!(iH^)k≡U^${\displaystyle \exp (i\hat{H})=\sum _{k=0}^{\mathrm{\infty }}\frac{1}{k!}(i\hat{H}{)}^k\equiv \hat{U}}$

とすると

U^†=∑k=0∞1k!(−iH^†)k=∑k=0∞1k!(−iH^)k=exp(−iH^)$\begin{array}{rl}{\hat{U}}^{†}& =\sum _{k=0}^{\mathrm{\infty }}\frac{1}{k!}(-i{\hat{H}}^{†}{)}^k\\ & =\sum _{k=0}^{\mathrm{\infty }}\frac{1}{k!}(-i\hat{H}{)}^k\\ & =\exp (-i\hat{H})\end{array}$

ユニタリー行列の固有値

これは数学の話なので、以下の記事での解説に譲ります。

risalc.info

eman-physics.net

エッセンスとして重要なことは以下です。

• ユニタリ行列の固有値は大きさが1の複素数であること
• ユニタリ行列の固有値はeiλ$e^{i\lambda }$で表される複素数であること

これらのことから、ユニタリ行列をその固有ベクトルである量子状態にかけた場合以下になります。

U^|ψ⟩=eiθ|ψ⟩${\displaystyle \hat{U}|\psi ⟩=e^{i\theta }|\psi ⟩}$

この固有値eiλ$e^{i\lambda }$の位相λ$\lambda$が第一量子ビットの確率分布から分かるというのがアダマールテストです。

部分トレース

今回は第一量子ビットとそれ以外の量子ビットがあり、第一量子ビットだけ測定するというスキームですので、このような場合は部分トレースをとる考え方が必要になります。 whyitsso.net

部分トレースについては前回の記事でも解説しています。

sakumatcho.hatenablog.com

入力状態がユニタリー行列の固有値の場合

入力状態|ψ⟩$|\psi ⟩$がユニタリ行列U^$\hat{U}$の固有ベクトルだと仮定します。
するとU^|ψ⟩=eiλ|ψ⟩$\hat{U}|\psi ⟩=e^{i\lambda }|\psi ⟩$が成り立ちますので、これを使って式変形していきます。

12–√(|0⟩+|1⟩)⊗|ψ⟩⟶Ctrl-U12–√(|0⟩⊗|ψ⟩+|1⟩⊗U^|ψ⟩)=12–√(|0⟩⊗|ψ⟩+eiλ|1⟩⊗|ψ⟩)=12–√(|0⟩+eiλ|1⟩)⊗|ψ⟩$\begin{array}{rl}\frac{1}{\sqrt{2}}(|0⟩+|1⟩)\otimes |\psi ⟩& \stackrel{\text{Ctrl-}U}{⟶}\frac{1}{\sqrt{2}}(|0⟩\otimes |\psi ⟩+|1⟩\otimes \hat{U}|\psi ⟩)\\ & =\frac{1}{\sqrt{2}}(|0⟩\otimes |\psi ⟩+e^{i\lambda }|1⟩\otimes |\psi ⟩)\\ & =\frac{1}{\sqrt{2}}(|0⟩+e^{i\lambda }|1⟩)\otimes |\psi ⟩\end{array}$

この後、第一量子ビットにアダマールゲートをかけて整理すると以下になります。

12–√{12–√(|0⟩+|1⟩)+eiλ12–√(|0⟩−|1⟩)}=1+eiλ2|0⟩+1−eiλ2|1⟩${\displaystyle \frac{1}{\sqrt{2}}\{\frac{1}{\sqrt{2}}(|0⟩+|1⟩)+e^{i\lambda }\frac{1}{\sqrt{2}}(|0⟩-|1⟩)\}=\frac{1+e^{i\lambda }}{2}|0⟩+\frac{1-e^{i\lambda }}{2}|1⟩}$

従って、第一量子ビットを測定してその結果が0又は1となる確率はそれぞれ以下になります。

p0=∣∣∣1+eiλ2∣∣∣2=1+cosλ2p1=∣∣∣1−eiλ2∣∣∣2=1−cosλ2$\begin{array}{rl}& p_0={\left|\frac{1+e^{i\lambda }}{2}\right|}^2=\frac{1+\cos \lambda }{2}\\ & p_1={\left|\frac{1-e^{i\lambda }}{2}\right|}^2=\frac{1-\cos \lambda }{2}\end{array}$

この式と実際に測定した結果を照らし合わせてλ$\lambda$を推定します、というか実際には既知のユニタリーを実行するので、ユニタリーを利用してこの分布をコントロールします。

入力状態がユニタリー行列の固有値ではない場合

入力状態|ψ⟩$|\psi ⟩$がユニタリー行列U^$\hat{U}$の固有状態では無い場合も

12–√(|0⟩+|1⟩)⊗|ψ⟩⟶Ctrl-U12–√(|0⟩⊗|ψ⟩+|1⟩⊗U^|ψ⟩)${\displaystyle \frac{1}{\sqrt{2}}(|0⟩+|1⟩)\otimes |\psi ⟩\stackrel{\text{Ctrl-}U}{⟶}\frac{1}{\sqrt{2}}(|0⟩\otimes |\psi ⟩+|1⟩\otimes \hat{U}|\psi ⟩)}$

ここまでは固有状態である場合と同じですが、この後アダマールゲートを第一量子ビットにかけると状況がちょっと変わります。

12–√{(|0⟩+|1⟩)⊗|ψ⟩+(|0⟩−|1⟩)⊗U^|ψ⟩}=12{|0⟩⊗(|ψ⟩+U^|ψ⟩)+|1⟩⊗(|ψ⟩−U^|ψ⟩)}${\displaystyle \frac{1}{\sqrt{2}}\{(|0⟩+|1⟩)\otimes |\psi ⟩+(|0⟩-|1⟩)\otimes \hat{U}|\psi ⟩\}=\frac{1}{2}\{|0⟩\otimes (|\psi ⟩+\hat{U}|\psi ⟩)+|1⟩\otimes (|\psi ⟩-\hat{U}|\psi ⟩)\}}$

この後、第一量子ビットだけを測定したときに第一量子ビットの|0⟩,|1⟩$|0⟩,|1⟩$状態それぞれが観測される確率を求めるため、密度演算子をこの状態から求めます。

ρ^AB=14{|0⟩A⊗(I^+U^)|ψ⟩B+|1⟩A⊗(I^−U^)|ψ⟩B}{⟨0|A⊗⟨ψ|B(I^+U^)†+⟨1|A⊗⟨ψ|B(I^−U^)†}$$\begin{gathered} {\displaystyle {\hat{\rho }}^{AB}=\frac{1}{4}\{|0{⟩}_A\otimes (\hat{I}+\hat{U})|\psi {⟩}_B+|1{⟩}_A\otimes (\hat{I}-\hat{U})|\psi {⟩}_B\}\{⟨0{|}_A\otimes ⟨\psi {|}_B{(\hat{I}+\hat{U})}^{†}+⟨1{|}_A\otimes ⟨\psi {|}_B{(\hat{I}-\hat{U})}^{†}\} \\ } \end{gathered}$$

から

ρ^A=======14trB{|0⟩A⊗(I^+U^)|ψ⟩B⟨0|A⊗⟨ψ|B(I^+U^)†}+14trB{|0⟩A⊗(I^+U^)|ψ⟩B⟨1|A⊗⟨ψ|B(I^−U^)†}+14trB{|1⟩A⊗(I^−U^)|ψ⟩B⟨0|A⊗⟨ψ|B(I^+U^)†}+14trB{|1⟩A⊗(I^−U^)|ψ⟩B⟨1|A⊗⟨ψ|B(I^−U^)†}14trB{|0⟩A⟨0|A⊗(I^+U^)|ψ⟩B⟨ψ|B(I^+U^)†}+14trB{|0⟩A⟨1|A⊗(I^+U^)|ψ⟩B⟨ψ|B(I^−U^)†}+14trB{|1⟩A⟨0|A⊗(I^−U^)|ψ⟩B⟨ψ|B(I^+U^)†}+14trB{|1⟩A⟨1|A⊗(I^−U^)|ψ⟩B⟨ψ|B(I^−U^)†}14|0⟩A⟨0|Atr{(I^+U^)|ψ⟩B⟨ψ|B(I^+U^)†}+14|0⟩A⟨1|Atr{(I^+U^)|ψ⟩B⟨ψ|B(I^−U^)†}+14|1⟩A⟨0|Atr{(I^−U^)|ψ⟩B⟨ψ|B(I^+U^)†}+14|1⟩A⟨1|Atr{(I^−U^)|ψ⟩B⟨ψ|B(I^−U^)†}14|0⟩A⟨0|A⟨ψ|B(I^+U^)(I^+U^)†|ψ⟩B+14|0⟩A⟨1|A⟨ψ|B(I^−U^)(I^+U^)†|ψ⟩B+14|1⟩A⟨0|A⟨ψ|B(I^+U^)(I^−U^)†|ψ⟩B+14|1⟩A⟨1|A⟨ψ|B(I^−U^)(I^−U^)†|ψ⟩B14|0⟩A⟨0|A(⟨ψ|B2I^|ψ⟩B+⟨ψ|BU^|ψ⟩B+⟨ψ|BU^†|ψ⟩B)+略+14|1⟩A⟨1|A(⟨ψ|B2I^|ψ⟩B−⟨ψ|BU^|ψ⟩B−⟨ψ|BU^†|ψ⟩B)14|0⟩A⟨0|A{2+⟨ψ|BU^|ψ⟩B+(⟨ψ|BU^|ψ⟩B)∗}+略+14|1⟩A⟨1|A{2−⟨ψ|BU^|ψ⟩B−(⟨ψ|BU^|ψ⟩B)∗}14|0⟩A⟨0|A{2+Re(⟨ψ|BU^|ψ⟩B)}+略+14|1⟩A⟨1|A{2−Re(⟨ψ|BU^|ψ⟩B)}$\begin{array}{rl}{\hat{\rho }}^A=& \frac{1}{4}{\mathrm{t}\mathrm{r}}_B\{|0{⟩}_A\otimes (\hat{I}+\hat{U})|\psi {⟩}_B⟨0{|}_A\otimes ⟨\psi {|}_B{(\hat{I}+\hat{U})}^{†}\}\\ & +\frac{1}{4}{\mathrm{t}\mathrm{r}}_B\{|0{⟩}_A\otimes (\hat{I}+\hat{U})|\psi {⟩}_B⟨1{|}_A\otimes ⟨\psi {|}_B{(\hat{I}-\hat{U})}^{†}\}\\ & +\frac{1}{4}{\mathrm{t}\mathrm{r}}_B\{|1{⟩}_A\otimes (\hat{I}-\hat{U})|\psi {⟩}_B⟨0{|}_A\otimes ⟨\psi {|}_B{(\hat{I}+\hat{U})}^{†}\}\\ & +\frac{1}{4}{\mathrm{t}\mathrm{r}}_B\{|1{⟩}_A\otimes (\hat{I}-\hat{U})|\psi {⟩}_B⟨1{|}_A\otimes ⟨\psi {|}_B{(\hat{I}-\hat{U})}^{†}\}\\ =& \frac{1}{4}{\mathrm{t}\mathrm{r}}_B\{|0{⟩}_A⟨0{|}_A\otimes (\hat{I}+\hat{U})|\psi {⟩}_B⟨\psi {|}_B{(\hat{I}+\hat{U})}^{†}\}\\ & +\frac{1}{4}{\mathrm{t}\mathrm{r}}_B\{|0{⟩}_A⟨1{|}_A\otimes (\hat{I}+\hat{U})|\psi {⟩}_B⟨\psi {|}_B{(\hat{I}-\hat{U})}^{†}\}\\ & +\frac{1}{4}{\mathrm{t}\mathrm{r}}_B\{|1{⟩}_A⟨0{|}_A\otimes (\hat{I}-\hat{U})|\psi {⟩}_B⟨\psi {|}_B{(\hat{I}+\hat{U})}^{†}\}\\ & +\frac{1}{4}{\mathrm{t}\mathrm{r}}_B\{|1{⟩}_A⟨1{|}_A\otimes (\hat{I}-\hat{U})|\psi {⟩}_B⟨\psi {|}_B{(\hat{I}-\hat{U})}^{†}\}\\ =& \frac{1}{4}|0{⟩}_A⟨0{|}_A\mathrm{t}\mathrm{r}\{(\hat{I}+\hat{U})|\psi {⟩}_B⟨\psi {|}_B{(\hat{I}+\hat{U})}^{†}\}\\ & +\frac{1}{4}|0{⟩}_A⟨1{|}_A\mathrm{t}\mathrm{r}\{(\hat{I}+\hat{U})|\psi {⟩}_B⟨\psi {|}_B{(\hat{I}-\hat{U})}^{†}\}\\ & +\frac{1}{4}|1{⟩}_A⟨0{|}_A\mathrm{t}\mathrm{r}\{(\hat{I}-\hat{U})|\psi {⟩}_B⟨\psi {|}_B{(\hat{I}+\hat{U})}^{†}\}\\ & +\frac{1}{4}|1{⟩}_A⟨1{|}_A\mathrm{t}\mathrm{r}\{(\hat{I}-\hat{U})|\psi {⟩}_B⟨\psi {|}_B{(\hat{I}-\hat{U})}^{†}\}\\ =& \frac{1}{4}|0{⟩}_A⟨0{|}_A⟨\psi {|}_B(\hat{I}+\hat{U}){(\hat{I}+\hat{U})}^{†}|\psi {⟩}_B\\ & +\frac{1}{4}|0{⟩}_A⟨1{|}_A⟨\psi {|}_B(\hat{I}-\hat{U}){(\hat{I}+\hat{U})}^{†}|\psi {⟩}_B\\ & +\frac{1}{4}|1{⟩}_A⟨0{|}_A⟨\psi {|}_B(\hat{I}+\hat{U}){(\hat{I}-\hat{U})}^{†}|\psi {⟩}_B\\ & +\frac{1}{4}|1{⟩}_A⟨1{|}_A⟨\psi {|}_B(\hat{I}-\hat{U}){(\hat{I}-\hat{U})}^{†}|\psi {⟩}_B\\ =& \frac{1}{4}|0{⟩}_A⟨0{|}_A(⟨\psi {|}_{B}2\hat{I}|\psi {⟩}_B+⟨\psi {|}_B\hat{U}|\psi {⟩}_B+⟨\psi {|}_B{\hat{U}}^{†}|\psi {⟩}_B)\\ & +略\\ & +\frac{1}{4}|1{⟩}_A⟨1{|}_A(⟨\psi {|}_{B}2\hat{I}|\psi {⟩}_B-⟨\psi {|}_B\hat{U}|\psi {⟩}_B-⟨\psi {|}_B{\hat{U}}^{†}|\psi {⟩}_B)\\ =& \frac{1}{4}|0{⟩}_A⟨0{|}_A\{2+⟨\psi {|}_B\hat{U}|\psi {⟩}_B+{(⟨\psi {|}_B\hat{U}|\psi {⟩}_B)}^{\ast }\}\\ & +略\\ & +\frac{1}{4}|1{⟩}_A⟨1{|}_A\{2-⟨\psi {|}_B\hat{U}|\psi {⟩}_B-{(⟨\psi {|}_B\hat{U}|\psi {⟩}_B)}^{\ast }\}\\ =& \frac{1}{4}|0{⟩}_A⟨0{|}_A\{2+\mathrm{R}\mathrm{e}(⟨\psi {|}_B\hat{U}|\psi {⟩}_B)\}\\ & +略\\ & +\frac{1}{4}|1{⟩}_A⟨1{|}_A\{2-\mathrm{R}\mathrm{e}(⟨\psi {|}_B\hat{U}|\psi {⟩}_B)\}\end{array}$

よって、第一量子ビットを測定したときに|0⟩A$|0{⟩}_A$が観測される確率p0$p_0$および|1⟩A$|1{⟩}_A$が観測される確率p1$p_1$はそれぞれ以下になります。

p0=1+Re⟨ψ|BU^|ψ⟩B2${\displaystyle p_0=\frac{1+\mathrm{R}\mathrm{e}⟨\psi {|}_B\hat{U}|\psi {⟩}_B}{2}}$

p1=1−Re⟨ψ|BU^|ψ⟩B2${\displaystyle p_1=\frac{1-\mathrm{R}\mathrm{e}⟨\psi {|}_B\hat{U}|\psi {⟩}_B}{2}}$

ということでやっと教科書に書いてあるような式が現れました。
教科書ではサラッと書かれていますが、密度演算子と行列計算を使いこなせる必要があります。

制御ユニタリーがかかった後の量子ビットの状態

入力する量子ビットの状態は固有状態ではないが、(固有値がλ1,λ2${\lambda }_1,{\lambda }_2$で、固有状態が|u1⟩,|u2⟩$|u_1⟩,|u_2⟩$と予めわかっているとすると)以下のように書き換えられます。

|ψ⟩=c1|u1⟩+c2|u2⟩${\displaystyle |\psi ⟩=c_1|u_1⟩+c_2|u_2⟩}$

そうすると、2回目のアダマールゲートをかけた後の標的量子ビットの状態は以下になります。

|ψ⟩±U^|ψ⟩=c1|u1⟩+c2|u2⟩±c1U^|u1⟩±c2U^|u2⟩=c1|u1⟩+c2|u2⟩±c1eiλ1|u1⟩±c2eiλ2|u2⟩=c1(1±eiλ1)|u1⟩+c2(1±eiλ2)|u2⟩$\begin{array}{rl}|\psi ⟩\pm \hat{U}|\psi ⟩& =c_1|u_1⟩+c_2|u_2⟩\pm c_1\hat{U}|u_1⟩\pm c_2\hat{U}|u_2⟩\\ & =c_1|u_1⟩+c_2|u_2⟩\pm c_1{e}^{i{\lambda }_1}|u_1⟩\pm c_2{e}^{i{\lambda }_2}|u_2⟩\\ & =c_1(1\pm e^{i{\lambda }_1})|u_1⟩+c_2(1\pm e^{i{\lambda }_2})|u_2⟩\end{array}$

上記のように任意のユニタリーとそれに対応する固有値で表現するとわかりにくいですが、よく使うX^,Y^,Z^$\hat{X},\hat{Y},\hat{Z}$ゲートであれば固有値はλ=0,π$\lambda =0,\pi$で固有ベクトルはそれぞれの軸上の単位ベクトルになりますので、上記の式で言えば|u1⟩,|u2⟩$|u_1⟩,|u_2⟩$のどちらか一方だけが残ります。
従って、もう一つアダマールテスト回路を作成すれば「入力状態がユニタリー行列の固有値の場合」に帰着します。

以下で具体的に計算してみます。

制御Xゲートの場合

U^=X^$\hat{U}=\hat{X}$の固有ベクトルは|+⟩,|−⟩$|+⟩,|-⟩$で、対応する固有値は+1,−1$+1,-1$です。 従って

|ψ⟩+X^|ψ⟩=2c1|+⟩${\displaystyle |\psi ⟩+\hat{X}|\psi ⟩=2c_1|+⟩}$

|ψ⟩−X^|ψ⟩=2c2|−⟩${\displaystyle |\psi ⟩-\hat{X}|\psi ⟩=2c_2|-⟩}$

となります。

制御Zゲートの場合

U^=Z^$\hat{U}=\hat{Z}$の固有ベクトルは|0⟩,|1⟩$|0⟩,|1⟩$で、対応する固有値は+1,−1$+1,-1$です。 従って

|ψ⟩+Z^|ψ⟩=2c1|0⟩${\displaystyle |\psi ⟩+\hat{Z}|\psi ⟩=2c_1|0⟩}$

|ψ⟩−Z^|ψ⟩=2c2|1⟩${\displaystyle |\psi ⟩-\hat{Z}|\psi ⟩=2c_2|1⟩}$

となります。

制御Yゲートの場合

U^=Y^$\hat{U}=\hat{Y}$の固有ベクトルは|+⟩+i|−⟩,|+⟩−i|−⟩$|+⟩+i|-⟩,|+⟩-i|-⟩$ (ポアンカレ球の斜め直線偏光に相当する状態)で、対応する固有値は+1,−1$+1,-1$です。 従って

|ψ⟩+Y^|ψ⟩=2c1(|+⟩+i|−⟩)${\displaystyle |\psi ⟩+\hat{Y}|\psi ⟩=2c_1(|+⟩+i|-⟩)}$

|ψ⟩−Y^|ψ⟩=2c2(|+⟩−i|−⟩)${\displaystyle |\psi ⟩-\hat{Y}|\psi ⟩=2c_2(|+⟩-i|-⟩)}$

となります。

量子コンピュータで計算させてみる

入力状態がユニタリー行列の固有値の場合

ユニタリー行列が制御Xゲートの場合

U^=X^$\hat{U}=\hat{X}$の固有ベクトルは|+⟩,|−⟩$|+⟩,|-⟩$で、対応する固有値は+1,−1$+1,-1$です。 従ってλ+=0, λ−=π${\lambda }_{+}=0,{\lambda }_{-}=\pi$のため、第一量子ビットの確率は以下となります。

i) 第二量子ビットの状態が|ψ⟩=|+⟩$|\psi ⟩=|+⟩$の時

p0=1+12=1$\begin{array}{rl}{p}_0& =\frac{1+1}{2}\\ & =1\end{array}$

ii) 第二量子ビットの状態が|ψ⟩=|−⟩$|\psi ⟩=|-⟩$の時

p0=1−12=0$\begin{array}{rl}{p}_0& =\frac{1-1}{2}\\ & =0\end{array}$

i)について実行する量子回路です。

図1. アダマールテスト：ユニタリーがXゲートで、標的第二量子ビットがその固有ベクトル| + >の場合

シミュレータで実行した結果です。

図2. 図1の量子回路の実行結果(シミュレータ)

実際にIBMQの実機の結果です。

図3. 図1の量子回路の実行結果(実機)

ということで、予め計算した通りになりました。

これらのソースコードです。

qr_control = QuantumRegister(1, 'control')
qr_target = QuantumRegister(1, 'target')
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr_control, qr_target, cr)

qc.h(qr_target)

qc.barrier()

qc.h(qr_control)
qc.cx(qr_control, qr_target)
qc.h(qr_control)

qc.barrier()
qc.measure(qr_control, cr[0])

qc.draw(output='mpl')

## シミュレータ

backend = Aer.get_backend('qasm_simulator')
job = execute(qc, backend, shot=1024)

result = job.result()
plot_histogram(result.get_counts())

## 実機(マシンはその時空いているもので良いです)
backend = provider.get_backend('ibmq_belem')
job = execute(qc, backend, shots=1024)

result = job.result()
plot_histogram(result.get_counts())

とはいえ、i)の場合アダマール２回やって|0⟩$|0⟩$に戻ってきただけ、と思われても嫌なのでii)も計算しておきます。

ii)について実行する量子回路です。

図4. アダマールテスト：ユニタリーがXゲートで、標的量子ビットがその固有ベクトルの| - >の場合

シミュレータで実行した結果です。

図5. 図4の量子回路の実行結果(シミュレータ)

実機の結果です。

図6. 図4の量子回路の実行結果(実機)

ということで、第一量子ビットからしたらアダマール２回実行したようにしか見えませんがちゃんと第二量子ビットの影響を受けていることがわかります。

この時のソースコードです。

qr_control = QuantumRegister(1, 'control')
qr_target = QuantumRegister(1, 'target')
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr_control, qr_target, cr)

qc.x(qr_target)
qc.h(qr_target)

qc.barrier()

qc.h(qr_control)
qc.cx(qr_control, qr_target)
qc.h(qr_control)

qc.barrier()
qc.measure(qr_control, cr[0])

qc.draw(output='mpl')

### シミュレータ
backend = Aer.get_backend('qasm_simulator')
job = execute(qc, backend, shot=1024)

result = job.result()
plot_histogram(result.get_counts())

### 実機(その時に空いているマシンなら何でも可)
backend = provider.get_backend('ibmq_belem')
job = execute(qc, backend, shots=1024)

result = job.result()
plot_histogram(result.get_counts())

ユニタリー行列がRxゲートの場合

同じように任意の回転ゲートでも確率分布が変化することを確認します。

例えばRx(π2)$R_x\left(\frac{\pi }{2}\right)$の場合、固有値は(1±i)/2–√$(1\pm i)/\sqrt{2}$に対して固有ベクトルは|−⟩,|+⟩$|-⟩,|+⟩$です。

従ってλ+=π/4, λ−=−π/4${\lambda }_{+}=\pi /4,{\lambda }_{-}=-\pi /4$のため、第一量子ビットの確率は以下となります。

i) 第二量子ビットの状態が|ψ⟩=|−⟩$|\psi ⟩=|-⟩$の時、λ=λ+$\lambda ={\lambda }_{+}$のため

p0=1+2–√/22=2+2–√4∼0.85$\begin{array}{rl}{p}_0& =\frac{1+\sqrt{2}/2}{2}\\ & =\frac{2+\sqrt{2}}{4}\\ & \sim 0.85\end{array}$

p1=1−2–√/22=2−2–√4∼0.15$\begin{array}{rl}{p}_1& =\frac{1-\sqrt{2}/2}{2}\\ & =\frac{2-\sqrt{2}}{4}\\ & \sim 0.15\end{array}$

ii) 第二量子ビットの状態が|ψ⟩=|+⟩$|\psi ⟩=|+⟩$の時、λ=λ−$\lambda ={\lambda }_{-}$のため

p0=1+2–√/22=2+2–√4∼0.15$\begin{array}{rl}{p}_0& =\frac{1+\sqrt{2}/2}{2}\\ & =\frac{2+\sqrt{2}}{4}\\ & \sim 0.15\end{array}$

p1=1−2–√/22=2−2–√4∼0.15$\begin{array}{rl}{p}_1& =\frac{1-\sqrt{2}/2}{2}\\ & =\frac{2-\sqrt{2}}{4}\\ & \sim 0.15\end{array}$

i)について実行する量子回路は以下です。

図7. アダマールテスト：ユニタリーがRx(π/2)ゲートで、標的量子ビットがその固有ベクトルの| - >の場合

シミュレータで実行した結果です。

図8. 図7の量子回路の実行結果(シミュレータ)

実機の結果です。

図9. 図7の量子回路の実行結果(実機)

この時のソースコードです。

qr_control = QuantumRegister(1, 'control')
qr_target = QuantumRegister(1, 'target')
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr_control, qr_target, cr)

qc.h(qr_target)

qc.barrier()

qc.h(qr_control)
qc.crx(np.pi / 2, qr_control, qr_target)
qc.h(qr_control)

qc.barrier()
qc.measure(qr_control, cr[0])

qc.draw(output='mpl')

## 実行(シミュレーション)
backend = Aer.get_backend('qasm_simulator')
job = execute(qc, backend, shot=1024)

result = job.result()
plot_histogram(result.get_counts())

## 実行(実機)
backend = provider.get_backend('ibmq_quito')
job = execute(qc, backend, shots=1024)

result = job.result()
plot_histogram(result.get_counts())

ii)について実行する量子回路は以下です。

図10. アダマールテスト：ユニタリーがRx(π/2)ゲートで、標的量子ビットがその固有ベクトルの| + >の場合

図11. 図10の量子回路の実行結果(シミュレータ)

図12. 図10の量子回路の実行結果(実機)

この時のソースコードです。

qr_control = QuantumRegister(1, 'control')
qr_target = QuantumRegister(1, 'target')
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr_control, qr_target, cr)

qc.x(qr_target)
qc.h(qr_target)

qc.barrier()

qc.h(qr_control)
qc.crx(np.pi / 2, qr_control, qr_target)
qc.h(qr_control)

qc.barrier()
qc.measure(qr_control, cr[0])

qc.draw(output='mpl')

## 実行(シミュレータ)
backend = Aer.get_backend('qasm_simulator')
job = execute(qc, backend, shot=1024)

result = job.result()
plot_histogram(result.get_counts())

## 実行(実機)
backend = provider.get_backend('ibmq_quito')
job = execute(qc, backend, shots=1024)

result = job.result()
plot_histogram(result.get_counts())

ということで、予め求めた固有値の通りに確率分布が得られました。

入力状態がユニタリー行列の固有値ではない場合

入力状態は(|0⟩+i|1⟩)/2–√$(|0⟩+i|1⟩)/\sqrt{2}$でユニタリーはXゲートの場合

⟨ψ|X^|ψ⟩$⟨\psi |\hat{X}|\psi ⟩$が複素数になる場合です。

⟨ψ|X^|ψ⟩=12(1i)(0110)(1i)=12(i1)(1i)=i$\begin{array}{rl}⟨\psi |\hat{X}|\psi ⟩& =\frac{1}{2}\left(\begin{array}{cc}1& i\end{array}\right)\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}1\\ i\end{array}\right)\\ & =\frac{1}{2}\left(\begin{array}{cc}i& 1\end{array}\right)\left(\begin{array}{c}1\\ i\end{array}\right)\\ & =i\end{array}$

従って

p0=1+02=0.5$\begin{array}{rl}{p}_0& =\frac{1+0}{2}\\ & =0.5\end{array}$

p1=1−02=0.5$\begin{array}{rl}{p}_1& =\frac{1-0}{2}\\ & =0.5\end{array}$

この時の量子回路と実行回路です。

図13. アダマールテスト：ユニタリーがXゲートで、標的量子ビットがその固有ベクトルでない場合(1)

図14. 図13の量子回路の実行結果(シミュレータ)

図15. 図13の量子回路の実行結果(実機)

この時のソースコードです。

qr_control = QuantumRegister(1, 'control')
qr_target = QuantumRegister(1, 'target')
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr_control, qr_target, cr)

qc.rx(-np.pi / 3, qr_target)

qc.barrier()

qc.h(qr_control)
qc.cx(qr_control, qr_target)
qc.h(qr_control)

qc.barrier()
qc.measure(qr_control, cr[0])

qc.draw(output='mpl')

## 実行(シミュレータ)
backend = Aer.get_backend('qasm_simulator')
job = execute(qc, backend, shot=1024)

result = job.result()
plot_histogram(result.get_counts())

## 実行(実機)
backend = provider.get_backend('ibmq_quito')
job = execute(qc, backend, shots=1024)

result = job.result()
plot_histogram(result.get_counts())

入力状態は(3–√|0⟩+|1⟩)/2$(\sqrt{3}|0⟩+|1⟩)/2$でユニタリーはXゲートの場合

⟨ψ|X^|ψ⟩$⟨\psi |\hat{X}|\psi ⟩$が半端な実数になる場合です。

⟨ψ|X^|ψ⟩=14(3–√1)(0110)(3–√1)=14(13–√)(3–√1)=3–√2$\begin{array}{rl}⟨\psi |\hat{X}|\psi ⟩& =\frac{1}{4}\left(\begin{array}{cc}\sqrt{3}& 1\end{array}\right)\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}\sqrt{3}\\ 1\end{array}\right)\\ & =\frac{1}{4}\left(\begin{array}{cc}1& \sqrt{3}\end{array}\right)\left(\begin{array}{c}\sqrt{3}\\ 1\end{array}\right)\\ & =\frac{\sqrt{3}}{2}\end{array}$

従って

p0=1+3–√/22∼0.93$\begin{array}{rl}{p}_0& =\frac{1+\sqrt{3}/2}{2}\\ & \sim 0.93\end{array}$

p1=1−3–√/22∼0.07$\begin{array}{rl}{p}_1& =\frac{1-\sqrt{3}/2}{2}\\ & \sim 0.07\end{array}$

この時の量子回路と実行結果です。

図16. アダマールテスト：ユニタリーがXゲートで、標的量子ビットがその固有ベクトルでない場合(2)

図17. 図16の量子回路の実行結果(シミュレータ)

図18. 図16の量子回路の実行結果(実機)

この時のソースコードです。

qr_control = QuantumRegister(1, 'control')
qr_target = QuantumRegister(1, 'target')
cr = ClassicalRegister(1)
qc = QuantumCircuit(qr_control, qr_target, cr)

qc.ry(np.pi / 3, qr_target)

qc.barrier()

qc.h(qr_control)
qc.cx(qr_control, qr_target)
qc.h(qr_control)

qc.barrier()
qc.measure(qr_control, cr[0])

qc.draw(output='mpl')

## 実行(シミュレータ)
backend = Aer.get_backend('qasm_simulator')
job = execute(qc, backend, shot=1024)

result = job.result()
plot_histogram(result.get_counts())

## 実行(実機)
backend = provider.get_backend('ibmq_manila')
job = execute(qc, backend, shots=1024)

result = job.result()
plot_histogram(result.get_counts())

ということで、第一量子ビットの確率はRe⟨ψ|X^|ψ⟩$\mathrm{R}\mathrm{e}⟨\psi |\hat{X}|\psi ⟩$に従う分布になりました。

まとめ

• アダマールテストを用いると、制御ユニタリーに用いたユニタリー演算子の固有値がわかる
• アダマールテストによって第一量子ビットの確率分布がユニタリー演算子の固有値に応じた位相の影響を受ける、「位相キックバック」が起こる
• かけるユニタリーと入力する第二量子ビットの状態によって第一量子ビットの確率分布を制御できる

参考

• エルミート行列 - Wikipedia
• 行列指数関数 - Wikipedia

量子コンピューティング 基本アルゴリズムから量子機械学習まで [ 情報処理学会出版委員会 ]

密度演算子・密度行列

お久しぶりです。
色々忙しくて前回から随分間が空いてしまいました。

当初アダマールテストとかのアルゴリズムの記事を書こうとしましたが、それらを勉強したり説明するうえで、部分系や測定や確率分布の話が避けて通れませんでした。
先に密度演算子について解説して、部分系の取り扱い方法を整理しようと思います。

対象とする読者

• 量子プログラミングに興味を持ち始めたばかりのエンジニア
• 量子計算に興味を持ち始めたばかりの高校生・大学生

など

前提とする知識

• ブラケットを使った線形代数がわかること
• 量子状態をブラケットないしBloch球で表現できること
• ユニタリー発展がわかること
• 射影測定がわかること

バージョン情報

• Python 3.9.13
• Qiskit 0.20.2

目次

定義

密度演算子・密度行列ρ$\rho$は以下で定義されます。

ρ=∑ipi|ψi⟩⟨ψi|${\displaystyle \rho =\sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|}$

ここで、|ψi⟩$|{\psi }_i⟩$は量子系のとりうる状態で、pi$p_i$はそれに対する確率(実数)です。

性質

トレースは1

トレースとは行列の対角成分だけの和を取ることで、以下となります。

tr(ρ)=tr(∑ipi|ψi⟩⟨ψi|)=∑ipitr(|ψi⟩⟨ψi|)=∑ipi=1$\begin{array}{rl}\mathrm{t}\mathrm{r}(\rho )& =\mathrm{t}\mathrm{r}\left(\sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|\right)\\ & =\sum _i{p}_i\mathrm{t}\mathrm{r}(|{\psi }_i⟩⟨{\psi }_i|)\\ & =\sum _i{p}_i\\ & =1\end{array}$

つまり、トレースは確率の和を示しています。

なお、3行目から4行目の計算において、

|ψi⟩=∑jαj|ψij⟩(∑j|αj|2=1)${\displaystyle |{\psi }_i⟩=\sum _j{\alpha }_j|{\psi }_{ij}⟩(\sum _j|{\alpha }_j{|}^2=1)}$

で|ψij⟩$|{\psi }_{ij}⟩$は正規直交基底のため、tr(|ψi⟩⟨ψi|)=∑j|αj|2=1$\mathrm{t}\mathrm{r}(|{\psi }_i⟩⟨{\psi }_i|)=\sum _j|{\alpha }_j{|}^2=1$となることを利用しています。

対角成分はゼロ以上

対角成分は確率を示すのでそれはそうなのですが、ちゃんと式の上で確かめられます。 ある量子状態を|ϕ⟩$|\varphi ⟩$とした場合以下となります。

⟨ϕ|ρ|ϕ⟩=⟨ϕ|(∑ipi|ψi⟩⟨ψi|)|ϕ⟩=∑ipi⟨ϕ|ψi⟩⟨ψi|ϕ⟩=∑ipi|⟨ϕ|ψi⟩|2≥0 (∵pi≥0,|⟨ϕ|ψi⟩|2≥0)$\begin{array}{rl}⟨\varphi |\rho |\varphi ⟩& =⟨\varphi |\left(\sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|\right)|\varphi ⟩\\ & =\sum _i{p}_i⟨\varphi |{\psi }_i⟩⟨{\psi }_i|\varphi ⟩\\ & =\sum _i{p}_i|⟨\varphi |{\psi }_i⟩{|}^2\\ & \ge 0(\because p_i\ge 0,|⟨\varphi |{\psi }_i⟩{|}^2\ge 0)\end{array}$

確率を示す物理量とも言えるので、その期待値がちゃんと確率になっているというのがわかります。

定理

密度演算子のユニタリー発展

状態|ψi⟩$|{\psi }_i⟩$がユニタリー発展してU^|ψi⟩$\hat{U}|{\psi }_i⟩$となるので

ρinit=∑ipi|ψi⟩⟨ψi|−→U^∑ipiU^|ψi⟩⟨ψi|U^†=U^(∑ipi|ψi⟩⟨ψi|)U^†=U^ρinitU^†≡ρfinal$\begin{array}{rl}{\rho }_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}=\sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|& \stackrel{\hat{U}}{\to }\sum _i{p}_i\hat{U}|{\psi }_i⟩⟨{\psi }_i|{\hat{U}}^{†}\\ & =\hat{U}\left(\sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|\right){\hat{U}}^{†}\\ & =\hat{U}{\rho }_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}{\hat{U}}^{†}\equiv {\rho }_{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}\end{array}$

ということで、元の状態ρinit${\rho }_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}$とユニタリー発展U^$\hat{U}$が既知であれば計算から求められます。

測定後の密度演算子

測定値m$m$が得られる測定演算子M^m${\hat{M}}_m$(射影演算子とは限らない)を使って量子状態|ψ⟩$|\psi ⟩$を測定した時に、測定値m$m$が得られる確率p(m)$p(m)$はこちらで説明したように、以下の式で求まります。

p(m)=⟨ψ|M^†mM^m|ψ⟩${\displaystyle p(m)=⟨\psi |{\hat{M}}_m^{†}{\hat{M}}_m|\psi ⟩}$

i$i$で示されるいくつかの量子状態|ψi⟩$|{\psi }_i⟩$がある場合も同様に、そのi$i$に対する条件付き確率は以下で求まります。

p(m|i)=⟨ψi|M^†mM^m|ψi⟩=tr(M^†mM^m|ψi⟩⟨ψi|)$\begin{array}{rl}p(m|i)& =⟨{\psi }_i|{\hat{M}}_m^{†}{\hat{M}}_m|{\psi }_i⟩\\ & =\mathrm{t}\mathrm{r}({\hat{M}}_m^{†}{\hat{M}}_m|{\psi }_i⟩⟨{\psi }_i|)\end{array}$

ここでの式変形では次の行列の性質を利用しています。

tr(A^|ψ⟩⟨ψ|)=⟨ψ|A^|ψ⟩${\displaystyle \mathrm{t}\mathrm{r}(\hat{A}|\psi ⟩⟨\psi |)=⟨\psi |\hat{A}|\psi ⟩}$

さまざまなi$i$についてのp(m|i)$p(m|i)$を合計することでp(m)$p(m)$が得られます。 (確率論の言葉では「周辺化」と言います。)

p(m)=∑ip(m|i)⋅pi=∑ipi⋅tr(M^†mM^m|ψi⟩⟨ψj|)=tr(M^†mM^m∑ipi|ψi⟩⟨ψi|)$\begin{array}{rl}p(m)& =\sum _{i}p(m|i)\cdot p_i\\ & =\sum _i{p}_i\cdot \mathrm{t}\mathrm{r}({\hat{M}}_m^{†}{\hat{M}}_m|{\psi }_i⟩⟨{\psi }_j|)\\ & =\mathrm{t}\mathrm{r}\left({\hat{M}}_m^{†}{\hat{M}}_m\sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|\right)\end{array}$

以上を用いて測定後の密度演算子を求めるが、そのために測定後の量子状態を確認するとこちらで説明したように

|ψmi⟩=M^m|ψi⟩⟨ψi|M^†mM^m|ψi⟩−−−−−−−−−−−−√${\displaystyle |{\psi }_i^m⟩=\frac{{\hat{M}}_m|{\psi }_i⟩}{\sqrt{⟨{\psi }_i|{\hat{M}}_m^{†}{\hat{M}}_m|{\psi }_i⟩}}}$

となることを用いて、測定後の密度演算子ρm${\rho }_m$は以下になる

ρm=∑ip(i|m)|ψmi⟩⟨ψmi|=∑ip(i|m)M^m|ψi⟩⟨ψi|M^†m⟨ψi|M^†mM^m|ψi⟩$\begin{array}{rl}{\rho }_m& =\sum _{i}p(i|m)|{\psi }_i^m⟩⟨{\psi }_i^m|\\ & =\sum _{i}p(i|m)\frac{{\hat{M}}_m|{\psi }_i⟩⟨{\psi }_i|{\hat{M}}_m^{†}}{⟨{\psi }_i|{\hat{M}}_m^{†}{\hat{M}}_m|{\psi }_i⟩}\end{array}$

ここで、p(i|m)$p(i|m)$については確率の乗法法則

p(i|m)=p(m,i)p(m)=p(m|i)⋅pip(m)${\displaystyle p(i|m)=\frac{p(m,i)}{p(m)}=\frac{p(m|i)\cdot p_i}{p(m)}}$

より

ρm=∑ip(m|i)⋅pip(m)⋅M^m|ψi⟩⟨ψi|M^†m⟨ψi|M^†mM^m|ψi⟩=∑i⟨ψi|M^†mM^m|ψi⟩tr(M^†mM^mρ)⋅pi⋅M^m|ψi⟩⟨ψi|M^†m⟨ψi|M^†mM^m|ψi⟩=∑ipi⋅M^m|ψi⟩⟨ψi|M^†mtr(M^†mM^mρ)=M^m∑ipi|ψi⟩⟨ψi|M^†mtr(M^†mM^mρ)=M^mρM^†mtr(M^†mM^m)$\begin{array}{rl}{\rho }_m& =\sum _i\frac{p(m|i)\cdot p_i}{p(m)}\cdot \frac{{\hat{M}}_m|{\psi }_i⟩⟨{\psi }_i|{\hat{M}}_m^{†}}{⟨{\psi }_i|{\hat{M}}_m^{†}{\hat{M}}_m|{\psi }_i⟩}\\ & =\sum _i\frac{⟨{\psi }_i|{\hat{M}}_m^{†}{\hat{M}}_m|{\psi }_i⟩}{\mathrm{t}\mathrm{r}\left({\hat{M}}_m^{†}{\hat{M}}_m\rho \right)}\cdot p_i\cdot \frac{{\hat{M}}_m|{\psi }_i⟩⟨{\psi }_i|{\hat{M}}_m^{†}}{⟨{\psi }_i|{\hat{M}}_m^{†}{\hat{M}}_m|{\psi }_i⟩}\\ & =\sum _i{p}_i\cdot \frac{{\hat{M}}_m|{\psi }_i⟩⟨{\psi }_i|{\hat{M}}_m^{†}}{\mathrm{t}\mathrm{r}\left({\hat{M}}_m^{†}{\hat{M}}_m\rho \right)}\\ & =\frac{{\hat{M}}_m\sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|{\hat{M}}_m^{†}}{\mathrm{t}\mathrm{r}\left({\hat{M}}_m^{†}{\hat{M}}_m\rho \right)}\\ & =\frac{{\hat{M}}_m\rho {\hat{M}}_m^{†}}{\mathrm{t}\mathrm{r}\left({\hat{M}}_m^{†}{\hat{M}}_m\right)}\end{array}$

確率的に分布する量子状態

ある状態ρi=∑jpj|ψj⟩⟨ψj|${\rho }_i=\sum _j{p}_j|{\psi }_j⟩⟨{\psi }_j|$が確率pi$p_i$で存在するとき、全体の密度演算子は

ρ=∑i,jpipj|ψj⟩⟨ψj|=∑ipiρi$\begin{array}{rl}\rho & =\sum _{i,j}{p}_i{p}_j|{\psi }_j⟩⟨{\psi }_j|\\ & =\sum _i{p}_i{\rho }_i\end{array}$

なぜ密度演算子が必要なのか

量子状態(波動関数・状態ベクトル)はまさに状態そのものを示していますがその係数は確率振幅であるので測定に直接結びつきません。 一方で、密度演算子は確率そのものでそのもので表されているため、測定と直接結びつけることができます。

また、行列成分の現れ方でノイズの影響下にあるのか、エンタングルメント状態なのかなどを判断できます。 測定結果を理論と照らし合わせるために非常に重要です。

一方で、元々使っていた量子状態(波動関数・状態ベクトル)はどういう時に使うかというと位相を解析したい時ですかね、、。
詳しい方いらしゃったらコメントいただけますと嬉しいです。

密度演算子の具体例

よくある例ですけど、せっかくなのでQiskit使って測定まで含めて確認してみます。

純粋状態

測定対象となる状態は

|ψ⟩=12–√(|0⟩+|1⟩)${\displaystyle |\psi ⟩=\frac{1}{\sqrt{2}}(|0⟩+|1⟩)}$

です。 この1つの状態が単一で存在します。

|0⟩≡(10)${\displaystyle |0⟩\equiv \left(\begin{array}{c}1\\ 0\end{array}\right)}$

|1⟩≡(01)${\displaystyle |1⟩\equiv \left(\begin{array}{c}0\\ 1\end{array}\right)}$

とすると、密度演算子は

ρ=|ψ⟩⟨ψ|=12(|0⟩⟨0|+|0⟩⟨1|+|1⟩⟨0|+|1⟩⟨1|)=12(1111)$\begin{array}{rl}\rho & =|\psi ⟩⟨\psi |\\ & =\frac{1}{2}(|0⟩⟨0|+|0⟩⟨1|+|1⟩⟨0|+|1⟩⟨1|)\\ & =\frac{1}{2}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)\end{array}$

対角要素が各状態の確率に等しく、|0⟩$|0⟩$と|1⟩$|1⟩$の確率はそれぞれ1/2$1/2$です。

また、

tr(ρ2)=14(2222)=1$\begin{array}{rl}\mathrm{t}\mathrm{r}\left({\rho }^2\right)& =\frac{1}{4}\left(\begin{array}{cc}2& 2\\ 2& 2\end{array}\right)\\ & =1\end{array}$

となるのですが、このように

tr(ρ2)=1${\displaystyle \mathrm{t}\mathrm{r}\left({\rho }^2\right)=1}$

となる状態を純粋状態(Pure State)と言います。

これをZ測定します。 この時の量子回路は以下となります。

図1. 重ね合わせ状態のZ測定

この測定結果は下図です。 1024回測定しました。

図2. 図1の回路の測定結果

なお、この時のソースです。

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(0)
circuit.barrier()

circuit.measure(qr, cr)

backend = Aer.get_backend('qasm_simulator')
job = execute(circuit, backend=backend, shots=1024)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

混合状態

測定対象となる状態は

|ψ0⟩=|0⟩${\displaystyle |{\psi }_0⟩=|0⟩}$

か

|ψ1⟩=|1⟩${\displaystyle |{\psi }_1⟩=|1⟩}$

がそれぞれ確率50%です。 従って密度演算子は

ρ=12|ψ0⟩⟨ψ0|+12|ψ1⟩⟨ψ1|=12(|0⟩⟨0|+|1⟩⟨1|)=12(1001)$\begin{array}{rl}\rho & =\frac{1}{2}|{\psi }_0⟩⟨{\psi }_0|+\frac{1}{2}|{\psi }_1⟩⟨{\psi }_1|\\ & =\frac{1}{2}(|0⟩⟨0|+|1⟩⟨1|)\\ & =\frac{1}{2}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\end{array}$

また、

tr(ρ2)=14(1001)=12$\begin{array}{rl}\mathrm{t}\mathrm{r}\left({\rho }^2\right)& =\frac{1}{4}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\\ & =\frac{1}{2}\end{array}$

となるのですが、このように

tr(ρ2)<1${\displaystyle \mathrm{t}\mathrm{r}\left({\rho }^2\right)<1}$

となる状態を混合状態(Mixed State)と言います。

これをZ測定します。
こちらも対角要素が各状態の確率に等しく、|0⟩$|0⟩$と|1⟩$|1⟩$の確率はそれぞれ1/2$1/2$です。

ただし、量子回路は以下のどちらかを50%のランダムで実行することになります。

図3. 混合状態をZ測定する量子回路

量子ビットにXゲートがかかるかかからないかを乱数で50%の確率で決めます。

この測定結果は下図です。

図4. 図3の量子回路の測定結果

なお、この時のソースです。(※実行するには時間がかかります)

gates = [1 if val >= 0.5 else 0 for val in np.random.rand(1024)]

results = list()
for idx, gate in enumerate(gates):
    qr = QuantumRegister(1)
    cr = ClassicalRegister(1)
    circuit = QuantumCircuit(qr, cr)

    if gate:
        circuit.x(0)
    circuit.barrier()

    circuit.measure(qr, cr)
    
    backend = Aer.get_backend('qasm_simulator')
    job = execute(circuit, backend=backend, shots=1)
    job_monitor(job, interval = 2)

    result = job.result()
    count = result.get_counts()
    
    results.append(count)

counts_classic = {'0': 0, '1': 0}
for result in results:
    for bit, count in result.items():
        counts_classic[bit] += count

plot_histogram(counts_classic)

重ね合わせ状態と混合状態の見分け方

重ね合わせ状態も古典的混合も似たような結果になりました。
数式だとはっきり違いがわかりますが、測定だと違いがわかりません。
測定対象の量子系が重ね合わせ状態か混合状態かどちらなのか予めわからない場合、またより一般的に対象の系がどのような理論やモデルに基づくのかがわからない場合にどうやって測定から見分けたらいいでしょうか？
これが密度演算子で理解できます。

先ほどの例ではZ測定するつもりで、密度演算子もZ基底で記述しました。
これをX測定してみます。

重ね合わせ状態

測定対象となる状態は

|ψ⟩=12–√(|0⟩+|1⟩)≡|+⟩${\displaystyle |\psi ⟩=\frac{1}{\sqrt{2}}(|0⟩+|1⟩)\equiv |+⟩}$

です。 この1つの状態が単一で存在します。

|+⟩≡(10)${\displaystyle |+⟩\equiv \left(\begin{array}{c}1\\ 0\end{array}\right)}$

|−⟩≡(01)${\displaystyle |-⟩\equiv \left(\begin{array}{c}0\\ 1\end{array}\right)}$

とすると、密度演算子は

ρ=|ψ⟩⟨ψ|=|+⟩⟨+|=(1000)$\begin{array}{rl}\rho & =|\psi ⟩⟨\psi |\\ & =|+⟩⟨+|\\ & =\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)\end{array}$

|+⟩$|+⟩$状態である確率100%です。

今度の量子回路はこちらになります。

図5. 重ね合わせ状態のX測定

この測定結果は下図です。 1024回測定しました。

図6. 図5の量子回路の測定結果

予定通り100%で観測されました。(アダマールゲートで基底変換しているので、この図の|0⟩$|0⟩$が|+⟩$|+⟩$に相当します。)

この時のソースは以下です。

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(0)
circuit.barrier()

circuit.h(0)
circuit.measure(qr, cr)

backend = Aer.get_backend('qasm_simulator')
job = execute(circuit, backend=backend, shots=1024)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()
plot_histogram(count)

混合状態

一方で混合状態をX基底で記述すると以下になります。

|ψ0⟩=|0⟩=12–√(|+⟩+|−⟩)${\displaystyle |{\psi }_0⟩=|0⟩=\frac{1}{\sqrt{2}}(|+⟩+|-⟩)}$

か

|ψ1⟩=|1⟩=12–√(|+⟩−|−⟩)${\displaystyle |{\psi }_1⟩=|1⟩=\frac{1}{\sqrt{2}}(|+⟩-|-⟩)}$

がそれぞれ確率50%です。 従って密度演算子は

ρ=12|ψ0⟩⟨ψ0|+12|ψ1⟩⟨ψ1|=12{12(|+⟩+|−⟩)(⟨+|+⟨−|)+12(|+⟩−|−⟩)(⟨+|−⟨−|)}=14{(1000)+(0010)+(0100)+(0001)}+14{(1000)−(0010)−(0100)+(0001)}=12(1001)$\begin{array}{rl}\rho & =\frac{1}{2}|{\psi }_0⟩⟨{\psi }_0|+\frac{1}{2}|{\psi }_1⟩⟨{\psi }_1|\\ & =\frac{1}{2}\{\frac{1}{2}(|+⟩+|-⟩)(⟨+|+⟨-|)+\frac{1}{2}(|+⟩-|-⟩)(⟨+|-⟨-|)\}\\ & =\frac{1}{4}\{\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)+\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)+\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right)+\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)\}\\ & +\frac{1}{4}\{\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)-\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)-\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right)+\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)\}\\ & =\frac{1}{2}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\end{array}$

ということで、Z基底で記述した時と同じく、|+⟩$|+⟩$と|−⟩$|-⟩$状態が50%ずつの確率となります。

この時の量子回路は以下のどちらかとなります。

図7. |0⟩$|0⟩$と|1⟩$|1⟩$の混合状態をX測定する場合

この測定結果は以下となります。

図8. 図7の量子回路の測定結果

ということで、密度演算子から考えた通り、|+⟩$|+⟩$と|−⟩$|-⟩$状態が約50%ずつで観測されました。

このように未知の量子系であっても、適切な測定方法を用いることで得られる確率分布と理論を照らし合わせて、モデル・量子系を推定することができます。

この時のソースは以下です。

gates = [1 if val >= 0.5 else 0 for val in np.random.rand(1024)]

results = list()
for idx, gate in enumerate(gates):
    qr = QuantumRegister(1)
    cr = ClassicalRegister(1)
    circuit = QuantumCircuit(qr, cr)

    if gate:
        circuit.x(0)
    circuit.barrier()

    circuit.h(0)
    circuit.measure(qr, cr)

    if idx < 10:
        fig, ax = plt.subplots(1, 1)
        circuit.draw(output='mpl', ax=ax)
        plt.show()
        plt.close()
    
    backend = Aer.get_backend('qasm_simulator')
    job = execute(circuit, backend=backend, shots=1)
    job_monitor(job, interval = 2)

    result = job.result()
    count = result.get_counts()
    
    results.append(count)
    
counts_classic = {'0': 0, '1': 0}
for result in results:
    for bit, count in result.items():
        counts_classic[bit] += count

plot_histogram(counts_classic)

量子状態トモグラフィー

上記では簡単な例で実験結果と理論を比較して状態を推定するということを紹介しました。 実際にはより多量子ビット系になりますが要領は同じで、さまざまな測定から密度演算子の要素を推定し、理論と照らし合わせて量子状態を推定します。 このような手法を「量子状態トモグラフィー」と言います。 「トモグラフィー」は物理学に限らず一般的な言葉で、日常で触れるところではMRIの画像がそうです。 磁気共鳴の空間分布を解析することによって、脳や臓器の構造と状態を指定しそして可視化しています。

参考までに量子状態トモグラフィーについて以下の記事を紹介します。

qiita.com

arxiv.org

いずれこのブログでも解説してみたいですが、いつになることやら、、。

部分系・部分トレース

量子系Aと量子系Bの複合系がある時の全体系の密度演算子をρAB${\rho }^{AB}$とします。
例えば、二つの量子ビットからなる一つの量子系でも良いです。

この時、二つの量子ビットをどちらも測定すれば|00⟩,|01⟩,|10⟩,|11⟩$|00⟩,|01⟩,|10⟩,|11⟩$の4状態の確率は密度演算子の対角成分に示された通りになるでしょう。

では、量子ビットAだけを測定した場合、量子ビットAの状態 |0⟩A,|1⟩A$|0{⟩}_A,|1{⟩}_A$が測定される確率、またAの密度演算子ρA${\rho }^A$はどのようになるでしょうか？という問題です。

結論としては次の式です。

ρA=trB(ρAB)${\displaystyle {\rho }^A={\mathrm{t}\mathrm{r}}_B\left({\rho }^{AB}\right)}$

Bについてだけトレースを取ります。(トレースアウトする、とよく言います。)

実際にはより具体的にρAB${\rho }^{AB}$の成分を計算する際に以下のように計算します。

trB(|ai⟩⟨aj|⊗|bk⟩⟨bl|)=|ai⟩⟨aj|tr(|bk⟩⟨bl|)${\displaystyle {\mathrm{t}\mathrm{r}}_B(|a_i⟩⟨a_j|\otimes |b_k⟩⟨b_l|)=|a_i⟩⟨a_j|\mathrm{t}\mathrm{r}(|b_k⟩⟨b_l|)}$

もしくはρAB=ρA⊗ρB${\rho }^{AB}={\rho }^A\otimes {\rho }^B$のように簡単に表せる場合は以下で計算できます。

ρA=trB(ρA⊗ρB)=ρAtr(ρB)=ρA(∵(tr)(ρB)=1)$\begin{array}{rl}{\rho }^A& ={\mathrm{t}\mathrm{r}}_B({\rho }^A\otimes {\rho }^B)\\ & ={\rho }^A\mathrm{t}\mathrm{r}\left({\rho }^B\right)\\ & ={\rho }^A(\because (tr)\left({\rho }^B\right)=1)\end{array}$

Bell状態を使った具体例

Bell状態の一つ

|ψ⟩=12–√(|00⟩+|11⟩)${\displaystyle |\psi ⟩=\frac{1}{\sqrt{2}}(|00⟩+|11⟩)}$

の密度演算子は以下です。

ρAB=12(|00⟩+|11⟩)(⟨00|+⟨11|)=12(|00⟩⟨00|+|11⟩⟨00|+|00⟩⟨11|+|11⟩⟨11|)$\begin{array}{rl}{\rho }^{AB}& =\frac{1}{2}(|00⟩+|11⟩)(⟨00|+⟨11|)\\ & =\frac{1}{2}(|00⟩⟨00|+|11⟩⟨00|+|00⟩⟨11|+|11⟩⟨11|)\end{array}$

これを実際にB$B$についてトレースアウトすると以下です。

ρA=trB(ρAB)=12{trB(|00⟩⟨00|)+trB(|11⟩⟨00|)+trB(|00⟩⟨11|)+trB(|11⟩⟨11|)}=12{trB(|0⟩⟨0|⊗|0⟩⟨0|)+trB(|1⟩⟨0|⊗|1⟩⟨0|)+trB(|0⟩⟨1|⊗|0⟩⟨1|)+trB(|1⟩⟨1|⊗|1⟩⟨1|)}=12{|0⟩⟨0|tr(|0⟩⟨0|)+|1⟩⟨0|tr(|1⟩⟨0|)+|0⟩⟨1|tr(|0⟩⟨1|)+|1⟩⟨1|tr(|1⟩⟨1|)}=12(|0⟩⟨0|⟨0|0⟩+|1⟩⟨0|⟨0|1⟩+|0⟩⟨1|⟨1|0⟩+|1⟩⟨1|⟨1|1⟩)=12(|0⟩⟨0|+|1⟩⟨1|)$\begin{array}{rl}{\rho }^A& ={\mathrm{t}\mathrm{r}}_B\left({\rho }^{AB}\right)\\ & =\frac{1}{2}\{{\mathrm{t}\mathrm{r}}_B(|00⟩⟨00|)+{\mathrm{t}\mathrm{r}}_B(|11⟩⟨00|)+{\mathrm{t}\mathrm{r}}_B(|00⟩⟨11|)+{\mathrm{t}\mathrm{r}}_B(|11⟩⟨11|)\}\\ & =\frac{1}{2}\{{\mathrm{t}\mathrm{r}}_B(|0⟩⟨0|\otimes |0⟩⟨0|)+{\mathrm{t}\mathrm{r}}_B(|1⟩⟨0|\otimes |1⟩⟨0|)+{\mathrm{t}\mathrm{r}}_B(|0⟩⟨1|\otimes |0⟩⟨1|)+{\mathrm{t}\mathrm{r}}_B(|1⟩⟨1|\otimes |1⟩⟨1|)\}\\ & =\frac{1}{2}\{|0⟩⟨0|\mathrm{t}\mathrm{r}(|0⟩⟨0|)+|1⟩⟨0|\mathrm{t}\mathrm{r}(|1⟩⟨0|)+|0⟩⟨1|\mathrm{t}\mathrm{r}(|0⟩⟨1|)+|1⟩⟨1|\mathrm{t}\mathrm{r}(|1⟩⟨1|)\}\\ & =\frac{1}{2}(|0⟩⟨0|⟨0|0⟩+|1⟩⟨0|⟨0|1⟩+|0⟩⟨1|⟨1|0⟩+|1⟩⟨1|⟨1|1⟩)\\ & =\frac{1}{2}(|0⟩⟨0|+|1⟩⟨1|)\end{array}$

ということで、量子ビットAは|0⟩$|0⟩$と|1⟩$|1⟩$がそれぞれ確率1/2$1/2$で観測されます。 また

tr{(ρA)2}=12${\displaystyle \mathrm{t}\mathrm{r}\left\{{\left({\rho }^A\right)}^2\right\}=\frac{1}{2}}$

から、この状態は混合状態だということがわかりました。

このBell状態を使った例は計算結果と実際に測定される状態のイメージがつきやすいわかりやすい例でしたが、この部分トレースという方法は一般の自明でない状態についての確率を予め計算から見積もりたい時に強力なツールとなります。

まとめ

今回は密度演算子の基礎について説明しました。

• 密度演算子は各状態の確率を示す演算子
• 密度演算子のトレース(全状態の確率の和)は1
• 密度演算子の対角成分は非負
• 密度演算子の自乗のトレースを計算することでその状態が純粋状態か混合状態かを見分けることが可能
• 部分的に量子ビットを測定する場合の密度演算子は測定しない量子ビットの分を部分トレースすることで計算可能

次回はこれらを用いてアダマールテストについて解説する予定です。

また、今年、2022年4月よりフリーランスのデータサイエンティストとして活動することになりました。
量子アプリケーション開発というものはまだ市場がほとんどありませんので、主な事業は通常のシステム開発や人工知能アプリケーション開発、またデータ分析の業務委託です。
このブログでの機械学習や量子計算・量子コンピュータについての解説や情報発信は引き続き行いますので、今後とも読んでいただけますと嬉しいです。

私の事業用のホームページがこちらになります。

sakuma.wraptas.site

参考文献

• 量子情報科学ウィンタースクール2010 2日目：量子情報理論の基礎Ｉの14ページ以降
• エンタングル状態
• 密度演算子 - 純粋状態・混合状態 - 量子コンピュータ入門 - 物理とか
• 量子力学的測定の定式化と部分系・部分トレース - 物理とかの4.部分トレース
• 量子計算量理論入門の2.4混合状態(36ページ)
• 量子力学II講義ノートの1.5 密度演算子(16ページ)

物理量の期待値

前回の射影測定で状態がわかると物理量(演算子)の値がわかります。
また、一般に状態は重ね合わせのため、多数回測定したときには期待値が得られますのでその計算方法を紹介します。

量子力学では頻繁に登場する概念ですが、量子計算は状態を測定できればいいので物理量を意識することはほとんどありません。
ただ、VQEなど一部のアルゴリズムでは物理量(といってもσx,σy,σz${\sigma }_x,{\sigma }_y,{\sigma }_z$)の期待値が必要になる場合ばありますので量子コンピュータを用いた場合の期待値計算もご紹介します。

対象とする読者

• 量子プログラミングに興味を持ち始めたばかりのエンジニア
• 量子計算に興味を持ち始めたばかりの高校生・大学生 など

対象でない読者

• 量子情報処理の研究室ご出身の方
  など

※もちろん読んでいただける分には嬉しいですが、
特に前半部分は電子回路の授業で扱われていたり教科書にも載っている内容となりますので
いまさらと感じさせてしまうかもしれません。

前提とする知識

• 簡単な線形代数が分かること
• 固有値・固有ベクトルがわかること

バージョン情報

• Python 3.7.3
• Qiskit 0.23.1

目次

物理量の演算子

物理量

物理量とは可観測量(Observable)とも呼ばれ、量子力学で実際に何らかの値として測定されうる量のことです。
例えば下記があります。(太字はベクトル量です。)

• 位置 x$\mathbf{x}$
• 速度 v$\mathbf{v}$ (もしくは運動量 p$\mathbf{p}$ )
• 粒子数 n$n$
• 位相 θ$\theta$
• エネルギー H$H$
• スピン σ$\sigma$

量子計算で一番馴染みがありそうなのはエネルギーだと思います。
エネルギーの演算子H$H$がハミルトニアンで、エネルギー固有値がその状態が持つエネルギーで、その一番小さいエネルギーを持つ状態を|0⟩$|0⟩$とし、その1つ上のエネルギーを持つ状態を|1⟩$|1⟩$としたのが超伝導量子ビットです。

下図は初期の超伝導磁束量子ビットでRabi振動(量子状態のマイクロ波制御)を観測した論文に掲載されている図です。
超伝導磁束量子ビットのエネルギーは超伝導リングを貫く磁束(磁場)で変化します(右図A)。
一番下の曲線が基底状態(量子計算でいうところの|0⟩$|0⟩$状態)のエネルギーの変化を、またその一つ上の曲線が第一励起状態(|1⟩$|1⟩$)のエネルギーの変化の理論計算値を示しています。
図Bは飛ばして図Cが基底状態に置いた超伝導量子ビットにある磁場(横軸)をかけたあと、様々な周波数のマイクロ波を照射した結果、吸収された周波数(縦軸)を示します。
図Cの実線が基底状態と第一励起状態のエネルギーの理論値の差で、点が測定された周波数です。
マイクロ波フォトンのエネルギーはhν$h\nu$ (h$h$はPlank定数で、ν$\nu$が周波数)で計算されるので、周波数を測定することでその状態における系のエネルギーを測定できたことに相当します。

図1. 超伝導磁束量子ビットのエネルギー測定

引用した論文はこちら arxiv.org

演算子の定義

物理量の概念や物理学における実体が以上で理解いただけたとして、数学上の物理量A^$\hat{A}$の定義は一般に下記です。

A^=∑iai|ai⟩⟨ai|${\displaystyle \hat{A}=\sum _i{a}_i|a_i⟩⟨a_i|}$

ここで|ai⟩$|a_i⟩$は物理量A^$\hat{A}$の固有ベクトル(状態ベクトル)でai$a_i$がその固有ベクトルに対応する固有値です。

量子計算では状態ベクトルは基本的に|0⟩,|1⟩$|0⟩,|1⟩$ (とその組み合わせ)しか用いません。
これはスピンが上(|0⟩$|0⟩$)を向いているか、下(|1⟩$|1⟩$)を向いているかに相当します。
パウリ行列σz^$\widehat{{\sigma }_z}$の固有ベクトルを量子計算で使う|0⟩,|1⟩$|0⟩,|1⟩$に対応しており、下記のように表せます。

σz^=|0⟩⟨0|−|1⟩⟨1|${\displaystyle \widehat{{\sigma }_z}=|0⟩⟨0|-|1⟩⟨1|}$

量子ビットが|0⟩$|0⟩$なら「1」でZ軸上向き、そして|1⟩$|1⟩$なら「-1」でZ軸下向きを表します。

測定値に中間の値はなく、同一の状態|ψ⟩=α|0⟩+β|1⟩$|\psi ⟩=\alpha |0⟩+\beta |1⟩$として生成した量子ビットを複数回測定した場合に、σx${\sigma }_x$の測定値は|α|2$|\alpha {|}^2$の確率で「1」となり|β|2$|\beta {|}^2$の確率で「-1」となります。
また後述しますが、この時の期待値は通常の統計における期待値と同様に±1にそれぞれの重みを乗じたものの和になります。

パウリ行列

代表的な物理量としてパウリ行列をブラケット記法でご紹介します。

σz=|0⟩⟨0|−|1⟩⟨1|=(100−1)(Z基底の場合)$\begin{array}{rl}{\sigma }_z& =|0⟩⟨0|-|1⟩⟨1|=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)(Z基底の場合)\end{array}$

σx=(0110)=|0⟩⟨1|+|1⟩⟨0|(Z基底の場合)=|+⟩⟨+|−|−⟩⟨−|(X基底の場合)$\begin{array}{rl}{\sigma }_x& =\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)=|0⟩⟨1|+|1⟩⟨0|(Z基底の場合)\\ & =|+⟩⟨+|-|-⟩⟨-|(X基底の場合)\end{array}$

ここで|±⟩=(1/2–√)(|0⟩±|1⟩)$|\pm ⟩=(1/\sqrt{2})(|0⟩\pm |1⟩)$

σy=(0−ii0)=i|0⟩⟨1|−i|1⟩⟨0|(Z基底の場合)=|π2⟩⟨π2|−|−π2⟩⟨−π2|(Y基底の場合)$\begin{array}{rl}{\sigma }_y& =\left(\begin{array}{cc}0& i\\ -i& 0\end{array}\right)=i|0⟩⟨1|-i|1⟩⟨0|(Z基底の場合)\\ & =|\frac{\pi }{2}⟩⟨\frac{\pi }{2}|-|-\frac{\pi }{2}⟩⟨-\frac{\pi }{2}|(Y基底の場合)\end{array}$

ここで|±π/2⟩=(1/2–√)(|0⟩+e±iπ/2|1⟩)$|\pm \pi /2⟩=(1/\sqrt{2})(|0⟩+e^{\pm i\pi /2}|1⟩)$

それぞれのパウリ行列に対応する基底に変換すると、いずれについても±1$\pm 1$になるのが面白いところですね。
それぞれの軸について量子ビットがどちらを向いているかが確かに表されています。

余談ですが、物理学ではパウリ行列単独で物理量となることはなく、スピン(s^=(1/2)σ^$\hat{\mathbf{s}}=(1/2)\hat{\sigma }$)が物理量となります。

物理量の期待値の計算方法

物理量A^$\hat{A}$の期待値⟨A^⟩$⟨\hat{A}⟩$は任意の量子状態|ψ⟩$|\psi ⟩$を用いて

⟨A^⟩=⟨ψ|A^|ψ⟩${\displaystyle ⟨\hat{A}⟩=⟨\psi |\hat{A}|\psi ⟩}$

例えば、σz${\sigma }_z$について、状態が|ψ⟩=α|0⟩+β|1⟩$|\psi ⟩=\alpha |0⟩+\beta |1⟩$の場合における期待値⟨σz⟩$⟨{\sigma }_z⟩$は下記です。

⟨σz⟩=(α∗⟨0|+β∗⟨1|)(|0⟩⟨0|−|1⟩⟨1|)(α|0⟩+β|1⟩)=|α|2−|β|2$\begin{array}{rl}⟨{\sigma }_z⟩& =({\alpha }^{\ast }⟨0|+{\beta }^{\ast }⟨1|)(|0⟩⟨0|-|1⟩⟨1|)(\alpha |0⟩+\beta |1⟩)\\ & =|\alpha {|}^2-|\beta {|}^2\end{array}$

期待値計算の例題

単一の物理量

⟨σx⟩$⟨{\sigma }_x⟩$

• 量子状態が|ψ⟩=|0⟩$|\psi ⟩=|0⟩$の場合

⟨σx⟩=⟨0|(|+⟩⟨+|−|−⟩⟨−|)|0⟩=12(⟨+|+⟨−|)(|+⟩⟨+|−|−⟩⟨−|)(|+⟩+|−⟩)=0$\begin{array}{rl}⟨{\sigma }_x⟩& =⟨0|(|+⟩⟨+|-|-⟩⟨-|)|0⟩\\ & =\frac{1}{2}(⟨+|+⟨-|)(|+⟩⟨+|-|-⟩⟨-|)(|+⟩+|-⟩)\\ & =0\end{array}$

量子コンピュータ(IBMQ)上では測定の基底変換はできないので等価になるように状態を回転させます。

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.barrier()

circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

plot_histogram(count)

図2. |0⟩$|0⟩$をX基底で測定した結果 (左) 量子回路 (右) 測定結果

ヒストグラム上では|0⟩,|1⟩$|0⟩,|1⟩$と表記されているが、|+⟩,|−⟩$|+⟩,|-⟩$の分布を観察しているのと等価である。

そして得られた分布から

num_sampling = sum([dist for dist in count.values()])
print(f'サンプリング数: {num_sampling}')
>> サンプリング数: 8192

observable = [1, -1]
exp_val = sum([(cnt / num_sampling) * val for cnt, val in zip(observable, count.values())])
print(f'期待値: {exp_val}')
>> 期待値: 0.01953125

となり、代数計算した場合とほぼ同じくゼロとなった。

• 量子状態が|ψ⟩=|±⟩$|\psi ⟩=|\pm ⟩$の場合

⟨σx⟩=⟨±|(|+⟩⟨+|−|−⟩⟨−|)|±⟩=±1$\begin{array}{rl}⟨{\sigma }_x⟩& =⟨\pm |(|+⟩⟨+|-|-⟩⟨-|)|\pm ⟩\\ & =\pm 1\end{array}$

で量子コンピュータ(IBMQ)上では先ほどと同じ計算過程によって

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.barrier()

circuit.h(qr)

circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

plot_histogram(count)

図3. |+⟩$|+⟩$ をX基底で測定した結果

num_sampling = sum([dist for dist in count.values()])
print(f'サンプリング数: {num_sampling}')
>> サンプリング数: 8192

observable = [1, -1]
exp_val = sum([(cnt / num_sampling) * val for cnt, val in zip(observable, count.values())])
print(f'期待値: {exp_val}')
>> 期待値: 0.818603515625

である。

物理量の積

量子ビットが2つの場合の⟨σ1zσ2y⟩$⟨{\sigma }_{1z}{\sigma }_{2y}⟩$
(1つ目の量子ビットのσz${\sigma }_z$と2つ目の量子ビットのσy${\sigma }_y$の積の期待値)

全体の物理量は2状態かける2状態の総当たり(直積)となり、

• 量子状態が|ψ⟩=|0π2⟩$|\psi ⟩=|0\frac{\pi }{2}⟩$の場合

⟨σ1zσ2y⟩=⟨0π2|(|0π2⟩⟨0π2|+|1−π2⟩⟨1−π2)|−|0−π2⟩⟨0−π2|−|1π2⟩⟨1π2|)|0π2⟩=1$\begin{array}{rl}⟨{\sigma }_{1z}{\sigma }_{2y}⟩& =⟨0\frac{\pi }{2}|(|0\frac{\pi }{2}⟩⟨0\frac{\pi }{2}|+|1\frac{-\pi }{2}⟩⟨1\frac{-\pi }{2})|-|0\frac{-\pi }{2}⟩⟨0\frac{-\pi }{2}|-|1\frac{\pi }{2}⟩⟨1\frac{\pi }{2}|)|0\frac{\pi }{2}⟩\\ & =1\end{array}$

で、量子コンピュータ(IBMQ)を用いた場合もこれまで通り下記です。

qr = QuantumRegister(2)
cr = ClassicalRegister(2)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr[1])
circuit.rz(np.pi / 2, qr[1])
circuit.barrier()

circuit.rz(-np.pi / 2, qr[1])
circuit.h(qr[1])
circuit.barrier()

circuit.measure(qr, cr)

circuit.draw(output='mpl')

backend = provider.get_backend('ibmq_quito')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

num_sampling = sum([dist for dist in count.values()])
print(f'サンプリング数: {num_sampling}')
>> サンプリング数: 8192

observable = [1, -1, -1, 1]
exp_val = sum([(cnt / num_sampling) * val for cnt, val in zip(observable, count.values())])
print(f'期待値: {exp_val}')
>> 期待値: 0.9560546875

図4. |ψ⟩=|0π2⟩$|\psi ⟩=|0\frac{\pi }{2}⟩$をそれぞれZ基底とY基底で測定した結果

• 量子状態が|ψ⟩=12√(|0π2⟩+|1−π2⟩$|\psi ⟩=\frac{1}{\sqrt{2}}(|0\frac{\pi }{2}⟩+|1\frac{-\pi }{2}⟩$の場合

⟨σ1zσ2y⟩=12(⟨0π2|+⟨1−π2|)(|0π2⟩⟨0π2|+|1−π2⟩⟨1−π2)|−|0−π2⟩⟨0−π2|−|1π2⟩⟨1π2|)(|0π2⟩+|1−π2⟩)=1$\begin{array}{rl}⟨{\sigma }_{1z}{\sigma }_{2y}⟩& =\frac{1}{2}(⟨0\frac{\pi }{2}|+⟨1\frac{-\pi }{2}|)(|0\frac{\pi }{2}⟩⟨0\frac{\pi }{2}|+|1\frac{-\pi }{2}⟩⟨1\frac{-\pi }{2})|-|0\frac{-\pi }{2}⟩⟨0\frac{-\pi }{2}|-|1\frac{\pi }{2}⟩⟨1\frac{\pi }{2}|)(|0\frac{\pi }{2}⟩+|1\frac{-\pi }{2}⟩)\\ & =1\end{array}$

で、量子コンピュータ(IBMQ)を用いた場合もこれまで通り下記です。

qr = QuantumRegister(2)
cr = ClassicalRegister(2)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.rz(np.pi / 2, qr[1])
circuit.cz(qr[0], qr[1])
circuit.barrier()

circuit.rz(-np.pi / 2, qr[1])
circuit.h(qr[1])
circuit.barrier()

circuit.measure(qr, cr)

circuit.draw(output='mpl')

backend = provider.get_backend('ibmq_quito')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

plot_histogram(count)

num_sampling = sum([dist for dist in count.values()])
print(f'サンプリング数: {num_sampling}')
>> サンプリング数: 8192

observable = [1, -1, -1, 1]
exp_val = sum([(cnt / num_sampling) * val for cnt, val in zip(observable, count.values())])
print(f'期待値: {exp_val}')
>> 期待値: 0.74365234375

図5. 12√(|0π2⟩+|1−π2⟩$\frac{1}{\sqrt{2}}(|0\frac{\pi }{2}⟩+|1\frac{-\pi }{2}⟩$を測定した結果

わかりにくいですが、この重ね合わせ状態はグラフの両端(1/2–√)(|00⟩+|11⟩)$(1/\sqrt{2})(|00⟩+|11⟩)$に対応します。

まとめ

• 物理量とは可観測量(Observable)とも呼ばれ、何らかの値として測定されうる量のこと
• 物理量の固有ベクトル(状態ベクトル)に対応する固有値が測定値になる
• 期待値はそれぞれの測定値を測定された状態の確率を乗じて求まる

量子コンピューティング 基本アルゴリズムから量子機械学習まで [ 情報処理学会出版委員会 ]

射影測定

変分分類器の解説をした時に計算結果を取り出すのに最後測定をするという話をしました。
量子機械学習に限った話ではありませんが、量子回路での計算(演算)結果を取り出すには適した量子ビットを適した方法で測定することが必要になります。
古典コンピュータではANDとかOR回路で演算した結果の電圧値を読み出すような行為に相当します。
量子ビットの場合はどの量子ビットをどのような方法で読み出すかで取り出せる情報が異なります。
将来的にこういうハードウェアを理解しないといけないようなことはソフトウェア工学的には無くなるべきと思いますが、現状ではまだこのような原始的なレベルの理解が必要なので勉強しておきましょう。

対象とする読者

• 量子プログラミングに興味を持ち始めたばかりのエンジニア
• 量子計算に興味を持ち始めたばかりの高校生・大学生
• 数学・論理演算・電子回路に慣れていない方
  など

前提とする知識

• 2～3次元のベクトル・行列を計算できること
• Dirac記法(ブラ⟨⋅|$⟨\cdot |$とケット|⋅⟩$|\cdot ⟩$を使った量子状態・演算子の表現)を用いた計算式を読めること
• 基本的な量子ゲート(X, Y, Z)が分かること

バージョン情報

• Python 3.9.5
• Qiskit 0.29.0

目次

測定理論

量子測定の仮定

量子力学にはいくつかの仮定(postulate)があり、その一つに測定についての項目があります。
様々な本で解説されており、本によって具体的な文は異なりますが、主旨は同一です。
ここではNielsen-Chuangの2.2.3節にある文を引用します。

Postulate3: Quantum measurements are described by a collection {Mm}$\left\{M_m\right\}$ of measurement operators. These are operators acting on the state space of the system being measured. The index m refers to the measurement outcomes that may occur in the experiment. If the state of the quantum system is |ψ⟩$|\psi ⟩$ immediately before the measurement then the probability that result m occurs is given by

p(m)=⟨ψ|M†mMm|ψ⟩,(2.92)$\begin{array}{}\text{(}2.92\text{)}& {\displaystyle p(m)=⟨\psi |M_m^{†}{M}_m|\psi ⟩,}\end{array}$

and the state of the system after the measurement is

Mm|ψ⟩⟨ψ|M†mMm|ψ⟩−−−−−−−−−−−√.(2.93)$\begin{array}{}\text{(}2.93\text{)}& {\displaystyle \frac{M_m|\psi ⟩}{\sqrt{⟨\psi |M_m^{†}{M}_m|\psi ⟩}}.}\end{array}$

The measurement operators satisfy the completeness equation,

∑mM†mMm=I.(2.94)$\begin{array}{}\text{(}2.94\text{)}& {\displaystyle \sum _m{M}_m^{†}{M}_m=I.}\end{array}$

手元に英語版しかなくすみません、、

これはコペンハーゲン解釈ではいわゆる「波束の収縮」に当たる文言ですが、その哲学的な話には立ち入らず、量子計算に用いることができる数学的なエッセンスだけ取り出して解釈します。

この仮定から量子状態とその測定について理解できることは次の3点です。

• ある測定値(に対応するインデックス)m$m$が得られる測定演算子(行列)をMm$M_m$とする
• 測定前の量子状態|ψ⟩$|\psi ⟩$において測定値m$m$が得られる確率p(m)$p(m)$は(2.92)式で得られる
• 測定後の量子状態は(2.93)式で得られる状態に変化する

測定というのは対象となる量子系(ここでは量子ビット)に対して測定器と結合させることによって測定値が得られるのですが、この際には何らかの物理現象を介しています。
超伝導量子ビットでは初期では電荷計や磁気センサーを、また最近ではマイクロ波の透過・吸収・反射などを利用して読み出しています。
これらの測定器との結合時に起こる物理現象によってそれまで純粋に他と相互作用することなくユニタリー発展していた量子系はその測定後に何らかの状態に遷移してしまうということです。
この遷移するプロセスも厳密には何らかの時間発展なのだと思いますが、ここでは不連続な変化が起きることとします。

一般の量子測定

今回のテーマは「射影測定」なのですが、それ以外の測定もありますので「量子測定」の全体像を整理しますと下記のベン図になります。

図1. 量子測定の全体像

この図の通り、まず一般的な測定があって、その特別な場合として射影測定があります。
まずは一般的な測定から解説します。

測定とは

具体的な方法としては先の電荷計や磁気センサーなどで、量子系から物理量を得て、それに対応する情報を得るということです。

測定演算子

ある測定値 m$m$が得られる測定演算子は M^m${\hat{M}}_m$ で、その集合を{M^m}$\{{\hat{M}}_m\}$と表します。
この時に具体的にどのような測定器でどのような測定をするかは考慮していません。
先の(2.94)式

∑mM†mMm=I(2.94)$\begin{array}{}\text{(}2.94\text{)}& {\displaystyle \sum _m{M}_m^{†}{M}_m=I}\end{array}$

を満たしていればよいです。

測定前の状態を|ψ⟩$|\psi ⟩$とし、測定後の状態を|ψ′⟩$|{\psi }^{\prime }⟩$とした場合に、測定後の状態は測定前の状態と測定演算子を用いて下記で表します。

|ψ′⟩=Mm|ψ⟩⟨ψ|M†mMm|ψ⟩−−−−−−−−−−−√${\displaystyle |{\psi }^{\prime }⟩=\frac{M_m|\psi ⟩}{\sqrt{⟨\psi |M_m^{†}{M}_m|\psi ⟩}}}$

射影測定

量子計算でよく使うのは射影測定で、主に3種類あります。

• X測定(σx${\sigma }_x$測定)
  • 量子状態が|+⟩$|+⟩$か|−⟩$|-⟩$かを測定
• Y測定(σy${\sigma }_y$測定)
  • 光の偏光が|↗⟩$|\nearrow ⟩$か|↙⟩$|\swarrow ⟩$かを測定(本当は両矢印なのですが、tex記号に無い、、)
    (一般の量子状態で簡潔に表せる状態はなく、光の偏光状態くらいが特別な表現があります。)
• Z測定(σz${\sigma }_z$測定)
  • 量子状態が|0⟩$|0⟩$か|1⟩$|1⟩$かを測定

通常は|0⟩$|0⟩$か|1⟩$|1⟩$かを求めるZ測定がほとんどだと思います。
ただ、論文とかで式を追っていったときに⟨σx⟩,⟨σy⟩,⟨σz⟩$⟨{\sigma }_x⟩,⟨{\sigma }_y⟩,⟨{\sigma }_z⟩$が出てきたらやってることは射影測定だということを思い出していただければと思います。

ここでは簡単のために3種類紹介しましたが、一般に正規直交基底となる軸に射影する演算子を用いていれば、軸はX, Y, Zでなくとも良いですが量子計算ではまず使われないと思います。

それぞれの測定がハードウェア的に何をしているかは物理系によります。
偏光であれば|0⟩$|0⟩$は右回り円偏光、|1⟩$|1⟩$は左回り円偏光に対応し、偏光板などを用いて光が通過する/しないを観測します。
超伝導量子ビット(トランズモン)の場合は|0⟩$|0⟩$は基底状態に、|1⟩$|1⟩$は第一励起状態に対応しており、マイクロ波の反射・透過・吸収を計測して、どちらの量子状態にいるかを観測します。

arxiv.org

arxiv.org

その他色々ありますが、量子回路に落とし込んでいるため、プログラミングする際にどのような物理系かは意識しなくて良くなっています。

計算基底(Computational Basis)

1量子ビットの場合、具体的には|0⟩,|1⟩,|+⟩,|−⟩$|0⟩,|1⟩,|+⟩,|-⟩$など、x,y,z$x,y,z$軸上の基底そのものです。
複数量子ビットの場合、具体的には|00⟩,|01⟩,|10⟩,|11⟩$|00⟩,|01⟩,|10⟩,|11⟩$など、1量子ビットの計算基底の組み合わせです。

射影測定とは

基底の軸上で測定で、射影演算子P^m${\hat{P}}_m$を用いた測定です。
具体的にP^m${\hat{P}}_m$が満たすべき条件は下記です。

1. ∑mP^m=I^$\sum _m{\hat{P}}_m=\hat{I}$
2. P^mP^m′=δmm′P^m${\hat{P}}_m{\hat{P}}_{m^{\prime }}={\delta }_{m{m}^{\prime }}{\hat{P}}_m$ (δmm′${\delta }_{m{m}^{\prime }}$はKroneckerのdelta)

これらを満たせれば何でも良いのですが、通常は測定値m$m$に対応した状態|m⟩$|m⟩$を用いて

P^m=|m⟩⟨m|${\displaystyle {\hat{P}}_m=|m⟩⟨m|}$

が用いられます。

Z測定であればm=0,1$m=0,1$(エネルギー固有値Em$E_m$のインデックスm$m$に対応)であり、対応する状態は|0⟩=( 1 0 )T,|1⟩=( 0 1 )T$|0⟩=(10{)}^T,|1⟩=(01{)}^T$で、射影演算子は|0⟩⟨0|,|1⟩⟨1|$|0⟩⟨0|,|1⟩⟨1|$です。
ちゃんと1, 2を満たしていますので行列計算してみてください。

量子回路での対応

実際にそれぞれの基底で測定するにはどのようにプログラミングしてどのような量子回路を構成するのかを確認します。
以下のコードでcircuit.barrier()の上までが状態生成で、circuit.barrier()の下で測定方法を指定します。

Z測定

|0⟩,|1⟩,|+⟩,|−⟩$|0⟩,|1⟩,|+⟩,|-⟩$の4状態を実機を使ってZ基底で測定してみます。

| 0 >を測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.barrier()

circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図2. |0⟩$|0⟩$をZ基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

| 1 >を測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.x(qr)
circuit.barrier()

circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図3. |1⟩$|1⟩$をZ基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

Z基底による測定によって|0⟩$|0⟩$と|1⟩$|1⟩$を見分けられています。

| + >を測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.barrier()

circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図4. |+⟩$|+⟩$をZ基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

図4右のグラフで|0⟩$|0⟩$と表示されている状態が|+⟩$|+⟩$に相当し、|1⟩$|1⟩$が|−⟩$|-⟩$に相当します。

| - >を測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.z(qr)
circuit.barrier()

circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図5. |−⟩$|-⟩$をZ基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

|+⟩,|−⟩$|+⟩,|-⟩$状態はBloch球上では明らか別の状態なのですが、Z基底による測定上では見分けることができません。
下記のX測定で改めて確認しますが、X測定では見分けることが可能です。

X測定

|0⟩,|1⟩,|+⟩,|−⟩$|0⟩,|1⟩,|+⟩,|-⟩$の4状態を実機を使ってX基底で測定してみます。

| 0 >を測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.barrier()

circuit.h(qr)
circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図6. |0⟩$|0⟩$をX基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

| 1 >を測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.x(qr)
circuit.barrier()

circuit.h(qr)
circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図7. |1⟩$|1⟩$をX基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

X基底の場合は|0⟩$|0⟩$状態と|1⟩$|1⟩$状態の見分けがつかなくなりました。

| + >を測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.barrier()

circuit.h(qr)
circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図8. |+⟩$|+⟩$をX基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

| - >を測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.z(qr)
circuit.barrier()

circuit.h(qr)
circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図9. |−⟩$|-⟩$をX基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

X測定の場合は|+⟩$|+⟩$と|−⟩$|-⟩$を見分けることができています。

測定したい状態に合わせて基底を変えましょう！

また、ここでは測定したい状態の方を回転させましたが、偏光を量子ビットとした場合は量子状態はそのままに、偏光板の向きや位相など測定の方を変えることができます。
超伝導量子ビットの場合は状態の方を変えるしかない、、、と思います。
間違ってたらすみません、どなたか教えてください。

Y測定

|+⟩+|−⟩,|+⟩−|−⟩$|+⟩+|-⟩,|+⟩-|-⟩$の2状態を実機を使ってY基底で測定してみます。

Y軸上を向いた状態を見分けられ、X測定・Z測定では見分けられず、Y測定で見分けられます。 まずはX, Z測定してみます。

| + > + | - >をX基底で測定

import numpy as np

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.rz(np.pi / 2, qr)

circuit.barrier()

circuit.h(qr)
circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図10. |+⟩+|−⟩$|+⟩+|-⟩$をX基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

| + > + | - >をZ基底で測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.rz(np.pi / 2, qr)
circuit.barrier()

circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図11. |+⟩+|−⟩$|+⟩+|-⟩$をZ基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

| + > + | - >をY基底で測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.rz(np.pi / 2, qr)
circuit.barrier()

circuit.rz(-np.pi / 2, qr)
circuit.h(qr)
circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図12. |+⟩+|−⟩$|+⟩+|-⟩$をY基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

| + > - | - >をY基底で測定

qr = QuantumRegister(1)
cr = ClassicalRegister(1)
circuit = QuantumCircuit(qr, cr)

circuit.h(qr)
circuit.rz(-np.pi / 2, qr)
circuit.barrier()

circuit.rz(-np.pi / 2, qr)
circuit.h(qr)
circuit.measure(qr, cr)

backend = provider.get_backend('ibmq_armonk')
job = execute(circuit, backend=backend, shots=8192)
job_monitor(job, interval = 2)

result = job.result()
count = result.get_counts()

図13. |+⟩−|−⟩$|+⟩-|-⟩$をY基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

量子測定をもっと知りたい方へ

量子計算する上では計測・観測というのは状態を読み出して終わりですが、量子力学としては重要な理論の一つです。
理論をちゃんと勉強しようとするのは大学・大学院での研究になってしまいますが、一般向けの読み物として下記のようなものがあります。

• 量子と情報
• ハイゼンベルクの顕微鏡
• 量子力学の奇妙なところが思ったほど奇妙でないわけ

　 　

まとめ

量子測定とは何か、数学的な定義、物理系との対応、そして実際に実機で測定する際の量子回路とコーディングをまとめました。

• 量子計算で主に使うのはX測定・Y測定・Z測定
• Z測定はそのままで、X測定はHゲートの後にZ測定、Y測定はRz(−π/2)$R_z(-\pi /2)$ゲート・Hゲートの後にZ測定

今回扱ったのは基本的な状態の観測だけなので、任意の状態をそれぞの基底で観測した時にどういう確率分布になるか試してみてください。

量子計測理論についてもNielsen & Chuangにあります。

PennyLaneを用いた変分分類器(VariationalClassifier)

前回は実用的でなかったので、今回は分類器を作る例をPennyLaneのチュートリアル2つを通して解説します。

今回対象のチュートリアル

pennylane.ai

こちらのチュートリアルではirisデータセットを使って分類アルゴリズムを実装しています。
なお、多クラス分類のチュートリアルも

pennylane.ai

にありますので、機会があればやってみてください。

どのような問題をときたいときに、どのような量子回路を定義すればいいのかなど細かい解説は今後に譲りますが、 全体の流れとしては量子回路の部分は通常の機械学習のネットワークのようなものだと感じてもらえればいいなと思います。

対象とする読者

• 量子計算の応用に興味のある方
• 量子機械学習に興味がある機械学習エンジニア・データサイエンティスト

前提とする知識

• 基本的な機械学習アルゴリズムがわかること
• 簡単な量子ゲート(X, Y, Z, Rx, Ry, Rz, CNOT)がわかること

バージョン情報

• Python 3.9.5
• PennyLane 0.15.1

目次

2クラス変分分類器

変分分類器とは

そもそも「変分分類器(Variational Classifier)」とはなにか？というのが気になると思います。
変分分類器とは、変分アルゴリズムによって最適化した量子回路によって作成された分類器です。

変分アルゴリズムと量子機械学習、量子回路学習

「変分」というのは「変分法」のことです。
機械学習の文脈で「変分法」というと「変分推論」「変分ベイズ学習」などが当てはまると思います。
適用例としては変分オートエンコーダ(Variational Auto Encoder; VAE)がそうです。
求めたい潜在変数の確率分布を何らか(VAEの場合は正規分布)を仮定して、学習で得られる分布との誤差(変分)をなるべく減らして近づける、というものだと思います。

量子計算での変分も似たようなものです。
真に知りたい量子状態|ψ⟩$|\psi ⟩$に対してある時得られた量子状態|ψ′⟩$|{\psi }^{\prime }⟩$との誤差(変分)を誤差関数で評価し、その誤差関数が最小なものが真に知りたかった状態に近いだろう、という考えのもとその量子状態を作り出す量子回路をデータから求めます。
(厳密には物理学の「変分法」は解析力学由来の言葉なので、機械学習文脈での「変分法」とは若干ニュアンスが違いますが、結局の使われどころは同じです。)

回路のパラメータを機械学習で求めることになるので「量子回路学習」とも呼ばれます。
またここで用いられる量子回路は「変分量子回路」と呼ばれます。 この辺りの詳細はまた別に解説することにして、変分アルゴリズムという難しい言い方をしましたが、機械学習と同じようにデータから量子回路のパラメータを最適化してその量子回路を使って分類問題を解きましょうということです。

下記で物理学会誌の記事で京都大学の藤井先生と大阪大学の御手洗先生が解説されています。 https://www.jps.or.jp/books/gakkaishi/2019/09/74-09seriesAIphys1.pdf

その中で示されている量子回路学習の手順の図が機械学習の手順と似ていてわかりやすいと思いました。

図1. 量子回路学習の概念図(こちらより転載)

V(x)$V(\mathbf{x})$で特徴量を量子ビット上にエンコーディングし、U(θ)$U(\theta )$が計算を行う操作で機械学習でいうところのニューラルネットワークなどの機械学習アルゴリズムに相当します。
これらの「量子コンピュータ」と書かれた枠の中だけが機械学習の手順と異なっており、そこを差し替えれば機械学習とやることは同じなのだと思います。

また、この辺りの説明は嶋田さんの本でも解説されています。

量子コンピューティング 基本アルゴリズムから量子機械学習まで [ 情報処理学会出版委員会 ]

2クラス変分分類器の量子回路

分類器の箇所のコードはこちらです。

import pennylane as qml
from pennylane import numpy as np
from pennylane.optimize import NesterovMomentumOptimizer


dev = qml.device("default.qubit", wires=2)


def layer(W):
    qml.Rot(W[0, 0], W[0, 1], W[0, 2], wires=0)
    qml.Rot(W[1, 0], W[1, 1], W[1, 2], wires=1)
    qml.CNOT(wires=[0, 1])


@qml.qnode(dev)
def circuit(weights, angles):
    statepreparation(angles)

    for W in weights:
        layer(W)

    return qml.expval(qml.PauliZ(0))

それを図にするとこちらです。
今回のレイヤー数は6になっています。
(statepreparationはデータセットの特徴量を量子ビットにエンコーディングする箇所(V(x)$V(\mathbf{x})$のこと)なので一旦説明を省きます。)

図2. 2クラス分類の変分量子回路

最後の第一量子ビットを測定することで、2クラスのどちらに分類するかの確率が求まります。
ソースコードではqml.expval(qml.PauliZ(0))となっているので求めているのは⟨σ^z⟩$⟨{\hat{\sigma }}_z⟩$です。
その結果が1,−1$1,-1$どちらであるかで2クラス分類しています。
PennyLaneではどの基底で測定するかを選ぶことができ、今回はZ$Z$基底(|0⟩,|1⟩$|0⟩,|1⟩$の2状態)で測定することが2クラス分類に対応しています。

pennylane.readthedocs.io

ただ、すみません、どうしてこの量子回路だと分類問題が解けるのかというのはまだ理解が足りていません、、
詳細についてはまたの機会に解説解説させていただければと思います。

また、現状の量子回路学習の課題として、解きたい問題に対してどのような回路を構成すべきかという定石(機械学習でいうところの画像だったらCNN, 時系列だったらRNNなど)が無いようです。
これからいろんな問題について適用されて知見が溜まることを期待します。

なお、測定についてはNielsen-Chuangに詳しく書かれています。

【新品】量子コンピュータと量子通信　2　量子コンピュータとアルゴリズム　Michael　A．Nielsen/共著　Isaac　L．Chuang/共著　木村達也/訳

NIIの根本先生の本も物理現象に即して解説されていてわかりやすいです。

www.saiensu.co.jp

データ読み込み〜モデル学習までの基本的な流れ

振幅エンコーディング

まずデータ(問題設定)を量子回路に置き換える方法は主に次の2通りがあります。

• 基底エンコーディング：量子ビットの |0⟩,|1⟩$|0⟩,|1⟩$をフラグのように対応させます。
  以前解説させていただいた、Groverのアルゴリズムの入力がこちらに当たります。
• 振幅エンコーディング：量子ビットの確率振幅(の位相)をデータに対応させます。
  今回用いるのがこちらの手法です。

チュートリアルでは下記になっています。

def get_angles(x):

    beta0 = 2 * np.arcsin(np.sqrt(x[1] ** 2) / np.sqrt(x[0] ** 2 + x[1] ** 2 + 1e-12))
    beta1 = 2 * np.arcsin(np.sqrt(x[3] ** 2) / np.sqrt(x[2] ** 2 + x[3] ** 2 + 1e-12))
    beta2 = 2 * np.arcsin(
        np.sqrt(x[2] ** 2 + x[3] ** 2)
        / np.sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2 + x[3] ** 2)
    )
    return np.array([beta2, -beta1 / 2, beta1 / 2, -beta0 / 2, beta0 / 2])


def statepreparation(a):
    qml.RY(a[0], wires=0)

    qml.CNOT(wires=[0, 1])
    qml.RY(a[1], wires=1)
    qml.CNOT(wires=[0, 1])
    qml.RY(a[2], wires=1)

    qml.PauliX(wires=0)
    qml.CNOT(wires=[0, 1])
    qml.RY(a[3], wires=1)
    qml.CNOT(wires=[0, 1])
    qml.RY(a[4], wires=1)
    qml.PauliX(wires=0)


@qml.qnode(dev)
def circuit(weights, angles):
    statepreparation(angles)
    for W in weights:
        layer(W)
    return qml.expval(qml.PauliZ(0))


def variational_classifier(var, angles):
    weights = var[0]
    bias = var[1]
    return circuit(weights, angles) + bias


def cost(weights, features, labels):
    predictions = [variational_classifier(weights, f) for f in features]
    return square_loss(labels, predictions)


# irisデータセットの読み込み
data = np.loadtxt("variational_classifier/data/iris_classes1and2_scaled.txt")
X = data[:, 0:2]

# パディング
padding = 0.3 * np.ones((len(X), 1))
X_pad = np.c_[np.c_[X, padding], np.zeros((len(X), 1))]

# 正規化(標準化)
normalization = np.sqrt(np.sum(X_pad ** 2, -1))
X_norm = (X_pad.T / normalization).T

# 角度(位相)に変換
features = np.array([get_angles(x) for x in X_norm])

今回のirisもそうですが一般に特徴量は実数です。
対して回転ゲートU^(θ)$\hat{U}(\theta )$は2π$2\pi$の周期性があります。
特徴量の値が大きいと何周回ったのかわからなくなります。
なのでできるだけ1周程度に収めるために正規化(ここでは標準化)をしています。

そしてそのあと、学習処理

var = opt.step(lambda v: cost(v, feats_train_batch, Y_train_batch), var)

にて、実際に振幅エンコーディングします。
(cost関数にて量子回路を定義しています。必要な関数は並べてあるので、お手数ですが少し式を追ってください。)

誤差関数

測定結果⟨σ^z⟩$⟨{\hat{\sigma }}_z⟩$の値に対して、教師データとの精度を誤差関数で評価します。

def square_loss(labels, predictions):
    loss = 0
    for l, p in zip(labels, predictions):
        loss = loss + (l - p) ** 2

    loss = loss / len(labels)
    return loss


def cost(weights, features, labels):
    predictions = [variational_classifier(weights, f) for f in features]
    return square_loss(labels, predictions)

予測結果の正解ラベルへの当てはまり具合を平均自乗誤差で評価します。
この辺りはもう機械学習と同じ流れです。

最適化アルゴリズム

誤差関数での評価結果に基づいて、パラメータ更新をします。

opt = NesterovMomentumOptimizer(0.01)

var = opt.step(lambda v: cost(v, feats_train_batch, Y_train_batch), var)

varが変分量子回路U^(w1,w2,w3),U^(w4,w5,w6)$\hat{U}(w_1,w_2,w_3),\hat{U}(w_4,w_5,w_6)$ x 6層分ののパラメータw1,w2,w3,…w36$w_1,w_2,w_3,\dots w_{36}$です。

分類結果

学習済みの量子回路を使って特徴量空間をクラスで色分けしたのが下図です。

図3. 変分量子回路を使ってirisデータセットを2クラス分類した結果(こちらより転載)

まとめ

• 量子回路(変分量子回路)を用いた2クラス分類ができることを確かめた
• 量子回路学習とは量子回路のパラメータを機械学習的に最適化することで、この時の量子回路を変分量子回路という
• どのような問題設定の時にどのような変分量子回路にすればいいかという定石は定まってなく、個別に計算式・量子回路を考えて実装する必要がある