矩陣特征值分解的算法研究與FPGA 實現

黃家俊，王建新

（南京理工大學電子工程與光電技術學院，江蘇南京 210094）

陣列信號處理中的一類典型問題是如何從接收到的多路觀測數據中解析出被混疊著的原始信號[1]。常用的陣列信號分離算法有波達方向（Direction of Arrival，DOA）估計[2-4]、盲源分離[5]、快速獨立成分分析（Fast Independent Component Analysis，FastICA）[6-9]等。這些算法都需要涉及大量的矩陣運算，而矩陣特征值分解就是其中的關鍵步驟之一。

常規的特征值分解算法主要針對實對稱矩陣，如Jacobi 旋轉法、QR 迭代法[10]，其中Jacobi 法計算簡單、精度高，應用更為廣泛。傳感器陣列觀測到的信號經數字正交下變頻后變成I/Q 兩路的復信號，其協方差矩陣是一個復Hermitian 矩陣，無法直接應用傳統的算法。所以需要通過一定的運算先將復Hermitian 矩陣轉化為實對稱矩陣。

FPGA 硬件相較于DSP 硬件，運行速度快，并行運算能力強，更適合解決陣列信號處理中的矩陣運算。因此，該文研究陣列信號處理中典型的復Hermitian 矩陣特征值分解問題，在FPGA 中實現矩陣轉化和特征值分解的運算。

1 矩陣轉化

假設由M個天線組成的天線陣列觀測到的信號經正交下變頻后得到的復觀測數據為x(t)，如式（1）所示。

其中，N為x(t)的采樣長度。顯然矩陣Rx是一個M×M階的復Hermitian 矩陣。

復Hermitian 矩陣向實對稱矩陣的轉化有兩種處理方法。第一種轉化方法是將M×M階的復Hermitian 矩陣轉換為2M×2M階的實對稱矩陣[11]，首先用實對稱矩陣Rr和反對稱實矩陣Ri作為實部和虛部來表示矩陣Rx。假設矩陣Rx的特征值σ對應的特征向量實部為u、虛部為v，可以得到下式：

易證式（3）中的分塊矩陣是實對稱的。綜上，M×M階的復Hermitian 矩陣Rx的特征值分解等價于求解由實矩陣Rr、Ri組成的2M×2M階分塊矩陣的特征值分解。這種轉化方法是嚴格等價的，沒有誤差。但是這種方法擴大了待求解矩陣的階數，運算量大。

為了減少運算量，第二種轉化方法將M×M階復Hermitian 矩陣轉化為M×M階的實對稱矩陣。定義一個M×M階分塊矩陣U,M為偶數時，結構如下：

定義滿足如下等式的復Hermitian 矩陣P為中心對稱的復Hermitian 矩陣：

不難證明可以通過下式將矩陣Rx轉化為一個中心對稱的復Hermitian 矩陣R[2,4]：

定義一個由矩陣R變換而來的Hermitian矩陣K：

綜上，矩陣K是一個實對稱矩陣。定義UK為矩陣K的特征向量矩陣，∑K為矩陣K的特征值矩陣，那么矩陣R的特征值分解可以表示為：

根據式（9）可以得到矩陣R的特征向量矩陣為UHUK，矩陣R的特征值矩陣就是∑K。

如果M為奇數，此時，矩陣U結構類似，矩陣I、J都是階的，上述證明依舊成立。

綜上，可以通過一定的運算將原本的復Hermitian 矩陣Rx先變換為一個中心對稱的復Hermitian 矩陣R，再將矩陣R等價變換為一個實對稱矩陣K。在求解K的特征值分解后得到R的特征值分解，從而替代Rx的特征值分解。第二種轉化方法不擴大原矩陣階數，計算量相較于第一種方法是大幅度減小的。當然這種轉化并不是嚴格等價的，對于矩陣R和矩陣Rx之間的轉化誤差，可以證明在所有的復Hermitian中心對稱矩陣中，矩陣R是對矩陣Rx在歐式距離（Euclidean Distance）下的最佳估計[12]。

2 Jacobi旋轉法

Jacobi 法的基本原理為如果一個方陣是三角矩陣或對角陣，則該方陣主對角線上的元素就是該方陣的特征值。因此可以利用旋轉矩陣不斷地對實對稱矩陣K進行正交相似變換，通過變換將方陣K中所有非主對角線元素都消去，在精度足夠高的情況下，最后得到的近似的對角陣就是原矩陣K的特征值矩陣。矩陣K的特征向量矩陣也可以通過類似的運算來求得[13-14]。

首先基于單位矩陣I定義一個M×M階旋轉矩陣Qij(1 ≤i

易證矩陣Qij是一個正交矩陣，令K0=K，利用旋轉矩陣Qij對矩陣K0進行一次正交相似變換，得到矩陣K1：

變換后的矩陣K1仍然是一個實對稱矩陣，設矩陣K1中第p行第q列（1 ≤p

1）當p≠i或j且q≠i或j時，變換前后的矩陣元素相等，不發生變化。

2）當p=i,q=j時：

3）當p=i或j，q≠i或j時：

為了提高效率，每次正交相似變化都應該消去矩陣上三角部分中絕對值最大的元素。定義實對稱矩陣G的特征值矩陣CV、特征向量矩陣EV，利用Jacobi 旋轉法求解矩陣G的特征值分解，一般步驟如下：

1）給矩陣CV和矩陣EV賦初值（矩陣CV初值為矩陣G，矩陣EV初值為單位矩陣I）；

2）尋找矩陣CV上三角部分中絕對值最大的元素cvij，將其與允許誤差ej進行比較，小于ej進行步驟6），否則進行步驟3）；

3）根據元素cvij、cvii、cvjj的值計算旋轉角θ；

4）計算旋轉角θ的正余弦值，構造旋轉矩陣Qij；

5）計算變換后的矩陣CV、矩陣EV，然后進行步驟2）：

6）輸出矩陣CV和矩陣EV。

Jacobi 旋轉法的計算時間取決于兩點：一是ej的大小，ej越小表明計算出來的特征值矩陣越接近對角陣，算法誤差越小，需要的迭代次數就越多。二是待分解矩陣G的大小，G越大，待消去的元素數量越多，迭代次數也就越多，并且單次迭代時矩陣乘法的運算量也越大。

3 Jacobi旋轉法的FPGA實現

上一節的算法流程中，Jacobi 法的迭代次數不固定。在FPGA 中迭代次數不固定會導致時延的不確定，因此需要根據待分解矩陣的階數以及數據的量化精度來固定迭代次數。對于ej而言，在FPGA 中會受到量化精度和計算精度的影響，過小的ej是無法達到的。根據仿真，16 位量化精度時，4 階方陣需要15 次迭代；24 位量化精度時，6 階方陣需要25 次迭代。

利用FPGA 實現Jacobi 法求解特征值矩陣時，需要解決兩個問題：一是三角函數的計算，二是矩陣乘法的實現。

Jacobi 迭代法中三角函數的計算主要包括兩個部分：一是利用反正切arctan 函數計算旋轉角θ的值，二是計算旋轉角θ的正、余弦值。這兩步計算都可以采用ISE 軟件中自帶的Cordic IP 核來完成，計算精度很高。

Jacobi 法中每次迭代需要進行3 次矩陣乘法運算。雖然可以用式(12)～(17)來直接計算每個元素值，但每次消去的元素位置不固定，導致每個元素的表達式也不固定，不如直接通過矩陣乘法來計算旋轉后的矩陣。在FPGA 中，可以采用脈動陣列（Systolic Array）的結構來構造矩陣乘法器[16-17]。脈動陣列由多個相同的處理單元（Processing Elements，PE）按照特定的規則互聯而成。

首先對矩陣乘法進行拆分，設m×n階矩陣C是由m×l階矩陣A和l×n階矩陣B相乘得到的：

其中，i=1，2，…，m；j=1，2，…，l。

不難發現每一個PE 需要一個乘法器和一個累加器來完成對應元素的計算。為了使整個運算過程流水線化，輸入的數據需要在PE 間“流動”起來，所以每個PE 都要對兩路輸入數據做一個單位的延時輸出來提供給下一級PE 進行運算[18]。綜上，矩陣乘法器中PE 的框圖如圖1 所示。

圖1 PE結構圖

為了實現m×l階矩陣A與l×n階矩陣B相乘，需要m×n個這樣的PE 來組成脈動陣列。每一個PE 的輸出對應矩陣C中的一個元素，用脈動陣列實現的矩陣乘法器結構如圖2 所示。

圖2 脈動陣列矩陣乘法器結構圖

上圖輸入數據中的0 表示對輸入數據做一個時鐘的延遲以匹配上一級PE 的延時輸出。如果預先對輸入數據流進行控制可以提高PE 的利用率，減少PE 的數量，從而節省乘法器和加法器資源。

在矩陣乘法運算后，要對結果進行截位，截位策略是保證未參與運算的元素值不發生變化。根據仿真結果，特征值矩陣需要提前進行位擴展，以保證矩陣對角線上的特征值不溢出[19]。

4 仿真結果

以16 位定點數量化后的4 階復Hermitian 矩陣Rx為例，其實部為：

使用第二種矩陣轉化方法后用Jacobi 旋轉法求特征值分解，其中Jacobi 法的迭代次數為15 次，Matlab 軟件中計算出的特征值矩陣CV如下：

在ISE 軟件中編寫好程序后，用Modelsim 仿真得到的特征值矩陣CVM和特征向量矩陣EVM，如圖3所示。模塊端口圖如圖4 所示。

圖3 Modelsim仿真結果圖

圖4 模塊端口圖

其中，clk 為輸入時鐘端口，en 為模塊輸入使能端口，rst_n 為模塊異步清零端口，Snr_in1～Snr_in4 是矩陣Rx實部輸入端口，Sni_in1～Sni_in4 是矩陣Rx虛部輸入端口，每一個端口輸入待求解矩陣的一行數據；valid 端口為輸出有效端口，計算結果輸出時同步置高，CV_out1～CV_out4 是矩陣CVM的輸出端口（特征值矩陣在迭代開始時就已經做了防溢出位擴展），EVr_out1～EVr_out4 是矩陣EVM實部輸出端口，EVi_out1～EVi_out4 是矩陣EVM虛部輸出端口。整個模塊從數據輸入到結果輸出的延時為1 074 個時鐘周期。由于Cordic IP 核的計算存在細小的誤差，導致仿真結果與Matlab 定點仿真結果有一定的誤差，不過對整體精度影響非常小。‖ ‖?2表示對矩陣取二范數，Modelsim 仿真結果與Matlab 仿真結果的歸一化誤差為：

整個模塊的資源消耗情況如表1 所示。

表1 FPGA模塊資源消耗

5 結束語

該文對矩陣特征值分解算法進行了研究并給出了求解復Hermitian 矩陣的特征值分解的一般方案。該方案利用FPGA 完成復Hermitian 矩陣到實對稱矩陣的轉化后，再運用Jacobi 旋轉法求解轉化后的實對稱矩陣特征值分解。整個方案具有結構簡單、計算精度高、時延小、資源消耗少的特點。