基于FPGA 的復數(shù)矩陣求逆設計

周金強，鳳繼鋒

（四川大學電子信息學院，成都610065）

0 引言

矩陣運算[1]在現(xiàn)代科學工程領(lǐng)域的基本問題，如何更快，效率更高地進行矩陣運算是科學計算研究的重點。在信號處理和圖像識別領(lǐng)域中，矩陣乘法、求逆是最基本的數(shù)學工具，大量算法都涉及到這些基本運算，因此如何更快、更節(jié)省運算資源來進行矩陣的求逆運算是現(xiàn)代工程計算領(lǐng)域一直在研究的重點。但是，當今信號處理算法大部分是在PC 端的軟件上進行實施的，在一些講究實時性的場景，如無人機電子偵查，電磁頻譜感知[2]等領(lǐng)域，這些場景對實時性要求高，并且受到環(huán)境影響，無法時刻裝配大型計算機來進行信號算法運行，如調(diào)制寬帶轉(zhuǎn)換器（Modulated Wideband Converter，MWC）[3]，其混頻、濾波、采樣都在硬件上執(zhí)行，但恢復算法則是在MATLAB 端進行，無法實時恢復出輸入信號的頻譜，這對該系統(tǒng)的實際應用形成了限制，使其無法滿足特定場景下實時恢復出信號的頻譜圖，無法滿足實時性要求，而FPGA 平臺則可以滿足這些特點，其功耗，運行速度都有優(yōu)勢，并且能夠適應更為復雜的環(huán)境因素，通過將其使用的OMP 算法進行FPGA 平臺移植可使該系統(tǒng)滿足實時性場景，而OMP算法的硬件化的關(guān)鍵問題在于最小二乘部分，該部分在數(shù)學上表示為矩陣求逆計算，若用傳統(tǒng)定義法或者SVD 分解法之間進行求偽逆需要消耗大量資源，使得硬件運算速度下降，因此如何選擇合適的分解求逆法，并在此基礎(chǔ)上設計高速，邏輯資源使用合理的求逆模塊是本文的重點。

1 矩陣分解求逆數(shù)學表示

1.1 改進型Cholesky分解

OMP 算法中所遇到的最小二乘問題矩陣為非方陣，可以通過Moore-Penrose 求解方法將偽逆問題轉(zhuǎn)化為求解共軛對稱方陣的逆，而對于該類型的矩陣，改進型Cholesky 方法[4]有著較好的分解效果，相比較傳統(tǒng)Cholesky 分解法，改進型Cholesky 分解規(guī)避了開方計算，而開方計算在硬件執(zhí)行上需要消耗大量資源，復雜運算只涉及到倒數(shù)計算，并且由于矩陣對角線都為正實數(shù)，因此所消耗的除法計算與矩陣維度一致，相比較其他分解法如QR 分解[5]，SVD 分解，復雜計算較少，因此本文選擇改進型Cholesky 分解法，下面介紹其數(shù)學原理：

假設可逆方陣為A，其元素以aij來表示，將矩陣進行分解為A=LDLH的形式：

其中L 為對角線元素全為1 的下三角矩陣，D 為對角矩陣，LH為L 矩陣的共軛轉(zhuǎn)置，通過矩陣計算，可得計算結(jié)果為：

將ljj=1代入，之后進行化解，移項，可得：

為方便計算，引入輔助矩陣G 和C，代入公式可得：

此為改進型Cholesky 分解的數(shù)學表達。

1.2 下三角矩陣求逆

根據(jù)分解公式，對A=LDLH進行兩邊求逆，可得A-1=(L-1)HD-1L-1，其中D 的逆即為輔助矩陣C，因此求逆關(guān)鍵在于下三角矩陣L 的求逆[6]，考慮矩陣B 為L 的逆，根據(jù)單位矩陣定義，可得：

令P=LB，根據(jù)乘法運算：

又因pi,j=0(i ≠j),pii=1，最后將對角線1 代入，解得遞推關(guān)系：

由此可得求逆的遞推公式，根據(jù)公式，該計算方法規(guī)避了除法計算，在硬件上可節(jié)省資源，使運算速度得以提升。

2 模塊設計

根據(jù)先前的計算公式，本文設計了2 個基本運算模塊[7]，一是改進型Cholesky 分解模塊，二是下三角矩陣求逆遞推模塊，最后的三角矩陣相乘則通過矩陣之間的數(shù)值關(guān)聯(lián)信息，考慮到最終矩陣也為共軛對稱陣，計算時只需計算上三角及對角線即可，本文所使用的數(shù)據(jù)類型為（1,7,19）的浮點數(shù)，因此基本運算模塊為ISE 14.7 的浮點數(shù)計算IP 核，并且由于是復數(shù)，考慮到資源消耗，復數(shù)乘法計算使用兩個浮點數(shù)乘法器，一個加法器一個減法器，消耗兩個計算時鐘來進行乘法計算，下面介紹子模塊的硬件設計圖。

圖1 為分解模塊框圖，該模塊的主要運算單元為3個復數(shù)浮點數(shù)乘法器，一組減法器（分別計算實部虛部）和一個浮點數(shù)除法器（用于計算倒數(shù)），三組RAM用于存儲迭代運算過程中的中間變量，RAM_L 存儲計算所得的L 矩陣的下三角元素，而RAM_C 則是存儲對角矩陣D 的逆矩陣元素。三個Sel 為數(shù)據(jù)選擇器，分別用于選擇乘法器、減法器、除法器的輸入數(shù)據(jù)以及將輸出數(shù)據(jù)寫入到對應RAM 的地址控制操作，根據(jù)所得的迭代順序，三個乘法器用于并行計算同列標的中間變量值以及L 矩陣元素，減法器則是將計算所得的中間變量進行累減操作，求得下一列的G 矩陣元素，新一列的G 矩陣元素根據(jù)行標大小進行有序計算，當求得第一個元素即gj,j后，倒數(shù)模塊啟動，除法器一端輸入固定參數(shù)(1f80000)16，即十進制1，作為被除數(shù)，另一端則輸入gj,j進行倒數(shù)計算。為了方便地址控制，便于后續(xù)更高維矩陣的分解計算，將每次計算的時鐘周期設置為35 個，保證浮點數(shù)乘除減法模塊的運算正確性，以4 維復數(shù)矩陣為例，列數(shù)為1 的中間變量元素有6 個，3 個復數(shù)乘法器用2 個計算時鐘即可求得，行標小的元素先進行計算，之后啟動累減模塊，根據(jù)減法次數(shù)消耗相應的計算時鐘，對角倒數(shù)操作則是在減法器求第二個gi,j時并行進行，新的對角倒數(shù)求得即可進行新一列的L 矩陣元素計算，此過程只涉及到乘法，此時啟動乘法器，根據(jù)減法器求得的元素順序進行下三角及中間變量計算，以此類推，由于只使用一組減法器，時鐘消耗主要在于減法計算次數(shù)。

圖1 分解模塊圖

分解后L 矩陣輸入到下三角求逆模塊，根據(jù)遞推公式，該模塊只涉及到乘加計算，因此結(jié)構(gòu)較為簡單，只需單個復數(shù)乘法器和一組加法器同時計算實部虛部即可，圖2 為模塊圖。

圖2 下三角求逆模塊圖

最后矩陣相乘為基本的點積計算，設計時考慮到數(shù)據(jù)的關(guān)聯(lián)性[8]，只需計算上三角部分及對角線部分，下三角部分元素只需對對應的上三角部分虛部數(shù)據(jù)進行取反操作即可，對于本文所使用的浮點數(shù)來說，將二進制數(shù)據(jù)前3 位于100 進行異或計算即可得出相反數(shù)，硬件上極易實現(xiàn)：

3 仿真結(jié)果分析

本文所用的編譯平臺為Xilinx 的ISE 14.7，使用的仿真軟件為ModelSim 10.1c，分解模塊的仿真時序圖如圖3。

下三角矩陣求逆的時序圖如圖4。

圖3 分解時序圖

圖4 求逆時序圖

最終相乘求逆結(jié)果的時序圖如圖5。

圖5 求逆時序圖

根據(jù)時序圖的計算結(jié)果，可見其都是按一定的遞推順序求解所得，滿足先前所推導的計算公式。將最終計算結(jié)果進行數(shù)值轉(zhuǎn)換，與MATLAB 端計算結(jié)果相比，進行相對誤差分析，所設計的模塊計算精度在保留5 位小數(shù)的情況下，與MATLAB 端計算結(jié)果一致，綜合相對誤差為百分之6.82×10-5，表1 為計算結(jié)果對比，考慮輸入矩陣維度為4 維方陣，取第一列數(shù)據(jù)，W 為分解過后求逆過程中最后矩陣相乘的結(jié)果，即最終的求逆結(jié)果。

表1 誤差分析（取第一列數(shù)據(jù)）

分解模塊的邏輯資源消耗情況如如圖6。

圖6 分解模塊邏輯資源

根據(jù)時序分析報告，分解模塊的運行頻率最高可以達到326.535MHz，而下三角矩陣求逆模塊的運行頻率最高為454.277MHz，圖7-圖8 為分析時序報告圖。

圖8 求逆時序報告

根據(jù)模塊設計的計算時鐘數(shù)，可計算出最終的硬件模塊矩陣求逆運行時間為4.322 μs，本文所使用的PC 平臺為CPU Intel i5-7300H，內(nèi)存為16GB，使用MATLAB 對相同四維復數(shù)矩陣進行求逆，即使用inv（）函數(shù)，通過tic，toc 測得最短運行時間為37 μs，經(jīng)計算，硬件模塊所消耗的計算時間比MATLAB 端快了8.5 倍，節(jié)省了88.32%的計算時間，同時計算精度高，滿足一般的工程需求，該模塊可以運用到涉及最小二乘法的信號處理硬件算法中，提高效率。

4 結(jié)語

本文根據(jù)最小二乘類問題的硬件實施上的難點，通過公式簡化將偽逆計算難點轉(zhuǎn)化為求解特殊方陣的逆，并選擇針對該類共軛對稱矩陣分解效果較好，硬件實施上較為合理的改進型Cholesky 法作為分解方法，通過分解完的求逆計算，將下三角矩陣的求逆問題通過數(shù)學遞推推導轉(zhuǎn)化為只涉及乘加計算的求逆計算，最后的矩陣相乘由于數(shù)據(jù)關(guān)聯(lián)性及共軛對稱的關(guān)系，節(jié)省了多余的計算步驟，節(jié)約了邏輯資源消耗，并根據(jù)所推導的數(shù)學關(guān)系設計了對應的硬件模塊，給出了其模塊原理圖，并通過ModelSim 仿真得出時序圖及計算結(jié)果，根據(jù)時序分析圖計算出模塊的理論運行時間，兩者與MATLAB 計算結(jié)果及時間消耗進行比對，可見所設計的硬件模塊滿足一般工程精度要求，且計算時間也遠低于MATLAB 端，能夠很好的運用在實時性強的場景中。