0 引言

    在圖像視頻處理和機器學習領域，矩陣運算規模達數千維甚至上兆維^[1]。矩陣運算，尤其是矩陣乘法O(n³)成為影響上述應用實時性的關鍵。研究低成本和低時間開銷的大規模矩陣乘法求解方法具有極強的工程實用價值。

1 相關工作

    C=A×B(A、B和C矩陣均為M階)的矩陣乘法的結果數據之間無關，故C陣的元素可同時計算。利用這一特性，多種改進矩陣乘法器方法被提出獲得低時間復雜度，較為典型的有：STRASSEN V等人提出了利用張量代數實現的公式化矩陣乘法，降低矩陣乘的時間復雜度的Strassen算法^[2]；CANNON L E等人介紹了一種網格并行的Cannon算法^[3]，通過矩陣分塊和地址循環的方法實現矩陣乘，最快執行時間為O(n)；此后，FOX G C等人提出一種基于超立方體結構的Fox算法^[4]；KUNG H T等人提出了一種陣列結構的脈動算法^[5]；沈俊忠、田翔等人提出了兩種基于脈動的二維陣列實現方案^[6-7]。

    理論上采用上述矩陣陣列乘法器陣列規模越大，加速效果越顯著，但是實際硬件實現中陣列規模受兩個因素制約。

    首先，理論上脈動陣列的IO帶寬應滿足3fMB(這里f為工作頻率，M為陣列規模，B為數據位寬，當f和B一定時，可記為O(3M)），受存儲墻影響，M無法得到K級；

    其次，過高的M限制了陣列IO帶寬利用率。設脈動陣列輸入帶寬利用率B_ave定義為：

    式(1)的函數曲線如圖1所示，圖1表明隨陣列規模M的增長，B_ave值逐漸收斂為1/2。這就意味著單純增加M值會造成近一半IO帶寬的浪費，這又進一步加劇了IO帶寬的壓力。

    為克服經典脈動陣列結構的缺陷，KUMAR V K P等人提出一種IO帶寬固定的鏈式乘法器^[8-9]，該乘法器中每個運算單元內存儲容量固定，通過將矩陣A元素和矩陣B元素送到“低速通道”和“快速通道”來實現Block內和Block間的數據傳遞，這種數據組織形式將IO帶寬限制為9fB，避免了隨著矩陣規模增大帶來的帶寬問題，該結構另外的優點在于運算單元內部的存儲容量固定，但容量和運算單元數量呈反比；對于大規模矩陣乘，容量越小意味著運算單元數量越多，造成運算單元數量的浪費，并且由于線性乘法器僅有單一的數據輸入端，這會造成從數據輸入到所有運算單元都處于工作狀態的延時更大，降低了大規模矩陣運算效率。

    ALOGEELY M A等人提出的新型脈動陣列結構^[10]改進了KUMAR V K P等人提出的乘法器^[8-9]，通過對每個Block內的運算單元數量做了裁剪，構成梯度式增長的鏈式乘法器，這種結構解決了運算啟動的延時問題，但數據組織較為復雜，不具備良好的拓展性。

    針對上述問題，本文提出了一種低IO帶寬矩陣陣列乘法器——鏈式矩陣乘法器，主要特點是用列向量乘行向量和分時累加操作替代行向量乘列向量運算，使得乘法器內所有運算單元都參與每條向量計算過程，充分發揮了數據和運算單元可復用特性。本文給出了所設計鏈式矩陣乘法器的工作原理和對應硬件電路架構，并在FPGA芯片上完成乘法器硬件實現和性能測試。多種規模的矩陣乘法實測結果證實，本設計在極低IO帶寬下達到了經典陣列乘法器計算性能。

2 鏈式乘法器工作原理

    矩陣乘法運算C=A×B，(A、B和C均為M維矩陣)的偽碼如圖2所示。

    令C[0][i][j]=0；1≤i≤M；1≤j≤M，設M=3，可得如圖3左側所示的三階矩陣乘的DAG，圖中每個頂點代表一個乘累加運算，無括號坐標表征圖2中數據移動路線，該DAG的硬件實現如圖3右側所示。

    當IO帶寬滿足，規模為n×n的脈動乘法器可得計算時間下界O(n)。但受IO帶寬制約，典型脈動陣列規模有限，需要將大維度矩陣分解若干子陣分別處理，執行復雜的數據緩沖或訪存操作才能得到理想的加速效果。

    將圖3左側圖沿d₁=(0，0，1)方向投影可得圖4。

    圖2中沿著i坐標軸正方向存在三個與坐標軸jk平面平行的平面，第二個平面的輸入是行向量a₂和矩陣B，第三個平面輸入是行向量a₃和矩陣B。故其運算過程可描述為如下公式：

    式(2)中，矩陣A、B分別以列主式和行主方式展開乘法運算。由此可構造一種無源數據緩存的3個數據通道(分別對應源矩陣A和B，結果矩陣C)的鏈式乘法器，其計算過程的偽代碼如圖5所示，其運算路線如圖2中帶括號的坐標。

    綜上，本文提出了一種IO帶寬恒為3fB的鏈式乘法器，矩陣A、B中元素分別按列和行輸入乘法器內運算單元完成確定操作。該鏈式乘法器支持在線配置，能夠根據待處理矩陣規模被配置多鏈結構，適應不同的帶寬條件。

3 鏈式矩陣乘器硬件實現

3.1 鏈式矩陣乘法器架構設計

    由圖4得到行向量a₁與矩陣B相乘的數據組織形式如圖6(a)所示。圖中大方框中是用來執行乘累加操作并緩存結果矩陣元素的2組中間值，進而輸出最終結果。由此類推，行向量a₂和a₃與矩陣B相乘的數據組織形式如圖6(b)和圖6(c)所示。

    由圖6可知，M階(此時M=3)矩陣A在與M階矩陣B進行矩陣乘時，矩陣A中的行向量的每個元素都要保持M個周期才能向對應乘累加器中輸入下一個元素，即輸入帶寬利用率僅為矩陣A輸入總帶寬的1/M；矩陣B按行向每個乘累加器輸入相同的元素，其帶寬利用率也僅為矩陣B輸入總帶寬的1/M。觀察矩陣A的行向量的輸入規律，在行向量的每個元素的保持周期內，將矩陣A中的元素按列進行輸入(例如，在矩陣A第一個行向量的第一個元素a₁₁的保持周期內，將矩陣A中的元素a₁₁、a₂₁、a₃₁依次輸入到每一個乘累加器中)，可以得到三階矩陣列乘行的數據組織形式如圖7所示。圖7中矩陣A按列依次進行輸入，矩陣B按行依次進行輸入，但要使得矩陣B輸入端的帶寬得到充分的利用，矩陣B需要在每個相鄰乘累加器相差1cycle，故而可以將矩陣B也按照行輸入到每個乘累加器，通過使用寄存器將矩陣B的輸入延遲一個周期輸入到下一個乘累加器中，可以得到圖8的數據組織形式。

3.2 鏈式乘法器PE設計

    如圖8所示，為了將鏈式乘法器的輸出帶寬利用率達到最高，本文設計了如圖9所示PE的結構，包含一個浮點數乘法器和一個浮點數加法器，一對用于脈動輸入的寄存器，以及一組深度為m的乒乓存儲器，分別用于緩存加法器結果和輸出矩陣運算結果。

    當源數據srcA和srcB進入PE時，首先被送入浮點數乘法器中進行乘法運算，乘法結果會被送入浮點數加法器，與運算存儲器(用于緩存加法器結果的存儲器)中頂層(dstA0)數據相加，相加的結果進入運算存儲器底層(dstA(m-1))；加法器每進行一次運算，所有數據被刷新一次，即向前移動一個地址，直到處于頂層(dstA0)；反復執行上述操作，直到運算完畢。

    當前矩陣運算結束后，需要將結果輸出，此時存儲器執行乒乓操作，將運算存儲器和原先用于輸出矩陣運算結果的傳輸存儲器互相切換。

    這種結構將累加操作解耦，只將乘法和加法操作結合，使得PE更加獨立，同時避免了硬件實現中累加器的流水線級數造成的結果數據延遲效應。

3.3 鏈式乘法器硬件設計

    以M階矩陣為例，鏈式乘法器需要M個PE。在進行C=A×B的矩陣乘法時，輸入矩陣A和B分別按列和按行輸入到鏈式乘法器結構中。圖10所示為鏈式乘法器結構中M=3的具體結構，矩陣B的數據元素每周期依次進入PE₀的1端口，而矩陣A中的元素每周期依次對各個PE的0端口輪轉輸入。

    鏈式乘法器的工作流程如圖11所示，具體對應為一個3階方陣相乘的流程。可以看到，在初始M(M=3)個周期的啟動時間過后，所有PE處于滿載狀態。對于PE₀，直到第8個周期才輸出結果c11，c11是由第1、4、7周期的元素乘加得到的，三個階段的結果分別對應圖中c11₁、c11₂、c11₃，最終求和得到的結果為矩陣A的第1行向量與矩陣B的第1列向量的乘累加結果。圖中灰色部分表示新的矩陣A′和B′輸入并相乘，此時的結果傳輸時間被隱藏。

    根據這種矩陣乘法的結構特性，可以實現多個矩陣的連乘操作。由于沒有輸入緩沖，一組矩陣完成計算并切換到下一矩陣時沒有存儲上的時間消耗，結果存儲器采用乒乓操作，一部分用于參與當前矩陣的實時運算，另外一部分存儲了上次矩陣乘結果并回傳，隱藏了大部分數據回寫時間。

    在面對大規模矩陣乘的時候，由于硬件資源受限，難以一次性完成整個矩陣的運算，這時候必須將矩陣分成多個小塊，對各個子塊以及塊間進行乘加運算，鏈式乘法器能夠很好地處理分塊矩陣乘。

4 鏈式乘法器性能分析

    下面以在FPGA上實現的鏈式乘法器來對上述設計的性能進行分析，在本文中，設計在XC7V2000T芯片上進行驗證，實驗中采取單精度浮點數作為數據元素，經過實驗，片內最大可集成800個PE。本設計中所使用的FPGA開發環境和仿真環境為Xilinx Vivado 2016.3及Mentor Graphics Modelsim SE-64 10.2c。

4.1 峰值計算性能

    理想情況下，鏈式乘法器工作時所有PE均滿載，每個PE在單時鐘周期內消耗兩個浮點數，得到一個結果數據。實際上，由于矩陣類型、流水線延遲、存儲策略加上啟動延遲(數據從輸入開始到所有PE工作)等因素的影響，在某些周期內仍會存在部分PE處于空閑的情況，設PERF_peak表示運算的峰值性能，N表示PE的個數，f表示工作頻率，S表示整個大規模矩陣總運算次數，則有式(3)表示如下：

    本文主要圍繞大規模矩陣儲乘法進行硬件設計，通過探究N、S、f等因素對運算器的性能影響，進行了以下幾組實驗，實驗以200 MHz時鐘作為工作時鐘，以2片DDR3作為主存，存儲源數據和結果數據，2片DDR3峰值帶寬為204.8 Gb/s，考慮到多路并行輸入輸出時的帶寬限制問題，在實驗中矩陣元素采用32位浮點數表示，所以單條鏈式乘法器的峰值帶寬為18.75 Gb/s(2組源數據，一組結果數據)。

    在xc7v2000tflg1925-1 FPGA上實現該設計，FPGA內部資源使用情況如表1所示，根據布線后仿真的結果，該矩陣乘法器在未做優化的情況下工作頻率可達到100 MHz。由此得出，該設計的峰值單精度浮點計算性能可以達到156.25 GFLOPS。

    當PE資源不足以支持一次性運算整個矩陣，需要分塊運算，表2中陰影部分表示運算不需分塊，非陰影部分表示需要分塊運算。可以看出，隨矩陣規模增大，不同鏈數對運算性能的影響逐漸減小，圖中，在處理128×128及更大規模的矩陣乘時，所有鏈式乘法器運算周期基本持平，這是因為在處理這些矩陣分塊運算時，所有PE均處于滿負荷狀態，而矩陣分塊僅僅改變的是同一數據的讀取次數，對總的運算量并不影響。

4.2 不同鏈數對運算的影響

    分析表2中所列的數據，可以得到，當運算規模小的時候，鏈式乘法器的整體運算時間相較于數據運算量來說較大，這是因為數據回傳時間對于總運算周期占比較大，隨著矩陣規模的增大，數據回傳時間占總周期比重越來越小，計算效率逐漸提高。

    由于表2中數據范圍相差較大，為了直觀地得到各組乘法器的性能對比，對每組鏈的運算周期與理想脈動結構的周期做比值,得到結果如圖12所示。

    表2中同樣列舉了規模為8×8的脈動陣列的運算周期，200 MHz工作頻率下8×8脈動陣列的峰值帶寬為150 Gb/s，而本文設計的單鏈乘法器峰值帶寬僅18.75 Gb/s，由圖12可知，單鏈乘法器在處理小規模矩陣時與脈動相比性能有所欠缺，但在處理大規模矩陣乘時兩者表現相當，而單鏈乘法器在帶寬方面優勢非常明顯。

    設PPB表示鏈式結構單位帶寬的運算能力，并以脈動陣列作為衡量標準，PPB為脈動陣列所有帶寬的運算能力與鏈式結構的運算能力的歸一化比值，滿足如下關系：

    式(4)中，鏈數L₁表示脈動陣列規模(本實驗中L₁=8)，L₂表示鏈式乘法器的鏈數，T_systolic和T_chain分別表示脈動和鏈式結構的運算周期，將表2中數據帶入式(2)，得到鏈式乘法器的PPB如圖13所示。

    當運算小規模矩陣時，鏈式乘法器單位帶寬運算能力較低，然而當計算大規模矩陣乘時，由于鏈式乘法器能夠很好地處理矩陣分塊造成的運算效率問題，因此單位帶寬的運算能力逐漸提高，并逐漸趨近于L₁與L₂的比值。例如本實驗中，單鏈模式下，PPB最終趨近于8，表示鏈式乘法器單位帶寬的運算能力與脈動陣列的8帶寬的運算能力相當。由圖13可知，本文所設計的鏈式乘法器在單鏈(L=1)下工作性能最佳，同比與其他模式下的運算器，單鏈乘法器在帶寬方面有著明顯的優勢。

4.3 不同數量運算單元對運算的影響

    不同規模的鏈式乘法器對于矩陣乘法也有不同的加速效果，本文分別以PE集成度為16、32、64、128的單鏈乘法器為實驗對象，計算不同規模的矩陣乘，見表3。

    表3中的數字表示總運算周期，陰影部分表示不分塊運算結果，非陰影部分表示分塊運算結果。可以看出，在處理小規模矩陣時，所有運算器性能相當，這是因為PE資源足夠支撐一次性運算；在處理大規模矩陣時，PE數量成為限制運算器性能的主要因素，在處理128×128及更大規模的矩陣乘時，運算周期隨PE數量呈遞減趨勢。

    分塊運算時，受PE數量N限制，每次運算量為N×N，若將矩陣分割為K×K個N×N的矩陣塊，則總的運算周期近似N²×K³。設其他結構的乘法器所分塊得到的子陣規模為N′×N′，其分塊得到的子陣數量為K′×K′，可以得到，相同PE數量下，單鏈結構所對應的N值大于N′，K小于K′，所以總運算周期最優。

    綜上所述，本文提出的鏈式乘法器可以處理不同規模的矩陣乘法，并且運算性能卓越。

5 結論

    本文設計了一種鏈式并行浮點數矩陣乘法器，并在Xilinx的XC7V2000T芯片進行了原型驗證，通過實驗分析，該矩陣乘法器優點在于對IO帶寬的要求非常低，結構靈活，能夠適用于不同類型的矩陣乘法。該設計主要面向大規模矩陣的分塊運算，由于數據采取非緩沖的組織形式，運算器能夠實現數據流入的同時開始計算，數據流入完畢結束運算，極大地提高了整體運算的吞吐率；由于PE相對獨立，支持根據不同的矩陣規模根據配置信息在線調整PE的使用情況，滿足了靈活性與通用性的要求，同時具備了良好的拓展性，鏈式乘法器在運算大規模矩陣乘法時表現突出。

參考文獻

[1] 馮子勇.基于深度學習的圖像特征學習和分類方法的研究及應用[D].廣州：華南理工大學，2016.

[2] STRASSEN V.Gaussian elimination is not optimal[J].Numerische Mathematik，1969，13(4)：354-356.

[3] CANNON L E.A cellular computer implement the Kalman filter algorithm[D].Bozeman US：Montana State University，1969.

[4] FOX G C，OTTO S W，HEY A J G.Matrix algorithm on a hypercube I: matrix multiplication[J].Parallel Computing，1987，4(1)：17-31.

[5] KUNG H T，LEISERSON C E.Systolic arrays(for VLSI)[J].Proceeding Sparse Matrix Conference，1978：256-282.

[6] 沈俊忠，肖濤，喬寓然，等.一種支持優化分塊策略的矩陣乘加速器設計[J].計算機工程與科學，2016，38(9)：1748-1754.

[7] 田翔，周凡，陳耀武，等.基于FPGA的實時雙精度浮點矩陣乘法器設計[J].浙江大學學報(工學版)，2008，42(9)：1611-1615.

[8] KUMAR V K P，TSAI Y C.On synthesizing optimal family of linear systolic arrays for matrix multiplication[J].IEEE Transactions on Computers，1991，40(6)：770-774.

[9] KUMAR V K P，TSAI Y C.On Mapping algorithms to linear and fault-tolerant systolic arrays[J].IEEE Transactions on Computers，1989，38(3)：470-478.

[10] ALOGEELY M A，Al-Turaigi M A，ALSHEBEILI S A.A new approach for the design of linear systolic arrays for computing third-order cumulants[J].Integration the Vlsi Journal，1997，24(1)：1-17.

作者信息:

宋宇鯤，鄭強強，王澤中，張多利

(合肥工業大學 電子科學與應用物理學院，安徽 合肥230009)