基于改進QR算法的矩陣分解器設(shè)計

陳文杰，宋宇鯤，張多利

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

矩陣分解是矩陣運算中重要的運算之一，在神經(jīng)網(wǎng)絡(luò)加速器[1-2]、多進多出(Multiple-Input Multiple-Output，MIMO)[3-5]、自適應(yīng)濾波器[6]等領(lǐng)域有廣泛的應(yīng)用。常見的矩陣分解法主要有LU分解、QR分解等。LU分解包括Cholesky分解法、Dolittle分解法等。QR分解包括Gram-Schmidt正交分解法(MGS)、Givens-Rotation分解法(GR)、Householder分解法(HT)等。不同的分解法適用于不同的應(yīng)用場景，在高精度、高穩(wěn)定性的設(shè)計前提下，QR分解憑借其高數(shù)值穩(wěn)定性和對輸入矩陣無要求等特點得到了廣泛的關(guān)注。

在QR分解算法中，Gram-Schmidt正交分解算法的運算并行度高，但數(shù)值穩(wěn)定性較差[7]；Householder分解法運算復(fù)雜度過高[8]，一般不適用于工程應(yīng)用中；Givens-Rotation分解法在高精度計算時數(shù)值穩(wěn)定性較好，且算法本身存在并行性，易于硬件實現(xiàn)[9]。GR算法最早在文獻[10]中被提出，之后大量研究學(xué)者對其進行了改進。文獻[11～13]采用牛頓迭代法將GR算法中的開方操作省略，降低了算法的資源消耗。文獻[14]采用二維陣列結(jié)構(gòu)實現(xiàn)QR算法，提高運算并行度，縮短了分解時間。文獻[15]采用細粒度并行的Givens分解算法，縮短了分解的運算周期，使其具有良好的可擴展性，結(jié)構(gòu)也更加簡單。現(xiàn)有GR算法的硬件實現(xiàn)研究大多基于低維度、低精度的硬件實現(xiàn)，有關(guān)高維度、高精度的矩陣分解電路的研究較少。

在信息化高度發(fā)展的時代，信號處理系統(tǒng)需要實時輸出高維度、高精度的矩陣運算結(jié)果[16-17]。因此，本文在深入了解GR-QR算法的基礎(chǔ)上，通過查閱矩陣分解算法資料，找到一種基于按列Givens旋轉(zhuǎn)(Column-wise Givens Rotation，CGR-QR)的高并行度算法[18]，并基于該算法設(shè)計了一種高性能的CGR-QR分解電路結(jié)構(gòu)。該結(jié)構(gòu)消去了QR分解中多次迭代的運算過程，降低了消元過程的運算時間，提高了資源利用率。

1 算法原理

設(shè)A為m×n階非奇異矩陣，Q為m×m階正交矩陣，R為上三角矩陣，若滿足下式

A=QR

(1)

則稱其為矩陣A的QR分解[19]。

QR分解中，通過左乘Givens矩陣Gij，更新原矩陣的第i行第j列元素，其他行不變。通過不斷左乘Givens矩陣，可將矩陣的一列元素消元。

(2)

其中G矩陣的表達式如下

(3)

c與s為更新因子，表達式如下

(4)

將G的表達式帶入式(2)中，得到R(1)的計算式

(5)

式中，bij為需要消元的矩陣的第i行第j列的元素。式(5)中參數(shù)K、S、l、P的計算式如式(6)所示。

(6)

根據(jù)數(shù)學(xué)定義，可知矩陣Q為正交矩陣。經(jīng)過n-1次Givens旋轉(zhuǎn)，原矩陣被轉(zhuǎn)化為上三角矩陣R，如式(7)所示。

R=Qn…Q1A=

(7)

Q=QnQn-1…Q1Q=QnQn-1…Q1也為正交矩陣，R為上三角矩陣，上述過程即為矩陣A的CGR-QR分解。

2 QR分解的相關(guān)研究

2.1 QR-RD的細粒度并行結(jié)構(gòu)

文獻[15]提出一種并行的QR算法。該算法由一個外循環(huán)和兩個同級內(nèi)循環(huán)構(gòu)成，且兩個內(nèi)循環(huán)同時進行工作，其更新過程如圖1所示。

圖1 算法更新過程[15]

算法中，B(i)為外部循環(huán)中的參數(shù)計算過程，該過程負責(zé)QR分解過程中更新因子的計算。B1(i，j)為行更新過程，B2(i，j)為更新過程中的參數(shù)計算過程。

QR算法分為4個步驟：初始化階段、計算階段、同步階段、存儲階段和待定階段。每個處理元件(Processing Element，PE)在初始階段獲得其PE標(biāo)識符(pid)。計算階段是算法的主要階段，主要包括接收數(shù)據(jù)、更新數(shù)據(jù)和發(fā)送數(shù)據(jù)等階段，這些階段并行執(zhí)行。首個Loop從外部存儲器加載數(shù)據(jù)，之后的Loop從上一級Loop中獲得數(shù)據(jù)。最后一個Loop的結(jié)果將更新到存儲器中，其余的Loop結(jié)果傳到下一級Loop中。

算法細粒度并行QR分解算法。

輸入：源矩陣Ai。

輸出：分解得到的上三角矩陣Ri，共軛轉(zhuǎn)置矩陣Qi。

1.p:the number of PE pid: identifier of the current PE

2.Initial Phase:

3.SetPID[pid]；//設(shè)置當(dāng)前標(biāo)志符

4.Calculating://開始計算

5.if(pid==0)Load(Ai)；

6.else Rcv(pid-1，Ai)；//判斷并加載輸入矩陣

7.forj=(i+1)，…，n

8.Parallel do 1，2 and 3；//并行計算步驟1、2、3

9.1:if(pid==0)Load(Aj)；

10.esle Rcv(pid-1，Aj)；//加載當(dāng)前計算矩陣

11.2:Updata(Ai，Aj)；//更新矩陣數(shù)據(jù)

12.3:if(pid==p-1) Store(Aj)；

13.else Send(pid+1，Aj)；//判斷計算結(jié)果輸入到下一級還是存儲到存儲器

14.end

15.wait()；

16.Store1(pid，Ri，Qi)；//存儲分解得到的矩陣Ri和Qi

17.Ready Next Section：

18.Judge(pid)；//判斷當(dāng)前標(biāo)志符是否已經(jīng)完成分解。

主要消息傳遞基元的定義為：

(1)Load(X)指從外部存儲器中加載數(shù)據(jù)；

(2)Store(X)指將矩陣X存入外部存儲器中；

(3)Store 1(pid，XR，XQ)代表將PE編號為pid的矩陣XR和XQ存儲到外部存儲器中；

(4)Send(pid，Y)代表將矩陣Y發(fā)送到編號為pid的PE中；

(5)Rcv(pid，Y)代表從編號為pid的PE中接收Y矩陣數(shù)據(jù)。

細粒度并行分解硬件結(jié)構(gòu)如圖2所示。

圖2 細粒度并行分解硬件結(jié)構(gòu)

細粒度并行的QR分解算法將QR分解為更新因子計算、參數(shù)就算、行列更新3個部分，提高了分解過程中的運算并行度。由于該算法延續(xù)傳統(tǒng)的Gives-QR分解，每一次分解只能更新兩行行向量，雖然通過多組PE進行流水操作，但迭代次數(shù)并無減少。另外，該結(jié)構(gòu)采用大量的PE堆疊，運算過程中的中間數(shù)據(jù)需要額外的存儲資源寄存，在并行度較高時所需的運算資源較多。對于高速、高精度矩陣求逆器設(shè)計而言，需要在算法層面進行更新，消去矩陣元素的數(shù)據(jù)相關(guān)性，從而獲得高并行度的實現(xiàn)結(jié)構(gòu)，最大程度挖掘算法的運算性能。綜上所述，盡管細粒度并行的分解結(jié)構(gòu)性能優(yōu)異，但是不適用于高速、高精度大維度矩陣求逆器設(shè)計。

2.2 二維陣列結(jié)構(gòu)

在資源充足的情況下，二維脈動陣列分解結(jié)構(gòu)是性能最高的結(jié)構(gòu)[14]。二維陣列采用按堆積資源的方式處理數(shù)據(jù)更新，數(shù)據(jù)吞吐量和運算性能較高，其具體結(jié)構(gòu)如圖3所示。

圖3 二維陣列結(jié)構(gòu)圖

二維陣列由對角線單元和非對角線單元組成。對角線單元由一個Givens生成(Givens Generation，GG)單元和一個Givens旋轉(zhuǎn)(Givens Rotation，GR)單元組成，用于計算更新因子和行列更新。非對角單元由一個GR單元組成，僅用于行列數(shù)據(jù)更新。對角線處理單元按式(3)計算每一組Givens的旋轉(zhuǎn)因子并向右傳遞，同時更新該列的數(shù)據(jù)，并將數(shù)據(jù)回存至存儲器中。非對角線處理單元使用左側(cè)處理單元傳遞的旋轉(zhuǎn)因子，將該列的當(dāng)前行元素更新并向下傳遞。

二維脈沖陣列結(jié)構(gòu)通過堆積資源的方式實現(xiàn)性能提升，在處理一些低階矩陣時擁有較強的運算性能。但是隨著矩陣階數(shù)的增加，運算資源也呈現(xiàn)指數(shù)增長的趨勢。在高精度計算情況下，更新因子的計算涉及開方和除法運算，硬件消耗巨大。綜上所述，二維結(jié)構(gòu)雖然自身性能優(yōu)異，但不適用于處理高精度高緯度的矩陣求逆運算。

3 優(yōu)化的CGR-RD分解結(jié)構(gòu)

3.1 基于CGR-RD算法的一維分解結(jié)構(gòu)

由上文可知，二維陣列每次消元時，數(shù)據(jù)從上一級向下傳輸。當(dāng)所有矩陣輸入到陣列時，運算并行度達到最高峰，此時參與計算的對角運算單元數(shù)量為n/2，非對角單元數(shù)量為(3n3-2n)/4。在浮點運算過程中，二維陣列的硬件資源消耗巨大，且資源閑置率較高，因此不適用于實際應(yīng)用?；谝陨显颍疚奶岢鲆环N適用于2～32階高精度矩陣分解的一維線性結(jié)構(gòu)，具體結(jié)構(gòu)如圖4所示。

圖4 CGR-QR分解一維線性結(jié)構(gòu)

CGR-QR分解結(jié)構(gòu)僅使用一組更新因子計算單元，數(shù)據(jù)以流水形式傳入更新因子計算單元。其采用8組行列更新計算單元并行計算，同時對更新因子計算與行列更新計算的時間進行折疊處理?；贑GR-RD算法的GG單元和GR單元的結(jié)構(gòu)如圖5所示。

(a)

由式(6)可知，更新因子的計算僅需矩陣中的一列元素，因此當(dāng)矩陣更新完首列元素后，便可并行計算下一次消元的更新因子。如圖6所示，更新因子計算部分分為一個GG單元和兩組寄存器組。當(dāng)GG單元完成第一組更新因子的計算后，將數(shù)據(jù)寫入第一組寄存器組中，第二組寄存器組從上一個寄存器組中讀取更新因子，并將其傳入GR單元進行行列更新。當(dāng)GR單元完成首列元素的更新操作后，GG單元從存儲器讀取更新后的列元素，并進行下一次行列更新的更新因子計算。當(dāng)計算完成后，將數(shù)據(jù)寫入第1組寄存器。對于n階矩陣，GG單元和GR單元運行的時間對比如圖6所示。

圖6 n階矩陣GG、GR單元運算時間對比圖

由圖6可知，基于CGR-QR的一維線性結(jié)構(gòu)將GG運算單元與GR運算單元的運算時間折疊，將原本需要n-1次GG單元和n-1次GR單元的運算壓縮到僅需1次GG單元和n-1次GR單元運算。

相比于二維陣列結(jié)構(gòu)，該結(jié)構(gòu)的性能損失主要為：當(dāng)矩陣階數(shù)變大時，由于GR單元的數(shù)目要小于矩陣階數(shù)，一維線性結(jié)構(gòu)在行列向量消元的過程中需要[N/4]個運算周期進行消元運算，[i]為向上取整函數(shù)。

GR單元的數(shù)量決定了該一維線性結(jié)構(gòu)的運算性能，因此加大GR單元的數(shù)量可以提高分解的效率。然而當(dāng)GR單元數(shù)量過多時，當(dāng)前階段的數(shù)據(jù)更新會快于下一階段更新因子的計算，GR單元仍然需要等待更新因子計算，并且當(dāng)GR數(shù)量增大時，硬件資源消耗也將逐漸加大。因此，需要根據(jù)更新因子的計算時間和數(shù)據(jù)更新的計算時間等條件計算出最合適GR單元的數(shù)目，計算式為

(8)

式中，TGG為更新因子計算單元的計算時間；TGR為行列更新計算單元的計算時間；M為矩陣運算器兼容的最大階數(shù)；N為GR單元的數(shù)量，由于本文的Q和R矩陣是同時得到的，因此N為偶數(shù)，其計算式如式(9)所示。

(9)

根據(jù)設(shè)計指標(biāo)和實現(xiàn)工藝的情況，本文GR單元的數(shù)量N被選定為8，即8組GR單元。

3.2 一維線性結(jié)構(gòu)分析

基于CGR-QD算法的一維線性結(jié)構(gòu)與二維脈沖陣列的理論運算周期數(shù)和資源消耗分別如表1、表2所示。

表1 兩種結(jié)構(gòu)的理想運算周期數(shù)

表2 兩種結(jié)構(gòu)的浮點運算器資源理論消耗情況

由于一維陣列的運算資源數(shù)和二維陣列的運算資源數(shù)不同，故不能單純地比較兩者的運算周期數(shù)。本文在一維陣列與二維陣列的運算資源數(shù)對等的情況下，選取兩者計算單個數(shù)據(jù)的時間進行比對，結(jié)果如圖7所示。

圖7 二維陣列結(jié)構(gòu)和一維線性結(jié)構(gòu)性能對比圖

由圖7可知，在矩陣階數(shù)小于6階時，二維陣列的性能高于一維陣列結(jié)構(gòu)。由于二維陣列的運算資源與矩陣階數(shù)的平方成正比，因此當(dāng)矩陣階數(shù)逐漸增加時，二維陣列處理單個數(shù)據(jù)的時間逐漸多于一維線性結(jié)構(gòu)。根據(jù)表1可知，在矩陣階數(shù)小于6階時，一維陣列與二維陣列的運算周期數(shù)相差較少。綜上所述，在高速高精度矩陣求逆器設(shè)計中采用基于CGR的一維線性結(jié)構(gòu)所耗費的運算資源少，運算性能較高，能更好地滿足設(shè)計需求。

4 硬件驗證與性能分析

4.1 驗證方案

本文采用Verilog語言實現(xiàn)硬件設(shè)計，在TSMC 28 nm工藝下，使用Verdi_VL-2016.06-1與VCS-VL-2016.06對本文設(shè)計進行仿真驗證，并通過Design Complier軟件完成綜合實現(xiàn)，最終通過IC-Complier后端流程得到矩陣分解器的專用集成電路(Application Specific Integrated Circuit，ASIC)版圖。矩陣分解器的ASIC版圖如圖8所示。

圖8 矩陣分解器的ASIC版圖

本文設(shè)計采用Synopsys公司DesignWare庫中的浮點流水IP，包括浮點加、減、乘、除和開方。經(jīng)過綜合和后端時序調(diào)整，設(shè)置浮點加法器、浮點乘法器為4級流水，浮點除法器為32級流水，浮點開方器為25級流水，綜合考慮運算性能得到GR單元數(shù)目為4。本文在分解過程中同時進行Q和R陣的運算，因此GR單元總數(shù)為8，雙口RAM的數(shù)量為8。

4.2 資源和誤差分析

本文設(shè)計的基于CGR-QR的分解器在現(xiàn)場可編程門陣列(Field Programmable Gate Array，F(xiàn)PGA)上進行了驗證，F(xiàn)PGA開發(fā)板型號為Xilinx XCV440T，表3為FPGA資源使用情況。

表3 FPGA的硬件資源占用表

本設(shè)計通過MATLAB軟件在[-100，100]、[-1010，1010]、[-1020，1020]3個范圍內(nèi)隨機取數(shù)，數(shù)值精度為雙精度浮點。本文分別測試了4、8、16、32階非奇異矩陣。MATLAB求逆結(jié)果與設(shè)計結(jié)果誤差如表4所示。

表4 設(shè)計結(jié)果與MATLAB結(jié)果誤差

由表4可知，設(shè)計結(jié)果和MATLAB計算結(jié)果的相對誤差穩(wěn)定在10-12～10-15數(shù)量級，滿足高精度的設(shè)計要求。浮點計算器由于精度的限制，在進行計算時會發(fā)生截斷誤差，這是誤差的主要來源。由于矩陣階數(shù)相對較小，因此累計產(chǎn)生的誤差不大。

4.3 性能分析

不同結(jié)構(gòu)下，矩陣從開始分解到分解結(jié)束的運算周期數(shù)如表5所示。

表5 3種結(jié)構(gòu)的運算周期比較

相比于傳統(tǒng)一維結(jié)構(gòu)，本文采用基于CGR-RD的一維線性結(jié)構(gòu)，將分解過程中的消元步驟并行化，降低了計算時間。相比二維結(jié)構(gòu)，本文結(jié)構(gòu)的運算周期略大于二維結(jié)構(gòu)，但本文設(shè)計的分解器在處理高階矩陣時運算資源消耗較小，且可兼容2～32階矩陣分解，適用性更好。

另外，本文將所設(shè)計的分解器與NVIDIA RTX2070 GPU的CUDA軟件進行了對比，兩者的矩陣分解周期如表6所示。

表6 本文分解器與NVDIA GPU的矩陣分解周期對比

表6表明，矩陣規(guī)模在32階以下時，分解器加速效果較好，當(dāng)矩陣階數(shù)達到32階時，和RTX 2070 GPU相比，本文設(shè)計的分解器相對于RTX 2070的分解周期具有22.8倍加速比。加速的主要原因有兩點：(1)高度并行的CGR-QR分解算法；(2)硬件實現(xiàn)中采用一維線性結(jié)構(gòu)，縮短了運算時間。

5 結(jié)束語

本文在Givens-QR算法的基礎(chǔ)上，根據(jù)文獻[18]的分解過程，推導(dǎo)出CGR-QR算法中R矩陣和Q矩陣的計算式，并設(shè)計了基于CGR-QR算法的矩陣分解器。該分解器采用4路并行的流水結(jié)構(gòu)進行消元計算，減少了運算資源，降低了矩陣分解的運算周期。此外，本文對設(shè)計的矩陣分解器進行FPGA和ASIC驗證。結(jié)果表明，本文的矩陣分解器可以在TSMC 28 nm工藝下，以700 MHz主頻工作，其運算結(jié)果與MATLAB運算結(jié)果的相對誤差在10-12以下。當(dāng)矩陣階數(shù)達到12階時，該分解器相比于傳統(tǒng)一維分解結(jié)構(gòu)達到了2.8倍的加速比。當(dāng)矩陣階數(shù)為32階時，分解器的周期數(shù)相比于RTX2070達到22.8倍的加速比。綜上，本文設(shè)計的矩陣分解器滿足高速、高精度的設(shè)計需求，滿足現(xiàn)代信號處理的要求。