線性方程組(5)-克雷洛夫子空間與伽遼金原理

02-28

上一節講到了變分原理和共軛梯度法。在共軛梯度法的推導中，我們提出了一條引理

$spanleft{vec{p_1},vec{p_2},dots,vec{p_{k}} ight}=spanleft{vec{r}^{left(0 ight)},vec{r}^{left(1 ight)},dots,vec{r}^{left(k-1 ight)} ight}$

$=spanleft{vec{b},Avec{b},dots,A^{k-1}vec{b} ight}$

這個子空間稱為 $k$ 階克雷洛夫子空間（Krylov subspace）

$K_k=spanleft{vec{b},Avec{b},dots,A^{k-1}vec{b} ight}$

共軛梯度法中的 $vec{r}$ 便是這個子空間的一組正交基。

對一般的矩陣，能否求克雷洛夫子空間的一組單位正交基呢？

阿諾爾迪過程（Arnoldi process）

容易發現克雷洛夫子空間的遞推性質

$K_k=spanleft{K_{k-1},AK_{k-1} ight}$

假設已經求出了 $K_{k-2}$ 的一組單位正交基 $vec{v_1},vec{v_2},dots,vec{v_{k-2}}$ 和 $K_{k-1}$ 的一組單位正交基 $vec{v_1},vec{v_2},dots,vec{v_{k-2}},vec{v_{k-1}}$ ，那麼

$K_{k}=spanleft{ K_{k-1},Avec{v_1},Avec{v_2},dots,Avec{v_{k-2}},Avec{v_{k-1}} ight}$

又因為

$spanleft{Avec{v_1},Avec{v_2},dots,Avec{v_{k-2}} ight}=AK_{k-2}subset K_{k-1}$

所以

$K_{k}=spanleft{ K_{k-1},Avec{v_{k-1}} ight}=spanleft{ vec{v_1},vec{v_2},dots,vec{v_{k-2}},vec{v_{k-1}},Avec{v_{k-1}} ight}$

施密特正交化後得到

$vec{r_k}=Avec{v_{k-1}}-sum_{i=1}^{k-1}left(vec{v_i}cdot Avec{v_{k-1}} ight)vec{v_i}$

$vec{v_k}=frac{vec{r_k}}{lVertvec{r_k} Vert}$

這就是阿諾爾迪過程。

令

$h_{i,k-1}=vec{v_i}cdot Avec{v_{k-1}}$

$h_{k,k-1}=lVertvec{r_k} Vert$

可將等式改寫為

$Avec{v_{k-1}}=sum_{i=1}^{k}h_{i,k}vec{v_i}$

將第 $1,2,dots,k,k+1$ 步的等式都寫出來，合併成矩陣形式

$AV_k=V_{k+1}ar{H_{k}}$

其中 $ar{H_{k}}$ 是一個 $left(k+1 ight) imes k$ 的矩陣。當 $i>j+1$ 時 $h_{i,j}=0$ ，也就是說， $ar{H_{k}}$ 在對角線下方只有一條次對角線非零，其餘部分相當於上三角矩陣，這樣的矩陣稱為上海森堡陣（upper Hessenburg matrix）。

阿諾爾迪過程可寫為

$h_{j,i}=vec{v_j}cdot Avec{v_{i}}$

$vec{r_{i}}=Avec{v_{i-1}}-sum_{j=1}^{i-1}h_{j,i}vec{v_j}$

$h_{i,i-1}=lVertvec{r_i} Vert$

$vec{v_i}=frac{vec{r_i}}{h_{i,i-1}}$

正交化的數值穩定性

求一組標準正交基這個運算我們非常熟悉，這就是QR分解。QR分解除了施密特正交化以外，還有豪斯霍爾德變換等方法，施密特正交化還有改進版本（modified Gram-Schmit）。由於浮點誤差，原始版本的施密特正交化得到的向量的正交性不好，在QR分解時通常還是用豪斯霍爾德變換。

阿諾爾迪過程也可以通過改進施密特正交化或豪斯霍爾德變換來實現。但是在之後由阿諾爾迪過程推導出的演算法中，豪斯霍爾德變換可能導致計算量和存儲空間消耗過大的問題，所以除非特殊說明，一般還是用改進施密特正交化。

這裡用MATLAB代碼演示原始版和改進版施密特正交化的區別

function v = gs_orthogonalize(V, u) v = u; for vk = V h = dot(u, vk); % 用最初的u v = v - h * vk; endend

function v = mgs_orthogonalize(V, u) v = u; for vk = V h = dot(v, vk); % 用最新的v v = v - h * vk; endend

蘭佐斯過程（Lanczos process）

將等式

$V_{k+1}^T AV_k=V_{k+1}^T V_{k+1}ar{H_{k}}=ar{H_{k}}$

兩邊去掉最後一行，得到

$V_{k}^T AV_k=H_{k}$

其中 $H_{k}$ 是 $k imes k$ 的上海森堡陣。

容易發現，當 $A$ 對稱時， $H_{k}$ 也對稱，這時 $H_{k}$ 是三對角矩陣（tridiagonal matrix），我們為了強調這一點，將它記作 $T_k$ 。

這是阿諾爾迪過程的特殊情況，稱為蘭佐斯過程。

$t_{i,i}=vec{v_i}cdot Avec{v_{i}}$

$vec{r_{i}}=Avec{v_{i-1}}-t_{i-1,i-1}vec{v_{i-1}}-t_{i-1,i-2}vec{v_{i-2}}$

$t_{i,i-1}=lVertvec{r_i} Vert$

$vec{v_i}=frac{vec{r_i}}{t_{i,i-1}}$

可以發現，在阿諾爾迪過程中，求一個新的 $vec{v_k}$ 需要前面所有的 $vec{v_1},dots,vec{v_{k-1}}$ ；而在蘭佐斯過程中，求一個新的 $vec{v_k}$ 只需要 $vec{v_{k-2}}$ 和 $vec{v_{k-1}}$ 。正因為這個特點，用共軛梯度法解對稱矩陣的時候，不需要存儲 $vec{p_1},dots,vec{p_{k-1}}$ 。

伽遼金原理（Galerkin principle）

在共軛梯度法中，第 $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$ 維子空間 $K=spanleft{vec{v_1},dots,vec{v_k} ight}$ 內找到 $vec{x}^{left(k ight)}$ 使得殘差正交於某個 $k$ 維子空間 $L=spanleft{vec{w_1},dots,vec{w_k} ight}$ ，作為 $Avec{x}=vec{b}$ 的近似解。假設

$vec{x}^{left(k ight)}=V_{k} vec{y}$

則方程可寫成

$W_{k}^T left(vec{b}-AV_{k} vec{y} ight)=0$

即

$W_{k}^T AV_{k} vec{y} = W_{k}^T vec{b}$

這是一個在 $mathbb{R}^k$ 空間的線性方程組，在 $k$ 較小時，其求解比 $mathbb{R}^n$ 空間的線性方程組要容易得多。