科学しよう

量子計算のプログラミングの解説をメインに、データサイエンス・機械学習について勉強したことをご紹介します

MENU

スワップテスト

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

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

前提とする知識

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

バージョン情報

目次

SWAPテスト回路

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

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

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

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

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

\displaystyle{
p_0 = \frac{1 + \mathrm{Re} \langle \psi |_B \hat{U} | \psi \rangle_B}{2}
}
\displaystyle{
p_1 = \frac{1 - \mathrm{Re} \langle \psi |_B \hat{U} | \psi \rangle_B}{2}
}

SWAPテスト回路の理論

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

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

\displaystyle{
| a \rangle = \sum_{i=0}^{N-1} {a_i | i \rangle}
}
\displaystyle{
| b \rangle = \sum_{i=0}^{N-1} {b_i | i \rangle}
}

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

\displaystyle{
| \langle a | b \rangle|^2 = | \sum_{i=0}^{N-1} {a_i b_i}|^2
}

を評価できます。

図2. SWAPテスト回路

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

\displaystyle{
p_0 = \frac{1 + | \langle a | b \rangle |^2}{2}
}

SWAPテストの仕組み導出

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

\displaystyle{
| 0 \rangle | a \rangle | b \rangle \overset{H_1}\longrightarrow \frac{1}{\sqrt{2}} \left( | 0 \rangle | a \rangle | b \rangle + | 1 \rangle | a \rangle | b \rangle \right)
}

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

\displaystyle{
\frac{1}{\sqrt{2}} \left( | 0 \rangle | a \rangle | b \rangle + | 1 \rangle | b \rangle | a \rangle \right)
}

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


\begin{aligned}
| \psi \rangle &= \frac{1}{2} \left\{ \left( | 0 \rangle + | 1 \rangle \right) | a \rangle | b \rangle + \left( | 0 \rangle - | 1 \rangle \right) | b \rangle | a \rangle \right\} \\

&= \frac{1}{2} \left\{ | 0 \rangle \left( | a \rangle | b \rangle + | b \rangle | a \rangle \right) + | 1 \rangle \left( | a \rangle | b \rangle - | b \rangle | a \rangle \right) \right\}
\end{aligned}

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


\begin{aligned}
p_0 &= | \langle 0 | \psi \rangle |^2 \\

&= ( \langle 0 | \psi \rangle)^{\dagger}  \langle 0 | \psi \rangle \\

&= \langle \psi | 0 \rangle \langle 0 | \psi \rangle \\

&= \frac{1}{4} (\langle a | \langle b | + \langle b | \langle a | ) (| a \rangle | b \rangle + | b \rangle | a \rangle) \\

&= \frac{1}{4} (\langle a | \langle b | a \rangle | b \rangle + \langle a | \langle b | b \rangle | a \rangle + \langle b | \langle a | a \rangle | b \rangle + \langle b | \langle a | b \rangle | a \rangle)
\end{aligned}

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


\begin{aligned}
p_0 &= \frac{1}{4} ( \langle a | b \rangle \langle b | a \rangle + \langle a | a \rangle \langle b | b \rangle + \langle b | b \rangle \langle a | a \rangle + \langle b | a \rangle \langle a | b \rangle ) \\

&= \frac{1}{4} (| \langle a | b \rangle|^2 + 1 + 1 + |\langle a | b \rangle|^2 ) \\

&= \frac{1}{4} \left( 2 + 2 |\langle a | b \rangle|^2 \right) \\

&= \frac{1 + |\langle a | b \rangle|^2}{2}
\end{aligned}

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

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

SWAPゲートを示すユニタリー \hat{U}_{\mathrm{SWAP}}は以下です。

\displaystyle{
\hat{U}_{\mathrm{SWAP}} = \begin{pmatrix}
1 & 0 &0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}
}

この行列の固有ベクトル


\begin{aligned}
| x_+ \rangle &= \begin{pmatrix}
0 \\
1 \\
1 \\
0
\end{pmatrix} \\

&= | 1 \rangle | 0 \rangle + | 0 \rangle | 1 \rangle

\end{aligned}

\begin{aligned}
| x_- \rangle &= \begin{pmatrix}
0 \\
1 \\
-1 \\
0
\end{pmatrix} \\

&= | 1 \rangle | 0 \rangle - | 0 \rangle | 1 \rangle

\end{aligned}

のため、一般の入力状態


\begin{aligned}
| \psi \rangle &= | a \rangle | b \rangle \\

&= \left( \alpha | 0 \rangle + \beta | 1 \rangle \right) \left( u | 0 \rangle + v | 1 \rangle \right) \\

&= \alpha u | 0 \rangle | 0 \rangle + \alpha v | 0 \rangle | 1 \rangle + \beta u | 1 \rangle | 0 \rangle + \beta v | 1 \rangle | 1 \rangle \\

&= \begin{pmatrix}
\alpha u \\
\alpha v \\
\beta u \\
\beta v
\end{pmatrix}

\end{aligned}

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

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


\begin{aligned}
p_0 &= \frac{1 + \mathrm{Re} \langle \psi | \hat{U}_{\mathrm{SWAP}} | \psi \rangle}{2} \\

&= \frac{1 + \alpha^2 u^2 + 2 \alpha \beta u v + \beta^2 v^2}{2} \\

&= \frac{1 + (\alpha u + \beta v)^2}{2} \\

&= \frac{1 + \left| \begin{pmatrix} \alpha & \beta \end{pmatrix} \begin{pmatrix} u \\ v \end{pmatrix} \right|^2}{2} \\

&= \frac{1 + \left| \langle a | b \rangle \right|^2}{2}

\end{aligned}

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

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

SWAPテスト回路の実装

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

単純なベクトルの計算

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

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

わかりやすい例だと | a \rangle = | 0 \rangle, | b \rangle = | 0 \rangleなら \langle a | b \rangle = 1 | a \rangle = | 0 \rangle, | b \rangle = | 1 \rangleなら \langle a | b \rangle = 0です。

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


\begin{aligned}
\langle a | b \rangle &= \langle 0 | \left( \cos \frac{\theta}{2} | 0 \rangle + \sin \frac{\theta}{2} | 1 \rangle \right) \\

&= \cos \frac{\theta}{2}

\end{aligned}

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

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

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

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

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

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

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

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

\displaystyle{
p_0 = \frac{1 + \left| \langle a | b \rangle \right|^2}{2}
}

から \left| \langle a | b \rangle \right|^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

まとめ

参考


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

SWAPゲート・制御SWAPゲート

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

対象とする読者

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

前提とする知識

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

バージョン情報

目次

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

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

\displaystyle{
( \alpha | 0 \rangle + \beta | 1 \rangle)_1 \otimes (u | 0 \rangle + v | 1 \rangle)_2 \longleftrightarrow (u | 0 \rangle + v | 1 \rangle)_1 \otimes ( \alpha | 0 \rangle + \beta | 1 \rangle)_2
}

ができます。

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

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

\displaystyle{
\left( \frac{1}{2} | 0 \rangle + \frac{\sqrt{3}}{2} | 1 \rangle \right)_1 \otimes \left( \frac{\sqrt{3}}{2} | 0 \rangle + \frac{1}{2} | 1 \rangle \right)_2 = \frac{\sqrt{3}}{4} | 00 \rangle + \frac{1}{4} | 01 \rangle + \frac{3}{4} | 10 \rangle + \frac{\sqrt{3}}{4} | 11 \rangle
}

ビット列の読み順は式とグラフとで逆順なのでグラフ上での | 01 \rangleが式上の | 10 \rangleに対応していて確率振幅0.75の自乗で約0.56、グラフ上の | 10 \rangleが式上の | 01 \rangleに対応していて確率振幅0.25の自乗で約0.05、 | 00 \rangle | 11 \rangleはそのままで確率振幅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ゲート実装

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

\displaystyle{
| 000 \rangle_1 \otimes | 011 \rangle_2 \longleftrightarrow | 011 \rangle_1 \otimes | 000 \rangle_2
}

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

図4. 3量子ビットのSWAP

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

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

第一量子ビットの状態が | 011 \rangleで第二量子ビットの状態が | 000 \rangleになっていることを示しています。 よってちゃんと思った通りに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 \rangle (第一量子ビットが制御量子ビット)に対して、測定結果は | 010 \rangleのため、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 \rangle (第一量子ビットが制御量子ビット)に対して、測定結果は | 101 \rangleのため、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_quito1量子ビットを想定して、ancilaとして使われる想定です。 ただ、実行したjobを確認すると思ったようには動いてなさそうなので、transplieの仕組みを理解する必要があります。

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

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

まとめ

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

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

アダマールテスト

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

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

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

対象とする読者

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

前提とする知識

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

バージョン情報

目次

アダマールテストとは

量子回路

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

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

用途

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

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

代数的に解いてみる

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

ユニタリー行列とは

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

ユニタリ行列 - Wikipedia

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

\displaystyle{
\hat{U} \hat{U}^{\dagger} = \hat{U}^{\dagger} \hat{U} = \hat{I}
}

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

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

\displaystyle{
i \hbar \frac{d}{dt} | \psi (t) \rangle = \hat{H} | \psi (t) \rangle
}

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

\displaystyle{
| \psi (t) \rangle = \exp \left \{ \frac{1}{i \hbar} \hat{H}(t - t_0) \right \} | \psi (t_0) \rangle
}

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

\displaystyle{
\exp \left \{ \frac{1}{i \hbar} \hat{H} (t - t_0) \right \} \longrightarrow \exp \left \{ \frac{-1}{i \hbar} \hat{H} (t - t_0) \right \}
}

となるだけで、これは

\displaystyle{
\exp i \Delta \theta \longrightarrow \exp ( - i \Delta \theta)
}

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

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

\displaystyle{
\exp (i \hat{H}) = \sum_{k=0}^{\infty}{\frac{1}{k!}(i \hat{H})^k} \equiv \hat{U}
}

とすると


\begin{aligned}
\hat{U}^{\dagger} &= \sum_{k=0}^{\infty}{\frac{1}{k!}(-i \hat{H}^{\dagger})^k} \\
&= \sum_{k=0}^{\infty}{\frac{1}{k!}(-i \hat{H})^k} \\
&= \exp (-i \hat{H})
\end{aligned}

ユニタリー行列の固有値

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

risalc.info

eman-physics.net

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

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

\displaystyle{
\hat{U} | \psi \rangle = e^{i \theta} | \psi \rangle
}

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

部分トレース

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

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

sakumatcho.hatenablog.com

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

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


\begin{aligned}
\frac{1}{\sqrt{2}} \left( | 0 \rangle + | 1 \rangle \right) \otimes | \psi \rangle &\overset{\text{Ctrl-}U}\longrightarrow \frac{1}{\sqrt{2}} \left( | 0 \rangle \otimes | \psi \rangle + | 1 \rangle \otimes  \hat{U} | \psi \rangle \right) \\

&= \frac{1}{\sqrt{2}} \left( | 0 \rangle \otimes | \psi \rangle + e^{i \lambda }| 1 \rangle \otimes  | \psi \rangle \right) \\

&= \frac{1}{\sqrt{2}} \left( | 0 \rangle + e^{i \lambda } | 1 \rangle \right) \otimes | \psi \rangle \\
\end{aligned}

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

\displaystyle{
\frac{1}{\sqrt{2}} \left\{ \frac{1}{\sqrt{2}} \left( | 0 \rangle + | 1 \rangle \right) + e^{i \lambda } \frac{1}{\sqrt{2}} \left( | 0 \rangle - | 1 \rangle \right) \right\}
=
\frac{1 + e^{i \lambda}}{2} | 0 \rangle + \frac{1 - e^{i \lambda}}{2} | 1 \rangle
}

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


\begin{aligned}
&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{aligned}

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

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

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

\displaystyle{
\frac{1}{\sqrt{2}} \left( | 0 \rangle + | 1 \rangle \right) \otimes | \psi \rangle \overset{\text{Ctrl-}U}\longrightarrow \frac{1}{\sqrt{2}} \left( | 0 \rangle \otimes | \psi \rangle + | 1 \rangle \otimes  \hat{U} | \psi \rangle \right)
}

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

\displaystyle{
\frac{1}{\sqrt{2}} \left\{ \left( | 0 \rangle + | 1 \rangle \right) \otimes | \psi \rangle + \left( | 0 \rangle - | 1 \rangle \right) \otimes  \hat{U} | \psi \rangle \right\}
= \frac{1}{2} \left\{ | 0 \rangle \otimes \left( | \psi \rangle + \hat{U} | \psi \rangle \right) + | 1 \rangle \otimes \left( | \psi \rangle - \hat{U} | \psi \rangle \right) \right\}
}

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

\displaystyle{
\hat{\rho}^{AB} = \frac{1}{4} \left\{ | 0 \rangle_{A} \otimes \left( \hat{I} + \hat{U} \right) | \psi \rangle_{B} + | 1 \rangle_{A} \otimes \left( \hat{I} - \hat{U} \right) | \psi \rangle_{B} \right\} \left\{ \langle 0 |_{A} \otimes \langle \psi |_{B} \left( \hat{I} + \hat{U} \right)^{\dagger} + \langle 1 |_{A} \otimes \langle \psi |_{B} \left( \hat{I} - \hat{U} \right)^{\dagger} \right\} \\
}

から


\begin{aligned}
\hat{\rho}^{A}

= &\frac{1}{4} 
\mathrm{tr}_B \left\{ | 0 \rangle_A \otimes \left( \hat{I} + \hat{U} \right) | \psi \rangle_B \langle 0 |_A \otimes \langle \psi | _B \left( \hat{I} + \hat{U} \right)^{\dagger} \right\} \\

&+ \frac{1}{4} \mathrm{tr}_B \left\{ | 0 \rangle_A \otimes \left( \hat{I} + \hat{U} \right) | \psi \rangle_B \langle 1 |_A \otimes \langle \psi | _B \left( \hat{I} - \hat{U} \right)^{\dagger} \right\} \\

&+ \frac{1}{4} \mathrm{tr}_B \left\{ | 1 \rangle_A \otimes \left( \hat{I} - \hat{U} \right) | \psi \rangle_B \langle 0 |_A \otimes \langle \psi | _B \left( \hat{I} + \hat{U} \right)^{\dagger} \right\} \\

&+ \frac{1}{4} \mathrm{tr}_B \left\{ | 1 \rangle_A \otimes \left( \hat{I} - \hat{U} \right) | \psi \rangle_B \langle 1 |_A \otimes \langle \psi | _B \left( \hat{I} - \hat{U} \right)^{\dagger} \right\} \\


= &\frac{1}{4} 
\mathrm{tr}_B \left\{ | 0 \rangle_A \langle 0 |_A \otimes \left( \hat{I} + \hat{U} \right) | \psi \rangle_B  \langle \psi | _B \left( \hat{I} + \hat{U} \right)^{\dagger} \right\} \\

&+ \frac{1}{4} \mathrm{tr}_B \left\{ | 0 \rangle_A \langle 1 |_A \otimes \left( \hat{I} + \hat{U} \right) | \psi \rangle_B  \langle \psi | _B \left( \hat{I} - \hat{U} \right)^{\dagger} \right\} \\

&+ \frac{1}{4} \mathrm{tr}_B \left\{ | 1 \rangle_A \langle 0 |_A \otimes \left( \hat{I} - \hat{U} \right) | \psi \rangle_B  \langle \psi | _B \left( \hat{I} + \hat{U} \right)^{\dagger} \right\} \\

&+ \frac{1}{4} \mathrm{tr}_B \left\{ | 1 \rangle_A \langle 1 |_A \otimes \left( \hat{I} - \hat{U} \right) | \psi \rangle_B  \langle \psi | _B \left( \hat{I} - \hat{U} \right)^{\dagger} \right\} \\


= &\frac{1}{4} | 0 \rangle_A \langle 0 |_A \mathrm{tr} \left\{ \left( \hat{I} + \hat{U} \right) | \psi \rangle_B  \langle \psi | _B \left( \hat{I} + \hat{U} \right)^{\dagger} \right\} \\

&+ \frac{1}{4} | 0 \rangle_A \langle 1 |_A \mathrm{tr} \left\{ \left( \hat{I} + \hat{U} \right) | \psi \rangle_B  \langle \psi | _B \left( \hat{I} - \hat{U} \right)^{\dagger} \right\} \\

&+ \frac{1}{4} | 1 \rangle_A \langle 0 |_A \mathrm{tr} \left\{ \left( \hat{I} - \hat{U} \right) | \psi \rangle_B  \langle \psi | _B \left( \hat{I} + \hat{U} \right)^{\dagger} \right\} \\

&+ \frac{1}{4} | 1 \rangle_A \langle 1 |_A \mathrm{tr} \left\{ \left( \hat{I} - \hat{U} \right) | \psi \rangle_B  \langle \psi | _B \left( \hat{I} - \hat{U} \right)^{\dagger} \right\} \\


= &\frac{1}{4} | 0 \rangle_A \langle 0 |_A \langle \psi | _B \left( \hat{I} + \hat{U} \right)  \left( \hat{I} + \hat{U} \right)^{\dagger} | \psi \rangle_B  \\

&+ \frac{1}{4} | 0 \rangle_A \langle 1 |_A \langle \psi | _B \left( \hat{I} - \hat{U} \right)  \left( \hat{I} + \hat{U} \right)^{\dagger} | \psi \rangle_B  \\

&+ \frac{1}{4} | 1 \rangle_A \langle 0 |_A \langle \psi | _B \left( \hat{I} + \hat{U} \right)  \left( \hat{I} - \hat{U} \right)^{\dagger} | \psi \rangle_B  \\

&+ \frac{1}{4} | 1 \rangle_A \langle 1 |_A \langle \psi | _B \left( \hat{I} - \hat{U} \right)  \left( \hat{I} - \hat{U} \right)^{\dagger} | \psi \rangle_B  \\


= &\frac{1}{4} | 0 \rangle_A \langle 0 |_A \left( \langle \psi |_B 2 \hat{I} | \psi \rangle_B + \langle \psi |_B \hat{U} | \psi \rangle_B + \langle \psi |_B \hat{U}^{\dagger} | \psi \rangle_B \right) \\

&+ \mathrm{略} \\

&+ \frac{1}{4} | 1 \rangle_A \langle 1 |_A \left( \langle \psi |_B 2 \hat{I} | \psi \rangle_B - \langle \psi |_B \hat{U} | \psi \rangle_B - \langle \psi |_B \hat{U}^{\dagger} | \psi \rangle_B \right) \\


= &\frac{1}{4} | 0 \rangle_A \langle 0 |_A \left\{ 2 + \langle \psi |_B \hat{U} | \psi \rangle_B + \left( \langle \psi |_B \hat{U} | \psi \rangle_B \right)^{*} \right\} \\

&+ \mathrm{略} \\

&+ \frac{1}{4} | 1 \rangle_A \langle 1 |_A \left\{ 2 - \langle \psi |_B \hat{U} | \psi \rangle_B - \left( \langle \psi |_B \hat{U} | \psi \rangle_B \right)^{*} \right\} \\


= &\frac{1}{4} | 0 \rangle_A \langle 0 |_A \left\{ 2 + \mathrm{Re} \left( \langle \psi |_B \hat{U} | \psi \rangle_B \right) \right\} \\

&+ \mathrm{略} \\

&+ \frac{1}{4} | 1 \rangle_A \langle 1 |_A \left\{ 2 - \mathrm{Re} \left( \langle \psi |_B \hat{U} | \psi \rangle_B \right) \right\}

\end{aligned}

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

\displaystyle{
p_0 = \frac{1 + \mathrm{Re} \langle \psi |_B \hat{U} | \psi \rangle_B}{2}
}
\displaystyle{
p_1 = \frac{1 - \mathrm{Re} \langle \psi |_B \hat{U} | \psi \rangle_B}{2}
}

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

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

入力する量子ビットの状態は固有状態ではないが、(固有値 \lambda_1, \lambda_2で、固有状態が | u_1 \rangle, | u_2 \rangleと予めわかっているとすると)以下のように書き換えられます。

\displaystyle{
| \psi \rangle = c_1 | u_1 \rangle + c_2 | u_2 \rangle
}

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


\begin{aligned}
| \psi \rangle \pm \hat{U} | \psi \rangle &= c_1 | u_1 \rangle + c_2 | u_2 \rangle \pm c_1 \hat{U} |u_1 \rangle \pm c_2 \hat{U} | u_2 \rangle \\

&= c_1 | u_1 \rangle + c_2 | u_2 \rangle \pm c_1 e^{i \lambda_1} | u_1 \rangle \pm c_2 e^{i \lambda_2} | u_2 \rangle \\

&= c_1 \left( 1 \pm e^{i \lambda_1} \right) | u_1 \rangle + c_2 \left( 1 \pm e^{i \lambda_2} \right) | u_2 \rangle

\end{aligned}

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

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

制御Xゲートの場合

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

\displaystyle{
| \psi \rangle + \hat{X} | \psi \rangle = 2 c_1 | + \rangle
}
\displaystyle{
| \psi \rangle - \hat{X} | \psi \rangle = 2 c_2 | - \rangle
}

となります。

制御Zゲートの場合

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

\displaystyle{
| \psi \rangle + \hat{Z} | \psi \rangle = 2 c_1 | 0 \rangle
}
\displaystyle{
| \psi \rangle - \hat{Z} | \psi \rangle = 2 c_2 | 1 \rangle
}

となります。

制御Yゲートの場合

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

\displaystyle{
| \psi \rangle + \hat{Y} | \psi \rangle = 2 c_1 \left( | + \rangle + i | - \rangle \right)
}
\displaystyle{
| \psi \rangle - \hat{Y} | \psi \rangle = 2 c_2 \left( | + \rangle - i | - \rangle \right)
}

となります。

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

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

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

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

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


\begin{aligned}
p_0 &= \frac{1 + 1}{2}  \\

&= 1
\end{aligned}

ii) 第二量子ビットの状態が | \psi \rangle = | - \rangleの時


\begin{aligned}
p_0 &= \frac{1 - 1}{2} \\

&= 0
\end{aligned}

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)の場合アダマール2回やって | 0 \rangleに戻ってきただけ、と思われても嫌なのでii)も計算しておきます。

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

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

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

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

実機の結果です。

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

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

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

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ゲートの場合

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

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

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

i) 第二量子ビットの状態が | \psi \rangle = | - \rangleの時、 \lambda = \lambda_{+}のため


\begin{aligned}
p_0 &= \frac{1 + \sqrt{2} / 2}{2} \\

&= \frac{2 + \sqrt{2}}{4} \\

&\sim 0.85
\end{aligned}

\begin{aligned}
p_1 &= \frac{1 - \sqrt{2} / 2}{2} \\

&= \frac{2 - \sqrt{2}}{4} \\

&\sim 0.15
\end{aligned}

ii) 第二量子ビットの状態が | \psi \rangle = | + \rangleの時、 \lambda = \lambda_{-}のため


\begin{aligned}
p_0 &= \frac{1 + \sqrt{2} / 2}{2} \\

&= \frac{2 + \sqrt{2}}{4} \\

&\sim 0.15

\end{aligned}

\begin{aligned}
p_1 &= \frac{1 - \sqrt{2} / 2}{2} \\

&= \frac{2 - \sqrt{2}}{4} \\

&\sim 0.15

\end{aligned}

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 \rangle + i | 1 \rangle) / \sqrt{2}でユニタリーはXゲートの場合

 \langle \psi | \hat{X} | \psi \rangle複素数になる場合です。


\begin{aligned}
\langle \psi | \hat{X} | \psi \rangle &= \frac{1}{2} \begin{pmatrix} 1&i \end{pmatrix} \begin{pmatrix} 0&1\\ 1&0 \end{pmatrix} \begin{pmatrix} 1\\ i \end{pmatrix} \\

&= \frac{1}{2} \begin{pmatrix} i&1 \end{pmatrix} \begin{pmatrix} 1\\ i \end{pmatrix} \\

&= i
\end{aligned}

従って


\begin{aligned}
p_0 &= \frac{1 + 0}{2} \\

&= 0.5

\end{aligned}

\begin{aligned}
p_1 &= \frac{1 - 0}{2} \\

&= 0.5

\end{aligned}

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

図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())

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

 \langle \psi | \hat{X} | \psi \rangleが半端な実数になる場合です。


\begin{aligned}
\langle \psi | \hat{X} | \psi \rangle &= \frac{1}{4} \begin{pmatrix} \sqrt{3}&1 \end{pmatrix} \begin{pmatrix} 0&1\\ 1&0 \end{pmatrix} \begin{pmatrix} \sqrt{3}\\ 1 \end{pmatrix} \\

&= \frac{1}{4} \begin{pmatrix} 1&\sqrt{3} \end{pmatrix} \begin{pmatrix} \sqrt{3}\\ 1 \end{pmatrix} \\

&= \frac{\sqrt{3}}{2}
\end{aligned}

従って


\begin{aligned}
p_0 &= \frac{1 + \sqrt{3} / 2}{2} \\

&\sim 0.93

\end{aligned}

\begin{aligned}
p_1 &= \frac{1 - \sqrt{3} / 2}{2} \\

&\sim 0.07

\end{aligned}

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

図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())

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

まとめ

参考


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

密度演算子・密度行列

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

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

対象とする読者

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

など

前提とする知識

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

バージョン情報

目次

定義

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

\displaystyle{
\rho = \sum_i{p_i | \psi_i \rangle \langle \psi_i |}
}

ここで、 | \psi_i \rangleは量子系のとりうる状態で、 p_iはそれに対する確率(実数)です。

性質

トレースは1

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


\begin{aligned}
\mathrm{tr} (\rho) &= \mathrm{tr} \left( \sum_i {p_i | \psi_i \rangle \langle \psi_i |} \right) \\
&= \sum_i {p_i \mathrm{tr} \left(| \psi_i \rangle \langle \psi_i | \right)} \\
&= \sum_i {p_i} \\
&= 1
\end{aligned}

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

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

\displaystyle{
| \psi_i \rangle = \sum_j{\alpha_j | \psi_{ij} \rangle} \left( \sum_j{| \alpha_j |^2} = 1 \right)
}

 | \psi_{ij} \rangleは正規直交基底のため、 \mathrm{tr}(| \psi_i \rangle \langle \psi_i |) = \sum_j{| \alpha_j |^2} = 1となることを利用しています。

対角成分はゼロ以上

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


\begin{aligned}
\langle \phi | \rho | \phi \rangle &= \langle \phi | \left( \sum_i{p_i | \psi_i \rangle \langle \psi_i |} \right) | \phi \rangle \\
&= \sum_i{p_i \langle \phi | \psi_i \rangle \langle \psi_i | \phi \rangle} \\
&= \sum_i{p_i | \langle \phi | \psi_i \rangle |^2} \\
&\geq 0 \ (\because p_i \geq 0, | \langle \phi | \psi_i \rangle |^2 \geq 0)
\end{aligned}

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

定理

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

状態 | \psi_i \rangleがユニタリー発展して \hat{U} | \psi_i \rangleとなるので


\begin{aligned}
\rho_{\mathrm{init}} = \sum_i{p_i | \psi_i \rangle \langle \psi_i |} &\xrightarrow{\hat{U}} \sum_i{p_i \hat{U} | \psi_i \rangle \langle \psi_i | \hat{U}^{\dagger}}\\
&= \hat{U} \left( \sum_i{p_i | \psi_i \rangle \langle \psi_i |} \right)\hat{U}^{\dagger} \\
&= \hat{U} \rho_{\mathrm{init}} \hat{U}^{\dagger} \equiv \rho_{\mathrm{final}}
\end{aligned}

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

測定後の密度演算子

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

\displaystyle{
p(m) = \langle \psi | \hat{M}^{\dagger}_m \hat{M}_m | \psi \rangle 
}

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


\begin{aligned}
p (m | i) &= \langle \psi_i | \hat{M}^{\dagger}_m \hat{M}_m | \psi_i \rangle \\
&= \mathrm{tr} \left( \hat{M}^{\dagger}_m \hat{M}_m | \psi_i \rangle \langle \psi_i | \right)
\end{aligned}

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

\displaystyle{
\mathrm{tr} \left( \hat{A} | \psi \rangle \langle \psi | \right) = \langle \psi | \hat{A} | \psi \rangle
}

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


\begin{aligned}
p(m) &= \sum_i{p(m | i) \cdot p_i} \\
&= \sum_i{p_i \cdot \mathrm{tr}\left( \hat{M}^{\dagger}_m \hat{M}_m | \psi_i \rangle \langle \psi_j | \right)} \\
&= \mathrm{tr} \left( \hat{M}^{\dagger}_m \hat{M}_m \sum_i{p_i | \psi_i \rangle \langle \psi_i |} \right)
\end{aligned}

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

\displaystyle{
| \psi_i^m \rangle = \frac{\hat{M}_m | \psi_i \rangle}{\sqrt{\langle \psi_i | \hat{M}_m^{\dagger} \hat{M}_m | \psi_i \rangle}}
}

となることを用いて、測定後の密度演算子 \rho_mは以下になる


\begin{aligned}
\rho_m &= \sum_i{p(i | m) | \psi_i^m \rangle \langle \psi_i^m |} \\
&= \sum_i{p(i | m) \frac{\hat{M}_m | \psi_i \rangle \langle \psi_i | \hat{M}^{\dagger}_m}{\langle \psi_i | \hat{M}^{\dagger}_m \hat{M}_m | \psi_i \rangle}}
\end{aligned}

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

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

より


\begin{aligned}
\rho_m &= \sum_i{\frac{p(m | i) \cdot p_i}{p(m)} \cdot \frac{\hat{M}_m | \psi_i \rangle \langle \psi_i | \hat{M}^{\dagger}_m}{\langle \psi_i | \hat{M}^{\dagger}_m \hat{M}_m | \psi_i \rangle}} \\
&= \sum_i{\frac{\langle \psi_i | \hat{M}^{\dagger}_m \hat{M}_m | \psi_i \rangle}{\mathrm{tr} \left( \hat{M}^{\dagger}_m \hat{M}_m \rho \right)} \cdot p_i \cdot \frac{\hat{M}_m | \psi_i \rangle \langle \psi_i | \hat{M}^{\dagger}_m}{\langle \psi_i | \hat{M}^{\dagger}_m \hat{M}_m | \psi_i \rangle}} \\
&= \sum_i{p_i \cdot \frac{\hat{M}_m | \psi_i \rangle \langle \psi_i | \hat{M}^{\dagger}_m}{\mathrm{tr} \left( \hat{M}^{\dagger}_m \hat{M}_m \rho \right)}} \\
&= \frac{\hat{M}_m \sum_i{p_i | \psi_i \rangle \langle \psi_i |} \hat{M}^{\dagger}_m}{\mathrm{tr} \left( \hat{M}^{\dagger}_m \hat{M}_m \rho \right)} \\
&= \frac{\hat{M}_m \rho \hat{M}^{\dagger}_m}{\mathrm{tr} \left( \hat{M}^{\dagger}_m \hat{M}_m \right)}
\end{aligned}

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

ある状態 \rho_i = \sum_j{p_j | \psi_j \rangle \langle \psi_j |}が確率 p_iで存在するとき、全体の密度演算子


\begin{aligned}
\rho &= \sum_{i, j}{p_i p_j | \psi_j \rangle \langle \psi_j |} \\
&= \sum_i{p_i \rho_i}
\end{aligned}

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

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

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

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

密度演算子の具体例

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

純粋状態

測定対象となる状態は

\displaystyle{
| \psi \rangle = \frac{1}{\sqrt{2}} \left( | 0 \rangle + | 1 \rangle \right)
}

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

\displaystyle{
| 0 \rangle \equiv 
\begin{pmatrix} 1 \\ 0 \end{pmatrix}
}
\displaystyle{
| 1 \rangle \equiv 
\begin{pmatrix} 0 \\ 1 \end{pmatrix}
}

とすると、密度演算子


\begin{aligned}
\rho &= | \psi \rangle \langle \psi | \\
&= \frac{1}{2} \left( | 0 \rangle \langle 0| + | 0 \rangle \langle 1 | + | 1 \rangle \langle 0 | + | 1 \rangle \langle 1 | \right) \\
&= \frac{1}{2} \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}
\end{aligned}

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

また、


\begin{aligned}
\mathrm{tr} \left( \rho^{2} \right) &= \frac{1}{4} \begin{pmatrix} 2 & 2 \\ 2 & 2 \end{pmatrix} \\
&= 1
\end{aligned}

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

\displaystyle{
\mathrm{tr} \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()

混合状態

測定対象となる状態は

\displaystyle{
| \psi_0 \rangle = | 0 \rangle
}

\displaystyle{
| \psi_1 \rangle = | 1 \rangle
}

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


\begin{aligned}
\rho &=  \frac{1}{2} | \psi_0 \rangle \langle \psi_0 | + \frac{1}{2} | \psi_1 \rangle \langle \psi_1 | \\
&= \frac{1}{2} \left( | 0 \rangle \langle 0| + | 1 \rangle \langle 1 | \right) \\
&= \frac{1}{2} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}
\end{aligned}

また、


\begin{aligned}
\mathrm{tr} \left( \rho^{2} \right) &= \frac{1}{4} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \\
&= \frac{1}{2}
\end{aligned}

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

\displaystyle{
\mathrm{tr} \left( \rho^{2} \right) < 1
}

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

これをZ測定します。
こちらも対角要素が各状態の確率に等しく、 |0 \rangle | 1 \rangleの確率はそれぞれ 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測定してみます。

重ね合わせ状態

測定対象となる状態は

\displaystyle{
| \psi \rangle = \frac{1}{\sqrt{2}} \left( | 0 \rangle + | 1 \rangle \right) \equiv | + \rangle
}

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

\displaystyle{
| + \rangle \equiv 
\begin{pmatrix} 1 \\ 0 \end{pmatrix}
}
\displaystyle{
| - \rangle \equiv 
\begin{pmatrix} 0 \\ 1 \end{pmatrix}
}

とすると、密度演算子


\begin{aligned}
\rho &= | \psi \rangle \langle \psi | \\
&= | + \rangle \langle + | \\
&= \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}
\end{aligned}

 | + \rangle状態である確率100%です。

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

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

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

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

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

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

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基底で記述すると以下になります。

\displaystyle{
| \psi_0 \rangle = | 0 \rangle = \frac{1}{\sqrt{2}} \left( | + \rangle + | - \rangle \right)
}

\displaystyle{
| \psi_1 \rangle = | 1 \rangle = \frac{1}{\sqrt{2}} \left( | + \rangle - | - \rangle \right)
}

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


\begin{aligned}
\rho &=  \frac{1}{2} | \psi_0 \rangle \langle \psi_0 | + \frac{1}{2} | \psi_1 \rangle \langle \psi_1 | \\
&= \frac{1}{2} \left\{ \frac{1}{2} \left( | + \rangle + | - \rangle \right) \left( \langle + | + \langle - | \right) + \frac{1}{2} \left( | + \rangle - | - \rangle \right) \left( \langle + | - \langle - | \right) \right\} \\
&= \frac{1}{4} \left\{ \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} + \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} + \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix} \right\} \\
&+  \frac{1}{4} \left\{ \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} - \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} - \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} + \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix} \right\} \\
&= \frac{1}{2} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}
\end{aligned}

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

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

図7.  | 0 \rangle | 1 \rangleの混合状態をX測定する場合

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

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

ということで、密度演算子から考えた通り、 | + \rangle | - \rangle状態が約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の複合系がある時の全体系の密度演算子 \rho^{AB}とします。
例えば、二つの量子ビットからなる一つの量子系でも良いです。

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

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

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

\displaystyle{
\rho^{A} = \mathrm{tr}_B \left( \rho^{AB} \right)
}

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

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

\displaystyle{
\mathrm{tr}_B \left( | a_i \rangle \langle a_j | \otimes | b_k \rangle \langle b_l | \right) = | a_i \rangle \langle a_j | \mathrm{tr} \left( | b_k \rangle \langle b_l | \right)
}

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


\begin{aligned}
\rho^{A} &= \mathrm{tr}_B \left( \rho^{A} \otimes \rho^{B} \right) \\
&= \rho^{A} \mathrm{tr} \left( \rho^{B} \right) \\
&= \rho^{A} \left( \because \mathrm(tr) \left( \rho^{B} \right) = 1 \right)
\end{aligned}

Bell状態を使った具体例

Bell状態の一つ

\displaystyle{
| \psi \rangle = \frac{1}{\sqrt{2}} \left( | 00 \rangle + | 11 \rangle \right)
}

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


\begin{aligned}
\rho^{AB} &= \frac{1}{2} \left( | 00 \rangle + | 11 \rangle \right) \left( \langle 00 | + \langle 11 | \right) \\
&= \frac{1}{2} \left( | 00 \rangle \langle 00 | + | 11 \rangle \langle 00 | + | 00 \rangle \langle 11 | + | 11 \rangle \langle 11 | \right) \\
\end{aligned}

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


\begin{aligned}
\rho^{A} &= \mathrm{tr}_B \left( \rho^{AB} \right) \\
&= \frac{1}{2} \left\{ \mathrm{tr}_B \left( | 00 \rangle \langle 00 |  \right) + \mathrm{tr}_B \left( | 11 \rangle \langle 00 | \right) + \mathrm{tr}_B \left( | 00 \rangle \langle 11 | \right) + \mathrm{tr}_B \left( | 11 \rangle \langle 11 | \right) \right\} \\
&= \frac{1}{2} \left\{ \mathrm{tr}_B \left( | 0 \rangle \langle 0 | \otimes | 0 \rangle \langle 0 | \right) + \mathrm{tr}_B \left( | 1 \rangle \langle 0 | \otimes | 1 \rangle \langle 0 | \right) + \mathrm{tr}_B \left( | 0 \rangle \langle 1 | \otimes | 0 \rangle \langle 1 | \right)  + \mathrm{tr}_B \left( | 1 \rangle \langle 1 | \otimes | 1 \rangle \langle 1 | \right)\right\} \\
&= \frac{1}{2} \left\{ | 0 \rangle \langle 0 | \mathrm{tr} \left( | 0 \rangle \langle 0 | \right) + | 1 \rangle \langle 0 | \mathrm{tr} \left( | 1 \rangle \langle 0 | \right) + | 0 \rangle \langle 1 | \mathrm{tr} \left( | 0 \rangle \langle 1 | \right) + | 1 \rangle \langle 1 | \mathrm{tr} \left( | 1 \rangle \langle 1 | \right) \right\} \\
&= \frac{1}{2} \left( | 0 \rangle \langle 0 | \langle 0 | 0 \rangle + | 1 \rangle \langle 0 | \langle 0 | 1 \rangle + | 0 \rangle \langle 1 | \langle 1 | 0 \rangle + | 1 \rangle \langle 1 | \langle 1 | 1 \rangle \right) \\
&= \frac{1}{2} \left( | 0 \rangle \langle 0 | + | 1 \rangle \langle 1 | \right)
\end{aligned}

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

\displaystyle{
\mathrm{tr} \left\{ \left( \rho^{A} \right)^2 \right\} = \frac{1}{2}
}

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

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

まとめ

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

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

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

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

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

sakuma.wraptas.site

参考文献

物理量の期待値

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

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

対象とする読者

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

対象でない読者

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

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

前提とする知識

バージョン情報

目次

物理量の演算子

物理量

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

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

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

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

f:id:sakumadaisuke32:20211012001717p:plain
図1. 超伝導磁束量子ビットのエネルギー測定

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

演算子の定義

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

\displaystyle{
\hat{A} = \sum_{i}{a_i | a_i \rangle \langle a_i |}
}

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

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

\displaystyle{
\hat{\sigma_z} = | 0 \rangle \langle 0 | - | 1 \rangle \langle 1 |
}

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

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

パウリ行列

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


\begin{aligned}
\sigma_z &=  | 0 \rangle \langle 0 | - | 1 \rangle \langle 1 | = 
\begin{pmatrix}
1 & 0 \\
0 & -1 \\
\end{pmatrix} (Z基底の場合) \\
\end{aligned}

\begin{aligned}
\sigma_x &= 
\begin{pmatrix}
0 & 1 \\
1 & 0 \\
\end{pmatrix} 
= | 0 \rangle \langle 1 | + | 1 \rangle \langle 0 | (Z基底の場合) \\
&= | + \rangle \langle + | - | - \rangle \langle - | (X基底の場合) \\
\end{aligned}

ここで | \pm \rangle = (1 / \sqrt{2})(| 0 \rangle \pm | 1 \rangle)


\begin{aligned}
\sigma_y &= 
\begin{pmatrix}
0 & i \\
-i & 0 \\
\end{pmatrix} 
= i | 0 \rangle \langle 1 | -i | 1 \rangle \langle 0 | (Z基底の場合) \\
&= | \frac{\pi}{2} \rangle \langle \frac{\pi}{2} | - | - \frac{\pi}{2} \rangle \langle - \frac{\pi}{2} | (Y基底の場合) \\
\end{aligned}

ここで | \pm \pi / 2 \rangle = (1 / \sqrt{2})(| 0 \rangle + e^{\pm i \pi / 2} | 1 \rangle)

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

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

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

物理量 \hat{A}の期待値 \langle \hat{A} \rangleは任意の量子状態 | \psi \rangleを用いて

\displaystyle{
\langle \hat{A} \rangle = \langle \psi | \hat{A} | \psi \rangle
}

例えば、 \sigma_zについて、状態が | \psi \rangle = \alpha | 0 \rangle + \beta | 1 \rangleの場合における期待値 \langle \sigma_z \rangleは下記です。


\begin{aligned}
\langle \sigma_z \rangle &= ( \alpha^{*} \langle 0 | + \beta^{*} \langle 1 |) (| 0 \rangle \langle 0 | - | 1 \rangle \langle 1 |) (\alpha | 0 \rangle + \beta | 1 \rangle) \\
&= | \alpha | ^2 - | \beta | ^2
\end{aligned}

期待値計算の例題

単一の物理量

 \langle \sigma_x \rangle

  • 量子状態が | \psi \rangle = | 0 \rangleの場合

\begin{aligned}
\langle \sigma_x \rangle &= \langle 0 | (| + \rangle \langle + | - | - \rangle \langle - |) | 0 \rangle \\
&= \frac{1}{2} (\langle + | + \langle - |) (| + \rangle \langle + | - | - \rangle \langle - |) (| + \rangle + | - \rangle ) \\
&= 0
\end{aligned}

量子コンピュータ(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)

f:id:sakumadaisuke32:20211011091603p:plain
図2.  | 0 \rangleをX基底で測定した結果 (左) 量子回路 (右) 測定結果
ヒストグラム上では | 0 \rangle, | 1 \rangleと表記されているが、 | + \rangle, | - \rangleの分布を観察しているのと等価である。

そして得られた分布から

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 \rangle = | \pm \rangleの場合

\begin{aligned}
\langle \sigma_x \rangle &= \langle \pm | (| + \rangle \langle + | - | - \rangle \langle - |) | \pm \rangle \\
&= \pm 1
\end{aligned}

量子コンピュータ(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)

f:id:sakumadaisuke32:20211011092522p:plain
図3.  | + \rangle を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つの場合の \langle \sigma_{1z} \sigma_{2y} \rangle
(1つ目の量子ビット \sigma_zと2つ目の量子ビット \sigma_yの積の期待値)

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

  • 量子状態が | \psi \rangle = | 0 \frac{\pi}{2} \rangleの場合

\begin{aligned}
\langle \sigma_{1z} \sigma_{2y} \rangle &= \langle 0 \frac{\pi}{2} | (| 0 \frac{\pi}{2} \rangle \langle 0 \frac{\pi}{2} | + | 1 \frac{- \pi}{2} \rangle \langle 1 \frac{- \pi}{2}) | - | 0 \frac{- \pi}{2} \rangle \langle 0 \frac{ - \pi}{2} | - | 1 \frac{\pi}{2} \rangle \langle 1 \frac{\pi}{2} |) | 0 \frac{\pi}{2} \rangle \\
&= 1
\end{aligned}

で、量子コンピュータ(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

f:id:sakumadaisuke32:20211011230631p:plain
図4.  | \psi \rangle = | 0 \frac{\pi}{2} \rangleをそれぞれZ基底とY基底で測定した結果

  • 量子状態が | \psi \rangle = \frac{1}{\sqrt{2}}(| 0 \frac{\pi}{2} \rangle + | 1 \frac{- \pi}{2} \rangleの場合

\begin{aligned}
\langle \sigma_{1z} \sigma_{2y} \rangle &= \frac{1}{2} (\langle 0 \frac{\pi}{2} | + \langle 1 \frac{- \pi}{2} | ) (| 0 \frac{\pi}{2} \rangle \langle 0 \frac{\pi}{2} | + | 1 \frac{- \pi}{2} \rangle \langle 1 \frac{- \pi}{2}) | - | 0 \frac{- \pi}{2} \rangle \langle 0 \frac{ - \pi}{2} | - | 1 \frac{\pi}{2} \rangle \langle 1 \frac{\pi}{2} |) (| 0 \frac{\pi}{2} \rangle + | 1 \frac{- \pi}{2} \rangle ) \\
&= 1
\end{aligned}

で、量子コンピュータ(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

f:id:sakumadaisuke32:20211011234600p:plain
図5.  \frac{1}{\sqrt{2}}(| 0 \frac{\pi}{2} \rangle + | 1 \frac{- \pi}{2} \rangleを測定した結果

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

まとめ

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


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

射影測定

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

対象とする読者

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

前提とする知識

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

バージョン情報

目次

測定理論

量子測定の仮定

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

Postulate3: Quantum measurements are described by a collection  \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 \rangle immediately before the measurement then the probability that result m occurs is given by
\displaystyle{
p(m) = \langle \psi | M^{\dagger}_m M_m | \psi \rangle , \tag{$2.92$}
}
and the state of the system after the measurement is
\displaystyle{
\frac{M_m | \psi \rangle}{\sqrt{\langle \psi | M^{\dagger}_m M_m | \psi \rangle}} . \tag{$2.93$}
}
The measurement operators satisfy the completeness equation,
\displaystyle{
\sum_m{M^{\dagger}_m M_m} = I. \tag{$2.94$}
}

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

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

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

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

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

一般の量子測定

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

f:id:sakumadaisuke32:20210711102958p:plain
図1. 量子測定の全体像

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

測定とは

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

測定演算子

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

\displaystyle{
\sum_m{M^{\dagger}_m M_m} = I \tag{$2.94$}
}

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

測定前の状態を | \psi \rangleとし、測定後の状態を | \psi ' \rangleとした場合に、測定後の状態は測定前の状態と測定演算子を用いて下記で表します。

\displaystyle{
| \psi ' \rangle = \frac{M_m | \psi \rangle}{\sqrt{\langle \psi | M^{\dagger}_m M_m | \psi \rangle}}
}

射影測定

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

  • X測定( \sigma_x測定)
    • 量子状態が | + \rangle | - \rangleかを測定
  • Y測定( \sigma_y測定)
    • 光の偏光が | \nearrow \rangle | \swarrow \rangleかを測定(本当は両矢印なのですが、tex記号に無い、、)
      (一般の量子状態で簡潔に表せる状態はなく、光の偏光状態くらいが特別な表現があります。)
  • Z測定( \sigma_z測定)
    • 量子状態が | 0 \rangle | 1 \rangleかを測定

通常は | 0 \rangle | 1 \rangleかを求めるZ測定がほとんどだと思います。
ただ、論文とかで式を追っていったときに \langle \sigma_x \rangle, \langle \sigma_y \rangle, \langle \sigma_z \rangleが出てきたらやってることは射影測定だということを思い出していただければと思います。

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

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

arxiv.org

arxiv.org

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

計算基底(Computational Basis)

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

射影測定とは

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

  1.  \sum_m \hat{P}_m = \hat{I}
  2.  \hat{P}_m \hat{P}_{m'} = \delta_{mm'} \hat{P}_m ( \delta_{mm'}はKroneckerのdelta)

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

\displaystyle{
\hat{P}_m = | m \rangle \langle m |
}

が用いられます。

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

量子回路での対応

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

Z測定

 | 0 \rangle, | 1 \rangle, | + \rangle, | - \rangleの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()

f:id:sakumadaisuke32:20210822231441p:plain
図2.  | 0 \rangleを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()

f:id:sakumadaisuke32:20210822232001p:plain
図3.  | 1 \rangleをZ基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

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

| + >を測定

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()

f:id:sakumadaisuke32:20210822232415p:plain
図4.  | + \rangleをZ基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

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

| - >を測定

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()

f:id:sakumadaisuke32:20210822232846p:plain
図5.  | - \rangleをZ基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

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

X測定

 | 0 \rangle, | 1 \rangle, | + \rangle, | - \rangleの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()

f:id:sakumadaisuke32:20210825212015p:plain
図6.  | 0 \rangleを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()

f:id:sakumadaisuke32:20210825212543p:plain
図7.  | 1 \rangleをX基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

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

| + >を測定

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()

f:id:sakumadaisuke32:20210825213059p:plain
図8.  | + \rangleを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()

f:id:sakumadaisuke32:20210825213504p:plain
図9.  | - \rangleをX基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

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

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

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

Y測定

 | + \rangle + | - \rangle, | + \rangle - | - \rangleの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()

f:id:sakumadaisuke32:20210825223634p:plain
図10.  | + \rangle + | - \rangleを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()

f:id:sakumadaisuke32:20210825224737p:plain
図11.  | + \rangle + | - \rangleを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()

f:id:sakumadaisuke32:20210825225441p:plain
図12.  | + \rangle + | - \rangleを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()

f:id:sakumadaisuke32:20210825225645p:plain
図13.  | + \rangle - | - \rangleをY基底で測定する場合.
(左) 測定対象の状態を示すBloch球. (中) 測定まで含めた量子回路. (右) 測定結果.

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

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

   

まとめ

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

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

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

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

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

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

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

pennylane.ai

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

pennylane.ai

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

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

対象とする読者

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

前提とする知識

バージョン情報

  • Python 3.9.5
  • PennyLane 0.15.1

目次

2クラス変分分類器

変分分類器とは

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

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

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

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

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

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

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

f:id:sakumadaisuke32:20210704102716p:plain
図1. 量子回路学習の概念図(こちらより転載)

 V(\mathbf{x})で特徴量を量子ビット上にエンコーディングし、 U(\mathbf{\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(\mathbf{x})のこと)なので一旦説明を省きます。)

f:id:sakumadaisuke32:20210704110604p:plain
図2. 2クラス分類の変分量子回路

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

pennylane.readthedocs.io

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

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

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

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

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

www.saiensu.co.jp

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

振幅エンコーディング

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

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

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

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

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

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

誤差関数

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

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が変分量子回路 \hat{U}(w_1, w_2, w_3), \hat{U}(w_4, w_5, w_6) x 6層分ののパラメータ w_1, w_2, w_3, \dots w_{36}です。

分類結果

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

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

まとめ

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