Q1. The Fast Fourier Transform¶
In this exercise, we derive the "most important numerical algorithm of our lifetime" (G. Strang): the Fast Fourier Transform (FFT). The FFT is a fast algorithm to compute the Discrete Fourier Transform (DFT). The DFT of a vector $x \in \mathbb{C}^n$ is the vector $y\in \mathbb{C}^n$ defined by:
$$ y_k = \sum_{j=0}^{n-1} x_j \exp(-i2\pi k j / n ) \ . $$(Here, vectors are 0-indexed). We restrict our attention to the case where the size $n$ is a power of $2$, although the ideas extend to the general case. The DFT is the result of the matrix-vector multiply $Fx$ with the matrix
$$ F = \begin{bmatrix} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega & \omega^2 & \omega^3 & \cdots & \omega^{n-1} \\ 1 & \omega^2 & \omega^4 & \omega^6 & \cdots & \omega^{2(n-1)} \\ 1 & \omega^3 & \omega^6 & \omega^9 & \cdots & \omega^{2(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots \\ 1 & \omega^{n-1} & \omega^{2(n-1)} & \omega^{3(n-1)} & \cdots & \omega^{(n-1)(n-1)} \\ \end{bmatrix} $$where $\omega = \exp(-i 2 \pi /n) $ (note that $\omega$ is a function of n).
To illustrate the FFT, we take $n = 8$ (this is the presentation in Golub & Van Loan).
$$ F_8 = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & \omega & \omega^2 & \omega^3 & \omega^4 & \omega^5 & \omega^6 & \omega^7 \\ 1 & \omega^2 & \omega^4 & \omega^6 & \omega^1 & \omega^2 & \omega^4 & \omega^6 \\ 1 & \omega^3 & \omega^6 & \omega & \omega^4 & \omega^7 & \omega^2 & \omega^5 \\ 1 & \omega^4 & 1 & \omega^4 & 1 & \omega^4 & 1 & \omega^4 \\ 1 & \omega^5 & \omega^2 & \omega^7 & \omega^4 & \omega & \omega^6 & \omega^3 \\ 1 & \omega^6 & \omega^4 & \omega^2 & 1 & \omega^6 & \omega^4 & \omega^2 \\ 1 & \omega^7 & \omega^6 & \omega^5 & \omega^4 & \omega^3 & \omega^2 & \omega \\ \end{bmatrix} $$where we have used the fact that $$\omega^{n+k} = \omega^{n}\omega^{k} = \exp(-i 2 \pi / n)^{n} \omega^k = \exp(-i 2 \pi n/n)\omega^k = 1 \omega^k $$
Next, we reorder the columns so that the odd-indexed columns come first (they are now ordered as [1,3,5,7,2,4,6,8]), which is equivalent to multiplying $F_8$ by the permutation matrix on the right:
$$ P_8 = \begin{bmatrix} 1 & & & & & & & \\ & & & &1 & & & \\ & 1& & & & & & \\ & & & & &1 & & \\ & &1 & & & & & \\ & & & & & & 1 & \\ & & &1 & & & & \\ & & & & & & &1 \\ \end{bmatrix} $$Thus, the matrix $F_8P_8$ is:
$$F_8P_8 = \left[ \begin{array}{cccc|cccc} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & \omega^2 & \omega^4 & \omega^6 & \omega & \omega^3 & \omega^5 & \omega^7 \\ 1 & \omega^4 & 1 & \omega^4 & \omega^2 & \omega^6 & \omega^2 & \omega^6 \\ 1 & \omega^6 & \omega^4 & \omega^2 & \omega^4 & \omega^7 & \omega^2 & \omega^5 \\ \hline 1 & 1 & 1 & 1 & -1 & -1 & -1 & -1 \\ 1 & \omega^2 & \omega^4 & \omega^6 & -\omega & -\omega^3 & -\omega^5 & -\omega^7 \\ 1 & \omega^4 & 1 & \omega^4 & \omega^2 & -\omega^6 & -\omega^2 & -\omega^6 \\ 1 & \omega^6 & \omega^4 & \omega^2 & -\omega^4 & -\omega^7 & -\omega^2 & -\omega^5 \\ \end{array} \right] = \left[ \begin{array}{c|c} F_4 & D_4 F_4 \\ \hline F_4 & - D_4 F_4 \\ \end{array} \right] $$where $D_4$ is a diagonal matrix with entries $\text{diag}(D_4) = [1, \omega_8, \omega_8^2, \omega_8^3] $ (the subscript on $\omega_8$ is to emphasize that it is $\omega_8 = \exp(-i 2 \pi /8) $) and $F_4$ is the DFT matrix with $n= 4$:
$$ F_{4} = \left[ \begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & \omega^2 & \omega^4 & \omega^6 \\ 1 & \omega^4 & 1 & \omega^4 \\ 1 & \omega^6 & \omega^4 & \omega^2 \\ \end{array} \right] $$That is, up to a re-ordering of columns and a scaling matrix $D_4$, $F_8$ is made out of four copies of $F_4$. The matrix can be further factorized as:
$$ F_8P_8 = \left[ \begin{array}{c|c} F_4 & D_4 F_4 \\ \hline F_4 & - D_4 F_4 \\ \end{array} \right] = \left[ \begin{array}{c|c} I & D_4 \\ \hline I & - D_4 \\ \end{array} \right] \begin{bmatrix} F_4 & 0 \\ 0 & F_4 \end{bmatrix} $$
It follows that:
$$F_8 x = F_8P_8 P_8^{-1} x = \left[ \begin{array}{c|c} F_4 & D_4 F_4 \\ \hline F_4 & - D_4 F_4 \\ \end{array} \right] P_8^{T} x = \left[ \begin{array}{c|c} F_4 & D_4 F_4 \\ \hline F_4 & - D_4 F_4 \\ \end{array} \right] \begin{bmatrix} x_{0:2:7} \\ x_{1:2:7} \end{bmatrix} $$where $x_{0:2:7} = [x_0, x_2, x_4, x_6]^T$, $ x_{1:2:7} = [x_1, x_3, x_5, x_7]^T$. Thus
$$ F_8 x= \left[ \begin{array}{c|c} F_4 & D_4 F_4 \\ \hline F_4 & - D_4 F_4 \\ \end{array} \right] \begin{bmatrix} x_{1:2:8} \\ x_{2:2:8} \end{bmatrix} = \left[ \begin{array}{c} (F_4 x_{0:2:7}) + D_4 ( F_4 x_{1:2:7} ) \\ (F_4 x_{0:2:7}) - D_4 ( F_4 x_{1:2:7} ) \end{array} \right] $$The argument generalizes to any n (divisible by 2), and thus we have:
$$ F_n P_n = \left[ \begin{array}{c|c} I & D_{n/2} \\ \hline I & - D_{n/2} \\ \end{array} \right] \begin{bmatrix} F_{n/2} & 0 \\ 0 & F_{n/2} \end{bmatrix} $$where the diagonal entries of $D_n$ are : $[ 1, \omega_n, \omega_n^2, \omega_n^3 \dots \omega_n^{n/2-1}]$. The idea of the FFT is to use this formula recursively to divide-and-conquer the computation.
At the top of the tree, we seek a DFT of size $n$ which we "divide-and-conquer" into two DFTs of size $n/2$, we then scale and add those vectors which takes $8n$ operations. At next level, where we compute two $F_{n/2}$, we do $8(n/2)$ work to scale/add vectors, thus $8 (n/2) \times 2 = 8n$ total work. Recursively, we see that we do $8n$ work at each level of the tree. After $\log_2(n) = d$ rounds of divide and conquer, we hit the base case of $n =1$, for which $F_1x = x$. Summing up the work of each level, we thus do $ ~ 8n d = 8n \log_2 n$ work.
Write a recursive algorithm for the DFT assuming $n = 2^d$.