$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$
はじめに
自作の量子計算シミュレータqlazyで、量子アルゴリズムの基本を確認するよ。今回は「アダマールテスト」です。以下を参考にさせていただきました。
アダマールテスト
ユニタリ演算子の期待値推定
あるユニタリ演算子$U$と量子状態$\ket{\psi}$があったとき、期待値$\bra{\psi}U\ket{\psi}$を求めたい場合があります。単なる行列計算なので、そんなもんやりゃいいじゃん、と思われるかもしれませんが、その次元は量子ビット数に対して指数的に増大するので、甘く見ていると痛い目を見ます。しかし、以下のような回路で定義される「アダマールテスト」によって、現実的に推定可能となることがわかっています。つまり、$U$の演算を直接相手にするのではなく、補助量子ビットを制御ビットとした制御ユニタリゲートを用意しておいて、間接的に測定することで、上手に推定できるようになるというわけです。ということで、まずはそれについて見ていきます。
|0> ---H---*---H---M---
|psi> -/-----U-----------
(ここで、|psi>の部分は1本しか線が書かれていませんが、複数本のつもりです。"-/--"は複数あるということを表しています、と読んでください。また、"U"は一番上の量子ビットを制御ビットにした何らかの制御ユニタリゲートです、というつもりです。)
上の回路図を素直に左から右に計算していきます。
\begin{align}
&\ket{0} \otimes \ket{\psi} \\
&\Rightarrow \frac{1}{\sqrt{2}} (\ket{0}+\ket{1}) \otimes \ket{\psi} \\
&\Rightarrow \frac{1}{\sqrt{2}} (\ket{0} \otimes \ket{\psi} + \ket{1} \otimes U\ket{\psi}) \\
&\Rightarrow \frac{1}{\sqrt{2}} (\ket{+} \otimes \ket{\psi} + \ket{-} \otimes U\ket{\psi}) \\
&= \frac{1}{2} \{ (\ket{0}+\ket{1}) \otimes \ket{\psi} + (\ket{0}-\ket{1}) \otimes U\ket{\psi} \} \\
&= \frac{1}{2} \{ \ket{0} \otimes (I+U)\ket{\psi} + \ket{1} \otimes (I-U)\ket{\psi} \}
\end{align}
最後の行が測定直前の状態です。測定した結果、0を得る確率は、
\begin{align}
p_{0} &= \frac{1}{4} ||(I+U)\ket{\psi}||^2 \\
&= \frac{1}{4} \bra{\psi} (I+U^{\dagger})(I+U) \ket{\psi} \\
&= \frac{1}{2} + \frac{1}{4} \bra{\psi} (U^{\dagger} + U) \ket{\psi} \\
&= \frac{1}{2} + \frac{1}{4} (\bra{\psi}U\ket{\psi}^{*} + \bra{\psi}U\ket{\psi}) \\
&= \frac{1}{2} (1 + Re \bra{\psi} U \ket{\psi}) \\
\end{align}
となります。同様に、1を得る確率は、
p_{1} = \frac{1}{2} (1 - Re \bra{\psi} U \ket{\psi})
となります。したがって、$U$の期待値の実部は、
Re \bra{\psi} U \ket{\psi}) = p_{0} - p_{1}
と計算できます。一方、虚部については、補助量子ビットに、位相ゲート
S =
\begin{pmatrix}
1 & 0 \\
0 & i
\end{pmatrix}
を付加して、以下のような回路で推定できます。
|0> ---H--S--*---H---M---
|psi> -/-------U-----------
上と同様な計算により、0,1を得る確率$p_{0}^{\prime},p_{1}^{\prime}$は、
\begin{align}
p_{0}^{\prime} &= \frac{1}{2} (1 - Im \bra{\psi} U \ket{\psi}) \\
p_{1}^{\prime} &= \frac{1}{2} (1 + Im \bra{\psi} U \ket{\psi})
\end{align}
となるので、$U$の期待値の虚部は、
Im \bra{\psi} U \ket{\psi} = p_{1}^{\prime} - p_{0}^{\prime}
となります。実部と虚部がわかったので、期待値が求まったことになります。
というわけで、本当にこれで推定できるか、qlazyで確認します。コードは以下です。対象としたユニタリ演算子は2量子ビット演算で、$U=H_{0}RX_{1}(\pi/4)$としてみました(何でも良いので適当に決めました)。
from qlazypy import QState
def hadamard_test_real():
print("== expectaion value prediction (real part) == ")
# prepare initial state
qs_0 = QState(1)
qs_psi = QState(2)
qs = qs_0.tenspro(qs_psi)
# circuit for hadamard test
qs.h(0)
qs.ch(0,1).crx(0,2,phase=0.25)
qs.h(0)
shots = 1000
md = qs.m(id=[0], shots=shots)
p0 = md.frq[0]/shots
p1 = md.frq[1]/shots
exp_pred = p0-p1
print("expectation value (predict) = {0:.3f}".format(exp_pred))
# theoretical expactation value
qs_op = qs_psi.clone()
qs_op.h(0).rx(1,phase=0.25)
exp_theo = qs_psi.inpro(qs_op).real
print("expectation value (theoretical) = {0:.3f}".format(exp_theo))
# free memory
qs_0.free()
qs_psi.free()
qs_op.free()
qs.free()
def hadamard_test_imag():
print("== expectaion value prediction (imaginary part) == ")
# prepare initial state
qs_0 = QState(1)
qs_psi = QState(2)
qs = qs_0.tenspro(qs_psi)
# circuit for hadamard test
qs.h(0).s(0)
qs.ch(0,1).crx(0,2,phase=0.25)
qs.h(0)
shots = 1000
md = qs.m(id=[0], shots=shots)
p0 = md.frq[0]/shots
p1 = md.frq[1]/shots
# exp_pred = p0-p1
exp_pred = p1-p0
print("expectation value (predict) = {0:.3f}".format(exp_pred))
# theoretical expactation value
qs_op = qs_psi.clone()
qs_op.h(0).rx(1,phase=0.25)
exp_theo = qs_psi.inpro(qs_op).imag
print("expectation value (theoretical) = {0:.3f}".format(exp_theo))
# free memory
qs_0.free()
qs_psi.free()
qs_op.free()
qs.free()
if __name__ == '__main__':
hadamard_test_real()
hadamard_test_imag()
実部と虚部を別の関数で定義していますが、ほぼ同じなので、実部の関数hadamard_test_real()についてのみ簡単に説明します。
qs_0 = QState(1)
qs_psi = QState(2)
qs = qs_0.tenspro(qs_psi)
で、補助量子ビットを$\ket{0}$、任意に用意する量子状態$\ket{\psi}$を$\ket{00}$とし、そのテンソル積の状態$\ket{0}\otimes\ket{00}$を変数qsに格納します(v0.0.14でテンソル積を実行するtensproメソッドを追加しました)。
qs.h(0)
qs.ch(0,1).crx(0,2,phase=0.25)
qs.h(0)
で、アダマールテストの回路を適用します。ここでch(0,1)は、量子ビット(0,1)に対する制御アダマールゲートで、crx(0,2,phase=0.25)は、量子ビット(0,2)に対する(回転角$\pi/4$の)制御X軸回転ゲートを表しています。これで、ユニタリゲート$U=H_{1}RX_{2}(\pi/4)$に対する制御ユニタリが実現できます(v0.014でch,crxゲートを追加しました)。
shots = 1000
md = qs.m(id=[0], shots=shots)
p0 = md.frq[0]/shots
p1 = md.frq[1]/shots
exp_pred = p0-p1
print("expectation value (predict) = {0:.3f}".format(exp_pred))
で、100O回測定したときの結果から、上に示した推定式によって期待値(実部)を計算し表示します。
qs_op = qs_psi.clone()
qs_op.h(0).rx(1,phase=0.25)
exp_theo = qs_psi.inpro(qs_op).real
print("expectation value (theoretical) = {0:.3f}".format(exp_theo))
最後のこの部分で、期待値を真面目に計算しています。対象となっている2量子状態$\ket{00}$に$H_{0}RX_{1}(\pi/4)$を適用した状態と元々の2量子状態$\ket{00}$との内積によって、期待値を得ています(cloneメソッドは状態のコピー、inproメソッドは内積です)。
実行結果は、以下です。
== expectaion value prediction (real part) ==
expectation value (predict) = 0.658
expectation value (theoretical) = 0.653
== expectaion value prediction (imaginary part) ==
expectation value (predict) = -0.064
expectation value (theoretical) = 0.000
というわけで、だいたい推定できているかと思います。
ユニタリ演算子の固有状態への収束
先程は任意の状態に対するユニタリ演算子の期待値を推定しましたが、アダマールテストの結果の状態を再度入力する操作を繰り返すことで、任意の状態からユニタリ演算子の固有状態に徐々に収束させることができるようなので、次にそれを確認します。
先程の回路を繰り返すだけなので、いきなりそのコードを示します(対象としたユニタリ演算子は同じものです)。
from qlazypy import QState
def main():
print("== convergence on engenstate ==")
# prepare initial state
qs_0 = QState(1)
# qs_psi = QState(2)
qs_psi = QState(2).h(0).h(1)
qs = qs_0.tenspro(qs_psi)
# circuit for hadamard test
iter = 100
for _ in range(iter):
qs.h(0)
qs.ch(0,1).crx(0,2,phase=0.25)
qs.h(0)
qs.m(id=[0],shots=1)
print("- final |psi>")
qs.show(id=[1,2])
print("- unitary operated |psi>")
qs.h(1).rx(2,phase=0.25).show(id=[1,2])
qs_0.free()
qs_psi.free()
qs.free()
if __name__ == '__main__':
main()
何をやっているかというと、まず状態の初期化は同じなので省略して、
iter = 100
for _ in range(iter):
qs.h(0)
qs.ch(0,1).crx(0,2,phase=0.25)
qs.h(0)
qs.m(id=[0],shots=1)
で、アダマールテストを100回繰り返します。
print("- final |psi>")
qs.show(id=[1,2])
で、(補助量子ビットを除いた)最終状態を表示します。もし、これがユニタリ演算子$H_{1}RX_{2}(\pi/4)$の固有状態であるなら、$H_{1}RX_{2}(\pi/4)$を最終状態に適用した結果は、固有値分を除き一致しているはずです。
print("- unitary operated |psi>")
qs.h(1).rx(2,phase=0.25).show(id=[1,2])
で、それを確認しています。
実行結果は、以下です。
== convergence on engenstate ==
- final |psi>
c[00] = +0.6533+0.0000*i : 0.4268 |+++++
c[01] = +0.6533+0.0000*i : 0.4268 |+++++
c[10] = +0.2706+0.0000*i : 0.0732 |++
c[11] = +0.2706+0.0000*i : 0.0732 |++
- unitary operated |psi>
c[00] = +0.6533-0.0000*i : 0.4268 |+++++
c[01] = +0.6533-0.0000*i : 0.4268 |+++++
c[10] = +0.2706+0.0000*i : 0.0732 |++
c[11] = +0.2706+0.0000*i : 0.0732 |++
というわけで、きれいに一致しました。つまり、$U$の固有状態に収束しました。qlazyのshowメソッドは状態の位相項を無視して正規化までやってしまうので、固有値がどうなっているか残念ながら表示できませんが、、。 また、何度か実行すればわかりますが、固有状態は２種類あるので、もうひとつの また、qs_psiで表される初期状態をいろいろ変えて実行すればわかりますが、固有状態は４種類あるので、上記以外の固有状態に収束することもあります。
ユニタリ演算子の固有値の位相推定
ユニタリ演算子の固有状態が用意できたところで、今度はその固有値の位相が推定できることを示してみます。ユニタリ演算子$U$に対する一つの固有状態を$\ket{\psi}$とすると、
U \ket{\psi} = e^{i\gamma} \ket{\psi}
が成り立ちます(ユニタリなので固有値は必ずこの形になります)。この前提でアダマールテストをやってみると、測定前の状態は、先程の計算の$U$を単純に$e^{i\gamma}$に置き換えれば良くて、結局、
(\frac{1+e^{i\gamma}}{2} \ket{0} + \frac{1-e^{i\gamma}}{2} \ket{1}) \otimes \ket{\psi}
となります。測定をして$m=0,1$を得る確率$p_{m}$は、
p_{m} = |\frac{1+(-1)^m e^{i\gamma}}{2}|^2 = \frac{1+(-1)^{m}cos\gamma}{2}
となります。この式から、
cos\gamma = p_{0} - p_{1}
なので、これで位相$\gamma$が求まります。
では、本当にこれで推定できるか、qlazyで試してみます。コードは以下です。
import math
from qlazypy import QState
def main():
print("== phase prediction == ")
# prepare initial state
qs_0 = QState(1)
qs_psi = QState(1).h(0)
qs = qs_0.tenspro(qs_psi)
# circuit for hadamard test
qs.h(0)
qs.crx(0,1,phase=0.5)
qs.h(0)
shots = 1000
md = qs.m(id=[0], shots=shots)
p0 = md.frq[0]/shots
p1 = md.frq[1]/shots
print("[predicted] cos(gamma) = {0:.3f}".format(p0-p1))
print("[theoretical] cos(gamma) = {0:.3f}".format(math.cos(-0.25*math.pi)))
qs_0.free()
qs_psi.free()
qs.free()
if __name__ == '__main__':
main()
わかりやすく確認するため、今度は、ユニタリ演算子$U$としてX軸周りの角度$\pi/2$の回転演算子を採用しました。その固有状態は$\ket{+},\ket{-}$ですが、とりあえず$\ket{+}$でやってみました。
\begin{align}
&U = RX(\pi/2) =
\begin{pmatrix}
cos \frac{\pi}{4} & -i sin \frac{\pi}{4} \\
-i sin \frac{\pi}{4} & cos \frac{\pi}{4}
\end{pmatrix} \\
&\ket{\psi} = \ket{+} = \frac{1}{\sqrt{2}}
\begin{pmatrix}
1 \\
1
\end{pmatrix}
\end{align}
このとき、固有値は$e^{-i\pi/4}$となるので、アダマールテストで、$cos\gamma$の推定値が$cos(-\pi/4) = 1/\sqrt{2}$になれば推定成功です。では、上のコードを簡単に説明します(と言っても、ほぼ見た通りですが)。
qs_0 = QState(1)
qs_psi = QState(1).h(0)
qs = qs_0.tenspro(qs_psi)
で、初期状態を用意します。$\ket{\psi} = \ket{+}$にするため、アダマールゲートをかけています。
qs.h(0)
qs.crx(0,1,phase=0.5)
qs.h(0)
対象としているユニタリ演算子は$\pi/2$のX軸周りの回転演算子なので、それに対応した制御ユニタリゲートをかけています(v0.0.14でcrxメソッドを追加しています)。
shots = 1000
md = qs.m(id=[0], shots=shots)
p0 = md.frq[0]/shots
p1 = md.frq[1]/shots
print("[predicted] cos(gamma) = {0:.3f}".format(p0-p1))
print("[theoretical] cos(gamma) = {0:.3f}".format(math.cos(-0.25*math.pi)))
で、1000回測定した結果に基づき、位相の$cos$の推定値を計算し、さらに$cos(-\pi/4)$を合わせて、表示しています。
実行結果は、以下です。
== phase prediction ==
[predicted] cos(gamma) = 0.710
[theoretical] cos(gamma) = 0.707
というわけで、だいたい推定できているかと思います。
おわりに
アダマールテストが一体どんなところに役に立つのかという説明もなしに、淡々と確認してみましたが、実はいろんな量子アルゴリズムのサブルーチンとして使われているとのことです。例えば、補助系をもう少し工夫した「位相推定アルゴリズム」を用いることで、素因数分解アルゴリズムができたり、量子化学計算等で登場する時間発展演算子の計算ができたりします。次回以降、このあたりの理解を進めてみたいと思います。
補足
更新履歴
2020.1.25
「ユニタリ演算子の固有状態への収束」のシミュレーション結果が間違っていたので、訂正しました。
初期状態を準備するところで、
qs_psi = QState(2)
としていましたが、これだとうまく収束してくれません。
qs_psi = QState(2).h(0).h(1)
という具合に、アダマールをかけて2量子ビットの均等重ね合わせ状態を初期状態とすることで、うまく収束してくれます。アダマールをかけなくても収束しているように見えていたのは、シミュレータのバグでした、汗。v0.0.32で修正しました。
2020.2.3
「ユニタリ演算子の期待値推定」で、関数hadamard_test_imag()中のexp_pred=p0-p1はexp_perd=p1-p0の間違いだったので、訂正しました。
2020.2.4
「ユニタリ演算子の固有状態への収束」の記述が間違っていたので訂正しました。固有状態は「２種類」ではなく「４種類」でした。また、「何度か実行すれば」どれかの固有状態に収束するという言い方をしていましたが、この例で実際にやってみいると、他の固有状態に収束するというのは、かなり稀な現象のようです。それよりも、初期状態をいろいろいじったほうが、他の固有状態への収束が容易に観測できるようなので、そのように本文の記述を変えました。
以上