深度学习算子优化-FFT

作者：严健文 | 旷视 MegEngine 架构师

背景

在数字信号和数字图像领域，对频域的研究是一个重要分支。
我们日常“加工”的图像都是像素级，被称为是图像的空域数据。空域数据表征我们“可读”的细节。如果我们将同一张图像视为信号，进行频谱分析，可以得到图像的频域数据。观察下面这组图 (来源)，频域图中的亮点为低频信号，代表图像的大部分能量，也就是图像的主体信息。暗点为高频信号，代表图像的边缘和噪声。从组图可以看出，Degraded Goofy 与 Goofy 相比，近似的低频信号保留住了 Goofy 的“轮廓”，而其高频信号的增加使得背景噪点更加明显。频域分析使我们可以了解图像的组成，进而做更多的抽象分析和细节处理。

Goofy and Degraded Goofy

实现图像空域和频域转换的工具，就是傅立叶变换。由于图像数据在空间上是离散的，我们使用傅立叶变换的离散形式 DFT（Discrete Fourier Transform）及其逆变换 IDFT（Inverse Discrete Fourier Transform)。Cooley-Tuckey 在 DFT 的基础上，开发了更快的算法 FFT（Fast Fourier Transform）。

DFT/FFT 在数字图像领域还有一些延伸应用。比如基于 DFT 的 DCT（Discrete Cosine Transform，离散余弦变换）就用在了图像压缩 JPEG 算法 (来源) 和图像水印算法（来源）。

JPEG 编码是通过色彩空间转换、抽样分块、DCT 变换、量化编码实现的。其中 DCT 变换的使用将图像低频信息和高频信息区分开，在量化编码过程中压缩了少量低频信息、大量高频信息从而获得尺寸上压缩。从猫脸图上可看出随着压缩比增大画质会变差，但是主体信息还是得以保留。

猫脸不同 jpeg 画质（压缩比）

图像水印算法通过 DCT 将原图转换至频域，选取合适的位置嵌入水印图像信息，并通过 IDCT 转换回原图。这样对原图像的改变较小不易察觉，且水印通过操作可以被提取。

DFT/FFT 在深度学习领域也有延伸应用。比如利用 FFT 可以降低卷积计算量的特点，FFT_Conv 算法也成为常见的深度学习卷积算法。本文我们就来探究一下频域算法的原理和优化策略。

DFT 的原理及优化

公式

无论是多维的 DFT 运算，还是有基于 DFT 的 DCT/FFT_Conv, 底层的计算单元都是 DFT_1D。因此，DFT_1D 的优化是整个 FFT 类算子优化的基础。
DFT_1D 的计算公式：
$X_{k}=\sum_{n=0}^{\mathrm{N}-1} x_{n} e^{-j 2 \pi k \frac{n}{N}} \quad k=0, \ldots, N-1$

其中 $x_{n}$ 为长度为 N 的输入信号， $e^{-j 2 \pi k \frac{n}{N}}$ 是 1 的 N 次根， $X_{k}$ 为长度为 N 的输出信号。
该公式的矩阵形式为：

$\left[\begin{array}{c}X(0) \\ X(1) \\ \vdots \\ X(N-1)\end{array}\right]=\left[W_{N}^{n k}\right]\left[\begin{array}{c} \left.x(0\right) \\ x(1) \\ \vdots \\ x(N-1)\end{array}\right]$

单位复根的性质

DFT_1D 中的 $W_{N}^{nk} = e^{-j 2 \pi k \frac{n}{N}}$ 是 1 的单位复根。直观地看，就是将复平面划分为 N 份，根据 k * n 的值逆时针扫过复平面的圆周。

单位复根有着周期性和对称性，我们依据这两个性质可以对 W 矩阵做大量的简化，构成 DFT_1D 的快速算法的基础。
周期性： $W_{N}^{k +N}=W_{N}^{k}$
对称性： $W_{N}^{k+N / 2}=-W_{N}^{k}$

Cooley-Tuckey FFT 算法

DFT_1D 的多种快速算法中，使用最频繁的是 Cooley-Tuckey FFT 算法。算法采用用分治的思想，将输入尺寸为 N 的序列，按照不同的基 radix，分解为 N/radix 个子序列，并对每个子序列再划分，直到不能再被划分为止。每一次划分都可以得到一级 stage，将所有的级自下而上组合在一起，计算得到最后的输出序列。
这里以 N = 8， radix=2 为例展示推理过程。
其中 $x(k)$ 为 N=8 的序列， $X^{F}(k)$ 为 DFT 输出序列。
根据 DFT 的计算公式
$X^{F}(k)=W_{8}^{0} x_{0}+W_{8}^{k} x_{1}+W_{8}^{2 k} x_{2}+W_{8}^{3k} x_{3}+W_{8}^{4k} x_{4} + W_{8}^{5k} x_{5}+W_{8}^{6k} x_{6} +W_{8}^{7k} x_{7}$

根据奇偶项拆开，分成两个长度为 4 的序列 $G(k)$ , $H(k)$ 。

$X^{F}(k)= W_{8}^{0} x_{0}+W_{8}^{2 k} x_{2}+W_{8}^{4 k} x_{4}+W_{8}^{6 k} x_{6}$

$+W_{8}^{k}\left(W_{8}^{0} x_{1}+W_{8}^{2 k} x_{3}+W_{8}^{4 k} x_{5}+W_{8}^{6 k} x_{7}\right)$
$=G^{F}(k)+W_{8}^{k} H^{F}(k)$

$X^{F}(k+4)=W_{8}^{0} x_{0}+W_{8}^{2(k+4)} x_{2}+W_{8}^{4(k+4)} x_{4}+W_{8}^{6(k+4)} x_{6}$
$+W_{8}^{(k+4)}\left(W_{8}^{0} x_{1}+W_{8}^{2(k+4)} x_{3}+W_{8}^{4(k+4)} x_{5}+W_{8}^{6(k+4)} x_{7}\right)$
$=G^{F}(k)+W_{8}^{k+4} H^{F}(k)$
$=G^{F}(k)-W_{8}^{k} H^{F}(k)$

$G^{F}(k)$ 和 $H^{F}(k)$ 为 $G(k)$ 和 $H(k)$ 的 DFT 结果。 $G^{F}(k)$ 和 $H^{F}(k)$ 乘以对应的旋转因子 $W_{8}^{k}$ ，进行简单的加减运算可以得到输出 $X^{F}(k)$ 。
同理，对 $G(k)$ 和 $H(k)$ 也做一样的迭代， $A(k)$ ， $B(k)$ , $C(k)$ , $D(k)$ 都是 N=2 的序列，用他们的 DFT 结果进行组合运算可以得到 $G^{F}(k)$ 和 $H^{F}(k)$ 。

$\begin{aligned} &G^{F}(k)=A^{F}(k) + W_{4}^{k}B^{F}(k)\\ \end{aligned}$
$\begin{aligned} &G^{F}(k+2)=A^{F}(k)-W_{4}^{k}B^{F}(k)\\ \end{aligned}$
$\begin{aligned} &H^{F}(k)=C^{F}(k)+W_{4}^{k}D^{F}(k)\\ \end{aligned}$
$\begin{aligned} &H^{F}(k+2)=C^{F}(k)-W_{4}^{k}D^{F}(k)\\ \end{aligned}$

计算 N=2 的序列 $A^{F}(k)$ , $B^{F}(k)$ , $C^{F}(k)$ , $D^{F}(k)$ , 因为 $k=0$ ，旋转因子 $W_{2}^{0}$ = 1。只要进行加减运算得到结果。
$\left[\begin{array}{l} A^{F}(0) \\ A^{F}(1) \end{array}\right]=\left[\begin{array}{ll} 1 & 1 \\ 1 & -1 \end{array}\right]\left[\begin{array}{l} x_{0} \\ x_{4} \\ \end{array}\right]$

$\left[\begin{array}{l} B^{F}(0) \\ B^{F}(1) \end{array}\right]=\left[\begin{array}{ll} 1 & 1 \\ 1 & -1 \end{array}\right]\left[\begin{array}{l} x_{2} \\ x_{6} \\ \end{array}\right]$

$\left[\begin{array}{l} C^{F}(0) \\ C^{F}(1) \end{array}\right]=\left[\begin{array}{ll} 1 & 1 \\ 1 & -1 \end{array}\right]\left[\begin{array}{l} x_{1} \\ x_{5} \\ \end{array}\right]$

$\left[\begin{array}{l} D^{F}(0) \\ D^{F}(1) \end{array}\right]=\left[\begin{array}{ll} 1 & 1 \\ 1 & -1 \end{array}\right]\left[\begin{array}{l} x_{3} \\ x_{7} \\ \end{array}\right]$

用算法图形表示，每一层的计算会产生多个蝶形，因此该算法又被称为蝶形算法。
这里我们要介绍碟形网络的基本组成，对下文的分析有所帮助。

N=8 碟形算法图

N=8 的计算序列被分成了 3 级，每一级 (stage) 有一个或多个块 (section)，每个块中包含了一个或者多个蝶形（butterfly), 蝶形的计算就是 DFT 运算的 kernel。
每一个 stage 的计算顺序：

取输入
乘以转换因子
for section_num, for butterfly_num，执行 radixN_kernel
写入输出。

看 N=8 的蝶形算法图，stage = 1 时，运算被分成了 4 个 section，每个 section 的 butterfly_num = 1。stage = 2 时，section_num = 2，butterfly_num = 2。 stage = 3 时，section_num = 1，butterfly_num = 4。
可以观察到，从左到右过程中 section_num 不断减少，butterfly_num 不断增加，蝶形群在“变大变密”，然而每一级总的碟形次数是不变的。
实际上，对于长度为 N，radix = r 的算法，我们可以推得到：

$\text { Sec_num }=N / r^{S}$
$\text { Butterfly_num }= r^{S-1}$
$\text { Sec_stride }=r^{S}$
$\text { Butterfly_stride }=1$

S 为当前的 stage，sec/butterfly_stride 是每个 section/butterfly 的间隔。

这个算法可以将复杂度从 O(n^2) 下降到 O(nlogn)，显得高效而优雅。我们基于蝶形算法，对于不同的 radix 进行算法的进一步划分和优化，主要分为 radix - 2 的幂次的和 radix – 非 2 的幂次两类。

radix-2 的幂次优化

DFT_1D 的 kernel 即为矩阵形式中的 $W_{N}^{nk}$ 矩阵，我们对 radix_2^n 的 kernel 进行分析。

背景里提到， DFT 公式的矩阵形式为：
$\left[\begin{array}{c}X(0) \\ X(1) \\ \vdots \\ X(N-1)\end{array}\right]=\left[W_{N}^{n k}\right]\left[\begin{array}{c} \left.x(0\right) \\ x(1) \\ \vdots \\ x(N-1)\end{array}\right]$
其中 $x(0)$ ~ $x(N-1)$ 为乘以旋转因子 $W_{N}^{kn}$ 后的输入

当 radix = 2 时，由于 $W_{2}^1$ = -1, $W_{2}^2$ = 1, radix_2 的 DFT 矩阵形式可以写为：

$\left[\begin{array}{c}\mathrm{X}_{\mathrm{k}} \\ \mathrm{X}_{\mathrm{k}+\mathrm{N} / 2}\end{array}\right]$ $=\left[\begin{array}{cc}1 & 1 \\ 1 & -1\end{array}\right]\left[\begin{array}{l}\mathrm{W}_{\mathrm{N}}^{0} \mathrm{A}_{\mathrm{k}} \\ \mathrm{W}_{\mathrm{N}}^{\mathrm{k}} \mathrm{B}_{\mathrm{k}}\end{array}\right]$

当 radix = 4 时，由于 $W_{4}^1$ = -j, $W_{4}^2$ = -1, $W_{4}^3$ = j, $W_{4}^4$ = 1，radix_4 的 DFT 矩阵形式可以写为：

$\left[\begin{array}{c}\mathrm{X}_{\mathrm{k}} \\ \mathrm{X}_{\mathrm{k}+\mathrm{N} / 4} \\ \mathrm{X}_{\mathrm{k}+\mathrm{N} / 2} \\ \mathrm{X}_{\mathrm{k}+3 \mathrm{~N} / 4}\end{array}\right]=\left[\begin{array}{cccc}1 & 1 & 1 & 1 \\ 1 & -\mathrm{j} & -1 & \mathrm{j} \\ 1 & -1 & 1 & -1 \\ 1 & \mathrm{j} & -1 & -\mathrm{j}\end{array}\right]\left[\begin{array}{c}\mathrm{W}_{\mathrm{N}}^{0} \mathrm{A}_{\mathrm{k}} \\ \mathrm{W}_{\mathrm{N}}^{\mathrm{k}} \mathrm{B}_{\mathrm{k}} \\ \mathrm{W}_{\mathrm{N}}^{2 \mathrm{k}} \mathrm{C}_{\mathrm{k}} \\ \mathrm{W}_{\mathrm{N}}^{3 \mathrm{k}} \mathrm{D}_{\mathrm{k}}\end{array}\right]$

同理推得到 radix_8 的 kernel 为：

$\left[\begin{array}{cccccccc}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & \mathrm{~W}_{8}^{1} & -j & \mathrm{~W}_{8}^{3} & -1 & -\mathrm{W}_{8}^{1} & j & -\mathrm{W}_{8}^{3} \\ 1 & -j & -1 & j & 1 & -j & -1 & j \\ 1 & \mathrm{~W}_{8}^{3} & j & \mathrm{~W}_{8}^{1} & -1 & -\mathrm{W}_{8}^{3} & -j & -\mathrm{W}_{8}^{1} \\ 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 \\ 1 & -\mathrm{W}_{8}^{1} & -j & -\mathrm{W}_{8}^{3} & -1 & \mathrm{~W}_{8}^{1} & j & \mathrm{~W}_{8}^{3} \\ 1 & j & -1 & -j & 1 & j & -1 & -j \\ 1 & -\mathrm{W}_{8}^{3} & j & -\mathrm{W}_{8}^{1} & -1 & \mathrm{~W}_{8}^{3} & -j & \mathrm{~W}_{8}^{1}\end{array}\right]$

我们先来看访存，现代处理器对于计算性能的优化要优于对于访存的优化，在计算和访存相近的场景下，访存通常是性能瓶颈。

DFT1D 中，对于不同基底的算法 r-2/r-4/r-8, 每一个 stage 有着相等的存取量：2 * butterfly_num * radix = 2N, 而不同的基底对应的 stage 数有着明显差异（ $\log_2N$ vs $\log_4N$ vs $\log_8N$ )。

因此对于 DFT, 在不显著增加计算量的条件下，选用较大的 kernel 会在访存上取得明显的优势。观察推导的 kernel 图， r-2 的 kernel 每个蝶形对应 4 次访存操作和，2 次复数浮点加减运算。r-4 的 kernel 每个蝶形算法对应 8 次 load/store、8 次复数浮点加减操作（合并相同的运算），在计算量略增加的同时 stage 由 $\log_2N$ 下降到 $\log_4N$ , 降低了总访存的次数，因此会有性能的提升。r-8 的 kernel 每个蝶形对应 16 次 load/store、24 次复数浮点加法和 8 次浮点乘法。浮点乘法的存在使得计算代价有所上升， stage 由 $\log_4N$ 进一步下降到 $\log_8N$ ，但由于 N 日常并不会太大， r-4 到 r-8 的 stage 减少不算明显，所以优化有限

我们再来看计算的开销。减少计算的开销通常有两种办法：减少多余的运算、并行化。

以 r-4 算法为例，kernel 部分的计算为：

radix_4_first_stage(src, dst, sec_num, butterfly_num)
radix_4_other_stage(src, dst, sec_num, butterfly_num)
- for Sec_num
  - for butterfly_num
    - raidx_4_kernel

radix4_first_stage 的数据由于 k=0, 旋转因子都为 1，可以省去这部分复数乘法运算，单独优化。 radix4_other_stage 部分，从第 2 个 stage 往后， butterfly_num = 4^(s-1) 都为 4 的倍数，而每个 butterfly 数组读取/存储都是间隔的。可以对最里层的循环做循环展开加向量化，实现 4 个或更多 butterfly 并行运算。循环展开和 SIMD 指令的使用不仅可以提高并行性，也可以提升 cacheline 利用的效率，可以带来较大的性能提升。以 SM8150(armv8) 为例，r-4 的并行优化可以达到 r2 的 1.6x 的性能。

尺寸：1 * 2048（r2c) 环境：SM8150 大核

总之，对于 radix-2^n 的优化，选用合适的 radix 以减少多 stage 带来的访存开销，并且利用单位复根性质以及并行化降低计算的开销，可以带来较大的性能提升。

radix-非 2 的幂次优化

当输入长度 N = radix1^m1 * radix2^m2... 且 radix 都不为 2 的幂次时，如果使用 naive 的 O(n^2) 算法，性能就会急剧下降。常见的解决办法对原长补 0、使用 radix_N 算法、特殊的 radix_N 算法 (chirp-z transform)。补 0 至 2 的幂次方法对于大尺寸的输入要增加很多运算量和存储量，而 chirp-z transform 是用卷积计算 DFT，算法过于复杂。因此对非 2 的幂次 radix-N 的优化也是必要的。

radix-N 计算流程和 radix-2 幂次一样，我们同样可以利用单位复根的周期性和对称性，对 kernel 进行计算的简化。以 radix-5 为例，radix-5 的 DFT_kernel 为：
$\left[\begin{array}{cccc} 1&1&1&1&1\\ 1 &\mathrm{W}_{\mathrm{5}}^{1} & \mathrm{W}_{\mathrm{5}}^{2} & \mathrm{W}_{\mathrm{5}}^{-2} & \mathrm{W}_{\mathrm{5}}^{-1} \\ 1 &\mathrm{W}_{\mathrm{5}}^{2} & \mathrm{W}_{\mathrm{5}}^{-1} & \mathrm{W}_{\mathrm{5}}^{1} & \mathrm{W}_{\mathrm{5}}^{-2} \\ 1 &\mathrm{W}_{\mathrm{5}}^{-2} & \mathrm{W}_{\mathrm{5}}^{1} & \mathrm{W}_{\mathrm{5}}^{-1} & \mathrm{W}_{\mathrm{5}}^{2} \\ 1 &\mathrm{W}_{\mathrm{5}}^{-1} & \mathrm{W}_{\mathrm{5}}^{-2} & \mathrm{W}_{\mathrm{5}}^{2} & \mathrm{W}_{\mathrm{5}}^{1} \\ \end{array}\right]$

$W_5^k$ 和 $W_{5}^{-k}$ 在复平面上根据 x 轴对称，有相同的实部和相反的虚部。根据这个性质。如下图所示，对于每一个 stage，可以合并公共项 A，B，C，D，再根据公共项计算出该 stage 的输出。

$\begin{array}{l} A=\left(x_{1}+x_{4}\right) * W_{5}^{1} \cdot r+\left(x_{2}+x_{3}\right) * W_{5}^{2} \cdot r\\\end{array}$

$B=(-j) * \left[\left(x_{1}-x_{4}\right) * W_{5}^{1} \cdot i+\left(x_{2}-x_{3}\right) * W_{5}^{2} \cdot i\right] \$

$C=\left(x_{1}+x_{4}\right) * W_{5}^{2} \cdot r+\left(x_{2}+x_{3}\right) * W_{5}^{1} \cdot r\$

$D=j * \left[\left(x_{1}-x_{4}\right) * W_{5}^{2} \cdot i-\left(x_{2}-x_{3}\right) * W_{5}^{1} \cdot i\right] \$

$\begin{array}{l} X(k)=x_{0}+\left(x_{1}+x_{4}\right)+\left(x_{2}+x_{3}\right)\\ \end{array}$
$\begin{array}{l} X(k+N/5)=x_{0}+\mathrm{A}-\mathrm{B}\\ X(k+2N/5)=x_{0}+\mathrm{C}+\mathrm{D}\\ X(k+3N/5)=x_{0}+C-D\\ X(k+4N/5)=x_{0}+\mathrm{A}+\mathrm{B}\\ \end{array}$

这种算法减少了很多重复的运算。同时，在 stage>=2 的时候，同样对 butterfly 做循环展开加并行化，进一步减少计算的开销。
radix-5 的优化思想可以外推至 radix-N。对于 radix_N 的每一个 stage, 计算流程为：

取输入
乘以对应的转换因子
计算公共项， radix_N 有 N-1 个公共项
执行并行化的 radix_N_kernel
写入输出

其他优化

上述两个章节描述的是 DFT_1D 的通用优化，在此基础上还可以做更细致的优化，可以参考本文引用的论文。

对于全实数输入的，由于输入的虚部为 0，进行旋转因子以及 radix_N_kernel 的复数运算时会有多余的运算和多余的存储，可以利用 split r2c 算法，视为长度为 N/2 的复数序列，计算 DFT 结果并进行 split 操作得到 N 长实数序列的结果。
对于 radix-2 的幂次算法，重新计算每个 stage 的输入/输出 stride 以取消第一级的位元翻转可以进一步减少访存的开销。
对于 radix-N 算法，在混合基框架下 N = radix1^m1 * radix2^m2，合并较小的 radix 为大的 radix 以减少 stage。

DFT 延展算法的原理及优化

DCT 和 FFT_conv 两个典型的基于 DFT 延展的算法，DFT_1D/2D 的优化可以很好的用在这类算法中。

DCT

DCT 算法（Discrete Cosine Transform，离散余弦变换）可以看作是 DFT 取其正弦分量并经过工业校正的算法。DFT_1D 的计算公式为：

$\begin{aligned} X[k] &=C(k) \sum_{n=0}^{N-1} x[n] \cos \left(\frac{(2 n+1) \pi k}{2 N}\right) \\ &C(k)=\sqrt{\frac{1}{n}} \\&k=1 \\ &C(k)=\sqrt{\frac{2}{n}} \\&k!=1 \\ \end{aligned}$

该算法 naive 实现是 O(n^2) 的，而我们将其转换成 DFT_1D 算法，可以将算法复杂度降至 O(nlogn)。
基于 DFT 的 DCT 算法流程为：

对于 DCT 的输入序列 x[n], 创建长为 2N 的输入序列 y[n] 满足 y[n] = x[n] + x[2N-n-1], 即做一个镜像对称。
对输入序列 y[n] 进行 DFT 运算，得到输出序列 Y[K]。
由 Y[K] 计算得到原输入序列的输出 X[K] 。

我们尝试推导一下这个算法：

${l} y[n]=x[n]+x [2 N-1-n] \$

${l} Y[k]=\sum_{n=0}^{N-1} x[n]\cdot e^{-j \frac{2 \pi k n}{2 N}} +\sum_{n=N}^{2 N-1} x[2 N-1-n] \cdot e^{-j \frac{2 \pi k n}{2 N}}$

$=\sum_{n=0}^{N-1} x[n]\cdot e^{-j \frac{2 \pi k n}{2 N}} +\sum_{n=0}^{N-1} x[n] \cdot e^{-j \frac{2 \pi k (2N-1-n)}{2 N}}$
$=e^{-j \frac{2 \pi k }{2 N}} \cdot \sum_{n=0}^{N-1} x[n] (e^{-j \frac{2\pi}{2 N} kn} \cdot e^{-j \frac{\pi}{2 N}k}+e^{j \frac{2\pi}{2 N} kn} \cdot e^{j \frac{\pi}{2 N}k})$
$=e^{-j \frac{2 \pi k }{2 N}} \cdot \sum_{n=0}^{N-1} x[n] \cdot 2\cdot\cos(\frac{2n+1}{2N} k\pi)$
$=e^{-j \frac{2 \pi k }{2 N}} \cdot C(u) \cdot X[k]$

对 y[n] 依照 DFT 公式展开，整理展开的两项并提取公共项 $e^{-j \frac{2 \pi k }{2 N}}$ , 根据欧拉公式和诱导函数，整理非公共项 $(e^{-j \frac{2\pi}{2 N} kn} \cdot e^{-j \frac{\pi}{2 N}k}+e^{j \frac{2\pi}{2 N} kn} \cdot e^{j \frac{\pi}{2 N}k})$ 。可以看出得到的结果正是 x[k] 和与 k 有关的系数的乘积。这样就可以通过先计算 $Y[k]$ 得到 x[n] 的 DCT 输出 $X[k]$ 。

在理解算法的基础上，我们对 DFT_1D 的优化可以完整地应用到 DCT 上。DCT_2D 的计算过程是依次对行、列做 DCT_1D, 我们用多线程对 DCT_1D 进行并行，可以进一步优化算法。

FFT_conv

Conv 是深度学习最常见的运算，计算 conv 常用的方法有 IMG2COL+GEMM， Winograd， FFT_conv。三种算法都有各自的使用场景。

FFT_conv 的数学原理是时域中的循环卷积对应于其离散傅里叶变换的乘积。如下图所示， f 和 g 的卷积等同于将 f 和 g 各自做傅立叶变幻 F，进行点乘并通过傅立叶逆变换计算后的结果。
$f \underset{\text { Circ }}{*} g=\mathcal{F}^{-1}(\mathcal{F}(f) \cdot \mathcal{F}(g))$

直观的理论证明可下图（来源)。

$\mathcal{F}[f * g] \$

$=\int_{-\infty}^{\infty}\left[\left(\int_{-\infty}^{\infty}g(z)f(x-z)dz\right)e^{-i k x}\right]dx$

$=\int_{-\infty}^{\infty} g(z)\left[\int_{-\infty}^{\infty} f(x-z) e^{-i k x} d x\right] d z$
$=\int_{-\infty}^{\infty} g(z)\left[\int_{-\infty}^{\infty} f(y) e^{-i k(y+z)} d y\right] d z$
$=\left[\int_{-\infty}^{\infty} g(z) e^{-i k z} d z\right]\left[\int_{-\infty}^{\infty} f(y) e^{-i k y} d y\right]$
$=\mathcal{F}[f] \cdot \mathcal{F}[g]$

将卷积公式和离散傅立叶变换展开，改变积分的顺序并且替换变量，可以证明结论。
注意这里的卷积是循环卷积，和我们深度学习中常用的线性卷积是有区别的。利用循环卷积计算线性卷积的条件为循环卷积长度 L⩾| f |+| g |−1。因此我们要对 Feature Map 和 Kernel 做 zero-padding，并从最终结果中取有效的线性计算结果。

FFT_conv 算法的流程：

将 Feature Map 和 Kernel 都 zero-pad 到同一个尺寸，进行 DFT 转换。
矩阵点乘
将计算结果通过 IDFT 计算出结果。

该算法将卷积转换成点乘，算法复杂度是 O(nlogn), 小于卷积的 O(n^2), 在输入的尺寸比较大时可以减少运算量，适用于大 kernel 的 conv 算法。

深度学习计算中， Kernel 的尺寸要远小于 Feature Map，因此 FFT_conv 第一步的 zero-padding 会有很大的开销，参考论文 2 里提到可以通过对 Feature map 进行分块，分块后的 Feature Map 和 Kernel 需要 padding 到的尺寸较小，可以大幅减小这一部分的开销。优化后 fft_conv 的计算流程为：

合理安排缓存计算出合适的 tile 尺寸，对原图进行分块
分块后的小图和 kernel 进行 zero-padding, 并进行 DFT 运算
小图矩阵点乘
进行逆运算并组合成大图。

同时我们可以观察到，FFT_conv 的核心计算模块还是针对小图的 DFT 运算，因此我们可以将前一章节对 DFT 的优化代入此处，辅以多线程，进一步提升 FFT_Conv 的计算效率。

参考资料

陈暾，李志豪，贾海鹏，张云泉。基于 ARMV8 平台的多维 FFT 实现与优化研究
Qinglin Wang,Dongsheng Li. Optimizing FFT-Based Convolutionon ARMv8 Multi-core CPUs
Aleksandar Zlateski, Zhen Jia, Kai Li, Fredo Durand. FFT Convolutions are Faster than Winograd onModern CPUs, Here’s Why

附：

GitHub：MegEngine 天元

官网：MegEngine-深度学习，简单开发

欢迎加入 MegEngine 技术交流 QQ 群：1029741705

最后编辑于：2021.07.20 18:50:47

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 204,732评论 6赞 478
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 87,496评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 151,264评论 0赞 338
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,807评论 1赞 277
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,806评论 5赞 368
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,675评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 38,029评论 3赞 399
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,683评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 41,704评论 1赞 299
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,666评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,773评论 1赞 332
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,413评论 4赞 321
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 39,016评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,978评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,204评论 1赞 260
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 45,083评论 2赞 350
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,503评论 2赞 343