一種實(shí)對(duì)稱矩陣特征值快速求解方法研究

彭永健 廖曉慧 張鵬飛

摘要：實(shí)對(duì)稱矩陣特征分解在工程項(xiàng)目中經(jīng)常遇到需要快速求解特征值和特征向量，文章提出了一種基于FPGA的高效并行實(shí)對(duì)稱矩陣特征值分解方法。通過(guò)構(gòu)造旋轉(zhuǎn)矩陣一次消除矩陣多個(gè)元素，讓矩陣快速收斂為對(duì)角陣，從而得到原實(shí)對(duì)稱矩陣的特征值和特征向量。并根據(jù)工程誤差需求通過(guò)Matlab仿真確定了算法循環(huán)次數(shù)，在FPGA上驗(yàn)證了算法的可行性和快速性。

關(guān)鍵詞：實(shí)對(duì)稱矩陣；特征值分解；快速收斂；旋轉(zhuǎn)矩陣

中圖分類號(hào)：TP3? 文獻(xiàn)標(biāo)志碼：A

0 引言

在工程應(yīng)用中，經(jīng)常需要計(jì)算實(shí)對(duì)稱矩陣的特征值和特征向量，特別是信號(hào)子空間相關(guān)的陣列算法中，特征分解是其中的關(guān)鍵步驟［1］。針對(duì)實(shí)對(duì)稱矩陣，通常思路是采用正交變換法計(jì)算得到全部特征值和特征向量［2］，即通過(guò)多次正交相似變換將實(shí)對(duì)稱矩陣轉(zhuǎn)化為相似矩陣，根據(jù)計(jì)算過(guò)程可以得出原矩陣的特征值和特征向量。

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

1 特征值分解計(jì)算方法

1.1 單點(diǎn)消除的雅克比分解

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

針對(duì)N階實(shí)對(duì)稱矩陣A構(gòu)造旋轉(zhuǎn)矩陣，對(duì)角線上的元素只有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）］。

通過(guò)矩陣計(jì)算A1=PTAP，可知該運(yùn)算只改變?cè)仃嘇的第i行、第i列、第j行、第j列，并且aij和aji變?yōu)?。繼續(xù)該過(guò)程，依次將上三角中的元素消為0，共需要N×（N-1）/2次運(yùn)算可以遍歷一次所有非對(duì)角線元素，該過(guò)程記為一個(gè)循環(huán)LOOP。在一個(gè)LOOP中，將后面的元素消為0的過(guò)程會(huì)將前面已經(jīng)消為0的元素改變，但整個(gè)LOOP結(jié)束后，對(duì)角線元素的平方和增大，非對(duì)角線元素的平方和減小。經(jīng)過(guò)若干個(gè)LOOP后，矩陣逼近于一個(gè)對(duì)角陣。此時(shí)，對(duì)角線上的元素就是原矩陣A的特征值，所有的旋轉(zhuǎn)矩陣P相乘就是特征向量矩陣。

1.2 多點(diǎn)消除的雅克比分解

上述串行分解運(yùn)算過(guò)程中，一個(gè)LOOP中每一步僅消去上三角中的一個(gè)元素，PTAP涉及兩次矩陣乘法運(yùn)算，共需要N×（N-1）次矩陣乘法，在FPGA中實(shí)現(xiàn)該過(guò)程耗時(shí)較長(zhǎng)。

消去上三角中元素aij時(shí)，只改變?cè)仃嘇的第i行、第i列、第j行、第j列，而消去上三角中元素akl計(jì)算旋轉(zhuǎn)矩陣時(shí)，僅需要用到akk，all，akl，因此可以同時(shí)消aij和akl，只需要滿足i，j，k，l互不相等即可。對(duì)于N階矩陣A，可以構(gòu)造一個(gè)矩陣P，一次將上三角中的N/2個(gè)元素消為0，需要（N-1）次運(yùn)算即可將非對(duì)角線元素遍歷一次，矩陣乘法的次數(shù)降為2（N-1）。

1.3 實(shí)對(duì)稱矩陣特征值和特征向量求解方法

根據(jù)上文所述，以16階實(shí)對(duì)稱矩陣為例，并行雅可比分解法計(jì)算特征值和特征向量的實(shí)現(xiàn)步驟如下。

（1）根據(jù)公式計(jì)算旋轉(zhuǎn)矩陣P11。

（2）計(jì)算A11=PT11AP11，將8個(gè)元素消為0。

（3）重新選定8個(gè)元素位置，計(jì)算旋轉(zhuǎn)矩陣P12。

（4）計(jì)算A12=PT12A11P12。

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

（6）重復(fù)進(jìn)行M個(gè)LOOP，直到對(duì)角線元素趨于0。

（7）對(duì)角線元素即為特征值，所有的旋轉(zhuǎn)矩陣相乘，最終得到的矩陣每一列都是對(duì)角陣相應(yīng)列特征值的特征向量。

2 Matlab仿真

針對(duì)以上算法，編寫(xiě)Matlab程序，仿真每次消單點(diǎn)和消多點(diǎn)在相同的循環(huán)迭代次數(shù)情況下，收斂速度是否有區(qū)別，并以此確定循環(huán)迭代次數(shù)LOOP的值。收斂速度按照絕對(duì)誤差為判斷標(biāo)準(zhǔn)，將Matlab系統(tǒng)函數(shù)Eig的計(jì)算結(jié)果作為理論值，雅可比分解算法得到的特征值與理論值的差值作為絕對(duì)誤差。針對(duì)不同的LOOP值，隨機(jī)生成16階實(shí)對(duì)稱矩陣，仿真次數(shù)為10 000。

由仿真結(jié)果可知，當(dāng)LOOP值為6時(shí)，消單點(diǎn)的計(jì)算結(jié)果誤差為10^（-6），消多點(diǎn)的計(jì)算結(jié)果誤差為10^（-9），可以看出消多點(diǎn)雅可比算法能夠更快的收斂到對(duì)角陣，得到與理論值誤差更小的特征值。并且，10^（-9）的精度能滿足大部分工程應(yīng)用場(chǎng)景，后續(xù)的FPGA實(shí)現(xiàn)按照LOOP為6實(shí)施。

3 FPGA實(shí)現(xiàn)

在FPGA中實(shí)現(xiàn)并行雅可比分解算法如圖1所示。每次計(jì)算旋轉(zhuǎn)矩陣P后，需要計(jì)算正交變換A1=PTAP，并且將矩陣P累乘。正交變換和矩陣?yán)鄢硕夹枰玫骄仃嚦朔ú僮?，設(shè)計(jì)時(shí)將矩陣?yán)鄢撕陀?jì)算旋轉(zhuǎn)矩陣并行，能夠節(jié)省整個(gè)循環(huán)迭代的時(shí)間，即計(jì)算出旋轉(zhuǎn)矩陣后，先進(jìn)行正交變換，再同時(shí)開(kāi)展矩陣?yán)鄢撕陀?jì)算下一個(gè)旋轉(zhuǎn)矩陣的工作。并且，這樣操作能夠復(fù)用矩陣乘法資源，僅例化一套矩陣乘法運(yùn)算資源，交替進(jìn)行正交變換運(yùn)算和變換矩陣?yán)鄢诉\(yùn)算。

整個(gè)實(shí)現(xiàn)過(guò)程中，資源消耗量大的是矩陣乘法運(yùn)算。以16階矩陣乘法B=A×P為例，計(jì)算B中的1個(gè)點(diǎn)需要16個(gè)乘法器、15個(gè)加法器，計(jì)算所有點(diǎn)共需要4 096個(gè)乘法器、3 840個(gè)加法器，DSP資源消耗巨大，無(wú)法在一片F(xiàn)PGA上實(shí)現(xiàn)。因此，需要做適當(dāng)?shù)拇谢?，降低資源消耗。經(jīng)過(guò)時(shí)間和資源的折中考慮，在本設(shè)計(jì)中例化16套資源，并行計(jì)算B中的16個(gè)點(diǎn)，再重復(fù)使用16次完成B中所有點(diǎn)的計(jì)算。資源消耗為1 248個(gè)DSP，計(jì)算時(shí)間為75個(gè)時(shí)鐘周期。

在Modelsim中對(duì)開(kāi)發(fā)的算法程序進(jìn)行仿真。從輸入矩陣A，到計(jì)算出特征值和特征向量（LOOP次數(shù)為6），共需要29 000個(gè)時(shí)鐘周期，使用300 MHz的計(jì)算時(shí)鐘，所需時(shí)間為96.6 us，絕對(duì)誤差最大僅10^（-6），滿足工程應(yīng)用需求。

4 結(jié)語(yǔ)

本文提出了一種高效并行實(shí)對(duì)稱矩陣特征值分解方法，通過(guò)Matlab仿真確定了循環(huán)次數(shù)，并在FPGA上驗(yàn)證了算法的有效性。

參考文獻(xiàn)

［1］陳建華.線性代數(shù)［M］.4版.北京：機(jī)械工業(yè)出版?? 社，2017.

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

［3］郝英軍.MUSIC算法中特征值分解及信源數(shù)估計(jì)的FPGA實(shí)現(xiàn)［D］.哈爾濱：哈爾濱工程大學(xué)，2020.

［4］胡樂(lè)宇，蔡邢菊.低計(jì)算精度下實(shí)對(duì)稱矩陣的特征值分解［J］.高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào)，2021（2）：117-133.

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

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

（編輯 沈 強(qiáng)）

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