基于擴展FEAST的大規模特征值求解問題研究

李玲玲 李 華

(河南城建學院數理學院 河南 平頂山 467036)

0 引 言

非線性特征值問題[1](Nonlinear Eigenvalue Pro-blems,NLEVP)是計算科學的熱門研究課題。其矩陣值函數的有用性，都將受到問題求解難度的限制，伴隨著一些特有難題，例如其特征向量通常并不構成n維向量空間的基礎，且矩陣值函數的特定形式會影響到求解方法的效率。然而，非線性特征值問題與線性特征值問題有著相同的核心難題：在n的數值較小時有效的求解方法并不一定能很好地擴展到高維問題[2]。此類問題通常側重于尋找少量特定的、有重要物理學意義的“特征值-特征向量”對。

大部分求解此類高維問題的方法需要生成一個與期望特征對足夠接近的初始猜測，以確保收斂性。同時，需要使用能夠計算特征值接近的連續特征對，但不會向同一個特征對反復收斂的方法[3]。尤其在期望得到大量特征對的情況下，以上問題更加難以解決。NLEVP算法進一步加大了計算大量特征對的難度，導致不能很好地適應在子空間內求解的降維后的特征值問題。而求解非線性特征值問題所需的計算量至少與求解單個線性特征值問題相等。因此，保持較小的維數至關重要，以確保能夠適應較大規模的問題。

圍道積分法是解決上述難題的一類算法[4-5]。圍道積分法利用柯西積分公式，從而僅計算其特征值位于復平面的用戶定義特定區域中的特征對。這樣就不需要計算初始猜測，并對維數為n的原始問題中距離很近的特征值進行了分離。此外，圍道積分法支持大量特征對的高效并行計算，通過將復平面中的感興趣區域分割為不相交的子區域，在每個子區域中獨立求解特征值，且與其他子區域的特征值無關。

而線性FEAST算法[6-7]是求解線性特征值問題的圍道積分法樣例。該算法能夠求解大規模廣義線性特征值問題，具有較好的魯棒性和優秀的并行可擴展性。本文對線性FEAST算法進行推廣，使其能夠求解非線性特征值問題，并通過求解一些物理模型問題來解釋該算法的特性。

1 非線性特征值問題

給定一個非空的開集D?C和一個矩陣值函數T:D→Cn×n，考慮一個非線性特征值問題：找到標量λ∈D和非零向量x∈Cn、y∈Cn，使得：

T(λ)x=0 且y·T(λ)=0

(1)

本文將(λ,x,y)稱為T的特征三元組，(λ,x)或(λ,y)均被稱為特征對。對于T(λ)=λI-A和T(λ)=λB-A，通過式(1)分別將其簡化為標準特征值問題Ax=λx和廣義特征值問題Ax=λBx。廣義非線性特征值問題的特性與線性特征值問題有著很大的不同，例如T(即使是正則值，即det(T(λ))≠0)可能有無窮多個特征值，屬于不同特征值的特征向量不一定是線性無關的，且一個孤立特征值的代數重數雖然是有限的，但不受問題大小n的限制[8]。所有這些特性，使得廣義非線性特征值問題非常難以求解。關于該問題更詳細的描述，可參考二次特征值問題的相關文獻或廣義非線性特征值問題[9]。

為便于討論，本文側重于研究方陣類多項式特征值問題。找到標量λ∈D?C和非零向量x,y∈Cn，以使得：

P(λ)x=0y·P(λ)=0

(2)

式中：P(λ)為n×n的矩陣多項式。

(3)

若P為正則，即det(P)≠0，則式(2)有r個有限特征值(det(P)=0的根)和kn-r個無限特征值(復特征值)。通過線性化，可知k次的n×n矩陣多項式P(λ)有著kn個特征值(有限或無限)，最多kn個右相關特征向量，及最多kn個左相關特征向量。

2 擴展FEAST算法

2.1 FEAST算法

原始FEAST算法被用于求解廣義特征值問題，即T(λ)=λB-A，其是一個子空間迭代方法，使用Rayleigh-Ritz投影和近似光譜投影器作為濾波器。FEAST計算與位于復平面中用戶定義的特定區域中的特征值相對應的特征對。因此，該算法可被用于并行計算大量特征對，其中將復平面分割為不相交的區域集合，并獨立求解每個區域中的特征對(與其他區域中的特征對無關)。

FEAST首先將任意隨機初始子空間X(0)投影到感興趣特征向量所覆蓋的子空間來選擇感興趣的特征向量和特征值，然后使用Rayleigh-Ritz程序[10]在該子空間中提取特征值/特征向量的近似。通過使用復圍道積分來完成向感興趣子空間的投影：

(4)

式中：Q為新的投影子空間，在包圍著感興趣特征值的圍道中完成積分。若精確完成積分，則濾波器ρ(A,B)成為光譜投影器，ρ(A,B)=XYHB，其中X和Y分別是與感興趣特征值(即位于閉合曲線內的特征值)相關的精確的右特征向量矩陣和左特征向量矩陣。然而在實踐中，無法在式(4)中精確完成積分，因此使用求積規則對其進行逼近。則FEAST算法實際生成的子空間為：

(5)

式中：zj和wj分別為積分節點和權值，通過求解線性系統(zjB-A)Qj=BX(0)來完成求和。由于每個線性系統均獨立于其他系統，因此可以并行地同時求解。FEAST通過反復應用求和，將得到的子空間Q正交化，并將其用于Rayleigh-Ritz程序中，對感興趣特征向量的估計進行迭代式改進。

本文旨在利用FEAST框架的優點來求解NLEVP,使用圍道積分對固定維數的子空間進行迭代優化，以求解特征值位于復平面的用戶定義區域中的特征對。

2.2 NLFEAST算法

線性FEAST算法可以解釋為在復平面上使用多個移位的平移-逆迭代的推廣。利用相似的推理，假定非線性FEAST(NLFEAST)算法是使用多個移位的非線性移位-逆迭代的推廣：

x(i+1)=T(z)-1T′(λ(i))x(i)

(6)

式中：x(i)和λ(i)為特征向量和特征值在迭代i處的估計;T′(·)為T(·)的導數;z為復平面中位于感興趣特征值附近的移位。由此預期非線性FEAST圍道積分(單個向量)將類似于：

(7)

積分圍道包圍著感興趣特征值。但在實踐中，式(7)中的更新程序未能向正確的特征值-特征向量對收斂。這是因為在式(7)中使用恒定移位z會造成非線性逆迭代向錯誤的解收斂。

原始FEAST算法可以解釋為移位-逆迭代的推廣，使用圍道積分以高效地同時使用多個移位。而非線性FEAST算法則可以解釋為殘差逆迭代的推廣，使用圍道積分以高效地同時使用多個移位。本文給出了非線性FEAST算法的算法框架，如算法1所示。

算法1T(λ)x=0時的NLFEAST算法

步驟1設子空間Q=X(0)，對Q的列向量進行正交歸一化。

步驟2求解投影后的非線性特征值問題:

QHT(λ)Qy=0

(8)

以得到近似特征對(λ,Qy)。

Λ=diag(λ1,λ2,…,λm0)X=[Qy1,Qy2,…,Qym0]

(9)

步驟5通過執行數值積分對搜索子空間進行更新：

(10)

使用求積規則并求解線性系統：

(11)

步驟6對新的搜索子空間Q的列向量進行正交歸一化(例如使用QR分解)，并回到步驟2。

非線性FEAST算法中Rayleigh-Ritz步驟與線性算法中的工作原理相同：將殘差函數T(λ)投影到子空間Q上，然后選擇合適的方法(本文使用了簡單線性化)對維數降低后的非線性特征值問題進行求解。得到的Ritz值及Ritz向量被用作期望特征對的新估計。

用于Rayleigh-Ritz程序的子空間，即通過執行圍道積分所生成的子空間，應該是正交歸一化的(如使用QR分解)。正交歸一化能夠防止期望特征向量的范數出現發散，并防止Rayleigh-Ritz程序生成偽特征對(即與全尺寸特征值問題的任意特征對均不對應的Ritz對)，由此提升了FEAST迭代的數值穩定性。

在m0=1和nc=1的情況下，算法1大致等價于殘差逆迭代。這兩種算法的區別在于非線性特征值和特征向量的估計方式不同，以及步驟5中更新向量的縮放差別。算法1通過求解步驟2中的投影非線性特征值問題來計算近似特征向量和近似特征值，而殘差逆迭代則使用步驟5中的子空間更新程序對近似特征向量進行更新。

3 數值分析實驗

將本文的NLFEAST算法在3個有實踐意義的多項式特征值問題上做數值分析實驗。

3.1 樣例1：非過阻尼質點彈簧系統

考慮一個與質點彈簧系統[11]相關的自共軛(Hermitian)二次特征值問題：

T(λ)x=(λ2A2+λA1+A0)x=0

(12)

式中:

參照文獻[12]，首先分析非過阻尼系統。本文設n=1 000，并選擇τ=0.620 2，k=0.480 7。圖1展示了Q的所有特征值。要考慮計算的特征向量的特征值沒有虛部；質點彈簧系統原運動公式的解的指數增長或衰減是這些特征向量的線性組合。

圖1 非過阻尼質點彈簧系統的所有特征值

使用算法1中的非線性FEAST方法，計算這20個實特征值的近似。選擇的圍道為一個橢圓，其中心為區間(-1.6,-1.5)的中點，使用nc=16個求積節點，實軸方向的半徑ra=0.05，虛軸方向的半徑rb=0.003 5，殘值的收斂標準為tol=10-10。在維數m0=22的子空間中，非線性FEAST算法在三次迭代中找到了區間內的全部m=20個特征值。表1給出了使用線性化方法和非線性FEAST算法得到的特征值近似。為完備性起見，本文還列舉了所有的最終殘值，即‖ri‖2=‖T(λi)xi‖2。

表1 非過阻尼質點彈簧系統在的所有實特征值列表

續表1

表2給出了當搜索區間(a,b)增大時，收斂所需的非線性FEAST迭代次數。當求積節點nc的數量保持在16個不變時，擴大搜索區間的影響等同于將求積節點放置在遠離感興趣特征值的地方，由此使得圍道積分的準確度降低。區間大小的適度變化僅會造成收斂所需的工作量稍微上升；而區間大小的較大變化往往會要求非線性FEAST迭代次數隨之按比例增加。實踐中，與期望特征值分布相關的搜索區間的大小而造成的性能降低，可以通過簡單地使用大量求積節點而緩解，這種情況下需要更多的并行處理能力來同時求解更多的線性系統。

表2 非過阻尼質點彈簧系統的特征值數量

3.2 樣例2：一維散射共振

非線性特征值問題的另一個應用是開邊界量子透射問題。在開邊界量子投射問題中，需要計算量子勢的準束縛態，其中粒子可以進入或離開系統，最好不要對外部源進行顯式建模。這可以通過求解非線性特征值問題來實現?？紤]一個開邊界量子透射問題的樣例，即一維散射共振問題，其中包含以下緊支有限二乘模型：

(13)

在n=302個點的網格上的有限元離散化得到大小為302的帶散射共振的二次特征值，其特征值結果如圖2所示?？梢钥闯?，r=15.5為半徑的圓內共計22個復散射共振。非線性FEAST使用nc=16個求積節點，求解子空間維數m0=30、V0=10的散射共振問題。4次迭代計算這22個散射共振的近似值(收斂容差10-10)，其結果如圖3所示。

圖2 散射共振的所有特征值

圖3 一維散射共振的NLFEAST算法特征值估計

3.3 樣例3：軌道振動

較長軌道的振動可建模為以下形式的二次特征值問題[13]：

T(λ)=λ2I+λ(I+A2)+A2+A+I

(14)

式中：A為一個n×n的轉換矩陣。對于任意維數n，通過求解以下二次方程給出式(14)的二次特征值問題的特征值：

(15)

式中：μk=-4sin2((k-1)π/n)(1≤k≤n)為矩陣A的特征值，特征值以復共軛對的形式出現。

考慮A的維數為n=50 000，使用NLFEAST算法計算與區間[-7.119 2,-6.965 0]內的特征值相對應的250個特征向量，其中圓形圍道的圓心坐標為(-7.042，0)，半徑0.077 1。圖4展示了其特征值位置和積分圍道。

圖4 軌道振動特征圖示

為進一步比較，本文還使用Beyn法利用相同的積分圍道計算相同的250個特征值。表3給出了Beyn法在使用不同的子空間維數和各種數量的求積節點時，與積分圍道內的特征值相對應的特征向量的最終特征向量殘值。與NLFEAST算法一樣，僅考慮殘值低于10-2的特征向量為正確特征向量，由此排除偽特征值。

表3 使用Beyn法計算特征值時的向量殘值

表3的結果表明，Beyn法和NLFEAST算法在計算一個解時需要的計算量大致相當。例如，當m0=500，nc=8時，NLFEAST算法需要6次迭代以收斂至約10-11，即總計500×8×6=24 000個線性系統解;Beyn法在m0=500時需要nc=64個求積節點以收斂至殘差10-12的特征向量，即對應于500×8×8=28 500個線性系統解。

以上兩種方法的差異源于計算方式的不同。Beyn法理論上可以并行求解所有相關的線性系統，而NLFEAST算法則可以并行求解最多m0×nc個線性系統。Beyn法的并行性較高，但需要計算更多數量的矩陣分解。在上述案例中，Beyn法需要64次矩陣分解，而NLFEAST算法僅需8次。

4 結 語

本文提出了線性FEAST特征值算法的一個擴展算法，以求解非線性特征值問題。提出的NLFEAST算法使用與線性FEAST算法相同的一系列運算，但修改了圍道積分，支持在求解非線性特征值問題時使用固定移位集合和固定子空間維數。線性FEAST算法可被解釋為移位和逆迭代的使用多次移位的推廣，而NLFEAST算法則可被解釋為殘差逆迭代的使用多次移位的推廣。本文改進數值圍道積分提升了近似特征向量子空間的維數，從而系統地提高了算法的收斂速度。