一種實對稱矩陣特征值快速求解方法研究

彭永健 廖曉慧 張鵬飛

摘要：實對稱矩陣特征分解在工程項目中經常遇到需要快速求解特征值和特征向量，文章提出了一種基于FPGA的高效并行實對稱矩陣特征值分解方法。通過構造旋轉矩陣一次消除矩陣多個元素，讓矩陣快速收斂為對角陣，從而得到原實對稱矩陣的特征值和特征向量。并根據工程誤差需求通過Matlab仿真確定了算法循環次數，在FPGA上驗證了算法的可行性和快速性。

關鍵詞：實對稱矩陣；特征值分解；快速收斂；旋轉矩陣

中圖分類號：TP3? 文獻標志碼：A

0 引言

在工程應用中，經常需要計算實對稱矩陣的特征值和特征向量，特別是信號子空間相關的陣列算法中，特征分解是其中的關鍵步驟［1］。針對實對稱矩陣，通常思路是采用正交變換法計算得到全部特征值和特征向量［2］，即通過多次正交相似變換將實對稱矩陣轉化為相似矩陣，根據計算過程可以得出原矩陣的特征值和特征向量。

QR分解算法無需將矩陣完全化簡，利用Householder變換，將原矩陣化簡為對稱三對角陣，再通過QR分解迭代，求得特征向量以及特征值［3］。雅可比分解算法通過一系列平面旋轉矩陣迭代，最終將原矩陣轉化為對角陣［4］。大部分實現算法都是非線性運算，而FPGA的優勢主要是并行線性運算帶來運算速度的提升。本文介紹一種基于FPGA的高效并行實對稱矩陣特征值分解方法。

1 特征值分解計算方法

1.1 單點消除的雅克比分解

兩個相似矩陣具有相同的特征值和特征向量［5］。根據這一基本原理，雅可比分解算法通過構造旋轉矩陣，多次正交相似變換后，將N階實對稱矩陣A的非對角線元素逐漸變換為零，最終趨近于對角陣，此時對角線元素即為特征值，所有旋轉矩陣可計算得到特征向量［6］。

針對N階實對稱矩陣A構造旋轉矩陣，對角線上的元素只有pii和pjj為cosθ，其他位置為1；上三角的元素中，只有pij為-sinθ，其他位置為0；下三角的元素中，只有pji為sinθ，其他位置為0；m=-apq，n=0.5×（aqq-app），sinθ=sign（n）×sign（m）×m2/（m2+n2）2［1+n2/（m2+n2）］，cosθ=1-m2/（m2+n2）2［1+n2/（m2+n2）］。

通過矩陣計算A1=PTAP，可知該運算只改變原矩陣A的第i行、第i列、第j行、第j列，并且aij和aji變為0。繼續該過程，依次將上三角中的元素消為0，共需要N×（N-1）/2次運算可以遍歷一次所有非對角線元素，該過程記為一個循環LOOP。在一個LOOP中，將后面的元素消為0的過程會將前面已經消為0的元素改變，但整個LOOP結束后，對角線元素的平方和增大，非對角線元素的平方和減小。經過若干個LOOP后，矩陣逼近于一個對角陣。此時，對角線上的元素就是原矩陣A的特征值，所有的旋轉矩陣P相乘就是特征向量矩陣。

1.2 多點消除的雅克比分解

上述串行分解運算過程中，一個LOOP中每一步僅消去上三角中的一個元素，PTAP涉及兩次矩陣乘法運算，共需要N×（N-1）次矩陣乘法，在FPGA中實現該過程耗時較長。

消去上三角中元素aij時，只改變原矩陣A的第i行、第i列、第j行、第j列，而消去上三角中元素akl計算旋轉矩陣時，僅需要用到akk，all，akl，因此可以同時消aij和akl，只需要滿足i，j，k，l互不相等即可。對于N階矩陣A，可以構造一個矩陣P，一次將上三角中的N/2個元素消為0，需要（N-1）次運算即可將非對角線元素遍歷一次，矩陣乘法的次數降為2（N-1）。

1.3 實對稱矩陣特征值和特征向量求解方法

根據上文所述，以16階實對稱矩陣為例，并行雅可比分解法計算特征值和特征向量的實現步驟如下。

（1）根據公式計算旋轉矩陣P11。

（2）計算A11=PT11AP11，將8個元素消為0。

（3）重新選定8個元素位置，計算旋轉矩陣P12。

（4）計算A12=PT12A11P12。

（5）重復步驟（3）和（4），直到上三角中的120個元素全都完成一次消除，一個循環LOOP內的消除順序如下表1所示，表中（x，y）代表第x行、第y列。

（6）重復進行M個LOOP，直到對角線元素趨于0。

（7）對角線元素即為特征值，所有的旋轉矩陣相乘，最終得到的矩陣每一列都是對角陣相應列特征值的特征向量。

2 Matlab仿真

針對以上算法，編寫Matlab程序，仿真每次消單點和消多點在相同的循環迭代次數情況下，收斂速度是否有區別，并以此確定循環迭代次數LOOP的值。收斂速度按照絕對誤差為判斷標準，將Matlab系統函數Eig的計算結果作為理論值，雅可比分解算法得到的特征值與理論值的差值作為絕對誤差。針對不同的LOOP值，隨機生成16階實對稱矩陣，仿真次數為10 000。

由仿真結果可知，當LOOP值為6時，消單點的計算結果誤差為10^（-6），消多點的計算結果誤差為10^（-9），可以看出消多點雅可比算法能夠更快的收斂到對角陣，得到與理論值誤差更小的特征值。并且，10^（-9）的精度能滿足大部分工程應用場景，后續的FPGA實現按照LOOP為6實施。

3 FPGA實現

在FPGA中實現并行雅可比分解算法如圖1所示。每次計算旋轉矩陣P后，需要計算正交變換A1=PTAP，并且將矩陣P累乘。正交變換和矩陣累乘都需要用到矩陣乘法操作，設計時將矩陣累乘和計算旋轉矩陣并行，能夠節省整個循環迭代的時間，即計算出旋轉矩陣后，先進行正交變換，再同時開展矩陣累乘和計算下一個旋轉矩陣的工作。并且，這樣操作能夠復用矩陣乘法資源，僅例化一套矩陣乘法運算資源，交替進行正交變換運算和變換矩陣累乘運算。

整個實現過程中，資源消耗量大的是矩陣乘法運算。以16階矩陣乘法B=A×P為例，計算B中的1個點需要16個乘法器、15個加法器，計算所有點共需要4 096個乘法器、3 840個加法器，DSP資源消耗巨大，無法在一片FPGA上實現。因此，需要做適當的串行化，降低資源消耗。經過時間和資源的折中考慮，在本設計中例化16套資源，并行計算B中的16個點，再重復使用16次完成B中所有點的計算。資源消耗為1 248個DSP，計算時間為75個時鐘周期。

在Modelsim中對開發的算法程序進行仿真。從輸入矩陣A，到計算出特征值和特征向量（LOOP次數為6），共需要29 000個時鐘周期，使用300 MHz的計算時鐘，所需時間為96.6 us，絕對誤差最大僅10^（-6），滿足工程應用需求。

4 結語

本文提出了一種高效并行實對稱矩陣特征值分解方法，通過Matlab仿真確定了循環次數，并在FPGA上驗證了算法的有效性。

參考文獻

［1］陳建華.線性代數［M］.4版.北京：機械工業出版?? 社，2017.

［2］徐士良，馬爾妮.常用算法程序集（C/C++描述）［M］.北京：清華大學出版社，2013.

［3］郝英軍.MUSIC算法中特征值分解及信源數估計的FPGA實現［D］.哈爾濱：哈爾濱工程大學，2020.

［4］胡樂宇，蔡邢菊.低計算精度下實對稱矩陣的特征值分解［J］.高等學校計算數學學報，2021（2）：117-133.

［5］胡樂宇.線性方程組的隨機求解算法以及低精度矩陣的特征值分解［D］.南京：南京師范大學，2020.

［6］唐浩偉，羅林，王濤.基于特征值分解的時間反轉超聲成像技術［J］.信息技術，2019（9）：52-55.

（編輯 沈 強）

Research on a fast method for solving eigenvalues of real symmetric matrix

Peng? Yongjian， Liao? Xiaohui， Zhang? Pengfei

（The 29th Research Institute of China Electronics Technology Group Corporation， Chengdu 610029， China）

Abstract：? Eigendecomposition of real symmetric matrices is often encountered in engineering projects， and it is necessary to quickly solve eigenvalues and eigenvectors. This paper proposes an efficient parallel real symmetric matrix eigenvalue decomposition method based on FPGA. By constructing a rotation matrix to eliminate multiple elements of the matrix at one time， the matrix can quickly converge to a diagonal matrix， and the eigenvalues and eigenvectors of the original real symmetric matrix can be obtained. According to the engineering error requirements， the cycle times of the algorithm are determined by Matlab simulation， and the feasibility and rapidity of the algorithm are verified on FPGA.

Key words： real symmetric matrices; eigenvalue decomposition; quickly converges; rotation matrix