量子振幅推定(QAE)とは？
量子振幅推定(Quantum Amplitude Estimation :QAE)は、所望の量子状態の振幅を効率的に推定できるアルゴリズムです。
ふつう、量子状態の振幅（の絶対値）を知りたいならば、その量子状態をサンプリングしまくればよいです。
（振幅と観測確率は比例するため）
ただ、この方法だと精度よく推定するのにはサンプリング回数がたくさんいるので、QAEで賢く推定するのが良いわけです。
先の記事からQAEについて引用します。
QAEは、QAAで使うゲート$Q$とQPE(Quantum Phase Estimation) を組み合わせた発展的アルゴリズムです。
QAEは、ゲート$Q$を用いますが、振幅増幅が目的ではありません。
QAEでは、初期状態の偏角$\theta$、つまり初期状態に含まれる解状態の重み$\sin\theta$が未知である状況を扱います。
$Q$とは、初期状態の偏角を$\theta$としたときに、その偏角を$2\theta$回す回転ゲートのようなイメージでした。
実は、ここから$Q$の固有値の位相は、$2\theta$に等しいという性質が導かれます。
よって、ゲート$Q$が用意できれば、そこへQPEを適応して固有値の位相を読み出すことで$\theta$が分かります。
これで、初期状態に含まれていた解状態の重み$\sin\theta$が推定できました。
振幅を推定するので、amplitude estimation と呼ばれます。
ただ、よく考えると
自分で準備した初期状態なのに、解状態の重み(係数$\sin\theta$)は知らないという、一見して不思議な設定になっています。
不思議ですが、量子数値積分アルゴリズムではこのような状況が出現します。
$f(x)$はしっている(ゲートで書ける)が$ \int f(x) dx$は知らない ということがこれに対応します。
QAEの例題
論文の問題設定を流用します。
Practical numerical integration on NISQ devices
量子状態
$|\psi \rangle = \sqrt{\frac{3}{4}}|\psi_{0} \rangle + \sqrt{\frac{1}{4}}|\psi_{1} \rangle$
があるとします。
ここで
$|\psi_{0} \rangle = \frac{1}{\sqrt{3}}(|00 \rangle + |10 \rangle + |11 \rangle ) \otimes |0 \rangle$
$|\psi_{1} \rangle = |01 \rangle \otimes |1 \rangle$
であり、$|\psi_{0} \rangle $と$|\psi_{1} \rangle$ は正規直交基底になっています。1
最も右の$1$量子ビットは補助量子ビットです。
我々はこの$|\psi \rangle$を$|000 \rangle$から効率的に作るゲート$U$は手に入れているとします。
しかし、$|\psi_{1} \rangle$の”重み”である$\sqrt{\frac{1}{4}}$はわからないという状況を考えます。2
QAEを利用して、”重み”を推定していきます。
増幅演算子$Q$を作って、$Q$へQPEをします。
ここで、事前準備として
$|\psi \rangle = \cos\theta|\psi_{0} \rangle + \sin\theta|\psi_{1} \rangle$
とおきます。我々の目的は、$\sin\theta = \sqrt{\frac{1}{4}} = 1/2$ あるいは$\theta = \pi/6$を知ることです。
増幅演算子Qの構成
QAEをするために、対象となる増幅演算子$Q$を構成します。以下が参考になります。
https://whyitsso.net/physics/quantum_mechanics/AA_AE_Grover.html
$Q$は、次のように構成されます。
$Q = -S_{\psi}S_{f}$
ここで、
$S_{f} = |\psi_{0} \rangle \langle \psi_{0}| - |\psi_{1} \rangle \langle \psi_{1}|$
は、補助量子ビットが$1$であるような状態(今は$|01 \rangle \otimes |1 \rangle$)だけ符号反転させるゲートです。
また、
$-S_{\psi} = 2|\psi \rangle \langle \psi | - I$
は、状態を初期状態$| \psi \rangle$を軸として反転させるゲートです。
参考サイト(https://whyitsso.net/physics/quantum_mechanics/AA_AE_Grover.html)によれば、更にシンプルな書き方ができて、
$-S_{\psi} = -US_{0}U^{\dagger}$
$S_{0} = I - 2|0 \rangle \langle 0|$
$|\psi \rangle = U|000 \rangle$
となります。問題説明のところで断ったように、$U$は既知です。
よって、$S_{0}$を構築すれば$-S_{\psi}$は作れます。
Uを作る
一般の場合に対応する$U$を作ることは難しいでしょうが、いまは比較的に$|\psi \rangle$が単純なので簡単に出来ます。
nqubit = 3
state_preparation = QuantumCircuit(nqubit, name='State preparation')
state_preparation.h([1,2])
state_preparation.x(2)
state_preparation.ccx(2,1,0)
state_preparation.x(2)
このようにすればよいです。なおブラケット表示で最も右の量子ビット(つまり補助量子ビット)が回路図中では$q_{0}$に相当します。
このゲートを通せば
$|\psi \rangle = \sqrt{\frac{3}{4}}|\psi_{0} \rangle + \sqrt{\frac{1}{4}}|\psi_{1} \rangle$
が作れることは、手計算でも確認できます。
Sfを作る
$S_{f}$は簡単です。補助量子ビットが$1$である状態の符号反転させればいいので、
単に補助量子ビットに$Z$ゲートを掛ければよいです。
oracle = QuantumCircuit(nqubit)
oracle.z(0)
S0を作る
$S_{0}$を作る一般的な手続きはわからないのですが、論文によると以下のような回路でいいそうです。3
Sψを作る
$-S_{\psi} = -US_{0}U^{\dagger}$
ただ、残念ながら qiskitでは "-1倍" というゲートがないようなので、実際には
$+S_{\psi} = +US_{0}U^{\dagger}$ を作っています。
追記(2021/2/11)
グローバル位相を与える機能を使うと -1倍 が一応できるようでした。後述。
追記(2021/02/25)
$Q$の符号の話は別記事を立てました。
Groverの増幅演算子における反転操作S0の符号の問題
組み合わせてQにする
$Q = -S_{\psi}S_{f}$
circ_Q = QuantumCircuit(nqubit)
circ_Q.z(0)
circ_Q.barrier()
circ_Q.x(2)
circ_Q.ccx(2,1,0)
circ_Q.x(2)
circ_Q.h([1,2])
circ_Q.barrier()
circ_Q.x([0,1,2])
circ_Q.h(0)
circ_Q.ccx(2,1,0)
circ_Q.h(0)
circ_Q.x([0,1,2])
circ_Q.barrier()
circ_Q.h([1,2])
circ_Q.x(2)
circ_Q.ccx(2,1,0)
circ_Q.x(2)
circ_Q.barrier()
やはりqiskitでは$Q = +S_{\psi}S_{f}$ を作っています。
Qの固有値
あとはQPEで$Q$の固有値を読み出すだけですが、参考までに古典ソルバーで$Q$の固有値を直接出してみましょう。
$Q$の行列表示を取るには、以下のコマンドが簡単です。
from qiskit.quantum_info.operators import Operator
Unitary = Operator(circ_Q).data
QAEの原理によれば、求めたい位相角$\theta$が、増幅演算子$Q$の固有値$e^{\pm j 2\theta}$として現れているはずです。
いま、$\theta = \pi/6$のはずですから、$Q$の固有値は
$$e^{\pm j \frac{\pi}{3}} = \frac{1}{2} \pm j\frac{\sqrt{3}}{2} = 0.5 \pm j 0.866$$
となっているはずです。
固有値をみてみると、
a_eig = np.linalg.eig(+1*Unitary)
print("固有値\n {}\n".format(a_eig[0]))
となります。それらしい固有値が第2,3固有値として出てきていますが、符号が反転しています！
この理由は、先に断った
やはりqiskitでは$Q = +S_{\psi}S_{f}$ を作っています。
です。本来の増幅演算子は$-S_{\psi}S_{f}$なのですが、qiskitでは$+S_{\psi}S_{f}$しか作れないので、
行列が符号反転しており、固有値も符号反転しているわけです。
そのため、$2\theta = \pm\frac{1}{3}\pi$ではなく$2\theta = \frac{2}{3}\pi, \frac{4}{3}\pi$として求まっているのです。4
2021/02/25追記
符号反転は$e^{j\pi}$、即ち位相を$\pi$足すことです。$2\theta = \pm\frac{1}{3}\pi + \pi = \frac{4}{3}\pi, \frac{2}{3}\pi$ となります。
とにかく、$\theta$が固有値から分かりました。
2021/02/11追記
-1倍の演算は、qiskitではグローバル位相$\pi$を与えることで表現できるようです。
例えばqiskit.aquaのgroverライブラリで増幅演算子を自動生成させた場合
このように表示され、左上に Globai Phase : $\pi$ と入ります。
この演算子の固有値を計算すると、$2\theta = \pm \frac{1}{3}\pi$が出てきます。
なお論文（および当記事）中とは$S_{0}$の部分が上下反転しているのですが、これは影響を及ぼしません。
上下反転はエルミートを取るに相当するのですが、$S_{0}$は自己共役演算子$S_{0}^{\dagger}=S_{0}$なので何も起きないのです。
他方、全く同じ回路を手打ちで作った場合は、
このようになり、ゲートの並びは全く同じですが、この演算子の固有値は$2\theta = \frac{2}{3}\pi, \frac{4}{3}\pi$と出ます。
2021/02/25追記
この回路にglobal phase を追加して符号を直したい場合は
circ.global_phase = math.pi とします。するとaquaのライブラリ回路と一致します。
QへのQPE
では古典的な固有値計算の代わりに、QPEで固有値推定を行います。
正確に言えばQPEでは固有値ではなく固有値の位相部分が求まりますので、直接$2\theta = \pm\pi/3$が出てきます。
ただ、$Q$の符号がひっくり返ってますので、$2\theta = \frac{2}{3}\pi, \frac{4}{3}\pi$ として出てきます。
また、QPEでは位相そのものではなく位相を$2\pi$で規格化したものが2進小数表示で出てきます。
規格した位相を先に計算しておくと、$2\theta_{norm.} = 0.333 , 0.666$ となります。
QPEのコードは以下のようになります。qiskitのチュートリアルを改造すればすぐ作れます。
import math
n_ancilla = 3
qpe = QuantumCircuit(nqubit+n_ancilla, n_ancilla) #
qpe.append(state_preparation,[n_ancilla,n_ancilla+1,n_ancilla+2])
for qubit in range(n_ancilla):
qpe.h(qubit)
repetitions = 1
for counting_qubit in range(n_ancilla):
for i in range(repetitions):
Arb_gate = UnitaryGate(Unitary, label='Q')
Arb_gate_ctrl = add_control(Arb_gate, num_ctrl_qubits=1, label='CQ', ctrl_state='1')
qpe.append(Arb_gate_ctrl,[counting_qubit,n_ancilla,n_ancilla+1,n_ancilla+2])
repetitions *= 2
def qft_dagger(circ, n):
"""n-qubit QFTdagger the first n qubits in circ"""
# Don't forget the Swaps!
for qubit in range(n//2):
circ.swap(qubit, n-qubit-1)
for j in range(n):
for m in range(j):
circ.cu1(-math.pi/float(2**(j-m)), m, j)
circ.h(j)
qpe.barrier()
# Apply inverse QFT
qft_dagger(qpe, n_ancilla)
# Measure
qpe.barrier()
for n in range(n_ancilla):
qpe.measure(n,n)
測定します。
shots = 2048
results = execute(qpe, backend=Aer.get_backend('qasm_simulator'), shots=shots).result()
answer = results.get_counts()
plot_histogram(answer)
ピークが2つ見えます。それぞれ$2\theta_{norm.} = 0.333 , 0.666$であろうと予測されます。
2進小数表示を10進数にしてみると(ここで、右端のビットが$2^{-3}$です)、
$ 010 \to 1/4 = 0.25$
$ 101 \to 1/2+1/8 = 0.625$
ということで、まぁまぁ合っていそうです。
推定の精度は、QPEで使う補助量子ビット数で制限されそうだということがわかります。
以上で、誤差はあるものの確かに$\theta$が推定でき、
$|\psi \rangle = \cos\theta|\psi_{0} \rangle + \sin\theta|\psi_{1} \rangle$
の振幅である$\sin\theta$がわかりました。
これがQAEです。
2021/02/25追記
符号が正しい$+Q$に対する結果も載せておきます。
$2\theta = \pm\frac{1}{3}\pi$、即ち正の偏角としては $2\theta = \frac{1}{3}\pi, \frac{5}{3}\pi$です。
$2\theta_{norm.} = 0.166 , 0.8333$ です。
$ 001 \to 1/8 = 0.125$
$ 110 \to 1/2+1/4 = 0.75$
となりますので、やはり正しく動いています。
結論
- QAEは、準備した状態に解状態が含まれている場合、その重み（振幅）を推定するアルゴリズム
- 振幅を固有値に持つ増幅演算子$Q$をうまく準備して、QPEを行うことにより抽出する
- 推定精度はQPEの読み出し精度（補助量子ビット数）の影響を受ける
なお、論文中ではQPEを回避したQAEも紹介されています。
QPEは回路が長くなるし、補助量子ビットを多数使うためです。