LoginSignup
3
4

More than 3 years have passed since last update.

@notori48(Kuma)

QiskitでHHLアルゴリズム(線形方程式を解く)

Last updated at Posted at 2021-01-11

HHLアルゴリズム

HHLアルゴリズムは、量子アルゴリズムであって線形方程式 $Ax=b$を解くものです。
先人の解説は、

あたりが良いと思います。現論文は、

(参考) 発展形として変分法による解法もあります。
量子コンピューターで連立一次方程式を変分法で解く量子アルゴリズム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$

image.png

に対して、$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)

回路は、
image.png

上の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)

image.png

わけわからないですが、実はここから解が掘り出せます。
そのために、補助量子ビットが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)

image.png

補助量子ビットが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))

image.png

ちゃんと一致しています。
規格化を元に戻す方法もあるみたいですが、勉強中です。
（またしてもqiskitのチュートリアルが間違ってるように見えるので、元の論文に立ち返っています）

結論

HHLが動かせた。固有制御ゲートやスケーリングがめんどくさい。

3
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
4