有點感覺了，趕緊記錄。
量子傅立葉轉換 (Quantum Fourier Transform, QFT) 是在量子運算上實現離散傅立葉轉換 (Discrete Fourier Transform, DFT)。QFT 的存在是為了之後 shor 演算法1，不過本篇先專注討論 QFT。
有錯再麻煩指正，感謝。
另外，本文強調使用實際例子，所以用的規模比較小，一般化的公式及程式再請大家參考所附的資料或等待之後的文章。
常用符號
- $\log$: 資訊領域多用以 $2$ 為底，也就是 $\log_2$，不過本文使用 $\log$。
- $i$: 虛數單位，也就是 $\sqrt{-1}$。
Discrete Fourier Transform
雖然小時候大一普物 (一) 有使用 Matlab 實作傅立葉級數，傅立葉轉換的部分在大二下修電機系的通訊導論也有遇到，不過還是有些不懂，所以就順便記錄。
定義
一般 DFT (classical Discrete Fourier Transform) 是針對一組長度為 $N = 2^n$ 的複數向量 $(x_0, x_1, …, x_{N-1})\in\mathbb{C}^N$ 轉換，其對應的結果也是一組複數向量 $(y_0, y_1, …, y_{N-1})\in\mathbb{C}^N$。函數對應如下:
$$
\begin{equation*}
\begin{split}
y_k =& \frac{1}{\sqrt{N}}\sum^{N-1}_{n = 0}{x_n{\omega_N}^{-kn}}, k = 0, 1, 2, …, N-1
\newline
\newline
\omega_N =& e^{\frac{2 \pi i}{N}}
\newline
{\omega_N}^k =& e^{{\frac{2 \pi i}{N}}k}, \text{Root of unity}
\end{split}
\end{equation*}
$$
根據尤拉公式2，可以將公式表示為:
$$
y_k = \frac{1}{\sqrt{N}}\sum^{N-1}_{n = 0}{x_n\big(\cos(\frac{2\pi}{N}kn)-i\sin(\frac{2\pi}{N}kn)\big)}
$$
對 DFT 來說，正規化並不重要，不過為了方便跟量子形式，$\frac{1}{\sqrt{N}}$ 做的事就是使此函數為么正變換3。一般 DFT 的物理意義如下圖所示，將時間軸的資料作為輸入，接著輸出不同頻率下的發生率。
與單位根的關係
傅立葉轉換利用的複數平面的旋轉特性。其中與 Root of unity 單位根4有關，其定義為 $z^n = 1$ 在複數平面上有 n 次單位根，這裡有一個旋轉及取模的概念也呈現在下圖中。同時注意若 $kn$ 前有負號的話就會順時針取根如下圖；反之，沒有負號就會逆時針，如稍後會提到的 IDFT。
MarekSchmidt, Public domain, via Wikimedia Commons link
離散傅里葉變換矩陣 (DFT Matrix)
$\omega_N$ 所構成的運算可以看作一個矩陣，離散傅里葉變換矩陣5，$N = 4, (\omega = e^{-\frac{2\pi i}{4}} = -i)$ 時如下:
$$
M_{DFT} = {\frac{1}{\sqrt{N}}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & \omega & \omega^2 & \omega^3\\ 1 & \omega^2 & \omega^4 & \omega^6\\ 1 & \omega^3 & \omega^6 & \omega^9\end{pmatrix}} = {\frac{1}{\sqrt{4}}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & \omega & \omega^2 & \omega^3\\ 1 & \omega^2 & 1 & \omega^2\\ 1 & \omega^3 & \omega^2 & \omega\end{pmatrix}} = {\frac{1}{2}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & -i & -1 & i\\ 1 & -1 & 1 & -1\\ 1 & i & -1 & -i\end{pmatrix}}
$$
$$\small
DFT: \begin{pmatrix}y_0 \\ y_1 \\ y_2 \\ y_3\end{pmatrix} = M_{DFT} \begin{pmatrix}x_0 \\ x_1 \\ x_2 \\ x_3\end{pmatrix}
$$
實際運算把輸入乘上該矩陣即可得到結果 (範例來自: The DFT Output and Its Dimensions，為了方便把向量轉置):
$$
\begin{equation*}
\begin{split}
\begin{pmatrix}0 & 1 & 0 & -1\end{pmatrix}^T&\xrightarrow{DFT}\begin{pmatrix}0 & -2i & 0 & 2i\end{pmatrix}^T
\newline
\newline
\begin{pmatrix}1 & 0 & -1 & 0\end{pmatrix}^T&\xrightarrow{DFT}\begin{pmatrix}0 & 2 & 0 & 2\end{pmatrix}^T
\end{split}
\end{equation*}
$$
DFT Matrix 視覺化呈現如下圖:
John Ringland, CC BY-SA 4.0, via Wikimedia Commons link
此外，我這邊極度推薦看 3Blue1Brown “利用旋轉找出質心” 的高品質視覺化解釋影片: But what is the Fourier Transform? A visual introduction. - youtube 和續集 The more general uncertainty principle, beyond quantum - youtube (就不附在本篇，自行前往觀看)。
IDFT
DFT (時間到頻率) 的逆轉換是 IDFT (Inverse DFT)，也就是頻率到時間。其定義跟 DFT 在次方上差一個負號:
$$
\begin{equation*}
\begin{split}
y_k =& \frac{1}{\sqrt{N}}\sum^{N-1}_{n = 0}{x_n{\omega_N}^{kn}}, k = 0, 1, 2, …, N-1
\newline
\newline
\omega_N =& e^{\frac{2 \pi i}{N}}
\newline
{\omega_N}^k =& e^{{\frac{2 \pi i}{N}}k}, \text{Root of unity}
\end{split}
\end{equation*}
$$
$N = 4, (\omega = e^{\frac{2\pi i}{4}} = i)$ 時 IDFT Matrix 如下:
$$
M_{IDFT} = {\frac{1}{\sqrt{N}}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & \omega & \omega^2 & \omega^3\\ 1 & \omega^2 & \omega^4 & \omega^6\\ 1 & \omega^3 & \omega^6 & \omega^9\end{pmatrix}} = {\frac{1}{\sqrt{4}}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & \omega & \omega^2 & \omega^3\\ 1 & \omega^2 & 1 & \omega^2\\ 1 & \omega^3 & \omega^2 & \omega\end{pmatrix}} = {\frac{1}{2}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & i & -1 & -i\\ 1 & -1 & 1 & -1\\ 1 & -i & -1 & i\end{pmatrix}}
$$
以上是大部分人使用的版本，我取自下方來源:
- Discrete Fourier transform - wiki
- Quantum Fourier transform - wiki
- Introduction to the Fourier Transform (Part 1) - youtube
- Discrete Fourier Transform (numpy.fft) - numpy
- 2.161 Signal Processing: Continuous and Discrete Fall 2008 LECTURE 10 - MIT OpenCourseWare
另外有發現有些文章會將 DFT, IQFT 的運算對調以便說明 QFT。例如 IBM 教學中: Quantum Fourier Transform 就把 DFT 的負號移除。
Quantum Fourier Transform
量子傅立葉轉換也是一個對應關係，原先對資料旋轉的操作在量子上改成對 Z 軸來旋轉
- 使用 $n = \log N$ 個量子位元的疊加態，來疊加出 $2^n = N$ 種可能，由這 $N$ 種可能當作輸入 $\lvert x\rangle = \sum^{N-1}_{j = 0}x_j\lvert j\rangle$，$x_j \in \mathbb{C}$ 且 $\lvert {x_j}\rvert^2$ 代表第 j 個 state $\lvert j\rangle$ 出現的機率。
- 輸出為 $\lvert y\rangle = \sum^{N-1}_{j = 0}y_j\lvert j\rangle$，與輸入的設定相同。
- 其量子形式的函數定義如下:
$$
\begin{equation*}
\begin{split}
y_k =& \frac{1}{\sqrt{N}}\sum^{N-1}_{n = 0}{x_n{\omega_N}^{kn}}, k = 0, 1, 2, …, N-1
\newline
\newline
\omega_N =& e^{\frac{2 \pi i}{N}}
\newline
{\omega_N}^k =& e^{{\frac{2 \pi i}{N}}k}, \text{Root of unity}
\end{split}
\end{equation*}
$$
根據 IBM 教學抽象一點可以看作下方的形式:
$$
\begin{equation*}
\begin{split}
\lvert \text{State in Computational Basis} \rangle &\xrightarrow{QFT} \lvert \text{State in Fourier Basis}\rangle
\newline
\newline
\lvert x \rangle &\xrightarrow{QFT} \lvert y\rangle = \lvert \tilde{x}\rangle
\end{split}
\end{equation*}
$$
其中 $\omega$ 的定義一樣是來自單位根:
$N = 4$ 的情況，驗證6。CC BY 4.0 Huang, Po-Hsuan
可以注意到其實 QFT 是跟 IDFT 比較像，同樣是逆時針旋轉。這裡一開始可能會比較困惑，因為也有些人也是反過來看，兩種都有人用。常見的原因是使用慣例 (convention)。這邊可以參考 What are the input and output of QFT and IQFT, respectively? 中的討論。
$N = 8 = 2^n (n = 3)$。來源: link
QFT Matrix
$N = 4, (\omega = e^{\frac{2\pi i}{4}} = i)$ 時 QFT Matrix 如下 (跟 IDFT Matrix 一樣):
$$
M_{QFT} = {\frac{1}{\sqrt{N}}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & \omega & \omega^2 & \omega^3\\ 1 & \omega^2 & \omega^4 & \omega^6\\ 1 & \omega^3 & \omega^6 & \omega^9\end{pmatrix}} = {\frac{1}{\sqrt{4}}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & \omega & \omega^2 & \omega^3\\ 1 & \omega^2 & 1 & \omega^2\\ 1 & \omega^3 & \omega^2 & \omega\end{pmatrix}} = {\frac{1}{2}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & i & -1 & -i\\ 1 & -1 & 1 & -1\\ 1 & -i & -1 & i\end{pmatrix}}
$$
IQFT Matrix
在公式補上負號，就變成 Inverse Quantum Fourier Transform。
$N = 4, (\omega = e^{-\frac{2\pi i}{4}} = -i)$ 時 IQFT Matrix 如下 (跟 DFT Matrix 一樣):
$$
M_{IQFT} = {\frac{1}{\sqrt{N}}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & \omega & \omega^2 & \omega^3\\ 1 & \omega^2 & \omega^4 & \omega^6\\ 1 & \omega^3 & \omega^6 & \omega^9\end{pmatrix}} = {\frac{1}{\sqrt{4}}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & \omega & \omega^2 & \omega^3\\ 1 & \omega^2 & 1 & \omega^2\\ 1 & \omega^3 & \omega^2 & \omega\end{pmatrix}} = {\frac{1}{2}\small\begin{pmatrix}1 & 1 & 1 & 1\\ 1 & -i & -1 & i\\ 1 & -1 & 1 & -1\\ 1 & i & -1 & -i\end{pmatrix}}
$$
對我來說很特別的是，$M_{QFT}$ 跟 $M_{IQFT}$ 剛好互為共軛轉置7，所以 ${M_{QFT}}^\dagger = M_{IQFT}$ 且 $M_{QFT}\times M_{IQFT} = I$。
量子邏輯閘實作
Controlled rotation
一般用 CU1 實作對 Z 軸的旋轉，且控制控制位元為 $\lvert 1 \rangle$ 才會讓目標位元轉。$k$ 為控制位元與目標位元距離再加 1，例如第一個位元被第二個位元控制時，$k = (2-1)+1 = 2$ 則 $e^{i\lambda} = e^{\frac{\pi i}{2}} = i$，而也就是利用 $k$ 來實現不同角度的旋轉。
$$
CU1(\lambda) = {\small\begin{pmatrix}1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & e^{i\lambda}\end{pmatrix}}, \lambda = \frac{2\pi}{2^k}
$$
對 Z 軸旋轉不會影響 quantum state 出現的機率，只會變更 phase。來源: link
SWAP
最後輸出會用到交換 (SWAP)，讓輸出可以跟公式相同。Qiskit QFT 文件中8有提到，若是 SWAP 是整個電路的末端的話，可以將 SWAP 的工作移動到古典運算實作就可以。
QFT 邏輯閘
透過 IBM Quantum Composer 實作 $N = 4$ ($n = 2$) 的 QFT，輸入為 $[0 1 0 0]$ 為例。可透過 NOT 對兩個原為 $\lvert 0 \rangle$ 的值產生 $\lvert 1 \rangle$，因此有四種輸入可能:
- $\lvert 00\rangle$ ($\lvert 0 \rangle= \begin{pmatrix}1&0&0&0\end{pmatrix}^T$)
- $\lvert 01\rangle$ ($\lvert 1 \rangle= \begin{pmatrix}0&1&0&0\end{pmatrix}^T$): 下圖中使用的
- $\lvert 10\rangle$ ($\lvert 2 \rangle= \begin{pmatrix}0&0&1&0\end{pmatrix}^T$)
- $\lvert 11\rangle$ ($\lvert 3 \rangle= \begin{pmatrix}0&0&0&1\end{pmatrix}^T$)
2-bit QFT with $\lvert 01 \rangle$ as input
輸出說明:
q[0]: $\frac{1}{\sqrt{2}}(\lvert 0\rangle+e^{2\pi i\cdot\frac{1}{2}}\lvert 1\rangle) = \frac{1}{\sqrt{2}}(\lvert 0\rangle - \lvert 1\rangle)$q[1]: $\frac{1}{\sqrt{2}}(\lvert 0\rangle+e^{2\pi i\cdot\frac{1}{4}}\lvert 1\rangle) = \frac{1}{\sqrt{2}}(\lvert 0\rangle + i \lvert 1\rangle)$- $\Rightarrow \frac{1}{2}(\lvert00\rangle+i\lvert01\rangle-\lvert10\rangle-i\lvert11\rangle)$
OPENQASM 2.0;
include "qelib1.inc";
qreg q[2];
creg c[2];
x q[0];
barrier q[0],q[1];
h q[1];
cu1(pi/2) q[0],q[1];
h q[0];
swap q[0],q[1];
$\lvert 01 \rangle$ 為輸入的輸出模擬結果
驗證實作與公式是否對應
這邊為了減少手殘，所以使用 numpy 幫忙驗證。$N = 4 (n = 2)$，的公式驗證如下，CU1 中有一個 j 是因為 $e^{\frac{2\pi i}{4}}$ 是 $i$。 si 為輸入，可自行更改。這裡感謝 Small Quantum Fourier transforms | QuTech Academy - youtube 影片示範推導過程 (雖然影片中的邏輯閘位置有錯誤)。
import numpy as np
si = '0100'
Input = np.matrix(f'{si[0]};{si[1]};{si[2]};{si[3]}')
I = np.matrix('\
1 0; \
0 1')
H = np.matrix('\
1 1; \
1 -1')*(1/np.sqrt(2))
IH = np.kron(I, H) # tensor product
HI = np.kron(H, I) # tensor product
SWAP = np.matrix('\
1 0 0 0; \
0 0 1 0; \
0 1 0 0; \
0 0 0 1')
CU1 = np.matrix('\
1 0 0 0; \
0 1 0 0; \
0 0 1 0; \
0 0 0 1j')
QFT = SWAP.dot(IH).dot(CU1).dot(HI)
print(f'QFT:\n{QFT}\n')
print(f'QFT(|{si}>):\n{QFT.dot(Input)}')
驗證執行結果 (numpy 用 j 表示虛數單位 $i$):
QFT:
[[ 0.5+0.j 0.5+0.j 0.5+0.j 0.5+0.j ]
[ 0.5+0.j 0. +0.5j -0.5+0.j 0. -0.5j]
[ 0.5+0.j -0.5+0.j 0.5+0.j -0.5+0.j ]
[ 0.5+0.j 0. -0.5j -0.5+0.j 0. +0.5j]]
QFT(|0100>):
[[ 0.5+0.j ]
[ 0. +0.5j]
[-0.5+0.j ]
[ 0. -0.5j]]
複雜度比較
參考 QC Theory Lecture 20 Shor’s algorithm part III Quantum Fourier Transform QFT - youtube 內的說法得出下列複雜度比較 ($N = 2^n$)。最主要提升的原因是量子版本只需要 $\log N$ 位元即可運算 $N$ 種狀態。另外 Quantum Fourier Transform and Inverse QFT in Python - QISKit - youtube 提到 QFT 無法加速 DFT 在古典資料上的運算，原因是輸入輸出資料的表達。不過也有相關研究是在用量子版本取代古典版本，例如: QFFT9。
- DFT: $O(N^2) = O(2^{2n})$
- FFT: $O(N\log N) = O(n2^n)$
- QFT: $O(\log^2N) = O(n^2)$
結語
本篇利用從 DFT 的原理介紹到實際用小規模的範例完整示範 QFT 的原理及實作過程，應該可以稱得上是目前最完整的文章 (除非是一些學校超長的課程)，若有其他更完整資料再請提供給我參考。接下來 “有空的話” 會拓展到一般化的形式，到時候會有更多公式。另外，本文章也沒實作 IQFT，雖然是正負號換一下就好，不過因為篇幅已經過長所以就保留到下次。