0 引言

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

1 相關(guān)工作

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

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

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

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

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

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

    ALOGEELY M A等人提出的新型脈動陣列結(jié)構(gòu)^[10]改進(jìn)了KUMAR V K P等人提出的乘法器^[8-9]，通過對每個Block內(nèi)的運算單元數(shù)量做了裁剪，構(gòu)成梯度式增長的鏈?zhǔn)匠朔ㄆ?，這種結(jié)構(gòu)解決了運算啟動的延時問題，但數(shù)據(jù)組織較為復(fù)雜，不具備良好的拓展性。

    針對上述問題，本文提出了一種低IO帶寬矩陣陣列乘法器——鏈?zhǔn)骄仃嚦朔ㄆ?，主要特點是用列向量乘行向量和分時累加操作替代行向量乘列向量運算，使得乘法器內(nèi)所有運算單元都參與每條向量計算過程，充分發(fā)揮了數(shù)據(jù)和運算單元可復(fù)用特性。本文給出了所設(shè)計鏈?zhǔn)骄仃嚦朔ㄆ鞯墓ぷ髟砗蛯?yīng)硬件電路架構(gòu)，并在FPGA芯片上完成乘法器硬件實現(xiàn)和性能測試。多種規(guī)模的矩陣乘法實測結(jié)果證實，本設(shè)計在極低IO帶寬下達(dá)到了經(jīng)典陣列乘法器計算性能。

2 鏈?zhǔn)匠朔ㄆ鞴ぷ髟?/strong>

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

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

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

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

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

    式(2)中，矩陣A、B分別以列主式和行主方式展開乘法運算。由此可構(gòu)造一種無源數(shù)據(jù)緩存的3個數(shù)據(jù)通道(分別對應(yīng)源矩陣A和B，結(jié)果矩陣C)的鏈?zhǔn)匠朔ㄆ?，其計算過程的偽代碼如圖5所示，其運算路線如圖2中帶括號的坐標(biāo)。

    綜上，本文提出了一種IO帶寬恒為3fB的鏈?zhǔn)匠朔ㄆ鳎仃嘇、B中元素分別按列和行輸入乘法器內(nèi)運算單元完成確定操作。該鏈?zhǔn)匠朔ㄆ髦С衷诰€配置，能夠根據(jù)待處理矩陣規(guī)模被配置多鏈結(jié)構(gòu)，適應(yīng)不同的帶寬條件。

3 鏈?zhǔn)骄仃嚦似饔布崿F(xiàn)

3.1 鏈?zhǔn)骄仃嚦朔ㄆ骷軜?gòu)設(shè)計

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

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

3.2 鏈?zhǔn)匠朔ㄆ鱌E設(shè)計

    如圖8所示，為了將鏈?zhǔn)匠朔ㄆ鞯妮敵鰩捓寐蔬_(dá)到最高，本文設(shè)計了如圖9所示PE的結(jié)構(gòu)，包含一個浮點數(shù)乘法器和一個浮點數(shù)加法器，一對用于脈動輸入的寄存器，以及一組深度為m的乒乓存儲器，分別用于緩存加法器結(jié)果和輸出矩陣運算結(jié)果。

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

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

    這種結(jié)構(gòu)將累加操作解耦，只將乘法和加法操作結(jié)合，使得PE更加獨立，同時避免了硬件實現(xiàn)中累加器的流水線級數(shù)造成的結(jié)果數(shù)據(jù)延遲效應(yīng)。

3.3 鏈?zhǔn)匠朔ㄆ饔布O(shè)計

    以M階矩陣為例，鏈?zhǔn)匠朔ㄆ餍枰狹個PE。在進(jìn)行C=A×B的矩陣乘法時，輸入矩陣A和B分別按列和按行輸入到鏈?zhǔn)匠朔ㄆ鹘Y(jié)構(gòu)中。圖10所示為鏈?zhǔn)匠朔ㄆ鹘Y(jié)構(gòu)中M=3的具體結(jié)構(gòu)，矩陣B的數(shù)據(jù)元素每周期依次進(jìn)入PE₀的1端口，而矩陣A中的元素每周期依次對各個PE的0端口輪轉(zhuǎn)輸入。

    鏈?zhǔn)匠朔ㄆ鞯墓ぷ髁鞒倘鐖D11所示，具體對應(yīng)為一個3階方陣相乘的流程?？梢钥吹?，在初始M(M=3)個周期的啟動時間過后，所有PE處于滿載狀態(tài)。對于PE₀，直到第8個周期才輸出結(jié)果c11，c11是由第1、4、7周期的元素乘加得到的，三個階段的結(jié)果分別對應(yīng)圖中c11₁、c11₂、c11₃，最終求和得到的結(jié)果為矩陣A的第1行向量與矩陣B的第1列向量的乘累加結(jié)果。圖中灰色部分表示新的矩陣A′和B′輸入并相乘，此時的結(jié)果傳輸時間被隱藏。

    根據(jù)這種矩陣乘法的結(jié)構(gòu)特性，可以實現(xiàn)多個矩陣的連乘操作。由于沒有輸入緩沖，一組矩陣完成計算并切換到下一矩陣時沒有存儲上的時間消耗，結(jié)果存儲器采用乒乓操作，一部分用于參與當(dāng)前矩陣的實時運算，另外一部分存儲了上次矩陣乘結(jié)果并回傳，隱藏了大部分?jǐn)?shù)據(jù)回寫時間。

    在面對大規(guī)模矩陣乘的時候，由于硬件資源受限，難以一次性完成整個矩陣的運算，這時候必須將矩陣分成多個小塊，對各個子塊以及塊間進(jìn)行乘加運算，鏈?zhǔn)匠朔ㄆ髂軌蚝芎玫靥幚矸謮K矩陣乘。

4 鏈?zhǔn)匠朔ㄆ餍阅芊治?/strong>

    下面以在FPGA上實現(xiàn)的鏈?zhǔn)匠朔ㄆ鱽韺ι鲜鲈O(shè)計的性能進(jìn)行分析，在本文中，設(shè)計在XC7V2000T芯片上進(jìn)行驗證，實驗中采取單精度浮點數(shù)作為數(shù)據(jù)元素，經(jīng)過實驗，片內(nèi)最大可集成800個PE。本設(shè)計中所使用的FPGA開發(fā)環(huán)境和仿真環(huán)境為Xilinx Vivado 2016.3及Mentor Graphics Modelsim SE-64 10.2c。

4.1 峰值計算性能

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

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

    在xc7v2000tflg1925-1 FPGA上實現(xiàn)該設(shè)計，F(xiàn)PGA內(nèi)部資源使用情況如表1所示，根據(jù)布線后仿真的結(jié)果，該矩陣乘法器在未做優(yōu)化的情況下工作頻率可達(dá)到100 MHz。由此得出，該設(shè)計的峰值單精度浮點計算性能可以達(dá)到156.25 GFLOPS。

    當(dāng)PE資源不足以支持一次性運算整個矩陣，需要分塊運算，表2中陰影部分表示運算不需分塊，非陰影部分表示需要分塊運算。可以看出，隨矩陣規(guī)模增大，不同鏈數(shù)對運算性能的影響逐漸減小，圖中，在處理128×128及更大規(guī)模的矩陣乘時，所有鏈?zhǔn)匠朔ㄆ鬟\算周期基本持平，這是因為在處理這些矩陣分塊運算時，所有PE均處于滿負(fù)荷狀態(tài)，而矩陣分塊僅僅改變的是同一數(shù)據(jù)的讀取次數(shù)，對總的運算量并不影響。

4.2 不同鏈數(shù)對運算的影響

    分析表2中所列的數(shù)據(jù)，可以得到，當(dāng)運算規(guī)模小的時候，鏈?zhǔn)匠朔ㄆ鞯恼w運算時間相較于數(shù)據(jù)運算量來說較大，這是因為數(shù)據(jù)回傳時間對于總運算周期占比較大，隨著矩陣規(guī)模的增大，數(shù)據(jù)回傳時間占總周期比重越來越小，計算效率逐漸提高。

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

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

    設(shè)PPB表示鏈?zhǔn)浇Y(jié)構(gòu)單位帶寬的運算能力，并以脈動陣列作為衡量標(biāo)準(zhǔn)，PPB為脈動陣列所有帶寬的運算能力與鏈?zhǔn)浇Y(jié)構(gòu)的運算能力的歸一化比值，滿足如下關(guān)系：

    式(4)中，鏈數(shù)L₁表示脈動陣列規(guī)模(本實驗中L₁=8)，L₂表示鏈?zhǔn)匠朔ㄆ鞯逆湐?shù)，T_systolic和T_chain分別表示脈動和鏈?zhǔn)浇Y(jié)構(gòu)的運算周期，將表2中數(shù)據(jù)帶入式(2)，得到鏈?zhǔn)匠朔ㄆ鞯腜PB如圖13所示。

    當(dāng)運算小規(guī)模矩陣時，鏈?zhǔn)匠朔ㄆ鲉挝粠掃\算能力較低，然而當(dāng)計算大規(guī)模矩陣乘時，由于鏈?zhǔn)匠朔ㄆ髂軌蚝芎玫靥幚砭仃嚪謮K造成的運算效率問題，因此單位帶寬的運算能力逐漸提高，并逐漸趨近于L₁與L₂的比值。例如本實驗中，單鏈模式下，PPB最終趨近于8，表示鏈?zhǔn)匠朔ㄆ鲉挝粠挼倪\算能力與脈動陣列的8帶寬的運算能力相當(dāng)。由圖13可知，本文所設(shè)計的鏈?zhǔn)匠朔ㄆ髟趩捂?L=1)下工作性能最佳，同比與其他模式下的運算器，單鏈乘法器在帶寬方面有著明顯的優(yōu)勢。

4.3 不同數(shù)量運算單元對運算的影響

    不同規(guī)模的鏈?zhǔn)匠朔ㄆ鲗τ诰仃嚦朔ㄒ灿胁煌募铀傩Ч疚姆謩e以PE集成度為16、32、64、128的單鏈乘法器為實驗對象，計算不同規(guī)模的矩陣乘，見表3。

    表3中的數(shù)字表示總運算周期，陰影部分表示不分塊運算結(jié)果，非陰影部分表示分塊運算結(jié)果?？梢钥闯?，在處理小規(guī)模矩陣時，所有運算器性能相當(dāng)，這是因為PE資源足夠支撐一次性運算；在處理大規(guī)模矩陣時，PE數(shù)量成為限制運算器性能的主要因素，在處理128×128及更大規(guī)模的矩陣乘時，運算周期隨PE數(shù)量呈遞減趨勢。

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

    綜上所述，本文提出的鏈?zhǔn)匠朔ㄆ骺梢蕴幚聿煌?guī)模的矩陣乘法，并且運算性能卓越。

5 結(jié)論

    本文設(shè)計了一種鏈?zhǔn)讲⑿懈↑c數(shù)矩陣乘法器，并在Xilinx的XC7V2000T芯片進(jìn)行了原型驗證，通過實驗分析，該矩陣乘法器優(yōu)點在于對IO帶寬的要求非常低，結(jié)構(gòu)靈活，能夠適用于不同類型的矩陣乘法。該設(shè)計主要面向大規(guī)模矩陣的分塊運算，由于數(shù)據(jù)采取非緩沖的組織形式，運算器能夠?qū)崿F(xiàn)數(shù)據(jù)流入的同時開始計算，數(shù)據(jù)流入完畢結(jié)束運算，極大地提高了整體運算的吞吐率；由于PE相對獨立，支持根據(jù)不同的矩陣規(guī)模根據(jù)配置信息在線調(diào)整PE的使用情況，滿足了靈活性與通用性的要求，同時具備了良好的拓展性，鏈?zhǔn)匠朔ㄆ髟谶\算大規(guī)模矩陣乘法時表現(xiàn)突出。

參考文獻(xiàn)

[1] 馮子勇.基于深度學(xué)習(xí)的圖像特征學(xué)習(xí)和分類方法的研究及應(yīng)用[D].廣州：華南理工大學(xué)，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] 沈俊忠，肖濤，喬寓然，等.一種支持優(yōu)化分塊策略的矩陣乘加速器設(shè)計[J].計算機工程與科學(xué)，2016，38(9)：1748-1754.

[7] 田翔，周凡，陳耀武，等.基于FPGA的實時雙精度浮點矩陣乘法器設(shè)計[J].浙江大學(xué)學(xué)報(工學(xué)版)，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.

作者信息:

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

(合肥工業(yè)大學(xué) 電子科學(xué)與應(yīng)用物理學(xué)院，安徽 合肥230009)