線性方程組(7)-最小殘差法

02-28

在上一節中，從共軛梯度法得到啟發，取 $K=L=K_k$ ，我們得到了阿諾爾迪演算法。通過證明蘭佐斯演算法與共軛梯度法的等價性，找到了變分原理與伽遼金原理之間的某種聯繫。

我們重新考慮變分原理並選取函數

$Psileft(vec{x} ight)=lVertvec{r} Vert^2=lVertvec{b}-Avec{x} Vert^2$

廣義最小殘差法（Generalized minimal residual method, GMRES）

若

$vec{x}^{left(k ight)}=V_kvec{y_k}$

由阿諾爾迪過程

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

則有

$lVert vec{r}^{left(k ight)} Vert=lVert vec{b}-AV_kvec{y_k} Vert=lVert vec{b}-V_{k+1}ar{H_k}vec{y_k} Vert=lVert V_{k+1}^T vec{b}-ar{H_k}vec{y_k} Vert$

其中

$V_{k+1}^T vec{b}=lVertvec{b} Vertvec{e_1}$

所以，求 $lVert vec{r}^{left(k ight)} Vert$ 的最小值等價於用最小二乘法求解超定方程

$ar{H_k}vec{y_k}=lVertvec{b} Vertvec{e_1}$

設 $ar{H_{k-1}}$ 和 $ar{H_{k}}$ 有QR分解

$ar{H_{k-1}}=Q_{k-1}ar{R_{k-1}}$

$ar{H_k}=Q_kar{R_k}$

若

$ar{H_k}=left(egin{matrix} ar{H_{k-1}} & vec{h_k}\ 0 & h_{k+1,k} end{matrix} ight)$

$ar{R_k}=left(egin{matrix} ar{R_{k-1}} & vec{r_k}\ 0& 0end{matrix} ight)$

矩陣分塊乘法可得

$Q_{k}left(egin{matrix} vec{r_k}\ 0end{matrix} ight)=left(egin{matrix} vec{h_k}\ h_{k+1,k} end{matrix} ight)$

令

$vec{h_k^prime}=Q_{k-1}^T vec{h_k}=left(egin{matrix} h_{k,1}^prime\ vdots\ h_{k,k}^prime end{matrix} ight)$

再進行一次吉文斯旋轉（Givens rotation）即可完成QR分解

$c_k=frac{h_{k,k}^prime}{sqrt{h_{k,k}^{^prime 2}+h_{k+1,k}^2}}$

$s_k=frac{h_{k+1,k}}{sqrt{h_{k,k}^{prime 2}+h_{k+1,k}^2}}$

$G_k=left(egin{matrix}I &0&0\ 0&c_k&s_k\ 0&-s_k&c_k end{matrix} ight)$

則有

$left(egin{matrix} vec{r_k}\ 0end{matrix} ight)=G_kleft(egin{matrix} vec{h_k^prime}\ h_{k+1,k} end{matrix} ight)$

$Q_k=left(egin{matrix} Q_{k-1}&0\ 0&1 end{matrix} ight)G_k^T$

方程變為

$ar{R_k}vec{y_k}=Q_k^TlVertvec{b} Vertvec{e_1}$

最小二乘解相當於解前 $k$ 行。

令 $ar{R_k}$ 的前 $k$ 行為 $R_k$ ，和阿諾爾迪演算法一樣，可以令

$P_k=V_k R_k^{-1}$

$vec{xi_k}=Q_k^TlVertvec{b} Vertvec{e_1}=left(egin{matrix} vec{xi_k^prime}\ eta_{k} end{matrix} ight)=left(egin{matrix} xi_{1}\ vdots\ xi_{k}\ eta_{k} end{matrix} ight)$

得到的解和殘差

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

$lVert vec{r}^{left(k ight)} Vert=lvert eta_k vert$

同樣得到

$vec{p_k}=frac{1}{r_{k,k}}left(vec{v_k}-P_{k-1}vec{r_k^prime} ight)$

$xi_{k}=c_keta_{k-1}$

$eta_{k}=-s_keta_{k-1}$

初值

$xi_1=c_1lVert vec{b} Vert$

$eta_{1}=-s_1lVert vec{b} Vert$

$xi$ 和 $eta$ 的初值可統一為

$eta_{0}=lVert vec{b} Vert$

有限存儲空間的廣義最小殘差法

和阿諾爾迪演算法一樣，GMRES也需要保存前面所有的 $vec{v_1},dots,vec{v_{k-1}}$ 。為了限制存儲空間的增長，也有重啟法和截斷法兩種變種。

這裡再說說截斷法。當 $k>m$ 時，

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

有

$r_{1,k}=cdots=r_{k-m-2,k}=0$

即 $ar{R_k}$ 的最後一列比 $ar{H_k}$ 的最後一列多一個非零元素 $r_{k-m-1,k}$ ，對 $Q_{k-1}$ 也只需要存儲

$c_{k-m-1},dots,c_{k-1}$ 和 $s_{k-m-1},dots,s_{k-1}$ 。

最小殘差法（Minimal residual method, MINRES）

在這裡我們先講廣義最小殘差法，然後將狹義的最小殘差法作為特殊情況。和蘭佐斯演算法一樣，對於對稱的 $A$ ， $ar{H_k}$ 變為三對角矩陣 $ar{T_k}$ ，有

$r_{1,k}=cdots=r_{k-3,k}=0$

令

$alpha_i=t_{i,i}$

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

有

$r_{k-2,k}=s_{k-2}eta_k$

$r_{k-1,k}=c_{k-1}c_{k-2}eta_k+s_{k-1}alpha_k$

$r_{k,k}=sqrt{left(-s_{k-1}c_{k-2}eta_k+c_{k-1}alpha_k ight)^2+eta_{k+1}^2}$

$s_k=frac{eta_{k+1}}{sqrt{left(-s_{k-1}c_{k-2}eta_k+c_{k-1}alpha_k ight)^2+eta_{k+1}^2}}$

$c_k=frac{-s_{k-1}c_{k-2}eta_k+c_{k-1}alpha_k}{sqrt{left(-s_{k-1}c_{k-2}eta_k+c_{k-1}alpha_k ight)^2+eta_{k+1}^2}}$

以及前面以及推導的

$eta_{k+1}vec{v_{k+1}}=Avec{v_k}-alpha_kvec{v_k}-eta_{k}vec{v_{k-1}}$

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

$xi_{k}=c_keta_{k-1}$

$eta_{k}=-s_keta_{k-1}$

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

$lVert vec{r}^{left(k ight)} Vert=lvert eta_k vert$

用MATLAB代碼描述

function [x, iter, e] = minres_solve(A, b, x0, epsilon, iter_max) iter = 0; n = size(A, 1); % 循環隊列 p = zeros(n, 2); c = ones(2, 1); s = zeros(2, 1); v = zeros(n, 2); beta = zeros(2, 1); x = x0; w = b - A * x0; eta = norm(w); v(:, 1) = w / eta; while abs(eta) >= epsilon && iter < iter_max iter = iter + 1; e(iter) = abs(eta); i_0 = mod(iter + 1, 2) + 1; % k i_1 = mod(iter, 2) + 1; % k - 1, k + 1 i_2 = i_0; % k - 2 w = A * v(:, i_0) - beta(i_0) * v(:, i_1); alpha = dot(w, v(:, i_0)); w = w - alpha * v(:, i_0); beta(i_1) = norm(w); % eta_{k+1} v(:, i_1) = w / beta(i_1); % v_{k+1} r_2 = s(i_2) * beta(i_0); t_1 = c(i_2) * beta(i_0); r_1 = c(i_1) * t_1 + s(i_1) * alpha; t_0 = - s(i_1) * t_1 + c(i_1) * alpha; r_0 = sqrt(t_0 ^ 2 + beta(i_1) ^ 2); c(i_0) = t_0 / r_0; s(i_0) = beta(i_1) / r_0; p(:, i_0) = (v(:, i_0) - r_2 * p(:, i_2) - r_1 * p(:, i_1)) / r_0; xi = c(i_0) * eta; eta = - s(i_0) * eta; x = x + xi * p(:, i_0); endend