一種基于Cholesky分解的快速矩陣求逆方法設計

魏嬋娟，張春水，劉 健

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

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

1 Cholesky分解求逆方法及其流水設計

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

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

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

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

在FPGA實現赫米特矩陣求逆的關鍵在于實現方法的流水設計，合理高效的流水方法可以最大程度地節約資源提高實現速度，以下分3部分實現赫米特矩陣Cholesky分解求逆的實現方法設計[2-4]。

1.1 矩陣的Cholesky分解

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

將矩陣L及A進行分塊表示如下：

由A=LLH推導可得：

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

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

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

對以上的Cholesky矩陣分解方法進行FPGA實現的流水設計，以4×4矩陣為例，可分為四步分析各個數據之間的依賴關系并計算L各列的值。

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

步驟1：計算L矩陣的第1列；之后更新A矩陣（圖中，左側的矩陣為A矩陣，右側為L矩陣，箭頭所指向的方向為計算當前值所依賴的數據，例如將頭由a指向b，表示計算a時需要用到b，需在b計算完之后才可進行a的計算）。

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

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

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

對于4×4矩陣，在經過4步之后可獲得L矩陣。

1.2 下三角矩陣求逆

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

L矩陣求逆數據依賴關系及流水設計如圖2所示。

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

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

1.3 矩陣相乘

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

矩陣相乘流水設計數據間的依賴關系（圖中左側矩陣為L矩陣，右側矩陣為IA矩陣）。

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

圖2 下三角矩陣求逆數據依賴關系圖Fig.2 Data dependency relationship of lower triangular matrix inversion

圖3 矩陣相乘的數據依賴關系圖Fig.3 Data dependency relationship of matrix vector multiplication

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

2 FPGA實現

2.1 狀態控制及功能描述

與1節Cholesky分解求逆流水方法對應，FPGA硬件實現電路包括3個主要模塊：Cholesky分解、下三角矩陣求逆和矩陣相乘。其在FPGA中實現的狀態控制及存儲如下圖所示：

圖4 FPGA實現Cholesky矩陣分解求逆狀態控制圖Fig.4 FPGA implementation of Cholesky decomposition and inversion

功能描述：

1）Start_inv為矩陣分解求逆開始標志，save狀態下，將生成的協方差矩陣存入RAM_A；

2）Chol狀態下完成協方差矩陣的Cholesky分解，得到的下三角矩陣L存入RAM_B；

3）Inv狀態下完成下三角矩陣的求逆計算，將得到的逆矩陣IL存回RAM_A；

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

5）在Inv_mult狀態結束后進入Output狀態將RAM_B中的逆矩陣輸出。

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

2.2 用時測試

跟據以上的流水設計進行FPGA代碼的編寫實現，矩陣求逆的整體用時如下表所示：

表1 FPGA實現Cholesky矩陣求逆用時Tab.1 FPGA implementation time table

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

3 結 論

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

[1]張賢達.矩陣分析與應用[M].北京:清華大學出版社, 2004.

[2]郭磊.矩陣運算的硬件加速技術研究[D].長沙:國防科學技術大學, 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]魯翠仙.分塊矩陣在矩陣求逆中的應用[D].昆明:云南大學,2009.

[6]彭玲,彭大芹. 一種下三角復矩陣求逆方法的IP設計與實現[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.