卷积

卷积的概念

回顾一下上周介绍过的卷积的概念：

对于两个序列
A→=(A0,A1,…,Ap−1,Ap)、B→=(B0,B1,…,Bm−1,Bm)\\overrightarrow {A} = (A_0,A_1,…,A_{p-1},A_p)、\\overrightarrow {B} = (B_0,B_1,…,B_{m-1},B_m)A=(A0,A1,…,Ap−1,Ap)、B=(B0,B1,…,Bm−1,Bm)

定义

C→=(A0B0,A1B0+A0B1,…,ApBm)\\overrightarrow{C} = (A_0B_0,A_1B_0+A_0B_1,…,A_pB_m)C=(A0B0,A1B0+A0B1,…,ApBm)

写成求和的形式：

Cj=∑i+k=jAiBk,0≤j≤p+mC_j=\\sum _{i+k=j}A_iB_k,\\ \\ \\ \\ \\ \\ 0\\leq j \\leq p+mCj=i+k=j∑AiBk,0≤j≤p+m

则 CjC_jCj 所组成的序列 C→=(C0,C1,…,Cp+m−1,Cp+m)\\overrightarrow {C} = (C_0,C_1,…,C_{p+m-1},C_{p+m})C=(C0,C1,…,Cp+m−1,Cp+m) 就是序列A→、B→\\overrightarrow {A}、\\overrightarrow {B}A、B 线性卷积的结果。记为
C→=A→∗B→\\overrightarrow {C}=\\overrightarrow {A}* \\overrightarrow {B}C=A∗B

刚开始光看这个定义容易让人摸不着头脑，我将从3个角度来解释卷积。

卷积的理解

1 多项式法

多项式法，也就是我们在 COMP9101 (Algorithm) Week 2 主定理和大数乘法这篇文章中提到的，建立序列 A→、B→\\overrightarrow {A}、\\overrightarrow {B}A、B 对应的多项式
PA(x)=A0+A1x+…+Ap−1xp−1+ApxpP_A(x) = A_0+A_1x+…+A_{p-1}x^{p-1}+A_px^p PA(x)=A0+A1x+…+Ap−1xp−1+Apxp

PB(x)=B0+B1x+…+Bm−1xm−1+BmxmP_B(x) = B_0+B_1x+…+ B_{m-1}x^{m-1}+B_mx^mPB(x)=B0+B1x+…+Bm−1xm−1+Bmxm
令 PC(x)=PA(x)⋅PB(x)P_C(x) = P_A(x)·P_B(x)PC(x)=PA(x)⋅PB(x) 我们可以得到

PC(x)=C0+…+Cp+m−1xp+m−1+Cp+m⋅xp+mP_C(x) = C_0+…+ C_{p+m-1}{x^{p+m-1}}+C_{p+m}·{x^{p+m}} PC(x)=C0+…+Cp+m−1xp+m−1+Cp+m⋅xp+m
写成加和的形式如下：
PC(x)=∑j=0p+mCjxjP_C(x) = \\sum ^{p+m}_{j=0} C_jx^jPC(x)=j=0∑p+mCjxj

其中
Cj=∑i+k=jAiBk,0≤j≤p+mC_j=\\sum _{i+k=j}A_iB_k,\\ \\ \\ \\ \\ \\ 0\\leq j \\leq p+mCj=i+k=j∑AiBk,0≤j≤p+m

也就是说由多项式PC(x)P_C(x)PC(x) 的系数组成的序列 C→=(C0,C1,…,Cp+m−1,Cp+m)\\overrightarrow {C} = (C_0,C_1,…,C_{p+m-1},C_{p+m})C=(C0,C1,…,Cp+m−1,Cp+m) 就是序列A→、B→\\overrightarrow {A}、\\overrightarrow {B}A、B 线性卷积的结果。

2 矩阵法

建立一个 (p+1)⋅(m+1)(p+1)·(m+1)(p+1)⋅(m+1) 的矩阵 MMM，其中 Mij=aibjM_{ij} =a_ib_jMij=aibj
如下图所示。
[a0b0a0b1…a0bm−1a0bma1b0a1b1…a1bm−1a1bm……………apb0apb0…apbm−1apbm]\\begin{bmatrix} a_{0}b_{0} & a_{0}b_{1} & \\ldots & a_{0}b_{m-1} & a_{0}b_{m} \\\\ a_{1}b_{0} & a_{1}b_{1} &\\ldots & a_{1}b_{m-1}& a_{1}b_{m}& \\\\ \\ldots & \\ldots & \\ldots &\\ldots &\\ldots \\\\ a_{p}b_{0} & a_{p}b_{0} &\\ldots & a_{p}b_{m-1} & a_{p}b_{m} \\end{bmatrix}⎣⎢⎢⎡a0b0a1b0…apb0a0b1a1b1…apb0…………a0bm−1a1bm−1…apbm−1a0bma1bm…apbm⎦⎥⎥⎤
那么序列 CCC 中的每一项都等于该矩阵对应对角线元素之和。

3 实例图示法

对于序列 A→=(A0,A1,…,An−1,An)、B→=(B0,B1,…,Bm−1,Bm)\\overrightarrow {A} = (A_0,A_1,…,A_{n-1},A_n)、\\overrightarrow {B} = (B_0,B_1,…,B_{m-1},B_m)A=(A0,A1,…,An−1,An)、B=(B0,B1,…,Bm−1,Bm)

先将序列 B→\\overrightarrow BB 进行翻转得到 B′→\\overrightarrow {B\’}B′，那么序列A→、B→\\overrightarrow {A}、\\overrightarrow {B}A、B 的卷积就相当于对序列A→\\overrightarrow {A}A中的元素以序列 B′→\\overrightarrow {B\’}B′ 中的元素为权重进行加权求和。

如下图所示。

知道了卷积的定义，我们将介绍求卷积的方法。根据 9101 week2的内容，我们可以用以下方法：

求卷积A→∗B→\\overrightarrow {A} * \\overrightarrow {B}A∗B：
step 1:
构造相应的多项式
PA(x)=A0+A1x+…+Ap−1xp−1+ApxpP_A(x) = A_0+A_1x+…+A_{p-1}x^{p-1}+A_px^p PA(x)=A0+A1x+…+Ap−1xp−1+Apxp

PB(x)=B0+B1x+…+Bm−1xm−1+BmxmP_B(x) = B_0+B_1x+…+ B_{m-1}x^{m-1}+B_mx^mPB(x)=B0+B1x+…+Bm−1xm−1+Bmxm

step 2:
代入 p+m+1p+m+1p+m+1 个 xxx 的值，得到相应的 PA(x)P_A(x)PA(x) 和 PB(x)P_B(x)PB(x) ，对应项相乘得到 p+m+1p+m+1p+m+1 个 PC(x)P_C(x)PC(x) 的值

step 3:
建立关于多项式PC(x)P_C(x)PC(x) 系数的线性方程，对线性方程进行求解得到 PC(x)P_C(x)PC(x) 系数的值

【在 step 2 赋值的过程中，我们需要计算 xix_ixi 的 ppp 次幂，这就造成了数值大小的激增，所以我们引入单位 1 的复数根，来简化这种幂运算。】

单位 1 的复数根

在复平面上，复数 zzz 的表示形式：

z=a+biz = a + biz=a+bi

z=∣z∣eiarg(z)=∣z∣(cosarg(z)+isinarg(z))z = |z|e^{i\\ arg(z)} =|z|(cos\\ arg(z) + i\\ sin\\ arg(z))z=∣z∣eiarg(z)=∣z∣(cosarg(z)+isinarg(z))

其中 ∣z∣=a2+b2|z| =\\sqrt {a^{2}+b^{2}}∣z∣=a2+b2， arg(z)∈(−π,π]arg(z) \\in (-\\pi,\\pi ]arg(z)∈(−π,π]

如下图所示

那么
z2=(∣z∣eiarg(z))2=∣z∣2(cos2arg(z)+isin2arg(z))z^2 = (|z|e^{i\\ arg(z)})^2 = |z|^2(cos\\ 2\\ arg(z) + i\\ sin\\ 2\\ arg(z))z2=(∣z∣eiarg(z))2=∣z∣2(cos2arg(z)+isin2arg(z))

zn=(∣z∣eiarg(z))n=∣z∣n(cosnarg(z)+isinnarg(z))z^n = (|z|e^{i\\ arg(z)})^n = |z|^n(cos\\ n\\ arg(z) + i\\ sin\\ n\\ arg(z))zn=(∣z∣eiarg(z))n=∣z∣n(cosnarg(z)+isinnarg(z))

单位 1 的 nnn 次复数根就是满足 zn=1z^n = 1zn=1 的 zzz，表示为：

∣z∣=1，narg(z)=2kπ，k∈Z|z| = 1， \\ \\ n\\ arg(z) = 2k\\pi，k\\in Z∣z∣=1，narg(z)=2kπ，k∈Z

arg(z)=2kπn，k∈Zarg(z) = \\dfrac {2k\\pi}{n}，k\\in Zarg(z)=n2kπ，k∈Z

所以单位 111 的所有n次复数根可以表示为：
wnk=(ei2π/n)k,k∈[0,n−1]w_n^k = (e^{i\\ 2\\pi /n})^k,k\\in [0,n-1]wnk=(ei2π/n)k,k∈[0,n−1]

假设 n=16n =16n=16 ，那么单位 1 的16次复数根如下图所示。

复数根的性质：
1 任意 2 个相同阶次的复数根相乘得到的结果仍然是该阶次的复数根。
2 对于所有的整数 kmnk\\ m\\ nkmn，wknkm=wnmw_{kn}^{km}=w_n^mwknkm=wnm

所以，我们获取 p+m+1p+m+1p+m+1 个 PC(x)P_C(x)PC(x) 值的方式如下：

根据序列 A=⟨A0,A1,…,Ap⟩A =\\ \\langle A_0,A_1,…,A_{p}\\rangleA=⟨A0,A1,…,Ap⟩ 代入 111 的 p+m+1p+m+1p+m+1 次复数根，我们就能得到对应的PA(wp+m+10)P_A(w_{p+m+1}^0)PA(wp+m+10), PA(wp+m+11)P_A(w_{p+m+1}^1)PA(wp+m+11), PA(wp+m+12)P_A(w_{p+m+1}^2)PA(wp+m+12)…

同理我们也可以得到PB(wp+m+10)P_B(w_{p+m+1}^0)PB(wp+m+10), PB(wp+m+11)P_B(w_{p+m+1}^1)PB(wp+m+11), PB(wp+m+12)P_B(w_{p+m+1}^2)PB(wp+m+12)…进而求出相应的PC(x)P_C(x)PC(x) 的值

【事实上，这是对序列 AAA 进行离散傅立叶变换的过程。为了凑齐 p+m+1p+m+1p+m+1 个值，我们需要在序列 AAA 后面补上mmm 个 000 。也就是说，为了求序列AAA , BBB 的卷积，我们需要将序列 AAA , BBB 补 000 后求二者的离散傅立叶变换】

离散傅立叶变换 DFT

令 AAA 为 nnn 个复数组成的序列 A=⟨A0,A1,…,An−1⟩A =\\ \\langle A_0,A_1,…,A_{n-1}\\rangleA=⟨A0,A1,…,An−1⟩

AAA 对应的多项式为 PA(x)=∑j=0n−1AjxjP_A(x) = \\sum _{j=0}^{n-1}A_jx^jPA(x)=j=0∑n−1Ajxj
将单位 1 的 nnn 个复数根代入到这个多项式中得到的新的序列: ⟨PA(1),PA(wn),PA(wn2),…,PA(wnn−1)⟩\\langle P_A(1),P_A(w_n),P_A(w_n^2),…,P_A(w_n^{n-1})\\rangle⟨PA(1),PA(wn),PA(wn2),…,PA(wnn−1)⟩ 就叫做序列 AAA 的离散傅立叶变换 Discrete Fourier Transform （DFT），记作
A^=⟨A^0,A^1,…,A^n−1⟩\\widehat {A} =\\langle \\widehat {A}_0,\\widehat {A}_1,…,\\widehat {A}_{n-1}\\rangleA=⟨A0,A1,…,An−1⟩

对于每一项，我们都需要进行 nnn 次系数AjA_jAj与wnjw_n^jwnj的乘法运算，所以计算离散傅立叶变换所需的时间为 O(n2)O(n^2)O(n2)

快速傅立叶变换 FFT

【快速傅立叶变换是一种利用分治法求离散傅立叶变换的方法，该方法的时间复杂度为 O(nlogn)O(nlogn)O(nlogn)】

为了求给定序列 A=⟨A0,A1,…,An−1⟩A =\\ \\langle A_0,A_1,…,A_{n-1}\\rangleA=⟨A0,A1,…,An−1⟩ 的傅立叶变换
我们考虑把多项式 PA(x)P_A(x)PA(x) 按照 xxx 的指数的奇偶性分成两个部分

A0+A2x2+A4x4+…+An−2xn−2A_0+A_2x^2+A_4x^4 + …+ A_{n-2}x^{n-2}A0+A2x2+A4x4+…+An−2xn−2

和

A1x+A3x3+A5x5+…+An−1xn−1 A_1x+A_3x^3 +A_5x^5+ …+ A_{n-1}x^{n-1}A1x+A3x3+A5x5+…+An−1xn−1
定义

PAeven(x)=A0+A2x+A4x2+…+An−2x(n−2)/2P_A^{even}(x)=A_0+A_2x+A_4x^2 + …+ A_{n-2}x^{(n-2)/2}PAeven(x)=A0+A2x+A4x2+…+An−2x(n−2)/2

PAodd(x)=A1+A3x+A5x2+…+An−1x(n−2)/2P_A^{odd}(x) = A_1+A_3x +A_5x^2+ …+ A_{n-1}x^{(n-2)/2}PAodd(x)=A1+A3x+A5x2+…+An−1x(n−2)/2
那么 PA(x)=PAeven(x2)+x⋅PAodd(x2)P_A(x) = P_A^{even}(x^2) + x·P_A^{odd}(x^2)PA(x)=PAeven(x2)+x⋅PAodd(x2)

所以，序列 AAA 的 DFTDFTDFT 的结果为：

PA(wnk)=PAeven((wnk)2)+wnk⋅PAodd((wnk)2)P_A(w_n^k) = P_A^{even}(({w_n^k})^2)+w_n^k·P_A^{odd}(({w_n^k})^2)PA(wnk)=PAeven((wnk)2)+wnk⋅PAodd((wnk)2)

也就是
PA(wnk)=PAeven(wn/2k)+wnk⋅PAodd(wn/2k),k∈[0,n−1]P_A(w_n^k) = P_A^{even}({w_{n/2}^k})+w_n^k·P_A^{odd}({w_{n/2}^k}), \\ \\ k\\in [0,n-1]PA(wnk)=PAeven(wn/2k)+wnk⋅PAodd(wn/2k),k∈[0,n−1]

等同于当 k∈[0,n/2−1]k\\in [0,n/2-1]k∈[0,n/2−1] 时：
PA(wnk)=PAeven(wn/2k)+wnk⋅PAodd(wn/2k)P_A(w_n^k) = P_A^{even}({w_{n/2}^k})+w_n^k·P_A^{odd}({w_{n/2}^k})PA(wnk)=PAeven(wn/2k)+wnk⋅PAodd(wn/2k)

PA(wnn/2+k)=PAeven(wn/2k)−wnk⋅PAodd(wn/2k)P_A(w_n^{n/2+k}) = P_A^{even}({w_{n/2}^k})-w_n^k·P_A^{odd}({w_{n/2}^k})PA(wnn/2+k)=PAeven(wn/2k)−wnk⋅PAodd(wn/2k)
完整步骤如下：

所以 FFTFFTFFT 时间复杂度为：
T(n)=T(n2)+O(n)T(n) = T(\\dfrac{n}{2}) + O(n)T(n)=T(2n)+O(n)

T(n)∈Θ(nlogn)T(n)\\in \\Theta(nlogn)T(n)∈Θ(nlogn)

【快速傅立叶逆变换】：
我们知道已知序列 A=⟨A0,A1,…,An−1⟩A =\\ \\langle A_0,A_1,…,A_{n-1}\\rangleA=⟨A0,A1,…,An−1⟩ 可以对其进行快速傅立叶变换来获取对应多项式的值 PA(wnk)P_A(w_n^k)PA(wnk)。那么我们也可以通过 PC(wnk)P_C(w_n^k)PC(wnk) 来获取对应的系数序列 C=⟨C0,C1,…,Cn−1⟩C =\\ \\langle C_0,C_1,…,C_{n-1}\\rangleC=⟨C0,C1,…,Cn−1⟩, 这个过程可以采用快速傅立叶逆变换。
具体步骤如下：