線性方程組(6)-阿諾爾迪演算法

02-28

上一節講到了克雷洛夫子空間單位正交基的計算方法，和伽遼金原理。根據伽遼金原理，把共軛梯度法推廣到一般情況。對一般的矩陣 $A$ ，使得第 $k$ 步的解

$vec{x}^{left(k ight)}in K_k$

其殘差

$vec{r}^{left(k ight)}=vec{b}-Avec{x}^{left(k ight)}ot K_k$

即 $K=L=K_k$ ，這就是阿諾爾迪演算法（Arnoldi algorithm）。

根據伽遼金原理，要解得方程是

$V_k^TAV_k vec{y_k}=V_k^Tb$

其中 $V_k$ 得列是 $K_k$ 的一組單位正交基，且 $vec{v_1}=frac{vec{b}}{lVert vec{b} Vert}$ 。則有

$H_kvec{y_k}=lVert vec{b} Vertvec{e_1}$

若 $H_{k-1}$ 和 $H_k$ 有LU分解

$H_{k-1}=L_{k-1}U_{k-1}$

$H_{k}=L_{k}U_{k}$

由於 $H$ 是上海森堡陣，所以這裡的 $L$ 只有主對角線和下次對角線非零。

$H_k=left(egin{matrix} H_{k-1} & vec{h_k}\ h_{k,k-1}vec{e_{k-1}}^T & h_{k,k} end{matrix} ight)$

$L_k=left(egin{matrix} L_{k-1} & 0\ l_{k}vec{e_{k-1}}^T & 1 end{matrix} ight)$

$U_k=left(egin{matrix} U_{k-1} & vec{u_k}\ 0 & u_{k,k} end{matrix} ight)$

矩陣分塊乘法可得

$L_{k-1}vec{u_k}=vec{h_k}$

$l_k u_{k-1,k}+u_{k,k}=h_{k,k}$

$l_k u_{k-1,k-1}=h_{k,k-1}$

則有

$vec{u_k}=L_{k-1}^{-1}vec{h_k}$

$u_{k,k}=h_{k,k}-l_k u_{k-1,k}$

$l_k =frac{h_{k,k-1}}{u_{k-1,k-1}}$

最終得到的解

$vec{x}^{left(k ight)}=V_k vec{y_k}=V_k U_k^{-1} L_k^{-1}lVertvec{b} Vertvec{e_1}$

令

$V_k U_k^{-1}= P_k$

即

$V_k = P_k U_k$

取最後一列，得到

$vec{p_k}=frac{1}{u_{k,k}}left(vec{v_k}-P_{k-1}vec{u_k} ight)$

令

$L_k^{-1}lVertvec{b} Vertvec{e_1}=vec{xi_k}$

即

$lVertvec{b} Vertvec{e_1}=L_kvec{xi_k}$

取最後一行，得到

$xi_k=-l_kxi_{k-1}$

最終得到

$vec{x}^{left(k ight)}=P_k vec{xi_k}=P_{k-1} vec{xi_{k-1}}+xi_{k}vec{p_k}=vec{x}^{left(k-1 ight)}+xi_{k}vec{p_k}$

初值

$vec{p_1}=frac{vec{v_1}}{u_{1,1}}=frac{vec{v_1}}{h_{1,1}}$

$xi_1=lVert vec{b} Vert$

由阿諾爾迪過程的矩陣形式

$AV_{k}=V_{k+1} ar{H_{k}}=V_{k} H_{k}+h_{k+1,k}vec{v_{k+1}}vec{e_{k}}^T$

右乘向量 $vec{y_k}$ 得到

$Avec{x}^{left(k ight)}=vec{b}+h_{k+1,k}y_{k}vec{v_{k+1}}$

而

$U_k vec{y_k}=vec{xi_k}$

即

$y_k=frac{xi_k}{u_{k,k}}$

最後得到

$vec{r}^{left(k ight)}=-y_{k}vec{v_{k+1}}=-frac{ h_{k+1,k}xi_k}{u_{k,k}}vec{v_{k+1}}$

$lVertvec{r}^{left(k ight)} Vert= vertfrac{ h_{k+1,k}xi_k}{u_{k,k}} vert$

與共軛梯度法一樣，由於 $vec{x}$ 只涉及加法， $vec{x}^{left(0 ight)}=vec{x_0}$ 非零不影響結果，只要用 $vec{r_0}=vec{b}-Avec{x_0}$ 替換 $vec{b}$ 即可。

有限存儲空間的阿諾爾迪演算法

阿諾爾迪演算法需要保存前面所有的向量 $vec{v_1},dots,vec{v_{k-1}}$ ，隨著迭代步數的增加，消耗的存儲空間不斷增加，每步的計算量也不斷增加。為了限制其增長，阿諾爾迪演算法有重啟法和截斷法兩種變種。

重啟法（restarted method）就是當進行到第 $m$ 步時，將得到的 $vec{x}^{left(m ight)}$ 作為初值重新開始阿諾爾迪演算法。重啟法可以不用每一步都求出解 $vec{x}^{left(k ight)}$ ，而是直接算出 $V_m$ 和 $H_m$ 並求解得到 $vec{x}^{left(m ight)}$ ，而且這裡可以用豪斯霍爾德變換法求 $V_m$ 和 $H_m$ 。

截斷法（truncated method）又稱為不完全正交法（incomplete orthogonal method）。忽略 $vec{v_k}$ 與 $vec{v_1},dots,vec{v_{k-m-1}}$ 之間的正交性要求，使得當 $k>m$ 時

$h_{1,k}=cdots=h_{k-m-1,k}=0$

同時

$u_{1,k}=cdots=u_{k-m-1,k}=0$

若

$H_k=left(egin{matrix} *&* & 0\ *&H_{k-1}^prime & vec{h_k^prime}\ 0 & h_{k,k-1}vec{e_{m}}^T & h_{k,k} end{matrix} ight)$

$L_k=left(egin{matrix} * & 0 & 0 \ * & L_{k-1}^prime & 0\ 0 & l_{k}vec{e_{m}}^T & 1 end{matrix} ight)$

$U_k=left(egin{matrix} * & * & 0 \ 0 & U_{k-1}^prime & vec{u_k^prime}\ 0& 0 & u_{k,k} end{matrix} ight)$

則

$vec{u_k^prime}=L_{k-1}^{prime -1}vec{h_k^prime}$

$u_{k,k}=h_{k,k}-l_k u_{k-1,k}$

$l_k =frac{h_{k,k-1}}{u_{k-1,k-1}}$

以及前面已經推導的

$h_{k,k-1} vec{v_k} = A vec{v_{k-1}}-sum_{i=k-m}^{k-1}h_{i,k}vec{v_{i}}$

$vec{p_k}=frac{1}{u_{k,k}}left(vec{v_k}-sum_{i=k-m}^{k-1}u_{i,k}vec{p_{i}} ight)$

$xi_k=-l_kxi_{k-1}$

$vec{x}^{left(k ight)}=vec{x}^{left(k-1 ight)}+xi_{k}vec{p_k}$

$lVertvec{r}^{left(k ight)} Vert= vertfrac{ h_{k+1,k}xi_k}{u_{k,k}} vert$

在實現中需要用到循環隊列（circular queue），有興趣的讀者可以嘗試。

蘭佐斯演算法（Lanczos algorithm）

當 $A$ 對稱時， $H_k$ 變為對稱三對角矩陣 $T_k$ ，於是演算法可以得到簡化。有

$u_{1,k}=cdots=u_{k-2,k}=0$

令

$alpha_i = t_{i,i}$

$eta_i=t_{i-1,i}=t_{i,i-1}$

則

$T_k=left(egin{matrix} T_{k-1} & eta_kvec{e_{k-1}}\ eta_kvec{e_{k-1}}^T & alpha_k end{matrix} ight)$

$L_k=left(egin{matrix} L_{k-1} & 0\ l_{k}vec{e_{k-1}}^T & 1 end{matrix} ight)$

$U_k=left(egin{matrix} U_{k-1} & u_{k-1,k}vec{e_{k-1}}\ 0 & u_{k,k} end{matrix} ight)$

分塊乘法可得

$u_{k-1,k}=eta_k$

$l_k u_{k-1,k-1}=eta_k$

$l_k u_{k-1,k}+u_{k,k}=alpha_k$

即

$l_k =frac{eta_k}{u_{k-1,k-1}}$

$u_{k,k}=alpha_k-l_k eta_k$

以及前面已經推導的

$eta_k vec{v_k} = A vec{v_{k-1}}-alpha_{k-1}vec{v_{k-1}}-eta_{k-1}vec{v_{k-2}}$

$vec{p_k}=frac{1}{u_{k,k}}left(vec{v_k}-eta_kvec{p_{k-1}} ight)$

$xi_k=-l_kxi_{k-1}$

$vec{x}^{left(k ight)}=vec{x}^{left(k-1 ight)}+xi_{k}vec{p_k}$

$Vertvec{r}^{left(k ight)} Vert=lvertfrac{ eta_{k+1}xi_k}{u_{k,k}} vert$

用MATLAB代碼描述（寫代碼的時候發現l不好看，把 $l$ 換成了 $lambda$ ）

function [x, iter] = lanczos_solve(A, b, x0, epsilon, iter_max) iter = 0; n = size(A, 1); x = x0; w = b - A * x0; r_norm = norm(w); v = w / r_norm; beta = 0; p = zeros(n, 1); while r_norm >= epsilon && iter < iter_max iter = iter + 1; if iter == 1 w = A * v; lambda = 0; xi = r_norm; else w = A * v - beta * v_prev; lambda = beta / u; % u == u_{k-1,k-1} xi = - lambda * xi; end v_prev = v; alpha = dot(w, v); % Modified Gram-Schmit u = alpha - lambda * beta; % u := u_{k,k} p = (v - beta * p) / u; x = x + xi * p; w = w - alpha * v; beta = norm(w); % beta := beta_{k+1} v = w / beta; % v := v_{k+1} r_norm = abs(beta * xi / u); endend

每一步仍然時一次矩陣向量乘和兩次內積，但有五次向量的線性組合（算上了 $vec{v_k} =frac{ vec{w_k}}{eta_k}$ ），相比共軛梯度法多了兩次。

不考慮浮點誤差，蘭佐斯演算法與共軛梯度法得到的迭代結果是相等的。在推導蘭佐斯演算法的過程中，沒有預設 $A$ 的正定性，這也說明共軛梯度法可以不局限於正定矩陣。但是當 $A$ 不定時，蘭佐斯演算法與共軛梯度法都有數值不穩定的問題，收斂速度較慢。

下面證明蘭佐斯演算法與共軛梯度法的等價性。

首先，在蘭佐斯演算法中

$D_k=P_k^T AP_k=U_k^{-T}V_k^T AV_k U_k^{-1}=U_k^{-T}T_k U_k^{-1}=U_k^{-T}L_k$

注意到 $U_k^{-T}$ 和 $L_k$ 均為下三角陣，所以 $D_k$ 是下三角陣；又因為 $A$ 對稱，所以 $P_k^T AP_k$ 對稱。也就是說 $D_k$ 是對角陣，這就證明了對於 $i ot= j$ ，有

$vec{p_i}cdot Avec{p_j}=0$

又由

$U_k^T D_k=L_k$

取最後一行最後一列，有

$u_{k,k}vec{p_k}cdot Avec{p_k}=1$

將共軛梯度法第 $k$ 步的解、殘差和搜索方向分別記為 $vec{x}^{primeleft(k ight)}$ 、 $vec{r}^{primeleft(k ight)}$ 和 $vec{p_k^prime}$ 。

在第一步，初值相等， $vec{p_{1}}$ 與 $vec{p_{1}^prime}$ 均與 $vec{r}^{left(0 ight)}$ 共線。

假設在蘭佐斯演算法的第 $k-1$ 步得到了與共軛梯度法中相同的解，且 $vec{p_{k-1}}$ 與 $vec{p_{k-1}^prime}$ 共線。

由

$vec{p_k}in spanleft{vec{v_k},vec{p_{k-1}} ight}=spanleft{vec{r}^{left(k-1 ight)},vec{p_{k-1}} ight}=spanleft{vec{r}^{primeleft(k-1 ight)},vec{p_{k-1}^prime} ight}$

和

$vec{p_k}cdot Avec{p_{k-1}^prime}=vec{p_k}cdot Avec{p_{k-1}}=0$

可知 $vec{p_k}$ 與 $vec{p_k^prime}$ 共線。

同時

$xi_k=-l_kxi_{k-1}=-frac{eta_k xi_{k-1}}{u_{k-1,k-1}}=-frac{eta_k xi_{k-1}}{u_{k-1,k-1}u_{k,k}vec{p_k}cdot Avec{p_k}}$

$=-frac{frac{eta_k xi_{k-1}}{u_{k-1,k-1}}vec{v_k}cdot frac{1}{u_{k,k}}left(vec{v_k}-eta_kvec{p_{k-1}} ight)}{vec{p_k}cdot Avec{p_k}}=frac{vec{r}^{left(k-1 ight)}cdot vec{p_k}}{vec{p_k}cdot Avec{p_k}}$