漫談CUDA優化

CV技術指南（公眾號）發表於2021-08-03

作者：Lawliet

翻譯：彷彿若有光

前言：

幾個月前，我根據 Simoncelli 2016 年的論文編寫了自己的自動編碼器，用於研究目的。一開始，我想使用一些流行的深度學習框架（例如 Tensor Flow、Caffe2 或 MXNet）來做我的實驗。然而，在對所有這些框架進行了幾周的調查之後，我發現了一個非常令人頭疼的問題——可擴充套件性。我不是說這些框架設計得不好，而是不允許使用者開發第三方運算元，就像寫一個外掛一樣，你給我一個沒有任何引數的函式。那麼改變函式行為的唯一方法就是修改原始碼，由於文件組織不善，這無疑是一個巨大的工程。（這似乎是開源軟體的通病。）因此，由於不常見的運算元 GDN 並未包含在所有這些框架中，因此設計一個新框架似乎是唯一的解決方案。

點個關注，專注於計算機視覺的技術總結和分享

GDN

這個運算元是這個理論中的核心非線性函式，表示式如下（公式不重要，如果你不喜歡這些該死的符號，你可以直接跳過這一節。）：

上標(k)和(k+1)表示層數，w和u是多通道影像的輸入和輸出，下標i是通道數。β 和 γ 是我要訓練的引數。假設我們有 N 個通道，那麼 γ 是一個 N × N 矩陣，β 是一個 N × 1 向量。乍一看，這個功能與 cudnn 和所有深度學習框架都很好地支援的批量歸一化 (BN) 或區域性響應歸一化 (LRN) 非常相似。但相信我，不要讓你的眼睛欺騙你。這是非常不同的。（注意大除法是元素除法。）

前向不會消耗太多計算能力，而後向會消耗我 GPU 的大部分能量。現在讓我們看看後面。我需要計算 3 個梯度，∇β、∇γ 和 ∇u。

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

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

Fig 1. Computation of γ

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

與一些高階演算法不同，這種計算對大多數人來說非常直觀，我們可以輕鬆編寫 CPU 程式來處理它。只要稍微瞭解一下 CUDA，每個人都可以將他們的 CPU 程式碼移植到 GPU。但是，如果您可以選擇不同的組織來啟動核心，則速度會有很大的不同。

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

我稱這種方法“不只是天真”是因為這是我用過的第一種方法。即使使用小尺寸影像作為輸入，它也幾乎耗盡了我所有的 GPU 記憶體，並實現了最慢的效能。沒有利用任何記憶體重用，我只是垂直和水平復制所有這些小矩形以獲得更大的矩陣，如下圖所示，並啟動許多一維組織的核心。然後將它們相加。

Fig 2. Less than naive Algo.

該演算法唯一的優點是不需要在每個CUDA執行緒中計算索引，因為執行緒id只是唯一對應的記憶體索引。所以你需要做的就是一些乘法，然後使用 cublas 將每個小彩色矩形與 1 向量（一個充滿所有 1 的向量）的點積相加。但是正如你所看到的，矩形的大小並不像我這裡畫的那麼小，大小和影像一樣。對於這張圖片中的每個向量，大小將為 N x N x imageSize x batchSize。很明顯，我們浪費了 (N-1) x N x imageSize x batchSize x 4 個位元組，更不用說浪費在訪問所有這些冗餘全域性記憶體上的時間了。

2. 樸素演算法。

對於第一種演算法，我每次迭代只能在我的網路中訓練不到 4 張大小為 128 x 128 的影像，時間幾乎為 2 秒。（我的 GPU 是 GTX 1080。）這個現實迫使我改進我的演算法，否則，我必須等待近 2 個月才能得到我的結果。

因為我需要啟動的核心數量肯定比我GPU中的CUDA核心多很多，所以不管我用什麼方法，cuda驅動都會把這些任務序列化。然後我決定不復制所有這些記憶。相反，我將啟動 N x 一維組織的 N x imageSize 核心 N 次（N 是通道總數）。

Fig 3. Without memory replication

可以看出，改進是顯而易見的。因為，我們不再需要大量複製資料。 GPU 中的全域性記憶體訪問非常昂貴。記憶體訪問模式也很簡單，因為當您獲得執行緒 id 時，只需使用一個 mod 操作就可以獲得記憶體索引（記憶體索引 = 執行緒 id % imageSize）。但是，在這種方法中，由於核心仍然是一維組織的，並且我們使用for迴圈來啟動所有這些核心，那麼我們可能無法從GPU更智慧的排程演算法中受益，儘管我已經嚐到了血的滋味.現在，通過這個小小的改變，2 個月的訓練時間可以縮短到將近 2 周。

3. 更智慧的組織演算法。

到目前為止，我還沒有考慮過共享記憶體的威力，因為對我來說，通常設計一個好的核心模式是枯燥和頭痛的。顯然，一維核心模式是最容易編寫的程式碼。然而，更好的效能值得更仔細的設計。令我驚訝的是，本節中的演算法實現了第二個演算法的 3 倍速度。

回到圖 1，可以看到前 3 個右側矩陣的第一行 δ0、w0 和 D0 是相同的。因此，我們可以在一個塊中計算一行 γ，對於每個塊我們可以啟動 imageSize 個執行緒，並且對於每個執行緒我們可以使用 for 迴圈計算所有通道。

Fig 5. Computation in one block

所以從圖 5 來看，將 δ0、w0 和 D0 放在共享記憶體中是非常直觀的，而對於執行緒 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 選擇行主計算而不是列主計算是因為在一個網格中計算一行，我們可以共享 3 個向量 δ0、w0 和 D0。但是如果我們像在 Algo 中那樣計算一列，我們只能共享 1 個向量 w0。（再次參見圖 1。）。

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

無索引計算也是一個優勢。通過設計一維模式，我們必須使用執行緒id來計算記憶體索引，但這裡不需要將blockId和threadId轉換為一維記憶體索引來訪問資料。

最後，因為我的資料儲存在列major中，這意味著，像向量δ0一樣，這個向量中的所有元素都是連續儲存的。所以它受益於全域性記憶體合併機制。全域性記憶體也是cuda中的一個重要概念。

在硬體方面，16個cuda核心被組織在一個warp中。當其中一個執行緒訪問資料時，例如上圖中的 a1，資料匯流排不僅會傳輸 a1，還會將 a1~a32 傳輸到快取中，以加速其他 15 個核心的資料訪問。因此，當我讀取全域性資料以共享記憶體時，每 32 個位元組我只讀取一次，所有其他位元組都從快取中讀取，速度快了數百。多虧了時空局域性理論。

4. 多一點改進

今天突然發現其實我不需要共享記憶體，但是可以使用const記憶體。因為對於向量δ0、w0和D0，一個block中的每個執行緒只需要訪問一次。所以在for迴圈之前，我們實際上可以將元素快取在const記憶體中。另一個糖是因為每個執行緒只訪問一個元素，不需要執行緒同步。

程式碼如下：

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次，本地記憶體總是儲存在本地暫存器中。因此，頻寬大於共享記憶體。

Reduce Operation

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

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

N 個矩陣向量乘法 VS 批處理矩陣向量乘法的 for 迴圈。

結果表明for迴圈要快得多。但是我不知道原因，也許是因為我這裡的 N 太小（N = 256）。

我不會展示如何計算 ∇β 和 ∇u，因為它們類似於 ∇γ。我知道必須有比我更進一步的優化或更好的設計。CUDA 優化對於不深入瞭解 GPU 組織的人來說通常是困難的。熟悉 CPU 的程式設計師總是受益於現代作業系統和強大的編譯器。然而，GPU 在編寫足夠的程式碼方面與 CPU 有很大不同和複雜性，儘管它比以前使用圖形著色器進行計算要方便得多。生態環境的完善還需要幾年時間。

原文連結：

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

本文來源於公眾號 CV技術指南的論文分享系列。

歡迎關注公眾號 CV技術指南，專注於計算機視覺的技術總結、最新技術跟蹤、經典論文解讀。

在公眾號中回覆關鍵字 “技術總結” 可獲取以下文章的彙總pdf。