量子フーリエ変換

はじめに

 先の回では「量子もつれ」という現象を紹介し、blueqatという名前のフレームワークを用いて量子もつれをシミュレートした。今回は、離散フーリエ変換を量子コンピュータ(量子回路)で行うアルゴリズム「量子フーリエ変換」を紹介し、先と同じようにblueqatでシミュレートする。

離散フーリエ変換

 離散値x_j\;(j=0,\cdots,N-1)の離散フーリエ変換は次式で定義される。

(1)    \begin{equation*} y_k=\frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}x_j\exp{\left(i\frac{2\pi kj}{N}\right)} \end{equation*}

ここで、k=0,\cdots,N-1であり、iは虚数単位である。いま

(2)    \begin{equation*} w_{kj}=\frac{1}{\sqrt{N}}\exp{\left(i\frac{2\pi kj}{N}\right)} \end{equation*}

を定義すると、式(1)は

(3)    \begin{equation*} y_k=\sum_{j=0}^{N-1}w_{kj}x_j \end{equation*}

と書くことができる。w_{kj}を成分に持つ行列Wを導入すると、式(3)は

(4)    \begin{equation*} y=Wx \end{equation*}

と書ける。ここで、y\in\mathbb{R}^N,x\in\mathbb{R}^N,W\in\mathbb{R}^{N\times N}である。行列Wはユニタリー行列になっていることに注意する。すなわち

(5)    \begin{equation*} W^{\dagger}W=WW^{\dagger}=I \end{equation*}

が成り立つ。ここで、W^{\dagger}Wのエルミート共役(転置して複素共役をとる)を表す。Iは単位行列である。ユニタリー性を満たさない演算は量子回路上に載せることができないので、離散フーリエ変換のユニタリー性(式(5))は重要である。

振幅エンコーディング

 xyのようなベクトルをどのように量子コンピュータに載せるのかを考える。これまでのブログ(1,2)で1量子ビットと2量子ビットを扱った。任意の量子状態は、1量子ビットの場合は|0\rangle,|1\rangleの2つの基底状態の重ね合わせで、2量子ビットの場合は|00\rangle,|01\rangle,|10\rangle,|11\rangleの4つの基底状態の重ね合わせで表現されることを見た。ここで、0や1の並びは0以上の整数の2進数表記に相当し、2量子ビットの場合の00,01,10,11は、それぞれ、10進数表記で0,1,2,3に相当する。いま、4次元の単位ベクトルx=(x_0,x_1,x_2,x_3)が与えられたとき、各成分を係数に持つ次の量子状態を考える。

(6)    \begin{eqnarray*} |x\rangle &\equiv& \sum_{j=0}^{3}x_j|j\rangle \\ &=&x_0|0\rangle+x_1|1\rangle+x_2|2\rangle+x_3|3\rangle \\ &=&x_0|00\rangle+x_1|01\rangle+x_2|10\rangle+x_3|11\rangle \end{eqnarray*}

xは単位ベクトルとしたから量子状態が満たすべき条件(規格化条件:\sum_{j=0}^{3}|x_j|^2=1)は満たされている。つまり、4次元ベクトルは2量子ビットを用いて上のように表すことができる。これをnビットに拡張する。

(7)    \begin{equation*} |x\rangle =\sum_{j=0}^{N-1}x_j|j\rangle \end{equation*}

ただし、N=2^nである。式(7)はN次元ベクトルを表す量子状態である。n=3なら

(8)    \begin{eqnarray*} |x\rangle &=& \sum_{j=0}^{2^3-1}x_j|j\rangle \\ &=& x_0|0\rangle+x_1|1\rangle+\cdots+x_7|7\rangle \\ &=& x_0|000\rangle+x_0|001\rangle+\cdots+x_7|111\rangle \\ \end{eqnarray*}

となり、8(=2^3)次元ベクトルを表すことができる。このようにベクトルの成分を量子状態の係数(確率振幅)に埋め込むことを振幅エンコーディングと呼ぶ。どのようにしてこの状態を作り出すのかという話は少々複雑なので今回は言及しない。

量子フーリエ変換

 振幅エンコーディングにより、ベクトルyを量子状態に変換する。

(9)    \begin{equation*} |y\rangle =\sum_{k=0}^{N-1}y_k|k\rangle \end{equation*}

ここに式(3)を代入すると

(10)    \begin{eqnarray*} |y\rangle &=& \sum_{k=0}^{N-1}y_k|k\rangle \\ &=& \sum_{k=0}^{N-1}\left[\sum_{j=0}^{N-1}w_{kj}x_j\right]|k\rangle \\ &=& \sum_{j=0}^{N-1} x_j \left[\sum_{k=0}^{N-1}w_{kj}|k\rangle\right] \end{eqnarray*}

を得る。ところでxの量子状態は

(11)    \begin{equation*} |x\rangle =\sum_{j=0}^{N-1}x_j|j\rangle \end{equation*}

であった。この2式を見比べると、|x\rangleから|y\rangleに変換するには次の変換を行えばよいことが分かる。

(12)    \begin{equation*} |j\rangle \rightarrow \sum_{k=0}^{N-1}w_{kj}|k\rangle \end{equation*}

すなわち、この変換をF_nと書くことにすると

(13)    \begin{equation*} F_n|j\rangle =\sum_{k=0}^{N-1}w_{kj}|k\rangle \end{equation*}

であるから

(14)    \begin{eqnarray*} F_n|x\rangle &=& F_n\sum_{j=0}^{N-1}x_j|j\rangle \\ &=& \sum_{j=0}^{N-1}x_j\left(F_n|j\rangle\right) \\ &=& \sum_{j=0}^{N-1}x_j\left[\sum_{k=0}^{N-1}w_{kj}|k\rangle\right] \\ &=& \sum_{k=0}^{N-1}\left[\sum_{j=0}^{N-1}w_{kj}x_j\right]|k\rangle \\ &=& \sum_{k=0}^{N-1}y_k|k\rangle \\ &=& |y\rangle \end{eqnarray*}

となる。変換F_nを量子フーリエ変換(Quantum Fourier Transform:QFT)と呼ぶ。上式を見ると、F_n|x\rangleの結果からフーリエ係数y_kを求めるには|k\rangleの係数を取り出す必要があることに注意しなければならない。これを行う手間についてはまた後で触れることにする。

QFTアルゴリズム(ステップ1)

 QFTを量子回路に載せる手順(QFTアルゴリズム)を説明する。まず最初に、QFTの式を量子回路に載せやすい形に変形する必要がある。ステップ1ではその変形の過程を示す。古典回路と同様に量子回路の場合も2進数の言葉で書き直す必要がある

 式(12)の右辺に式(2)を代入し、N2^nに戻すと

(15)    \begin{equation*} \sum_{k=0}^{2^n-1}\frac{1}{\sqrt{2^n}}\exp{\left(i\frac{2\pi kj}{2^n}\right)} |k\rangle \end{equation*}

となる。ところで、knビット2進数表記がk_{0}k_{1}\cdots k_{n-1}となるとき(例えば、3の2進数表記は11なのでk_0=1,k_1=1)、次式が成り立つ。

(16)    \begin{equation*} k=\sum_{l=0}^{n-1}2^{n-l-1} k_{l} \end{equation*}

k_lは0か1を取る変数である。状態|k\rangleを2進数で書き直すと

(17)    \begin{eqnarray*} |k\rangle &=& |k_{0}k_{1}\cdots k_{n-1}\rangle \\ &=& |k_{0}\rangle|k_{1}\rangle\cdots|k_{n-1}\rangle \end{eqnarray*}

となる。表記を簡単にするため直積記号\otimesは省略した。さらにkについての和は次のように書き直すことができる。

(18)    \begin{equation*} \sum_{k=0}^{2^n-1}\rightarrow\sum_{k_{0}=0}^{1}\;\sum_{k_{1}=0}^{1}\cdots\sum_{k_{n-1}=0}^{1}  \end{equation*}

以上のことを考慮すると式(15)は

(19)    \begin{eqnarray*} &&\sum_{k=0}^{2^n-1}\frac{1}{\sqrt{2^n}}\exp{\left(i\frac{2\pi kj}{2^n}\right)} |k\rangle \\  &=&\frac{1}{\sqrt{2^n}} \sum_{k_{0}=0}^{1}\;\sum_{k_{1}=0}^{1}\cdots\sum_{k_{n-1}=0}^{1}\exp{\left(i\frac{2\pi j\sum_{l=0}^{n-1}2^{n-l-1} k_l}{2^n}\right)} |k_{0}\rangle|k_{1}\rangle\cdots|k_{n-1}\rangle \\ &=&\frac{1}{\sqrt{2^n}} \sum_{k_{0}=0}^{1}\;\sum_{k_{1}=0}^{1}\cdots\sum_{k_{n-1}=0}^{1}\exp{\left(i2\pi j\sum_{l=0}^{n-1}2^{-l-1} k_l\right)} |k_{0}\rangle|k_{1}\rangle\cdots|k_{n-1}\rangle \\ &=&\frac{1}{\sqrt{2^n}} \prod_{l=0}^{n-1}\left[\sum_{k_{l}=0}^{1}\exp{\left(i2\pi j2^{-l-1} k_l\right)} |k_{l}\rangle\right] \\ &=&\frac{1}{\sqrt{2^n}} \prod_{l=0}^{n-1}\left[|0\rangle+\exp{\left(i2\pi j2^{-l-1}\right)} |1\rangle\right] \end{eqnarray*}

となる。ここで、乗積を展開するときはl=0の項が一番左側に来ると約束する。jについても2進数表記を行うと

(20)    \begin{equation*} j=\sum_{l=0}^{n-1}2^{n-l-1} j_l \end{equation*}

従って

(21)    \begin{eqnarray*} j2^{-l-1} &=& 2^{-l-1}\left(2^{n-1}j_{0}+\cdots+2^{l+1}j_{n-l-2}+2^{l}j_{n-l-1}+\cdots+2^0j_{n-1}\right) \\ &=& 2^{n-l-2}j_{0}+\cdots+2^{0}j_{n-l-2}+2^{-1}j_{n-l-1}+\cdots+2^{-l-1}j_{n-1} \end{eqnarray*}

となる。右辺の先頭のn-l-1個の和は正の整数である。これを整数M=2^{n-l-2}j_{0}+\cdots+2^{0}j_{n-l-2}に置き換え、式(21)を指数関数の肩に載せると

(22)    \begin{eqnarray*} \exp{\left(i2\pi j2^{-l-1}\right)} &=& \exp{\left(i2\pi\left[M+2^{-1}j_{n-l-1}+\cdots+2^{-l-1}j_{n-1}\right]\right)} \\ &=& \exp{\left(i2\pi\left[2^{-1}j_{n-l-1}+\cdots+2^{-l-1}j_{n-1}\right]\right)} \end{eqnarray*}

となる。ここで、\exp{\left(i2\pi M\right)}=1を用いた。浮動小数点の2進数表記を用いると

(23)    \begin{equation*} 2^{-1}j_{n-l-1}+\cdots+2^{-l-1}j_{n-1}=0.j_{n-l-1}j_{n-l}\cdots j_{n-1} \end{equation*}

であるから、結局

(24)    \begin{equation*} \exp{\left(i2\pi j2^{-l-1}\right)}= \exp{\left(i2\pi\;0.j_{n-l-1}j_{n-l}\cdots j_{n-1}\right)} \end{equation*}

を得る。以上の結果をまとめると、式(19)は最終的に次のように変形される。

(25)    \begin{equation*} &&\sum_{k=0}^{2^n-1}\frac{1}{\sqrt{2^n}}\exp{\left(i\frac{2\pi kj}{2^n}\right)} |k\rangle= \frac{1}{\sqrt{2^n}} \prod_{l=0}^{n-1}\left[|0\rangle+\exp{\left(i2\pi0.j_{n-l-1}j_{n-l}\cdots j_{n-1} \right)} |1\rangle\right] \end{equation*}

改めて量子フーリエ変換を書くと

(26)    \begin{eqnarray*} F_n|j\rangle &=& \frac{1}{\sqrt{2^n}} \prod_{l=0}^{n-1}\left[|0\rangle+\exp{\left(i2\pi0.j_{n-l-1}j_{n-l}\cdots j_{n-1} \right)}|1\rangle\right] \\ &=& \frac{1}{\sqrt{2^n}}\left(|0\rangle+e^{i2\pi0.j_{n-1}}|1\rangle\right)\left(|0\rangle+e^{i2\pi0.j_{n-2}j_{n-1}}|1\rangle\right)\cdots\left(|0\rangle+e^{i2\pi 0.j_{0}j_{1}\cdots j_{n-1}}|1\rangle\right) \end{eqnarray*}

となる。乗積を展開するときはl=0の項が一番左側に来るように配置する。あとはこの式を量子回路で表現すれば良い。その手順を次に示す。

QFTアルゴリズム(ステップ2)

 いきなりn量子ビットの場合を考えると難しいので、最初に1量子ビット、2量子ビットの場合を考える。

QFTアルゴリズム:1量子ビットの場合

 n=1のときQFTは次式になる。

(27)    \begin{equation*} F_1|j\rangle &=& \frac{1}{\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_0}|1\rangle\right) \end{equation*}

jの2進数表記はj_0なので

(28)    \begin{equation*} F_1|j_0\rangle = \frac{1}{\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_0}|1\rangle\right) \end{equation*}

である。j_0=0,1を代入すると

(29)    \begin{eqnarray*} F_1|0\rangle &=& \frac{1}{\sqrt{2}}\left(|0\rangle+|1\rangle\right) \\ F_1|1\rangle &=& \frac{1}{\sqrt{2}}\left(|0\rangle-|1\rangle\right) \\ \end{eqnarray*}

を得る。ここで、e^{i2\pi0.1}=e^{i2\pi/2}=-1を用いた。2進数0.1は10進数で表すと1/2であることに注意する。式(29)は先のブログで紹介したアダマールゲートHで実現できる変換である。すなわち、F_1=Hである。n=1の場合のQFTを実現する量子回路は下図である。

アダマールゲートHの行列表記は

(30)    \begin{equation*} H=\frac{1}{\sqrt{2}}\left(     \begin{array}{cc}       1 & 1 \\       1 & -1 \\     \end{array}   \right) \end{equation*}

である。H|0\rangle=(1,0)|1\rangle=(0,1)に作用させると式(29)が得られることを確認してほしい。

QFTアルゴリズム:2量子ビットの場合

 天下り的であるが先に2量子ビットのQFTを実現する量子回路を示す。

最初の状態は

(31)    \begin{equation*} |\psi_0\rangle\equiv|j_0j_1\rangle=|j_0\rangle|j_1\rangle \end{equation*}

である。先頭の量子ビット|j_0\rangleにアダマールゲートHをかけると、n=1の場合のQFTの式を参考にして

(32)    \begin{eqnarray*} |\psi_0\rangle &\rightarrow&  \left(H|j_0\rangle\right)|j_1\rangle \\ &=&  \frac{1}{\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_0}|1\rangle\right)|j_1\rangle\equiv|\psi_1\rangle \end{eqnarray*}

を得る。次にj_1=1のときだけ先頭の状態|j_0\rangleに位相ゲートR_1を作用させる(このとき量子ビット|j_1\rangleは制御ビットと呼ばれる)。R_1は次式で定義される行列である。

(33)    \begin{equation*} R_1=\left(     \begin{array}{cc}       1 & 0 \\       0 & e^{i2\pi/4} \\     \end{array}   \right) \end{equation*}

この変換を丁寧に見ていく。j_1=0のときは何もしないので状態は|\psi_1\rangleのままである。j_1=1のときは、量子ビット|j_0\rangleR_1を作用させるので

(34)    \begin{eqnarray*} |\psi_1\rangle &\rightarrow& \left(R_1\frac{1} {\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_0}|1\rangle\right)\right)|j_1=1\rangle \\ &=&\frac{1} {\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_01}|1\rangle\right)|j_1=1\rangle \end{eqnarray*}

となる。j_1=0,1をまとめて書くと

(35)    \begin{equation*} |\psi_1\rangle\rightarrow \frac{1} {\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_0j_1}|1\rangle\right)|j_1\rangle\equiv|\psi_2\rangle \end{equation*}

と変化することになる。この状態の後ろの量子ビット|j_1\rangleにアダマールゲートHをかけると

(36)    \begin{eqnarray*} |\psi_2\rangle &\rightarrow&  \frac{1}{\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_0j_1}|1\rangle\right) \left(H|j_1\rangle\right) \\ &\rightarrow&  \frac{1}{\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_0j_1}|1\rangle\right) \frac{1}{\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_1}|1\rangle\right)\equiv|\psi_3\rangle \end{eqnarray*}

最後に、2つの量子ビットの順番を交換するSWAPゲートSを作用させる(今回はSWAPゲートの詳細は省略する)。

(37)    \begin{equation*} |\psi_3\rangle\rightarrow  \frac{1}{\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_1}|1\rangle\right) \frac{1}{\sqrt{2}}\left(|0\rangle+e^{i2\pi0.j_0j_1}|1\rangle\right)\equiv|\psi_4\rangle \end{equation*}

最終的に得られた状態|\psi_4\rangleは式(26)でn=2としたものである。すなわち

(38)    \begin{equation*} F_2|\psi_0\rangle=|\psi_4\rangle \end{equation*}

である。

QFTアルゴリズム:n量子ビットの場合

 n=2の量子回路を拡張してn=3の量子回路を得る。

さらに、拡張すれば、一般のnの場合の量子回路を得る。

ただし

(39)    \begin{equation*} R_k=\left(     \begin{array}{cc}       1 & 0 \\       0 & e^{i2\pi/2^{k+1}} \\     \end{array}   \right) \end{equation*}

である。n\geq3の場合、SWAPゲートSによりビット列は逆順に並べ替えられる。以上で量子フーリエ変換を量子回路で表現できた。

古典アルゴリズム(FFT)との比較

 ここでは、QFTアルゴリズムと古典アルゴリズム(FFT)の計算量の見積もりを行う。まず最初にQFTアルゴリズムから考える。量子ゲートであるHR_kの個数の和は、回路内の横線上にある四角の数を数えればよいから

(40)    \begin{equation*} \sum_{k=1}^{n}k=\frac{1}{2}n(n+1) \end{equation*}

である。図のSで示したSWAPゲートについては今回はその詳細を説明しないが、3n/2個のCNOTゲートから構成される(CNOTゲートは先のブログで紹介した)。従って、QFTアルゴリズムで使う総ゲート数は

(41)    \begin{equation*} \frac{1}{2}n(n+1)+\frac{3n}{2} \end{equation*}

となる。これはnが大きくなるとn^2のオーダーで大きくなる量\mathcal{O}\left(n^2\right)である。フーリエ係数の個数N=2^nに換算すると\mathcal{O}\left((\log_2{N})^2\right)となる。これが、QFTアルゴリズムの計算量である。一方、FFTの計算量は、\mathcal{O}\left(N\log_2{N}\right)である。下図に見る通り、QFTの方が高速であることが分かる。

式(14)の後の文章で触れたとおり、QFTの場合、計算後にN個のフーリエ係数を取り出す必要がある。この作業は指数関数的な時間を必要とするため、単にフーリエ係数を知りたいだけならFFTの方が速いというオチが付く。

blueqatによる実装

 今回作成したソースコードはここにある。まず最初に2量子ビットの場合のQFTを考える。このときの量子回路は以下であった。

初期状態が|j_0j_1\rangle=|00\rangleの場合のQFTを実現するコードは以下の通り(main_2.py)。

コード内のコメントで示した通り、量子ゲートをひとつひとつ再現しているだけである。ただし、SWAPゲートは省略した。出力は以下の通り。

各基底状態(|00\rangle,|01\rangle,|10\rangle,|11\rangle)の係数(複素振幅)が出力される。これらは複素数である(Pythonでは虚数単位はjと表記される)。上の回路図に従って手計算を行うとすべての係数が0.5となるので(試してほしい)上の結果と一致している。

 次のコードは同じ2量子ビットのQFTであるが、初期状態が|j_0j_1\rangle=|11\rangleの場合である。

6行目で|j_0\rangle|j_1\rangleの両ビットにXゲートを作用させている。Xゲートはビットを反転させる量子ゲートである。

(42)    \begin{eqnarray*} X|0\rangle&=&|1\rangle \\ X|1\rangle&=&|0\rangle \end{eqnarray*}

従って、6行目で状態は|00\rangleから|11\rangleに変換される。以降の処理は先のコードと同じである。出力は以下の通り。

ここまでのコードでは、F_2|j\rangleまでの計算しかしていない。フーリエ係数まで求めるコードは以下の通りである。

このコードは一般のn量子ビットの場合のQFTである。入力がx=(x_0,\cdots,x_{N-1}))、出力がy=(y_0,\cdots,y_{N-1}))である。ここで、N=2^nである。7行目から11行目までのコードで量子回路Fを構築している。14行目から28行目までのコードで|j\rangle=|j_0\cdots j_{n-1}\rangleの取り得るすべての状態cを再現し、その各々をFに接続している(26行目:d=c+F)。31行目から35行目までの計算でフーリエ係数y_kを計算している。n量子ビットの場合の量子回路を再掲しておく。


上で実装した関数qftを用いて次式を離散フーリエ変換してみる。

(43)    \begin{equation*} y=5\sin{2x} \end{equation*}

これを離散化するため、0から2\piまでをN=2^4=16等分する。

(44)    \begin{equation*} x_j=5\sin{2\left(\frac{2\pi j}{16}\right)}, j=0,\cdots,15 \end{equation*}

このとき、フーリエ係数は

(45)    \begin{equation*} y_k=\frac{1}{4}\sum_{j=0}^{N-1}x_j\exp{\left(i\frac{2\pi kj}{16}\right)} \end{equation*}

となる。この結果をqftで用いて再現してみる。

1行目の関数は式(44)を、5行目の関数は式(45)を表す。これらの式を用いて計算したフーリエ係数と、上で実装したqftによる結果とを比較する。

ここでは、フーリエ係数の大きさ(|y_k|,k=0,\cdots,15)を出力している。結果は以下の通り。

どちらの結果も、0010(=2)番目と1110(=14)番目の項が10、その他の項は0となる。y_2y_{14}だけが値を持つのは元の式y=5\sin{2x}が周波数2の正弦波のためである(y_{14}y_{2}の複素共役である)。

まとめ

 今回は、量子コンピュータを用いた計算の具体例としてフーリエ変換を取り上げ、かなり詳しく解説した(かなり長い記事になってしまった)。量子コンピュータで計算するとは、アルゴリズムを量子回路で表現することである。次回は、回帰を量子コンピュータで計算する方法を説明する予定である。

参考文献

  • 量子コンピューティング 基本アルゴリズムから量子機械学習まで
  • 量子コンピュータ入門[第2版]
  • 量子フーリエ変換とは
  • blueqatで量子フーリエ変換
  • Kumada Seiya

    Kumada Seiya

    仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

    最近の記事

    • 関連記事
    • おすすめ記事
    • 特集記事

    アーカイブ

    カテゴリー

    PAGE TOP