嚴新文 彭永健 張續(xù)瑩

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

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

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

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

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

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

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

1.1 實數(shù)化

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

1.2 并行雅克比分解

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

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

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

1.3 特征值排序

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

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

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

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

（2）根據(jù)雅可比基本原理的公式計算旋轉(zhuǎn)矩陣R_temp。

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

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

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

2 仿真

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

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

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

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

3 結(jié)語

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

參考文獻

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

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

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