MIT 6.1220 — Lectures 22–23

Fast Fourier Transform

Definition (Convolution). Given sequences {aj}j=0m\{a_j\}_{j = 0}^m and {bk}k=0n\{b_k\}_{k = 0}^n of numbers, their convolution a∗ba \ast b is:

{cℓ}ℓ=0n+m     given by     cℓ=(a∗b)ℓ:=∑j=0ℓajbℓ−j.\{c_\ell\}_{\ell = 0}^{n + m} ~~~~ \text{ given by } ~~~~ c_{\ell} = (a \ast b)_{\ell} := \sum_{j = 0}^{\ell} a_j b_{\ell - j}.

This is just polynomial multiplication, or the PMF of the sum of two independent discrete random variables.

Problem. (FFT) Given {aj}j=0m\{a_j\}_{j = 0}^m and {bk}k=0n\{b_k\}_{k = 0}^n, compute their convolution a∗ba \ast b in O((m+n)log⁡(m+n))O((m + n) \log (m + n)) time.

Remark. The naive brute-force algorithm runs in O(mn)O(mn) time.

Let's black-box a solution to FFT for now. How might it be useful?

Problem. (String Matching) Fix a string s=(s0,…,sn)∈{0,1}n+1s = (s_0, \dots, s_n) \in \{0, 1\}^{n + 1} and a pattern p=(p0,…,pm)∈{0,1}m+1p = (p_0, \dots, p_m) \in \{0, 1\}^{m + 1}. Find the set HH of all indices kk for which sk+i=pis_{k + i} = p_i for all i=0,…,mi = 0, \dots, m.

Solution. (String Matching) Switch to a ±1\pm 1 encoding, reverse pp, and compute s∗ps \ast p. This is O(nlog⁡n)O(n \log n).   ■\blacksquare

Now let's un-black-box the solution to FFT. How is FFT implemented in O((m+n)log⁡(m+n))O((m + n) \log (m + n)) time?

Solution. (FFT) It's equivalent to compute p(x)⋅q(x)p(x) \cdot q(x) for polynomials satisfying deg⁡(p)=m\deg (p) = m and deg⁡(q)=n\deg(q) = n.

The driving principle behind our solution is the following.

Fact. Any degree-nn polynomial is uniquely determined by its evaluations on any n+1n + 1 distinct values.

A good choice of values to evaluate our polynomials on, it turns out, is the full set of (2t)th(2^t)^{\text{th}} roots of unity.

Claim (Discrete Fourier Transform). Let nn be a power of two, and let ω:=e2πin\omega := e^{\frac{2\pi i}{n}}. For any polynomial pp with degree (n−1)(n - 1) or lower, we can compute {p(ωk)}k=0n−1\{p(\omega^{k}) \}_{k = 0}^{n - 1} in O(nlog⁡n)O(n \log n) time.

Proof: Divide and conquer. Say the runtime is T(n)T(n).

  • Denote peven(x)p_{\text{even}}(x) and podd(x)p_{\text{odd}}(x) so that p(x)=peven(x2)+x⋅podd(x2)p(x) = p_{\text{even}}(x^2) + x \cdot p_{\text{odd}}(x^2).

  • In 2T(n2)2T\left(\frac{n}{2}\right) time, compute {peven(ω2k)}k=0n/2−1\{p_{\text{even}}(\omega^{2k})\}_{k = 0}^{n / 2 - 1} and {podd(ω2k)}k=0n/2−1\{p_{\text{odd}}(\omega^{2k})\}_{k = 0}^{n / 2 - 1}.

  • In O(n)O(n) time, we can then compute:

    p(ωk)={peven(ω2k)+ωk⋅podd(ω2k)if k<n2.peven(ω2k−n)−ωk−n/2⋅podd(ω2k−n)if k≥n2.p(\omega^k) = \begin{cases} p_{\text{even}}\left(\omega^{2k}\right) + \omega^k \cdot p_{\text{odd}}\left(\omega^{2k}\right) & \text{if $k < \frac{n}{2}$.} \\ p_{\text{even}}\left(\omega^{2k - n}\right) - \omega^{k - n/2} \cdot p_{\text{odd}}\left(\omega^{2k - n}\right) & \text{if $k \geq \frac{n}{2}$.} \end{cases}

This yields T(n)=2T(n2)+O(n)T(n) = 2T\left(\frac{n}{2}\right) + O(n), which solves to T(n)=O(nlog⁡n)T(n) = O(n \log n).   □\square

Thus, we select ω=e2πiN\omega = e^{\frac{2\pi i}{N}} for some power-of-two N>m+nN > m + n, then compute {p(ωk)⋅q(ωk)}k=0N−1\{p(\omega^k) \cdot q(\omega^k)\}_{k = 0}^{N - 1} in O((m+n)log⁡(m+n))O((m + n) \log (m + n)) time. It remains to invert the DFT process above to extract p(x)⋅q(x)p(x) \cdot q(x) from here.

Claim (Inverse Discrete Fourier Transform). Given the values {p(ωk)}k=0n−1\{p(\omega^k)\}_{k = 0}^{n - 1} for some polynomial pp with degree (n−1)(n - 1) or lower, we can compute p(x)p(x) in O(nlog⁡n)O(n \log n) time.

Proof: The point is that the Inverse DFT is no harder than the DFT.

Say p(x)=a0+a1x+a2x2+⋯+an−1xn−1p(x) = a_0 + a_1x + a_2x^2 + \dots + a_{n - 1}x^{n - 1}. Then the DFT computes a matrix product:

[p(ω0)p(ω1)p(ω2)⋮p(ωn−1)]=[111⋯11ωω2⋯ωn−11ω2ω4⋯ω2(n−1)⋮⋮⋮⋱⋮1ωn−1ω2(n−1)⋯ω(n−1)(n−1)]⏟W[a0a1a2⋮an−1].\begin{bmatrix} p(\omega^0) \\ p(\omega^1) \\ p(\omega^2) \\ \vdots \\ p(\omega^{n - 1}) \end{bmatrix} = \underbrace{\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega & \omega^2 & \cdots & \omega^{n - 1} \\ 1 & \omega^2 & \omega^4 & \cdots & \omega^{2(n - 1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{n - 1} & \omega^{2(n - 1)} & \cdots & \omega^{(n - 1)(n - 1)} \end{bmatrix}}_{W} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n - 1} \end{bmatrix}.

We effectively wish to find an efficient way to multiply the LHS by W−1W^{-1}. The trick is that WW satisfies a really nice property:

(W‾W)s,t=∑k=0n−1W‾s,k⋅Wk,t=∑k=0n−1ωs⋅k‾⋅ωk⋅t=∑k=0n−1ωk⋅(t−s)={n if s=t.0 otherwise.\left(\overline{W} W\right)_{s, t} = \sum_{k = 0}^{n - 1} \overline{W}_{s, k} \cdot W_{k, t} = \sum_{k = 0}^{n - 1} \overline{\omega^{s \cdot k}} \cdot \omega^{k \cdot t} = \sum_{k = 0}^{n - 1} \omega^{k \cdot (t - s)} = \begin{cases} n & \text{ if } s = t. \\ 0 & \text{ otherwise.} \end{cases}

Therefore, W‾W=nI\overline{W}W = nI, so all we need to do is perform the DFT (with ω‾\overline{\omega} in place of ω\omega) on the LHS and divide by nn.   □\square

And that completes the implementation of FFT!   ■\blacksquare

Remark. In 18.701 terms, the DFT matrix WW satisfies the property that 1nW\frac{1}{\sqrt{n}}W is unitary.