邱俊豪,宋宇鯤,陳文杰,侯 寧

(1.合肥工業(yè)大學(xué) 微電子設(shè)計研究所,安徽 合肥 230601; 2.合肥工業(yè)大學(xué) 教育部IC設(shè)計網(wǎng)上合作研究中心,安徽 合肥 230601;3.河南城建學(xué)院 電氣與控制工程學(xué)院,河南 平頂山 467000)

0 引 言

矩陣分解是線性代數(shù)中最重要的運算之一,在諸如MIMO系統(tǒng)[1-2]、自適應(yīng)濾波[3]和智能機器人[4]等領(lǐng)域有著廣泛的應(yīng)用。QR分解是矩陣分解的常用方法,典型的QR分解[5]有Gram-Schmidt正交法、Householder變換(Householder rotation,HT)和Givens旋轉(zhuǎn)法(Givens rotation,GR)。由于HT變換的運算復(fù)雜度相對較高,工程實踐中通常使用另外2類算法。一類是基于修正的Gram-Schmidt正交[6]的矩陣分解算法MGS,該算法的運算并行度較高,但數(shù)值穩(wěn)定性相對較差。另一類是基于GR旋轉(zhuǎn)[7-11]的矩陣分解算法,該算法具有多粒度并行的特點,易于硬件實現(xiàn),且在高精度計算時,GR算法的數(shù)值穩(wěn)定性優(yōu)于MGS算法。

自文獻[7]提出GR算法以來,國內(nèi)外學(xué)者進行了一系列的改進。文獻[8]使用牛頓迭代法和查找表法替代GR算法中的平方根倒數(shù),實現(xiàn)高速近似運算;文獻[9]使用二維陣列結(jié)構(gòu)實現(xiàn)單精度浮點矩陣的GR-QR分解;文獻[10]采用GR分解的改進算法減少算法運算量,并采用二維陣列結(jié)構(gòu)實現(xiàn);文獻[11]使用基于CORDIC算法的流水線結(jié)構(gòu)實現(xiàn)高速的GR-QR分解。已知的GR算法的硬件實現(xiàn)集中在非高精度矩陣的高速實現(xiàn)上,關(guān)于高精度矩陣分解電路實現(xiàn)的研究較少。

在大數(shù)據(jù)時代,復(fù)雜信號處理要求能實時提供高維度、高精度的矩陣求解結(jié)果[12]。因此,本文在深入研究GR-QR算法的基礎(chǔ)上,提出一種高性能GR-QR分解電路結(jié)構(gòu),該結(jié)構(gòu)充分利用GR-QR分解算法中存在的向量計算和元素計算并行執(zhí)行特征,設(shè)計多路并行計算通道配套的地址無沖突存儲策略,實現(xiàn)高計算速度、高運算資源利用率。

1 算法原理

設(shè)A為m×n階非奇異矩陣,Q為m×m階的正交矩陣,R為m×n階上三角矩陣,若滿足下式則稱其為矩陣A的正交三角(QR)分解。

A=QR

(1)

QR分解中,Givens旋轉(zhuǎn)通過旋轉(zhuǎn)矩陣G左乘原矩陣A,更新原矩陣第i行以及第j行數(shù)據(jù),其他行不變[13],即

(2)

其中,c、s為Givens旋轉(zhuǎn)的旋轉(zhuǎn)因子,表達式如下:

j=i+1,i=1,2,…,n

(3)

(4)

由(4)式可得:

Q=(Gn-1,n,…,G12)-1,

由文獻[13]可知,Q為正交矩陣,R為上三角矩陣,上述過程稱為矩陣A的GR-QR分解。

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

2.1 GR-QR分解的混合粒度并行

假設(shè)A為n階矩陣,其行向量可記為a1、a2、…、an。若對矩陣A進行GR-QR分解,相當于對相鄰的向量依次進行數(shù)據(jù)消元。一組相鄰的行向量ai、ai+1的數(shù)據(jù)消元記作一次消元任務(wù),分2步進行,首先計算旋轉(zhuǎn)因子c、s,然后根據(jù)旋轉(zhuǎn)因子更新向量ai、ai+1。GR-QR分解過程中的消元任務(wù)可進行粗細粒度并行化處理。其中：細粒度并行是指行向量的各列元素相互獨立,可進行同時計算;粗粒度并行是指多組不同的相鄰向量可同時進行。結(jié)合粗細粒度并行,4階矩陣的GR-QR分解過程如圖1所示。

細粒度并行如圖1階段2所示,3組列向量元素a32、a42,a33、a43,a34、a44同時進行數(shù)據(jù)更新;粗粒度并行如圖1階段3所示,a1、a2,a3、a42組不相關(guān)的相鄰行向量同時開始旋轉(zhuǎn)因子計算。除了顯而易見的粗細粒度并行,GR-QR分解在相鄰任務(wù)的計算上也存在并行化可能。盡管由于數(shù)據(jù)相關(guān)性,相鄰任務(wù)的旋轉(zhuǎn)因子計算不能并行,但在一定條件下,相鄰任務(wù)的數(shù)據(jù)更新與旋轉(zhuǎn)因子計算可同時執(zhí)行。

圖1 粗細粒度并行的GR-QR分解步驟

如圖1階段3所示,a2、a3的旋轉(zhuǎn)因子計算完成并更新對應(yīng)行向量的前2列a21、a31,a22、a32后,a1、a2,a3、a4的旋轉(zhuǎn)因子計算開始運行,無需等待a23、a33,a24、a34更新完畢。上述過程稱為GR-QR分解的混合粒度并行過程[14],具體的任務(wù)調(diào)度如圖2所示。

圖2 混合粒度并行的4階矩陣GR-QR分解任務(wù)調(diào)度

理想情況下,每個階段的GR-0任務(wù)內(nèi)完成該階段的所有數(shù)據(jù)更新,相當于每個階段僅需1次旋轉(zhuǎn)因子計算時間和1次數(shù)據(jù)更新時間,由此得到n階矩陣的理論運算時間,即

Ttotal=(Tupdate+Tgivens)(2n-3)

(5)

其中:Ttotal為總運算周期;Tupdata為矩陣更新時間;Tgivens為旋轉(zhuǎn)因子計算時間。

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

理想硬件資源下,二維陣列結(jié)構(gòu)實現(xiàn)GR-QR分解的效率最高[9-10]。分解器有極高的數(shù)據(jù)吞吐量和運算性能,如圖3所示。

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

圖3中的二維陣列結(jié)構(gòu)包括對角單元與非對角單元。對角單元(圖3中顯示為橢圓)按(3)式計算輸入矩陣對應(yīng)列的旋轉(zhuǎn)因子c、s,生成的旋轉(zhuǎn)因子被廣播至同一行的所有非對角單元(圖3中顯示為矩陣)。非對角單元根據(jù)旋轉(zhuǎn)因子對輸入的原矩陣數(shù)據(jù)進行更新,即每行非對角單元根據(jù)旋轉(zhuǎn)因子計算結(jié)果對涉及的兩行矩陣進行消元操作。

二維陣列結(jié)構(gòu)的工作方式如下:原矩陣A=[a1，a2,…,an]T的行向量從an開始依次從頂部輸入二維陣列。當接收到第1組數(shù)據(jù)an、an-1時,對角單元PE1,1計算該行的旋轉(zhuǎn)因子,非對角單元PE1,2,PE1,3,…,PE1,n根據(jù)PE1,1的計算結(jié)果更新矩陣數(shù)據(jù),an的更新結(jié)果發(fā)送給下一層的對角單元與非對角單元,an-1留在當前位置,待an-2發(fā)送過來后開始下一階段數(shù)據(jù)消元。

依次類推,數(shù)據(jù)流從上往下傳遞,當所有矩陣數(shù)據(jù)均通過二維陣列,每個PE單元保存的值即為三角矩陣R的輸出值。上述過程嚴格遵守混合粒度并行的GR-QR分解的計算流程。圖2每個階段的每個旋轉(zhuǎn)因子計算和數(shù)據(jù)更新在不同的運算單元內(nèi)執(zhí)行,保證每個階段在1次旋轉(zhuǎn)因子計算時間和1次數(shù)據(jù)更新時間內(nèi)完成,運算性能達到理論運算周期。

二維結(jié)構(gòu)以資源全展開的方式獲得理論性能,其所需資源消耗是矩陣階數(shù)的2次冪,且在高精度計算場景,除法、開方運算無法轉(zhuǎn)換為近似運算[8,11]形式,硬件開銷過于巨大。綜上所述,受資源規(guī)模限制,二維結(jié)構(gòu)不適用于高精度矩陣分解實現(xiàn)。

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

3.1 基于原位替換的一維線性結(jié)構(gòu)

由上文可知,二維陣列的第k層先后通過向量ai、ai+1和ai-1、ai進行2次消元,向第k+1層傳遞更新后的向量ai+1和ai,接著第k層計算ai-2、ai-1的消元,第k+1層計算ai+1、ai的消元。原矩陣數(shù)據(jù)按上述規(guī)律向下傳遞,當所有矩陣均輸入至二維陣列時,運算并行度達到峰值,此時僅有n/2個對角單元以及3/4n2-n/2個非對角單元處于工作狀態(tài)。

在GR-QR分解的多數(shù)時間,二維陣列的大部分處理單元處于閑置狀態(tài),資源利用率極低。在浮點矩陣計算時,浮點運算器的資源消耗巨大,高閑置率是硬件實現(xiàn)時無法接受的。

針對浮點矩陣分解需求,本文從2個方面減少運算資源消耗,提高資源利用率。對于旋轉(zhuǎn)因子計算,利用流水線結(jié)構(gòu)特征,使用一個旋轉(zhuǎn)因子計算單元代替二維結(jié)構(gòu)的n個對角單元,矩陣數(shù)據(jù)依次流水輸入至旋轉(zhuǎn)因子計算單元,計算時間僅增加少量流水線填充時間。

對于矩陣數(shù)據(jù)更新,按圖1階段3所示,a23、a33和a24、a34的數(shù)據(jù)更新與向量a1、a2和a3、a4的旋轉(zhuǎn)因子計算并行,不必進行二維陣列類似的飽和計算,僅需設(shè)置合適的數(shù)據(jù)更新單元數(shù)目,保證第n階段的數(shù)據(jù)更新在第n+1階段的旋轉(zhuǎn)因子計算結(jié)束前完成。綜上所述,本文設(shè)計了基于原位替換的GR-QR分解一維線性結(jié)構(gòu),如圖4所示。

圖4 基于原位替換的一維線性結(jié)構(gòu)

圖4中,GG單元用于旋轉(zhuǎn)因子計算,GR單元用于矩陣數(shù)據(jù)更新。在旋轉(zhuǎn)因子計算時,當前階段所需的矩陣數(shù)據(jù)存儲于數(shù)據(jù)寄存器中,依次流水填充至GG單元內(nèi),計算結(jié)果存儲于旋轉(zhuǎn)因子寄存器中;在數(shù)據(jù)更新時,GR單元讀取旋轉(zhuǎn)因子寄存器中的旋轉(zhuǎn)因子,對本地存儲的矩陣數(shù)據(jù)進行數(shù)據(jù)更新,計算結(jié)果替換原矩陣數(shù)據(jù)。

一維線性結(jié)構(gòu)的執(zhí)行過程如下:

(1) 數(shù)據(jù)寄存器從本地存儲Mem中讀取初始數(shù)據(jù)an-1,1、an,1,發(fā)送給GG單元,計算an-1、an的第1組旋轉(zhuǎn)因子;

(3) 數(shù)據(jù)寄存器接收到矩陣數(shù)據(jù)后,GG單元激活,開始計算an-1、an-2的旋轉(zhuǎn)因子,此時an-1、an-2的旋轉(zhuǎn)因子計算與an-1、an的數(shù)據(jù)更新同時進行。

(4) 依次類推,一維線性結(jié)構(gòu)的GR-QR分解分為2個過程交替進行,分別為GG單元與GR單元均激活時的旋轉(zhuǎn)因子與矩陣數(shù)據(jù)并行計算過程以及僅GR單元激活的矩陣數(shù)據(jù)更新過程。此外,在并行計算過程,旋轉(zhuǎn)因子寄存器需存儲2個階段的旋轉(zhuǎn)因子,因此需設(shè)置2組寄存器,以乒乓操作形式執(zhí)行運算。

本文設(shè)計需設(shè)置合適的GR單元數(shù)目,使得矩陣數(shù)據(jù)更新在旋轉(zhuǎn)因子計算結(jié)束前完成。GR單元數(shù)目的計算公式為:

(6)

其中:TGG為GG單元的流水級數(shù);TGR為GR單元的流水級數(shù);M為矩陣最大兼容階數(shù);N為GR單元數(shù)目。由(6)式可推得：

(7)

3.2 基于運算并行度的無沖突存儲規(guī)則

(7)式成立的前提條件是所有GR單元在數(shù)據(jù)更新完成前一直處于工作狀態(tài),因此本文采用集中資源逐行計算的矩陣數(shù)據(jù)更新方式。單個GR單元對應(yīng)雙讀雙寫,n個GR單元最少需要2n個簡單雙口RAM。在更新過程中,數(shù)據(jù)讀寫應(yīng)充分利用存儲端口同時避免地址沖突。若數(shù)據(jù)存儲按行存儲,則每1行的端口數(shù)僅能支持1個GR單元運算;若按列存儲,GR單元數(shù)據(jù)更新需要的數(shù)據(jù)位于同1片RAM內(nèi),產(chǎn)生讀寫地址沖突。

表1 運算并行度為n的M階矩陣存儲規(guī)則

4 硬件驗證與性能分析

4.1 驗證方案

本文采用VerilogHDL語言實現(xiàn)硬件設(shè)計,在TSMC28 nm工藝下,使用verdi-vL-2016.06-1與vcs-vL-2016.06對設(shè)計進行仿真驗證,然后通過Design Complier軟件完成綜合實現(xiàn),通過ICC后端流程最終得到矩陣分解器的ASIC版圖,如圖5所示。

圖5 矩陣分解器的ASIC版圖

設(shè)計采用Synopsys公司的DesignWare庫中的浮點流水IP,包括了浮點加、減、除以及開方。經(jīng)過綜合與后端的時序調(diào)整,設(shè)置浮點加法為12級流水、浮點乘法為8級流水、浮點除法為27級流水以及浮點開方為45級流水,綜合考慮運算性能與運算面積結(jié)合(7)式,得到GR單元數(shù)目為4。本文的矩陣分解器支持R矩陣和Q矩陣同時求解,因此GR單元的總數(shù)為8,雙口RAM的數(shù)量為16。

后端設(shè)計結(jié)果表明,在TSMC28 nm工藝下,本文設(shè)計在700 MHz主頻下完成2～32階雙精度浮點矩陣的矩陣分解,芯片面積為2 mm2,利用率為63%,折合實際面積為1.26 mm2。

4.2 精度分析

除了ASIC流程外,使用FPGA進行基于硬件的原型驗證是ASIC設(shè)計的重要一環(huán)。本文在Xlinx XC7V2000T上進行了分解器的硬件性能測試,FPGA的硬件資源使用情況,見表2所列。

表2 FPGA資源使用情況

通過Matlab在[-10,10]、[-10100,10100]、[-10200,10200]3個取值范圍內(nèi)隨機取數(shù),數(shù)值精度為雙精度浮點,分別構(gòu)成4、8、16、32階非奇異矩陣。每個取值范圍內(nèi)取1 000組矩陣,對ASIC仿真結(jié)果,FPGA驗證結(jié)果與Matlab軟件結(jié)果的誤差分析，見表3所列,以Matlab結(jié)果作為參考值。

表3 不同階數(shù)下ASIC仿真與FPGA驗證結(jié)果誤差

由表3可知,輸入數(shù)據(jù)的數(shù)量級影響輸出結(jié)果的數(shù)量級,若ASIC與FPGA的均方差按照浮點數(shù)的尾數(shù)計算,在10-14～10-15數(shù)量級,滿足高精度的設(shè)計要求。

4.3 性能分析

與文獻[9]單精度浮點數(shù)的二維陣列結(jié)構(gòu)相比,計算單元歸一化后的運算周期,見表4所列。

表4 與其他結(jié)構(gòu)的運算周期比較

相比于串行執(zhí)行的傳統(tǒng)一維結(jié)構(gòu),本文優(yōu)化的一維結(jié)構(gòu)基于混合粒度并行算法,具有較好的加速。相比于二維陣列結(jié)構(gòu),本文的運算周期稍大,但本文的資源消耗顯然小于二維結(jié)構(gòu)。

由于工藝、適用范圍等因素不同,同類型設(shè)計的直接比較較為困難,本文將不同階的矩陣分解分別在硬件設(shè)計和NVIDIA GPU的CUDA軟件實現(xiàn)進行對比測試,得到單位面積的性能比如圖6所示。

圖6 硬件設(shè)計與GPU單位面積性能比

圖6表明,矩陣規(guī)模越大,加速越明顯,當矩陣階數(shù)達到32階時,和RTX2070 GPU相比,本文設(shè)計的單位性能比可達到95倍。加速的主要原因如下:

(1) GR-QR分解的混合粒度并行。

(2) 硬件實現(xiàn)中采用一維線性結(jié)構(gòu),減少運算資源消耗,提高運算資源利用率。

5 結(jié) 論

本文基于GR-QR分解的混合粒度并行算法設(shè)計實現(xiàn)用于高精度矩陣分解的硬件電路,提出基于原位替換的一維線性結(jié)構(gòu),以流水線形式運行GR-QR分解,減少運算資源開銷,提高資源利用率,并設(shè)計基于運算并行度的無沖突存儲規(guī)則,充分利用存儲端口,提高運算并行度。在ASIC平臺完成硬件實現(xiàn),實驗結(jié)果表明,在TSMC 28 nm工藝下,本文設(shè)計可在700 MHz主頻下完成2～32階雙精度浮點矩陣的矩陣分解,芯片面積為2 mm2,運算精度可達10-15,單位面積性能是RTX 2070的95倍,達到高速、高精度和低資源消耗的目標,滿足目標SoC芯片的要求。