フーリエさん、ひさしぶり？
このワークシートはMath by Codeの一部です。
今回は、量子フーリエ変換QFTにとりくもう。
その前に、
フーリエ級数からフーリエ変換への流れを確認しておこう。
フーリエ級数というのは、周期関数を三角関数の和に展開するものだったね。
＜フーリエ級数を思い出そう＞
一般に、周期Tのｆ(t)のフーリエ級数は
f(t)=a0+ Σ( an cos[ 2π n t/T ] + bn sin[ 2π n t/T] )
a0 = 1/T ∫ f(t) dt (-T/2, T/2)
an = 2/T ∫ f(t) cos[ 2π n t/T ] dt ( -T/2, T/2)
bn = 2/T ∫ f(t) sin[ 2π n t/T ] dt (-T/2, T/2)
とされていて、音の波形を作るのに役立ったね。
純音はｙ＝a sin[ 2π f t] （f=440でA４の音）だった。
フーリエ級数展開では、周期関数f(t)と倍音周波数の純音に和分解していることになる。
ノコギリ波は[ーπ、＋π]の区間で、f(t)=t, T=2πのときのフーリエ級数展開は
f(t) =2 Σ(-1)(n+1)/n sin(nt)
矩形波は[0, π]区間で作れたら奇関数で、f(t)=1, T=2πのときのフーリエ級数展開は
f(t) = Σ(1-(-1)n) /nπ sin(nt)
三角波はf(t) = t (0からπ/2) ；π-t (π/2, π）を偶関数化したフーリエ級数展開は
f(t)=π/4- 2/π{ cos 2t + 1/32 cos 6t +1/ 52 cos 10t + ......1/ (2m+1)2 cos2(2m+1)+.......}
基本周波数をｆとすると、(周波数の倍率,振幅の倍率)のペアのリストができる。
純音 ；(１，１）
ノコギリ波；(1,1),(2,-1/2), (3, 1/3),(4,-1/4),(5,1/5),........ ( k, (-1)^(k+1)/k ),.....
矩形波 ；(1,1) , (3, 1/3) ,(5,1/5),....., ( k, (1-(-1)^k)/k ),.....
三角波 ；(1,1), ,(3, 1/32) ,(5,1/52),,.... (2k-1, 1/(2k-1)^2),.....
これらを決めるには、部分積分と三角関数の０、πの値が0や-1になることを利用したね。
＜複素数化と区間の無限化でフーリエ変換＞
また、オイラーの公式により、複素数を使うとフーリエ級数が複素フーリエ級数になった。
式が見やすくなり、f(t)に指数関数e^ixをかけた形でシンプルになったね。
複素フーリエ級数f(t) =...+c_-1e^(-ik-1t)+c0+c_1 e^(ik1t)+......= Σc_n e^(iknt)(n from ‐∞ to ∞)
複素フーリエ係数cn=1/T ∫ f(t) e^(-iknt) dt( -T/2, T/2)
さらに、複素フーリエ級数の周期T→∞にしたものがフーリエ変換と逆フーリエ変換だ。
f(t)のフーリエ変換はF(k)=∫ f(t) e^(-ikt) dt (‐∞,∞)
F(k)の逆フーリエ変換はf(t)= 1/2π∫F(k) e^(ikt) dk (‐∞,∞)
関数を波のベクトルとみなすと、
フーリエ級数が直交基底ベクトルの和でできているといえる。
{Un(x)}=１/√π{１/√2, cosx ,sinx, cos2x, sin2x, ........} も正規直交する基底関数だ。（ｘが-π以上とπ未満）
f(x)=Σ am | Um＞ ( m=[0、∞] )のように2π周期関数f(x)がUnの級数として表示できる。
これがフーリエ級数だった。
フーリエ級数が正規直交関数系（直交多項式）でできてることが、
フーリエ変換や微分方程式の解法につながった。
（例）
フーリエ級数y=1・sinx+1/3 ・sin3x+ 1/5・sin5x は、
３つの正規直交基底Fm={sin x, sin3x, sin5x}の一次結合。
（正規直交関数である理由）
(sinx・sin3x)=(sin3x・sin5x)=(sin5x・sinx)=0
(sinx・sinx)=(sin3x・sin3x)=(sin5x・sin5x)=π。
∫XYdx[ -π, π] は奇関数X×奇関数Y＝奇関数の積分は原点対称で０になる。
X=mxとして、∫(sin X)2 dx=1/2 ∫(1- cos 2X) dx =1/2 ∫1 dx-1/2 ∫cos 2Xdx=1/2[x][-π,π] -0=π。
関数を波のベクトルとみなせば、
f(t)のフーリエ変換はF(k)=∫ f(t) e^(-ikt) dt (‐∞,∞)
F(k)の逆フーリエ変換f(t)= 1/2π∫F(k) e^(ikt) dk (‐∞,∞)
この２式が内積に見えてくるね。
f(t)という基本の時間波に周波数kごとの回転をかけてたし算する内積だ。
波にその周波数があれば強められてなければ弱められます。波の干渉性ですね。
これって、量子論理とスゴク相性がいいですよね。
ヒルベルト空間を思い出すと内積がピンとくるでしょう。
逆変換は周波数ごとの波F(k)に時間tの逆回転をして、それを合計すると、もとの基本波ｆ（t）に戻るということだ。
だから、連続フーリエ変換はｆ(t)という時間領域の非周期関数から、F(k)という周波数領域の連続関数への変換をするので、フーリエ変換は周波数分解、スペクトル分解と言い換えることもできる。
「周波数を横軸にとると、特定の周波数ｋの成分がどれくらい含まれているかが視覚化」できます。
もっと、実用的にかみ砕くと、Waveファイルから音高データを導くので、Midiファイル作成にもつながるね。
＜DFTとFFT＞
この関数の内積（F=∫f(t)*e dt ）をデジタルな内積(y=Σx*e ）にしたものが
離散フーリエ変換（DFT)だ。
調べる時間がＴ秒間で、サンプリング間隔がＤ秒だとするとn=T/D。
時刻リストts={0,D,2D,....(n-1)D}に対するf(t)のデータx={x0, x1,x2,....x(n-1)}が得られたとしよう。
ｆ(t)が連続的なときと同様にxt =...+c-1e-ik-1t+c0+c1eik1t+......= Σc_n e^(iknt) (‐∞,∞)
４つの時刻に対する音や波の信号データから、フーリエ係数c0,c1,c2,c3を逆算して、
求めることができるはずだ。
ベクトルと行列で再表現してみよう。
この式を、ベクトルA＝[c0,c1,c2,c3]からベクトルx=[x0,x1,x2,x3]を求める
かけ算ｘ＝Ａｃとするとき、これが逆離散フーリエ変換（IDFT）で、
この逆算の c=A^-1 xが離散フーリエ変換(DFT)です。
そして、DFTを高速化するものがFFTだったね。
FFTのアルゴリズムの基本は、
ｘをLとRに分けて、対応する番号順の和と差にし、
差については、回転子にリストの番号乗をしてかけることです。
これをn=8,n=4,n=2に対して再帰的に実行しよう。名付けて、再帰分割（バタフライ演算）。
すると、
時間順に区切った音波の振幅データ
x=[x0,x1,x2,x3,x4,x5,x6,x7]に対して
X=[X0,X4,X2,X6,X1,X5,X3,X7]を順に求めることができますね。
このインデックスの入れ替わりは2進数にしたときに左右逆転する数同士になっていること（ビット反転）
を使うと、計算で求めることもできるはず。
このXをインデックス順に並べ替えることで、
[X0,X1,X2,X3,X4,X5,X6,X7]が周波数を小さい順にならべた周波数特性のデータになります。
たとえば、n=8の場合複素数のかけ算が8*8=64個出てきます。
FFTを使うと、n=4にsizeダウンしたことで、8/2=4個。n=2にサイズダウンしたことで、4/2=2個を2回。
結果として、4+2*2=4×2=8個ですみます。
一般に、サイズがnなら複素数のかけ算の回数がn^2回からFFTでn/2×(log2(n)-1)回に激減します。
N*NのオーダーからN*logNのオーダーに減るということですね。
＜FFTからQFT＞
フーリエ変換の関数の内積（F=∫f(t)*e dt ）をデジタルな内積(y=Σx*e ）するとどうなったでしょうか。
連続フーリエ変換はF(k)=∫ f(t) e^(-ikt) dt (‐∞,∞) は時間の波関数に周波数の波をかけて積分しました。
離散フーリエ変換は時間による波ベクトルｘjから周波数ベクトルｙkを求めるために
周波数の波eをかけて総和を求めます。ｘ、ｙともにサイズはNで添え字は０からN-1としましょう。
y_k=１/√N Σ x_j e^(2π ijk/N)
ここで、ω=e^2πi /N とすると、ω^jk= e^(2π i jk/N) となるね。
だから、FTはN行N列の行列でかくことができる。列ベクトル＝|m)=(1,ω,ω^2,....,, ω^N-1)とおくと、
FT=１/√N[ |m^0), |m^1), |m^2),......,|m^N-1))
これは基底１つの変換だったから、重なり状態の波については次のようになる。
FT: Σ xj |j) →Σ yk |k)
y=FT xのような行列方程式になるといえるね。
FTはユニタリー行列で、状態ベクトルの回転を表している。
この計算の工夫が、QFTのポイントになるということだね。
N=2^nとすれば、
ωn=e^2πi /N= e^(2πi /2^n)
ωl=e^(2πi /2^l)
ω=e^2πi
とおくこともできて、
ωnは複素平面上の単位円のN等分したベクトル
ωlは複素平面上の単位円の2^l等分したベクトルと見られる。
ωはベクトルにかけると回転子となり、波に干渉することができる。
Nをnビットの2進数とみることもできる。
時間をデジタルにN個に切った基底ベクトル|j)を
N個の周波数の基底ベクトル|k)に対して周波数{k)のベクトルをかけた内積を求めるのがQFTだ。
周波数の基底ベクトル|k)は、|k_1)から|k_n)までのｎ個のビットに区切って考えることができる。
だから、
FFTでN＝８なら8→４→2と
周波数基底による線形結合を入れ子にして、
ペアのペアのぺア
{[(x0-x4)+p2(x2-x6)]+[p1(x1-x5)+p3(x3-x7)]}
のように、２ビットまで2等分を繰り返して、
「多層的に表現」できる。
再帰的表現の代わりに指数でk1ビットが２を多くかけることになるね。
ただし、k1からknというのはｎビットのインデックスの意味であり、実際は０か１でしかない。
だから、最後はどのビットも０か１の2通りに分解できる。
代数式とみると、「入れ子の式は因数分解で2項式のかけ算」に直せる。
また、ｊ/2^nを2進小数j_1j_2....j_n-1j_nと表すことができる。
回転は1回転を超えると切り捨てられることを考えると、
２進小数0.j_1j_2....j_n-1j_nを2倍2倍していき、整数を切り捨てていくと、ｋ1の回転子の指数は0.j_nだけになる。
FT:
|j)=｜j_1,j_2,....,j_n-1,j_n)
->１/√N Σ ωn^jk |k)
=１/√N Σ_(k1) …Σ_(kn) ωn^j (k1*2^(n-1)+....kn*2^0) |k1...kn)
=１/√N (Σ_(k1)ω1^k1 |k1)) … ( Σ_(kn) ωn^kn |kn) )
=１/√N (|0) +ω^0.j_n|1)) (|0) +ω^0.j_n-1j_n |1)) .... (|0) +ω^0.j_1j_2....j_n-1j_n|1))
ビットごとに独立、並列して計算ができることがわかった。
あとは回転子の指数にある小数の計算を回路でどう実現するかという課題があるね。
＜QFTの回路＞
二進数はオンオフと親和性が高いですね。
しかも、回転するのは|1)です。
回転の最小単位はj/2^nです。これを位相ゲートRn=P(n)=[(1,0),)(0,e^2πi/2^n)]としましょう。
第kのビットがたっているときだけ、Rk=P(k)=[(1,0),)(0,e^2πi/2^k)]まわす制御位相ゲートをCRkとします。
第１ビットが１のときだけ、R1回すのは、Hでできますから、
ｊが３ビットのときは、
H, R2,R3を1ビット目にして、
2ビット目にはH,R2をして、
３ビット目にはHをします。
合計で３＋２＋１＝６個のゲートで３つのビットに回転をかけられますね。
古典論理のFFTのときに、インデックスの大小を最後に入れ替えて整理したように、
量子論理でも、インデックスのソートをします。ソートは大小のSWAPを使いましょう。
SWAPは前にやりましたが、入れ替える２つのワイヤでCNOT3連発すればよかった。
結局、ｎビットのQFTのゲート数は回転でn+n-1+....+1=n(n+1)/2個横に並べ、
その後、大小整理で３*floor(n/2) 個。
合計でn(n+1)/2＋3floor(n/2)個のゲートで実現できるね。
FFTがｎ*2^ｎのオーダーの計算量レベルだったのに対して、
QFTはｎ*ｎのオーダーに減る。
QFTから見れば、FFTは指数関数的な計算量だ。
しかし、QFTは理論的にはすばらしいが、量子コンピュータでは
実際に音波に直接使うには音の振幅を直接測定することが難しいようです。
＜QFTの振り返りとユースケース＞
QFTを振り返ると、
まず、時間の関数である波をデジタル化して、2進数として基底を渡します。
FTでは、周波数を2進化して、基底と回転をセットにして、波にかけます。
この内積によって、周波数のあるところは強め、ないところはゼロにすることで、
周波数スペクトルを抜き出したのです。
時間の波から周波数を抜き出したということですね。
これが、FFTよりも驚異的な高効率ということだった。
このQFTのアルゴリズムは、位相推定、素因数分解、位数発見、暗号解読、
隠れ部分群の発見など、高難度で計算量の多い課題に対して
その有効性が確認されているようです。
課題：DFTとFFTとQFTの計算量のオーダーで視覚化しよう。
N=２^nとすると、
DFTはN＊Nのオーダー、FFTはn×Nのオーダー、QFTはｎ×ｎのオーダーです。
N^2、NlogN, (logN)^2のグラフをかいてみよう。