$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

はじめに

自作の量子計算シミュレータqlazyで、量子アルゴリズムの基本を確認するよ。今回は「量子フーリエ変換」です。解説記事はここに書ききれないくらいたくさんありましたが、主に参考にさせていただいたものを、以下に挙げておきます。

量子フーリエ変換の導出

量子フーリエ変換というのは、離散フーリエ変換を量子回路で実現したものです。これからその導出を行います。と言っても、ほぼ上の参考記事の受け売りになってしまいますが、、(汗)。

まず、離散フーリエ変換ですが、以下のように、$N$個のデータ系列{$x_j$}から{$y_k$}への変換式として定義されます。

y_{k} = \frac{1}{\sqrt{N}} \sum_{j-0}^{N-1} exp(i\frac{2\pi kj}{N}) x_{j}  \tag{1}

これを量子計算として実行するために、データ系列を量子状態のN個の基底に対する係数(確率振幅)と見立てるようにします。つまり、

\begin{align}
\ket{\boldsymbol{x}} &= \sum_{j=0}^{N-1} x_{j} \ket{j} \\
\ket{\boldsymbol{y}} &= \sum_{j=0}^{N-1} y_{j} \ket{j}
\end{align}

です。この2番目の式に基づき、(1)式の$y_{k}$を係数とした重ね合わせ状態をつくると、以下のようになります。

\begin{align}
\ket{\boldsymbol{y}} &= \sum_{k=0}^{N-1} y_{k} \ket{k}
= \sum_{k=0}^{N-1} \sum_{j=0}^{N-1} W_{kj} x_{j} \ket{k} \\
&= \sum_{j=0}^{N-1} x_{j} \bigl[ \sum_{k=0}^{N-1} W_{kj} \ket{k} \bigr]
\end{align}

ここで、

W_{kj} = \frac{1}{\sqrt{N}}exp(i \frac{2\pi kj}{N})

と定義しました。ちょっと式を見比べてみればわかりますが、この変換は、

\ket{j} \rightarrow \sum_{k=0}^{N-1} W_{kj} \ket{k}
= \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} exp(i \frac{2\pi kj}{N}) \ket{k}

という基底変換を行ったことと等価です。なので、今後、この基底変換に注目していきます。

まず、量子状態の次元(基底の数)は常に2のべき乗と決まっているので、いままで$N$と書いていたものを$2^{n}$と書き直して、式を展開します。

\begin{align}
\ket{j}
&\rightarrow \frac{1}{\sqrt{2^{n}}} \sum_{k=0}^{2^{n}-1}exp(i\frac{2\pi kj}{2^{n}})\ket{k}\\
&= \frac{1}{\sqrt{2^{n}}} \sum_{k_{0}=0,1} \sum_{k_{1}=0,1} \cdots \sum_{k_{n-1}=0,1}
exp(i \frac{2\pi j}{2^{n}} \sum_{l=0}^{n-1} k_{l} 2^{l}) \ket{k_{n-1} \cdots k_{0}} \\
&= \frac{1}{\sqrt{2^{n}}} \sum_{k_{0}=0,1} \sum_{k_{1}=0,1} \cdots \sum_{k_{n-1}=0,1}
\prod_{l=0}^{n-1} exp(i2\pi jk_{l}2^{l-n}) \ket{k_{l}} \\
&= \frac{1}{\sqrt{2^{n}}} \prod_{l=0}^{n-1} (\sum_{k_{l}=0,1} exp(i 2\pi jk_{l}2^{l-n}) \ket{k_{l}}) \\
&= \frac{1}{\sqrt{2^{n}}} \prod_{l=0}^{n-1} (\ket{0} + exp(i 2\pi j2^{l-n}) \ket{1}) \\
&= \frac{1}{\sqrt{2^{n}}} \prod_{l=0}^{n-1} (\ket{0} + exp(i 2\pi 0.j_{n-l-1} \cdots j_{0}) \ket{1}) \\
&= \frac{1}{\sqrt{2^{n}}} (\ket{0}+e^{i2\pi 0.j_{0}}\ket{1})(\ket{0}+e^{i2\pi 0.j_{1}j_{0}}\ket{1}) \cdots (\ket{0}+e^{i2\pi 0.j_{n-1}j_{n-2} \cdots j_{0}}\ket{1}) \tag{2}
\end{align}

ここで、$j$を2進数で"$j_{n-1}j_{n-2} \cdots j_{0}$"と書くことにすると、

\begin{align}
j2^{l-n} &= 2^{l} \times 0.j_{n-1}j_{n-2} \cdots j_{0} \\
&= j_{n-1}j_{n-2} \cdots j_{n-l}.j_{n-l-1} \cdots j_{0}
\end{align}

となり、これを$exp$の肩にのせます。そうすると、

exp(i 2\pi j2^{l-n}) = exp(i 2\pi 0.j_{n-l-1} \cdots j_{0})

となります($2\pi$があるので小数部しか残らないことに注意)。上の(2)式を導出する途中で、これを使いました。

さて、ここまで来たところで、量子回路での実現を具体的に考えてみます。いきなり(2)式を相手にするのはしんどいので、まず、簡単なところから考えます。

n=1の場合

(2)式は、

\ket{j_{0}} \rightarrow \frac{1}{\sqrt{2}} (\ket{0} + e^{i2\pi 0.j_{0}} \ket{1})

です。つまり、

\begin{align}
\ket{0} \rightarrow \frac{1}{\sqrt{2}} (\ket{0} + \ket{1}) \\
\ket{1} \rightarrow \frac{1}{\sqrt{2}} (\ket{0} - \ket{1})
\end{align}

ということなので、量子回路はアダマールゲート１つです。

|j_0> ---H---

n=2の場合

(2)式は、

\ket{j_{1}j_{0}} \rightarrow \frac{1}{2} (\ket{0}+e^{i2\pi 0.j_{0}}\ket{1})
(\ket{0}+e^{i2\pi 0.j_{1}j_{0}}\ket{1})

という形になります。右辺の１番目のカッコは、先程と同じで$j_{0}$の量子ビットにアダマールをかけたものです。問題は２番目のカッコです。よぉーく考えればだんだんわかってきますが、$j_{1}$の量子ビットにアダマールをかけた上で、さらに$j_{0}$の量子ビットが1のときのみ$\ket{1}$に$e^{i2\pi/4}$をかけたものです。つまり、２番目のカッコは、$j_{0}$の量子ビットが1のときのみ、$j_{1}$の量子ビットに位相ゲート、

P =
\begin{pmatrix}
1 & 0 \\
0 & e^{i2\pi/4}
\end{pmatrix}

をかけたもの、つまり、制御位相ゲートをかけたものになります。したがって、量子回路は、

|j_1> ---H--P(1)-----
|j_0> -------*----H--

となります。ここでP(k)は、

P_{k} = P(2\pi/2^{k+1}) =
\begin{pmatrix}
1 & 0 \\
0 & e^{i2\pi/2^{k+1}}
\end{pmatrix}

を表し、回路中の"-*-"と合わせて制御位相ゲートを表すものとします。

n=3の場合

式の展開は省略します。同様の議論により、量子回路は、

|j_2> ---H-P(1)-P(2)--------------
|j_1> ------*----------H-P(1)-----
|j_0> -----------*--------*----H--

という形になります。

n=任意の場合

n=1,2,3の場合を順によぉーく見比べると、なんとなくパターンが見えてきます。結局、nが任意の場合は、

|j_n-1> --H-P(1)-...-P(n-2)-P(n-1)---------...---------------...-------------
|j_n-2> -----*---...----------------H-P(1)-...-P(n-3)-P(n-2)-...-------------
 ...             ...                       ...               ...               
 
|j_1>   ---------...--*--------------------...--*-------------...-H-P(1)-----
|j_0>   ---------------------*-------------------------*------...----*----H--

という形になります。

シミュレータで確認

それでは、qlazyで、量子フーリエ変換の動作を確認します。まず、n=2の場合です。コードは以下の通りです。

【2021.9.5追記】qlazy最新版でのソースコードはここに置いてあります。

from qlazypy import QState

def swap(self,q0,q1):

    self.cx(q0,q1).cx(q1,q0).cx(q0,q1)
    return self

def qft2(self,q1,q0):

    self.h(q1).cp(q0,q1,phase=0.5)
    self.h(q0)
    self.swap(q0,q1)
    return self
  
def main():

    QState.swap = swap
    QState.qft2 = qft2
    
    qs = QState(2)
    print("== before ==")
    qs.show()

    qs.qft2(0,1)
    print("== after ==")
    qs.show()

    qs.free()

if __name__ == '__main__':
    main()

何をやっているか、簡単に説明します。

def swap(self,q0,q1):

    self.cx(q0,q1).cx(q1,q0).cx(q0,q1)
    return self

で、swapゲートを定義しています。

def qft2(self,q1,q0):

    self.h(q1).cp(q0,q1,phase=0.5)
    self.h(q0)
    self.swap(q0,q1)
    return self

で、2量子ビットの量子フーリエ変換を定義しています(hはアダマール、cpは制御位相ゲートを表しています)。その上で、main()の中の、

QState.swap = swap
QState.qft2 = qft2

で、この関数swap(),qft2()をQStateクラスのメソッドとして定義しています。

(Pythonでは、ユーザー側で定義したメソッドを既存クラスのメソッドとして、こんな風に簡単に追加することができます。が、既存のメソッドと名前が衝突しないようにしないといけないので要注意です。本当はqlazyの側で、ユーザー定義ゲートを登録する手段を用意しておくべきだと思うのですが、、要検討です)

あとは、初期状態qsを用意して、

qs.qft2(0,1)

で、量子フーリエ変換ゲートをかけているだけです。

実行結果は以下の通りです。

== before ==
c[00] = +1.0000+0.0000*i : 1.0000 |+++++++++++
c[01] = +0.0000+0.0000*i : 0.0000 |
c[10] = +0.0000+0.0000*i : 0.0000 |
c[11] = +0.0000+0.0000*i : 0.0000 |
== after ==
c[00] = +0.5000+0.0000*i : 0.2500 |++++
c[01] = +0.5000+0.0000*i : 0.2500 |++++
c[10] = +0.5000+0.0000*i : 0.2500 |++++
c[11] = +0.5000+0.0000*i : 0.2500 |++++

というわけで、ちゃんとフーリエ変換が実行できていることがわかります。

次に、任意の量子ビットに対する量子フーリエ変換をやってみます。コードは以下です。

from qlazypy import QState

def swap(self,q0,q1):

    self.cx(q0,q1).cx(q1,q0).cx(q0,q1)
    return self

def qft(self,id=None):

    dim = len(id)
    iid = id[::-1]

    for i in range(dim):
        self.h(iid[dim-i-1])
        phase = 1.0
        for j in range(dim-i-1):
            phase /= 2.0
            self.cp(iid[dim-i-j-2],iid[dim-i-1],phase=phase)

    i = 0
    while i < dim-1-i:
        self.swap(iid[i], iid[dim-1-i])
        i += 1

    return self

def main():

    QState.swap = swap
    QState.qft = qft
    
    qs = QState(3)
    print("== before ==")
    qs.show()

    qs.qft(id=[0,1,2])
    print("== after ==")
    qs.show()

    qs.free()

if __name__ == '__main__':
    main()

関数qft()の引数として、量子ビット番号のリストidを指定するようにしています。中身は、先程の量子回路を素直に実装しただけです。これを先程と同じように、QStateクラスのメソッドとして追加して、3量子ビットの状態$\ket{000}$に対して、量子フーリエ変換するようにしました。

実行結果は以下の通りです。

== before ==
c[000] = +1.0000+0.0000*i : 1.0000 |+++++++++++
c[001] = +0.0000+0.0000*i : 0.0000 |
c[010] = +0.0000+0.0000*i : 0.0000 |
c[011] = +0.0000+0.0000*i : 0.0000 |
c[100] = +0.0000+0.0000*i : 0.0000 |
c[101] = +0.0000+0.0000*i : 0.0000 |
c[110] = +0.0000+0.0000*i : 0.0000 |
c[111] = +0.0000+0.0000*i : 0.0000 |
== after ==
c[000] = +0.3536+0.0000*i : 0.1250 |++
c[001] = +0.3536+0.0000*i : 0.1250 |++
c[010] = +0.3536+0.0000*i : 0.1250 |++
c[011] = +0.3536+0.0000*i : 0.1250 |++
c[100] = +0.3536+0.0000*i : 0.1250 |++
c[101] = +0.3536+0.0000*i : 0.1250 |++
c[110] = +0.3536+0.0000*i : 0.1250 |++
c[111] = +0.3536+0.0000*i : 0.1250 |++

というわけで、こちらもちゃんとフーリエ変換が実行できました。

おわりに

前回の「アダマールテスト」に続いて、今回「量子フーリエ変換」を確認しました。両方とも、これ単独で何か役に立つアルゴリズムというよりも、何かのサブルーチンとして使われるもののようです。特に、量子フーリエ変換は、直接観測することができない量子状態の係数(確率振幅)を変換するアルゴリズムです。量子フーリエ変換によって特定の係数のみピークを持つような形にしてから観測することで、何らかの意味のある結果を見出すような、そんな使い方がなされているようです。ショアのアルゴリズムやグローバーのアルゴリズムが、その代表例です。それら代表選手を確認していく前に、次回は、もう一つのサブルーチン的な基本アルゴリズムとして重要な「位相推定アルゴリズム」について、見ていこうと思います。

訂正(2019.6.11)

QFTのPythonコード(量子ビット数が2の場合と任意の場合の両方)が間違っていたので、訂正しました。引数として指定するビットの並びを逆順にして、さらにアプトプットの量子状態のビットの並びも逆転させないといけませんでした。そのため、swapゲートも追加で定義して、QFT関数の中で使うようにしました。QFT実装のときの「ハマりどころ」にきれいにハマってしまいました(汗)。皆様も１から実装する際にはご注意くださいまし。

以上