lanczos演算法及C 實現(一)框架及簡單實現

2022-09-03 18:18:07 字數 4187 閱讀 7857

其中$u$是列正交矩陣,即 $u^tu=i$,每一列為乙個特徵向量,$s$是對角陣,對角線上每個元素為特徵值。$r$為分解的秩

lanczos演算法分三步求解:

1) 對$x$進行正交變換得到乙個三對角陣$t$:$x=ptp^t$,其中$p$為正交陣

$t= \left[

\begin \alpha_0 & \beta_1 & 0 & 0 & 0& 0& \dots \\

\beta_1 & \alpha_1 & \beta_2 & 0 & 0& 0& \dots \\

0 & \beta_2 &\alpha_3 &\beta_3 & 0& 0 & \dots \\

0 & 0 & \beta_3 &\alpha_4 &\beta_4 & 0& \dots \\

0 & 0 & 0 & \beta_4 &\alpha_5 &\beta_5 & \dots \\

\dots & \dots &\dots &\dots &\dots &\dots &\dots\\

\end

\right]$

2) 對$t$進行奇異值分解:$t=w_sw^t_$

1) $t$矩陣是乙個三對角陣,很稀疏,因此對它的特徵值分解會非常快。

2) 如果$r$很小,可以不用求整個t,而是求其左上$1.5r$階的方陣即可得到很好的近似。這樣對t的特徵值分解會更快。( 另外一種方法動態決定求多少lanczos向量, 這裡簡單總結一下,初始選取$n \le m$個lanczos向量,然後求t的最大特徵值和最小特徵值,然後逐漸增加$n$,並更新這兩個特徵值,直到這種更新非常小為止。如果用在求$r$個主特徵向量,那麼可以更新第$1$個和第$r$個特徵值。

考慮冪迭代產生的一系列向量$p = [v, xv, x^2v, \dots x^v]$,其中$v$是任意向量。假設$x$的特徵值分解為$x = \sum_i \lambda_i \xi_i\xi^t_i$,並且特徵值按絕對值從大到小排序 $|\lambda_1| \ge |\lambda_2| \ge \dots$。則 $xv = \sum_i \lambda_i \xi_i \xi_i^tv = \sum_i \lambda_i v'_i \xi_i$,其中$v'_i = \xi_i^tv$是$v$在正交基$[\xi_1, \xi_2\dots ]$下的座標。 類似的,我們可以得到$x^kv = \sum_i \lambda^k_i v'_i \xi_i$。也就是說特徵值大的特徵向量在$x^kv$中比重越大(這也是冪迭代收斂到主特徵值的原因) 。換句話說,$p$中的每列$x^kv$都是成特徵向量$\xi_1, \xi_2, \dots$的線性組合,這些特徵向量構成了 $p$ 的一組正交基,並且特徵值越大的特徵向量所佔的比例越高。乙個自然的想法就是,對$p$的前$l$個向量所構成的矩陣進行特徵值分解,應該能得到$x$最大的幾個特徵向量的近似解,$l$越大,越近似。

然而,這種原始的方法並不穩定。因為較小的特徵值被快速的拋棄,即使是第二大特徵值,也是隨$l$增大,其比重指數速度縮小($|\lambda_1^l| >> |\lambda_2^l|$)。因此,為了能夠保留小特徵向量,需要作正交化,即$p$是列正交的。由此便產生了arnoldi's algorithm 使用施密特正交化:每次迭代所得向量都要去掉之前所有向量的投影:$p_=xp_k \ominus p_k \ominus p_ \ominus p_ \dots$

其中$a \ominus b$ 表示 a在b的正交空間的投影:

正交化以後,雖然$x$特徵向量仍舊是$p$的一組正交基,但大特徵向量在$p$中的比重就不一定高了。因此$p$的最大特徵向量並不是$x$特徵向量的近似。

注意到如果對$p$每列歸一化: $p_k = \frac$,則$p$就是乙個正交陣,因此$p^txp$和$x$有著相同的特徵值,這意味著如果能求得$p^txp$的主特徵向量$w_1, w_2 \dots$,則$pw_1, pw_2, \dots$就是$x$的主特徵向量。很重要的是$p^txp$是三對角陣! 證明可參見後文。而求解三對角陣的特徵向量會容易很多!

不過arnoldi's algorithm使用的是施密特正交化,要求$o(l^2)$次向量內積,相減。隨著$l$增大,該方法速度堪憂。

lanczos提出了更快的$o(l)$複雜度的正交化方法。由於$x$是實對稱陣,可以證明,每次迭代得到的向量只需要和前兩次得到的向量正交後,即可保證該向量和其他所有向量正交。在列歸一化後,滿足$x=ptp^t$, $t$為三對角陣。

$p_=\frac}\|}$

**  解釋了該方法實際上就是用迭代後的向量和前兩次向量做正交,但是並沒有證明這樣得到的向量和其他向量一定正交和為什麼得到的是三對角陣,本文在附錄a中,給出了lanczos方法的正交以及三對角性的證明。

三對角化後的矩陣t有三種常用的方法進行分解:

1) qr法,即對t進行qr分解,$t=qr$,其中$q$為正交陣,$r$為上三角陣,然後令$t'=rq$,如此反覆進行,直至收斂。其中產生的$q$矩陣的連乘就得到了$t$的特徵向量。具體演算法及實現請參考

文章第二頁給出了偽**。在2.2.2節中,還提到了針對上海森伯格陣的快速qr演算法,正好適用於這裡的$t$矩陣。這裡不再贅述。相比於另兩種方法,qr方法並不高效,但是在求最大的少數特徵向量時,$t$通常很小,因此為實現簡單,以及穩定性考慮,qr方法常被使用。

3) mrrr方法。這個方法是理論上最快的方法,$o(l^2)$複雜度,主要採用逆迭代來加快收斂。詳細討論見後續文章。該方法作者的博士** a new o(n2) algorithm for the symmetric tridiagonal eigenvalue/eigenvector problem

這裡討論大型、稀疏、非方陣a的奇異值分解的實現,並求其少數主奇異向量。主要從如下角度考慮:

1)對於非方陣的svd可以轉化為對稱方陣,本文採用了通過對$x=a*a^t$的特徵值分解從而得到$a$的奇異值分解的方法。$x$的特徵向量即為$a$的左奇異向量,而$a$的右奇異向量可以通過左奇異向量左乘$a$得到。

2)由於矩陣非常大,不能完全讀入記憶體,因此是存在硬碟上的。在計算lanczos向量時,每次計算$a*a^t*p$需要掃兩遍磁碟。$a$在磁碟上順序儲存,順序讀取,這樣可以利用磁碟的讀快取。

3)考慮到多執行緒平行計算$a*a^t*p$,需要開闢讀快取,每次裝載一批資料讀入記憶體。每個執行緒維護自己的乙個$p$向量,迭代完之後合併。

4)由於只求少數主奇異向量,因此lanczos向量也不會很多,因此可以假設lanczos向量可以完全放在記憶體中。

本文中求lanczos向量採用的是

附錄a:由lanczos方法得到的$p$矩陣是正交陣,$p^txp$是三對角陣。

首先引入乙個符號,對於兩個向量$a,b$, $a \oplus b=\lambda a + \theta b$,表示由$a$和$b$張成的空間中的某乙個向量

證明:採用數學歸納法。需要證明兩個命題

(1)正交性: $\forall i

(2) 三對角:$\forall i

先證明正交性。當$k\le2$時,結論顯然成立。當$k>2$時,假設結論成立,即對任意的$i < j \le k, p_i^tp_j=0$,現在要證明對所有$i \le k, p_i^tp_=0$

分兩種情況,如果$i < k - 1$,則有

$\begin\nonumber & p_^tp_i \\\nonumber =& (xp_k \ominus p_k \ominus p_)^tp_i \\\nonumber =& (xp_k)^tp_i \\\nonumber =& p_k^txp_i \\\nonumber =& p_k^t(p_\oplus p_i \oplus p_)\\\nonumber =& p_k^tp_ \\\nonumber = & 0\end$

否則如果$i = k - 1$或者$i = k$,根據$p_ =xp_k \ominus p_ \ominus p_$,直接可以得到$p_i^tp_=0$

再證明三對角性。當$k\le 2$時,結論顯然成立,當$k> 2$時,假設結論成立,則當$k+1$時,對於$i+1 < k+1$的任意$i$有

$ p_^txp_i = p_(p_\oplus p_i \oplus p_) =0$

最後乙個等號由正交性得出。證畢

目錄:lanczos演算法及c++實現(〇)c++**

lanczos演算法及c++實現(一)框架及簡單實現

lanczos演算法及c++實現(二)實對稱陣奇異值分解的qr演算法

lanczos演算法及c++實現(三)實對稱三對角陣特徵值分解的分治演算法

排序演算法分析及實現C

參考資料 經典排序演算法過程 經典排序演算法性質 王道資料結構 嚴蔚敏 資料結構 穩定 如果a原本在b前面,而a b,排序之後a仍然在b的前面。不穩定 如果a原本在b的前面,而a b,排序之後 a 可能會出現在 b 的後面。時間複雜度 對排序資料的總的操作次數。反映當n變化時,操作次數呈現什麼規律。...

插入演算法演算法詳解及實現 c語言

插入排序原理很簡單,將一組資料分成兩組,分別將其稱為有序組與待插入組。每次從待插入組中取出乙個元素,與有序組的元素進行比較,並找到合適的位置,將該元素插到有序組當中。就這樣,每次插入乙個元素,有序組增加,待插入組減少。直到待插入組元素個數為0。實現 include include include 插...

CRC16 演算法及c實現

標準crc生成多項式如下表 名稱 生成多項式 簡記式 標準引用 crc 4 x4 x 1 3 itu g.704 crc 8 x8 x5 x4 1 0x31 crc 8 x8 x2 x1 1 0x07 crc 8 x8 x6 x4 x3 x2 x1 0x5e crc 12 x12 x11 x3 x ...