[轉]矩陣運算優化問題

form https://www.zhihu.com/question/19706331

(一) MKL庫問題

作者：過擬合
鏈接：https://www.zhihu.com/question/19706331/answer/25255444
來源：知乎
著作權歸作者所有，轉載請聯系作者獲得授權。

MATLAB的矩陣計算使用的是Intel自己出的Math kernel library（MKL）,這個庫遠比其他的blas/lapack庫要快。C快在循環，要想矩陣計算也和MATLAB一樣快，那就得鏈接MKL，寫起來免不了各種折騰。而且，即使你鏈接上了，編譯時各種優化選項之類的還是比不上人家專業的設定，速度很難接近MATLAB。
我自己在Gentoo上試過源里的所有blas/lapack庫，無一能與MKL匹敵，而且連接近都不可能。甚至我把python的NumPy庫鏈接上MKL后，速度也只是勉強接近。由於Gentoo的MKL庫永遠是最新的，而每一個新版本的MKL庫對矩陣計算都有略微提升，導致可能暫時NumPy與MATLAB可以匹敵。但是一旦更新版本的MATLAB出來后，它會使用上更新的MKL庫，這種領先優勢就又喪失殆盡。你可以在MATLAB文檔搜索中輸入MKL，這樣會被定位到MATLAB release notes，而里面就會含有這么一句話“Upgrade to Intel Math Kernel Libraries”，這就是每一個版本MATLAB矩陣計算都越發變態快的原因。
當然，剛我提到的python，其矩陣計算速度雖然微微落后於MATLAB，但是在很多其他地方是可以大大強於MATLAB的。例如繪制大規模三維點雲，以及輕松調用gpu之類的。因此python在矩陣計算的微小速度劣勢完全可以忽略，可以考慮用於科學計算。關於NumPy鏈接MKL，我之前寫過一篇博文： Numpy使用MKL庫提升計算性能

 

（二）MATLAB的優化

作者：王洋子豪
鏈接：https://www.zhihu.com/question/19706331/answer/16134037
來源：知乎
著作權歸作者所有，轉載請聯系作者獲得授權。

矩陣乘法是一個相對成熟的問題，根據矩陣的稀疏程度有不同的優化算法。
不使用GPU加速的MATLAB版本采用的是BLAS中的General Matrix Multiplication[1]。學術界有各種矩陣乘法算法將其復雜度降低到O(n^2.x)，例如Strassen和Winograd算法，在BLAS中應該已經使用了Strassen算法。
如果你的MATLAB是安裝了Parallel Computing Toolbox的話，那么很可能它已經在使用GPU進行計算了。這種情況下采用的是MAGMA[2]。我沒有使用過MAGMA，但我猜測它應該使用了cuBLAS來計算矩陣乘法。
宏觀角度上對矩陣乘法的優化包括對局部內存使用的優化(Blocked/Tiled)以及對中間運算步驟的優化(Strassen/Winograd)，實現細節上的優化就非常繁多了。比如loop unrolling，多級的tiling，指令級並行等等。其中會牽扯到一些編譯器和體系結構的知識，似乎對僅僅希望使用矩陣乘法函數的用戶來講沒有什么太大必要去探究。所以我認為寫出比MATLAB更快的矩陣乘法是可行但困難的。
[1] General Matrix Multiply
[2] MAGMA

 

（三）矩陣乘法cache優化

作者：楊樹下的狐狸
鏈接：https://www.zhihu.com/question/19706331/answer/15392676
來源：知乎
著作權歸作者所有，轉載請聯系作者獲得授權。

偏個題。
說到這個想起之前朋友和我說到他最近在上一個課。
那個課上，教授要求他們寫Cache friendly code。
尤其像矩陣這種很大塊的東西，在運算時，會導致cache根本無法完全裝下需要使用的數據。
此刻，如果程序沒有設計得很有技巧，不斷地刷新cache，會需要浪費大量的時間。
所以，他們教了好幾種方法去計算矩陣，讓整個計算過程中盡量減少cache的重新載入：
以下是引用朋友給我的郵件，作者是 @Tian Tan :

先給你說個好玩的。這是我上的一門課的內容，叫high performance computing

這周在超級計算機上做了一個實驗，
實驗內容是想盡辦法優化很小一段代碼，比如矩陣乘法。

先說說cache的特點。
當訪問內存中一個element的時候，
cpu會把那個element放到cache里面，
同時還會把它臨近的elements放進cache。

衡量CPU快慢的一個標准是MFLOPS, 全稱為millions of floating point operations per second. 實驗的宗旨是寫cache friendly code. 直接上例子吧。

A, B, C都是浮點數矩陣
for (i=0; i<N; i++)
for (j=0; j<N; j++)
for (k=0; k<N; k++)
C[i][j] += A[i][k]*B[k][j];

這是個很簡單的矩陣乘法算法。但是這么寫效率是不高的，
原因在於當N很大的時候，比如2048， 4096，會產生cache conflict。要解釋這個術語比較麻煩，
想知道的話去看computer system: a programmer's perspective那本書的第六章。

但是矩陣乘法有別的寫法。比如

for (i=0; i<N; i++)
for (k=0; k<N; k++)
for (j=0; j<N; j++)
C[i][j] += A[i][k]*B[k][j];

這種寫法交換了k和j的位置，效率應該會比前面那個高些。（
術語是loop permutation）

還有的寫法叫loop tiling，
tiling的實質是將大矩陣的乘法變成小的分塊矩陣的乘法。

就用上一個例子吧。
for(it=0; it<N; it+=T)
for(jt=0; jt<N; jt+=T)
for(kt=0; kt<N; kt+=T)
for(i=it; i<it+T; i++)
for(j=jt; j<jt+T; j++)
for(k=kt; k<kt+T; k++)
C[i][j] += A[i][k]*B[k][j];

其中的T叫tiling size，能整除N。

這樣先算小矩陣的話，cache 就能裝下參與運算的elements, 對速度提升很大。

在實驗中有一道題，
經過優化之后把運行時間從49秒降到13秒了。

矩陣乘法只是最簡單例子，不同的code優化方式各異，
但是基本思想一樣。

另外，正如上面說的。
使用針對你自己的CPU的編譯器，編譯器有可能能夠識別到你的功能，做出相應的優化。