DC: FFT Part 2

PDF of Eric’s handwritten notes are here.

Complex Numbers Review

We can represent a complex number $z = a + b i$ as a point $(a,b)$ in the complex plane or in polar coordinates $(r,\theta)$ where:
$z = (a,b) = (r \cos{\theta}, r \sin{\theta})$
Euler’s Identity gives a convenient form for a complex number by noting that $\cos{\theta} + i \sin{\theta} = e^{i\theta}$ and thus $z = r e^{i\theta}$ .

A benefit of the representation in polar coordinates is that they are convenient for multiplication: $(r_1, \theta_{1}) . (r_2, \theta_{2}) = (r_1 r_2, \theta_{1} + \theta_{2})$
A few convenient basic facts:
Multiplication by $-1$ : $-(r, \theta) = (1, \pi) . (r, \theta) = (r, \theta + \pi)$
Exponentiation: ${(r, \theta)}^n = (r^n, n \theta)$
Hence, ${(1, \theta)}^n = (1, n \theta)$

We will use the n^th complex roots of unity:
A complex number that is a solution to the equation $z^n = 1$ is referred to as an $n^{\text{th}}$ root of unity. If we work in polar coordinates, then we can obtain a very nice classification of all $n^{\text{th}}$ roots of unity.

We are looking for $z$ such that $z^n = (r,\theta)^n = (r^n, n \theta) = (1,0)$ . This implies that $r = 1$ and $\theta = \frac{2 \pi j}{n}$ for some $j \in \{0,1,\dots, n-1\}$ . We denote:
$\omega_n = (1, \frac{2\pi}{n}) = e^{i\frac{2\pi}{n}}$ ,
and then all $n$ of the $n^{th}$ roots of unity are
$\omega_n^0, \omega_n^1, \omega_n^2, \dots, \omega_n^{n-1}$ .

Key Properties:

They satisfy the $\pm$ property when $n$ is even i.e.
$\omega_n^0 = - \omega_n^{\frac{n}{2}}$
$\omega_n^1 = - \omega_n^{\frac{n}{2}+1}$
$\vdots$
$\omega_n^{\frac{n}{2}-1} = - \omega_n^{n-1}$
To see why, observe that $\omega_n^j = (1, \frac{2\pi j}{n})$ and $\omega_n^{\frac{n}{2} + j} = (1, \frac{2\pi }{n} ( {\frac{n}{2} + j})) = (1, \pi + \frac{2\pi j}{n}) = - (1, \frac{2\pi j}{n}) = -\omega_n^j$
The squares of the $n^{th}$ roots of unity are the ${\frac{n}{2}}^{th}$ roots of unity when $n$ is a power of $2$ . To see why, observe that ${(\omega_n^j)}^2 = {(1, \frac{2\pi j}{n})}^2 = (1, \frac{4\pi j }{n}) = (1, \frac{2\pi j }{\frac{n}{2}}) = \omega_{\frac{n}{2}}^j$

FFT Algorithm

Given a polynomial of degree at most $n-1$ , $A(x) = a_0 + a_1 x + \dots + a_{n-1}x^{n-1}$ , we want to choose $n$ points $x_0, x_1, \dots, x_{n-1}$ and evaluate $A(x)$ at these points. (We can view $A(x)$ as a polynomial of degree $\leq 2n-1$ in order to get it at $2n$ points.) We require the points $x_0, x_1, \dots, x_{n-1}$ to satisfy the $\pm$ property, i.e., $x_j = - x_{\frac{n}{2} + j}$ for $j \in \{0,1,\dots, \frac{n}{2}-1\}$ . Hence we will choose the points to be the $n^{th}$ roots of unity, i.e., $x_j = \omega_n^j$ for $j \in \{0, 1, \dots, n - 1\}$ .

Define the two polynomials based on the even and odd terms of $A(x)$ . Let $A_{even}(y) = a_0 + a_2 y + \dots + a_{n-2} y^{\frac{n}{2} - 1}$ and $A_{odd}(y) = a_1 + a_3 y + \dots + a_{n-1} y^{\frac{n}{2} - 1}$ . We will recursively compute $A_{even}(y), A_{odd}(y)$ at $\frac{n}{2}$ points $y_0 = x_0^2 = x_{\frac{n}{2}}^2, y_1 = x_1^2 = x_{\frac{n}{2}+1}^2, \dots, y_{\frac{n}{2}-1} = x_{\frac{n}{2}-1}^2 = x_{n-1}^2$ . Property (2) above says that these $n/2$ points $y_0,\dots,y_{\frac{n}{2}-1}$ are the ${\frac{n}{2}}^{th}$ roots of unity. Then we obtain $A(x)$ at $n$ points by using the relation $A(x) = A_{even}(x^2) + x A_{odd}(x^2)$ .

Notice that our original problem was to evaluate a polynomial $A(x)$ of degree $\leq n-1$ at the $n^{th}$ roots of unity. We then have two recursive subproblems, namely evaluating the polynomials $A_{even}(y)$ and $A_{odd}(y)$ of degree $\leq \frac{n}{2}-1$ at the ${\frac{n}{2}}^{th}$ roots of unity. Thus we have 2 subproblems which are of the same form and of half the size. We are now ready to present the FFT algorithm.

Later (to do Inverse FFT) we will evaluate a polynomial again at the $n^{th}$ roots of unity, but in the reverse order $\omega_n^0, \omega_n^{n-1}, \omega_n^{n-2}, \dots, \omega_n^1$ . To make the procedure modular enough to fit both situations we take as an input parameter $\omega$ where $\omega$ is any $n^{th}$ root of unity. Our FFT procedure will follow by setting $\omega=\omega_n$ and we will later see that our inverse FFT procedure will follow by setting $\omega=\omega_n^{n-1}$ .

Input: Coefficients of $A(x)$ : $\vec{a} = (a_0, a_1, \dots, a_{n-1})$ and a $n^{th}$ root of unity $\omega$

Output: $A(\omega^0), A(\omega^1), \dots, A(\omega^{n-1})$

$FFT(\vec{a}, \omega)$ :
If $n = 1$ , Return $(A(1))$
Let $\vec{a}_{even} = (a_0, a_2, \dots, a_{n-2}), \vec{a}_{odd} = (a_1, a_3, \dots, a_{n-1})$
Call $FFT(\vec{a}_{even}, \omega^2)$ to get $A_{even}(\omega^0), A_{even}(\omega^2), \dots, A_{even}(\omega^{2(\frac{n}{2}-1)})$
Call $FFT(\vec{a}_{odd}, \omega^2)$ to get $A_{odd}(\omega^0), A_{odd}(\omega^2), \dots, A_{odd}(\omega^{2(\frac{n}{2}-1)})$
For $j = 0$ to $\frac{n}{2} -1$ :
$\indent A(\omega^j) = A_{even}(\omega^{2j}) + \omega^j A_{odd}(\omega^{2j})$
$\indent A(\omega^{\frac{n}{2} + j}) = A_{even}(\omega^{2j}) + \omega^{\frac{n}{2}+j} A_{odd}(\omega^{2j})$
Return $(A(\omega^0), A(\omega^1), \dots, A(\omega^{n-1}))$

This can be written more concisely as:

$FFT(\vec{a}, \omega)$ :
If $n = 1$ , Return $(A(1))$
Let $\vec{a}_{even} = (a_0, a_2, \dots, a_{n-2}), \vec{a}_{odd} = (a_1, a_3, \dots, a_{n-1})$
$(s_0, s_1, \dots, s_{\frac{n}{2} - 1}) = FFT(\vec{a}_{even}, \omega^2)$
$(t_0, t_1, \dots, t_{\frac{n}{2} - 1}) = FFT(\vec{a}_{odd}, \omega^2)$
For $j = 0$ to $\frac{n}{2} -1$ :
$r_j = s_j + \omega^j t_j$
$r_{\frac{n}{2} + j} = s_j - \omega^{j} t_j$
Return $(r_0, r_1, \dots, r_{n-1})$

Running Time:

To evaluate $A(x)$ at $n$ points, we evaluate $A_{even}(x)$ and $A_{odd}(x)$ at $\frac{n}{2}$ points, and combine them in $O(n)$ time. Thus, we have the recurrence $T(n) = 2T(\frac{n}{2}) + O(n)$ which solves to $T(n) = O(n \log{} n)$ .

Inverse FFT

Recall from the previous lecture that to multiply $A(x)$ and $B(x)$ , we evaluate these polynomials at $2n$ points and multiply to obtain $C(x)$ at $2n$ points. Finally, we need to convert back from point representation to coefficient representation to get the coefficients of $C(x)$ . In this section, we will see how to perform this key interpolation step.

For a polynomial $A(x)$ , when we apply $FFT(\vec{a}, \omega_n)$ , we get $A(\omega^0_n), A(\omega^1_n), \dots, A(\omega^{n-1}_n)$ .

In general, for points $x_0, x_1, \dots, x_{n-1}$ , we have:

$\begin{bmatrix} A(x_0)\\ A(x_1)\\ .\\ .\\ .\\ A(x_{n-1}) \end{bmatrix} = \begin{bmatrix} 1 & x_0 & x_0^2 & \dots & x_0^{n-1}\\ 1 & x_1 & x_1^2 & \dots & x_1^{n-1}\\ .\\ .\\ .\\ 1 & x_{n-1} & x_{n-1}^2 & \dots & x_{n-1}^{n-1}\\ \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ .\\ .\\ .\\ a_{n-1} \end{bmatrix}$

In particular, for points $x_j = w_n^j$ , we have:

$\begin{bmatrix} A(1)\\ A(\omega_n)\\ A(\omega_n^2)\\ .\\ .\\ A(\omega_n^{n-1}) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \dots & 1\\ 1 & \omega_n & \omega_n^2 & \dots & \omega_n^{n-1}\\ 1 & \omega_n^2 & \omega_n^4 & \dots & \omega_n^{2(n-1)}\\ . & . & . & \dots & . \\ . & . & . & \dots & . \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \dots & \omega_n^{(n-1)(n-1)}\\ \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ .\\ .\\ .\\ a_{n-1} \end{bmatrix}$

Let $\vec{A} = \begin{bmatrix} A(1)\\ A(\omega_n)\\ A(\omega_n^2)\\ .\\ .\\ A(\omega_n^{n-1}) \end{bmatrix}$ , $M_n(\omega_n) = \begin{bmatrix} 1 & 1 & 1 & \dots & 1\\ 1 & \omega_n & \omega_n^2 & \dots & \omega_n^{n-1}\\ 1 & \omega_n^2 & \omega_n^4 & \dots & \omega_n^{2(n-1)}\\ . & . & . & \dots & . \\ . & . & . & \dots & . \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \dots & \omega_n^{(n-1)(n-1)}\\ \end{bmatrix}$ , $\vec{a} = \begin{bmatrix} a_0\\ a_1\\ .\\ .\\ .\\ a_{n-1} \end{bmatrix}$

Thus, $\vec{A} = M_n(\omega_n) \vec{a}$

In order to compute the vector $\vec{a}$ from vector $\vec{A}$ , we need an efficient way to compute ${M_n(\omega_n)}^{-1}$ .

Lemma: ${M_n(\omega_n)}^{-1} = \frac{1}{n} M_n(\omega_n^{-1})$

Observe that $\omega_n^{-1} = \omega_n^{n-1}$ since $\omega_n . \omega_n^{n-1} = \omega_n^n = 1$ as $\omega_n$ is an $n^{th}$ root of unity.

So, $\vec{a} = {M_n(\omega_n)}^{-1} \vec{A} = \frac{1}{n} M_n(\omega_n^{n-1}) \vec{A}$

Then, to compute ${M_n(\omega_n)}^{-1}$ , we will apply FFT.

In particular, $FFT(\vec{A}, \omega_n^{n-1}) = M_n(\omega_n^{-1}) \vec{A} = n {M_n(\omega_n)}^{-1} \vec{A} = n \vec{a}$

Thus, $\vec{a} = \frac{1}{n} FFT(\vec{A}, \omega_n^{n-1})$

That gives our algorithm for inverse FFT. Given the values of $A(x)$ at the $n^{\text{th}}$ roots of unity, view those values as the coefficients for a new polynomial. Then run FFT for this new polynomial with parameter $\omega=\omega_n^{n-1})$ , scale the output by a factor $n$ and this gives the coefficients for $A(x)$ .

It remains to prove the lemma. We will need the following basic fact.

Claim: For any $n^{\text{th}}$ root of unity $\omega \neq 1$ , then $1 + \omega + \omega^2 + \ldots + \omega^{n-1} = 0$ .

Proof: For any number $z$ notice that

$(z-1)(1+z+z^2+\dots+z^{n-1})=(z+z^2+\dots+z^n) - (1+z+\dots + z^{n-1}) = z^n-1$ .

Let $z=\omega$ . Since $\omega$ is a $n^{\text{th}}$ root of unity we have that $\omega^n-1=0$ . Hence:

$(\omega -1)(1+\omega + \omega^2 + \ldots + \omega^{n-1})=0$ .

Thus, either $(\omega -1)=0$ and/or $(1+\omega + \omega^2 + \dots + \omega^{n-1})=0$ . We assumed that $\omega \neq 1$ and thus we have that $(1+\omega + \omega^2 + \dots + \omega^{n-1})=0$ , which proves the claim.

Proof of Key Lemma:

To prove the lemma, we need to show that

$\frac{1}{n} M_n(\omega_n) M_n(\omega_n^{-1}) = I$

where $I$ is the $n \times n$ identity matrix (this matrix has 1’s on the diagonal and 0 for all other entries).

Let’s first look at the diagonal entries of $M_n(\omega_n) M_n(\omega_n^{-1})$ . Thus, consider entry $(k,k)$ of $M_n(\omega_n)\times M_n(\omega_n^{-1})$ for some $k=0,\ldots,n-1$ :

$\begin{aligned} &= (1,\omega_n^k, \omega_n^{2k}, \ldots, \omega_n^{(n-1)k}) \cdot (1,\omega_n^{-k},\omega_n^{-2k},\ldots,\omega_n^{-(n-1)k}) \\ &= (1 + (\omega_n\omega_n^{-1})^k + (\omega_n\omega_n^{-1})^{2k} + \dots + (\omega_n\omega_n^{-1})^{(n-1)k} )\\ &=1+1+\ldots+1 \\ &=n \end{aligned}$

Now, consider entry $(k,j)$ where $k\neq j$ of $M_n(\omega_n)\times M_n(\omega_n^{-1})$ :

$\begin{aligned} &=(1,\omega_n^k,\omega_n^{2k},\ldots ,\omega_n^{(n-1)k})\cdot (1,\omega_n^{-j},\omega_n^{-2j},\dots ,\omega_n^{-(n-1)j}) \\ &=1+\omega_n^{k-j} + \omega_n^{2(k-j)} + \dots + \omega_n^{(n-1)(k-j)} \end{aligned}$

Let $\omega = \omega_n^{k-j}$ and note that $w \neq 1$ since $k \neq j$ and $0 \le k,j \le n-1$ . Then,

$\begin{aligned} &=1+\omega + \omega^2 + \ldots + \omega^{n-1} \\ &=0 \qquad \text{ from our earlier claim.} \end{aligned}$ .

This proves the lemma.

Algorithms

Lecture videos+notes

DC: FFT Part 2