前言：

幾個月前，我根據(jù) Simoncelli 2016 年的論文編寫了自己的自動編碼器，用于研究目的。一開始，我想使用一些流行的深度學(xué)習(xí)框架（例如 Tensor Flow、Caffe2 或 MXNet）來做我的實驗。然而，在對所有這些框架進行了幾周的調(diào)查之后，我發(fā)現(xiàn)了一個非常令人頭疼的問題——可擴展性。我不是說這些框架設(shè)計得不好，而是不允許用戶開發(fā)第三方算子，就像寫一個插件一樣，你給我一個沒有任何參數(shù)的函數(shù)。那么改變函數(shù)行為的唯一方法就是修改源代碼，由于文檔組織不善，這無疑是一個巨大的工程。（這似乎是開源軟件的通病。）因此，由于不常見的算子 GDN 并未包含在所有這些框架中，因此設(shè)計一個新框架似乎是唯一的解決方案。

GDN

這個算子是這個理論中的核心非線性函數(shù)，表達式如下（公式不重要，如果你不喜歡這些該死的符號，你可以直接跳過這一節(jié)。）：

上標(biāo)(k)和(k+1)表示層數(shù)，w和u是多通道圖像的輸入和輸出，下標(biāo)i是通道數(shù)。β 和 γ 是我要訓(xùn)練的參數(shù)。假設(shè)我們有 N 個通道，那么 γ 是一個 N × N 矩陣，β 是一個 N × 1 向量。乍一看，這個功能與 cudnn 和所有深度學(xué)習(xí)框架都很好地支持的批量歸一化 (BN) 或局部響應(yīng)歸一化 (LRN) 非常相似。但相信我，不要讓你的眼睛欺騙你。這是非常不同的。（注意大除法是元素除法。）

前向不會消耗太多計算能力，而后向會消耗我 GPU 的大部分能量?，F(xiàn)在讓我們看看后面。我需要計算 3 個梯度，?β、?γ 和 ?u。

我知道人們第一次看到這個的感覺，因為我第一次看到這個怪物時也想自殺。 但如果我能為所有這些狗屎畫一幅畫，你會感覺更舒服。

首先，我們可以很容易地注意到輸入可以看作是一個長度為 m x n 的向量。其次，(blabla...)^(-3/2) 出現(xiàn)在所有這些梯度中。這意味著我們可以只計算該術(shù)語 1 次，并將它們緩存以備后用。我們稱其為“(blabla...)^(-1/2)”矩陣 D 。最后，δ 是傳播到前一層的誤差。

Fig 1. Computation of γ

經(jīng)過一些簡化，它更清楚了，對吧？ 我知道仍然需要一些解釋。 對于等式的右側(cè)，每個矩形都是由我們上面提到的矩陣堆疊而成的向量。 D 是 GDN 公式中的分母項，還記得我們剛剛提到的“(blabla...)^(-1/2)”嗎？

與一些高級算法不同，這種計算對大多數(shù)人來說非常直觀，我們可以輕松編寫 CPU 程序來處理它。只要稍微了解一下 CUDA，每個人都可以將他們的 CPU 代碼移植到 GPU。但是，如果您可以選擇不同的組織來啟動內(nèi)核，則速度會有很大的不同。

1. 不僅僅是天真的算法。

我稱這種方法“不只是天真”是因為這是我用過的第一種方法。即使使用小尺寸圖像作為輸入，它也幾乎耗盡了我所有的 GPU 內(nèi)存，并實現(xiàn)了最慢的性能。沒有利用任何內(nèi)存重用，我只是垂直和水平復(fù)制所有這些小矩形以獲得更大的矩陣，如下圖所示，并啟動許多一維組織的內(nèi)核。然后將它們相加。

Fig 2. Less than naive Algo.

該算法唯一的優(yōu)點是不需要在每個CUDA線程中計算索引，因為線程id只是唯一對應(yīng)的內(nèi)存索引。所以你需要做的就是一些乘法，然后使用 cublas 將每個小彩色矩形與 1 向量（一個充滿所有 1 的向量）的點積相加。但是正如你所看到的，矩形的大小并不像我這里畫的那么小，大小和圖像一樣。對于這張圖片中的每個向量，大小將為 N x N x imageSize x batchSize。很明顯，我們浪費了 (N-1) x N x imageSize x batchSize x 4 個字節(jié)，更不用說浪費在訪問所有這些冗余全局內(nèi)存上的時間了。

2. 樸素算法。

對于第一種算法，我每次迭代只能在我的網(wǎng)絡(luò)中訓(xùn)練不到 4 張大小為 128 x 128 的圖像，時間幾乎為 2 秒。（我的 GPU 是 GTX 1080。）這個現(xiàn)實迫使我改進我的算法，否則，我必須等待近 2 個月才能得到我的結(jié)果。

因為我需要啟動的內(nèi)核數(shù)量肯定比我GPU中的CUDA內(nèi)核多很多，所以不管我用什么方法，cuda驅(qū)動都會把這些任務(wù)序列化。然后我決定不復(fù)制所有這些記憶。相反，我將啟動 N x 一維組織的 N x imageSize 內(nèi)核 N 次（N 是通道總數(shù)）。

Fig 3. Without memory replication

可以看出，改進是顯而易見的。因為，我們不再需要大量復(fù)制數(shù)據(jù)。 GPU 中的全局內(nèi)存訪問非常昂貴。內(nèi)存訪問模式也很簡單，因為當(dāng)您獲得線程 id 時，只需使用一個 mod 操作就可以獲得內(nèi)存索引（內(nèi)存索引 = 線程 id % imageSize）。但是，在這種方法中，由于內(nèi)核仍然是一維組織的，并且我們使用for循環(huán)來啟動所有這些內(nèi)核，那么我們可能無法從GPU更智能的調(diào)度算法中受益，盡管我已經(jīng)嘗到了血的滋味.現(xiàn)在，通過這個小小的改變，2 個月的訓(xùn)練時間可以縮短到將近 2 周。

3. 更智能的組織算法。

到目前為止，我還沒有考慮過共享內(nèi)存的威力，因為對我來說，通常設(shè)計一個好的內(nèi)核模式是枯燥和頭痛的。顯然，一維內(nèi)核模式是最容易編寫的代碼。然而，更好的性能值得更仔細(xì)的設(shè)計。令我驚訝的是，本節(jié)中的算法實現(xiàn)了第二個算法的 3 倍速度。

回到圖 1，可以看到前 3 個右側(cè)矩陣的第一行 δ0、w0 和 D0 是相同的。因此，我們可以在一個塊中計算一行 γ，對于每個塊我們可以啟動 imageSize 個線程，并且對于每個線程我們可以使用 for 循環(huán)計算所有通道。

Fig 5. Computation in one block

所以從圖 5 來看，將 δ0、w0 和 D0 放在共享內(nèi)存中是非常直觀的，而對于線程 i，它從 0 到 N-1 讀取 N 個通道中的一個像素與 δ0、w0 和 D0 相乘 分享回憶。偽代碼如下：

blockId = blockIdx.x; 
threadId = threadIdx.x;shareDelta <- delta[blockId];  
shareW <- W[blockId];
shareD <- D[blockId];
_synchronize();for(i = 0; i < N-1; i++)
{
   result[threadIdx i*imgSize] = shareDelta[threadId] *
                                 shareW[threadId] *
                                 shareD[threadId] * 
                                 W[threadId + i*imgSize];
}

Algo 2 選擇行主計算而不是列主計算是因為在一個網(wǎng)格中計算一行，我們可以共享 3 個向量 δ0、w0 和 D0。但是如果我們像在 Algo 中那樣計算一列，我們只能共享 1 個向量 w0。（再次參見圖 1。）。

在這段代碼片段中，沒有 if ... else ... 塊。這在并行計算中非常重要。因為所有線程都是并行運行的，理想的情況是所有這些線程同時完成它們的工作。但是如果有 if ... else ... 阻塞，分支會讓這些線程做不同的任務(wù)，以便它們在不同的時間完成。然后計算時間將由最慢的線程決定。

無索引計算也是一個優(yōu)勢。通過設(shè)計一維模式，我們必須使用線程id來計算內(nèi)存索引，但這里不需要將blockId和threadId轉(zhuǎn)換為一維內(nèi)存索引來訪問數(shù)據(jù)。

最后，因為我的數(shù)據(jù)存儲在列major中，這意味著，像向量δ0一樣，這個向量中的所有元素都是連續(xù)存儲的。所以它受益于全局內(nèi)存合并機制。全局內(nèi)存也是cuda中的一個重要概念。

在硬件方面，16個cuda內(nèi)核被組織在一個warp中。當(dāng)其中一個線程訪問數(shù)據(jù)時，例如上圖中的 a1，數(shù)據(jù)總線不僅會傳輸 a1，還會將 a1~a32 傳輸?shù)骄彺嬷校约铀倨渌?15 個內(nèi)核的數(shù)據(jù)訪問。因此，當(dāng)我讀取全局?jǐn)?shù)據(jù)以共享內(nèi)存時，每 32 個字節(jié)我只讀取一次，所有其他字節(jié)都從緩存中讀取，速度快了數(shù)百。多虧了時空局域性理論。

4. 多一點改進

今天突然發(fā)現(xiàn)其實我不需要共享內(nèi)存，但是可以使用const內(nèi)存。因為對于向量δ0、w0和D0，一個block中的每個線程只需要訪問一次。所以在for循環(huán)之前，我們實際上可以將元素緩存在const內(nèi)存中。另一個糖是因為每個線程只訪問一個元素，不需要線程同步。

代碼如下：

blockId = blockIdx.x; 
threadId = threadIdx.x;const float constDelta = delta[blockId * imgSize + threadId];  
const float constW = W[blockId * imgSize + threadId];
const float constD = D[blockId * imgSize + threadId];for(i = 0; i < N-1; i++)
{
   result[threadIdx + i*imgSize] = constDelta * constW *
                                   constD * 
                                   W[threadId + i*imgSize];
}

從上面的代碼可以看出，constDelta、constW、constD可以從本地內(nèi)存中重復(fù)使用N次，本地內(nèi)存總是存儲在本地寄存器中。因此，帶寬大于共享內(nèi)存。

Reduce Operation

我講的所有算法都沒有完成，因為我從上述算法中得到的實際上都是原始γ，如下所示：

我需要在左側(cè)累積每個向量以獲得一個元素。第一個選擇是 cublas API，cublasSsbmv。此函數(shù)將進行矩陣向量乘法。所以我們可以把左邊的向量看成一個矩陣，將它與一個全1向量相乘，得到γ的一行梯度。并重復(fù)N次以獲得最終結(jié)果。但我注意到還有其他 API cublasSgemmBatched。此函數(shù)可以進行批量矩陣向量乘法。然后我做了一個實驗來測試哪個更快：

N 個矩陣向量乘法 VS 批處理矩陣向量乘法的 for 循環(huán)。

結(jié)果表明for循環(huán)要快得多。但是我不知道原因，也許是因為我這里的 N 太?。?/span>N = 256）。

我不會展示如何計算 ?β 和 ?u，因為它們類似于 ?γ。我知道必須有比我更進一步的優(yōu)化或更好的設(shè)計。CUDA 優(yōu)化對于不深入了解 GPU 組織的人來說通常是困難的。熟悉 CPU 的程序員總是受益于現(xiàn)代操作系統(tǒng)和強大的編譯器。然而，GPU 在編寫足夠的代碼方面與 CPU 有很大不同和復(fù)雜性，盡管它比以前使用圖形著色器進行計算要方便得多。生態(tài)環(huán)境的完善還需要幾年時間。

原文鏈接：

https://medium.com/@Lawliet0320/ramble-in-cuda-optimization-8fbbcf81e7c5

整理不易，點贊三連↓