分治法(5)：快速傅立叶变换(FFT)

本文介绍了多项式乘法的快速傅立叶变换算法，分为插值、乘法、恢复三个步骤，将乘法的时间复杂度降到O(nlogn)。

我们之前提到，要想将大整数乘法的时间复杂度降低到$\bm{O(n \log n)}$，可以采用快速傅立叶变换（Fast Fourier Transform，FFT）。这次，我们从一个类似的问题——多项式乘法入手，学习这个神奇的算法。

多项式乘法

基本问题

给定两个$d - 1$阶多项式$p(x)$和$q(x)$，我们要计算他们的乘积$r(x) = p(x)q(x)$。

每个$d - 1$阶多项式都可以表示为一个$d$维向量。设

• $p(x) = \sum_{i=0}^{d-1} a_i x^i \rightarrow (a_0, a_1, \ldots, a_{d - 1})$
• $q(x) = \sum_{i=0}^{d-1} b_i x^i \rightarrow (b_0, b_1, \ldots, b_{d - 1})$

那么，我们就是要计算$r(x) = \sum_{i=0}^{2d-2} c_i x^i$的系数$c_i = \sum_{k=0}^{i} a_k b_{i-k}$。

朴素算法的时间复杂度为$O(d^2)$，因为我们要计算系数的两两乘积再相加。所以，我们应该比$\bm{O(d^2)}$做得更好。

回想Karatsuba算法

如何改进？我们可以回想大整数乘法的Karatsuba算法。在这里，我们尝试类似的算法。

假设$d$是2的整数次幂，我们先将$p(x)$写成

$$p(x) = p_1(x) + p_2(x) \cdot x^{\frac{d}{2}}$$

其中

$$p_1(x) = a_0 + a_1 x + \cdots + a_{\frac{d}{2} - 1}$$

$$p_2(x) = a_{\frac{d}{2}} + a_{\frac{d}{2} + 1} x + \cdots + a_{d - 1} x^{\frac{d}{2} - 1}$$

这里$p_1(x)$相当于大整数乘法中的低位，$p_2(x)$相当于大整数乘法中的高位。

类似地，我们将$q(x)$写成

$$q(x) = q_1(x) + q_2(x) \cdot x^{\frac{d}{2}}$$

从而有

$$r = p_1 q_1 + (p_1 q_2 + p_2 q_1) x^{\frac{d}{2}} + p_2 q_2 x^d$$

让我们应用Karatsuba算法！我们要计算$p_1 q_1$、$p_2 q_2$、$p_1 q_2 + p_2 q_1$，但并不需要计算四次乘法，因为$p_1 q_2 + p_2 q_1$可以由下式得到。

$$p_1 q_2 + p_2 q_1 = (p_1 + p_2)(q_1 + q_2) - p_1 q_1 - p_2 q_2$$

这样，我们只需要计算三次乘法，即计算

• $p_1 q_1$
• $p_2 q_2$
• $(p_1 + p_2)(q_1 + q_2)$

于是，我们把一个规模为$d$的乘法问题，分解为了三个规模$\dfrac{d}{2}$的乘法问题，也就是

$$T(n) = 3 T\left(\frac{d}{2}\right) + O(d) \Rightarrow T(n) = O(n^{\log_2 3})$$

与大整数乘法的联系和区别

你一定发现了，多项式乘法和大整数乘法联系紧密。例如，

可以看出，大整数乘法似乎就是$x=10$的多项式乘法。但是，他们也有区别！

主要的区别是进位，多项式乘法中不存在大整数乘法中的进位，这会简化我们的问题。

事实上，我们今天使用多项式乘法来讲解FFT，主要就是为了避免进位问题。

插值定理

我们先做一些准备工作。为了表示一个多项式，我们可以它的系数表示，也可以用它曲线上的几个点表示。这就是多项式插值的想法。

那么，要几个不同的点才能确定一个多项式呢？我们从简单情况考虑，确定一条直线需要2个点，确定一条抛物线需要3个点。由此我们猜测，需要 $\bm{d}$ 个不同的点 $(x_0, p(x_0)), (x_1, p(x_1)), \ldots, (x_{d - 1}, p(x_{d - 1}))$ 来唯一确定一个 $\bm{d - 1}$ 维多项式。

• 插值定理：给定$d$个点$(x_0, y_0), (x_1, y_1), \ldots, (x_{d - 1}, y_{d - 1})$，其中对任意$i \neq j$都有$x_i \neq x_j$，那么存在一个唯一的不超过$d - 1$阶的多项式$p(x)$，使得对每个$i$都有$p(x_i) = y_i$。

我们可以用线性代数的知识证明这个定理。设$p(x) = \sum_{t=0}^{d-1} a_t x^t$，那么我们有$y_i = \sum_{t=0}^{d-1} a_t x_i^t$，即

$$\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{d-1} \end{bmatrix} = \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^{d} \\ 1 & x_1 & x_1^2 & \cdots & x_1^{d} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{d-1} & x_{d-1}^2 & \cdots & x_{d-1}^{d} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{d-1} \end{bmatrix}$$

我们想要证明，满足等式的$\bm{a} = (a_0, a_1, \ldots, a_{d - 1})$是唯一的。

注意到，中间的方阵是一个范德蒙德矩阵，其行列式等于$\prod_{0 \leq i < j \leq d-1} (x_j - x_i)$。又因为$x_i \neq x_j$，因此行列式不为零，即该矩阵可逆，从而满足等式的$\bm{a}$有唯一解。

FFT的整体框架

FFT的核心思路就是我们刚才说的多项式插值，其算法主要分为以下三步：

1. 插值：选取$2d-1$个不同的数$x_0, x_1, \ldots, x_{2d-2}$，计算

   • $p(x_0), p(x_1), \ldots, p(x_{2d-2})$
   • $q(x_0), q(x_1), \ldots, q(x_{2d-2})$

2. 乘法：对于每个$i = 0, 1, \ldots, 2d-2$，计算 $r(x_i) = p(x_i) q(x_i)$。这其实获得了$r(x)$的插值$(x_0, r(x_0)), (x_1, r(x_1)), \ldots, (x_{2d-2}, r(x_{2d-2}))$。

3. 恢复：从前一步获得的插值中恢复出$r(x) = \sum_{i=0}^{2d-2}$的系数$(c_0, c_1, \ldots, c_{2d-2})$。

整个FFT的过程相当于“用插值绕了一圈”，主要的目的是将乘法的次数从 $\bm{d^2}$ 次减少到 $\bm{2d-1}$ 次，因为我们只需要完成$2d-1$个整数乘法。这样，我们期望能降低整体的时间复杂度。

步骤一：插值

要想在$O(d \log d)$时间内做完多项式乘法，就要求每一步的时间复杂度不超过$O(d \log d)$。我们先来看第一步——插值。

最简单的想法

最直接的做法是，随意找$2d-1$个点$x_0, x_1, \ldots, x_{2d-2}$，再分别算出$p(x_0), p(x_1), \ldots, p(x_{2d-2})$ 和 $q(x_0), q(x_1), \ldots, q(x_{2d-2})$

假设我们能在$O(1)$时间内算出$x^d$，那么算出每个$p(x_i)$或$q(x_i)$需要$O(d)$的时间，而我们一共要算$4d-2$个。因此，总时间复杂度为$O(d^2)$。这显然是不让人满意的，因为我们要求每一步都要在$O(d \log d)$以内完成。

那么，我们可以用分治法做得更好吗？

引入分治法

我们首先做一些规定和假设。令$D = 2d-1$，并假设$D$是2的整数次幂（因为我们将会对$D$分治）。假设我们能在$O(1)$的时间内算出$x^d$。

回到正题，我们如何能加快插值的速度呢？目前的时间花销来源是，要计算$O(d)$个$p(x_i)$或$q(x_i)$，而每个需要花费$O(d)$时间。所以，我们有两种优化思路：

• 内部优化：用分治法更快地算出单个$p(x_i)$
• 协同优化：在算不同的$p(x_i)$时分治

内部优化的尝试

首先，我们尝试内部优化的思路。为了用分治法更快算出$p(x_i)$，我们将$p(x_i)$分成高位和低位，即

$$p(x_i) = p_1(x_i) + p_2(x_i) \cdot x^{\frac{D}{2}}$$

然后，我们分别递归地求出$p_1(x_i)$和$p_2(x_i)$。

相信你已经意识到，这样做不可能降低时间复杂度。验证一下，我们有

$$T(D) = 2T\left(\frac{D}{2}\right) + O(1) \Rightarrow T(D) = O(D)$$

这和直接计算没有区别，因为我们只是改变了计算顺序，并没有减少任何计算。事实上，将$x_i$代入多项式求值的过程并没有什么优化空间。

协同优化的尝试

我们尝试在计算不同的$p(x_i)$时使用分治。如果简单地将$D$个$p(x_i)$分成两半，再递归地计算每一部分，这不会和之前有任何区别，因为不同的$p(x_i)$之间毫无关系，我们无法减少任何计算。

• 不同的$p(x_i)$之间没有关联，我们要单独算$O(D)$个$p(x_i)$；
• 单个$p(x_i)$的各项之间没有关联，我们要单独算$D$个单项式的值。

这么看，我们对这个$O(d^2)$的算法根本无法优化！

单次求多值

让我们再仔细思考一下。计算每个$p(x_i)$需要$O(D)$的时间，这确实无法优化，因为多项式的系数并不由我们决定，我们只能逐项代入计算。

但是计算不同的$p(x_i)$呢？刚才我们说，不能优化是因为不同的 $\bm{p(x_i)}$ 之间没有关联。但其实，这$D$个点只是要求出$r(x)$的系数，选择什么点完全取决于我们！

那么，能不能通过选取特殊的 $\bm{x_0, x_1, \ldots, x_{D-1}}$ 来减少计算呢？例如，选取$x_1 = 1$和$x_2 = -1$，我们能不能同时计算$p(x_1)$和$p(x_2)$？

奇偶分治

以$x_1 = 1$和$x_2 = -1$为例，我们观察到

• $p(1) = a_0 + a_1 + a_2 + a_3 + \cdots$
• $p(-1) = a_0 - a_1 + a_2 - a_3 + \cdots$

这么看，我们可以分开计算奇数项和偶数项，再分别相加和相减，就能同时得到$p(x_i)$和$p(-x_i)$！也就是说，我们用$d$次运算解决了原来$2d$次的计算！

更一般地，我们有如下奇偶分治的思路：设

$$p(x) = p_e(x^2) + x \cdot p_o(x^2)$$

其中

$$p_e(x) = a_0 + a_2 x + a_4 x^2 + \cdots + a_{D-2} x^{\frac{D-2}{2}}$$

$$p_o(x) = a_1 + a_3 x + a_5 x^2 + \cdots + a_{D-1} x^{\frac{D-2}{2}}$$

我们选择互为相反数的$x_1$和$x_2$，那么有

$$p_e(x_1^2) = p_e(x_2^2) \quad \mathrm{and} \quad p_o(x_1^2) = p_o(x_2^2)$$

进而有

$$p(x_1) = p_e(x_1^2) + x_1 \cdot p_o(x_1^2)$$

$$p(x_2) = p_e(x_2^2) + x_2 \cdot p_o(x_2^2)$$

原本我们计算$p(x_1)$和$p(x_2)$，需要计算$p_e(x_1^2), p_o(x_1^2), p_e(x_2^2), p_o(x_2^2)$这4项，但现在，我们只需要计算2项。也就是说，我们原先将2个规模为$D$的问题，分解为4个规模$\dfrac{D}{2}$的问题，而现在，我们只需要分解为2个规模$\dfrac{D}{2}$的问题！

分治法思路

这样，我们就设计出了一种基本的分治算法。

1. 选择 $x_0, x_1, \ldots, x_{D-1}$，使得 $x_0 = -x_1, x_2 = -x_3, \ldots, x_{D-2} = -x_{D-1}$。

2. 划分：分别计算奇部分和偶部分。

   • $p_e(x_0^2), p_e(x_2^2), \ldots, p_e(x_{D-2}^2)$

   • $p_o(x_0^2), p_o(x_2^2), \ldots, p_o(x_{D-2}^2)$

3. 合并：计算 $p(x_i) = p_e(x_i^2) + x_i \cdot p_o(x_i^2)$。

这样，我们每次就将计算$D$个$D$阶$p(x)$的问题，转换为了计算$D$个$\dfrac{D}{2}$阶$p(x)$的问题，即分解为了2个“计算$\dfrac{D}{2}$个$\dfrac{D}{2}$阶$p(x)$的问题”。如果$T(D)$表示计算$D$个$D$阶的$p(x)$的时间代价，那么我们有

$$T(D) = 2T\left(\frac{D}{2}\right) + O(D) \Rightarrow T(D) = O(D \log D)$$

这似乎非常巧妙。但这个算法中存在着很大的问题。

引入复数

分治法中的问题

请注意，我们的算法是否能够递归调用。在第一轮中，我们取了互为相反数的整数，如$x_0 = -x_1, x_2 = -x_3$等。但在第二次调用中，要代入的整数变为了$x_0^2, x_2^2, x_4^2, \ldots$，全部是正数，无法再次采用相反数的性质进行递归。

因此，为了让我们的递归进行下去，就要保证每次要代入的整数都两两互为相反数。而想要让$x_0^2 = -x_2^2$，我们必须引入复数。

复数基础复习

复数有三种基本的表示形式

1. 直角坐标：$z = a + bi$，其中$a$为实部，$b$为虚部，$i = \sqrt{-1}$为虚数单位。

2. 极坐标：$z = r(\cos \theta + i \sin \theta)$，其中$r$是向量$(a, b)$的长度，$\theta$是向量$(a, b)$与实轴的夹角。

3. 指数形式：$z = r \cdot e^{\theta i} = r(\cos \theta + i \sin \theta)$，由欧拉公式可得。

接下来，我们都采用指数形式表示复数，并且默认$r = 1$，因为我们在FFT中只会用到$r = 1$（即单位圆上）的复数。

1. 复数的平方：$e^{\theta i}$的平方是$e^{2 \theta i}$，只需要把$e^{\theta i}$旋转角度$\theta$。

2. 复数的相反数：互为相反数的复数关于原点中心对称（图中两个红色向量），并且平方相等。

$$\left(e^{(\theta + \pi)i}\right)^2 = e^{2 \theta i} \cdot e^{2 \pi i} = e ^{2 \theta i} = \left(e^{\theta i}\right)^2$$

3. 复数的平方根：$e^{\theta i}$的平方根是$e^{\frac{\theta}{2} i}$和$e^{(\frac{\theta}{2} + \pi) i}$。

用复数实现递归

有了复数，我们如何解决刚才的问题？也就是，如何让平方后的数还两两互为相反数？

我们先考虑原先有4个数的情况，即怎么样取这4个数，才能让他们平方后得到的2个数仍然是相反数？不难想到，可以取单位圆与坐标轴的四个交点（或者说从实轴正方向起取单位圆的四等分点），即

$$x_0 = 1, \quad x_1 = -1, \quad x_2 = i, \quad x_3 = -i$$

这样平方后，有

$$x_0^2 = 1, \quad x_2^2 = -1$$

仍为相反数。

类似地，考虑 $D=8$ 的情况，我们取从实轴正方向起取单位圆的八等分点，即

$$\omega_0 = 1, \quad \omega_1 = e^{\frac{\pi}{4}i}, \quad \omega_2 = e^{\frac{\pi}{2}i}, \quad \omega_3 = e^{\frac{3\pi}{4}i}$$

$$\omega_4 = e^{\pi i}, \quad \omega_5 = e^{\frac{5\pi}{4}i}, \quad \omega_6 = e^{\frac{3\pi}{2}i}, \quad \omega_7 = e^{\frac{7\pi}{4}i}$$

可以验证，这8个数平方后得到的4个数两两互为相反数，这4个数平方后得到的2个数仍两两互为相反数。这样，我们分治法中的递归就可以进行了。

那么一般地，对于任意的$D$，我们如何找到起始的$D$个点呢？

和之前一样，从实轴正方向起取单位圆的 $\bm{D}$ 等分点。设 $\omega = e^{\frac{2\pi}{D}i}$，那么我们就取$\omega^0, \omega^1, \omega^2, \ldots, \omega^{D-1}$。可以看到我们只需要一个参数 $\bm{\omega}$ 来表示 $\bm{D}$ 个点。

那么在下一阶段，要计算的点就变为了 $\omega^0, \omega^2, \omega^4, \ldots, \omega^{D-2}$。也就是说，我们可以用 $\omega^2$ 来表示下一阶段的所有点。

插值整体算法

至此，我们得到了第一步——插值的整体算法。

给定一个$D-1$阶多项式 $p$ 和 $\omega = e^{\frac{2\pi}{D}i}$。

• 基本情况：如果 $\omega = 1$，此时问题规模$D = 1$，直接返回 $p(1)$。

• 由 $p(x)$，构造偶多项式 $p_o(x)$ 和奇多项式 $p_e(x)$。

• 对 $p_o(x)$ 和 $p_e(x)$，分别在点 $\omega^0, \omega^2, \ldots, \omega^{D-2}$ 上递归计算。

• 依次合并 $D$ 个点的运算结果，即

$$p(\omega^t) = p_e(\omega^{2t}) + \omega^t \cdot p_o(\omega^{2t})$$

• 返回 $D$ 个点的计算结果。

设 $T(D)$ 表示计算 FFT($p$, $\omega$) 的时间复杂度，其中 $p$ 的阶数为 $D - 1$。之前已经分析过，我们有

$$T(D) = 2T\left(\frac{D}{2}\right) + O(D) \Rightarrow T(D) = O(D \log D)$$

所以，插值这一步的时间复杂度 $T(D) = \bm{O(d \log d)}$。

步骤二：乘法

在这一步中，我们要对每个 $x_i \,(0 \le i \le 2d-2)$，计算 $r(x_i) = p(x_i) q(x_i)$。

当我们完成了插值，这一步是简单的，我们只需要做 $2d-1$ 次数值乘法即可，故乘法这一步的时间复杂度为 $\bm{O(d)}$。

步骤三：恢复

第三步，我们需要从 $2d-1$ 个 $(x_i, r(x_i))$ 中，恢复出 $r(x)$ 的系数 $(c_0, c_1, \ldots, c_{2d-1})$。

最简单的想法

恢复其实是插值的逆过程，本质上是要解如下的方程。

$$\begin{bmatrix} r(1) \\ r(\omega) \\ r(\omega^2) \\ \vdots \\ r(\omega^{D-1}) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega & \omega^2 & \cdots & \omega^{D-1} \\ 1 & \omega^2 & \omega^4 & \cdots & \omega^{2(D-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{D-1} & \omega^{2(D-1)} & \cdots & \omega^{(D-1)(D-1)} \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{D-1} \end{bmatrix}$$

我们已知 $(1, r(1)), (\omega, r(\omega)), (\omega^2, r(\omega^2)), \ldots, (\omega^{D-1}, r(\omega^{D-1}))$，其中 $\omega = e^{\frac{2\pi}{D}i}$，现在要求系数 $(c_0, c_1, \ldots, c_D)$。

令

$$A = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega & \omega^2 & \cdots & \omega^{D-1} \\ 1 & \omega^2 & \omega^4 & \cdots & \omega^{2(D-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{D-1} & \omega^{2(D-1)} & \cdots & \omega^{(D-1)(D-1)} \end{bmatrix}$$

那么，最简单的想法是

$$\begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{D-1} \end{bmatrix} = A^{-1} \begin{bmatrix} r(1) \\ r(\omega) \\ r(\omega^2) \\ \vdots \\ r(\omega^{D-1}) \end{bmatrix}$$

即直接解矩阵方程，或者说求 $A^{-1}$。

但是，这种做法的时间复杂度是 $O(D^3)$，这是我们不可接受的。

因此，我们需要利用 $A$ 矩阵的特殊性质，来更快地解出方程。

复矩阵基础复习

1. 复数的共轭：$z = a + bi$ 的共轭复数是 $\bar{z} = a - bi$。

2. 复向量的内积：两个复向量 $\bm{a}, \bm{b} \in \mathbb{C}^n$，它们的内积定义为

$$\left<\bm{a}, \bm{b}\right> = \sum_{j=1}^{n} \overline{a_j} \cdot b_j$$

3. 复向量的正交：如果 $\left<\bm{a}, \bm{b}\right> = 0$，那么 $\bm{a}$ 和 $\bm{b}$ 是正交的。

4. 复向量的标准正交：如果 $\left<\bm{a}, \bm{b}\right> = 0$，且 $\left<\bm{a}, \bm{a}\right> = \left<\bm{b}, \bm{b}\right> = 1$，那么 $\bm{a}$ 和 $\bm{b}$ 是标准正交的。

5. 酉矩阵（标准正交复矩阵）：如果复方阵 $A$ 的任意两列都标准正交，那么 $A$ 是酉矩阵（标准正交复矩阵）。

6. 复矩阵的共轭转置：复矩阵 $A$ 的共轭转置 $A^*$，定义为 $(A^*)_{i, j} = \overline{A_{j, i}}$。

举个例子，比如

$$A = \begin{pmatrix} a & c \\ b & d \end{pmatrix} \quad A^* = \begin{pmatrix} \overline{a} & \overline{b} \\ \overline{c} & \overline{d} \end{pmatrix}$$

7. 酉矩阵的重要性质：如果 $A$ 是酉矩阵，那么 $A$ 可逆，且 $\bm{A^{-1} = A^*}$。

接着上面的例子，如果 $A$ 是酉矩阵，那么

$$A A^* = \begin{pmatrix} a & c \\ b & d \end{pmatrix} \begin{pmatrix} \overline{a} & \overline{b} \\ \overline{c} & \overline{d} \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$$

因此我们发现，如果 $A$ 是一个酉矩阵，那么求它的逆矩阵会非常方便。

A矩阵的性质

那么我们会想去验证，恢复步骤中的 $\bm{A(\omega)}$ 矩阵是一个酉矩阵吗？

还记得，

$$A = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega & \omega^2 & \cdots & \omega^{D-1} \\ 1 & \omega^2 & \omega^4 & \cdots & \omega^{2(D-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{D-1} & \omega^{2(D-1)} & \cdots & \omega^{(D-1)(D-1)} \end{bmatrix}$$

注意到，$\overline{\omega} \omega = 1$，$\omega^D = e^{2\pi i} = 1$。

首先验证正交性。挑选两列，

$$\bm{c_i} = (1, \omega, \omega^2, \ldots, \omega^{D-1})$$

$$\bm{c_j} = (1, \omega^2, \omega^4, \ldots, \omega^{2(D-1)})$$

那么有

$$\left<\bm{c_i}, \bm{c_j}\right> = 1 + \omega + \omega^2 + \cdots + \omega^{D-1} = 0$$

一般地，对于任意不同的两列 $\bm{c_i}, \bm{c_j}$，

$$\begin{aligned} \langle \bm{c}_i, \bm{c}_j \rangle &= \sum_{k=1}^{D} \overline{\omega^{(k-1)(i-1)}} \, \omega^{(k-1)(j-1)} \\ &= \sum_{k=1}^{D} \omega^{(k-1)(j-i)} = \frac{1 - \omega^{(j-i)D}}{1 - \omega^{j-i}} = 0 \end{aligned}$$

再验证标准性。但是，$A$ 矩阵明显不满足标准性（即每一列和自己的内积为1），因为第1列为全1。

事实上，对任意一列 $\bm{c_i}$，我们有

$$\langle \bm{c_i}, \bm{c_i} \rangle = \sum_{k=1}^{D} \overline{\omega^{(k-1)(i-1)}} \, \omega^{(k-1)(i-1)} = D$$

因此，$A$ 矩阵并不是酉矩阵！但是，我们可以通过缩放来使得它标准正交。

我们已经计算过，$A(\omega)$ 矩阵不同的两列内积为0，相同的两列内积为$D$。可以猜测，$\bm{\dfrac{1}{\sqrt{D}}A(\omega)}$ 应该是标准正交的！以下是对于这个命题的证明，可以跳过。

▶

命题及证明

命题：矩阵$\dfrac{1}{\sqrt{D}} A(\omega)$ 对于 $\omega = e^{\frac{2\pi i}{D}}$ 是标准正交的

证明：设 $\bm{c}_i, \bm{c}_j$ 是 $\dfrac{1}{\sqrt{D}} A(\omega)$ 的任意两列，那么有

$$\langle \bm{c}_i, \bm{c}_j \rangle = \sum_{k=1}^D \frac{1}{D} \overline{\omega^{(k-1)(i-1)}} \omega^{(k-1)(j-1)} = \frac{1}{D} \sum_{k=1}^D \omega^{(k-1)(j-i)}$$

• 如果 $i = j$，那么有

$$\langle \bm{c}_i, \bm{c}_j \rangle = \frac{1}{D} \sum_{k=1}^D \omega^0 = 1$$

• 如果 $i \neq j$，那么有

$$\langle \bm{c}_i, \bm{c}_j \rangle = \frac{1}{D} \sum_{k=1}^D \omega^{(k-1)(j-i)} = \frac{1}{D} \frac{1 - \omega^{(j-i)D}}{1 - \omega^{j-i}} = 0$$

因此， $\dfrac{1}{\sqrt{D}} A(\omega)$ 是标准正交的。

接下来，我们要求 $\bm{A(\omega)^{-1}}$。简单来说，由于

$$A(\omega) A(\omega)^* = D I$$

因此

$$A(\omega)^{-1} = \frac{1}{D} \cdot A(\omega)^{*} = \frac{1}{D} \cdot A(\omega^{-1})$$

由于 $A$ 对称，转置后不变。又 $\overline{\omega} = \omega^{-1}$，故 $A$ 的共轭转置 $A(\omega)^* = A(\omega^{-1})$。

当然，我们还可以严谨地证明，此部分可以跳过。

▶

求A的逆矩阵的严谨证明

$$\begin{aligned}A(\omega)^{-1} =& \left(\sqrt{D} \cdot \frac{1}{\sqrt{D}} \cdot A(\omega) \right)^{-1} = \sqrt{D} \left(\frac{1}{\sqrt{D}} \cdot A(\omega) \right)^{-1} \\=& \sqrt{D} \left(\frac{1}{\sqrt{D}} \cdot A(\omega) \right)^* = \frac{1}{D} A(\omega)^*\end{aligned}$$

因此，

$$(A(\omega)^{-1})_{i, j} = \frac{1}{D} \overline{(A(\omega)^{-1})_{j, i}} = \frac{1}{D} \cdot \omega^{-(i-1)(j-1)} = \frac{1}{D}(\omega^{-1})^{(i-1)(j-1)}$$

也就是说，

$$A(\omega)^{-1} = \frac{1}{D} \cdot A(\omega^{-1})$$

因此，

$$A^{-1} = \frac{1}{D} \cdot \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega^{-1} & \omega^{-2} & \cdots & \omega^{-(D-1)} \\ 1 & \omega^{-2} & \omega^{-4} & \cdots & \omega^{-2(D-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{-(D-1)} & \omega^{-2(D-1)} & \cdots & \omega^{-(D-1)(D-1)} \end{bmatrix}$$

恢复整体算法

现在，我们就得到了

$$\begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{D-1} \end{bmatrix} = \frac{1}{D} \cdot \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega^{-1} & \omega^{-2} & \cdots & \omega^{-(D-1)} \\ 1 & \omega^{-2} & \omega^{-4} & \cdots & \omega^{-2(D-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{-(D-1)} & \omega^{-2(D-1)} & \cdots & \omega^{-(D-1)(D-1)} \end{bmatrix} \begin{bmatrix} r(1) \\ r(\omega) \\ r(\omega^2) \\ \vdots \\ r(\omega^{D-1}) \end{bmatrix}$$

然而，直接计算这个矩阵乘法仍然需要 $O(D^2)$ 的时间，我们还是不能接受。

该如何改进呢？我们发现，现在要做的事情和第一步很类似。我们都是在求一个多项式在 $\bm{D}$ 个点处的值，只不过现在的多项式是

$$s(x) = r(1) + r(\omega) \cdot x + r(\omega^2) \cdot x^2 + \cdots + r(\omega^{D - 1}) \cdot x^{D - 1}$$

而现在的 $D$ 个数分别是 $\omega^0, \omega^{-1}, \omega^{-2}, \ldots, \omega^{-(D-1)}$。

这可以使用第一步的分治法吗？

现在和第一步的区别是：在第一步中，我们选择了 $D$ 个很好的数 $\omega^0, \omega^1, \ldots, \omega^{D-1}$，而在第三步中，我们无法选择数字，现在的数字就是 $\omega^0, \omega^{-1}, \omega^{-2}, \ldots, \omega^{-(D-1)}$。

那么，这 $D$ 个数是否也很好呢？也就是说，它们是否有平方后两两互为相反数的良好性质，使得我们的递归可以很好进行下去呢？

直观上看，是有的。因为我们在第一步在单位圆上取 $D$ 等分点，采取了顺时针方向；而这次，采用的是逆时针方向，当然能取到完全相同的点。也就是说，$(\omega^0, \omega^1, \ldots, \omega^{D-1})$ 和 $(\omega^0, \omega^{-1}, \omega^{-2}, \ldots, \omega^{-(D-1)})$ 只有旋转的方向不同。

所以，我们可以直接调用第一步中的函数 FFT$(s, \omega^{-1})$！从而，恢复步骤的时间复杂度也为 $\bm{O(d \log d)}$。

整体算法与时间复杂度分析

至此，快速傅立叶变换的三个步骤已经全部设计完成！

让我们把它们组装起来：

• 插值（第1～4步）：计算 $D$ 和 $\omega$，再计算出 $\omega^i$ 对应的函数值 $p(\omega^i)$ 和 $q(\omega^i)$。
• 乘法（第5步）：对每个点，计算 $r(\omega^i) = p(\omega^i)q(\omega^i)$。
• 恢复（第6～9步）：利用 $r(\omega_i)$ 结合 FFT 函数，解出 $r(x)$ 的系数。

总的时间复杂度为：

$$O(d \log d) + O(d) + O(d \log d) = \bm{O(d \log d)}$$

注：本文中所有图片均来自张宇昊老师的课程PPT。