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

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

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

中圖分類號：O129 文獻(xiàn)標(biāo)識碼：A 文章編號：1672-3791（2018）04（b）-0108-02

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

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

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

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

1.1 實(shí)數(shù)化

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

1.2 并行雅克比分解

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

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

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

1.3 特征值排序

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

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

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

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

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

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

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

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

2 仿真

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

設(shè)置雅可比sweep迭代次數(shù)為4，在此情況下，絕對誤差最大可達(dá)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)。為使絕對誤差達(dá)到10^（-6），在后續(xù)實(shí)現(xiàn)過程中迭代次數(shù)設(shè)置為5。

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

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

3 結(jié)語

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

參考文獻(xiàn)

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

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

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