線性方程組(4)-變分原理與共軛梯度法

02-27

上一節講到了靜態迭代法。因為收斂太慢，更常用的還是非靜態迭代法。

這裡先從變分原理入手。考慮函數

$Phileft(vec{x} ight)=frac{1}{2}vec{x}cdot Avec{x}-vec{x}cdot vec{b}$

若 $A$ 對稱，則函數的梯度

$abla Phileft(vec{x} ight)=Avec{x}- vec{b}$

若 $A$ 對稱正定，則解方程 $Avec{x}= vec{b}$ 等價於求函數 $Phileft(vec{x} ight)$ 的最小值。方程的殘差 $vec{r}=vec{b}-Avec{x}$ 即為函數 $Phileft(vec{x} ight)$ 的負梯度。

一維優化問題

假設已經得到了 $k-1$ 步迭代的解 $vec{x}^{left(k-1 ight)}$ ，在第 $k$ 步中，沿 $vec{p_k}$ 方向搜索最優解，即

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

代入函數表達式可得

$Phileft(vec{x}^{left(k ight)} ight)=frac{1}{2}vec{p_k}cdot Avec{p_k}alpha_k^2-left(vec{b}-Avec{x}^{left(k-1 ight)} ight)cdotvec{p_k}alpha_k+frac{1}{2}vec{x}^{left(k-1 ight)}cdot Avec{x}^{left(k-1 ight)}-vec{x}^{left(k-1 ight)}cdot vec{b}$

最優的 $alpha_k$ 取值應是

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

得到新的解後，新的殘差

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

最速下降法（Steepest descent method）

搜索方向的最簡單的取法就是取 $vec{p_k}=vec{r}^{left(k-1 ight)}$ ，即負梯度方向，這就是優化問題中很常用的最速下降法。用MATLAB代碼描述

function [x, iter] = sd_solve(A, b, x0, epsilon, iter_max) iter = 0; x = x0; r = b - A * x0; rr = dot(r, r); while sqrt(rr) >= epsilon && iter < iter_max iter = iter + 1; Ar = A * r; alpha = rr / dot(r, Ar); x = x + alpha * r; r = r - alpha * Ar; rr = dot(r, r); endend

注意提取出了 $vec{r}cdot vec{r}$ 和 $Avec{r}$ 兩個公共表達式，避免了重複計算。

每一步迭代需要計算一次矩陣向量乘，兩次向量內積，兩次向量線性組合。

注意到 $vec{r}^{left(k ight)}$ 的計算可能有誤差積累，干擾迭代終止條件的判斷，使得迭代在達到要求的精度之前終止。可能的解決方法有：用 $vec{r}^{left(k ight)}=vec{b}-Avec{x}^{left(k ight)}$ 來計算更精確的殘差，引入一次額外的矩陣向量乘；或用更小的 $varepsilon$ ；或用其它迭代終止條件。

假設精確解為 $vec{x}^*$ ，則有結論

$lVert vec{x}^{left(k ight)}-vec{x}^{*} Vert leqslant frac{lambda_{max}-lambda_{min}}{lambda_{max}+lambda_{min}}lVert vec{x}^{left(k-1 ight)}-vec{x}^{*} Vert=frac{condleft(A ight)-1}{condleft(A ight)+1}lVert vec{x}^{left(k-1 ight)}-vec{x}^{*} Vert$

共軛梯度法（Conjugate gradient method, CG）

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

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

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

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

滿足

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

將 $vec{x}^{left(k-1 ight)}$ 和 $vec{x}^{left(k ight)}$ 在子空間中展開

$vec{x}^{left(k-1 ight)}=sum_{i=1}^{k-1}alpha_ivec{p_{i}}$

$vec{x}^{left(k ight)}=sum_{i=1}^{k}alpha_ivec{p_{i}}$

有

$Phileft(vec{x}^{left(k ight)} ight)=frac{1}{2}vec{p_k}cdot Avec{p_k}alpha_k^2-vec{b}cdotvec{p_k}alpha_k-Avec{x}^{left(k ight)}cdotvec{p_k}alpha_k+ Phileft(vec{x}^{left(k-1 ight)} ight)$

對於 $i=1,2,cdots,k-1$

$frac{partial}{partialalpha_i} Phileft(vec{x}^{left(k ight)} ight)=-alpha_kvec{p_k}cdot Avec{p_i}+ frac{partial}{partialalpha_i}Phileft(vec{x}^{left(k-1 ight)} ight)$

由前面的條件

$frac{partial}{partialalpha_i} Phileft(vec{x}^{left(k ight)} ight)=frac{partial}{partialalpha_i} Phileft(vec{x}^{left(k-1 ight)} ight)=0$

得到

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

這就是共軛梯度法，即選擇兩兩 $A$ 共軛的一系列搜索方向。僅有這個約束還不足以確定 $vec{p_k}$ ，為了得到一個確定的搜索方向，我們選擇將 $vec{r}^{left(k-1 ight)}$ 對 $Avec{p_i}$ 正交化

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

其中 $gamma_i$ 是使得 $lVertvec{p_k} Vert$ 最小的係數（垂線最短）。

為了下面的推導，先提出一個引理，證明留作練習（提示：數學歸納法）

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

由

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

可知 $vec{r}^{left(k-1 ight)}=- abla Phileft(vec{x}^{left(k-1 ight)} ight)$ 正交於 $spanleft{vec{p_1},vec{p_2},dots,vec{p_{k-1}} ight}$ ，再根據引理可以推出下面的推導中所需要的所有正交關係。

因為

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

所以

$vec{p_k}=vec{r}^{left(k-1 ight)}-frac{gamma_{k-1}}{alpha_{k-1}}left(vec{r}^{left(k-1 ight)}-vec{r}^{left(k-2 ight)} ight)-sum_{i=1}^{k-2}gamma_i Avec{p_i}$

$=left(1-frac{gamma_{k-1}}{alpha_{k-1}} ight)vec{r}^{left(k-1 ight)}+frac{gamma_{k-1}}{alpha_{k-1}}left(vec{r}^{left(k-2 ight)}-sum_{i=1}^{k-2}frac{alpha_{k-1}}{gamma_{k-1}}gamma_i Avec{p_i} ight)$

由正交性

$lVertvec{p_k} Vert=left(1-frac{gamma_{k-1}}{alpha_{k-1}} ight)lVertvec{r}^{left(k-1 ight)} Vert+frac{gamma_{k-1}}{alpha_{k-1}}lVertvec{r}^{left(k-2 ight)}-sum_{i=1}^{k-2}frac{alpha_{k-1}}{gamma_{k-1}}gamma_i Avec{p_i} Vert$

第二項在 $vec{r}^{left(k-2 ight)}-sum_{i=1}^{k-2}frac{alpha_{k-1}}{gamma_{k-1}}gamma_i Avec{p_i}=vec{p_{k-1}}$ 時取最小值，所以

$vec{p_k}=left(1-frac{gamma_{k-1}}{alpha_{k-1}} ight)vec{r}^{left(k-1 ight)}+frac{gamma_{k-1}}{alpha_{k-1}}vec{p_{k-1}}$

又因為我們只關心 $vec{p_k}$ 的方向，不關心其大小，可以寫成

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

$alpha_k$ 的計算式可以繼續簡化

$alpha_k =frac{ vec{r}^{left(k-1 ight)}cdotvec{p_k}}{vec{p_k}cdot Avec{p_k}}=frac{ vec{r}^{left(k-1 ight)}cdotleft(vec{r}^{left(k-1 ight)}+eta_{k}vec{p_{k-1}} ight)}{vec{p_k}cdot Avec{p_k}}=frac{ vec{r}^{left(k-1 ight)}cdotvec{r}^{left(k-1 ight)}}{vec{p_k}cdot Avec{p_k}}$

因為

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

所以

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

經過上面一大段推導，共軛梯度法終於呈現在我們眼前，用MATLAB代碼描述

function [x, iter] = cg_solve(A, b, x0, epsilon, iter_max) iter = 0; x = x0; r = b - A * x0; p = r; rr_prev = dot(r, r); % 上一步得到的解的殘差平方和 while sqrt(rr_prev) >= epsilon && iter < iter_max iter = iter + 1; Ap = A * p; alpha = rr_prev / dot(p, Ap); x = x + alpha * p; r = r - alpha * Ap; rr_curr = dot(r, r); % 這一步得到的解的殘差平方和 beta = rr_curr / rr_prev; p = r + beta * p; % 下一步的p rr_prev = rr_curr; endend

每一步迭代需要計算一次矩陣向量乘，兩次向量內積，三次向量線性組合。

注意，因為要求

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

所以初值 $vec{x}^{left(0 ight)}=0$

如果能找到一個更接近精確解的非零初值 $vec{x_0}$ ，則等價於解 $Avec{z}=vec{b}-Avec{x_0}$ 中的 $vec{z}$ ，最後得到原方程的解 $vec{x}=vec{z}+vec{x_0}$ ；由於 $vec{x}^{left(k ight)}$ 只涉及加法，根據加法交換律和結合律，這個累加放在開頭還是結尾不影響結果，於是可以取 $vec{x}^{left(0 ight)}=0+vec{x_0}=vec{x_0}$ 。

按照前面的理論，共軛梯度法迭代到第 $n$ 步時

$Phileft(vec{x}^{left(n ight)} ight)=min_{vec{x}in spanleft{vec{p_1},vec{p_2},dots,vec{p_{n}} ight}} Phileft(vec{x} ight)=min_{vec{x}in mathbb{R}^n} Phileft(vec{x} ight)$

即得到精確解。也就是說，共軛梯度法是一種理論上的直接解法。

但是和LU分解法一樣，在實際計算中，由於浮點誤差的積累，理論上的精確解沒有意義。我們只需要把共軛梯度法看成迭代法即可，達到迭代終止條件之前的迭代步數小於 $n$ 或者大於 $n$ 都沒關係。這個現象在之後還會再遇到。

共軛梯度法的收斂性

$lVert vec{x}^{left(k ight)}-vec{x}^{*} Vert leqslant frac{sqrt{condleft(A ight)}-1}{sqrt{condleft(A ight)}+1}lVert vec{x}^{left(k-1 ight)}-vec{x}^{*} Vert$

對於條件數較大的矩陣，共軛梯度法的收斂速度比最速下降法快很多。

共軛梯度法要求係數矩陣正定，因為只有正定才存在最小值；實際上對於不定矩陣，共軛梯度法也很有可能收斂到梯度為零所對應的鞍點，不過不能從理論上保證收斂性。

共軛梯度法的變種

將共軛梯度法推廣到非對稱矩陣，最簡單的辦法是把方程

$Avec{x}=vec{b}$

轉換為

$A^T Avec{x}=A^T vec{b}$

由於 $A^T A$ 對稱正定，可以使用共軛梯度法解該方程。這就是正規方程的共軛梯度法。由於

$condleft(A^T A ight)=condleft(A ight)^2$

所以收斂較慢。

雙共軛梯度法（Bi-conjugate gradient method）也可以應用於非對稱矩陣。計算兩組 $vec{r}$ 和 $vec{p}$

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

$vec{s}^{left(k ight)}=vec{s}^{left(k-1 ight)}-alpha_k A^Tvec{q_k}$

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

$vec{q_k}=vec{s}^{left(k-1 ight)}+eta_{k}vec{q_{k-1}}$

其中

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

$eta_k=frac{vec{s}^{left(k-1 ight)}cdot vec{r}^{left(k-1 ight)}}{vec{s}^{left(k-2 ight)}cdotvec{r}^{left(k-2 ight)}}$

維護 $vec{r}$ 的雙正交性和 $vec{p}$ 的雙共軛性，對於 $i ot=j$

$vec{s}^{left(i ight)}cdot vec{r}^{left(j ight)}=vec{q_i}cdot Avec{p_j}=0$

雙共軛梯度法有數值不穩定的現象，在迭代一定步數後，由於浮點誤差積累，雙正交性和雙共軛性不再保持。穩定雙共軛梯度法（Bi-conjugate gradient stablilized method）通過將一維搜索擴增到高維來修正這一誤差，當搜索維數高時收斂速度提升，但單步迭代計算量和消耗的存儲空間也將提升。