这篇文章上次修改于 353 天前,可能其部分内容已经发生变化,如有疑问可询问作者。

算法学习笔记3——快速傅里叶变换(FFT)

问题:

内容:对于两个多项式快速取得两个多项式的乘积。

对于最基本的解法,便是通过分配率将两个多项式中的每一项都相乘,然后取得计算结果。

对于两个n项式,其复杂度为O(N2),然而通过快速傅里叶变换可以将复杂度降低为O(NlogN)

前置知识

系数表示法

对于一个多项式,我们可以用每一项的系数表示,对于n项式我们需要(n+1)个系数,成为一组系数,对于每一组系数,这个多项式都是唯一的。这种表示我们称为系数表示法

点值表示法

然而多项式也可以通过点值来表示,例如两个端点确定一条直线。其结论是对于n次多项式,都可以用n+1个点来表示一个唯一的多项式,只需这n+1个点的坐标,便可以确定这个多项式。这种表示我们称为点值表示法

证明

对于多项式

代入d + 1个点可得该矩阵

对于该矩阵,可见范德蒙行列式,易得该多项式唯一。

解决思路

因为点值表示法也能表达一个多项式,故我们可以尝试使用点值表达来获得结果的多项式。

对于两个d次多项式,其乘积为2d次,此时只需要从前两个多项式中找到2d + 1个点,再让这2d+1个点的值相乘,便是结果多项式的2d + 1个点,再通过2d + 1个点来解决问题。

此时我们的问题就被转化为了三个步骤。

  • 由系数表达法转化为点值表达法(转换两个要相乘的多项式,并有2d + 1个点)
  • 使两个多项式的点值表达相乘得到结果表达式的点值表达法
  • 结果表达式的点值表达法转化为表达式的点值表达法

这三个步骤看似简化了O(N2)的复杂计算,实则找到2d + 1更为复杂,每个点都需要被计算,结果值将会有更多计算。

系数表达法转化为点值表达法

我们把目光转移到第一步,如何快速的找到2d + 1个点

对于这种问题,我们肯定不能去死算2d + 1次,不妨先考虑一些特殊例子

F(x) = x^2

对于这个例子,当我计算出(1, 1)是它的一个点时,另一个点也就是(-1, 1),如果我需要计算2d + 1个点,我实际上只需要计算d个点。

对于偶函数的这种特性,我们不妨对原多项式进行以下的转化 —— 取出所有的奇数次项和偶数次项

此时,当我们计算到x=1时,当x=-1也就是计算两者相加还是相减了。

同样,计算到Pe(x2)和Po(x2)时,这就是一个同理的问题(递归的味道)。我们也只需要通过递归便可找到下面的值来求上面的解。同样一个子问题的解会变成其父问题的两个解。

然而事情并没有这么简单,我们之前的预设都是在相反数的基础上,例如1和-1,它们的平方都为1故能满足我们的需求。但例如下一层计算为1的值,上一层便可或得到1和-1,若我们需要第三层,也就是8个值我们就会遇到i^2 = -1,这会阻止递归的继续,我们需要找到一个恰到好处的值,便于我们多次的递归。或者说,我们设定的子问题参数都为x的平方,都是正,故继续使用相反数便不能成立。

但我们不妨作图,发现一些规律

graph TD
1((1))
2((1))
3(( -1))
4((1))
5(( -1))
6(( -i))
7(( i))

11((x1))
12((x2))
13((x3))
14((x4))
15((x5))
16((x6))
17((x7))
18((x8))


1 --> 2
1 --> 3
2 --> 4
2 --> 5
3 --> 6
3 --> 7

4 --> 11
4 --> 12
5 --> 13
5 --> 14
6 --> 15
6 --> 16
7 --> 17
7 --> 18

一个子问题的解会变成其父问题的两个解。故我们找到的点的数量必定是2的幂次

其中对于该问题,我们要求x的8次方等于1,对于这样的x我们需要使用复数来解决。

此时我们不妨使用复数坐标系,以0作为原点,向上为i,向右为1。以1作为半径做一个圆,圆上的所有点到原点的距离为1,设定任意一个点与正方向夹角为θ,该点坐标即为(cosθ, sinθ),代表的数就是cosθ + isinθ

此时我们引入欧拉公式

此时的w便是我们参考用值的目标,接下来便可进行FFT

FFT

根据以上思路转化为伪代码

def FFT(P):
    n = len(P) # n是2的幂
    if n == 1:
        return P
    w = e**(2*pi*i/n)
    Pe, Po = [p0, p2, ...pn-2], [p1, p3, ...pn-1]
    ye, yo = FFT(Pe), FFT(Po)
    y = [0] * n
    for j in range(n/2):
    	y[j] = ye[j] + w**j*yo[j]
    	y[j + n/2] = ye[j] - w**j*yo[j]
    return y

我们同时可以再次整理和理解当前算法的目的,由于我们实际上代入w,故x也可以替换掉

中间的这个矩阵我们也就称为离散傅里叶变换矩阵(DFT)

实际上我们就是通过代入p得到对应的P的函数值,是由右边一列得到左边一列。

通常这种向量乘上离散傅里叶变换矩阵的运算我们都称为快速傅里叶变换

点值表达法转化为系数表达法

IFFT

点值表达式代入多个点,以计算多项式的系数,我们同样可以用矩阵得到表达式

此时也就是基本的矩阵运算,将另一边移过去即可,离散傅里叶变换矩阵的逆矩阵也就是

此时和我们原来的FFT函数的运算几乎一样,只是修改了其中矩阵的值,并且修改也不大,便可简单的获得其伪函数:

def IFFT(P):
    n = len(P) # n是2的幂
    if n == 1:
        return P
    w = e**(-2*pi*i/n)
    Pe, Po = [p0, p2, ...pn-2], [p1, p3, ...pn-1]
    ye, yo = IFFT(Pe), IFFT(Po)
    y = [0] * n
    for j in range(n/2):
    	y[j] = ye[j] + w**j*yo[j]
    	y[j + n/2] = ye[j] - w**j*yo[j]
    return y

参考视频