改進的實對稱陣特征值分解并行化算法

金鷹翰， 趙 培， 趙景琰， 王進祥

（哈爾濱工業大學, 黑龍江 哈爾濱 150001）

0 引言

波達方向估計技術是陣列信號處理的關鍵技術之一，其中最具代表性的算法有MUSIC算法[1-2]和ESPRIT算法[3]。這兩種算法都需要對接收矢量的協方差矩陣進行特征值分解，且該部分占算法的主要運算量，通過實數化方法可以將協方差矩陣轉換為實對稱陣，因此采用并行算法實現實對稱陣特征值分解能顯著縮短算法時間。Jacobi算法是對稱陣特征值分解的常用算法，該算法為串行算法，需要時間隨矩陣階數增長迅速增長。David J. Kuck和Ahmed H. Sameh提出了一種基于Jacobi算法的并行算法[3]，但由于算法對矩陣元素重復操作，使得算法在實際實現后，效率不是很高。在分析了串行算法和David J. Kuck和Ahmed H. Sameh的并行算法的優缺點后，本文實現了一種更為高效的實對稱矩陣特征值分解的并行算法，該算法不增加計算復雜度，經過計算機仿真，效率有很大提高。

1 基于Jacobi算法的實對稱陣特征值分解并行算法

在實對稱矩陣特征值分解的串行算法中，通常采用Jacobi算法求矩陣的特征值和特征向量。David J. Kuck和Ahmed H. Sameh提出了一種基于Jacobi算法的實對稱陣特征值分解并行算法。串行Jacobi算法每次正交旋轉變換只消去將兩個非對角線元素，而通過如下正交旋轉變換可以消去n個非對角線元素（n為矩陣階數，這里只考慮偶數階矩陣，奇數階矩陣多出的一行和一列可以單獨處理）。

為了給下一次變換做準備，需要做以下變換：

其中：即將Ak+1的第2行移到最后一行，第2行以后各行依次向前移一行；將Ak+1的第2列移到最后一列，第2列以后各列依次向前移一列。

在實現中，發現此算法效率仍不是很高，主要原因是重復消去矩陣非對角線元素。以8階矩陣為例，第一次變換消去矩陣上三角元素（由于是對稱陣，只列出上三角元素）而第三次變換消去的元素對應于消去原矩陣元素與第一次變換存在大量重復消去元素。

2 改進的實對稱陣特征值分解并行化算法

針對David J. Kuck和Ahmed H. Sameh提出的并行算法的缺點，本文提出了一種改進的實對稱陣特征值分解并行化算法。保留了原并行算法中的并行正交旋轉變換而通過修改φ矩陣，減少重復消去元素的個數，并且只使用一種變換。

φ矩陣應為n階矩陣，每行每列有且僅有一個元素為1，其余元素均為0，可以用1到n的一個排列表示。φ矩陣第i行第j列元素為1，表示將的第i行移到第j行，第i列移到第j列。為記錄經過變換后的矩陣元素與原矩陣元素的位置對應關系，構造矩陣，對所有非對角線的上三角元素進行編號。記錄元素，即正交旋轉變換消去的元素在原矩陣中的對應位置，然后經過變換A'=φTAφ。重復上述操作n-1次，如果記錄的值包含1到n(n - 1 )/2的所有值，則說明經過 n -1次變換，可以將所有非對角線元素無重復無遺漏地消去一遍。

計算φ可以通過生成1到n的所有排列,然后驗證符合上述條件的排列。滿足條件的排列有多個，如下給出6階矩陣的一種變換矩陣中值為一的元素,其余元素均為零：

3 硬件結構

圖1為根據改進的特征值分解算法實現的模塊內部結構。實對稱矩陣和特征向量矩陣使用存儲器實現，變換 A =φTAφ并不需要對存儲器進行讀寫，只需要控制數據選擇器讀存儲器的地址即可實現。數據選擇器和數據分配器由一個狀態機控制，負責訪問存儲器，數據選擇器將數據傳給運算單元，數據分配器將結果存回存儲器。運算單元完成如下運算：

正交旋轉變換中的矩陣乘法都可以分解成這種運算。

圖1 硬件模塊結構

圖2為運算單元硬件結構，該結構將該運算過程分為計算θ和完成矩陣乘法兩部分，采用CORDIC算法[4]可將正弦、余弦和乘法運算轉化為加法和移位運算，便于硬件實現。

圖2 計算單元硬件結構

圖2（a）為計算CORDIC符號集的結構，kξ為CORDIC符號集寄存器，控制右邊兩個加法器進行加運算還是減運算。圖2（b）為完成矩陣乘法的結構，同樣由kξ控制加法器進行加或減運算。運算單元最多可以并行2/2n 個，可以根據設計要求選擇。

下面分析同樣根據以上硬件結構，串行 Jacobi算法、David J. Kuck和Ahmed H. Sameh算法和本文改進算法的運算量。假定CORDIC算法的階數為m，分解n階矩陣。

串行 Jacobi算法：單次迭代需要(6n+4)m次移位和(6n+4)m次加/減運算。完成一輪迭代（即將所有上三角元素消去一遍）需要n(n-1)(3n+2)m次移位和n(n-1)(3n+2)m次加/減運算。需要4輪迭代可以滿足要求，需要4n(n-1)(3n+2)m次移位和4n(n-1)(3n+2)m次加/減運算。

David J. Kuck和Ahmed H. Sameh算法：完成一輪迭代需要n(n-1)次單次迭代，一般需要2輪迭代，所以總共需要2n2(n-1)(3n+2)m次移位和2n2(n-1)(3n+2)m次加/減運算。

改進算法：改進算法需要4論迭代，單次迭代的運算量與David J. Kuck和Ahmed H. Sameh算法相同，總運算量與串行Jacobi算法相同。可見當n＞4時，改進算法的運算量要小于David J. Kuck和Ahmed H. Sameh的算法，且n越大，改進效果越明顯。

4 計算機仿真

在MUSIC算法仿真中，分別使用David J. Kuck和Ahmed H. Sameh提出的并行算法和本文提出的改進的并行算法進行特征值分解，兩者均能得到正確的結果。表1為不同天線陣元數MUSIC算法20次仿真的平均統計數據，天線采用均勻圓陣模型，信號源數均為 3，采用隨機正整數序列輸入，最大幅值分別為2、4、8，精度為10-4。

表1 仿真統計數據

5 結語

本文在實現波達方向估計MUSIC算法時，針對占主要運算量的實對稱陣特征值分解提出了一種改進的并行算法。該算法保留了David J. Kuck和Ahmed H. Sameh提出的并行算法中并行化的思想，同時解決了其重復消去矩陣元素的缺點。由于該算法增加的矩陣變換不需要對存儲器進行讀寫，只需要控制數據選擇器讀存儲器的地址即可實現，因此便于硬件實現，控制較簡單，不增加計算量，效率得到了明顯提高。由仿真數據可以看出，改進效率隨著矩陣階數增加而提高。

[1] Schmidt R O. Multiple Emitter Location and Signal Parameter Estimation[J].IEEE Transactions on Antennas and Propagation,1986,34(03):276-280．

[2] 張家平,劉青,張洪順.低信噪比中 MUSIC算法研究[J].通信技術,2009,42(01):87-89．

[3] 龍娟,周圍,吳江華.TD-SCDMA中DOA估計的一種ESPRIT算法[J].通信技術,2007,40(09):31-33．

[3] Kuck D J,Lawire D H,Sameh A H.High Speed Computer and Algorithm Organization[M].New York:Academic Press,1977:71-84．

[4] Kim M,Ichige K,Arai H.Design of JACOBI EVD Processor Basedon CORDIC for DOA Estimation with MUSIC Algorithm[C]//The 13th IEEE International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC 2002).Lisbon,Portugal:IEEE, 2002:120-124．