提起矩陣計算，學過《高等數學》的人可能都聽過，但若不是這個領域的研究者，恐怕也只停在“聽過”的程度。在矩陣計算領域，開源項目OpenBLAS影響巨大，除IBM、華為等巨頭公司在使用外，還吸引了全球的研究院校、開發者們關注。
雷鋒網 AI 研習社近日有幸邀請到了澎峰科技創始人、OpenBLAS專案創始人和主要維護者張先軼，他將為我們介紹OpenBLAS開源項目以及矩陣乘法的優化。
嘉賓介紹
張先軼，中國科學院博士，MIT博士後，OpenBLAS開源專案創始人和主要維護者，PerfXLab澎峰科技創始人。曾獲2016年中國電腦學會科技進步二等獎。
張先軼師從張雲泉教授，在中科院取得高性能計算博士學位。讀博期間，基於GotoBLAS的原有基礎，他創建了開源矩陣計算庫OpenBLAS，領導團隊不斷進行修補和維護，目前在矩陣計算的細分領域，成為影響力較大的開源項目。
他在MIT萌生創業想法，歸國之後，針對“深度學習”，創辦PerfXLab澎峰科技，為電腦視覺、語音辨識等公司提供一體化性能優化方案。雷鋒網【新智造】頻道此前曾採訪並報導了PerfXLab澎峰科技。
課程內容
OpenBLAS專案介紹
矩陣乘法優化演算法
一步步調優實現
以下為公開課完整視頻，共64分鐘：
以下為公開課內容的文字及 PPT 整理。
雷鋒網的朋友們大家好，我是張先軼，今天主要介紹一下我們的開源矩陣計算庫OpenBLAS以及矩陣乘法的優化。
首先，什麼是BLAS？
BLAS是 Basic Linear Algebra Subprograms （基本線性代數副程式）的首字母縮寫，主要用來做基礎的矩陣計算，或者是向量計算。它分為三級：
BLAS 1級，主要做向量與向量間的dot或乘加運算，對應元素的計算；
BLAS 2級，主要做矩陣和向量，就類似PPT中藍色部分所示，矩陣A*向量x， 得到一個向量y。除此之外，可能還會有對稱的矩陣變形；
BLAS 3級，主要是矩陣和矩陣的計算，最典型的是A矩陣*B矩陣，得到一個C矩陣。由矩陣的寬、高，得到一個m*n的C矩陣。
為什麼BLAS是一個非常重要的庫或者介面，是因為它是很多科學計算的核心之一。每年做超級電腦的排行榜，都要做LINPACK測試，該測試很多部分就是做BLAS 3級矩陣和矩陣的計算。此外，還有很多科學和工程的類比，在轉換後都變成了一種矩陣上的操作。如果你把矩陣優化的特別好的話，對整個應用的提升，都是非常有幫助的。
BLAS與 Deep Learning 的關係，深度學習這幾年火了之後，比如大家非常瞭解的Alexnet，如果做具體的性能劃分，PPT上的這張圖來源於某篇論文，cut down之後看每個部分花了多少時間，發現它大部分的時間花費在卷積層（Conv Layer），另外不少時間花在了全連接層（FC layer）。卷基層目前通用的實現是展成矩陣，變成矩陣與矩陣的乘法，就是BLAS 3級。而全連接層一般是變成一個矩陣和向量的乘法，也落成了BLAS操作。
也就是說，基於矩陣類學習的深度學習，有90%或者更多的時間是通過BLAS來操作的。當然，隨著新的演算法出現，卷積層對3*3的卷積核有專門的演算法，或者用FFT類類演算法也可以做，但是在通用上，展矩陣來做也非常廣泛。
目前，BLAS只是一個定義了的實現介面，但是具體的實現其實有很多種。從商業的角度來講，存在很多商業版本，比如說 Intel、AMD、NVIDIA、IBM公司等，基本上為了搭配自己的硬體，對其做了更優的優化，出了商業版本。
針對開源的而言，有如下幾種，之前比較知名的是GoToBLAS，和OpenBLAS有很深的淵源，但是在2010年終止開發了，有時間在給大家分析其背後的原因，主力的開發人員後藤，離開了UT Austin的研究組，進入了公司，就終止了開發。ATLAS是美國一個學校做的，OpenBLAS是我們基於GotoBLAS做的，BLIS是後藤走了之後，基於GotoBLAS擴展出來的一個項目，目前還處在相對早期的階段，成熟度還差一些。
OpenBLAS歷史已經有幾年了，從2011年初開始進入，最初的原因是GotoBLAS放棄了，我們重新fork了一個社區發行版本，繼續開發和維護，它的維護不是一個簡單修BUG的過程，如果想要獲得比較好的性能，需要不停跟著硬體走，比如說新出一種新的硬體架構，或者適配更廣的硬體架構，都要進行一定的優化，來獲得比較好的加速效果。OpenBLAS算是目前全球最好的開源矩陣計算庫，在去年的時候得到了中國電腦學會科技進步二等獎，同時也進入了很多主流的Linux安裝包，比如說Ubuntu裡面就有我們的OpenBLAS Package，你能想到的Linux發行版本幾乎都進去了，但這不是我們主動去做的。還有一個OpenHPC的套件，也是最近做高性能計算的一個源。
目前OpenBLAS的進展是，支援幾乎全部的主流CPU處理器，同時都能達到比較好的優化性能。從作業系統來說，基本上常見主流的OS都支援。整體上，從適配的處理器範圍和支援的作業系統，在開源庫中算是最廣的實現。
因此，OpenBLAS的用戶也是比較多的。比如有開源專案Julia語言、GNU octave等；深度學習方面有大家熟悉的mxnet、Caffe都可以選OpenBLAS，作為CPU端的計算backend; IBM公司、ARM公司也都在他們的產品裡邊使用了OpenBLAS，NVIDIA公司在做一些跟CPU的對比測試時，把OpenBLAS列為了一個基準。其他還有一些做編譯器的以色列創業公司，還有國內的一些互聯網公司，比如搜狗。前段時間還和搜狗公司的人聊過，我們的庫線上上已經穩定運行一年多以上的時間，所以說它的工程品質上還是還是可以的。
IBM前段時間，因為深度學習非常火，做了一個Power AI的軟體框架，可以看到，最上面一層是一些開源的框架，底層的計算中就有我們的OpenBLAS。當然是為了搭載他的伺服器。
簡要的看一下性能，BLAS庫的性能是越高越好。在Intel的 Sandy Bridge平臺上，相比MKL的性能，基本上是重合在一起的，達到了一個相當的性能。
這張圖展示了在龍芯上做的一個結果，測得比較全，整體的BLAS多執行緒的，性能全測試了，性能比較高的都是我們，提高了一倍到兩倍。這是因為我們針對龍芯3A做了優化，所以取得了非常好的效果。
剛才主要介紹了OpenBLAS的性能和效果，我們在GitHub上做了託管，歡迎對矩陣乘法或優化感興趣的同學能加入進來，貢獻代碼，我們公司願意拿出一筆錢來支持這個項目繼續往前走。接下來會開始一些技術類的乾貨，主要講一下大家對優化比較感興趣的部分，我參考了矩陣乘法的這幾篇教程，UT Austin Flame組做的教程。我把他的內容基本上是摳出來了，一步步帶著大家過一下，如果我們從最簡單的矩陣乘法實現，到一個高性能的矩陣乘法實現，大概是幾步，怎麼來的？或者是為什麼優化，每一步能獲得多少性能收益。這樣大家對於一些優化的其他程式，希望能提供一些幫助。
我們首先看一下基本實現。
我想只要學過《線性代數》之類的，這種矩陣乘法，是一個非常簡單的問題，如果轉換成C代碼來做的話，就是一個三重迴圈，我在這張圖裡列出了一個【i j k】的三重迴圈，這裡面矩陣乘法的代碼就已經是，它實現的功能就是矩陣A*矩陣B，加到矩陣C裡面，C是結果矩陣，這裡面C的代碼，和在課本上看到的累加的公式是一樣的。找到i行，對應這個位置的結果C，把i行的每個元素，都取出來乘以B列，對應的乘，然後加起來就可以得到這個結果。但是這種實現，如果你放到現在的處理器上跑性能的話，和優化後的BLAS庫的實現，性能會差很多倍，甚至會差10倍。
為什麼呢，我們就做了一下最後的性能測試
這張圖也是截自教程裡，代表了一個性能圖，越高越好。這裡的測試平臺是Intel core i5 ，只是測了單執行緒，沒管多執行緒的事情。這種初始實現可能是1 GFlop/s。隨著規模變大，矩陣的性能在下降是為什麼呢？因為在實現的過程中，沒有考慮到cache的原因，當矩陣比較小的時候，速度還能快一些，當矩陣大了的時候，一定會跌下去，所以圖裡就有一個下滑的過程。
這個是非常原始、基礎的實現。
再往上走一步，怎麼樣才能再稍微優化一下。我們需要把矩陣的乘法順序調一下，我們在這裡做了一個小的分塊，把p單獨提到了一個函數裡，以點乘的形式寫出來，每次做一個1*4的結果，單獨提出來變成一個函數。p的這一步，要把計算順序稍微換一下，把i放到裡面，j放到外面，這塊背景為什麼要換一下，實際上是因為我們假設矩陣在存儲的時候是以列優先存儲的，在列項的數值是連續存儲，行之間是有間隔的，這對於仿存更有優勢。變成這樣的實現之後，對整體的性能其實沒什麼説明，那為什麼換成這種形式來寫呢，是為了之後的優化，留下一定的空間。
當矩陣比較小，在cache裡面的時候，性能基本是沒什麼變化的，但是超過cache的時候，它的性能稍微好了一點，從剛才非常低的值，也達到了接近1GFlop/s主要的原因是對A(0,p)做了一定的複用，它省了一些cache。另外一方面，它本身在做迴圈的利用來說，相當於比這部分做了一定迴圈的展開，所以在效率上得到了一定的提升。
這塊的複用，只從記憶體讀取一次，所以對超過cache的情況有了一定改善。
在這個基礎上，我們就需要看一下有什麼更好的方法來做優化。我們的基準就是，AddDot1*4的基準上怎麼做，我們想到第一點做的是，我們可不可以用寄存器變數來做，而不是操作記憶體。我可以申請一堆C 00，01這樣的寄存器變數，在C語言中是register double，還有矩陣A的部分，也用寄存器變數。
剩下的操作都是一些寄存器的變數，當然B還是之前的方式，最後再寫回C裡面。它完成的流程基本跟與之前的實習一樣，只是我們引入了寄存器變數，讓更多的資料保存到寄存器裡，而不是放到cache緩存裡，來減輕cache的壓力，這也能獲得一部分性能的提升。
可以看到，紅色的線是我們優化後的性能，達到了2GFlop/s，藍色的線是之前的性能，這裡做的優化就是利用寄存器降低cache的消耗，取得了50%左右的性能提升。完成了這一步之後，我們還可以再怎麼樣做優化呢？
我們剛才在A、C的部分，已經用寄存器做了替換，那麼B仿存這部分，我們有沒有可能也做一些優化。在之前實現的時候，B部分每次的座標計算都需要算出來，B訪問每個位置都要算一次，這樣就引入了額外的開銷，其實是可以避免的。
我們使用指標的方式，一開始先把對應的指標位置指好，每次計算的時候只要指標連續移動就好，而不是每次讀一個位置重新算一遍，這樣速度就會快一些。
我們看一下，做完這個小優化之後，降低了B index的消耗之後，從剛才的2G F…達到了4G的性能。這塊的改善對於小矩陣有效果，大矩陣全都超出了cache範圍，就不會有效果的。所以假設矩陣都在cache裡面，這塊也獲得了不小的性能提升。
當你完成這一部分的時候，你可以想，我把A矩陣做了寄存器替換，B矩陣通過index改進，我們下一步還能怎麼做？
其實就是一個比較常用的方式，做展開。
在最裡層迴圈，是不是可以展開成4次，在做這個的時候，我們可以降低整個迴圈這部分的開銷，而且讓它流水的情況更好。
這部分也會對性能有一些改善。這張圖比較的當初在中間階段的時候比起開始階段，得到了多少提升。通過使用寄存器變數，使用了指標，在做了一定的底層迴圈展開之後，達到了紅色線的性能，已經比藍色線有了明顯的提升，但是這個還不算完，只是一個基礎。但是在1*4的層面，已經沒什麼油水可挖了，所以我們需要在更上層找一些方法。
在1*4的時候，A的值產生了一些重用，但是如果塊再放大一點，比如說變成4*4時，也就是說每次計算的時候算的是一個方塊，而不是列。這個對於整個的優化來說，可以複用你的訪存，而且能夠更充分的發揮計算能力。
當我們變成小的這種4*4的方塊時，AddDot函數也要寫成這個模式。當然，這部分也要用剛才做過的那些1*4的方法，A這邊之前是1個值，現在是4個值，用寄存器的變數，C部分已經是4*4共有16個，也全都是寄存器變數，B的部分全部用指標來優化。
但這樣做的話，對於整體的性能提升是比較有限的。因為這只是一個初始的結果，可以看到，對於小矩陣，在cache範圍內，沒有什麼起色。但是超過cache後，對於大規模的矩陣，是有了一定性能提升。
在以4*4的結果優化基礎上，我們可以再做進一步的優化和開發。為什麼要轉換成4*4的優化，是因為我們現在CPU的處理器，基本上想獲得高的性能，必須要用向量化指令，不管是老的SSE2，AVX或者AVX 2.0等，對於CPU的優化，如果想達到高性能，必須要用到單指令多資料的向量化指令。
做了4*4的分塊之後，在這種情況下，會更有利於向量指令。在這裡以向量指令重寫了這部分迴圈的內容，當然這部分指令我沒有任何的內嵌彙編或者純彙編的操作，我就直接用了Intel Intrinsic的形式來寫，可以看到這部分寫的就是一個128位的sse，這是做一個雙精度浮點double的一個矩陣，資料都是double的，從A裡把這兩個值load進來。後兩個load進這個向量寄存器裡，B部分每次都要用load複製的這種指令load進去，剩下的這塊基本都是用向量的Intrinsic的變數來做了操作，在這塊跟之前看起來差不多，所以在編譯的時候都變成了向量化的指令。這兩行就算前部分C的值，這部分就算後部分C的值。
通過這種向量寄存器的指令使用後，他的性能提升非常明顯，從剛才4G可以達到超過6G ，而且這一部分是一個整體的變化，cache內向量加速效果是非常明顯的，基本上是翻了一倍。
下一步需要解決的是這個cache的問題，問題是沒有做大的分塊，超過cache大小之後性能就會下滑，要解決這個問題的話，需要在更上一層做Blocking。
轉換成代碼的話，在這一層做一個K的切分，下面一層做一個m的切分，至於kc和mc都是參數。這些參數都是可調的，都要根據L2 cache的大小進行調整，然後每次做的是一個小塊c的矩陣乘，相當於一個子問題，這個子問題的實現基本和剛才4*4的實現是一樣的。
這一部分blocking做完的性能如圖所示，藍色的線是沒有做Blocking的性能，紅色線是做過之後。當問題規模在cache內，它的性能改善基本比較小，但是當大規模的矩陣，通過做了這次Blocking之後，可以看到獲得了非常明顯的提升，變快了2倍。
對於大矩陣，為了充分的利用cache，讓子問題變小，提升它的資料局部性，在做其他問題優化的時候也很有必要的。下一步當我們做到blocking的時候，如果只是代碼級別變化的時候，基本已經做完了。
此後再進一步優化的點，引入引入一些操作。當我們分析程式存在的 性能瓶頸，對於A的訪存和B的訪存是比較慢，很多訪存在矩陣中是不連續的，所以訪存性能就差了很多，一方面不能利用cache，一方面在TLB上也有影響，當然C部分也有一些影響，C矩陣往往很大，沒有辦法做packing，只能對A和B來做，packing的意思是說，我在這裡有一部分連續的記憶體空間，m*k，對應前面的mc和kc，在這塊記憶體空間，在每次做計算之前，我把所需要用到的A的矩陣，從原始矩陣讀取出來，存放到連續的一塊記憶體空間裡。 Packed Matrix A這個函數的具體實現非常簡單，基本上就是從對應的位置取出來，放在連續的記憶體位址就結束。
為什麼會做這步操作呢？這步操作的意義在於，通過pack之後，下一步AddDot4*4裡讀的元素就不是從A矩陣讀，而是從pack後緩存區的位置讀。一個好處是，A矩陣已經預熱，放進CPU的cache裡了；第二個好處是，你可以看到我在存儲的時候，這種連續性的存儲，讀的時候也是連續性讀取，效率會非常高，cache效率也非常高。加上通過pack這一步，對於性能的改善，是非常明顯的。
這張圖是上一步的操作，做了packing之後，除了極小矩陣規模沒什麼效果，或者引入了額外開銷，效果還變差之外，絕大部分的性能提升是非常明顯的，有50%以上，達到了9GFlop/s。
對於矩陣乘法實現的話，packing矩陣是一個非常重要的優化方式。再往後大家會想，對於A來說做了Packing，對於B是不是也能做Packing，同樣道理也是可以的，就是把它拷貝到一個連續空間。B部分的Packing操作和A部分類似，也是把它的資料從原始矩陣裡讀出來，然後放到一個連續空間裡，使它的記憶體訪問做連續的訪存。當然這部分，因為B訪存是個流式的訪存，所以在這部分的改進會稍微小一點，相比A只有大概10%左右的提升。
對於矩陣乘法實現的話，packing矩陣是一個非常重要的優化方式。再往後大家會想，對於A來說做了Packing，對於B是不是也能做Packing，同樣道理也是可以的，就是把它拷貝到一個連續空間。B部分的Packing操作和A部分類似，也是把它的資料從原始矩陣裡讀出來，然後放到一個連續空間裡，使它的記憶體訪問做連續的訪存。當然這部分，因為B訪存是個流式的訪存，所以在這部分的改進會稍微小一點，相比A只有大概10%左右的提升。
當你完成到這一步的時候，相比最開始三重迴圈的性能改進，你的矩陣乘法的性能已經有很明顯的提升了。你如果想做的更好的話，內部的核心可能不止要寫intrinsic的指令，還要寫內嵌彙編，重排流水線，使硬體資源能夠發揮更多，可能還會提升10%。當然這部分對實現BLAS比較重要，會摳的比較細。
我們再整體回顧一下矩陣乘法的演算法，我把演算法的這部分放到最後，從開始一步步實現之後，做到最後大家可能看的比較清楚一些。A矩陣*B矩陣得到C矩陣，對應的是最外層的迴圈，每一步往下的時候，其實都是在做分塊，做分塊的原因剛才有提到，就是為了配合硬體結構，因為memory、cache等都是分層的，它是越來越小的，做分塊實際上是提高了cache的各層的利用率。
今天就分享到這裡，謝謝大家。
問答解答
問題1：什麼是訪存優化？
張先軼：訪存優化解決的是處理器讀取資料的性能。從計算上來說，是相對好優化的，但是優化訪存會非常困難，稠密矩陣乘法的資料還是相對規整的，讀數據的順序是有規則的，更容易優化一些。但是我們也做過很多疏鬆陣列的優化，比如疏鬆陣列乘向量的優化，這個對訪存來說更困難一些，因為沒有辦法預測到下一次訪存在什麼位置，這造成了優化的困難。
問題2：OpenBLAS和其他矩陣庫有什麼關係？
張先軼：OpenBLAS和其他BLAS實現其實都是完成了介面，BLAS只是介面的定義，具體來說可以有多種實現。我們認為我們OpenBLAS在開源的實現是非常好的。如果是標準BLAS，有參考實現，只是一個非常簡單的Fortran實現，性能很差的，我們要比他們快很多。MKL是Intel公司自己做的BLAS，我們跟他們相當。Engien我們沒有完全測過，它號稱自己做的很好，但是他們的做法在X86的平臺可能有些效果，但是對其他平臺的效果我表示懷疑。不過我沒有具體做對比測試，所以發言權不大。
問題3：從入門到精通需要多久？
張先軼：如果我指導的話，幾個月時間就可以上手做一些事情。歡迎大家。
問題4：比起高通的庫表現如何？
張先軼：說實話高通的庫沒有測過，我覺得它號稱比較快，是因為在32位的ARM上，我們OpenBLAS沒有做向量化優化，高通的那個部分做了，所以它可能會比我們快一些，但是在我們公司內部的PerfBLAS是優化了的。
問題5：分塊的目的是什麼？
張先軼：就是優化訪存，通過分塊之後讓訪存都集中在一定區域，能夠提高了資料局部性，從而提高cache利用率，性能就會更好。
問題6：硬體不給力能玩神經網路麼？
張先軼：我們給出的一個資料是，我們在ARM CortexA57的平臺上，4核1.7GHz，跑AlexNet模型，一張圖是150ms，這個速度還算比較快。另一方面我們還在做一些其他的模型，做了精簡優化，再配合底層庫的優化。在某個小模型下，在跑一張小圖片的inference只用了50ms。所以，在ARM的處理器上，還是可以做到即時當地語系化的神經網路inference。
問題7：內部版本和開源版本差別大麼？
張先軼：內部版本是針對深度學習做了一些差異化處理，性能高的可能會到1倍多，這部分的優化，一部分是矩陣的規模，和剛才講的做方陣不一樣，深度學習的矩陣大部分是中小型，某個維度會比較小，要用到的優化策略，或者分塊策略會不一樣，甚至有時候特別小根本不用分塊或packing，直接做可能更好一些。