HHLアルゴリズム
HHLアルゴリズムは、量子アルゴリズムであって線形方程式 $Ax=b$を解くものです。
先人の解説は、
- QiitaのHHL アルゴリズム
- qiskitのチュートリアル Solving Linear Systems of Equations using HHL
- Dojoの記事 7-2. Harrow-Hassidim-Lloyd (HHL) アルゴリズム
あたりが良いと思います。現論文は、
(参考) 発展形として変分法による解法もあります。
量子コンピューターで連立一次方程式を変分法で解く量子アルゴリズムVQLSを調査した件
こちらは、動画を観るとよくわかります。動画version
今回はqiskitのチュートリアルのほうを見ていきます。
ただ、後述するように、このチュートリアルは具体的な計算例に間違いがかなりあります。
qiskit でHHLアルゴリズム
import
#initialization
import matplotlib.pyplot as plt
import numpy as np
import math
# importing Qiskit
from qiskit import IBMQ, Aer
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, execute
# import basic plot tools
from qiskit.visualization import plot_histogram
from qiskit.quantum_info.operators import Operator, Pauli
from qiskit.circuit.add_control import add_control
from qiskit.extensions import UnitaryGate
from scipy.linalg import expm
from qiskit.quantum_info import Statevector
from qiskit.circuit.library import RYGate
問題の定義
もともとの解きたい問題は、ユニタリ行列$A$
に対して、$Ax=b$を解くことです。
です。
$U = \exp(iA)$ を定義して、$U$に対して量子位相推定(QPE)を使います。
このとき$U$の固有値の位相、すなわちそれは$A$の固有値なのですが、これが2進小数の有限桁で求まります。
今回は簡単のために、$A$の固有値を2進小数2桁で書けるようにスケーリングしておきます。
もちろんこれは$A$の固有値を知っていないとできない調整なので、ずるですが、固有値が未知でも
QPEの桁数以上の精度が出ないだけで、アルゴリズムは動きますのでご心配なく。
これには、$A$に$2\pi(3/8)$をかければいいことが（手計算で）わかります。
そのうえで、スケーリングされた$A$に対して固有値と固有ベクトルを計算してみます。
（この操作は本来のHHLアルゴリズムでは、しません。今回は動作を追っていきたいのでやっているだけです）
すると
- 第一の固有値は$\lambda_{1} = \pi = 10_{2}$ , 固有ベクトルは $u_{1} = [1, -1]^{T}$
- 第二の固有値は$\lambda_{2} = \frac{\pi}{2} = 01_{2}$ , 固有ベクトルは $u_{2} = [1, 1]^{T}$
ということがわかります。
添字の $_{2}$ は、2進小数表示です。(つまり$2\pi$で割って$2^{-n}$で展開したときの係数)
$A$をスケーリングしたことで、2桁で書けています。
また、ついでに$Ax=b$の解$x$を規格化したものも求めます。
HHLアルゴリズムでは、解$x$が規格化されたものが量子状態の振幅に埋め込まれて出てきます。
規格化係数も後で測定で求められるようですが、完全に理解できていないので今回は扱いません。
コードで書くと
Scale = 2*np.pi*(3/8)
A_wo_scaling = np.array([[1, -1/3],[-1/3, 1]])
A = np.array([[1, -1/3],[-1/3, 1]])*Scale
b = np.array([1,0])
U = expm(1j*A)
print("A with scaling = \n {:}\n".format(A))
# calc. exact classical solution
A_inv = np.linalg.inv(A)
x = np.dot(A_inv, b)
norm = np.linalg.norm(x)
x_normalized = x / norm
print("Normalized solution, x_norm = {:}\n".format(x_normalized)) # normalized solution
# 2π規格化されたAの固有値。
lam01 = 1/4
lam10 = 1/2
print("scaled eigen values = {:} and {:}".format(lam01,lam10))
A with scaling =
[[ 2.35619449 -0.78539816]
[-0.78539816 2.35619449]]
Normalized solution, x_norm = [0.9486833 0.31622777]
scaled eigen values = 0.25 and 0.5
ということです。
qiskitのチュートリアルでは $u_{1}$と$u_{2}$がひっくり返っていますが、間違いです。
QPE
nqubit_ancilla_phase_kickback = 2
nqubit_eigen_state_input = 1
nqubit_ancilla_control = 1 # should be 1
nqubit_total = nqubit_eigen_state_input + nqubit_ancilla_phase_kickback + nqubit_ancilla_control
Ugate = UnitaryGate(U)
Ugate_ctrl = add_control(Ugate, num_ctrl_qubits=1, label='CU', ctrl_state='1')
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)
def qpe(circ, n):
for qubit in range(n):
circ.h(qubit)
repetitions = 1
for counting_qubit in range(n):
for i in range(repetitions):
circ.append(Ugate_ctrl, [counting_qubit, n]) # This is Controlled-U
repetitions *= 2
circ.barrier()
# Apply inverse QFT
qft_dagger(circ, n)
circ.barrier()
CUゲートは無理やり作ってます。
固有値逆回転
C = 1/8
cos_eff = np.sqrt(1-(C**2/lam01**2))
sin_eff = 1*C/lam01
theta_01 = np.arctan(sin_eff/cos_eff)
Ugate_01 = RYGate(2*theta_01)
Ugate_01_ctrl = add_control(Ugate_01, num_ctrl_qubits=2, label='CU_01', ctrl_state='01')
cos_eff = np.sqrt(1-(C**2/lam10**2))
sin_eff = 1*C/lam10
theta_10 = np.arctan(sin_eff/cos_eff)
Ugate_10 = RYGate(2*theta_10)
Ugate_10_ctrl = add_control(Ugate_10, num_ctrl_qubits=2, label='CU_10', ctrl_state='10')
def eigen_rot(circ):
circ.append(Ugate_01_ctrl,[0,1,3])
circ.append(Ugate_10_ctrl,[0,1,3])
QPEで読みだした固有値に応じて、補助量子ビットにその逆数を掛けます
これは位相回転ゲートで実現できます。
今回は固有値がわかっているとして、固有値のビットを制御ビットとみなしてCCUで実現しています。
無理やりです。一般にどうやったらこの操作が（固有値が未知で）簡単にできるのかは難しいらしいです。
なお、固有位相回転をユニタリ行列として実現するために、スケーリングパラメータ$C=1/8$をかけています。
qiskitのチュートリアルでは $ C = 3/8$ と書いてありますが、たぶん間違いです。
$C$は$A$の最小固有値の$2\pi$規格化されたもの(つまり1/4)よりも小さくなければならず、これより大きく取ってしまうとユニタリ行列にならなくなります。
Dojoの記事 7-2. Harrow-Hassidim-Lloyd (HHL) アルゴリズム にもこの旨記載があります。
$C$を小さくしすぎると、後で射影測定で解の状態を選び出すときに、その確率が下がります。
逆QPE
次にQPEの逆をやります。これが結構難しそうだったので、QPEを表現する行列を一度作ってしまって、
（ユニタリなので）逆行列とは共役転置であることを利用して
逆QPEの行列を出します。そしてゲートにしています。
# Define QPE and QPE_inv module
circ_qpe = QuantumCircuit(nqubit_total) # 1 bit for eigen state , 2 bits for phase kick back and 1 bit for controll
qpe(circ_qpe, nqubit_ancilla_phase_kickback)
simulator_unitary = Aer.get_backend('unitary_simulator')
result = execute(circ_qpe, simulator_unitary).result()
QPE_mat = result.get_unitary(circ_qpe)
from qiskit.quantum_info.operators import Operator
Gate_QPE = Operator(QPE_mat)
Gate_QPE_inv = Operator(QPE_mat.conj().T)
HHL
HHL本体を実行します。
HHL = QuantumCircuit(nqubit_total, nqubit_total) # 1 bit for eigen state , 2 bits for phase kick back and 1 bit for controll
HHL.initialize(b,nqubit_total-2)
HHL.unitary(Gate_QPE, range(4), label='QPE')
eigen_rot(HHL)
HHL.unitary(Gate_QPE_inv, range(4), label='QPE_inv')
HHL.draw('mpl', scale=0.5, fold=80)
回路は、
上の2量子ビットq0,q1がQPEの位相書き出し用の補助量子ビット(2桁なので2qubit)です。
その下q2がQPEの入力ベクトルです。
一番下q3が固有値逆回転で回転を食らうビットです。
実行していきます。
backend = Aer.get_backend('statevector_simulator')
results = execute(HHL, backend, shots=1024).result()
statevec = results.get_statevector()
print(statevec)
わけわからないですが、実はここから解が掘り出せます。
そのために、補助量子ビットがq3=1かつq0=q1=1 である状態を射影測定で求めます。
良い機能が見つからなかったので、バカみたいなコードを書きました。
def Sum_statevec_over_ancilla_bit(statevec):
Sum_statevec = []
for i in range(2**4):
bit_str = format(i,'04b')
if bit_str[0]=='1':
if bit_str[2]=='0':
if bit_str[3]=='0':
print(str(bit_str) + '->' + str(statevec[i]))
if bit_str[1]=='0':
Sum_statevec.append(statevec[i]*np.array([1,0]))
else:
Sum_statevec.append(statevec[i]*np.array([0,1]))
return Sum_statevec
Statevec_HHL_solution = sum(Sum_statevec_over_ancilla_bit(statevec))
print(Statevec_HHL_solution)
補助量子ビットがq3=1かつq0=q1=1 である状態 は2種類あります。1000と1100です。
なので、射影測定の結果としてはこれの和になります。それが$[0.375, 0.125]^{T}$だそうです。
しかしこれはノルム1ではない状態です。実際には射影測定後の状態は、これをノルム1となるように規格化したものです。
規格化は、
norm = np.linalg.norm(Statevec_HHL_solution)
Statevec_HHL_solution = Statevec_HHL_solution / norm
この状態ベクトルを古典の（規格化された）解$x$と比較してみると、
print("Exact = \n{:}\n".format(x_normalized))
print("HHL alg. = \n{:}\n".format(Statevec_HHL_solution))
ちゃんと一致しています。
規格化を元に戻す方法もあるみたいですが、勉強中です。
（またしてもqiskitのチュートリアルが間違ってるように見えるので、元の論文に立ち返っています）
結論
HHLが動かせた。固有制御ゲートやスケーリングがめんどくさい。