共軛梯度法及其淺顯分析

更多內(nèi)容請關(guān)注我的個人博客
浩源小站

背景

無處不在的線性方程組,昭示著一個優(yōu)秀的算法能帶來的巨大效益。對于簡單的低階方程組,只需要用普通的高斯消元法就能得到相當(dāng)不錯的結(jié)果,但是對于大型的方程組,常常與之相伴的就是稀疏性,如果還堅持使用高斯消元法,那么將帶來的是空間上的巨大浪費。為此,迭代法求解方程組便應(yīng)運而生。

共軛梯度法就是其中的佼佼者。

由來

二次型的最小值

我們熟知的線性方程組,常?？梢詫懗?img class="math-inline" src="https://math.jianshu.com/math?formula=Ax%3Db" alt="Ax=b" mathimg="1">的形式,實際上當(dāng) $A$ 是實對稱矩陣時,就是二次型 $f(x)=\frac{1}{2} x^TAx-b^Tx+c$ 對 $x$ 的導(dǎo)數(shù)為0時的表達式。

$f(x)=\frac{1}{2}x^TAx-b^Tx+c\\ f'(x)=\frac{1}{2}A^Tx+\frac{1}{2}Ax-b$

當(dāng) $A?$ 為實對稱矩陣的時候, $f'(x)=0?$ 就等價于 $Ax=b?$ .

而我們求解線性方程組的問題,也可以轉(zhuǎn)化成求解 $minf(x)$

舉個例子,當(dāng)

$A=\left[ \begin{matrix} 3&2\\ 2&6 \end{matrix} \right], b=\left[ \begin{matrix} 2\\8 \end{matrix} \right],c=0$

時, $f(x)?$ 的圖如下

image

等高線圖如下

image

看得出,此時的函數(shù)具有唯一的最小值。

由代數(shù)學(xué)的知識,我們知道,當(dāng)矩陣 $A?$ 為正定、半正定、負定、以及不定的時候,方程組 $Ax=b?$ 具有不同的解的情況。對應(yīng) $f‘(x)=0?$ 則是不同的最小值情況。

下面這幅圖說明了矩陣A的性質(zhì)對于 $f(x)$ 的影響，依次是正定、負定、半正定和不定的情況。

A的性質(zhì)對于f的影響

可以看到,當(dāng) $A$ 不定的時候, $f(x)$ 具有鞍點,這個時候無法通過導(dǎo)數(shù)為0來求得最值。

在接下來的共軛梯度方法中,我們都首先假定, $A$ 具有良好的性質(zhì),也就是對稱和正定。

最速下降法(梯度下降法)

如何得到 $f(x)$ 的最小值?

有一個形象的方法,從任意一點 $x_{(0)}$ 出發(fā),每次沿著當(dāng)前位置下滑最快的方向走,也就是該點處的梯度方向。

從而得到一系列的解序列: $x_{(1)},x_{(2)},\dots$ 直到兩次下降的差小于給定的誤差限為止。

由上面的結(jié)論, $f'(x_{(i)})=Ax_{(i)}-b$ ，下面記誤差(error)為 $e_{(i)}=x_{(i)}-x$ ，殘差(residual)為 $r_{(i)}=b-Ax_{(i)}$

于是有

$\begin{array}{**lr**} r_{(i)}=-f'(x_{(i)})\\ x_{(i+1)}=x_{(i)}+\alpha_{(i)}r_{(i)} \end{array}$

其中 $\alpha_{(i)}$ 為第 $i$ 步沿著方向 $r_{(i)}$ 的步長

如何選取步長 $\alpha_{(i)}$ 呢？樸素的想法是選取的步長盡量要讓 $f(x_{i+1})$ 的值最小,這樣才能更快的到達 $f(x)$ 的全局最小值。

從簡單的微積分知識中,我們把 $f(x_{(i+1)})$ 看作是 $\alpha_{(i)}$ 的函數(shù),要使得步長最合適,也就相當(dāng)于

$\fracu0z1t8os{d\alpha_{(i)}}f(x_{(i+1)})=0$

根據(jù)鏈式法則

$\begin{align*} \fracu0z1t8os{d\alpha_{(i)}}f(x_{(i+1)})&=f'(x_{(i+1)})^T\fracu0z1t8os{d\alpha_{(i)}}x_{(i+1)}\\ &=-r_{(i+1)}^Tr_{(i)} \end{align*}$

也就是說, $r_{(i)}$ 與 $r_{(i+1)}?$ 正交

$\begin{align*} r_{(i+1)}^Tr_{(i)}&=0\\ (b-Ax_{(i+1)})^Tr_{(i)}&=0\\ (b-A(x_{(i)}+\alpha_{(i)}r_{(i)}))^Tr_{(i)}&=0\\ (b-Ax_{(i)})^Tr_{(i)}-\alpha_{(i)}(Ar_{(i)})^Tr_{(i)}&=0\\ (b-Ax_{(i)})^Tr_{(i)}&=\alpha_{(i)}(Ar_{(i)})^Tr_{(i)}\\ \alpha_{(i)}&=\frac{r_{(i)}^Tr_{(i)}}{r_{(i)}^TAr_{(i)}}. \end{align*}$

也就是說,每一步的步長都可以根據(jù)當(dāng)前的殘差 $r_{(i)}$ 來確定

總結(jié)一下,最速下降法就是:

$\begin{align*} r_{(i)}&=b-Ax_{(i)}\\ \alpha_{(i)}&=\frac{r_{(i)}^Tr_{(i)}}{r_{(i)}^TAr_{(i)}}\\ x_{(i+1)}&=x_{(i)}+\alpha_{(i)}r_{(i)}. \end{align*}$

這樣迭代下去,直到達到事先給定的誤差限 $\epsilon$ 為止

當(dāng)然,有的時候,出于減小復(fù)雜度的考慮,我們會換用這種手段來求 $r_{(i)}?$

$r_{(i+1)}=r_{(i)}-\alpha_{(i)}Ar_{(i)}$

這個公式是在 $x_{(i+1)}=x_{(i)}+\alpha_{(i)}r_{(i)}$ 的兩邊左乘 $-A$ 再加上 $b$ 得到的。

這樣每次迭代過程我們只需要算一次矩陣乘法,可以減少運算。

雅可比迭代法

在開始共軛梯度法前,我們還需要再來看看它的同行們都是怎么干活的。

對于方程 $Ax=b$ ，我們把 $A$ 分成兩個矩陣 $D$ 和 $E$ ，其中 $D$ 是對角陣，而 $E$ 是剩下的部分，也就是說:

$\begin{align*} Ax&=b\\ (D+E)x&=b\\ Dx&=b-Ex\\ x&=D^{-1}(b-Ex). \end{align*}$

這里由于 $D$ 是對角陣,所以能保證它的可逆性。

這里為了便于分析,把上面的式子寫作

$x=Bx+z$

其中 $B=-D^{-1}E,z=D^{-1}b$

下面考察該方法迭代過程的收斂性,依舊把誤差(error)記作 $e_{(i)}$

$\begin{align*} x_{(i+1)}&=Bx_{(i)}+z\\ &=B(x+e_{(i)})+z\\ &=Bx+z+Be_{(i)}\\ &=x+Be_{(i)}\\ \end{align*}$

故而

$e_{(i+1)}=Be_{(i)}$

這里判斷迭代過程的斂散性主要看矩陣 $B$ 的性質(zhì),如果 $B$ 是實對稱矩陣,那么 $e_{(i)}$ 就可以用 $B$ 的特征向量線性表示,迭代過程就可以寫成:(下面以二階矩陣為例)

$\begin{align*} e_{(i+1)}&=Be_{(i)}\\ &=B^ie_{(0)}\\ &=B^i(x_1+x_2)\\ &=\lambda_1^ix_1+\lambda_2^ix_2 \end{align*}$

其中 $\lambda_1，\lambda_2$ 分別是矩陣 $B$ 的特征向量 $x_1,x_2$ 的特征值。

由此看得出,如果 $B$ 的最大的特征值小于1,這個迭代過程就收斂,不然迭代過程會發(fā)散。

而一個矩陣的最大特征值也就是它的譜半徑(spectral radius) $\rho(B)$

對于B不是實對稱矩陣的情況,也有類似的結(jié)論,只是說明過于復(fù)雜,在此不提。

總之，雅可比迭代法收斂的充分條件是迭代矩陣 $B$ 的譜半徑 $\rho(B)<1?$

下面還是以最初的方程舉例子，來分析雅可比迭代的具體過程，這有助于改進迭代過程，得到更好的算法

迭代過程現(xiàn)在是

$x_{(i+1)}=\left[ \begin{matrix} 0&-\frac{2}{3}\\ -\frac{1}{3}&0 \end{matrix} \right]x_{(i)}+ \left[ \begin{matrix} \frac{2}{3}\\ -\frac{4}{3} \end{matrix} \right]$

此時 $B$ 的特征向量以及與其對應(yīng)的特征值分別為

$$
x_1=\left[
\begin{matrix}
\sqrt2\1
\end{matrix}
\right],\lambda_1=-\frac{\sqrt2}{3}\
x_2=\left[
\begin{matrix}
-\sqrt2\1
\end{matrix}
\right],\lambda_1=\frac{\sqrt2}{3}\

下圖是迭代過程中,二次型函數(shù) $f(x)$ 的等高線圖，注意兩個箭頭分別表示特征向量

image

而下圖則是迭代過程

image

考慮每一次迭代與特征向量的關(guān)系

image

我們會發(fā)現(xiàn),每一次的迭代都朝著兩個特征向量的線性組合方向前進。

最速下降法的分析

前面我們得到了最速下降法的迭代過程，并且得到了一個十分有意思的結(jié)論

r_{(i+1)}^Tr_{(i)}=0

這意味著,最速下降法每一次的迭代過程,下降的方向都與上一次的方向正交。

用具體例子來表述是這樣的：

image

可以看到,(b)是最差的情況,需要蜿蜒前進很多次才接近解,而其他比較好的情況則兩次迭代就收斂了.

但即便迭代的次數(shù)不同,每一次迭代的方向都與上一次正交，這就是最速下降法的特點。

在(b)中還有一個情況值得注意,由于這是一個二維的圖形，故而如果每一次都和上一次正交,那么相鄰一次迭代的方向是平行的！

換句話來說,在最速下降法中,有一個很糟糕的現(xiàn)象,為了收斂到解附近,同樣的迭代方向可能走了不止一次。

這個現(xiàn)象不只是最速下降法中會出現(xiàn)，雅可比迭代法中也同樣的出現(xiàn)了，表現(xiàn)就是每一次的方向都是特征向量的線性組合，而且大多數(shù)情況下,前一次迭代過的特征向量方向上的分量,在下一次的迭代中還繼續(xù)存在。

共軛的想法

如果我能選取一系列線性無關(guān)的方向向量 $\{d_{(0)},d_{(1)},d_{(2)},\dots,d_{(n-1)}\}$ ,沿著每個方向只走一次，最后就能到達解 $x$ 處,是不是能很好的解決前面雅可比迭代和最速下降法的問題？

最直接的想法來源于直角坐標系，如果每一個方向都是正交的，既自然，有具有很好的性質(zhì)，是不是能得到很好的結(jié)果呢？

如果這樣選取了方向,那么應(yīng)該有誤差 $e_{(i+1)}$ 與方向 $d_{(i)}$ 正交

$\begin{align*} d_{(i)}^Te_{(i+1)}&=0\\ d_{(i)}^T(e_{(i)}+\alpha_{(i)}d_{(i)})&=0\\ \end{align*}$

所以

$\alpha_{(i)}=-\frac{d_{(i)}^Te_{(i)}}{d_{(i)}^Td_{(i)}}$

但是這里步長大小依賴于當(dāng)前的誤差 $e_{(i)}$ ，而不知道正確解 $x$ 是無法求出誤差來的,故而這種想法雖然很自然,但是無法得到有用的結(jié)果。

考慮到矩陣具有變換向量的作用,不妨記 $x^TAy=0$ 為向量 $x,y$ 關(guān)于矩陣 $A$ 共軛。

接下來選取的一系列方向向量都關(guān)于矩陣 $A$ 兩兩共軛, $\{ d_{(0)},d_{(1)},\dots,d_{(n-1)}\}$

我們從上一個想法得到靈感，要求這里的 $e_{(i+1)}$ 也與 $d_{(i)}$ 共軛，這是很自然的，只需要想象是對原本的空間做了一個變換,原本不正交的向量在這里正交,那么這里的要求就是合乎情理的了。

下面也一樣，對步長大小進行簡單的分析

$\begin{align*} \fracu0z1t8os{d\alpha_{(i)}}f(x_{(i+1)})&=f'(x_{(i+1)})^T\fracu0z1t8os{d\alpha_{(i)}}x_{(i+1)}\\ f'(x_{(i+1)})^Td_{(i)}&=0\\ d_{(i)}^Tr_{(i+1)}&=0\\ d_{(i+1)}^TA(e_{(i)}+\alpha_{(i)}d_{(i)})&=0\\ \end{align*}$

所以得

$\alpha_{(i)}=\frac{d_{i}^Tr_{(i)}}{d_{(i)}^TAd_{(i)}}$

這里得到的步長不僅可以算,而且可以證明,這樣迭代,對于n階的方程組,至多n步就可以收斂到正確解。

簡單的證明

首先，把誤差 $e_{(0)}$ 表示為方向向量的線性組合 $e_{(0)}=\sum_{j=0}^{n-1}\delta_jd_{(j)}$

兩側(cè)左乘 $d_{(k)}^TA$ 有：

$\begin{align*} d_{(k)}^TAe_{0}&=\sum_{j=0}^{n-1}\delta_{j}d_{(k)}^TAd_{(j)}\\ &=\delta_{k}d_{(k)}^TAd_{(k)} \end{align*}$

上式利用了之前提到的共軛性質(zhì),由此得

$\delta_{k}=\frac{d_{(k)}^TAe_{(k)}}{d_{(k)}^TAd_{(k)}}$

與上面的步長比較，有

$\alpha_{(k)}=-\delta_k$

于是對于第 $i$ 步迭代的誤差 $e_{(i)}$

$\begin{align*} e_{(i)}&=e_{(0)}+\sum_{j=0}^{i-1}\alpha_{(i)}d_{(j)}\\ &=\sum_{j=0}^{n-1}\delta_{(j)}d_{(j)}-\sum_{j=0}^{i-1}\delta_{(j)}d_{(j)}\\ &=\sum_{j=i}^{n-1}\delta_{(j)}d_{(j)} \end{align*}$

由此可知,當(dāng) $i=n$ 的時候, $e_{(n)}=0$

$Q.E.D$

并且我們能從上面的證明中得到這樣一個事實,第 $i$ 迭代時,有關(guān)前 $i-1$ 次迭代的方向上都已經(jīng)沒有誤差了,也就是這樣的方法,每次都完全消除某個方向上的誤差，這樣至多 $n$ 步得到精確解就是很自然的了。

共軛梯度法

上面的由來中，已經(jīng)充分做好了準備來得到共軛梯度法，現(xiàn)在只剩一個問題，那一組兩兩共軛的方向向量該如何選取?

選取了一組基后，只需要進行Gram-Schmidt正交化就可以得到一組正交基，那么同樣的也可以用這樣的方法來得到一組共軛的方向向量。

問題是，如果隨便選取一組線性無關(guān)的方向基，就進行Gram-Schmidt共軛化是否足夠好呢？

我們總希望這組方向基不僅容易得到，而且具有很好的性質(zhì)，減少計算量。

考慮到之前得到的結(jié)論，每一次迭代之間的殘差都是相互正交的，我們不妨就選殘差

$\{ r_{(0)},r_{(1)},r_{(2)},\dots,r_{(n-1)} \}$

作為共軛化之前的基。由于使用共軛化的方向向量來迭代至多只有n步，且每步都將該方向上的誤差消滅掉，故而這一組基不僅線性無關(guān)，而且由于它們是殘差，還具有正交的良好性質(zhì)。

推導(dǎo)

首先，對于第 $i?$ +1次迭代后的殘差 $r_{(i+1)}?$ 有

$\begin{align*} r_{(i+1)}&=-Ae_{(i+1)}\\ &=-A(e_{(i)}+\alpha_{(i)}d_{(i)})\\ &=r_{(i)}-\alpha_{(i)}Ad_{(i)}. \end{align*}$

接著我們來推導(dǎo)一下Gram-Schmidt共軛化的公式

第 $i$ 次的方向向量為 $d_{(i)}$ ，我們希望 $d_{(i)}$ 從 $r_{(i)}+\{ d_{(0)},d_{(1)},d_{(2)},\dots,d_{(i-1)}\}$ 中得到

,也就是

$d_{(i)}=r_{(i)}+\sum_{k=0}^{i-1}\beta_{ik}d_{(k)},(i>k)$

在上式的兩側(cè)同左乘上 $d_{(i)}^TA$ 有

$\begin{align*} d_{(i)}^TAd_{(j)}&=r_{(i)}^TAd_{(j)}+\sum_{k=0}^{i-1}\beta_{ik}d_{(k)}^TAd_{(j)}\\ 0&=r_{(i)}^TAd_{(j)}+\beta_{ij}d_{(j)}^TAd_{(j)},\quad i>j\\ \beta_{ij}&=-\frac{r_{(i)}^TAd_{(j)}}{d_{(j)}^TAd_{(j)}} \end{align*}$

下面來結(jié)合殘差 $r_{(j+1)}=r_{(j)}-\alpha_{(j)}d_{(j)}$ ,繼續(xù)做一些推演

在兩側(cè)同左乘 $r_{(i)}^T$

$$
\begin{align*}
r_{(i)}^{Tr_{(j+1)}&=r_{(i)}}Tr_{(j)}-\alpha_{(j)}r_{(i)}^TAd_{(j)}\
\alpha_{(j)}r_{(i)}^{TAd_{(j)}&=r_{(i)}}Tr_{(j)}-r_{(i)}^Tr_{(i)}\

\end{align*}
$$

故而

$r_{(i)}^TAd_{(j)}= \left\{ \begin{align*} \frac{1}{\alpha_{(i)}}r_{(i)}^Tr_{(i)},\quad i&=j\\ -\frac{1}{\alpha_{(i-1)}}r_{(i)}^Tr_{(i)},\quad i&=j+1\\ 0,\quad otherwise \end{align*} \right.$

結(jié)合 $\beta_{ij}$ 的值，有

$$
\beta_{ij}=
\left{
\begin{align}
\frac{1}{\alpha_{i-1}}\frac{r_{(i)}^{Tr_{(i)}}{d_{(i-1)}}TAd_{(i-1)}}&,\quad i=j+1\
0&,\quad i>j
\end{align}

\right.
$$

考慮到完全可以把 $\beta_{ij}$ 寫成 $\beta_i$ ，因為 $j$ 只能取固定的值。

神奇的事情發(fā)生了,原本需要 $i$ 個 $\beta?$ 才能確定的方向向量,在共軛化的條件下,只需要當(dāng)前的數(shù)據(jù)和前一步的數(shù)據(jù)就可以得到,而不必存儲之前所有走過的路徑信息。

結(jié)合前面得到的 $\alpha_(i)=\frac{d_{(i)}^Tr_{(i)}}{d_{(i)}^TAd_{(i)}}$

$\beta_{i}=\frac{r_{(i)}^Tr_{(i)}}{d_{(i-1)}^Tr_{(i-1)}}$

實際上，上面的式子還能夠?qū)懗筛悠恋臉幼?/p>

先由 $d_{(i)}=r_{(i)}+\sum_{k=0}^{i-1}\beta_{ik}d_{(k)},(i>k)$

做其與 $r_{(j)}$ 的內(nèi)積有

$\begin{align*} d_{(i)}^Tr_{(j)}&=r_{(i)}^Tr_{(j)}+\sum_{k=0}^{i-1}\beta_{ik}d_{(k)}^Tr_{(j)}\\ 0&=r_{(i)}^Tr_{(j)},\quad i<j \end{align*}$

令 $j=i$ 有

$d_{(i)}^Tr_{(i)}=r_{(i)}^Tr_{(i)}$

帶入 $\beta_i$ 中有

$\beta_i=\frac{r_{(i)}^Tr_{(i)}}{r_{(i-1)}^Tr_{(i-1)}}$

這樣，Gram-Schmidt共軛化就完美的實現(xiàn)了，不僅實現(xiàn)了每一個方向只迭代一次，而且需要存儲的數(shù)據(jù)只有上一步的殘差。

下面是共軛梯度法涉及到的所有公式

$\begin{align*} d_{(0)}&=r_{(0)}=b-Ax_{(0)}\\ \alpha_{(i)}&=\frac{r_{i}^Tr_{(i)}}{d_{(i)}^TAd_{(i)}}\\ x_{(i+1)}&=x_{(i)}+\alpha_{(i)}d_{(i)}\\ r_{(i+1)}&=r_{(i)}-\alpha_{(i)}Ad_{(i)}\\ \beta_{(i+1)}&=\frac{r_{(i+1)}^Tr_{(i+1)}}{r_{(i)}^Tr_{(i)}}\\ d_{(i+1)}&=r_{(i+1)}+\beta_{(i+1)}d_{(i)} \end{align*}$

最后給出共軛梯度法的Matlab實現(xiàn)

function x=Conjugate_gradient(A,b,x0)
x=x0;
d=b-A*x0;
r=d;
for k=0:size(A)-1
    if r==0
        break;
    end
    alpha=(r'*r)/(d'*A*d);
    x=x+alpha*d;
    temp=r'*r;
    r=r-alpha*A*d;
    belta=(r'*r)/(temp);
    d=r+belta*d;
end
end

分析

正交性是共軛梯度法成功的關(guān)鍵,注意到每一次迭代中，新的殘差 $r_{(i+1)}$ 與前面所有的 $r_{(i)}$ 正交，如果一個 $r_{(i)}=0$ 那么 $Ax_{(i)}=b$ 方程已經(jīng)解完了。

否則在n步迭代后， $r_{(n)}$ 和n個兩兩正交的向量 $\{ r_{(0)},r_{(1)},\dots,r_{(n-1)}\}$ 所張成的空間正交,而 $r_{(0)},r_{(1)},\dots,r_{(n-1)}$ 這n個向量是 $R^n$ 空間中的所有正交向量,由此 $r_{(n)}$ 必須是零向量,所以 $Ax_{(n)}=b$

比較

共軛梯度法從某種程度上要簡單于高斯消元法,不必考慮行和列的相消,而且代碼實現(xiàn)也十分簡潔。下面來比較一下共軛梯度法和高斯消元法在復(fù)雜度上的優(yōu)勢。

共軛梯度法的每一次迭代都要做一次矩陣向量乘法和一些向量內(nèi)積的計算，復(fù)雜度為 $n^2$ ，當(dāng)做完n次迭代后，復(fù)雜度變成了 $n^3$ ，而高斯消元法只是 $\frac{1}{3}n^3$ 左右，從這一點上，當(dāng)矩陣不是稀疏矩陣的時候，共軛梯度法沒有優(yōu)勢，而當(dāng)矩陣變得稀疏的時候，n大的驚人，高斯消元法如果要得到解，就必須要做完所有的運算才能得到解，這樣需要消耗大量的資源。而共軛梯度法不一樣，它每一次迭代都能在某個分量上得到解，并且可以通過殘差來度量解的精確情況，我們并不需要將算法進行到底，只需要解達到精度要求就可以退出。