blueqatで量子フーリエ変換

## この記事に関して
ここではpythonライブラリのblueqatを用いて量子プログラミングを行っていこうと思います。
内容は[公式のチュートリアル](https://github.com/Blueqat/Blueqat-tutorials/tree/master/tutorial-ja)を元に作成しています。
今回は量子フーリエ変換(QFT)を実装していきます。

## 離散フーリエ変換
関数を三角関数の和で書き直すためにフーリエ変換というものを使います。
具体的に言うと三角関数の和で書き直したとき、各三角関数の角振動数とその係数を取り出すということをします。

例えば sin2x なら 2, 1 を取り出すという操作です。
sinx + sin3x なら 1, 1 と 3, 1 を取り出す操作です。

このフーリエ変換は全ての点を考えるのに対し有限個の点のみで考えるフーリエ変換を離散フーリエ変換と言います。

例えば $5\sin 2x$ を考えます。これは離散フーリエ変換をすると 5, 2 が取り出せることがわかります。具体的に計算してみます。
このグラフは以下のようになります。

![img.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/517405/ecd93d6d-11d6-0418-2a35-2100383acd06.png)

これを16分割で離散フーリエ変換します。

![img2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/517405/c9fa7740-e50a-9661-979d-d9819b7776dc.png)

各点の y 座標を $x_j$ と置くと

``` math
x_j = 5\sin2\left(\frac{2\pi j}{16}\right)
```
となります。

これを離散フーリエ変換すると

``` math
y_k = \frac{1}{4}\sum_{j=0}^{15}e^{\frac{2\pi ikj}{16}}x_j = \frac{5}{4}\sum_{j=0}^{15}e^{\frac{2\pi ikj}{16}}\sin2\left(\frac{2\pi j}{16}\right)
```
となるので $\left|y_k\right|/2$ のグラフをかくと

![img3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/517405/803b4786-3773-5146-b1db-32a017ea43f2.png)

上のグラフは振動数 2、 振幅 5 なのできちんと表せることがわかります。
ちなみにコードは以下の通り

``` python:sample.ipynb
import numpy as np
import matplotlib.pyplot as plt

j = 1j
x = np.arange(0, 8, 1.0)
y = []
for k in range(8):
    s = 0
    for i in range(16):
        s += np.exp(2*np.pi*i*j*k/16)*np.sin(4*np.pi*i/16)
    else:
         y.append((np.abs(s))*5/4/2)

plt.plot(x, y)
plt.xlabel("frequency")
plt.ylabel("abs")
plt.show()
```

## 量子フーリエ変換
次に離散フーリエ変換から量子フーリエ変換を考えます。
上の各 $x_j$ を各ビットに埋め込むことを考えます。

すなわち入力する状態を以下のようにします。

``` math
\left|x\right> = \sum_{j=0}^{2^n-1}x_j\left|j\right>
```
この $x_j$ を離散フーリエ変換すると

``` math
y_k = \frac{1}{\sqrt{2^n}}\sum_{j=0}^{2^n-1}e^{\frac{2\pi ijk}{2^n}}x_j
```
となる。また、$\left|y\right>$ について

``` math
\left|y\right> = \sum_{k=0}^{2^n-1}y_k\left|k\right>
```
と置くと、

``` math
\left|y\right> = \frac{1}{\sqrt{2^n}}\sum_{j=0}^{2^n-1}\left(\sum_{k=0}^{2^n-1}e^{\frac{2\pi ijk}{2^n}}\left|k\right>\right)x_j
```
括弧の部分に注目して以下の変換

``` math
QFT: \left|j\right> → \frac{1}{\sqrt{2^n}}\sum_{k=0}^{2^n-1}e^{\frac{2\pi ijk}{2^n}}\left|k\right>
```
のことを量子フーリエ変換と言います。

### 回路の作成
上の式をゲートで実装するために以下のように式を変換します。

``` math
\begin{align}
\frac{1}{\sqrt{2^n}}\sum_{k=0}^{2^n-1}e^{\frac{2\pi ijk}{2^n}}\left|k\right> 
&=\frac{1}{\sqrt{2^n}}\sum_{k_1,\cdots,k_n=0}^{1}
e^{\frac{2\pi ij}{2^n}(2^{n-1}k_1+2^{n-2}k_2+\cdots+2^0k_n)}\left|k_1k_2\cdots k_{n-1}k_n\right> \\
&=\frac{1}{\sqrt{2^n}}\sum_{k_1,\cdots,k_n=0}^{1}
e^{2\pi ij(\frac{k_1}{2}+\frac{k_2}{2^2}+\cdots+\frac{k_n}{2^n})}\left|k_1k_2\cdots k_{n-1}k_n\right> \\
&=\frac{1}{\sqrt{2^n}}\sum_{k_1,\cdots,k_n=0}^{1}
e^{2\pi ij\frac{k_1}{2^n}}\left|k_1\right>\otimes e^{2\pi ij\frac{k_2}{2^{n-1}}}\left|k_2\right>\otimes \cdots\otimes e^{2\pi ij\frac{k_n}{2}}\left|k_n\right> \\
&=\frac{1}{\sqrt{2^n}}
(\left|0\right> + e^{2\pi i\frac{j}{2^n}}\left|1\right>)\otimes (\left|0\right> + e^{2\pi i\frac{j}{2^{n-1}}}\left|1\right>)\otimes \cdots\otimes (\left|0\right> + e^{2\pi i\frac{j}{2}}\left|1\right>) \\
&=\frac{1}{\sqrt2}\left(\left|0\right> + e^{2\pi i\frac{j}{2^n}}\left|1\right>\right)\otimes \frac{1}{\sqrt2}\left(\left|0\right> + e^{2\pi i\frac{j}{2^{n-1}}}\left|1\right>\right)\otimes \cdots\otimes \frac{1}{\sqrt2}\left(\left|0\right> + e^{2\pi i\frac{j}{2}}\left|1\right>\right) \\
&=\frac{1}{\sqrt2}\left(\left|0\right> + e^{2\pi i0.i_1i_2\cdots i_n}\left|1\right>\right)\otimes \frac{1}{\sqrt2}\left(\left|0\right> + e^{2\pi i0.i_{n-1}i_n}\left|1\right>\right)\otimes \cdots\otimes \frac{1}{\sqrt2}\left(\left|0\right> + e^{2\pi i0.i_n}\left|1\right>\right) 
\end{align} \\
(j = 2^{n-1}i_1 + 2^{n-2} + \cdots +2^0i_n)
```
回路は以下のようになります。

``` math
\begin{array}{cccc}
\left|i_1\right> -&H&R_2&\cdots&R_{n-1}&R_n&-&-&-&-&-&\cdots&-&-&-&-&-\\
\left|i_2\right> -&-&*&-&|&|&H&\cdots&R_{n-2}&R_{n-1}&-&\cdots&-&-&-&-&-\\
\vdots&&&&\vdots&\vdots&&&\vdots&\vdots \\
\left|i_{n-1}\right> -&-&-&-&*&|&-&-&*&|&-&\cdots&-&H&R_2&-&-\\
\left|i_n\right> -&-&-&-&-&*&-&-&-&*&-&\cdots&-&-&*&H&-
\end{array}
```
$\left|j\right>$ を量子フーリエ変換する場合は j を 2進数表示した $j = i_1\cdots i_n$ を入力とします。
また $R_n$ は以下の位相ゲートとします。

``` math
R_n = 
\left(\begin{array}{cc}
1 & 0 \\
0 & e^{\frac{2\pi i}{2^n}}
\end{array}\right)
```
この回路から量子フーリエ変換を実装します。

今回は $\left|00\right>, \left|000\right>$ の量子フーリエ変換をします。

#### 2量子ビット
式は以下のようになります。

``` math
\left|00\right> \xrightarrow{QFT} \frac{1}{2}\sum_{k=0}^{3}e^{\frac{2\pi i0k}{4}}\left|k\right> = \frac{1}{2}\sum_{k=0}^{3}\left|k\right>
```
回路は以下のようにできます。

``` math
\begin{array}{ccc}
\left|0\right> &-&H&-&R_2&-&-&-\\
&&&&|& \\
\left|0\right> &-&-&-&*&-&H&-
\end{array}
```
実装します。

``` python:sample.ipynb
from blueqat import Circuit
import math

c = Circuit().h[0].cu1(math.pi/2)[1,0].h[1]
c.run()

結果
array([0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])
```
上の式の通りになりました。

#### 3量子ビット
式は以下の通りになります。

``` math
\left|000\right> \xrightarrow{QFT} \frac{1}{\sqrt8}\sum_{k=0}^{7}e^{\frac{2\pi i0k}{8}}\left|k\right> = \frac{1}{\sqrt8}\sum_{k=0}^{7}\left|k\right>
```
回路は以下のようになります。

``` math
\begin{array}{ccc}
\left|0\right> &-&H&R_2&R_3&-&-&-&-\\
&&&|&|& \\
\left|0\right> &-&-&*&|&H&R_2&-&- \\
&&&&|&&| \\
\left|0\right> &-&-&-&*&-&*&H&-
\end{array}
```
上の回路を実装します。

``` python:sample.ipynb
from blueqat import Circuit
import math

c = Circuit().h[0].cu1(math.pi/2)[1,0].cu1(math.pi/4)[2,0] # 0番目
c.h[1].cu1(math.pi/2)[2,1] # 1番目
c.h[2] # 2番目
c.run()

結果
array([0.35355339+0.00000000e+00j, 0.35355339+1.38777878e-17j,
       0.35355339+0.00000000e+00j, 0.35355339+0.00000000e+00j,
       0.35355339+0.00000000e+00j, 0.35355339+1.38777878e-17j,
       0.35355339+0.00000000e+00j, 0.35355339+0.00000000e+00j])
```
上の式の通りになりました。

## 応用
最後に始めに例で示した離散フーリエ変換を量子フーリエ変換で求めてみようと思います。

とりあえずQFTを以下のように作りました。

``` python:qft.ipynb
from blueqat import Circuit
import numpy as np
import matplotlib.pyplot as plt
import math

def qft(ini):
    """ qftalgorithm"""
    l = len(ini)
    lb = format(l-1, 'b')

    # QFT circuit
    c_qft = Circuit()
    for i in range(len(lb)):
        c_qft.h[i]
        for j in range(1,len(lb)-i):
            c_qft.cu1(math.pi/(2**j))[i+j,i]

    # execute qft
    y_q = []
    for k in range(l):
        c = Circuit()
        b = format(k, 'b')
        bl = len(b)
        for a in range(bl):
            if b[bl-a-1] == '1':
                c = c.x[len(lb)-1-a]
        c = c+c_qft
        y_q.append(c.run())

    # change to y_k
    y = []
    for k in range(l):
        y_k = 0
        for j in range(l):
            y_k += y_q[j][k] * ini[j]
        y.append(y_k)

    return y
```
始めの ini は $x_j$ が入った配列でこれをQFTで $y_k$ の配列に変換しています。
これを使って上の関数を量子フーリエ変換します。

``` python:sample.ipynb
x = np.arange(0, 8, 1.0)
x_j = []
for j in range(16):
    x_j.append(5*np.sin(4*np.pi*j/16)) # x_jを格納

y_k = qft(x_j) # QFT実行

y = []
for j in range(int(len(y_k)/2)):
    y.append(np.abs(y_k[j])/2) # 絶対値を取得

plt.plot(x, y)
plt.xlabel("frequency")
plt.ylabel("abs")
plt.show()
```
この $\left|y_k\right|/2$ は以下のようになります。

![img4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/517405/a2010fde-ab4e-4b06-ae9b-a0f688f95aaf.png)


上手くできたことがわかります。

## まとめ
今回はblueqatを用いて量子フーリエ変換(QFT)を実装しました。
次回はグローバーのアルゴリズムをblueqatで実装してみます。

## 参考文献
宮野 健次郎・古澤 明 (2016) 『量子コンピュータ入門[第２版]』 日本評論社