一種用于矩陣求逆的原位替換算法及硬件實現(xiàn)

張多利, 蔣 雯, 葉紫燕, 宋宇鯤, 汪 健

(1.合肥工業(yè)大學 電子科學與應(yīng)用物理學院,安徽 合肥 230601; 2.中國兵器工業(yè)集團 第214研究所,安徽 蚌埠 233000)

0 引 言

隨著數(shù)字信號處理技術(shù)的發(fā)展,矩陣求逆[1]運算的應(yīng)用領(lǐng)域日益廣泛。人工智能技術(shù)[2]、圖形圖像處理、雷達技術(shù)以及實時通信技術(shù)等都涉及了大量的矩陣求逆運算[3],而這些工程領(lǐng)域都有高吞吐量、實時性運算要求,如何利用有效資源實現(xiàn)快速大規(guī)模矩陣求逆運算成為了硬件設(shè)計的重點和難點。

傳統(tǒng)的矩陣求逆方法,如初等變換法、伴隨矩陣法、Gauss-Jordan消去法[4]等,由于其固有的數(shù)學屬性而不適于硬件實現(xiàn)，工程上常用的矩陣求逆方法有Cholesky分解求逆法[5]和QR分解求逆法[6-7]。Cholesky分解求逆法計算量不大,但其是對對稱正定矩陣進行三角分解,適用矩陣范圍太小。QR分解求逆法運算量較大,中間有大量數(shù)據(jù)相關(guān)的迭代過程,還涉及大量開方和范數(shù)運算,運算時間較長,并行度不高,資源消耗量太大。

為了平衡計算效率和資源消耗,本文利用基于LU分解的按位替換算法[8]對矩陣進行求逆。該算法運算算式相對簡單,不需要大量的開方、求范數(shù)等復雜運算;且該算法與利用LU分解[9]進行矩陣求逆的常規(guī)算法相比,可以直接并行求三角逆矩陣,具有優(yōu)越的可并行性,降低了運算時間復雜度?；谠撍惴ǖ挠布軜?gòu)在源數(shù)據(jù)存儲空間上進行求逆運算和結(jié)果存儲,無需占用額外存儲資源,實現(xiàn)了運算資源和存儲資源的復用,具有優(yōu)越的可應(yīng)用性。

1 算法分析

設(shè)所有順序主子式不為0的M階方陣矩陣A為:

其中,aij為第i行第j列元素(i,j=1,2,3,…,M)。采用按位替換法進行矩陣求逆,需要先對待求逆矩陣A進行約化系數(shù)求解,得到M階約化系數(shù)方陣N。N的第1行和第1列元素為:

(1)

其中,i,j=1,2,3,…,M。N的其他元素為:

(2)

其中,i=2,3,…,M;j=i+1,i+2,…,M。最后得到約化系數(shù)矩陣N,即

由N可直接并行求得待求逆矩陣A分解所得的下三角陣L的逆矩陣L-1和上三角矩陣R的逆矩陣R-1。由L-1左乘R-1可得到矩陣A的逆矩陣,即A-1=R-1L-1。

L-1的元素lii為:

(3)

其中,i=1,2,3,…,M。L-1的元素lji為:

(4)

其中,i=1,2,3,…,M-1;j=i+1,i+2,…,M。則L-1為:

R-1的元素rii為:

(5)

其中,i=1,2,3,…,M。R-1的元素rji為:

(6)

其中,j=1,2,3,…,M-1;i=j+1,j+2,…,M。則R-1為:

按原位替換算法進行矩陣求逆的算法原理與傳統(tǒng)矩陣求逆算法相比,避免了傳統(tǒng)算法中大量的開方、平方、求范數(shù)等復雜運算;通過創(chuàng)建約化系數(shù)矩陣,沿對角線對稱的上三角矩陣元素和下三角矩陣元素可以并行計算,避免了傳統(tǒng)算法中大量的串行迭代過程,運算時間得到壓縮,并行度有所提高;采用按位替換的方式,在硬件設(shè)計中所有數(shù)據(jù)都存在原存儲空間中,減少了存儲資源消耗。

2 地址無沖突設(shè)計

基于原位替換的矩陣求逆算法,即整個運算過程的數(shù)據(jù)更替都在原矩陣所占的存儲空間內(nèi)進行。由于整個運算過程都是采用高效并行計算,矩陣源數(shù)據(jù)、中間數(shù)據(jù)和結(jié)果數(shù)據(jù)的讀取順序和存儲方式直接決定著整個設(shè)計的運算高效性、正確性及數(shù)據(jù)讀取的安全性。因此,讀取順序和存儲方式的設(shè)計以及無沖突地址的生成是整個矩陣求逆設(shè)計的重中之重。

本文先進行約化系數(shù)矩陣計算,再對矩陣進行數(shù)據(jù)預處理和三角逆矩陣計算,最后是上、下三角逆矩陣乘運算。其中數(shù)據(jù)的預處理過程為接下來的三角逆矩陣運算作了簡化,避免了三角逆過程中由于數(shù)據(jù)存取方式變化帶來的地址讀寫沖突。運算所需存儲器中,有1組塊存儲器Memory0由8小塊1 024×64 bit的片上隨機存取存儲器(random access memory,RAM)組成;另3組存儲器dist-mem0、dist-mem1和dist-mem2均由8小塊16×64 bit小的分布式RAM組成。運算數(shù)據(jù)采用32位浮點數(shù)形式,如此,對于1個小塊存儲器來說,讀取1個地址可以得到2個數(shù)據(jù);對1個大塊存儲器來說,讀取1個地址可以得到16個數(shù)據(jù)。這樣在整個運算過程中,多組源數(shù)據(jù)可從2組存儲器中并行讀取,避免了由于高并行度緊湊計算帶來的數(shù)據(jù)冒險。

2.1 存儲資源管理分配與約化系數(shù)矩陣計算

在源數(shù)據(jù)存儲模式下,將上、下三角元素沿對角線“對折”,第i行第j列元素和第j行第i列元素兩兩拼接,分別作為高32位和低32位,放在1個存儲地址空間內(nèi),每讀1個地址,能輸出2個元素,避免了地址沖突。其中對角元素存入dist-mem0,非對角元素存入Memory0。

對于前M/2行(列)非對角元素,每8行(列)為1個循環(huán),每行(列)的第1個元素按規(guī)律存入memx的某一地址,同行(列)的前8-x個元素按順序存入memx至mem7的同一地址,剩余元素每8個按順序存入mem0至mem7的同一地址,其中memx為Memory0的組內(nèi)存儲器號。

對于后M/2行(列)非對角元素,每8行(列)為1個循環(huán),每行(列)的第1個元素按順序存入memx至mem7的同一地址,同1行(列)的前8-x個元素按規(guī)律存入memx的某一地址,剩余元素每8個按順序存入mem0至mem7的某一地址。對于對角元素,每8行(列)為1個循環(huán),順次存儲在dist-mem0中。

以8階矩陣求逆存儲方式為例,8階矩陣非對角元素無沖突存儲方式見表1所列,矩陣對角元素無沖突存儲方式見表2所列。

表1 8階矩陣非對角元素無沖突數(shù)據(jù)存儲方式

表2 8階矩陣對角元素無沖突數(shù)據(jù)存儲方式

約化系數(shù)計算模式下數(shù)據(jù)存取步驟如下:

(1) Memory0輸出當前行需要循環(huán)取的數(shù)據(jù),存入dist-mem1。

(2) Memory0和dist-mem1分別提供約化系數(shù)計算單元(processing element,PE)PE1中乘法器的2組源數(shù)據(jù)。

(3) dist-mem0向PE1中的除法器提供源數(shù)據(jù)。

(4) Memory0提供減法器所需要的1組源數(shù)據(jù)。

(5) 將非對角元素結(jié)果回寫至Memory0中,仍保持上、下三角元素合為1個地址的存儲方式,其數(shù)據(jù)地址與源數(shù)據(jù)地址一致。

(6) 在約化系數(shù)對角元素結(jié)果計算結(jié)束后,由三角逆矩陣對角元素計算單元PE2隨后計算出逆矩陣的對角元素,并存入dist-mem2中。

以4階矩陣進行矩陣求逆為例,數(shù)據(jù)搬運原理如圖1所示。

圖1 約化系數(shù)矩陣數(shù)據(jù)搬運原理示意圖

圖1中,采用原位替換的矩陣求逆算法,約化系數(shù)矩陣沿對角線對稱的上三角矩陣元素和下三角矩陣元素可以并行計算;未用箭頭標出的元素與已標出元素執(zhí)行相同的操作。

2.2 預處理和三角逆矩陣計算

在求得約化系數(shù)矩陣后對數(shù)據(jù)進行2步預處理,避免了由于地址讀寫復雜帶來的沖突,為求上、下三角逆矩陣作了簡化。在數(shù)據(jù)預處理計算模式下的數(shù)據(jù)存取步驟如下:

(1) Memory0提供乘法器所需的一組數(shù)據(jù)njk·(k-1)或nki·(k-1),dist-mem2提供乘法器所需要的另一組數(shù)據(jù)lkk或rkk,并將結(jié)果回寫至Memory0中,其數(shù)據(jù)地址與源數(shù)據(jù)地址一致。

(2) Memory0提供乘法器所需的一組數(shù)據(jù)njk·(k-1)·lkk或者nki·(k-1)·rkk,dist-mem2提供乘法器所需的另一組數(shù)據(jù)ljj或rii,并將結(jié)果回寫至Memory0中,其數(shù)據(jù)地址與源數(shù)據(jù)地址一致。

在上、下三角逆矩陣計算模式下,Memory0提供上、下三角逆矩陣計算單元PE3中乘法器的2組源數(shù)據(jù),并作為三角逆矩陣結(jié)果存儲模塊,在相應(yīng)周期后將數(shù)據(jù)回寫至Memory0,期間為保證數(shù)據(jù)存取有效,將中間結(jié)果緩存至dist-mem0。三角逆矩陣數(shù)據(jù)搬運原理如圖2所示。圖2中,矩陣中每條虛線上的元素運算時間相同,運算順序為先計算靠近對角元的虛線,再沿著2個頂點逐漸向外擴散進行運算;未用箭頭標出的元素與已標出元素執(zhí)行相同的操作。

圖2 三角逆矩陣數(shù)據(jù)搬運原理示意圖

2.3 上、下三角逆矩陣乘運算

在矩陣乘計算模式下,采用直接矩陣乘的方式,Memory0組內(nèi)的8組mem采用一端讀數(shù)一端寫數(shù)的模式,設(shè)置2組乘加PE,每8個乘加器為1組PE,循環(huán)取數(shù)并行完成第i行(列)的數(shù)據(jù)運算。在整個矩陣乘運算過程中,除當前行自相乘得到對角元素數(shù)據(jù)外,其余各行交叉提供源數(shù)據(jù),得到非對角元素結(jié)果。

以前2行(列)為例,在計算第1行(列)時,隨著列計數(shù)的遞增,上三角逆矩陣第1行元素先依次和下三角逆矩陣第1列元素相乘,再依次和第2列、第3列……元素相乘。

當行計數(shù)與列計數(shù)相等時,首先利用第1組矩陣乘運算PE得到對角線元素。

當列計數(shù)大于行計數(shù)時,上三角逆矩陣第1行元素利用第1組PE依次和下三角逆矩陣第2列、第3列……元素相乘,上三角逆矩陣第2行、第3行……利用第2組PE與下三角逆矩陣第1列相乘,2組PE共同作用完成結(jié)果矩陣第1行和第1列的逆矩陣結(jié)果輸出。

在計算第2行(列)時,同樣首先上三角逆矩陣第2行元素利用第1組PE和下三角逆矩陣第2列元素相乘,之后第3行、第4行……元素利用第2組PE與第2列元素相乘,2組PE共同作用完成結(jié)果矩陣第2行和第2列的結(jié)果輸出。

此后的每行(列)依照上述相同規(guī)律相乘直到完成整個逆矩陣的結(jié)果輸出。

求4階三角逆矩陣的矩陣乘結(jié)果數(shù)據(jù)搬運流程如圖3所示。圖3中沿著對角線對稱的上三角逆矩陣元素和下三角逆矩陣元素對稱存儲在Memory0的相同地址空間,結(jié)果元素回寫至Memory0中;未用箭頭標出的元素與已標出元素執(zhí)行相同的操作。

各PE內(nèi)部互連結(jié)構(gòu)如圖4所示。圖4中,M0in表示存入Memory0的數(shù)據(jù);D1in表示存入dist-mem1的數(shù)據(jù);D0in表示存入dist-mem0的數(shù)據(jù);MUX表示運算中的符號判斷;reg=1表示寄存器數(shù)據(jù)為1;Rout0表示求得的三角逆矩陣對角元素;Rout1、Rout2表示求得的三角逆矩陣非對角元素。

約化系數(shù)和預處理運算模式下,PE內(nèi)各運算器互連結(jié)構(gòu)示意圖如圖5所示。共2組PE,每組PE內(nèi)各8路PE1,并行完成約化系數(shù)矩陣的上、下三角元素計算,PE2則計算完成三角逆矩陣的對角元素數(shù)據(jù)。圖5中,Control-1表示約化系數(shù)運算和預處理控制端;MUL-Group表示乘法運算塊;PE1-Group表示約化系數(shù)運算塊;Mem表示各存儲塊。

三角逆矩陣運算模式下,PE內(nèi)各運算器互連結(jié)構(gòu)示意圖如圖6所示。共2組PE,每組PE內(nèi)各8路PE3,并行完成上、下三角逆矩陣計算。圖6中,Control-2表示上、下三角逆矩陣計算控制端;PE3-Group表示上、下三角逆矩陣元素運算塊。

圖3 對三角逆矩陣進行矩陣乘的數(shù)據(jù)搬運流程

圖4 各PE內(nèi)部互連結(jié)構(gòu)

圖5 2×8路約化系數(shù)運算和預處理互連結(jié)構(gòu)

圖62×8路三角逆矩陣運算互連結(jié)構(gòu)

3 實現(xiàn)與對比

本文利用Matlab選取矩陣元素在[-1,1]和[-100,100]2種取值范圍內(nèi),各階順序主子式不為0的2n矩陣作為測試源矩陣,基于存儲空間的選取大小,驗證完成24(16)階、25(32)階、26(64)階和27(128)階這4種規(guī)模的矩陣求逆。

各階數(shù)矩陣在不同源數(shù)據(jù)范圍內(nèi)的求逆硬件測試結(jié)果和Matlab仿真結(jié)果對比,誤差見表3所列。

由表3可知,本文采用原位替換算法進行矩陣求逆與Matlab仿真實驗結(jié)果相比,硬件電路求逆的結(jié)果精度能達到10-6,具有良好的運算性能。經(jīng)測試所設(shè)計的硬件架構(gòu)綜合頻率在120 MHz以上,能夠滿足高速運算的要求。

表3 矩陣求逆誤差

所設(shè)計的硬件架構(gòu)完成矩陣求逆運算的理論計算時間與實際計算時間對比見表4所列。

表4 完成求逆所用時間與理論時間

注：運算用時數(shù)據(jù)為周期數(shù)。

與文獻[10]相對比,本文設(shè)計用時更短,32階矩陣求逆能達到近10倍的加速比;與采用常規(guī)LU分解[11]和QR分解進行矩陣求逆的硬件設(shè)計相比,本文設(shè)計采用按位替換法,通過創(chuàng)建約化系數(shù)矩陣,分解了串行迭代的過程,提高了算法的可并行度。采用原位替換和LU分解進行矩陣求逆的硬件加速比見表5所列。

表5 矩陣求逆硬件運算加速比

原位替換法進行矩陣求逆的整個計算過程中,各個模塊都根據(jù)相應(yīng)的對稱軸進行并行運算,各個模塊結(jié)果均儲存在原存儲空間,不占用額外資源,有效地壓縮了資源消耗。由于各模塊分步有序進行,不會產(chǎn)生由于數(shù)據(jù)運算不完全帶來的拖尾效應(yīng),極大地保證了結(jié)果的正確性和數(shù)據(jù)的安全性。原位替換法實現(xiàn)了運算過程間分步有序,同一運算過程內(nèi)上、下三角并行,同一三角內(nèi)PE并行的“三維運算”。

4 結(jié) 論

本文分析了面向復雜信號處理和高密度計算的算法類型特點,根據(jù)按位替換法進行矩陣求逆的算法原理,設(shè)計實現(xiàn)了一款能夠完成2n階矩陣求逆運算的硬件架構(gòu),并且提出了一種針對算法的地址無沖突設(shè)計,在保證多維、高并行度計算的情況下,平衡了存儲接口、數(shù)據(jù)安全、資源消耗及運算速度的問題。