Convolution

Consider two polynomials $F(x) = \sum_kf_kx^k$, $G(x) = \sum_kg_kx^k$ and their product $H(x) = (F \cdot G)(x) = \sum_kh_kx^k$, where

$$ h_k = \sum_{i + j = k}f_ig_j. \tag*{} $$

Toeplitz Matrix

We can describe convolution using the language of linear algebra.

Embed the Toeplitz matrix into a circulant matrix $C$:

$$ C_{ij} = t_{i - j} = \begin{bmatrix} t_0 & t_{-1} & t_{-2} & \cdots & \cdots & t_{-(n - 1)} \\ t_1 & t_0 & t_{-1} & \ddots & & \vdots \\ t_2 & t_1 & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & t_{-1} & t_{-2} \\ \vdots & & \ddots & t_1 & t_0 & t_{-1} \\ t_{n - 1} & \cdots & \cdots & t_2 & t_1 & t_0 \end{bmatrix}, \tag*{} $$

where $t_k = t_{k \bmod n}$, so that $C$ is circulant. Then consider $\mathbf y = \mathbf t * \mathbf x$:

$$ y_i = \sum_kt_kx_{i - k} = \sum_jt_{i - j}x_j = \sum_jC_{ij}x_j \tag*{}. $$

Therefore $\mathbf y = \mathbf t * \mathbf x = C\mathbf x$.

We want to find a way to compute $C\mathbf x$ efficiently, so we want to diagonalize $C$.

Diagonalization

Define cyclic shift operator $S$:

$$ S = \begin{bmatrix} 0 & 1 & 0 & \cdots \\ 0 & 0 & 1 & \cdots \\ \vdots & \vdots & \vdots & \ddots \\ 1 & 0 & 0 & \cdots \end{bmatrix}. \tag*{} $$

Then

$$ C = t_0I + t_1S + t_2S^2 + \cdots + t_{n - 1}S^{n - 1}. \tag*{} $$

Since each circulant matrix is a polynomial in $S$, diagonalizing $S$ simultaneously diagonalizes all circulant matrices. We want to find a basis $\{\mathbf v_k\}$ such that

$$ S\mathbf v_k = \lambda_k\mathbf v_k, \tag*{} $$

then

$$ C\mathbf v_k = \left(\sum_jt_j\lambda_k^j\right)\mathbf v_k. \tag*{} $$

The Fourier basis consists of eigenvectors of $S$:

$$ \mathbf v_k = (1, \omega^k_n, \omega^{2k}_n, \ldots, \omega^{(n - 1)k}_n)^{\intercal}, \tag*{} $$

Where $\omega_n = e^{2\pi i/n}$.

Discrete Fourier Transform

Define operator

$$ C_{\mathbf a}(\mathbf b) := \mathbf a * \mathbf b. \tag*{} $$

By assembling the Fourier basis vectors into the discrete Fourier transform matrix, we obtain a matrix that diagonalizes $S$.

$$ F = [\omega_n^{jk}]_{0 \leq j, k < n} = \begin{bmatrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0 \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n - 1}\\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n - 1)} \\ \omega_n^0 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n - 1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \omega_n^0 & \omega_n^{n - 1} & \omega_n^{2(n - 1)} & \omega_n^{3(n - 1)} & \cdots & \omega_n^{(n - 1)^2} \end{bmatrix}. \tag*{} $$

We have

$$ F^{-1}SF = \operatorname{diag}(\omega_n^0, \omega_n^1, \ldots, \omega_n^{n - 1}), \tag*{} $$

so

$$ F^{-1}C_{\mathbf a}F = \sum_ja_j(F^{-1}SF)^j = \operatorname{diag}\left(\sum_ja_j\omega_n^{jk}\right)_{k = 0}^{n - 1}. \tag*{} $$

Let

$$ \hat a_k = \sum_{j = 0}^{n - 1}a_j\omega_n^{jk}. \tag*{} $$

Let $A(x) = \sum_{j = 0}^{n - 1}a_jx^j$.

$$ F^{-1}C_{\mathbf a}F = \operatorname{diag}(\hat{\mathbf a}) = \operatorname{diag}\bigl(A(\omega^k_n)\bigr)_{k = 0}^{n - 1}. \tag*{} $$

Let $\mathcal F\{\mathbf a\} = F\mathbf a = \hat{\mathbf a}$.

Back to our original problem: polynomial multiplication.

Since $(A \cdot B)(x) = A(x) \cdot B(x)$, by multiplying each element of one vector by the corresponding element of the other vector:

$$ \mathcal F\{\mathbf a * \mathbf b\} = \mathcal F\{\mathbf a\} \cdot \mathcal F\{\mathbf b\}. \tag*{} $$

Finally, applying the inverse DFT, we obtain:

$$ A \cdot B = \mathcal F^{-1}\bigl\{\mathcal F\{\mathbf a\} \cdot \mathcal F\{\mathbf b\}\bigr\}. \tag*{} $$

At this point, we already know how to use DFT to compute polynomial multiplication. The remaining problem is to apply DFT efficiently.

Fast Fourier Transform

Assume $n$ is a power of $2$.

$$ F_n = \underbrace{ \begin{bmatrix} I & D \\ I & -D \end{bmatrix} }_{\text{butterfly}} \cdot \underbrace{ \begin{bmatrix} F_{n / 2} & 0 \\ 0 & F_{n / 2} \end{bmatrix}}_{\text{recursion}} \cdot \underbrace{P}_{\text{rearrange}}, \tag*{} $$

where

$$ D = \operatorname{diag}(1, \omega_n^1, \omega_n^2, \ldots, \omega_n^{n/2 - 1}), \tag*{} $$

and $P$ is the permutation matrix that maps

$$ (a_0, a_1, \ldots, a_{n - 1}) \mapsto (a_0, a_2, a_4, \ldots, a_1, a_3, a_5, \ldots). \tag*{} $$

Let

$$ A(x) = \sum_{k = 0}^{n - 1}a_kx^k. \tag*{} $$

Split it into two polynomials:

$$ A_0(x) = \sum_{k = 0}^{n/2 - 1}a_{2k}x^k = a_0 + a_2x + a_4x^2 + \cdots + a_{n - 2}x^{n/2 - 1}, \tag*{} $$

$$ A_1(x) = \sum_{k = 0}^{n/2 - 1}a_{2k + 1}x^k = a_1 + a_3x + a_5x^2 + \cdots + a_{n - 1}x^{n/2 - 1}, \tag*{} $$

so that

$$ A(x) = A_0(x^2) + xA_1(x^2). \tag*{} $$

Pick $x = \omega_n^k$ in $A(x)$:

$$ \begin{aligned} A(\omega_n^k) &= A_0\bigl((\omega_n^k)^2\bigr) + \omega_n^kA_1\bigl((\omega_n^k)^2\bigr) \\ &= A_0(\omega_{n/2}^k) + \omega_n^kA_1(\omega_{n/2}^k). \end{aligned} \tag{$\dagger$} $$

And

$$ \begin{aligned} A(\omega_n^{k + n/2}) = A(-\omega_n^k) &= A_0\bigl((-\omega_n^k)^2\bigr) - \omega_n^kA_1\bigl((-\omega_n^k)^2\bigr) \\ &= A_0(\omega_{n/2}^k) - \omega_n^kA_1(\omega_{n/2}^k). \end{aligned} \tag{$\ddagger$} $$

Therefore, after applying DFT to $A_0$ and $A_1$, use $\dagger$ to compute $A(\omega_n^k)$ for $0 \leq k < n/2$ and use $\ddagger$ for $n/2 \leq k < n$.

We have

$$ T_{\text{DFT}}(n) = 2T_{\text{DFT}}(n/2) + O(n). \tag*{} $$

Through master theorem, $T_{\text{DFT}}(n) = O(n \log n)$.

Inverse DFT

At this point, we have

$$ \hat{\mathbf a} = \mathcal F\{\mathbf a\} = F\mathbf a. \tag*{} $$

Hence

$$ \mathbf a = \mathcal F^{-1}\{\hat{\mathbf a}\} = F^{-1}\hat{\mathbf a}. \tag*{} $$

A quick check can verify that the inverse of the matrix $F$ has the following form:

$$ F^{-1} = \frac1n[\omega_n^{-jk}]_{0 \leq j, k < n} = \frac1n \begin{bmatrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0 \\ \omega_n^0 & \omega_n^{-1} & \omega_n^{-2} & \omega_n^{-3} & \cdots & \omega_n^{-(n - 1)}\\ \omega_n^0 & \omega_n^{-2} & \omega_n^{-4} & \omega_n^{-6} & \cdots & \omega_n^{-2(n - 1)} \\ \omega_n^0 & \omega_n^{-3} & \omega_n^{-6} & \omega_n^{-9} & \cdots & \omega_n^{-3(n - 1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \omega_n^0 & \omega_n^{-(n - 1)} & \omega_n^{-2(n - 1)} & \omega_n^{-3(n - 1)} & \cdots & \omega_n^{-(n - 1)^2} \end{bmatrix}. \tag*{} $$

Note that these problems are almost the same as DFT, so the coefficients $a_k$ can be found by the same divide and conquer algorithm, as well as the direct FFT, only instead of $\omega_n^k$, we have to use $\omega_n^{-k}$, and at the end we need to divide the resulting coefficients by $n$.