【学习笔记】快速傅里叶变换(FFT)

发布于 2019-03-12  93 次阅读


0.0 前言

0.1 关于本篇文章

本篇文章写作主要有两个目的,一是本人想对这一块的知识进行系统的梳理,还有就是将《算法导论》上的知识写成更适合OIer学习的格式(这本书由于各种原因是真的难啃)

为了方便阅读理解,在看本篇文章的时候我们规定一些标记:
1. 重点难点强调将以加粗形式标注
2. 长篇幅引用将以如下形式标注:

我是引用内容

  1. (括号斜体)中内容为原书没有,本人为了方便理解,在不改变意思的情况下加上的内容,如果有误欢迎联系我
  2. 循环变量将优先使用j表示,i要作为虚数单位。
  3. 楷体字为定理或引理

0.2 FFT是用来干什么的?

将多项式乘法从 $O(n^2)$ 加速到 $O(nlogn)$ 。

0.3 前置知识

有以下的只是基础会让你的学习变得更加轻松

  1. 复数
  2. 多项式基本概念
  3. 矩阵基本运算

0.4 基本原理框架

为了更方便地理解下面每个概念的意义,先给出算法实现的大体框架:

1.0 多项式

1.1 多项式基本概念

为了让我们对多项式的基本概念达成共识,这部分简单内容还是稍微提下:
1. 多项式表示:$A(x)=\displaystyle\sum_{j=0}^{n-1}a_jx^j$
2. 多项式的系数:其中, \(\{a_0,a_1,\cdots,a_{n-1} \}\) 被称为系数所有系数都属于复数集合C
3. 多项式的次数: $A(x)$ 的具有非零系数的最高次的次数,记为 $degree(A)$
4. 次数界:任何一个严格大于次数的整数
这就是一般的多项式,如果拿来相乘的话,时间复杂度为$O(n^2)$

1.2 多项式的表示

1.2.1 系数表达

多项式 $A(x)=\displaystyle\sum_{j=0}^{n-1}a_jx^j$ 的系数表达为一个由系数组成的向量 $a=(a_0,a_1,\cdots,a_{n-1})$
这也是我们平时最常用的表达方式。

1.2.2 点值表达

一个次数界为 $n$ 的多项式 $A(x)$的点值表达就是一个由 $n$ 个点值对所组成的集合

\(\left\{ (x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1}) \right\}\)

使 $k=0,\ 1,\ 2, \cdots, n-1$ 所有的 $x_k$ 各不相同,
$$y_k=A(x_k)$$
一个多项式可以有很多不同的点值表达。

在求点值表达的时候我们可以使用秦九韶算法,但这样的做法的时间复杂度任然是 $O(n^2)$ 之后我们将会讨论如何优化这一过程。
求值计算的逆运算称为插值

1.2.3 插值多项式的唯一性(可跳过)

对于任意 $n$ 个点值对组成的集合 \(\{(x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\}\) ,其中所有的 $x_k$ 都不同;那么存在唯一的次数界为 $n$ 的多项式 $A(x)$ ,满足 $y_k=A(x_k),k=0,1,\cdots,n-1$ 。

证明 证明主要是根据某个矩阵存在逆矩阵。
$$\begin{bmatrix}1&x_0&x_0^2&\cdots&x_0^{n-1}\1&x_1&x_1^2&\cdots&x_1^{n-1}\\vdots&\vdots&\vdots&\cdots&\vdots\1&x_{n-1}&x_{n-1}^2&\cdots&x_{n-1}^{n-1}\end{bmatrix}\begin{bmatrix}a_0\a_1\\vdots\a_{n-1}\end{bmatrix}=\begin{bmatrix}y_0\y_1\\vdots\y_{n-1}\end{bmatrix}$$
左边矩阵的行列式值为:
$$\displaystyle\prod_{0\le j< k\le n-1}(x_k-x_j)$$

下略(需要用到矩阵相关知识,这一部分内容超出本人理解范围,而且不影响我们的算法)

2.0 单位复数根

2.1 单位复数根的定义

单位负数根就是我们之前说的选的特殊值,它具有一些神奇的性质以至于我们可以分治的方法,来用它求点值表达。
我们定义n次单位复数根是满足 $\omega_n=1$ 的复数 $\omega$ 。
为什么这样定义呢?我们现在需要用到欧拉公式
$$e^{iu}=cos(u)+isin(u)$$

如下图( $\omega$ 表示的是 $\omega_8^1$ 其他的都是它的次幂),说明 $n$ 个单位复数根均匀地分布在以复平面的原点为圆心的单位半径的圆周
$$\omega_n=e^{\frac{2\pi i}{n}}$$

细心的读者可能已经发现,在这 $n$ 个单位复数根内,似乎它们总是成对存在(比如 $\omega_8^8$ 和 $\omega_8^4$ ),这是不是意味着我们能分治地去处理求值问题呢?

2.2.0 三大引理

在这一小节中,将介绍单位复数根的三大引理,这也是降低算法时间复杂度的关键。

2.2.1 消去引理

对于任何整数 $n\ge 0, k\ge 0, d>0$,总有
$$\omega_{dn}^{dk}=\omega_{n}^{k}$$

证明 首先我们可以形象地理解这个结论,假设我把单位圆分成 $n$ 份,从 $\omega_{n}^{0}$ 开始,取第 $k$ 份,和我把单位圆分成 $dn$ 份,从 $\omega_{dn}^{0}$ 开始,取第 $dk$ 份是不是取到的边是一样的? 如果你喜欢严谨一点的话,下面的推导式子也不难理解。
$$\omega_{dn}^{dk}=(e^{\frac{2\pi i}{dn}})^{dk}=(e^{\frac{2\pi i}{n}})^k=\omega_{n}^{k}$$

2.2.2 折半引理

如果 $n>0$ (且)为偶数,那么 $n$ 个 $n$ 次单位复数根的平方的集合就是 $n/2$ 个 $n/2$ 次单位复数根的集合

从上面的文字来理解其实是有点抽象的,其实根据图像这个结论是显然的,在图像上,这个引理就是(下面虽然是引用但是是本人总结的):

我们称单位复数根与实轴正半轴的角为夹角,那么一个单位复数根逆时针转动两倍的夹角会等价于,它先转过180°,再转过两倍的夹角

证明 $$(\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}\omega_n^n=\omega_n^{2k}=(\omega_n^{k})^2$$

有了它,就可以进行分治了。
只要运用以上两个引理,就可以进行快速求值,因此建议读者先跳过求和引理,等到需要的时候再回头看(其实不是很难)。

2.2.3 求和引理(建议先跳过)

对于任意整数 $n\ge1$ 和不能被 $n$ 整除的非负整数 $k$ ,有
$$\displaystyle\sum_{j=0}^{n-1}(\omega_n^k)^j=0$$

这些引理一个比一个抽象……但是不难发现,如果把图中复平面内的向量全部加起来,就是一个零向量,显然成立QWQ

好吧 还是给出严谨的证明

证明 等比数列求和公式不但适用于实数,也适用于复数,因此:
$$\displaystyle\sum_{j=0}^{n-1}(\omega_n^k)^j=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}=\frac{(1)^k-1}{\omega_n^k-1}=0$$

3.0 FFT的递归实现

3.1 DFT离散傅里叶变换

被上面的一团东西搞晕了?别忘了我们的目的是什么?对,求值(最终目标是实现 $O(nlogn)$ 的多项式乘法)

普通求值的复杂度是平方级别的,那是基于随机选点(或者叫做随意选点)。

现在,我们希望计算多项式
$$A(x)=\displaystyle\sum_{j=0}^{n-1}a_jx^j$$
在 $n$ 个 $n$ 次单位复数根处的值。
现在多项式还是以系数的形式给出 $a=(a_0,a_1,\cdots,a_{n-1})$ ,定义结果为 $y_k$:
$$y_k=A(\omega_n^k)=\sum_{j=0}^{n-1}a_j\omega_n^{kj}$$

我们称 $y$ 就是 $a$ 的离散傅里叶变换(DFT)
我们记作 $y=DFT_n(a)$

3.2 FFT能够快速实现的奥秘

假设 &n& 为2的整数幂


哇,你居然发现了这里!