共軛梯度法（一）：線性共軛梯度

本文開始介紹一種全新的最佳化方法——共軛梯度法。

最初，共軛梯度法是用來求解線性方程

的，稱為線性共軛梯度法。後來，有人把這種方法擴充套件到了非線性最佳化問題中，稱為非線性共軛梯度法。本文從基本的線性共軛梯度講起，介紹它的思想和演算法流程。

問題描述

共軛梯度法並不適用於任意線性方程，它要求係數矩陣對稱且正定。我們會在後面逐漸提到為什麼會有這樣苛刻的要求。

比較巧的是，在非線性最佳化問題中，所有的二階方法的每次迭代都需要求解一個形如

這樣的方程。比如，前面提到的牛頓系列方法，需要求解

$B_k p_k=\nabla f_k$

，SLAM中常用的高斯牛頓法，需要求解

$J^T J \Delta x = - J^T f$

，這些係數矩陣

和

通常都滿足對稱正定的要求（如果不正定，則想辦法讓它們正定，比如加上一個單位陣）。

之前，我們都預設直接對係數矩陣求逆，得到

$x = A^{-1} b$

，而沒有考慮求逆的效能代價。即便是透過矩陣分解的方式求逆，計算量仍然很大。線性共軛梯度法的提出，讓我們能夠透過迭代的方式求

，而不必對係數矩陣求逆。

具體如何實現呢？我們先額外構造一個二次函式

$\begin{equation} \phi \left( x \right) = \frac{1}{2} x^TA x - b^T x \tag{5.2} \end{equation}$

如果對該函式求導，並令導數為零，我們得到

$\begin{equation} \nabla \phi \left( x \right) = A x - b = 0 \tag{5.2} \end{equation}$

你會發現，使得該函式導數為0的

恰好是線性方程

的解。換句話說，

$\phi \left( x \right)$

的極小值點（

是正定矩陣，意味著

$\phi \left( x \right)$

有極小值）恰好是線性方程

的解。於是，我們可以把求解線性方程的問題轉化為求二次函式極值的問題。

尋找思路

如何求二次函式極值？可以用之前講過的最速下降法、牛頓法之類的方法嗎？

顯然，牛頓法我們不用考慮了，因為剛才提到，求解線性方程其實是牛頓法的每次迭代都要面對的問題，所以它們不是一個層次上的問題。最速下降法倒是可以考慮，我們只需要求解當前位置的導數

$\nabla \phi \left( x_k \right)$

，然後在該方向上確定一個最優的步長，不斷迭代即可。回顧一下線搜尋（二）那篇文章中給出的最速下降法的收斂路徑，如下圖所示。

對於一個簡單的二次函式，最速下降法竟然走出這麼蜿蜒曲折的路徑，這要用在求解線性方程上，恐怕比直接矩陣求逆還要慢。

我們仍然以上圖所示的二維目標函式為例，最好的結果是走出下面這條路徑。

我們這樣解讀這條新的路徑：整個路徑只有兩段，即兩次迭代，對應著目標函式只有兩維。第一次迭代從

到

，沿著座標軸

的方向前進，達到了目標函式在該方向上可能達到的最優值。第二次迭代從

到

，沿著座標軸

的方向前進，達到了目標函式在該方向上可能達到的最優值。兩次迭代後就達到了極小值點。

重新審視這一路線，可以發現，每次迭代的前進方向是互相正交的，因為都是沿著座標軸前進的。但最速下降法的搜尋路徑也是互相正交的，為什麼就不行呢？如果我們把目標函式稍微旋轉一下，就會發現問題的關鍵。

這是另一個二次函式，它與前面的二次函式很不一樣，橢圓的軸不與座標軸平行。（這裡補充一點，之所以我們畫出的二次函式的圖形都是橢圓，是因為我們只考察係數矩陣為對稱陣的情況。）此時，仍然沿著座標軸搜尋，這就是不折不扣的最速下降法了。這說明，能不能在兩次迭代內達到極小值，與是否沿著座標軸搜尋沒什麼關係，關鍵在於橢圓的軸是否與搜尋方向一致。在第一種情況中，橢圓的軸與座標軸平行，搜尋方向也和座標軸平行，所以兩次迭代即可。第二種情況中，橢圓的軸不與座標軸平行，搜尋方向依舊與座標軸平行，所以兩次迭代是不夠的。這給了我們一些提示，如果我們先算出橢圓長短軸的方向，然後沿著這些方向去搜索，就能在兩個迭代後達到極小值。

從上面的二維函式推廣到N維函式，想象N維空間中的橢圓，我們仍然可以計算出每一維的軸向，然後沿著這些軸搜尋，就能確保在N次迭代後達到極小值。那如何求橢圓的軸呢？很簡單，橢圓的軸其實就是係數矩陣的特徵向量，對

做特徵值分解就可以了。說到這裡，其實有點怪怪的，特徵值分解的計算量恐怕不小吧，為了避免求逆卻整出來個特徵值分解，豈不是捨近求遠。

那麼，有沒有不求特徵值，也能在N次迭代達到極小值點的方法呢？比如，能不能走出下面這樣的路徑。

第一次迭代前進的方向並不與橢圓的軸平行，在這種情況下，第二次迭代卻能一次達到最優點。異想天開？當然不是，共軛梯度法就是要走出這樣的路線。

我們用數學工具探究一下這種路徑需要滿足的條件。對於二次函式中的任意一個點來說，該點與函式極值點的連線正是牛頓法中每次迭代的線搜尋方向。回顧線搜尋與信賴域中的式（2。14），套用這裡的係數矩陣，可以得到圖中的向量

$\overrightarrow{x_1 x^*} = - A^{-1} \nabla \phi \left( x_1 \right)$

。如果假設起始點

的搜尋方向是最速下降方向，即

$\overrightarrow{x_0 x_1} = - \nabla \phi \left( x_0 \right)$

，於是有

$\begin{aligned} \overrightarrow{x_0 x_1}^T A \space \overrightarrow{x_1 x^*} &= \nabla \phi \left( x_0 \right) A A^{-1} \nabla \phi \left( x_1 \right) \\ &= \nabla \phi \left( x_0 \right) \nabla \phi \left( x_1 \right) \\ &= 0 \end{aligned} \\$

其中，最後兩個梯度的乘積等於0是因為最速下降法的每次迭代方向都是垂直的。我們從二維例子中推出了一個看起來有些意義的結果，這個結果表明，連續的兩次迭代的前進方法關於係數矩陣

正交。也就是說，其中的一個方向經過

的變換後與另一個方向正交。

接下來的問題是如何把這一結論推廣到N維目標函式中去，我們先不關注證明，直接看一下結論。

共軛梯度法

對於一組向量

$\left\{ p_0, p_1, ..., p_{n-1} \right\}$

，如果任意兩個向量間滿足

$\begin{equation} p_i^T A p_j = 0 \tag{5.4} \end{equation}$

則稱這組向量與對稱正定矩陣

共軛。

可見，剛才我們發現的正是共軛的性質，只不過，推廣到高維空間後，任意兩個向量之間都要關於

正交，不僅僅是相鄰的。

共軛向量的好處正如我們前面發現的，依次沿著每個向量最佳化，最終可以在N個迭代後達到極小值。書中有嚴格的證明，但我們無意在此著墨，我認為結合上面的例子便可以理解其中的思想。

現在的問題是，給定任意一個對稱正定矩陣，如何求取一組共軛向量？

顯然，

的特徵向量就是一組共軛向量，因為

$u_i^T A u_j = \lambda u_i u_j = 0$

對於

的任意兩個特徵向量

和

都成立，滿足式（5。4），這就從公式上驗證了我們前面的發現。不過，之前也說了，實際中是不能用特徵向量的，因為計算太複雜。（請注意，這裡我們用到了

這一性質，只有對稱矩陣有此性質，呼應了文章開頭提到的對矩陣

的特殊要求。這一性質貫穿了整個共軛梯度法的證明過程。）

實際中，需要一種不依賴於矩陣分解，可以在每次迭代時以較小的代價計算出

的方法。回顧上文的最後一張圖，當時我們假設起始點

的搜尋方向是最速下降方向。這其實不是隨便假設的，這樣做會使得後面的計算變得簡單。接下來，我們就從

$p_0 = - \nabla \phi \left( x_0 \right)$

開始，推匯出後續每次迭代所需的

。

滿足式（5。4）的共軛向量有無數個，而我們只需要取其中的一個。在

的時刻，我們已知

和

$\nabla \phi \left( x_1 \right)$

，注意到這兩個向量一定是線性無關的，因為

是該方向上使得目標函式最小的步長，所以在該方向上目標函式不可能繼續下降。那麼就可以在這兩個向量構成的平面上生成新的

，令

$p_1 = - \nabla \phi \left( x_0 \right) + \beta_1 p_0$

，將

和

帶入式（5。4），即可求得

$\beta_1 = \frac{p_0^T A \nabla \phi \left( x_1 \right)}{p_0^T A p_0}$

，從而得到

。

在此基礎上，結合數學歸納法，我們可以證明，令

$\begin{equation} p_k = - \nabla \phi \left( x_k \right) + \beta_k p_{k-1} \tag{5.12} \end{equation}$

其中

$\beta_k = \frac{p_{k-1}^T A \nabla \phi \left( x_k \right)}{p_{k-1}^T A p_{k-1}} \\$

即可使得生成的所有

滿足式（5。4）的要求，從而得到一組共軛向量。有了每次迭代的方向

，就可以計算出

$\phi \left(x_k \right)$

在該方向上的最優步長。計算方式很簡單，我們用

$\alpha_k$

替換

$\phi \left( x_{k+1} \right)$

中的

$x_{k+1}$

，得到

$\phi \left( x_k + \alpha_k p_k \right)$

，將其對

$\alpha_k$

求導，令導數等於0，即可得到最優步長

$\alpha_k = - \frac{\nabla \phi \left( x_k \right)^T p_k}{p_k^T A p_k} \tag{5.6}\\$

於是，整個共軛梯度法的流程整理如下

$\begin{aligned} \alpha_{k} & \leftarrow-\frac{r_{k}^{T} p_{k}}{p_{k}^{T} A p_{k}} \\ x_{k+1} & \leftarrow x_{k}+\alpha_{k} p_{k} \\ r_{k+1} & \leftarrow A x_{k+1}-b \\ \beta_{k+1} & \leftarrow \frac{r_{k+1}^{T} A p_{k}}{p_{k}^{T} A p_{k}} \\ p_{k+1} & \leftarrow-r_{k+1}+\beta_{k+1} p_{k} \\ k & \leftarrow k+1 \end{aligned} \tag{5.13}$

其中，用

替代了

$\nabla \phi \left( x_k \right)$

。從

開始，

$r_0 = A x_0 - b, \space p_0 = -r_0$

，依次迭代下去。

實際中，為了降低計算量，一般會對式（5。13）稍做改動，得到下面的演算法流程

$\begin{aligned} \alpha_{k} & \leftarrow \frac{r_{k}^{T} r_{k}}{p_{k}^{T} A p_{k}} \\ x_{k+1} & \leftarrow x_{k}+\alpha_{k} p_{k} \\ r_{k+1} & \leftarrow r_{k}+\alpha_{k} A p_{k} \\ \beta_{k+1} & \leftarrow \frac{r_{k+1}^{T} r_{k+1}}{r_k^T r_k} \\ p_{k+1} & \leftarrow-r_{k+1}+\beta_{k+1} p_{k} \\ k & \leftarrow k+1 \end{aligned} \tag{5.23}$

這一版本的共軛梯度演算法涉及的矩陣乘法和加法會更少一些。

收斂速度

最佳化演算法的收斂速度是個值得關心的問題。前面提到了，共軛梯度法一定會在N次迭代內收斂到極值。但在N次迭代的過程中，收斂速度怎樣呢？

我們用線搜尋（二）：收斂性和收斂速度中提到過的向量的矩陣範數來衡量共軛梯度法的收斂速度。比如，評估

$\| x_{k+1} - x^* \|_A^2$

隨著

的增加的變化率。這裡省略證明過程，直接給出結論

$\begin{aligned} \left\|x_{k+1}-x^{*}\right\|_{A}^{2} \le\min _{P_{k}} \max _{1 \leq i \leq n}\left[1+\lambda_{i} P_{k}\left(\lambda_{i}\right)\right]^{2}\left\|x_{0}-x^{*}\right\|_{A}^{2} \end{aligned} \tag{5.32}$

其中，

$\lambda_i$

是矩陣

的第

個特徵值，

$P_k\left( \lambda_i \right)$

是任意一個關於

$\lambda_i$

的

階多項式。這個等式可以理解為，演算法在第

次迭代後，殘差的最大值與一個

階多項式的值有關，這個

階多項式是所有

階多項式裡面結果最小的那個，但每個多項式的輸入必須是使其取得最大值的

$\lambda_i$

。單從這個不等式其實是看不出什麼東西的，但卻可以推匯出一個驚人的結論。

如果矩陣

的N個特徵值只能取

個互不相同的值，這

個值分別是

$\tau_1 \lt \tau_2 \lt \cdots \lt \tau_r$

。我們定義一個

階多項式

$Q_{r}(\lambda)=\frac{(-1)^{r}}{\tau_{1} \tau_{2} \cdots \tau_{r}}\left(\lambda-\tau_{1}\right)\left(\lambda-\tau_{2}\right) \cdots\left(\lambda-\tau_{r}\right) \\$

可以發現該多項式滿足

$Q_r\left(0\right) = 1$

。我們再定義另一個多項式

$\bar{P}_{r-1}(\lambda)=\left(Q_{r}(\lambda)-1\right) / \lambda \\$

此時有，

$\bar{P}_{r-1}(0) = 0$

，

是它的一個根。而且，

$\bar{P}_{r-1}(\lambda)$

是一個

次多項式。如果我們把這個多項式帶入到式（5。32）中，可以得到

$\begin{aligned}\left\|x_r-x^{*}\right\|_{A}^{2} &\le\min _{P_{r-1}} \max _{1 \leq i \leq n}\left[1+\lambda_{i} P_{r-1}\left(\lambda_{i}\right)\right]^{2} \\ &\le \max _{1 \leq i \leq n}\left[1+\lambda_{i} \bar{P}_{r-1}\left(\lambda_{i}\right)\right]^{2} \\ &= \max _{1 \leq i \leq n} Q_r\left(\lambda_i \right) \\ &= 0 \end{aligned} \\$

其中，第二行用特定的

階多項式

$\bar{P}_{r-1} \left(\lambda_i\right)$

代替了任意的

階多項式

$P_{r-1} \left(\lambda_i\right)$

，相當於放大了表示式的值。第四行，根據

$Q_r\left(\lambda\right)$

的定義，輸入任意的

$\lambda_i$

，

$Q_r\left(\lambda_i\right)$

都等於0。因此，我們最終發現，在第

次迭代後，

就已經收斂到了極小值

。

這是個非常好的性質，它說明，共軛梯度法在N次迭代後收斂其實是最壞情況下的結果，當特徵值存在重疊時，只需要更少的迭代次數就能收斂。

不過，以上還沒有回答我們想要了解的問題，在迭代過程中，演算法的收斂速度如何。下面的式子在一定程度上給出了回答（這裡仍然省略證明）

$\left\|x_{k+1}-x^{*}\right\|_{A}^{2} \leq\left(\frac{\lambda_{n-k}-\lambda_{1}}{\lambda_{n-k}+\lambda_{1}}\right)^{2}\left\|x_{0}-x^{*}\right\|_{A}^{2} \tag{5.34} \\$

在第

次迭代後，殘差不大於

$\left(\frac{\lambda_{n-k}-\lambda_{1}}{\lambda_{n-k}+\lambda_{1}}\right)^{2}$

。這個結論也是很有價值的，這裡的

$\lambda_i$

是從小到大排列的特徵值，

$\lambda_{n-k}$

越小，殘差就越小。對於同樣的

次迭代，如果矩陣

的第

大的特徵值比矩陣

的第

大的特徵值小，矩陣

就會收斂得更快。另一方面，也可以理解為，整個收斂過程中的收斂速度是變化的，與特徵值的大小分佈有關。如果特徵值的取值呈聚類狀分佈，比如有5個較大的特徵值，還有幾個接近於1的特徵值，那麼迭代過程很可能像下圖這樣

在7次迭代後就達到了比較小的值。這是因為除了5個較大的特徵值，其它特徵值都集中在一個聚類中。根據前面得出的結論，如果

只有

個不同大小的特徵值，

次迭代後就能收斂。現在的情況與之類似，雖然聚類中的特徵值並非完全一樣，但我們仍然可以將其當做近似相等的特徵值，所以最佳化在6次迭代後理論上就會接近於極值。當然，從圖中可以看到，實際情況下，第7次迭代才接近極值，這也是合理的。

如果我們進一步對式（5。34）放縮，可以得到一個更寬鬆的收斂性質

$\left\|x_{k}-x^{*}\right\|_{A} \leq 2\left(\frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa(A)}+1}\right)^{k}\left\|x_{0}-x^{*}\right\|_{A} \tag{5.35} \\$

在該式中，收斂速度受

的條件數制約，條件數越小，收斂速度越快。在矩陣的條件數中我們介紹了條件數的定義和計算方式，即

$\kappa(A)=\|A\|_{2}\left\|A^{-1}\right\|_{2}=\lambda_{n} / \lambda_{1}$

。這說明，特徵值分佈越集中，收斂速度越快，分佈越分散，收斂速度越慢。極端情況下，當

時，所有特徵值都相等，

的條件數為1，此時只需要一次迭代即可收斂到極值點。

既然共軛梯度法的收斂速度如此依賴於特徵值分佈，我們理應做點什麼，以改善那些特徵值分佈不佳的情況。

預處理策略

最容易想到，也是最有效的方法就是，對

做個變換，讓變換後的

$\hat{A}$

具有更為集中的特徵值。比如，令

$\hat{x} = Cx \tag{5.36} \\$

針對

$\hat{x}$

的最佳化目標函式為

$\hat{\phi}(\hat{x})=\frac{1}{2} \hat{x}^{T}\left(C^{-T} A C^{-1}\right) \hat{x}-\left(C^{-T} b\right)^{T} \hat{x} \tag{5.37} \\$

其中，新最佳化方程中的係數矩陣為

$\hat{A} = C^{-T}AC^{-1}$

。如果我們找到合適的

，使得

$\hat{A}$

具有比較集中的特徵值以及較小的條件數，那麼對

$\hat{\phi}(\hat{x})$

的最佳化就能很快收斂。

但問題的關鍵在於如何選擇合適的

。極端情況下，如果令

$\hat{A} = I$

，可以發現，

恰好是

的Cholesky分解，即

。但我們不能這樣做，因為對

做Cholesky分解無異於直接求

的逆，計算量太大。一種更快的方案是Incomplete Cholesky分解，該分解得到使得

$A \approx \tilde{C}^T\tilde{C}$

的稀疏矩陣

$\tilde{C}$

，該矩陣可保證

$\tilde{C}^{-T} A \tilde{C}^{-1} \approx I$

，從而得到良好的特徵值分佈。

考慮預處理之後，式（5。23）轉化為如下演算法

$\begin{aligned} \alpha_{k} & \leftarrow \frac{r_{k}^{T} y_{k}}{p_{k}^{T} A p_{k}} \\ x_{k+1} & \leftarrow x_{k}+\alpha_{k} p_{k} \\ r_{k+1} & \leftarrow r_{k}+\alpha_{k}A p_{k} \\ M y_{k+1} & \leftarrow r_{k+1} \\ \beta_{k+1} & \leftarrow-\frac{r_{k+1}^{T} y_{k+1}}{r_{k}^{T} y_{k}} \\ p_{k+1} & \leftarrow-y_{k+1}+\beta_{k+1} p_{k} \\ k & \leftarrow k+1 \end{aligned}\tag{5.38} \\$

從

開始，

$r_0 = A x_0 - b, \space y_0 = M^{-1}r_0, \space p_0 = -y_0$

，依次迭代下去。不過，預處理雖然改善了收斂速度，卻多了求解

$y_k = M^{-1} r_k$

的步驟，會增加一部分計算量。

總結

本文介紹了線性共軛梯度法，這是一種以迭代的方式求解線性方程的方法。這種方法的妙處在於，不需要矩陣分解和求逆，只需要上次迭代的步長和當前的梯度，即可計算下次的步長。計算量極低，且不佔用額外的記憶體空間，適用於求解大規模的線性方程。

參考資料

Chapter 16 Preconditioning Illinois math course

上一篇：信賴域（二）：確切方法

下一篇：共軛梯度法（二）：非線性共軛梯度

小蜜蜂問答

小蜜蜂問答

共軛梯度法（一）：線性共軛梯度

推薦文章

小蜜蜂問答

小蜜蜂問答

共軛梯度法（一）：線性共軛梯度

相關文章

求最大無關組（線性代數問題，要過程，急！）？

Arzela-Ascoli 定理的一個小小的推廣

牛頓法和擬牛頓法——（書中附錄B）

梯度下降和牛頓法、擬牛頓法

推薦文章