線性方程組(8)-共軛殘差法

03-03

到了廣義最小殘差法。和蘭佐斯演算法一樣，最小殘差法利用了矩陣的對稱性，從而將遞推關係簡化到只與前兩步有關。而共軛梯度法則是蘭佐斯演算法在變分原理的觀點下推導出的等價演算法，因此我們也可以嘗試用變分原理的觀點去推導GMRES的等價演算法，即考慮函數

$Psileft(vec{x} ight)=lVert vec{r} Vert^2=lVert vec{b}-Avec{x} Vert^2=Avec{x}cdot Avec{x}-2vec{b}cdot Avec{x} + vec{b}cdotvec{b}$

其梯度

$abla Psileft(vec{x} ight)=2A^T Avec{x}-2A^T vec{b}=2A^T vec{r}$

一維優化問題

若前一步的解為 $vec{x}^{left(k-1 ight)}$ ，沿 $vec{p_k}$ 方向搜索 $Psileft(vec{x} ight)$ 的最小值

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

代入得

$Psileft(vec{x}^{left(k ight)} ight)=Avec{p_k}cdot Avec{p_k}alpha_k^2-2vec{r}^{left(k-1 ight)}cdot Avec{p_k}alpha_k+Psileft(vec{x}^{left(k-1 ight)} ight)$

可知最佳的 $alpha_k$ 為

$alpha_k=frac{vec{r}^{left(k-1 ight)}cdot Avec{p_k}}{Avec{p_k}cdot Avec{p_k}}$

廣義共軛殘差法（Generalized conjugate residual method, GCR）

假設前 $k-1$ 步使用的搜索方向是線性無關的向量 $vec{p_1},vec{p_2},dots,vec{p_{k-1}}$ ，而且前 $k-1$ 步得到的 $vec{x}^{left(k-1 ight)}$ 滿足

$Psileft(vec{x}^{left(k-1 ight)} ight)=min_{vec{x}in spanleft{vec{p_1},vec{p_2},dots,vec{p_{k-1}} ight}} Psileft(vec{x} ight)$

且希望第 $k$ 步得到的結果

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

滿足

$Psileft(vec{x}^{left(k ight)} ight)=min_{vec{x}in spanleft{vec{p_1},vec{p_2},dots,vec{p_{k}} ight}} Psileft(vec{x} ight)$

將 $vec{x}^{left(k-1 ight)}$ 和 $vec{x}^{left(k-1 ight)}$ 在對應子空間的基上展開，並對展開係數求偏導可知，對於 $i=1,2,cdots,k-1$

$Avec{p_k}cdot Avec{p_i}=0$

設

$vec{p_k}=vec{r}^{left(k-1 ight)}-sum_{i=1}^{k-1}eta_{i,k} vec{p_i}$

則有

$Avec{p_k}=Avec{r}^{left(k-1 ight)}+sum_{i=1}^{k-1}eta_{i,k} Avec{p_i}$

也就是將 $Avec{r}^{left(k-1 ight)}$ 對 $Avec{p_i}$ 施密特正交化

$eta_{i,k}=-frac{Avec{r}^{left(k-1 ight)}cdot Avec{p_i}}{Avec{p_i}cdot Avec{p_i}}$

於是就得到了廣義共軛殘差法。

當 $K_1,dots,K_{k-1}$ 不是 $A$ 的特徵子空間時（如果是，那麼在第 $k-1$ 步可求得精確解），可以證明

$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}=K_k$

當 $k=1$ 和 $k=2$ 時容易驗證成立。假設當 $k=l-1$ 和 $k=l$ 時成立，那麼

$vec{r}^{left(l ight)}=vec{r}^{left(l-1 ight)}+alpha_l Avec{p_l}in spanleft{K_l,AK_l ight}=K_{l+1}$

即

$spanleft{vec{r}^{left(0 ight)},vec{r}^{left(1 ight)},dots,vec{r}^{left(l ight)} ight}subset K_{l+1}$

又因為 $vec{p_l}in K_{l}/K_{l-1}$ ，所以 $vec{p_l} otin K_{l-1}$ ；又因為 $K_l$ 不是 $A$ 的特徵子空間，所以 $Avec{p_l} otin K_l$ ，即

$vec{r}^{left(l ight)} otin K_l$

所以 $spanleft{vec{r}^{left(0 ight)},vec{r}^{left(1 ight)},dots,vec{r}^{left(l ight)} ight}$ 的維度為 $l+1$ ，所以

$spanleft{vec{r}^{left(0 ight)},vec{r}^{left(1 ight)},dots,vec{r}^{left(l ight)} ight}= K_{l+1}$

由

$vec{p_{l+1}}=vec{r}^{left(l ight)}-sum_{i=1}^{l}eta_{i,l+1} vec{p_i}$

可得

$vec{p_{l+1}}in K_{l+1}$

$vec{p_{l+1}} otin K_l$

同理有

$spanleft{vec{p_1},vec{p_2},dots,vec{p_{l+1}} ight}=K_{l+1}$

即當 $k=l+1$ 時也成立。歸納可知對任意的 $k$ 都成立。

由

$Psileft(vec{x}^{left(k ight)} ight)=min_{vec{x}in spanleft{vec{p_1},vec{p_2},dots,vec{p_{k}} ight}} Psileft(vec{x} ight)$

可知

$abla Psileft(vec{x}^{left(k ight)} ight)=2A^T vec{r}^{left(k ight)}ot spanleft{vec{p_1},vec{p_2},dots,vec{p_{k}} ight}$

即

$vec{r}^{left(k ight)}ot AK_k$

對 $i=1,2,dots,k-1$ 有

$vec{r}^{left(k ight)}ot Avec{r}^{left(i ight)}$

每一步的殘差與前面所有的殘差 $A$ 共軛，這就是「共軛殘差」名字的來源。注意當 $A$ 非對稱時，共軛關係也是非對稱的

$vec{r}^{left(i ight)}ot Avec{r}^{left(j ight)} otLeftrightarrow Avec{r}^{left(i ight)}ot vec{r}^{left(j ight)}$

$alpha_k$ 的計算式可以進一步簡化

$alpha_k=frac{vec{r}^{left(k-1 ight)}cdot Avec{p_k}}{Avec{p_k}cdot Avec{p_k}}=frac{vec{r}^{left(k-1 ight)}cdot Avec{r}^{left(k-1 ight)}-sum_{i=1}^{k-1}eta_{i,k} vec{r}^{left(k-1 ight)}cdot Avec{p_i}}{Avec{p_k}cdot Avec{p_k}}=frac{vec{r}^{left(k-1 ight)}cdot Avec{r}^{left(k-1 ight)}}{Avec{p_k}cdot Avec{p_k}}$

注意到與共軛梯度法不同的是，GCR中的 $A$ 可以是任意可逆矩陣，且運算過程中仍然需要保存 $vec{p_1},dots,vec{p_{k-1}}$ 。為了限制其增長，也有重啟和截斷兩種方法。

用MATLAB代碼描述截斷的GCR

function [x, iter] = gcr_solve(A, b, x0, m, epsilon, iter_max) iter = 0; n = size(A, 1); p = zeros(n, m); Ap = zeros(n, m); Ap_square = zeros(m, 1); x = x0; r = b - A * x0; Ar = A * r; p(:, 1) = r; Ap(:, 1) = Ar; Ap_square(:, 1) = dot(Ar, Ar); while norm(r) >= epsilon && iter < iter_max iter = iter + 1; i_curr = mod(iter - 1, m) + 1; i_next = mod(iter, m) + 1; alpha = dot(r, Ar) / Ap_square(i_curr); x = x + alpha * p(:, i_curr); r = r - alpha * Ap(:, i_curr); Ar = A * r; p_next = r; Ap_next = Ar; for prev = max(iter - m + 1, 1): iter i_prev = mod(prev - 1, m) + 1; beta = - dot(Ap_next, Ap(:, i_prev)) / Ap_square(i_prev); p_next = p_next + beta * p(:, i_prev); Ap_next = Ap_next + beta * Ap(:, i_prev); end p(:, i_next) = p_next; Ap(:, i_next) = Ap_next; Ap_square(i_next) = dot(Ap_next, Ap_next); endend

共軛殘差法（conjugate residual method, CR）

共軛殘差法可以視為廣義共軛殘差法在 $A$ 對稱時的特殊情況。

因為

$Avec{p_i}in AK_{i}subset K_{i+1}$

當 $i<k-1$ 時，對於對稱的 $A$ ，有

$Avec{r}^{left(k-1 ight)}cdot Avec{p_i}=vec{r}^{left(k-1 ight)}cdot Aleft(Avec{p_i} ight)=0$

於是遞推關係可簡化為

$vec{p_k}=vec{r}^{left(k-1 ight)}-eta_{k-1,k} vec{p_{k-1}}$

其中

$eta_{k-1,k}=-frac{Avec{r}^{left(k-1 ight)}cdot Avec{p_{k-1}}}{Avec{p_{k-1}}cdot Avec{p_{k-1}}}=-frac{Avec{r}^{left(k-1 ight)}cdot left(vec{r}^{left(k-2 ight)}-vec{r}^{left(k-1 ight)} ight)}{alpha_{k-1}Avec{p_{k-1}}cdot Avec{p_{k-1}}}=frac{vec{r}^{left(k-1 ight)}cdot Avec{r}^{left(k-1 ight)}}{vec{r}^{left(k-2 ight)}cdot Avec{r}^{left(k-2 ight)}}$

用MATLAB代碼描述

function [x, iter] = cr_solve(A, b, x0, epsilon, iter_max) iter = 0; x = x0; r = b - A * x0; p = r; Ap = A * p; rAr_prev = dot(r, A * r); while norm(r) >= epsilon && iter < iter_max iter = iter + 1; alpha = rAr_prev / dot(Ap, Ap); x = x + alpha * p; r = r - alpha * Ap; Ar = A * r; rAr_curr = dot(r, Ar); beta = rAr_curr / rAr_prev; p = r + beta * p; Ap = Ar + beta * Ap; rAr_prev = rAr_curr; endend

GCR與GMRES的等價性

GCR與GMRES都是在 $vec{x}^{left(k ight)}in K_k$ 中最小化 $lVert vec{r}^{left(k ight)} Vert$ ，所以理論上應是等價的。嚴格證明留作練習。

迭代停滯

對於一般的 $A$ ，在GCR中可能出現 $vec{r}^{left(k-1 ight)}cdot Avec{r}^{left(k-1 ight)}=0$ ，則 $alpha_k=0$ ，導致迭代停滯。

在GMRES中同樣會出現迭代停滯，即

$vec{h_k^prime}=Q_{k-1}^Tvec{h_k}$

的最後一個元素為零時，會出現

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