浮點厄米特矩陣特征值分解的FPGA實現

王傳根 呂樂群 林都督

摘要：為了提高厄米特矩陣特征值分解計算的速度和精度，基于QR迭代的方法，在FPGA上實現了一種并行可重構的特征值分解方案。設計在Intel公司型號為10AX066H1F341的FPGA上進行驗證，通過調用硬浮點DSP計算資源，實現了8×8階的浮點厄米特矩陣的特征值分解。驗證結果顯示：該方案耗時11425個時鐘周期，計算誤差小于3.74e-7，資源消耗小于總資源的32%，最高運行時鐘頻率可達311MHz，可支撐MUSIC等復雜算法的高速硬件工程實現。

關鍵詞：厄米特矩陣;特征值分解;可重構;FPGA

中圖分類號：TN911.7 文獻標識碼：A 文章編號：1007-9416（2020）02-0133-03

在陣列信號處理領域中，基于特征分解的超分辨率空間譜估計算法以其高效高精度的特點，一直以來都是研究的焦點，這其中最為經典的就是多重信號分類（MUSIC）算法[1-2]。MUSIC算法具有測向精度高、算法穩(wěn)定、優(yōu)異的多信號測向能力等眾多優(yōu)點，然而其運算卻尤為復雜，不利于硬件實現。這其中以矩陣特征值分解（EVD）的運算最為復雜，大約占到了MUSIC計算60%的計算量，所以研究EVD在FPGA上的加速運算有著十分重要的意義。當前研究EVD在FPGA上的實現多是針對實對稱矩陣，并且多是采用定點的方式實現[3-5]。然而在實際工程應用中，經常需要計算復數自共軛矩陣（即厄米特矩陣）的特征值分解，并且隨著矩陣階數增大，定點數據在計算過程中的截位誤差累積，造成最終結果精度大幅下降，難以滿足應用需求。

常用的EVD分解的方法主要有基于Jacobi旋轉的方法和基于QR迭代的方法，基于Jacobi的方法精度高、收斂速度慢，且不能直接計算復數矩陣的特征值分解。相比Jacobi方法，基于QR迭代的方法精度稍差，但是其收斂速度更快，并行度高，可以直接對復數矩陣進行計算，對于加速計算的效果更為明顯，特別適合于FPGA并行實現。

本文以8×8階的厄米特矩陣特征值分解為目標，采用并行QR迭代的方案實現FPGA設計。通過采取浮點運算機制和算法結構優(yōu)化，在保證計算精度的同時，加快計算速度，降低資源消耗。

1 算法原理

1.1 基于QR分解的EVD計算

對于非奇異方陣A，進行QR分解，可得A=QR。其中Q為正交矩陣（即QTQ=I），R為上三角矩陣，則：

（1）

也即是說，RQ與A是正交相似的，它們具有相同的特征值。因此可以采用如下的QR迭代格式。

（2）

取A1=A，經過式（2）迭代，可以證明Ak收斂到矩陣A的特征值，Uk=Q1Q2…Qk收斂至矩陣A的特征向量。

1.2 基于QR分解的EVD計算

根據正交理論的思想，矩陣QR分解可以采用Givens旋轉變換實現，對矩陣Am×m=（a1，a2，…，am）T，Givens旋轉矩陣定義為：

（3）

其中，，，，且Gi，j為酉矩陣。執(zhí)行Givens旋轉變換，可使矩陣按照一定角度旋轉，使矩陣i行j列元素置換為0。由此，可以重復執(zhí)行這個變換過程，生成一系列Givens矩陣Gm，m-1…Gm，2…G3，2Gm，1…G2，1，將對角線以下的元素全部置換為0，得到上三角矩陣R，于是有：

（4）

按照（4）式所示的方法，可以迭代完成矩陣A的QR分解。然而每次使用Givens矩陣置換元素后，新的Givens矩陣必須基于Givens變換后的矩陣生成。由于前后數據的相關性，造成FPGA實現時只能順序執(zhí)行，無法發(fā)揮FPGA并行執(zhí)行的優(yōu)勢，難以提升運算速度。由此，對Givens旋轉矩陣進行改進，將每列Givens矩陣的計算展開，每次取矩陣A第j列對角線以下元素，令

（5）

（6）

將（5）式和（6）式代入（3）式計算Givens矩陣，一次得到j列的Givens矩陣，計算可將矩陣A的j列對角線以下的元素置換為0。因此可以每次并行計算，一次迭代置換1列元素，經過m-1次迭代即可完成QR分解。這樣的計算步驟非常利于FPGA并行實現，從而極大地提高運算速度。

2 FPGA實現方案

為了提升運算效率，設計選用搭載硬浮點DSP計算單元的IntelArria10系列FPGA實現。該DSP單元采用芯片硬件定制浮點設計，遵循IEEE754標準，單個DSP即可實現一個浮點加法和一個浮點乘法，無需額外的邏輯資源。此外通過調整DSP內部流水寄存器級數，可以滿足不同的latency和Fmax的需求，無需在timing上做過多努力，就能工作在極高的工作頻率。在浮點計算密集的情況下，Arria10比起業(yè)界同等級的FPGA在資源、效率、運算速度方面均有明顯優(yōu)勢。

在FPGA上實現EVD分解的結構如圖1所示，設計包含4個模塊：矩陣數據緩存（BUF），Givens矩陣計算單元（GCU）、可重構計算單元（RCU）和計算控制單元（CCU）。

BUF模塊實現矩陣A_in輸入緩存、迭代結果緩存和EVD輸出結果緩存，為加快隨機訪問速度，BUF模塊緩存均使用寄存器實現。CCU模塊為整個EVD運算的控制樞紐，根據迭代次數iter_num、啟動信號start以及各個模塊的狀態(tài)，控制迭代運算過程，完成EVD的計算，CCU模塊主體為控制狀態(tài)機，其狀態(tài)轉換如圖2所示。按照相關性可將EVD計算分為4個計算步驟：

（a）計算Givens旋轉矩陣G;

（b）計算旋轉矩陣連乘G_prod;

（c）計算G=G_prod*G和A=G_prod*A;

（d）計算Ak+1=Rk*Qk和Uk+1=Uk*Qk。

這4個步驟輸入輸出上存在關聯，無法拆解開來并行實現，但工作時間上是分開的，且計算相似性極高。故在保證步驟內高并行度的情況下，采用分時重構計算資源的設計方法。

GCU和RCU相互協作，根據CCU設置的工作模式，重構完成不同的計算功能。GCU實現計算過程中開平方和除法等個性的計算，同時設置RCU重構參數，調度RCU協處理向量二范數、復矩陣乘法等計算資源密集的共性操作。如此通過循環(huán)調度計算資源，即可完成EVD計算的各個步驟。RCU在不同工作模式下，完成如下計算：

在步驟（a）時，取矩陣的列向量按（5）式和（6）式計算各個旋轉矩陣的c和s，此時需要求取列向量的二范數，重構RCU并行完成多個的計算;

在步驟（b）時，需要計算一系列Givens矩陣的連乘，對于j列生成的Givens矩陣Gm，j乘以矩陣B，根據（3）式矩陣的特點，可以簡化計算，僅需計算乘積的第m行和j行，其余的行等于B，此時RCU完成（7）式的運算過程。

（7）

在步驟（c）和（d）時，需要計算兩個8×8矩陣的乘法，此時RCU完成矩陣乘法的計算，即計算。

經過分析剝離，RCU中只保留復數乘和復數加兩種操作，并以此為最小重構顆粒，設計RCU模塊如圖3所示。RCU將重構參數解析為數據選通MUX的選通信號，去選通相應的復數乘法器和加法器，實現以上RCU要求的功能，其中“取反”、“求共軛”操作由格式轉換子模塊實現。經分析，計算資源消耗最大的步驟為矩陣乘法計算，并行計算一個矩陣乘積元素需要8個乘法器和7個加法器。設計例化64個乘法器和56個加法器，可以一次并行計算8個（一行）的矩陣乘積結果，僅需8+N（N為計算延遲）拍就能完成一次矩陣乘法。其余的運算操作通過重構RCU，選通部分乘法器和加法器即可實現。通過采用RCU重構的方法，在步驟（a）～（d）中資源重構，可以節(jié)約60%以上資源。

3 設計實現及結果

本設計選用的器件為Intel公司的10AX066H1F341FPGA，使用Verilog語言進行電路描述。在Quartus上進行編譯綜合，并使用Modelsim進行仿真。為驗證設計正確性，在matlab上生成好激勵，同時將激勵用于Modelsim仿真和matlab上eig函數計算，然后比對仿真和計算結果。

由于矩陣分解特征向量不唯一，FPGA實現方法和matlab中eig函數實現方法不一致，導致特征向量不同，不能直接數值比較。本文主要進行特征值數值比較，同時根據矩陣特征值分解的特點，驗證特征向量是否為酉矩陣、特征值是否為對角線實數矩陣、VDVH是否等于矩陣A。為了評估設計的計算精度，按（8）式對特征值誤差進行計算。其中rDi，j為仿真得到的特征值，gDi，j為eig計算得到的特征值。

（8）

通過matlab隨機產生21組數據，仿真設置迭代次數為12次，對數據進行了仿真。結果顯示，仿真得到的特征向量均為酉矩陣，特征值為對角線實數矩陣，符合矩陣特征值分解的要求。特征值誤差如圖4所示，計算誤差均小于3.74e-7。仿真顯示，一次計算需11425個時鐘周期，耗時約57us（工作頻率200MHz時）。在Quartus上綜合報告如表1所示。

4 結語

本文基于Givens變換和QR分解的方法，采用并行可重構的設計策略，在FPGA上實現了厄米特矩陣特征值分解的計算。設計結果表明：在資源節(jié)省60%的情況下，仍能保證極高的精度和運算速度。該設計可以滿足MUSIC算法實時信號處理的需求，也可以在人工智能、圖像處理、盲源分離、頻譜分析等有EVD需求的領域作為IP使用。

參考文獻

[1] Schmidt R O.Multiple emitter location and signal parameter estimation[J].Antennas and Propagation，1986，34（3）：276-280.

[2] Stoica P，Nehorai A.MUSIC，maximum likelihood，and Cramer-Rao bound[M].IEEE Transactions on Acoustics，Speech，and Signal Processing，1989.

[3] 袁生光，沈海斌.基于Jacobi算法對稱矩陣特征值計算的FPGA實現[J].機電工程，2008，25（10）：80-82.

[4] 劉永勤.對稱矩陣特征值分解的FPGA實現[J].現代電子技術，2017（12）：15-18.

[5] 王飛，王建業(yè)，張安堂，等.實對稱矩陣特征值分解高速并行算法的FPGA實現[J].空軍工程大學學報（自然科學版），2008，9（6）：71-74.

FPGA Implementation of? Eigenvalue Decomposition for Floating-point

Hermite Matrix

WANG Chuan-gen，LV Le-qun，LIN Du-du

（The 29th Research Institute of? CETC， Chengdu? Sichuan? 610000）

Abstract：In order to improve the speed and accuracy of eigenvalue decomposition of Hermite matrix， a parallel reconfigurable FPGA solution based on QR iteration was implemented. For verifying this design， eigenvalue decomposition of 8×8 Hermite matrix was realized by calling hard floating-point DSP on Intels 10AX066H1F341 FPGA. The verification results show that the calculation takes 11425 clock cycles， the error is less than 3.74e-7， the resource consumption is less than 32%， and the maximum running frequency can reach 311MHz， which can support the high-speed hardware engineering implementation of complex algorithms such as MUSIC.

Key words：Hermite matrix;eigenvalue decomposition;reconfigurable;FPGA