一種厄米特矩陣特征值并行求解方法

嚴新文 彭永健 張續瑩

摘 要：厄米特矩陣特征值特征向量快速求解是在實際工程中常常會碰到的問題，但大部分基于計算機實現的算法都是利用串行思想進行編程求解，不能滿足實時求解需求。本文提出基于FPGA的并行算法來求解本類問題。首先，將厄米特矩陣實數化。其次，通過并行雅克比分解和特征值快速排序來恢復厄米特矩陣的特征值特征向量。最后，通過MATLAB仿真驗證了算法的有效性，并通過蒙特卡洛仿真確定了雅克比sweep的迭代次數，同時通過FPGA實現驗證了算法的實時性。

關鍵詞：厄米特矩陣 特征值 并行雅克比分解 蒙特卡洛仿真

中圖分類號：O129 文獻標識碼：A 文章編號：1672-3791（2018）04（b）-0108-02

矩陣特征值特征向量求解一直是線性代數中最重要的操作之一，在科學計算中有著非常廣泛的應用[1]。工程技術和物理等領域中的很多問題在數學上都歸結為求矩陣的特征值特征向量問題。例如振動問題和物理學中某些臨界值的確定，均歸結為求矩陣的特征值問題。

矩陣特征值和特征向量的求解方法有很多，冪法可以求解一般矩陣絕對值最大的特征值，雅可比方法可以求解對稱矩陣的全部特征值，QR方法可以求解一般實矩陣的全部特征值[2]。在一般的工程當中，常用的矩陣特征分解方法大部分都是針對實矩陣的，而在實際工程中，我們常需要處理復矩陣，而且大部分情況下為一個厄米特矩陣，由于厄米特矩陣具有對稱性，能簡化其特征值特征向量求解過程。

但是上述算法的計算機實現大部分都是利用串行思想進行編程求解，并不適合用于并行計算，本文詳細介紹一種可用于FPGA并行高效實現的厄米特矩陣特征值特征向量求解方法[3]。

1 厄米特矩陣特征值求解方法

1.1 實數化

如果直接對厄米特矩陣進行特征分解顯然是不合適的，也是沒有效率的，因此我們在對其進行特征分解前先進行一項預處理過程，即將它實數化。

1.2 并行雅克比分解

對于任意的一個實對稱矩陣A，只要能夠求得一個正交矩陣U，使得 成為一個對角陣D，就得到了A的所有特征值和對應的特征向量。雅可比法就是基于這一思想，用一系列的初等正交變換逐步消去A的非對角線元素，從而使矩陣A對角化。

基本雅可比法需要逐次尋找非主對角元中絕對值最大元素，這種方法并不適合于并行運算，因此我們可以按照一定的順序將A中的所有上三角元素全部消去一遍，這樣的過程稱為一次掃描（sweep）。由于對稱性，在一次掃描中，A的下三角元素同時被消去了一遍，被消為0的元素在消去后面元素的旋轉過程中可能又變為非0。但經過若干輪之后，可使S收斂于一個對角陣。

并行雅克比分解，通過構造旋轉矩陣，可以一次將矩陣N階矩陣A上三角中的N/2個點消為0，經過（N-1）次旋轉變換可將上三角中的每個點都消一次，使矩陣A接近對角陣，該過程作為一個迭代過程。為滿足一定的精度要求，可循環迭代多次，迭代次數通過蒙特卡洛仿真確定。

1.3 特征值排序

排序邏輯是結合冒泡排序和歸并排序的分組處理思想而設計的。冒泡排序是每次選出當前序列中的最大或最小元素，重復這個過程，直到輸入序列中的元素全部有序輸出時結束。歸并排序算法首先將輸入序列分為若干個組，每組內排序，對這些已排序的有序子序列再次分組，同組內的有序序列歸并為一個有序序列，重復這個過程直到所有元素有序排列。

1.4 厄米特矩陣特征值求解方法

根據上文所述，按照并行雅可比法計算厄米特矩陣特征值和特征向量的整體實現流程步驟如下：

（1）將輸入的N階復數矩陣B升維為2*N階實數對稱矩陣A。

（2）根據雅可比基本原理的公式計算旋轉矩陣R_temp。

（3）計算特征值，即計算R_temp'*A*R_temp，該步驟可將矩陣A的N個點消為0，經過多次循環迭代后，矩陣A接近對角陣，對角陣上的元素即是特征值。

（4）計算特征向量，將每次迭代的旋轉矩陣R_temp累乘，最終得到特征向量。

（5）將特征值按照從小到大順序排序，選擇奇數位置的特征值，同時選擇對應的特征向量，恢復出N階復數特征向量。

2 仿真

本文使用matlab蒙特卡洛仿真來驗證并行雅可比算法的有效性。仿真中構造隨機8*8復數赫爾米特矩陣R，通過設置不同的迭代次數，分別使用上述并行雅可比算法進行特征值分解，以及matlab系統自帶函數eig特征值分解，計算兩者特征值的差值，并取差值的最大值作為絕對誤差。設置仿真循環次數10萬次，統計每次仿真的絕對誤差，繪制絕對誤差直方圖。

設置雅可比sweep迭代次數為4，在此情況下，絕對誤差最大可達30，大部分誤差都集中在2以內。設置雅可比sweep迭代次數為5，在此情況下，絕對誤差最大為5×10^（-6），大部分誤差都集中在10^（-6）以內。設置雅可比sweep迭代次數為6，在此情況下，絕對誤差最大為7×10^（-10），大部分誤差都集中在10^（-9）以內。為使絕對誤差達到10^（-6），在后續實現過程中迭代次數設置為5。

在Vivado 2016.2的環境下，使用VHDL按照單精浮點模型完成代碼開發，然后編寫測試激勵信號，在Modelsim環境下仿真如圖5。從輸入矩陣A，到完成特征值和特征向量計算，共需要11576個時鐘節拍。在系統時鐘為150MHz情況下，所需時間為11576×6.66ns=77.1us。

將Modelsim仿真得到的結果與matlab計算結果對比，可以看出結果基本一致，誤差為10^（-6）。

3 結語

本文提出了一種通過將厄米特矩陣實數化，然后再通過并行雅克比分解，特征值快速排序來恢復厄米特矩陣的特征值特征向量的方法。最后通過MATLAB仿真驗證了算法的有效性，并通過蒙特卡洛仿真確定了雅克比sweep的迭代次數，同時通過FPGA實現驗證了算法的實時性。

參考文獻

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

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

[3] 王濤.求矩陣特征值的GPU并行算法的研究[D].黑龍江：黑龍江大學，2012.