自己紹介
こんにちは。量子計算大好き大学生のdaveと申します。
今回は線型方程式を解くアルゴリズム[1]の解説をやっていきます。
線型方程式を解く量子アルゴリズムには、かの有名なHHLアルゴリズム[2]がありますが、説明・実装共に煩雑を極めているので、今回は量子計算の基本(量子ビット&量子ゲート)が分かれば理解できる簡単なアルゴリズムを紹介します。
また、IBMが開発しているPython量子計算フレームワークのQiskit[3]を使ったデモも行います。
連立１次方程式とは
行列Aとベクトルbの間で成立する式
$$Ax = b$$を満たす未知数ベクトルxを解に持つ方程式です。
量子アルゴリズムの解説
行列Aとベクトルbの値が以下であることを仮定します。
A = \begin{bmatrix}
a_{00} & a_{01} \\
a_{10} & a_{11}
\end{bmatrix}
x = \begin{bmatrix}
x_0 \\
x_1
\end{bmatrix}
b = \begin{bmatrix}
b_0 \\
b_1
\end{bmatrix}
また、以下の量子ゲートを定義します。
$$R(\alpha,\beta) = R_z(\beta)R_y(\alpha)R_z(- \beta)$$$$U_{ij}(\alpha,\beta) = CNOT_{ji} \cdot (R(\alpha, \beta) \otimes I)\cdot CNOT_{ij} \cdot (R(-\alpha, -\beta) \otimes I)\cdot CNOT_{ji}$$
それでは手順を説明します
1.ベクトルbを量子状態にエンコードする
式がN本あったら、量子ビットはN+1個用意します。
この場合、連立方程式は2本なので量子ビットは3個用意します。
初期状態は以下のようになります。
$$ |\Psi\rangle_b = \sqrt{1-{b_0}^2-{b_1}^2}|000\rangle + b_0|001\rangle + b_1|010\rangle \cdots (1)$$
これは以下のようにゲートを掛けると実装できます。
$$|\Psi\rangle_b = (U_{01}(\beta_2) \otimes I)(R_{y}(\beta_1) \otimes I \otimes I)|000\rangle \
= \cos\frac{\beta_1}{2}|000\rangle + \sin{\beta_2}\sin{\frac{\beta_1}{2}}|001\rangle - \cos{\beta_2}\sin{\frac{\beta_1}{2}}|010\rangle \cdots(2)$$
まずは(1)を満たすパラメータ$\beta_0、\beta_1$を考えましょう。
ここで、$\beta_1 = \pi$とすると$$\beta_1 = \arctan\Bigl({-\frac{b_0}{b_1}}\Bigr)$$になるので計算が楽になります。
を計算すれば良いことが分かります。
2. 初期状態にAの逆行列を掛ける
今回は以下のユニタリー変換を使います。
$$U_{012}(\alpha_1, \alpha_2) = U_{12}(\alpha_2, 0)U_{01}(\alpha_1,0)$$
これを初期状態に掛けると以下のような状態になります。
$$U_{012}|\Psi\rangle_b = (x_0D_{11} + x_1D_{12})|001\rangle + (x_0D_{21} + x_1D_{22})|010\rangle + (x_0D_{31} + x_1D_{32})|011\rangle $$
$$ D_1 = a_{00}\cos{\alpha_1} + a_{10}\sin{\alpha_1}$$$$ D_2 = a_{01}\cos{\alpha_1} + a_{11}\sin{\alpha_1}$$$$ D_{11} =a_{10}\cos{\alpha_1}-a_{00}\cos{\alpha_1}$$$$ D_{12} = a_{11}\cos{\alpha_1} - a_{01}\sin{\alpha_1}$$$$ D_{21} = -D_1\sin{\alpha_2}$$$$ D_{22} = -D_2\sin{\alpha_2}$$$$D_{31}=D_1\cos{\alpha_1}$$$$D_{32}=D_1\cos{\alpha_2}$$
(上の計算は現時点では手計算で追えていません 後日加筆します)
ここで、$b_i = a_{i0}x_0 + a_{i1}x_1$より、右辺第2項に着目すると、
$$D_{21} = 1, D_{22} = 0 \rightarrow \langle 010 |U_{012}|\Psi\rangle_b = x_0 \cdots(3)$$$$D_{21} = 0, D_{22} = 1 \rightarrow \langle 010 |U_{012}|\Psi\rangle_b = x_1 \cdots(4)$$
になります。
ベクトル$x$の値が求められました！
ちなみに、$(3)$の時は
$$-(a_{00}\cos{\alpha_1} + a_{10} \sin{\alpha_1}) \sin{\alpha_2}=1$$$$-(a_{01} \cos{\alpha_1} + a_{11}\sin{\alpha_1})\sin{\alpha_2}=0 $$
より
$$ \alpha_1 = \arctan \Bigl(- \frac{a_{01}}{a_{11}} \Bigr), \alpha2 = \arcsin \Bigl({- \frac{1}{a_{01}\cos{\alpha_1}+a_{11}\sin{\alpha_1}}}\Bigr)$$
$(4)$の時は
$$-(a_{00}\cos{\alpha_1} + a_{10} \sin{\alpha_1}) \sin{\alpha_2}=0$$$$-(a_{01} \cos{\alpha_1} + a_{11}\sin{\alpha_1})\sin{\alpha_2}=1 $$
より
$$ \alpha_1 = \arctan \Bigl(- \frac{a_{00}}{a_{10}} \Bigr), \alpha2 = \arcsin \Bigl({- \frac{1}{a_{00}\cos{\alpha_1}+a_{10}\sin{\alpha_1}}}\Bigr)$$
を計算すればいいです。
3. Qiskitを用いた実装例
ここまで紹介したアルゴリズムをIBMが開発した量子計算フレームワークQiskitで実装します。ここでは、論文[1]中の１番簡単な実装例をご紹介します。
A = \begin{bmatrix}
-1.8 & 0.6 \\
-0.4 & 1.4
\end{bmatrix}
b = \begin{bmatrix}
-0.6 \\
0.8
\end{bmatrix}
まずは$Ax=b$の解を古典的に求めてみましょう。
import numpy as np
A = np.array([[-1.8, 0.6], [-0.4, 1.4]])
b = np.array([-0.6, 0.8])
A_inv = np.linalg.inv(A)
x_classical = np.dot(A_inv, b)
print("x:", x_classical) # x: [0.57894737 0.73684211]
従って、求めるべきベクトルは
x = \begin{bmatrix}
0.57894737 \\
0.73684211
\end{bmatrix}
であることが分かります。
それでは、このベクトルを量子アルゴリズムを使って求めてみましょう。
from qiskit import QuantumCircuit, Aer, execute
import numpy as np
def solve_linear_equation(A, b):
# 答えのベクトル
x = np.array([])
# パラメータ計算
alpha0_for_x0 = np.arctan(- A[0][1]/A[1][1])
alpha1_for_x0 = np.arcsin(-1/(A[0][0]*np.cos(alpha0_for_x0)+A[1][0]*np.sin(alpha0_for_x0)))
alpha0_for_x1 = np.arctan(- A[0][0]/A[1][0])
alpha1_for_x1 = np.arcsin(-1/(A[0][1]*np.cos(alpha0_for_x1)+A[1][1]*np.sin(alpha0_for_x1)))
beta0 = np.pi
beta1 = np.arctan(-b[0]/b[1])
# Rゲートを定義
def R(qc, i, alpha, beta):
qc.rz(beta, i)
qc.ry(alpha, i)
qc.rz(-beta, i)
# Uゲートを定義
def U(qc, i, j, alpha, beta):
qc.cx(i, j)
R(qc, i, alpha, beta)
qc.cx(j, i)
R(qc, i, -alpha, -beta)
qc.cx(i, j)
# x0を求める量子回路
qc = QuantumCircuit(3)
qc.ry(beta0, 0)
U(qc, 0, 1, beta1, 0)
U(qc, 0, 1, alpha0_for_x0, 0)
U(qc, 1, 2, alpha1_for_x0, 0)
# バックエンドを指定 (今回は状態ベクトルを求め得るのでstatevector_simulator)
backend = Aer.get_backend('statevector_simulator')
# ジョブを実行
job = execute(qc, backend)
# 実行結果を取得
result = job.result()
# 状態ベクトルを取得
statevector = result.get_statevector(qc)
# x0の値を計算
x0 = np.dot(statevector, np.array([0,0,1,0,0,0,0,0]))
# 答えのベクトルに追加
x.append(x0.real)
# x1を求める量子回路
qc = QuantumCircuit(3)
qc.ry(beta0, 0)
U(qc, 0, 1, beta1, 0)
U(qc, 0, 1, alpha0_for_x1, 0)
U(qc, 1, 2, alpha1_for_x1, 0)
# バックエンドを指定 (今回は状態ベクトルを求め得るのでstatevector_simulator)
backend = Aer.get_backend('statevector_simulator')
# ジョブを実行
job = execute(qc, backend)
# 実行結果を取得
result = job.result()
# 状態ベクトルを取得
statevector = result.get_statevector(qc)
# x0の値を計算
x1 = np.dot(statevector, np.array([0,0,1,0,0,0,0,0]))
# 答えのベクトルに追加
x.append(x1.real)
return x
if __name__ == "__main__":
x_quantum = solve_linear_equation(A, b)
print("x:", x_quantum) # x:[0.5789473684210527, 0.7368421052631582]
それでは、数値解x_classicalと量子アルゴリズムの解x_quantumがどれくらい近いか調べましょう。ベクトルの距離はいくつかありますが、今回は2つのベクトルの差に対してL2ノルムを計算します。
任意のベクトル$a$の間のL2ノルムはベクトルの大きさのようなもので、以下のように計算されます。
\sqrt{\sum_{a_i \in a}|a_i|^2}
ノルムの厳密な定義に関してはこちらをご覧下さい。
L2ノルムはnumpyのnumpy.linalg.norm関数で計算できるので、その結果を見てみましょう。
l2norm = np.linalg.norm(x_classical-x_quantum)
print("l2norm is ", l2norm) # l2norm is 2.220446049250313e-16
よって、量子アルゴリズムで求めた解が数値解と非常に近いことが分かります。
さて、次は「量子コンピュータで線型方程式を解こう - 上級編 -」と題してHHLアルゴリズムの解説を行います。ご期待下さい！
[1]Doronin, Sergey I., E. B. Fel’dman, and Alexandre I. Zenchuk. "Solving systems of linear algebraic equations via unitary transformations on quantum processor of IBM Quantum Experience." Quantum Information Processing 19.2 (2020): 1-20.
[2]Harrow, Aram W., Avinatan Hassidim, and Seth Lloyd. "Quantum algorithm for linear systems of equations." Physical review letters 103.15 (2009): 150502.