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 nth 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:

  1. 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}
    \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
  2. 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})=0which 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.