[矩陣計算]Lanczos方法:求稀疏矩陣特征值


更新: 29 JUL 2016

QR方法知,求矩陣$A$的特征值,大多需要先將其三對角化(詳細方法見徐樹方先生的教材。此處外鏈一個例子),即

$$ T=Q^TAQ $$

即找到正交矩陣$Q$使得$T$成為三對角矩陣。然而若$A$為大型稀疏矩陣,常用的方法如Householder和Givens變換都無法充分利用$A$的稀疏性,因此考慮直接計算$T$和$Q$的矩陣元以利用$A$的稀疏性加速運算。

 

一、Lanczos方法基本原理

將以上分解式中的$Q$寫成

$$ Q=[q_1,q_2,\cdots,q_n] $$

其中$ q_i $為$Q$的列向量。$T$寫成

$$ T=\begin{bmatrix} \alpha_1 & \beta_1 & & & 0 \\ \beta_1 & \alpha_2 &\ddots & & \\ & \ddots & \ddots & \ddots &  \\ & & \ddots & \alpha_{n-1} & \beta_{n-1} \\ 0 & & & \beta_{n-1} & \alpha_n  \end{bmatrix} $$

比較

$$ AQ=QT $$

兩邊矩陣的每一列,得到

$$ Aq_i = \beta_{i-1}q_{i-1}+a_iq_i+\beta_iq_{i+1}, \quad i=1,2,\cdots, n$$

由於$ \beta_0, \beta_n, q_0, q_{n+1} $未定義,補充定義$ \beta_0q_0 = \beta_nq_{n+1} = 0 $,這樣上式兩邊乘$q_i^T$得到

$$ \alpha_i = q_i^TAq_i, \qquad \beta_i = q_{i+1}^TAq_i=||Aq_i-\beta_{i-1}q_{i-1}-\alpha_iq_i||_2$$

可以看出只要給定任意的$q_1 \in \mathbb{R}^n$且$||q_1||_2=1$,就能夠利用遞推關系得到全部$ q_i, \alpha_i, \beta_i $,迭代格式為

$ \alpha_1 = q_1^TAq_1, $       //初始值

$ r_i=Aq_i-\alpha_iq_i-\beta_{i-1}q_{i-1}, $

$ \beta_i=||r_i||_2, $         //由上一輪$\alpha$,$q$產生新的值

$ q_{i+1}=r_i/\beta_1 \ (\beta_i \neq 0)$

$ \alpha_{i+1} = q_{i+1}^TAq_{i+1}, \ i=1,2,\cdots, n-1$

此即Lanczos迭代。其中$q_i$稱為Lanczos向量。在第$j$步產生的矩陣$T_j$稱為j階Lanczos矩陣,其特征值可能是$A$的部分特征值的很好的近似。詳細內容參考Krylov子空間。

注意其中若某一步$\beta_i=0$,則此時得到的$T_j$的特征值將都是$A$的特征值。

 

二、具體算法

Lanczos算法求大型稀疏矩陣$A$特征值(三對角矩陣)

1. 輸入$A \in S\mathbb{R}^{n\times n}, q_1 \in \mathbb{R}^n\ (||q_1||_2=1)$

2. $u_1:=Aq_1,\ j:=1$

3. $a_j:=q_j^Tu_j,\ r_j:=u_j-a_jq_j,\ \beta_j:=||r_j||_2$

4. 若$\beta_j=0$則結束,否則

    $q_{j+1}:=r_j/\beta_j,\ u_{j+1}:=Aq_{j+1}-\beta_jq_j$

    $j:=j+1$,轉步3

這一算法的主要工作量集中在計算矩陣$A$與向量$v$的乘積$Av$上。在實際使用時,應根據$A$的具體特點,設計一個計算$Av$的子程序,使算法的運算量盡可能少。

 

三、Lanczos現象

如果Lanczos算法是在沒有攝入誤差的情況下執行的,則所得到的Lanczos向量$q_1,\cdots,q_j$是相互正交的,而且之多$n$步必然終止。但是,在誤差出現的情況下,計算得到的Lanczos向量的正交性很快就會失去,有時甚至還是線性相關的。C.C.Paige指出失去正交性恰與近似特征值的精度提高有關。

在Lanczos矩陣$T_j$中,其特征值$\mu_j$只要與其他$\mu_i$不要太靠近,特征向量$||z_j||_2$就和1很接近。而迭代次數$k$越大(可以超過$n$),則$T_k$含有越多的$A$的較好地近似特征值。當$k$充分大時,$T_k$就含有$A$的所有相異的特征值,此即Lanczos現象。

粗略地講,位於$A$譜區間兩端且分離較好地特征值$\lambda$,在$k\ll n$時$T_k$的特征值內就含有$\lambda$的很好的近似值;位於區間內部而又與其他特征值分離得不好的特征值$\lambda$,需$k\gg n$,$T_k$的特征值中才會有$\lambda$的較好地近似值。

因此,當我們只需計算大型稀疏對稱矩陣$A$的少數幾個兩端特征值時,通常只需迭代很少幾步($k\ll n$),$T_k$的兩端特征值就是$A$兩端特征值的很好的近似值。

一個直觀的例子是Parllet的數據,當$n=10^4$時,取$k=300$就可以求出10個$A$的兩端特征值和對應的特征向量的很好的近似值。

若求$A$全部的特征值時,一般$k$要遠大於$n$。Cullum和Willoughby的經驗是對於絕大多數矩陣來講,只需$k\leqslant 3n$就可以求出其幾乎全部的不同特征值達到機器精度的近似值。

顯然,對於$\textbf{HC}=E\textbf{C}$問題的精確基態和幾個低能態來說,Lanczos方法是有用而且可以進一步發展的。在這個問題上最高能級一般不要求求出,那么得到低能級的精確值的方法還可以進一步加速。具體的方法就是計算化學中常用的Davidson對角化。


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2024 CODEPRJ.COM