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


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

- [QiitaのHHL アルゴリズム](https://qiita.com/KeiichiroHiga/items/ee165892c9851913afca)
- [qiskitのチュートリアル Solving Linear Systems of Equations using HHL](https://qiskit.org/textbook/ch-applications/hhl_tutorial.html) 
- [Dojoの記事 7-2. Harrow-Hassidim-Lloyd (HHL) アルゴリズム](https://dojo.qulacs.org/ja/latest/notebooks/7.2_Harrow-Hassidim-Lloyd_algorithm.html)

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

- [Quantum algorithm for solving linear systems of equations](https://arxiv.org/abs/0811.3171)

(参考) 発展形として変分法による解法もあります。
[量子コンピューターで連立一次方程式を変分法で解く量子アルゴリズムVQLSを調査した件](https://qiita.com/Ayumu_walker/items/6b1131c6be64a8c5ad37)
こちらは、動画を観るとよくわかります。[動画version](https://www.youtube.com/watch?v=j3gYU4-7VFc&t=59s)

今回はqiskitのチュートリアルのほうを見ていきます。
ただ、後述するように、このチュートリアルは具体的な計算例に間違いがかなりあります。

# qiskit でHHLアルゴリズム

## import

```python
#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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/7664bd35-c849-9219-9f2b-95f33979c216.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$が規格化されたものが量子状態の振幅に埋め込まれて出てきます。
規格化係数も後で測定で求められるようですが、完全に理解できていないので今回は扱いません。

コードで書くと

```python
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

```python
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ゲートは無理やり作ってます。

### 固有値逆回転

```python
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) アルゴリズム](https://dojo.qulacs.org/ja/latest/notebooks/7.2_Harrow-Hassidim-Lloyd_algorithm.html) にもこの旨記載があります。
$C$を小さくしすぎると、後で射影測定で解の状態を選び出すときに、その確率が下がります。

### 逆QPE

次にQPEの逆をやります。これが結構難しそうだったので、QPEを表現する行列を一度作ってしまって、
（ユニタリなので）逆行列とは共役転置であることを利用して
逆QPEの行列を出します。そしてゲートにしています。

```python
# 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本体を実行します。

```python

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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/51eddfd8-4b62-3e69-61e7-f8515e724659.png)

上の2量子ビットq0,q1がQPEの位相書き出し用の補助量子ビット(2桁なので2qubit)です。
その下q2がQPEの入力ベクトルです。
一番下q3が固有値逆回転で回転を食らうビットです。

実行していきます。

```python
backend = Aer.get_backend('statevector_simulator')
results = execute(HHL, backend, shots=1024).result()
statevec = results.get_statevector()
print(statevec)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/2837d52b-3d8e-5b37-956d-72a9a6573e77.png)

わけわからないですが、実はここから解が掘り出せます。
そのために、補助量子ビットがq3=1かつq0=q1=1 である状態を射影測定で求めます。
良い機能が見つからなかったので、バカみたいなコードを書きました。

```python
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
```

```python
Statevec_HHL_solution = sum(Sum_statevec_over_ancilla_bit(statevec))
print(Statevec_HHL_solution)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/c03f854c-7f47-f44c-a469-93c14c582cad.png)


補助量子ビットがq3=1かつq0=q1=1 である状態 は2種類あります。1000と1100です。
なので、射影測定の結果としてはこれの和になります。それが$[0.375, 0.125]^{T}$だそうです。

しかしこれはノルム1ではない状態です。実際には射影測定後の状態は、これをノルム1となるように規格化したものです。
規格化は、

```python
norm = np.linalg.norm(Statevec_HHL_solution)
Statevec_HHL_solution = Statevec_HHL_solution / norm
```

この状態ベクトルを古典の（規格化された）解$x$と比較してみると、

```python
print("Exact = \n{:}\n".format(x_normalized))
print("HHL alg. = \n{:}\n".format(Statevec_HHL_solution))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/df33fe4f-5cf2-2aa3-933f-86d162615ff3.png)

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

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