深度學習算子優(yōu)化-FFT

作者：嚴健文 | 曠視 MegEngine 架構師

背景

在數(shù)字信號和數(shù)字圖像領域揩晴，對頻域的研究是一個重要分支筐钟。
我們日常“加工”的圖像都是像素級疹鳄，被稱為是圖像的空域數(shù)據(jù)拧略。空域數(shù)據(jù)表征我們“可讀”的細節(jié)瘪弓。如果我們將同一張圖像視為信號辑鲤，進行頻譜分析，可以得到圖像的頻域數(shù)據(jù)杠茬。觀察下面這組圖 (來源)，頻域圖中的亮點為低頻信號弛随，代表圖像的大部分能量瓢喉，也就是圖像的主體信息。暗點為高頻信號舀透，代表圖像的邊緣和噪聲栓票。從組圖可以看出，Degraded Goofy 與 Goofy 相比愕够，近似的低頻信號保留住了 Goofy 的“輪廓”走贪，而其高頻信號的增加使得背景噪點更加明顯。頻域分析使我們可以了解圖像的組成惑芭，進而做更多的抽象分析和細節(jié)處理坠狡。

Goofy and Degraded Goofy

實現(xiàn)圖像空域和頻域轉換的工具，就是傅立葉變換遂跟。由于圖像數(shù)據(jù)在空間上是離散的逃沿，我們使用傅立葉變換的離散形式 DFT（Discrete Fourier Transform）及其逆變換 IDFT（Inverse Discrete Fourier Transform)婴渡。Cooley-Tuckey 在 DFT 的基礎上，開發(fā)了更快的算法 FFT（Fast Fourier Transform）凯亮。

DFT/FFT 在數(shù)字圖像領域還有一些延伸應用边臼。比如基于 DFT 的 DCT（Discrete Cosine Transform，離散余弦變換）就用在了圖像壓縮 JPEG 算法 (來源) 和圖像水印算法（來源）假消。

JPEG 編碼是通過色彩空間轉換柠并、抽樣分塊、DCT 變換富拗、量化編碼實現(xiàn)的态贤。其中 DCT 變換的使用將圖像低頻信息和高頻信息區(qū)分開，在量化編碼過程中壓縮了少量低頻信息蜕青、大量高頻信息從而獲得尺寸上壓縮钧敞。從貓臉圖上可看出隨著壓縮比增大畫質會變差，但是主體信息還是得以保留谅阿。

貓臉不同 jpeg 畫質（壓縮比）

圖像水印算法通過 DCT 將原圖轉換至頻域半哟，選取合適的位置嵌入水印圖像信息，并通過 IDCT 轉換回原圖签餐。這樣對原圖像的改變較小不易察覺寓涨，且水印通過操作可以被提取。

DFT/FFT 在深度學習領域也有延伸應用氯檐。比如利用 FFT 可以降低卷積計算量的特點戒良，F(xiàn)FT_Conv 算法也成為常見的深度學習卷積算法。本文我們就來探究一下頻域算法的原理和優(yōu)化策略冠摄。

DFT 的原理及優(yōu)化

公式

無論是多維的 DFT 運算糯崎，還是有基于 DFT 的 DCT/FFT_Conv, 底層的計算單元都是 DFT_1D。因此河泳，DFT_1D 的優(yōu)化是整個 FFT 類算子優(yōu)化的基礎沃呢。
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 份惰瓜，根據(jù) k * n 的值逆時針掃過復平面的圓周。

單位復根有著周期性和對稱性汉矿，我們依據(jù)這兩個性質可以對 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 輸出序列第股。
根據(jù) 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}$

根據(jù)奇偶項拆開，分成兩個長度為 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朵耕，執(zhí)行 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 不斷增加墓猎，蝶形群在“變大變密”，然而每一級總的碟形次數(shù)是不變的赚楚。
實際上毙沾，對于長度為 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)勇皇，顯得高效而優(yōu)雅罩句。我們基于蝶形算法，對于不同的 radix 進行算法的進一步劃分和優(yōu)化敛摘，主要分為 radix - 2 的冪次的和 radix – 非 2 的冪次兩類门烂。

radix-2 的冪次優(yōu)化

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]$

我們先來看訪存，現(xiàn)代處理器對于計算性能的優(yōu)化要優(yōu)于對于訪存的優(yōu)化泄私，在計算和訪存相近的場景下房揭，訪存通常是性能瓶頸。

DFT1D 中晌端，對于不同基底的算法 r-2/r-4/r-8, 每一個 stage 有著相等的存取量：2 * butterfly_num * radix = 2N, 而不同的基底對應的 stage 數(shù)有著明顯差異（ $\log_2N$ vs $\log_4N$ vs $\log_8N$ )捅暴。

因此對于 DFT, 在不顯著增加計算量的條件下，選用較大的 kernel 會在訪存上取得明顯的優(yōu)勢咧纠。觀察推導的 kernel 圖蓬痒， r-2 的 kernel 每個蝶形對應 4 次訪存操作和，2 次復數(shù)浮點加減運算漆羔。r-4 的 kernel 每個蝶形算法對應 8 次 load/store梧奢、8 次復數(shù)浮點加減操作（合并相同的運算）狱掂，在計算量略增加的同時 stage 由 $\log_2N$ 下降到 $\log_4N$ , 降低了總訪存的次數(shù)，因此會有性能的提升亲轨。r-8 的 kernel 每個蝶形對應 16 次 load/store趋惨、24 次復數(shù)浮點加法和 8 次浮點乘法。浮點乘法的存在使得計算代價有所上升瓶埋， stage 由 $\log_4N$ 進一步下降到 $\log_8N$ 希柿，但由于 N 日常并不會太大， r-4 到 r-8 的 stage 減少不算明顯养筒，所以優(yōu)化有限

我們再來看計算的開銷曾撤。減少計算的開銷通常有兩種辦法：減少多余的運算、并行化晕粪。

以 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 的數(shù)據(jù)由于 k=0, 旋轉因子都為 1，可以省去這部分復數(shù)乘法運算巫湘，單獨優(yōu)化装悲。 radix4_other_stage 部分，從第 2 個 stage 往后尚氛， butterfly_num = 4^(s-1) 都為 4 的倍數(shù)诀诊，而每個 butterfly 數(shù)組讀取/存儲都是間隔的≡乃唬可以對最里層的循環(huán)做循環(huán)展開加向量化属瓣，實現(xiàn) 4 個或更多 butterfly 并行運算。循環(huán)展開和 SIMD 指令的使用不僅可以提高并行性讯柔，也可以提升 cacheline 利用的效率抡蛙，可以帶來較大的性能提升。以 SM8150(armv8) 為例魂迄，r-4 的并行優(yōu)化可以達到 r2 的 1.6x 的性能粗截。

尺寸：1 * 2048（r2c) 環(huán)境：SM8150 大核

總之，對于 radix-2^n 的優(yōu)化捣炬，選用合適的 radix 以減少多 stage 帶來的訪存開銷熊昌，并且利用單位復根性質以及并行化降低計算的開銷，可以帶來較大的性能提升湿酸。

radix-非 2 的冪次優(yōu)化

當輸入長度 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 的優(yōu)化也是必要的窝撵。

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}$ 在復平面上根據(jù) x 軸對稱寒砖，有相同的實部和相反的虛部赐劣。根據(jù)這個性質。如下圖所示哩都，對于每一個 stage魁兼，可以合并公共項 A，B漠嵌，C咐汞，D，再根據(jù)公共項計算出該 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 做循環(huán)展開加并行化植阴，進一步減少計算的開銷。
radix-5 的優(yōu)化思想可以外推至 radix-N章钾。對于 radix_N 的每一個 stage, 計算流程為：

取輸入
乘以對應的轉換因子
計算公共項墙贱， radix_N 有 N-1 個公共項
執(zhí)行并行化的 radix_N_kernel
寫入輸出

其他優(yōu)化

上述兩個章節(jié)描述的是 DFT_1D 的通用優(yōu)化，在此基礎上還可以做更細致的優(yōu)化贱傀，可以參考本文引用的論文惨撇。

對于全實數(shù)輸入的，由于輸入的虛部為 0府寒，進行旋轉因子以及 radix_N_kernel 的復數(shù)運算時會有多余的運算和多余的存儲魁衙，可以利用 split r2c 算法，視為長度為 N/2 的復數(shù)序列株搔，計算 DFT 結果并進行 split 操作得到 N 長實數(shù)序列的結果剖淀。
對于 radix-2 的冪次算法，重新計算每個 stage 的輸入/輸出 stride 以取消第一級的位元翻轉可以進一步減少訪存的開銷纤房。
對于 radix-N 算法纵隔，在混合基框架下 N = radix1^m1 * radix2^m2，合并較小的 radix 為大的 radix 以減少 stage。

DFT 延展算法的原理及優(yōu)化

DCT 和 FFT_conv 兩個典型的基于 DFT 延展的算法捌刮，DFT_1D/2D 的優(yōu)化可以很好的用在這類算法中碰煌。

DCT

DCT 算法（Discrete Cosine Transform，離散余弦變換）可以看作是 DFT 取其正弦分量并經過工業(yè)校正的算法绅作。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 實現(xiàn)是 O(n^2) 的芦圾，而我們將其轉換成 DFT_1D 算法，可以將算法復雜度降至 O(nlogn)俄认。
基于 DFT 的 DCT 算法流程為：

對于 DCT 的輸入序列 x[n], 創(chuàng)建長為 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}}$ , 根據(jù)歐拉公式和誘導函數(shù)岂贩，整理非公共項 $(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 有關的系數(shù)的乘積科平。這樣就可以通過先計算 $Y[k]$ 得到 x[n] 的 DCT 輸出 $X[k]$ 。

在理解算法的基礎上姜性，我們對 DFT_1D 的優(yōu)化可以完整地應用到 DCT 上瞪慧。DCT_2D 的計算過程是依次對行、列做 DCT_1D, 我們用多線程對 DCT_1D 進行并行部念，可以進一步優(yōu)化算法弃酌。

FFT_conv

Conv 是深度學習最常見的運算，計算 conv 常用的方法有 IMG2COL+GEMM儡炼， Winograd妓湘， FFT_conv。三種算法都有各自的使用場景乌询。

FFT_conv 的數(shù)學原理是時域中的循環(huán)卷積對應于其離散傅里葉變換的乘積榜贴。如下圖所示， 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]$

將卷積公式和離散傅立葉變換展開鬼佣，改變積分的順序并且替換變量驶拱，可以證明結論。
注意這里的卷積是循環(huán)卷積晶衷，和我們深度學習中常用的線性卷積是有區(qū)別的蓝纲。利用循環(huán)卷積計算線性卷積的條件為循環(huán)卷積長度 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 到的尺寸較小蝎抽，可以大幅減小這一部分的開銷路克。優(yōu)化后 fft_conv 的計算流程為：

合理安排緩存計算出合適的 tile 尺寸衷戈，對原圖進行分塊
分塊后的小圖和 kernel 進行 zero-padding, 并進行 DFT 運算
小圖矩陣點乘
進行逆運算并組合成大圖。

同時我們可以觀察到刁笙，F(xiàn)FT_conv 的核心計算模塊還是針對小圖的 DFT 運算，因此我們可以將前一章節(jié)對 DFT 的優(yōu)化代入此處座每，輔以多線程摘悴，進一步提升 FFT_Conv 的計算效率。

參考資料

陳暾蹂喻，李志豪口四，賈海鵬，張云泉治笨。基于 ARMV8 平臺的多維 FFT 實現(xiàn)與優(yōu)化研究
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-深度學習旷赖，簡單開發(fā)

歡迎加入 MegEngine 技術交流 QQ 群：1029741705

最后編輯于：2021.07.20 18:50:47

?著作權歸作者所有,轉載或內容合作請聯(lián)系作者

人面猴
序言：七十年代末杠愧，一起剝皮案震驚了整個濱河市逞壁，隨后出現(xiàn)的幾起案子腌闯，更是在濱河造成了極大的恐慌，老刑警劉巖糖声，帶你破解...
沈念sama閱讀 216,544評論 6贊 501
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件蘸泻，死亡現(xiàn)場離奇詭異嘲玫，居然都是意外死亡去团，警方通過查閱死者的電腦和手機穷蛹，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 92,430評論 3贊 392
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門肴熏，熙熙樓的掌柜王于貴愁眉苦臉地迎上來顷窒，“玉大人鞋吉，你說我怎么就攤上這事。” “怎么了漆魔？”我有些...
開封第一講書人閱讀 162,764評論 0贊 353
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵改抡，是天一觀的道長系瓢。經常有香客問我夷陋，道長，這世上最難降的妖魔是什么藐窄？我笑而不...
開封第一講書人閱讀 58,193評論 1贊 292
?港島之戀（遺憾婚禮）
正文為了忘掉前任荆忍，我火速辦了婚禮刹枉，結果婚禮上屈呕，老公的妹妹穿的比我還像新娘。我一直安慰自己芥吟，他們只是感情好钟鸵，可當我...
茶點故事閱讀 67,216評論 6贊 388
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布棺耍。她就那樣靜靜地躺著，像睡著了一般俊卤。火紅的嫁衣襯著肌膚如雪害幅。梳的紋絲不亂的頭發(fā)上以现，一...
開封第一講書人閱讀 51,182評論 1贊 299
城市分裂傳說
那天邑遏，我揣著相機與錄音，去河邊找鬼憎蛤。笑死俩檬，一個胖子當著我的面吹牛碾盟，可吹牛的內容都是我干的。我是一名探鬼主播晚胡，決...
沈念sama閱讀 40,063評論 3贊 418
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼估盘，長吁一口氣：“原來是場噩夢啊……” “哼遣妥！你這毒婦竟也來了攀细？” 一聲冷哼從身側響起，我...
開封第一講書人閱讀 38,917評論 0贊 274
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤锦担，失蹤者是張志新（化名）和其女友劉穎慨削，沒想到半個月后缚态，有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體玫芦，經...
沈念sama閱讀 45,329評論 1贊 310
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡桥帆，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內容為張勛視角年9月15日...
茶點故事閱讀 37,543評論 2贊 332
?白月光啟示錄
正文我和宋清朗相戀三年环葵，在試婚紗的時候發(fā)現(xiàn)自己被綠了宝冕。大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片地梨。...
茶點故事閱讀 39,722評論 1贊 348
活死人
序言：一個原本活蹦亂跳的男人離奇死亡宝剖，死狀恐怖，靈堂內的尸體忽然破棺而出扑眉，到底是詐尸還是另有隱情腰素，我是刑警寧澤雪营，帶...
沈念sama閱讀 35,425評論 5贊 343
?日本核電站爆炸內幕
正文年R本政府宣布献起，位于F島的核電站，受9級特大地震影響姻政，放射性物質發(fā)生泄漏扶歪。R本人自食惡果不足惜，卻給世界環(huán)境...
茶點故事閱讀 41,019評論 3贊 326
男人毒藥：我在死后第九天來索命
文/蒙蒙一妹萨、第九天我趴在偏房一處隱蔽的房頂上張望乎完。院中可真熱鬧树姨，春花似錦桥状、人聲如沸辅斟。這莊子的主人今日做“春日...
開封第一講書人閱讀 31,671評論 0贊 22
一樁弒父案查邢，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽扰藕。三九已至芳撒，卻和暖如春笔刹，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背门躯。一陣腳步聲響...
開封第一講書人閱讀 32,825評論 1贊 269
情欲美人皮
我被黑心中介騙來泰國打工讶凉，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留懂讯，地道東北人。一個月前我還...
沈念sama閱讀 47,729評論 2贊 368
代替公主和親
正文我出身青樓勒庄，卻偏偏與公主長得像实蔽，于是被迫代替她去往敵國和親谨读。傳聞我的和親對象是個殘疾皇子劳殖，可洞房花燭夜當晚...
茶點故事閱讀 44,614評論 2贊 353