如何寫出比 MATLAB 更快的矩陣運算程序？

11-24

MATLAB 似乎將矩陣運算優化得非常之好，比手寫的 C 程序運算 2D 數組還快。想知道它是怎麼做到的，有沒有辦法能比它做得更好。

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庫提升計算性能。

偏個題。
說到這個想起之前朋友和我說到他最近在上一個課。
那個課上，教授要求他們寫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& for (j=0; j& for (k=0; 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& for (k=0; k& for (j=0; j& C[i][j] += A[i][k]*B[k][j];

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

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

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

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

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

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

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

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

矩陣乘法是一個相對成熟的問題，根據矩陣的稀疏程度有不同的優化演算法。
不使用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

在多核的情況下優化稀疏矩陣的運算速度可以從下面幾個方面入手：
1. 更好的partition，每個core計算一部分的矩陣運算，減少false sharing的情況，
2. 設計內存中的數據組織方式減少cache miss和TLB miss
3. 更好的load-balance，使得在計算的時候沒有閑著的core
4. 如果是分散式的，那麼還要考慮盡量減少通信的開銷

另外，一般稀疏矩陣的運算可以抽象成為一個圖，對此，有很多圖計算的框架可以使用
1. GraphChi，CMU的一幫搞ML的人寫的單機圖計算框架
使用場景：圖比較大，無法都放在內存中，計算的速度要比Hadoop塊上很多倍，甚至和多個節點的Spark性能差不多

具體可以看發表在OSDI12上的文章 http://select.cs.cmu.edu/publications/paperdir/osdi2012-kyrola-blelloch-guestrin.pdf
2. GraphLab，CMU的一幫搞ML的人和搞OS的人寫的，分散式圖計算框架
使用場景：圖比較大，單個機器的內存放不下，那我就加機器，做分散式共享內存來計算，速度甩出Hadoop好幾條街，性能接近MPI，具體可以看發表在OSDI12上的 http://arxiv.org/abs/1006.4990 同時發表在OSDI12上針對power-law的圖優化的PowerGraph，合併到GraphLab的v2.2上了。

另外：這兩個框架都實現了專門針對機器學習的常用的庫，如果LZ是做機器學習的，可以參考。

如果數據量比較大的情況下，可能matlab性能就比較弱了，可以用上面的兩個框架

兩個框架的tutorial可以看官方網站：GraphLab.org
後來GraphChi和GraphLab的作者開公司了，GraphLab | Unleash Data Science

=======
補充：
很多同學都提到了用GPU和FPGA之類的做計算，和他們比起來，
1. GraphChi只需要一台普通的PC，就行，對硬體要求不大，編程模型也比較簡單，就是C++的代碼，不需要安裝任何庫
2. GraphLab需要幾台PC和一個交換機，需要安裝MPI，C++代碼，編程模型也比較簡單。

矩陣運算在圖像處理領域經常出現，遇到這類反覆迭代的運算，首先想到的是利用SIMD指令去做運算並行化，比如sse系列。配合前面仁兄說的cache命中手段，程序性能可以提高100%到300%左右。哦，多核機器可以考慮多線程並發處理，手機打字累

不要輕視Matlab的優化。

Matlab雖然界面上用的JAVA JVM給人一種很卡的感覺，但內部的矩陣運算全部都是用針對特定CPU在彙編級別優化過的矩陣運算庫實現的（如BLAS）。因此，如果你想寫出比它還要快的程序，這個問題要考慮這樣兩個方面：

（1）你用來對比的Matlab程序是優化的代碼嗎？雖然Matlab發展到最新的幾個版本已經能夠對循環做自動的優化，但畢竟還很依賴程序員是怎樣實現演算法的。譬如循環內外的矩陣定義，初始化等等都會影響到自動優化的過程。

（2）其實使用C/C++編程完全也能夠藉助這些加速庫。針對Intel的CPU，你可以看看Intel的編譯器，裡面也有BLAS方面的介紹。

BTW：上面有人說2^2.367的乘法，其實不是矩陣運算吧，如果是N^2.367，其中N是個m*n的矩陣，那問題又不一樣了。

作為MATLAB的核心功能，其對矩陣運算的優化，要超越恐怕沒那麼容易。首先，它用的演算法肯定不是直接按定義傻乘，而是用了幾乎是世界上最好的矩陣乘法演算法。其次，它對具有特殊性質的大矩陣可能還會有特別的優化。最後，它針對具體硬體條件會有不少很重要的微調(調整cache footprint，利用GPU加速，利用SSE等等……）。JVM的性能penalty可能多少會有一些，但是一般沒有那麼大。現在的java還是做得蠻厲害的。如果完全按照MATLAB同樣的演算法設計重新用C或者彙編寫一遍，樂觀估計也可能只會比MATLAB快很少的那麼一點點。想做得更快的話更簡單的辦法是換硬體，比如用FPGA。

GPU加速，cuda，不過matlab也支持cuda了

使用好的 blas 庫即可.

如果不是稀疏矩陣，可以用BLAS，要是有GPU，可以用CUBLAS，比BLAS塊5倍以上。

矩陣類計算2010以後的matlab基本是調用mkl，mkl底層應該是彙編實現的，而且針對不同CPU構架專業優化，你想超過mkl的速度不是沒可能，但是需要花費的精力會很大，需要你對特定計算的演算法和硬體特點以及優化技術非常熟悉，而且就算比mkl快應該也快不了多少。所以Matlab在intel平台上涉及矩陣的運算，你基本可以理解為平台速度的上限或接近上限。

當然有一點需要注意，由於寄存器等壓力，實際上密集矩陣計算需要禁用超線程，禁用超線程會使mkl不少計算密集的矩陣運算獲得顯著提速。

matlab自己並不是先矩陣乘法, 它是依賴底層的MKL做到的.

業界MKL確實是最快的矩陣乘法. 直接給一個參考鏈接吧.

How does BLAS get such extreme performance?

看第三個回答, 及其對應的鏈接, 尤其是第一個參考鏈接, 非常詳細.

Matlab 好歹是標杆，別說怎麼更快，看看怎麼差不多吧。
https://modelingguru.nasa.gov/docs/DOC-1762
這裡是 Matlab 和其它一些庫的對比，其中 NumPy 差的不多，所以至少可以參考 NumPy 中的矩陣運算實現方式。
簡單說演算法是王道。能有 GPU 或者多核優化的演算法更好。
寫個三層 for loop 然後指望什麼彙編實現 cache optimisation 就能趕上 Matlab 的趁早洗洗睡。

@劉賀答案是怎麼得了三票還被踩到最下面去的？

你好，這是一個很有意思的話題。首先，演算法上來講，矩陣運算的時間複雜度目前最好是o(n^2.7)左右，超過這個的新演算法可以去領圖靈獎了。
所以目前加速矩陣運算都在硬體加速這一塊。對於特別大的矩陣，一般進行分散式計算，比如經典的scalapck，和最新的dplasma。CUDA的出現是對中型矩陣運算很好的覆蓋，因為cluster會花很多時間在communication，而且對程序員要求比較高如果程序優化的很好。
我們最近開發一款計算kernel給Matlab，R，GSL, Coin-OR這些高層演算法軟體，來讓他們自動把計算放到GPU而且是out of memory的，這樣媽媽再也不用手動移動數據從cpu到gpu。未來的release會包含多GPU的支持。用我們的BlasX編譯出來的Lapack，雖然不如手動優化的Magam，但是不會差到哪裡。

這個Proj叫BlasX，我們馬上會放到網上，敬請期待。

估計用GPU運算可能會更快...
另外, intel CPU的話換intel C compiler試一下, 據說效果很好

少不了彙編和硬體優化
DDT/crc at master · DingSoung/DDT · GitHub

矩陣運算，再也沒有比matlab更快的了，以前沒有，未來估計也不會有

可以換fortran語言寫代碼，fortran計算大數組比matlab快許多

INTEL乾爹發布的KNL處理器可以通過MKL庫來進行加速，親測8000*8000的矩陣相乘可以達到4TFlops。

矩陣運算，matlab無比匹敵，超過它的似乎還沒誕生。數學演算法優化、GPU加速、並行運算等它都早用上了。其矩陣運算比嵌套循環要快好幾個數量級。