FFT – Maths to Share

The discrete Fourier transform (DFT) is a mapping from $\mathbb{C}^n$ to $\mathbb{C}^n$ . Given a vector $u \in \mathbb{C}^n$ , the discrete Fourier transform of $u$ is another vector $U \in \mathbb{C}^n$ given by,

$\displaystyle{U_k = \sum_{j=0}^{n-1} \exp\left(-2 \pi i kj/m\right)u_j}.$

If we let $\omega_n = \exp\left(-2 \pi i /m\right)$ , then $U$ can be written more compactly as

$\displaystyle{U_k = \sum_{j=0}^{n-1} \omega^{kj}_n u_j }.$

The map from $u$ to $U$ is linear and can therefor be represented as multiplication by an $n \times n$ matrix $D$ . Looking at the equation above, we can see that if we set $D_{kj} = \omega^{kj}_n$ and use $0$ -indexing, then,

$U = Du.$

This suggests that calculating the vector $U = Du$ would require roughly $n^2$ operations. However, the famous fast Fourier transform shows that $U$ can be computed in roughly $n\log_2(n)$ operations. There are many derivations of the fast Fourier transform but I wanted to focus on one which is based on factorizing the matrix $D$ into sparse matrices.

Matrix factorizations

A factorization of a matrix $A$ is a way of writing $A$ as a product of matrices $A = A_L\cdots A_2A_1$ where typically $L$ is $2$ or $3$ . Matrix factorizations can tell us information about the matrix $A$ . For example, the singular value decomposition tells us the spectral norm and singular vectors $A$ (and lots more). Other matrix factorizations can help us compute things related to the matrix $A$ . The QR decomposition helps us solve least-squares problems and the LU factorization let us invert $A$ in (almost) the same amount of time it takes to solve the equation $Ax = b$ .

Each of the factorizations above applied to arbitrary matrices. Every possible matrix has an singular value decomposition and a QR factorization. Likewise all square matrices have an LU factorization. The DFT matrix $D$ has a factorization which most matrices other matrices do not. This factorization lets us write $D = D_L\cdots D_2D_1$ where each matrix $D_l$ is sparse.

Sparse matrices

A matrix $A$ is sparse if many of the entries of $A$ are zero. There isn’t a precise definition of “many”. A good heuristic is that if $A$ has size $m$ by $n$ , then $A$ is sparse if the number of non-zero entries is of the order $m$ or $n$ . In contrast, a dense matrix is one where the number of non-zero entries is of the same order of $mn$ .

The structure of sparse matrices makes them especially nice to work with. They take up a lot less memory and they are quicker to compute with. If $A \in \mathbb{C}^{m \times n}$ has only $k$ non-zero elements, then for any vector $x \in \mathbb{C}^n$ , the matrix-vector product $Ax$ can be computed in order $k$ operations. If $k$ is much less than $mn$ , then this can be a huge speed up.

This shows why a sparse factorization of the DFT matrix $D$ will give a fast way of calculating $Du$ . Specifically, suppose that we have $D = D_L \cdots D_2D_1$ and each matrix $D_l$ has just $k_l$ non-zero entries. Then, for any vector $u \in \mathbb{C}^n$ , we can compute $U= Du=D_L \cdots D_2 D_1 u$ in stages. First we multiply $u$ by $D_1$ to produce $u^{(1)}$ . The output $u^{(1)}$ is then multiplied by $D_2$ to produce $u^{(2)}$ . We continue in this way until we finally multiply $u^{(L-1)}$ by $D_L$ to produce $u = D_L u^{(L-1)}$ . The total computation cost is of the order $k_1 + k_2 + \cdots + k_L$ which can be drastically less than $n^2$ .

The DFT factorization

We’ll now show that $D$ does indeed have a sparse matrix factorization. To avoid some of the technicalities, we’ll assume that $n$ is a power of $2$ . That is, $n = 2^m$ for some $m$ . Remember that the discrete Fourier transform of $u$ is the vector $U$ given by

$\displaystyle{U_k = \sum_{j=0}^{n-1} \omega^{jk}_n u_j},$

where $\omega_n = \exp(-2\pi i/n)$ . The idea behind the fast Fourier transform is to split the above sum into two terms. One where we add up the even values of $j$ and one where we add up the odd values of $j$ . Specifically,

$\displaystyle{U_k = \sum_{j=0, \text{j even }}^{n-1} \omega_n^{jk} u_j + \sum_{j = 0, j \text{ odd}}^{n-1} \omega_n^{jk}u_j = \sum_{j=0}^{n/2 - 1} \omega^{2jk}_n u_{2j} + \sum_{j=0}^{n/2-1} \omega^{(2j+1)k}_n u_{2j+1}}$

We’ll now look at each of these sums separately. In the first sum, we have terms of the form $\omega^{2 jk}_n$ , these can by rewritten as follows,

$\omega^{2 jk}_n = \exp(-2 \pi i /n)^{2 jk} = \exp(-2 \pi i 2/n)^{jk} = \exp(-2 \pi i /(n/2))^{jk} = \omega_{n/2}^{jk}$ .

We can also use this result to rewrite the terms $\omega^{(2j+1)k}_n$ ,

$\omega^{(2j+1)k}_n = \omega_n^k\omega^{2jk}_n = \omega_n^k \omega^{jk}_{n/2}.$

If we plug both of these back into the expression for $U = Du$ , we get

$\displaystyle{U_k = \sum_{j=0}^{n/2 - 1} \omega_{n/2}^{jk} u_{2j} + \sum_{j=0}^{n/2-1} \omega_n^k \omega_{n/2}^{jk} u_{2j+1} = \sum_{j=0}^{n/2-1} \omega_{n/2}^{jk} u_{2j} + \omega_n^k \sum_{j=0}^{n/2-1}\omega_{n/2}^{jk} u_{2j+1}. }$

This shows that the Fourier transform of the $n$ -dimensional vector $u$ can be done by computing and combining the Fourier transforms of the $n/2$ -dimensional vectors. One of these vectors contain the even indexed elements of $u$ and the others contain the odd indexed elements of $u$ .

There is one more observation that we need. If $k$ is between $n/2$ and $n-1$ , then we can simplify the above expression for $U_k$ . Specifically, if we set $k'=k-n/2$ , then

$\omega_{n/2}^{jk} = \omega_{n/2}^{j(n/2+k')} = \omega_{n/2}^{jn/2}\omega_{n/2}^{jk'} = \exp(-2\pi i (n/2)/(n/2))^j \omega_{n/2}^{jk'} = \omega_{n/2}^{jk'}.$

If we plug this into the above formula for $U_k$ , then we get

$\displaystyle{U_k = \sum_{j=0}^{n/2-1} \omega_{n/2}^{j(k-n/2)} u_{2j} + \omega_n^k \sum_{j=0}^{n/2-1}\omega_{n/2}^{j(k-n/2)} u_{2j+1}}$ for $k = n/2, n/2+1,\ldots, n-1$ .

This reduces the computation since there are some computations we only have to do once but can use twice. The algebra we did above shows that we can compute the discrete Fourier transform of $u \in \mathbb{C}^n$ in the following step.

Form the vectors $u^{\text{even}} = (u_0,u_2,\ldots,u_{n-2}) \in \mathbb{C}^{n/2}$ and $u^{\text{odd}} = (u_1,u_3,\ldots,u_{n-1}) \in \mathbb{C}^{n/2}$ .
Calculate the discrete Fourier transforms of $u^{\text{even}}$ and $u^{\text{odd}}$ . This gives two vector $U^{\text{even}}, U^{\text{odd}} \in \mathbb{C}^{n/2}$ .
Compute $U \in \mathbb{C}^n$ by,

$\displaystyle{U_k = \begin{cases} U^{\text{even}}_k + \omega_n^k U^{\text{odd}}_k & \text{for } k = 0,1,\ldots, n/2-1,\\ U^{\text{even}}_{k-n/2}+\omega_n^k U^{\text{odd}}_{k-n/2} &\text{for } k=n/2,n/2+1,\ldots,n-1 \end{cases}}$

The steps above break down computing $U=Du$ into three simpler steps. Each of these simpler steps are also linear maps and can therefore be written as multiplication by a matrix. If we call these matrices $D_1,D_2$ and $D_3$ , then

The matrix $D_1$ maps $u$ to $\begin{bmatrix} u^{\text{even}} \\ u^{\text{odd}}\end{bmatrix}$ . That is

$D_1 = \begin{bmatrix} 1&0&0&0&\cdots&0&0\\0&0&1&0&\cdots &0&0 \\\vdots& \vdots &\vdots & \vdots & \ddots & \vdots & \vdots \\0&0&0&0&\cdots&1&0 \\ 0&1&0&0&\cdots&0&0 \\ 0&0&0&1&\cdots&0&0 \\ \vdots&\vdots &\vdots & \vdots & \ddots & \vdots & \vdots \\ 0&0&0&0&\cdots&0&1\end{bmatrix} \in \mathbb{C}^{n \times n}.$

The matrix $D_2$ is block-diagonal and each block is the $n/2$ -dimensional discrete Fourier transformation matrix. That is

$D_2 = \begin{bmatrix}D^{(n/2)} & 0\\ 0 & D^{(n/2)}\end{bmatrix}\in \mathbb{C}^{n \times n},$

where $D^{(n/2)}$ is given by $D^{(n/2)}_{jk} = \omega_{n/2}^{jk}$ .

The matrix $D_3$ combines two $n/2$ -vectors according to the formula in step 3 above,

$D_3 = \begin{bmatrix} 1&0&\cdots & 0 & \omega_n^0 & 0 & \cdots & 0\\ 0&1&\cdots&0 &0& \omega_n^1&\cdots & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots& \vdots &\ddots&\vdots\\ 0&0&\cdots & 1 & 0 & 0 & \cdots & \omega_n^{n/2-1}\\1&0&\cdots & 0 & \omega_n^{n/2} & 0 & \cdots & 0\\ 0&1&\cdots&0 &0& \omega_n^{n/2+1}&\cdots & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots& \vdots &\ddots&\vdots\\ 0&0&\cdots & 1 & 0 & 0 & \cdots & \omega_n^{n-1}\end{bmatrix}\in \mathbb{C}^{n \times n}$

The three-step procedure above for calculating $U=Du$ is equivalent to iteratively multiplying by $D_1$ , then $D_2$ and the finally by $D_3$ . Another way of saying this is that $D=D_3D_2D_1$ is a matrix factorization!

From the definitions, we can see that each of the matrices $D_1,D_2$ and $D_3$ have a lot of zeros. Specifically, $D_1$ has just $n$ non-zero entries (one in each row), $D_2$ has $2(n/2)^2$ non-zero entries (two blocks of $n/2 \times n/2$ non-zeros) and $D_3$ has $2n$ non-zero entries (two in each row). This means that the product $U=D_3D_2D_1u$ can be computed in roughly

$n+2(n/2)^2 + 2n = 3n+2(n/2)^2,$

operations. This is far away from the $n\log_2(n)$ operations promised at the start of this blog post. However, the matrix $D_2$ contains two copies of the $n/2 \times n/2$ discrete Fourier transform matrix. This means we could apply the exact same argument to the $n/2 \times n/2$ matrix $D^{(n/2)}$ . This lets us replace $(n/2)^2$ with $\frac{3}{2}n + 2(n/4)^2$ . This brings the number of operations down to