はじめに
この記事では、前回と同じくQulacsを使って、量子位相推定(Quantum Phase Estimation, QPE)を試します。
QPEとは、ユニタリー行列の固有値を量子ビット上に得るアルゴリズムです。より正確には、ユニタリー行列の固有値は絶対値が1なので、その偏角(argument)が得られれば、固有値を得たことになります。QPEはその2進展開を量子ビット上に実現します。
今回はエルミート行列 $A\in\mathbb{R}^{d\times d}$ に対して、$e^{-itA}$ の固有値を求める、というのを試したいと思います。ここで、
u\colon \mathbb{R}\ni t\mapsto e^{-itA}u(0)\in\mathbb{R}^d
は、シュレディンガー方程式
i\frac{du}{dt}=Au, \quad t\in\mathbb{R}
を満たします。
さて、なぜこれをやるのかというと、$A$ の固有展開
A=\sum_j \lambda_j P_j
に対し、
e^{-itA}=\sum_j e^{-it\lambda_j}P_j
であり1、$e^{-itA}$ の固有値の偏角を量子ビット上に得るということは、それはすなわち $A$ の固有値そのものを得るということになるからです。このことは、線型方程式を量子的に解くアルゴリズムとして知られているHHLアルゴリズムにおいて重要な役割を果たします。それも見据えてこの部分を動かしてみようというのが目的です。
今回の設定
今回は簡単にして、2次元のエルミート行列
A=
\begin{pmatrix}
0 & -i\frac{3\pi}{4} \\
i\frac{3\pi}{4} & 0
\end{pmatrix}
で考えます。固有値と固有ベクトルは対応順にそれぞれ $\pm\frac{3\pi}{4}$、$\frac{1}{\sqrt{2}}(1,\pm i)^\text{T}=\frac{|0\rangle\pm i|1\rangle}{\sqrt{2}}$ です。ただし、通常よく用いられるように、
|0\rangle :=
\begin{pmatrix}
1 \\
0
\end{pmatrix}, \quad
|1\rangle :=
\begin{pmatrix}
0 \\
1
\end{pmatrix}
です。このとき、
U:=\left. e^{-itA} \right|_{t=1}=\frac{1}{\sqrt{2}}
\begin{pmatrix}
-1 & -1 \\
1 & -1
\end{pmatrix}
に対してQPEを行います。
$U$ の固有値は $e^{-i3\pi/4}=e^{i2\pi 0.101}$ と $e^{i3\pi/4}=e^{i2\pi 0.011}$ です。ただし、小数は2進小数です。よって、レジスタービット数を3とすれば、QPEでそれぞれ $|101\rangle$ と $|011\rangle$が(確率1で)得られるはずです。
Qulacsで使える基本的なゲートについて
ここでは、Qulacsの基本的なゲートについて、備忘録的に記載しておきたいと思います。
- アダマールゲート
\begin{align*}
H &:= \frac{1}{\sqrt{2}}
\begin{pmatrix}
1 & 1 \\
1 & -1
\end{pmatrix}
\end{align*}
- パウリゲート
\begin{align*}
X &:=
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix}, \quad
Y :=
\begin{pmatrix}
0 & -i \\
i & 0
\end{pmatrix}, \quad
Z :=
\begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix}
\end{align*}
- ブロッホ球上の回転
Qulacsの動作を見ていると、通常2の $R_X,R_Y,R_Z$ ゲートとは異なり、角度 $\theta$ に対して逆に定義されているようです。
\begin{align*}
R_X(\theta)
&:= e^{i\frac{\theta}{2}X} =
\begin{pmatrix}
\cos\frac{\theta}{2} & i\sin\frac{\theta}{2} \\
i\sin\frac{\theta}{2} & \cos\frac{\theta}{2}
\end{pmatrix}, \\[.5em]
R_Y(\theta)
&:= e^{i\frac{\theta}{2}Y} =
\begin{pmatrix}
\cos\frac{\theta}{2} & \sin\frac{\theta}{2} \\
-\sin\frac{\theta}{2} & \cos\frac{\theta}{2}
\end{pmatrix}, \\[.5em]
R_Z(\theta)
&:= e^{i\frac{\theta}{2}Z} =
\begin{pmatrix}
e^{i\frac{\theta}{2}} & 0 \\
0 & e^{-i\frac{\theta}{2}}
\end{pmatrix}
\end{align*}
QPE回路図
量子回路図は次のようになります。
I-QFTは、逆量子フーリエ変換 (Inverse Quantum Fourier Transform) です。
ただし、
\begin{align*}
|\psi\rangle_\pm = \frac{|0\rangle\pm i|1\rangle}{\sqrt{2}}
\end{align*}
であり、
\begin{align*}
R^*_1 &= e^{-i\frac{\pi}{4}}R_Z\left( \frac{\pi}{2} \right), \quad R^*_2 = e^{-i\frac{\pi}{8}}R_Z\left( \frac{\pi}{4} \right)
\end{align*}
です。
上の回路図では 初期状態として $|\psi\rangle_\pm$ を描きましたが、プログラムではこれを用意してやらねばなりません。$R_X(\pm\pi/2)|0\rangle=|\psi\rangle_\pm$なので、初期状態を $|0000\rangle$ として、$R_X(\pm\pi/2)$ ゲートを入れて $|000\rangle \otimes |\psi\rangle_\pm$ を作ります。
プログラム
勉強のため、for文等使わず、回路図の通り愚直に書いています。
import numpy as np
from qulacs import QuantumState
from qulacs import QuantumCircuit
from qulacs import Observable
from qulacs.gate import DenseMatrix
from qulacs.gate import H
from qulacs.gate import RX
from qulacs.gate import RZ
from qulacs.gate import merge
from qulacs.gate import to_matrix_gate
# 量子状態の設定
# 量子ビット数を3(固有値用)+1(固有ベクトル用)=4で
n_qbits = 4
state = QuantumState(n_qbits)
state.set_zero_state()
# 量子回路
circuit = QuantumCircuit(n_qbits)
# A = [[0,-i3π/4],[i3π/4,0]]
# からユニタリー U = e^{-iA}とする
sqr2 = np.sqrt(2)
U = [
[-1.0 / sqr2, -1.0 / sqr2],
[+1.0 / sqr2, -1.0 / sqr2]
]
# ========== 前半QPE部 ==========
# |000>⊗|ψ>を準備。|ψ>はUの固有状態。
# H⊗H⊗H⊗RX(±π/2)|0000>=(H⊗H⊗H|000>)⊗|ψ>
# *** U=e^{-iA}の固有値と固有ベクトル情報 ***
# 固有状態:RX(+π/2)|0>=(|0>+i|1>)/√2, 固有値:exp(-i3π/4)=exp(i5π/4)→0.101
# 固有状態:RX(-π/2)|0>=(|0>-i|1>)/√2, 固有値:exp(+i3π/4)→0.011
# * RX(Θ)=exp(iΘX/2)
gate_list = [H(1), H(2), H(3)]
angle = np.pi / 2.0 # (☆) ±π/2
gate_list.append(RX(0, angle))
# controlled-U達のゲート
ctrl_U1_gate = DenseMatrix(0, U)
ctrl_U1_gate.add_control_qubit(3, 1)
ctrl_U2_gate = DenseMatrix(0, U)
ctrl_U2_gate.add_control_qubit(2, 1)
ctrl_U3_gate = DenseMatrix(0, U)
ctrl_U3_gate.add_control_qubit(1, 1)
# controlled gatesをリストに登録
# 1回
gate_list.append(ctrl_U1_gate)
# 2回
gate_list.append(ctrl_U2_gate)
gate_list.append(ctrl_U2_gate)
# 4回
gate_list.append(ctrl_U3_gate)
gate_list.append(ctrl_U3_gate)
gate_list.append(ctrl_U3_gate)
gate_list.append(ctrl_U3_gate)
# ========== 後半、I-QFT部 ==========
# 量子ビット1にHゲート追加
gate_list.append(H(1))
# 量子ビット2に制御R_1^*を追加、Hゲート追加
ctrl_adjoint_R1_gate = RZ(2, np.pi / 2.0)
ctrl_adjoint_R1_gate = to_matrix_gate(ctrl_adjoint_R1_gate)
ctrl_adjoint_R1_gate.multiply_scalar(np.exp(-1j * np.pi / 4.0))
ctrl_adjoint_R1_gate.add_control_qubit(1, 1)
gate_list.append(ctrl_adjoint_R1_gate)
gate_list.append(H(2))
# 量子ビット3に制御R_2^*、制御R_3^*、Hゲート追加
ctrl_adjoint_R2_gate = RZ(3, np.pi / 4.0)
ctrl_adjoint_R2_gate = to_matrix_gate(ctrl_adjoint_R2_gate)
ctrl_adjoint_R2_gate.multiply_scalar(np.exp(-1j * np.pi / 8.0))
ctrl_adjoint_R2_gate.add_control_qubit(1, 1)
ctrl_adjoint_R3_gate = RZ(3, np.pi / 2.0)
ctrl_adjoint_R3_gate = to_matrix_gate(ctrl_adjoint_R3_gate)
ctrl_adjoint_R3_gate.multiply_scalar(np.exp(-1j * np.pi / 4.0))
ctrl_adjoint_R3_gate.add_control_qubit(2, 1)
gate_list.append(ctrl_adjoint_R2_gate)
gate_list.append(ctrl_adjoint_R3_gate)
gate_list.append(H(3))
# ゲートをマージして回路に追加
merged_gate = merge(gate_list)
circuit.add_gate(merged_gate)
# 状態の更新
circuit.update_quantum_state(state)
# 状態のベクトル表示
print(state.get_vector())
# 量子ビット1～3と固有ベクトル|ψ>=(|0>±i|1>)/2の測定
coeff = 1.0
observable3 = Observable(n_qbits)
observable2 = Observable(n_qbits)
observable1 = Observable(n_qbits)
observable_eigen = Observable(n_qbits)
observable3.add_operator(coeff, 'Z 3')
observable2.add_operator(coeff, 'Z 2')
observable1.add_operator(coeff, 'Z 1')
observable_eigen.add_operator(coeff, 'Y 0')
value3 = observable3.get_expectation_value(state)
value2 = observable2.get_expectation_value(state)
value1 = observable1.get_expectation_value(state)
value_eigen = observable_eigen.get_expectation_value(state)
print('3:', value3)
print('2:', value2)
print('1:', value1)
print('eigen state:', value_eigen)
結果
- 初期状態を $|000\rangle\otimes|\psi\rangle_+$ としたとき
固有値は $e^{-i3\pi/4}=e^{i5\pi/4}$ なので、0.101、つまり $|101\rangle$ が得られるはずです。
[ 2.77555756e-17-1.38777878e-17j 1.38777878e-17+2.77555756e-17j
0.00000000e+00+1.96261557e-17j -1.96261557e-17+0.00000000e+00j
2.77555756e-17+1.38777878e-17j -1.38777878e-17+2.77555756e-17j
1.96261557e-17+5.55111512e-17j -5.55111512e-17+1.96261557e-17j
2.77555756e-17+6.93889390e-17j -6.93889390e-17+2.77555756e-17j
7.07106781e-01-1.96261557e-17j 1.96261557e-17+7.07106781e-01j
2.77555756e-17-6.93889390e-17j 6.93889390e-17+2.77555756e-17j
1.96261557e-17-5.55111512e-17j 5.55111512e-17+1.96261557e-17j]
3: -0.9999999999999989
2: 0.9999999999999989
1: -0.9999999999999989
eigen state: 0.9999999999999989
この測定はオブザーバブルの固有値が得られるので、$1$ は $|0\rangle$、$-1$ は $|1\rangle$ に対応します。なので、ちゃんと $|101\rangle$ が得られてることが分かります。また、初期状態もちゃんと設定されていることが分かります。
- 初期状態を $|000\rangle\otimes|\psi\rangle_-$ としたとき
固有値は $e^{i3\pi/4}$ なので、0.011、つまり $|011\rangle$ が得られるはずです。計算はプログラムの(☆)部分に"-"を付ければできます。
[ 2.77555756e-17+1.38777878e-17j 1.38777878e-17-2.77555756e-17j
1.96261557e-17+5.55111512e-17j 5.55111512e-17-1.96261557e-17j
2.77555756e-17+6.93889390e-17j 6.93889390e-17-2.77555756e-17j
7.07106781e-01-1.96261557e-17j -1.96261557e-17-7.07106781e-01j
2.77555756e-17-6.93889390e-17j -6.93889390e-17-2.77555756e-17j
1.96261557e-17-5.55111512e-17j -5.55111512e-17-1.96261557e-17j
2.77555756e-17-1.38777878e-17j -1.38777878e-17-2.77555756e-17j
0.00000000e+00+1.96261557e-17j 1.96261557e-17+0.00000000e+00j]
3: 0.9999999999999996
2: -0.9999999999999996
1: -0.9999999999999996
eigen state: -0.9999999999999996
ちゃんと、$|011\rangle$ になっています。初期状態もちゃんと $|\psi\rangle_-$ ですね。
まとめ
今回は簡単な例でQPEを試してみました。結論としてはちゃんと理論通りに固有値の2進小数が得られました。もっと桁が欲しい場合はレジスター量子ビットを増やして同様にすれば、固有値を得られます。
さて、実際の問題において $A$ のエルミート性を仮定しても一般性を失わない状況を考えます。正定値あるいは負定値であることが分かっていないときは、固有値は正負入り乱れている状況です。今回の結果を見ると分かりますが、得られる固有値の情報は $2\pi$ で規格化した偏角 $\theta\in [0,1)$ です。このアルゴリズムだと元の $A$ の固有値の符号までは区別がつきません。今後調べる必要がありそうです。