一種基于Cholesky分解的快速矩陣求逆方法設(shè)計(jì)

魏嬋娟，張春水，劉 健

（中國空間技術(shù)研究院航天恒星科技有限公司 北京 100086）

抗干擾接收機(jī)由于增加了抗干擾處理的部分將會對導(dǎo)航信號的實(shí)時解算造成一定的影響，求取抗干擾濾波權(quán)值用時越長，對解算的延遲越大，造成的定位誤差也將越大。在權(quán)值更新以及抗干擾濾波的整個處理過程中，矩陣求逆的用時占用了90%以上的時間，矩陣求逆用時過長是導(dǎo)致接收機(jī)權(quán)值更新過慢的最主要因素，因此若采用直接矩陣求逆對抗干擾算法進(jìn)行實(shí)現(xiàn)，必須保證矩陣求逆的用時極短。本文針對大維數(shù)赫米特矩陣求逆實(shí)現(xiàn)方法進(jìn)行FPGA設(shè)計(jì)及實(shí)現(xiàn)。

1 Cholesky分解求逆方法及其流水設(shè)計(jì)

A為正定的赫米特矩陣，采用Cholesky分解矩陣求逆方法的實(shí)現(xiàn)過程為[1]:

1）由Cholesky分解得到下三角矩陣L及其共軛轉(zhuǎn)置矩陣LH，A =L×LH；

2） 對下三角矩陣L進(jìn)行矩陣求逆計(jì)算，得到其逆矩陣IL=L-1；

3） 由IL矩 陣 求 得A的 逆 矩 陣IA，IA =（LH）-1×L-1=LH×IL。

在FPGA實(shí)現(xiàn)赫米特矩陣求逆的關(guān)鍵在于實(shí)現(xiàn)方法的流水設(shè)計(jì)，合理高效的流水方法可以最大程度地節(jié)約資源提高實(shí)現(xiàn)速度，以下分3部分實(shí)現(xiàn)赫米特矩陣Cholesky分解求逆的實(shí)現(xiàn)方法設(shè)計(jì)[2-4]。

1.1 矩陣的Cholesky分解

對A矩陣進(jìn)行Cholesky分解，得到下三角矩陣L，A=LLH。采用分塊的方法計(jì)算L矩陣[5]。

將矩陣L及A進(jìn)行分塊表示如下：

由A=LLH推導(dǎo)可得：

根據(jù)如上公式表示，L的第一列的值可由A的第一列的值獲得，其計(jì)算方法為：

在計(jì)算完L的第一列的之后，之后需要計(jì)算L1，根據(jù)式2可得:

在將A矩陣的部分更新為A1-llH后，可采用相同的方法計(jì)算出L1矩陣的第一列的值，即L矩陣的第2列的值，依次類推可計(jì)算出L的所有值。

對以上的Cholesky矩陣分解方法進(jìn)行FPGA實(shí)現(xiàn)的流水設(shè)計(jì)，以4×4矩陣為例，可分為四步分析各個數(shù)據(jù)之間的依賴關(guān)系并計(jì)算L各列的值。

圖1 矩陣Cholisky分解步驟圖Fig.1 Four steps of Cholisky decomposition

步驟1：計(jì)算L矩陣的第1列；之后更新A矩陣（圖中，左側(cè)的矩陣為A矩陣，右側(cè)為L矩陣，箭頭所指向的方向?yàn)橛?jì)算當(dāng)前值所依賴的數(shù)據(jù)，例如將頭由a指向b，表示計(jì)算a時需要用到b，需在b計(jì)算完之后才可進(jìn)行a的計(jì)算）。

步驟2：計(jì)算L矩陣的第2列；之后更新A矩陣。

步驟3：計(jì)算L矩陣的第3列；之后更新A矩陣。

步驟4：計(jì)算L矩陣的第4列；A矩陣無需再更新。

對于4×4矩陣，在經(jīng)過4步之后可獲得L矩陣。

1.2 下三角矩陣求逆

對下三角矩陣L進(jìn)行求逆計(jì)算得到的逆矩陣IL[6]，根據(jù)L*IL=E，分別求得IL各列的值。

L矩陣求逆數(shù)據(jù)依賴關(guān)系及流水設(shè)計(jì)如圖2所示。

步驟1：計(jì)算IL對角線上的值（圖中，左側(cè)的矩陣為L矩陣，右側(cè)的矩陣為IL矩陣）。步驟2：計(jì)算IL第二對角線上的值。步驟3：計(jì)算IL第三對角線上的值。步驟4：計(jì)算IL第四對角線上的值。

通過對對角線上的值進(jìn)行同步計(jì)算，在經(jīng)過4個步驟之后可計(jì)算出IL矩陣的所有值。

1.3 矩陣相乘

上三角矩陣LH與下三角矩陣L相乘得到A的逆矩陣IA，IA=ILHIL。

矩陣相乘流水設(shè)計(jì)數(shù)據(jù)間的依賴關(guān)系（圖中左側(cè)矩陣為L矩陣，右側(cè)矩陣為IA矩陣）。

矩陣相乘較為簡單，其各個數(shù)據(jù)間的計(jì)算無需等待，僅僅依賴于已知的下三角矩陣IL。

圖2 下三角矩陣求逆數(shù)據(jù)依賴關(guān)系圖Fig.2 Data dependency relationship of lower triangular matrix inversion

圖3 矩陣相乘的數(shù)據(jù)依賴關(guān)系圖Fig.3 Data dependency relationship of matrix vector multiplication

矩陣相乘較為簡單，其各個數(shù)據(jù)間的計(jì)算無需等待，僅僅依賴于已知的下三角矩陣IL。

2 FPGA實(shí)現(xiàn)

2.1 狀態(tài)控制及功能描述

與1節(jié)Cholesky分解求逆流水方法對應(yīng)，F(xiàn)PGA硬件實(shí)現(xiàn)電路包括3個主要模塊：Cholesky分解、下三角矩陣求逆和矩陣相乘。其在FPGA中實(shí)現(xiàn)的狀態(tài)控制及存儲如下圖所示：

圖4 FPGA實(shí)現(xiàn)Cholesky矩陣分解求逆狀態(tài)控制圖Fig.4 FPGA implementation of Cholesky decomposition and inversion

功能描述：

1）Start_inv為矩陣分解求逆開始標(biāo)志，save狀態(tài)下，將生成的協(xié)方差矩陣存入RAM_A；

2）Chol狀態(tài)下完成協(xié)方差矩陣的Cholesky分解，得到的下三角矩陣L存入RAM_B；

3）Inv狀態(tài)下完成下三角矩陣的求逆計(jì)算，將得到的逆矩陣IL存回RAM_A；

4）Inv_mult狀態(tài)下完成ILH矩陣與IL矩陣的相乘，即上三角陣與下三角陣的相乘，其結(jié)果存入RAM_B。

5）在Inv_mult狀態(tài)結(jié)束后進(jìn)入Output狀態(tài)將RAM_B中的逆矩陣輸出。

6 ）在Ouptput狀態(tài)結(jié)束后重新返回start狀態(tài)，開始對下一組協(xié)方差矩陣的分解求逆。

2.2 用時測試

跟據(jù)以上的流水設(shè)計(jì)進(jìn)行FPGA代碼的編寫實(shí)現(xiàn)，矩陣求逆的整體用時如下表所示：

表1 FPGA實(shí)現(xiàn)Cholesky矩陣求逆用時Tab.1 FPGA implementation time table

由上表看出，28×28維矩陣的求逆用時小于0.4 ms，35×35維矩陣的求逆用時僅需不到0.8 ms。

3 結(jié) 論

文中針對大維數(shù)赫米特矩陣的求逆的問題展開，介紹了適用于正定赫米特矩陣的Cholesky矩陣分解求逆方法，并對該方法進(jìn)行了高效的流水設(shè)計(jì)，分3步實(shí)現(xiàn)矩陣的FPGA分解求逆，經(jīng)過試驗(yàn)測試，用時極短。

[1]張賢達(dá).矩陣分析與應(yīng)用[M].北京:清華大學(xué)出版社, 2004.

[2]郭磊.矩陣運(yùn)算的硬件加速技術(shù)研究[D].長沙:國防科學(xué)技術(shù)大學(xué), 2010.

[3]Thompson, J. Benkrid, K. Xuezheng Chu. Rapid Prototyping of an Improved Cholesky Decomposition Based MIMO Detector on FPGAs [J].Adaptive Hardware and Systems, 2009(4):369-375.

[4]Mach D M,Koshak W J.General matrix inversion for the calibration of electric field sensor arrays on aircraft platforms [J].Marshall Space Flight Center, 2006(6):1576-1587.

[5]魯翠仙.分塊矩陣在矩陣求逆中的應(yīng)用[D].昆明:云南大學(xué),2009.

[6]彭玲,彭大芹. 一種下三角復(fù)矩陣求逆方法的IP設(shè)計(jì)與實(shí)現(xiàn)[J].電子測試, 2011(10):9-12.

PENG Ling, PENG Da-qin.IP design and implementarion of lower trlangular matrx invertion[J].Electronic Test,2011(10):9-12.