量子傅立葉轉換(上)

作者:
林昱誠(Yu-Cheng Lin)、施宇宸、徐育兆
閱讀時間:
15
分鐘
# 量子傅立葉變換(上) 量子傅立葉轉換(Quantum Fourier Transform, 簡稱 QFT),QFT 在量子計算中扮演重大角色,在後面用於解 [eigenvalue](https://www.entangletech.tw/lesson/math-08) 的 QPE 演算法(可以應用於化學問題)以及用於質因數分解的 Shor 演算法(可以應用於破解密碼)都會見到 QFT 的身影。QFT 是從離散傅立葉轉換(Discrete Fourier Transform, DFT)演變過來,因此在這篇,我們會先花不少篇幅介紹 DFT,如此才能理解 QFT 在做什麼,在[下一篇](https://www.entangletech.tw/lesson/algorithm-06)我們會更深入了解 QFT 的運作。 ## 何謂傅立葉轉換? 傅立葉轉換(Fourier transform)是一種數學技巧,最常見的應用就是分析“波”,可以將一個看似毫無規律的波(或說函數)轉換成週期函數($\sin$ 函數、 $\cos$ 函數等等)的線性組合。 這段[影片](https://youtu.be/spUNpyF58BY?si=AqrYefY1PmMR-aWV)透過圖示化幫助讀者了解傅立葉轉換。
圖片內容

對於不同種類的函數(或說波)會有不同版本的傅立葉轉換,在這裡我們將針對離散函數(下面會有例子說明什麼叫「離散」),處理離散函數的傅立葉轉換稱作離散傅立葉轉換(Discrete Fourier Transform, DFT) ### 離散傅立葉轉換 DFT 最具體的應用就是分析聲音,將聲音的波形轉換成聲音的頻率。
圖片內容

訊號經 DFT 轉換後成頻率

下圖是以鋼琴彈出 C 大調 1 秒內的波形
C 大調

C 大調為期一秒的波形

由於電腦是透過數位方式紀錄聲音,因此上圖看似連續的波形其實是透過一筆一筆資料儲存起來(這就是"離散")
C 大調

C 大調為期 0.05 秒的波形

總共有 44100 筆資料,大概會長這樣: \begin{array}{|c|c|} \hline \text{時間} & \text{震幅} \\ \hline 0 & -0.46933 \\ \hline 0.00002 & -0.46011 \\ \hline 0.00005 & -0.44931 \\ \hline \vdots & \vdots \\ \hline 0.99998 & 0.13571 \\ \hline \end{array} 我們無法從上面這樣的波看出什麼有意義的資訊,看不出 C 大調是由哪些頻率的波組合而成,這時候我們就需要 DFT 幫我們做分析,容我們先上 DFT 的數學公式 \begin{align}\tag{1} y_k = \frac{1}{\sqrt{N}} \sum_{j=0}^{N-1} a_j e^{i \frac{2\pi}{N} jk} \end{align} 其中 $N$ 就是有幾筆資料,以我們這例子就是 $N=1440$, $a_j$ 就是第 $j$ 筆資料的數值(在這裡是波的振幅),比方說 $a_0=-0.46933$,接著我們從 $k=0$ 開始列出以下式子: \begin{split} y_0 &= \frac{1}{\sqrt{44100}} (-0.46933 e^{i\frac{2\pi}{1440}0\cdot0}-0.46011e^{i\frac{2\pi}{1440}1\cdot 0}+...+0.13571e^{i\frac{2\pi}{1440}44099\cdot 0})=−0.0973861 \\ y_1 &= \frac{1}{\sqrt{44100}} (-0.46933 e^{i\frac{2\pi}{1440}0\cdot1}-0.46011e^{i\frac{2\pi}{1440}1\cdot 1}+...+0.13571e^{i\frac{2\pi}{1440}44099\cdot 1})=−0.118737 + 0.136405i \\ y_2 &= \frac{1}{\sqrt{44100}} (-0.46933 e^{i\frac{2\pi}{1440}0\cdot2}-0.46011e^{i\frac{2\pi}{1440}1\cdot 2}+...+0.13571e^{i\frac{2\pi}{1440}44099\cdot 2})=−0.106039 + 0.0597867i \\ &\vdots \end{split} 我們將 $y_k$ 平方開根號,得出 \begin{split} |y_0|&=0.097386 \\ |y_1|&=0.180844 \\ |y_2|&=0.121732 \\ &\vdots \end{split} 以 $k$ 為橫軸,$|y_k|$ 為縱軸,繪製成以下稱作頻譜的圖
圖片內容
C 大調的頻譜分析

從這張圖我們可以清楚看到 C 大調其實是由非常多頻率的聲音組合而成,不過訊號最強的是 262 Hz、330 Hz 與 392 Hz,這三個頻率分別對應到 Middle C、E 與 G(音樂術語),其他明顯的波峰都是這三個主音的倍頻,像是 522 Hz 是 Middle C 頻率的 2 倍,784 Hz 則為其 3 倍,660 Hz 是 E 頻率的 2 倍,也是 G 的 2 倍,990 Hz 是 E 的 3 倍。 ## 經典電腦的解法 從 $(1)$ 式中我們知道計算每個 $y_k$ 需要加總 $N$ 項數值,而總共有 $N$ 個 $y_k$,因此全部總共有 $N^2$ 項,因此一項一項加起來的算法複雜度為 $O(N^2)$(以 C 大調這例子,前 11000 筆資料,筆者花了 1 小時做計算,程式碼見附錄)。 其實我們也能把 $(1)$ 式想成一個矩陣乘法,為了方便書寫,我們定義 $\omega=e^{i\frac{2\pi}{N}}$,那 $(1)$ 式就變成 \begin{split} y_k=\frac{1}{\sqrt N}\sum_{j=0}^{N-1} a_j\omega^{jk} \end{split} 將其展開變成 \begin{split} y_0 = \frac{1}{\sqrt{N}} (a_0 + a_1 + a_2 + \cdots + a_{N-1}) \end{split}
這就是所有筆資料的平均值
\begin{split} y_1 = \frac{1}{\sqrt{N}} (a_0+a_1\omega+a_2\omega^2 \cdots + a_{N-1}\omega^{N-1}) \end{split}
可以看出其 y0 與 y1 的 phase 項差一個 omega
\begin{split} y_2 = \frac{1}{\sqrt{N}} (a_0+a_1\omega^2+a_2\omega^4 \cdots + a_{N-1}\omega^{2(N-1)}) \end{split}
可以看出其 y2 與 y0 的 phase 項相差 2 倍
將上式寫成矩陣 \begin{align}\tag{2} \begin{bmatrix} y_0 \\ y_1 \\ \vdots\\ y_{N-1} \end{bmatrix}= \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega & \omega^2 & \cdots & \omega^{N-1} \\ 1 & \omega^2 & \omega^4 & \cdots & \omega^{2(N-1)} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & \omega^{N-1} & \omega^{2(N-1)} & \cdots & \omega^{(N-1)^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots\\ a_{N-1} \end{bmatrix} \end{align} 這種視覺上看起來很舒適的矩陣,對於電腦而言也需要 $O(N^2)$ 步,不過不必灰心,有種演算法稱作 fast Fourier transform(FFT)可以將複雜度降低到 $O(N\log N)$,不過這不是今天的重點,我們就不針對這演算法做介紹。 ## 量子電腦的解法 從 $(2)$ 式中我們看到中間 $N\times N$ 矩陣是一個 unitary 的矩陣,這代表我們可以用量子邏輯閘組合出這個矩陣,以做到加速的效果,這就是量子傅立葉轉換(QFT)。 對應到量子電路,$(2)$ 式告訴我們這電路的初始態為 \begin{align} |a\rangle= \begin{bmatrix} a_0 \\ a_1 \\ \vdots\\ a_{N-1} \end{bmatrix}= a_0|0\rangle+a_1|1\rangle+a_2|2\rangle+\cdots+a_{N-1}|N-1\rangle= \sum_{j=0}^{N-1} a_j|j\rangle \end{align} 其中 \begin{align} |0\rangle= \begin{bmatrix} 1 \\ 0 \\ 0 \\ \vdots \end{bmatrix} , |1\rangle= \begin{bmatrix} 0 \\ 1 \\ 0 \\ \vdots \end{bmatrix} , |2\rangle= \begin{bmatrix} 0 \\ 0 \\ 1 \\ \vdots \end{bmatrix} \end{align} 經過中間 $N\times N$ 矩陣或說量子邏輯閘操作後,電路輸出 \begin{align} |y\rangle= \begin{bmatrix} y_0 \\ y_1 \\ \vdots\\ y_{N-1} \end{bmatrix}= y_0|0\rangle+y_1|1\rangle+y_2|2\rangle+\cdots+y_{N-1}|N-1\rangle=\sum_{j=0}^{N-1} y_k|k\rangle \end{align} 這時我們就稱 $|y\rangle$ 是 $|\psi\rangle$ 的量子傅立葉轉換。結合 $(1)$ 式: \begin{align}\tag{1} y_k = \frac{1}{\sqrt{N}} \sum_{j=0}^{N-1} a_j e^{i \frac{2\pi}{N} jk} \end{align} 以下這條數學式就是 QFT 要達成的事情: \begin{align} QFT|a\rangle=|y\rangle=\sum_{k=0}^{N-1} y_k|k\rangle= \frac{1}{\sqrt N}\sum_{k=0}^{N-1}\sum_{j=0}^{N-1}a_je^{i\frac{2\pi}{N}jk}|k\rangle \end{align} 不難看出 \begin{align}\tag{3} |j\rangle \rightarrow\frac{1}{\sqrt N}\sum_{k=0}^{N-1}\omega^{jk}|k\rangle \end{align} \begin{align}\tag{4} y_k=\frac{1}{\sqrt 2}\sum_{j=0}^{N-1}a_j \omega^{jk} \end{align} 在 DFT 中,$N$ 對應到資料筆數,在量子電腦中,$N=2^n$,$n$ 是 qubit 數,接著我們來看 QFT 運作起來的效果是什麼: ### 一個 qubit 我們先從最簡單的系統開始,即只有一個 qubit,即 $n=1$,那 $N=2$,$\omega=e^{i\frac{2\pi}{2}}=e^{i\pi}=-1$,初始態 $|a\rangle=a_0|0\rangle+a_1|1\rangle$
倒數第二句話用到歐拉公式
因此從(4)式: \begin{align} y_k&=\frac{1}{\sqrt 2}\sum_{j=0}^{N-1}a_j \omega^{jk} \\ &=\frac{1}{\sqrt 2}[a_0(-1)^{0k}+a_1(-1)^{1k}] \\ &=\frac{1}{\sqrt 2}[a_0+a_1(-1)^k] \end{align} 當電路初始態為 $|a\rangle=|0\rangle$,即 $a_0=1, a_1=0$,上式則為 \begin{align} y_0=y_1=\frac{1}{\sqrt 2} \end{align} 當電路初始態為 $|a\rangle=|1\rangle$,即 $a_0=0, a_1=1$,則為 \begin{align} y_0=\frac{1}{\sqrt 2} \\ y_1=-\frac{1}{\sqrt 2} \end{align} 綜合以上: \begin{align} QFT|0\rangle=\frac{1}{\sqrt 2}(|0\rangle+|1\rangle) \end{align} \begin{align} QFT|1\rangle=\frac{1}{\sqrt 2}(|0\rangle-|1\rangle) \end{align} 有沒有覺得這個哪裡似曾相似,沒錯,就是 H gate 的作用,所以在只有一個 qubit 的時候,QFT 長這樣
H gate

n=1 的 QFT 電路圖

將這 QFT 的作用寫成矩陣就是: \begin{split} \frac{1}{\sqrt 2} \begin{bmatrix} 1 & 1 \\ 1 & \omega \end{bmatrix}= \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \end{split} 剛好是 H gate 對應的矩陣,也剛好對應(2)式中的矩陣 ### 兩個 qubits 此時 $n=2, N=3$,初始態與輸出態為 \begin{split} |a\rangle=a_0|0\rangle+a_1|1\rangle+a_2|2\rangle+a_3|3\rangle \\ |y\rangle=y_0|0\rangle+y_1|1\rangle+y_2|2\rangle+y_3|3\rangle \end{split} 由於兩個 qubits 的情況就開始變複雜,為了方便起見,我們將 $|a\rangle$ 寫成[二進位制](https://www.entangletech.tw/lesson/basic-algorithm-02) $|j_1j_0\rangle$,$j_1$ 與 $j_0$ 可以是 $0$ 也可以是 $1$,今天要轉回去十進制就 $2j_1+j_0$,可以自己試試看,同理輸出態 $|y\rangle = |k_1k_0\rangle$,然後從(3)式出發(免去係數 $a$ 的煩惱): \begin{align}\tag{3} |j\rangle \rightarrow\frac{1}{\sqrt N}\sum_{k=0}^{N-1}\omega^{jk}|k\rangle \end{align}
如果數學是你的罩門,中間推導可以跳過,直接看電路圖長什麼樣子
所以: \begin{align} QFT|j_1j_0\rangle=\frac{1}{\sqrt 4}\sum_{k=0}^3 e^{i2\pi\frac{1}{4}j(2k_1+k_0)}|k_1k_0\rangle \end{align}
因為四式中的 k 是十進制,要從二進制 k1k0 轉回去就是 2k1_k0,同理 j 也是
上式接續下去 \begin{align} &=\frac{1}{\sqrt 4}\sum_{k=0}^3 e^{i2\pi\frac{1}{4}j(2k_1+k_0)}|k_1k_0\rangle \\ &= \frac{1}{2}\sum_{k_1=0}^1 e^{i2\pi 2\frac{1}{4}jk_1}|k_1\rangle \otimes \sum_{k_0=0}^1 e^{i2\pi \frac{1}{4}jk_0}|k_0\rangle \\ &= \frac{1}{2}(|0\rangle+e^{i\pi j}|1\rangle) \otimes (|0\rangle+e^{i\pi \frac{1}{2} j}|1\rangle) \\ &=\frac{1}{2}(|0\rangle+e^{i\pi (2j_1+j_0)}|1\rangle) \otimes (|0\rangle+e^{i\pi \frac{1}{2} (2j_1+j_0)}|1\rangle) \\ \end{align} 最後一行第一項可以拆解成 $e^{2i\pi j_1}e^{i\pi j_0}$,當 $j_1=1$ 時 $e^{2i\pi}=1$(當 $j_0=0$ 時也是 $1$),所以可以簡化成 $e^{i\pi j_0}$ \begin{align} &=\frac{1}{2} (|0\rangle+e^{i\pi j_0}|1\rangle) \otimes (|0\rangle+e^{i\pi j_1}e^{\frac{1}{2}i\pi j_0}|1\rangle) \\ &=\frac{1}{2} (|0\rangle+e^{i\pi j_0}|1\rangle) \otimes e^{\frac{1}{2}i\pi j_0} (|0\rangle+e^{i\pi j_1}|1\rangle) \\ &=\frac{1}{2} (|0\rangle+(-1)^{j_0}|1\rangle) \otimes S^{j_0} (|0\rangle+(-1)^{j_1}|1\rangle) \end{align} 如果以上看得眼花撩亂,那就跳過,但注意這個,這邊出現一個新的代數 $S^{j_0}$,他的效果是: \begin{split} \text{當 }j_0=0,S^{j_0}=e^0=I \\ \text{當 }j_0=1,S^{j_0}=e^{\frac{1}{2}\pi i}=S \\ \end{split} 所以他是一個 controlled S gate(CS gate),因此上面結果可以寫成 \begin{split} &=\frac{1}{2} (|0\rangle+(-1)^{j_0}|1\rangle) \otimes S^{j_0} (|0\rangle+(-1)^{j_1}|1\rangle) \\ &= (H\otimes I)CS(I\otimes H)|j_0j_1\rangle \\ &=(H\otimes I)CS(I\otimes H)SWAP|j_1j_0\rangle \end{split} 所以 $n=2$ 的 QFT 電路圖為:
n=2 QFT

n=2 的 QFT 電路圖

### 三個 qubits 我們知道讀者累了,這邊就不推導,有興趣的可以看附錄,他的電路圖長為這樣:
n=2 QFT

n=3 的 QFT 電路圖

### 多個 qubits 以此類推到多個 qubits,我們可以從前面三個例子推估完整的 QFT 電路應該會長這樣:
QFT
QFT 電路圖

有些教材的電路圖剛好與本圖上下相反,就只要把輸入的順序上下相反就好

我們將在下下一篇解釋這電路怎麼來的,圖中 $R_k$(雖然圖中是 $R_n$) 是一種 controlled phase rotation gate,其矩陣是: \begin{align} R_k= \begin{bmatrix} 1 & 0\\ 0 & e^{\frac{2\pi i}{2^k}} \end{bmatrix} \end{align} 它產生的效果會是 \begin{split} R_k |0\rangle &= |0\rangle \\ R_k |1\rangle &= e^{\frac{2\pi i}{2^k}}|1\rangle \end{split} 相信你在以上公式中反覆看到 $e^{i\frac{2\pi}{N}}$ 的東西而感到厭煩,我們將在[下一篇](https://www.entangletech.tw/lesson/algorithm-06)介紹這指數因子的作用。 ## 量子加速 QFT 電路中第一個 qubits 總共用了 $n$ 個 gates,第二個則用了 $n-1$ 個 gates,以此類推,因此 $n$ 個 qubits 將會用到: \begin{split} n+(n-1)+(n-2)+\cdots+3+2+1=\frac{n(n+1)}{2} \end{split} 個 gates,另外會有 $\frac{n}{2}$ 個 SWAP gate,所以總共需要 \begin{split} \frac{n(n+1)}{2}+\frac{n}{2} \end{split} 個 gate,以複雜度來說是 $O(\log^2N)$,比傳統 DFT 的 $O(N\log N)$ 還要快,達到量子加速。 ## QFT 的作用 如同前面介紹的 DFT,可以把聲波轉換成頻譜,以幫助我們更了解聲音的組成,QFT 也是,不過今天它不是要轉成我們人類看得懂的“譜“,而是量子電腦看得懂的。我們知道量子電腦是透過「疊加態」(與糾纏態)來做到更快的運算,然而我們要處理的問題或數據都不是疊加態下的資料,都是我們熟悉的二進位制,因此我們需要透過 QFT 將我們人類看得懂的二進位制轉成量子電腦看得懂的資料,藉此作運算,我們將在[下一篇](https://www.entangletech.tw/lesson/algorithm-06)用 Bloch sphere 圖示化這點。 那量子電腦運算完之後的結果還是在疊加態,我們需要一個方法轉成人類看得懂的二進位制,此時就需要 inverse QFT(逆 QFT)。這就類似我們把聲音轉成眼睛看得懂的頻譜後,我們在頻譜上某個頻率添加一個 peak,但我們並不知道加一個 peak 後的這頻譜聽起來會是什麼聲音,這時我們就需要 inverse DFT 把頻譜轉回去聲音,並播放給耳朵聽。 ## Inverse QFT 有時候我們會需要把 $|y\rangle$ 轉回 $|a\rangle$,以前面 C 大調為例,就是我們想知道這樣的頻譜會產生什麼樣的聲音,這時候就需要逆轉換,即: \begin{split} \frac{1}{\sqrt N}\sum_{k=0}^{N-1}e^{2i\frac{\pi}{N}jk}|k\rangle\rightarrow |j\rangle \end{split} 就很簡單地把 QFT 電路左右倒過來就是了。 ## 附錄 ### FT 程式碼 ```python=+ import pandas as pd import numpy as np N = len(amplitude) # 資料點數44100 dft =[] for i in range(N): phi=0 for k in range(N): phi += amplitude[k]*np.exp(-2j*np.pi*i*k/N) phi = phi*(1/np.sqrt(N)) dft.append(phi) dft_norm = np.abs(dft) ``` ### n=3 QFT \begin{split} QFT|j_2j_1j_0\rangle&=\frac{1}{\sqrt 2^3}(|0\rangle+e^{i\pi j_0}|1\rangle)\otimes (|0\rangle+e^{2i\pi(\frac{j_1}{2}+\frac{j_0}{4})}|1\rangle) \otimes (|0\rangle+e^{2i\pi (\frac{j_2}{2}+\frac{j_1}{4}+\frac{j_0}{8})}|0\rangle) \\ &= \frac{1}{\sqrt 8}[|0\rangle+(-1)^{j_0}|1\rangle]\otimes S^{j_0}[|0\rangle+(-1)^{j_1}|1\rangle] \otimes T^{j_0}S^{j_1}[|0\rangle+(-1)^{j_2}|1\rangle] \\ &=(H\otimes I \otimes I)S(I\otimes H\otimes I)TS(I\otimes I\otimes H)|j_0j_1j_2\rangle \\ &=(H\otimes I \otimes I)S(I\otimes H\otimes I)TS(I\otimes I\otimes H)SWAP|j_2j_1j_0\rangle \end{split}
課程目錄