Krylov 子空間法──線性方程的數值解法 (三):共軛梯度法

本文的閱讀等級:高級

共軛梯度法 (conjugate gradient method) 是一個適用於實對稱正定矩陣的線性方程數值解法。顧名思義,共軛梯度法的核心是共軛 (conjugacy) 和梯度 (一階導數)。共軛能夠加快收斂,梯度則提供正交基底。因為這兩個特性,共軛梯度法的結構簡單優美,儲存量及運算量少,並且無須設定參數。對於大尺寸矩陣,我們往往無法使用直接法求解,譬如 Cholesky 分解 (見“Cholesky 分解”),這時候可以採用以迭代方式計算的共軛梯度法。此外,對於大型非線性最佳化問題,共軛梯度法也是最有效的數值算法之一。

 
A 為一個 n\times n 階實對稱矩陣,\mathbf{b}\in\mathbb{R}^n 為非零向量。假設 \{\mathbf{p}_1,\ldots,\mathbf{p}_n\}\mathbb{R}^n 的一組基底。令 P=\begin{bmatrix}  \mathbf{p}_1&\cdots&\mathbf{p}_n  \end{bmatrix}。明顯地,Pn\times n 階可逆矩陣。線性方程 A\mathbf{x}=\mathbf{b} 的精確解 \mathbf{x}_\ast=A^{-1}\mathbf{b} 可表示為

\displaystyle  \mathbf{x}_\ast=\alpha_1\mathbf{p}_1+\cdots+\alpha_n\mathbf{p}_n=P\boldsymbol{\alpha}

其中 \boldsymbol{\alpha}=(\alpha_1,\ldots,\alpha_n)^T。線性方程 A\mathbf{x}=\mathbf{b} 等號兩邊左乘 P^T,並代入上式可得等價式

P^TAP\boldsymbol{\alpha}=P^T\mathbf{b}

其中係數矩陣 P^TAP 仍為對稱矩陣,\boldsymbol{\alpha} 是待解的向量。若 P^TAP 是對角矩陣,(P^TAP)_{ij}=\mathbf{p}_i^TA\mathbf{p}_j=0i\neq j,則

\displaystyle  \begin{bmatrix}  \mathbf{p}_1^TA\mathbf{p}_1&&\\  &\ddots&\\  &&\mathbf{p}_n^TA\mathbf{p}_n  \end{bmatrix}\begin{bmatrix}  \alpha_1\\  \vdots\\  \alpha_n  \end{bmatrix}=\begin{bmatrix}  \mathbf{p}_1^T\mathbf{b}\\  \vdots\\  \mathbf{p}_n^T\mathbf{b}  \end{bmatrix}

因為 \mathbf{p}_i\neq\mathbf{0},若 A 是正定矩陣,則 \mathbf{p}_i^TA\mathbf{p}_i>0 (見“特殊矩陣 (6):正定矩陣”)。換句話說,正定矩陣保證可逆,故

\displaystyle  \alpha_i=\frac{\mathbf{p}_i^T\mathbf{b}}{\mathbf{p}_i^TA\mathbf{p}_i},~~i=1,\ldots,n

精確解可表示為

\displaystyle  \mathbf{x}_\ast=\sum_{i=1}^n\frac{\mathbf{p}_i^T\mathbf{b}}{\mathbf{p}_i^TA\mathbf{p}_i}\mathbf{p}_i

前述推演引申出下列定義:對於一個實對稱矩陣 A,若 \mathbf{p}_i^TA\mathbf{p}_j=0,我們說 \mathbf{p}_i\mathbf{p}_j 對於 A 共軛 (conjugate) 或 A─正交 (A-orthogonal)。本文考慮 A 是正定矩陣,但此性質並不納入共軛的基本定義。若 A=0,則任何兩個向量皆為共軛。若 A=I,則共軛等價於正交。我們稱 \{\mathbf{p}_1,\ldots,\mathbf{p}_k\} 對於 A 為共軛向量集,若 \mathbf{p}_i^TA\mathbf{p}_j=0i\neq j (在不造成混淆的情況下,以下簡稱共軛向量集)。對於實對稱正定矩陣 A,若 \{\mathbf{p}_1,\ldots,\mathbf{p}_k\} 為非零共軛向量集,則它們是線性獨立向量。證明於下:假設數組 c_1,\ldots,c_k 使得 c_1\mathbf{p}_1+\cdots+c_k\mathbf{p}_k=\mathbf{0}。等號兩邊左乘 \mathbf{p}_j^TA\sum_{i=1}^kc_i\mathbf{p}_j^TA\mathbf{p}_i=c_j\mathbf{p}_j^TA\mathbf{p}_j=0。因為 \mathbf{p}_j^TA\mathbf{p}_j>0,證明 c_j=0。在向量空間 \mathbb{R}^n,一旦獲得完整的非零共軛向量集 \{\mathbf{p}_1,\ldots,\mathbf{p}_n\},透過矩陣─向量乘法與內積運算即可求出精確解。

 
如欲採用迭代法計算近似解,我們必須改寫精確解的表達式。給定初始猜測解 \mathbf{x}_0,以及非零共軛向量集 \{\mathbf{p}_0,\mathbf{p}_1,\ldots,\mathbf{p}_{n-1}\},令

\displaystyle  \mathbf{x}_\ast-\mathbf{x}_0=\alpha_0\mathbf{p}_0+\alpha_1\mathbf{p}_1+\cdots+\alpha_{n-1}\mathbf{p}_{n-1}

\mathbf{x}_j-\mathbf{x}_0=\alpha_0\mathbf{p}_0+\alpha_1\mathbf{p}_1+\cdots+\alpha_{j-1}\mathbf{p}_{j-1}j>0

因此,\mathbf{x}_n=\mathbf{x}_\ast。下列遞迴式成立:

\mathbf{x}_{j+1}=\mathbf{x}_j+\alpha_j\mathbf{p}_j,~~0\le j\le n-1

利用共軛性很容易解出 \alpha_j。考慮

\displaystyle  \mathbf{p}_j^TA(\mathbf{x}_\ast-\mathbf{x}_j)=\mathbf{p}_j^TA(\alpha_j\mathbf{p}_j+\cdots+\alpha_{n-1}\mathbf{p}_{n-1})=\alpha_j\mathbf{p}_j^TA\mathbf{p}_j

即得

\displaystyle  \alpha_j=\frac{\mathbf{p}_j^TA(\mathbf{x}_\ast-\mathbf{x}_j)}{\mathbf{p}_j^TA\mathbf{p}_j}=\frac{\mathbf{p}_j^T(\mathbf{b}-A\mathbf{x}_j)}{\mathbf{p}_j^TA\mathbf{p}_j}=\frac{\mathbf{p}_j^T\mathbf{r}_j}{\mathbf{p}_j^TA\mathbf{p}_j}

其中 \mathbf{r}_j=\mathbf{b}-A\mathbf{x}_j,稱為殘差 (residual)。

 
下面這個性質說明共軛向量集的價值所在。令 \mathcal{P}_j=\hbox{span}\{\mathbf{p}_0,\mathbf{p}_1,\ldots,\mathbf{p}_{j-1}\}。考慮二次型

\displaystyle  E(\mathbf{x})=\frac{1}{2}(\mathbf{x}_\ast-\mathbf{x})^TA(\mathbf{x}_\ast-\mathbf{x})

\mathbf{x}=\mathbf{x}_\astE(\mathbf{x}) 有極小值 0。限定 \mathbf{x}-\mathbf{x}_0\in\mathcal{P}_j,我們可以證明遞迴式產生的 \mathbf{x}_j\mathbf{x}_j-\mathbf{x}_0\in\mathcal{P}_j,最小化 E(\mathbf{x})。令 n\times j 階矩陣 P_j=\begin{bmatrix}  \mathbf{p}_0&\mathbf{p}_1&\cdots&\mathbf{p}_{j-1}  \end{bmatrix}。設 \mathbf{x}=\mathbf{x}_0+P_j\mathbf{y},上面的約束優化問題等價於尋找 \mathbf{y}\in\mathbb{R}^j 使最小化

\displaystyle  E(\mathbf{y})=\frac{1}{2}(\mathbf{x}_\ast-(\mathbf{x}_0+P_j\mathbf{y}))^TA(\mathbf{x}_\ast-(\mathbf{x}_0+P_j\mathbf{y}))

因為 A^T=A,求導可得 (見“矩陣導數”,SV-10)

\displaystyle  \frac{\partial E}{\partial\mathbf{y}}=-P_j^TA(\mathbf{x}_\ast-(\mathbf{x}_0+P_j\mathbf{y}))=-P_j^T(\mathbf{b}-A(\mathbf{x}_0+P_j\mathbf{y}))

最小化 E(\mathbf{y}) 的充要條件為 P_j^T(\mathbf{b}-A(\mathbf{x}_0+P_j\mathbf{y}))=P_j^T(\mathbf{b}-A\mathbf{x})=0,即 \mathbf{b}-A\mathbf{x} 屬於 P_j^T 的零空間 (nullspace) N(P_j^T),也就是說,\mathbf{b}-A\mathbf{x}\perp\mathcal{P}_j,因為 N(P_j^T)=\mathcal{P}_j^\perp (見“線性代數基本定理 (二)”)。我們用歸納法來證明 \mathbf{r}_j\perp\mathcal{P}_j (即 \mathbf{b}-A\mathbf{x}_j\perp\mathcal{P}_j),0\le j\le n-1。若 j=0\mathcal{P}_0 是空集合,命題自然成立。假設 \mathbf{r}_j\perp\mathcal{P}_j。寫出 \mathbf{r}_j 的遞迴式,

\mathbf{r}_{j+1}=\mathbf{b}-A\mathbf{x}_{j+1}=\mathbf{b}-A(\mathbf{x}_j+\alpha_j\mathbf{p}_j)=\mathbf{r}_j-\alpha_jA\mathbf{p}_j

等號兩邊左乘 \mathbf{p}_j^T,再套用 \alpha_j 的公式,

\mathbf{p}_j^T\mathbf{r}_{j+1}=\mathbf{p}_j^T\mathbf{r}_j-\alpha_j\mathbf{p}_j^TA\mathbf{p}_j=0

對於 i<j\mathbf{p}_i\in\mathcal{P}_j,使用先前假設 \mathbf{r}_j\perp\mathbf{p}_i,以及共軛性,可得

\mathbf{p}_i^T\mathbf{r}_{j+1}=\mathbf{p}_i^T\mathbf{r}_j-\alpha_j\mathbf{p}_i^TA\mathbf{p}_j=0

合併以上結果,證明 \mathbf{r}_{j+1}\perp\mathcal{P}_{j+1}

 
因為 \mathcal{P}_1\subset\mathcal{P}_2\subset\cdots\subset\mathcal{P}_n=\mathbb{R}^n,由上述性質可推論 E(\mathbf{x}_0)\ge E(\mathbf{x}_1)\ge\cdots\ge E(\mathbf{x}_n)=0。理論上,使用至多 n 次計算可得精確解 \mathbf{x}_\ast。接下來的工作是找出適當的非零共軛向量集 \{\mathbf{p}_0,\mathbf{p}_1,\ldots,\mathbf{p}_{j-1}\} 使得子空間 \mathcal{P}_j 包含較接近 \mathbf{x}_\ast 的近似解。最佳化理論提供了一些線索。將目標函數改寫如下:

\displaystyle  \begin{aligned}  E(\mathbf{x})&=\frac{1}{2}(\mathbf{x}_\ast-\mathbf{x})^TA(\mathbf{x}_\ast-\mathbf{x})\\  &=\frac{1}{2}\mathbf{x}^TA\mathbf{x}-\mathbf{x}^TA\mathbf{x}_\ast+\frac{1}{2}\mathbf{x}_\ast^TA\mathbf{x}_\ast\\  &=\frac{1}{2}\mathbf{x}^TA\mathbf{x}-\mathbf{x}^T\mathbf{b}+\frac{1}{2}\mathbf{x}_\ast^TA\mathbf{x}_\ast.  \end{aligned}

最小化 E(\mathbf{x}) 等價於最小化

\displaystyle  f(\mathbf{x})=\frac{1}{2}\mathbf{x}^TA\mathbf{x}-\mathbf{x}^T\mathbf{b}

求導可得

\displaystyle  \nabla f=\frac{\partial f}{\partial \mathbf{x}}=A\mathbf{x}-\mathbf{b}

故線性方程 A\mathbf{x}=\mathbf{b} 的唯一解等於二次優化問題 \min_{\mathbf{x}} f(\mathbf{x}) 的唯一解。給定初始猜測解 \mathbf{x}_0,梯度下降法 (gradient descent) 使用下列遞迴公式更新近似解 (見“最佳化理論與正定矩陣”):

\displaystyle  \mathbf{x}_{j+1}=\mathbf{x}_j-\eta_j\nabla f(\mathbf{x}_j)=\mathbf{x}_j+\eta_j\mathbf{r}_j,~~j=0,1,2,\ldots

其中 -\nabla f(\mathbf{x}_j)=\mathbf{b}-A\mathbf{x}_j=\mathbf{r}_j 表示朝精確解的前進方向,\eta_j 是一個夠小的正數,稱為步長 (step size),使得 f(\mathbf{x}_0)\ge f(\mathbf{x}_1)\ge f(\mathbf{x}_2)\ge\cdots。若 \mathbf{r}_j=\mathbf{0},則 \mathbf{x}_j 即為線性方程的解。我們可以設 \eta_j=\arg\min_{\eta>0}f(\mathbf{x}_j+\eta\mathbf{r}_j)[1],這個方法稱為最陡下降法 (steepest descent)。然而,當 A 的條件數 (condition number) 大時,\mathbf{r}_j 未必指向靠近精確解的方向,最陡下降法的收斂速度因此變得很慢 (見“條件數”)。為了克服這個問題,共軛梯度法提出一個非常精巧的想法:以共軛向量 \mathbf{p}_1,\mathbf{p}_2,\ldots 取代殘差 (負梯度) \mathbf{r}_1,\mathbf{r}_2,\ldots,這麼做的原因是限定 \mathbf{x}-\mathbf{x}_0\in\mathcal{P}_j,遞迴式 \mathbf{x}_j=\mathbf{x}_{j-1}+\alpha_{j-1}\mathbf{p}_{j-1} 得到的 \mathbf{x}_j 可使 f(\mathbf{x}) 有最小值。

 
下面介紹共軛梯度法。給定初始猜測解 \mathbf{x}_0,令 \mathbf{p}_0=\mathbf{r}_0=\mathbf{b}-A\mathbf{x}_0。重抄近似解和殘差的遞迴公式:

\begin{aligned}  \mathbf{x}_{j+1}&=\mathbf{x}_j+\alpha_j\mathbf{p}_j\\  \mathbf{r}_{j+1}&=\mathbf{r}_j-\alpha_jA\mathbf{p}_j,\end{aligned}

其中 \alpha_j={\mathbf{p}_j^T\mathbf{r}_j}/{\mathbf{p}_j^TA\mathbf{p}_j}。設想共軛向量 \mathbf{p}_{j+1} 可表示為 \mathbf{r}_0,\mathbf{r}_1,\ldots,\mathbf{r}_{j+1} 的線性組合,但為了減少計算,我們採用下列遞迴算式:

\mathbf{p}_{j+1}=\mathbf{r}_{j+1}+\beta_j\mathbf{p}_j

其中 \beta_j 稱為動量 (momentum),滿足共軛性 0=\mathbf{p}_{j+1}^TA\mathbf{p}_j=\mathbf{r}_{j+1}^TA\mathbf{p}_j+\beta_j\mathbf{p}_j^TA\mathbf{p}_j,由此解出 \beta_j=-{\mathbf{r}_{j+1}^TA\mathbf{p}_j}/{\mathbf{p}_j^TA\mathbf{p}_j}。如果 \beta_j=0,則共軛梯度法退化為最陡下降法。

 
定義 Krylov 子空間

\displaystyle  \mathcal{K}_j=\hbox{span}\{\mathbf{r}_0,A\mathbf{r}_0,\ldots,A^{j-1}\mathbf{r}_0\}

A(\mathcal{K}_j)=\{A\mathbf{y}\vert\mathbf{y}\in\mathcal{K}_j\}。若共軛梯度法於 \mathbf{x}_{j-1} 尚未停止,下列性質成立:

  1. \hbox{span}\{\mathbf{r}_0,\mathbf{r}_1,\ldots,\mathbf{r}_{j-1}\}=\mathcal{K}_j
  2. \mathcal{P}_j=\hbox{span}\{\mathbf{p}_0,\mathbf{p}_1,\ldots,\mathbf{p}_{j-1}\}=\mathcal{K}_j
  3. \mathbf{p}_i^TA\mathbf{p}_{j-1}=0i<j-1

根據性質(1)和(2),\mathbf{r}_j\perp\mathcal{P}_j 可推得 \mathbf{r}_j\perp\mathbf{r}_ii<j,也就是說,殘差 \{\mathbf{r}_0,\mathbf{r}_1,\ldots,\mathbf{r}_{j-1}\} 為正交集。性質(3)表明 \{\mathbf{p}_0,\mathbf{p}_1,\ldots,\mathbf{p}_{j-1}\} 為共軛向量集。

 
我們用歸納法同時證明這三個性質。當 j=1,性質(1-3)顯然成立。假設對於 j,性質(1-3)成立。考慮 \mathbf{r}_j=\mathbf{r}_{j-1}+\alpha_{j-1}A\mathbf{p}_{j-1}。根據假設,性質(1)給出 \mathbf{r}_{j-1}\in\mathcal{K}_j\subseteq\mathcal{K}_{j+1},由性質(2)可知 A\mathbf{p}_{j-1}\in A(\mathcal{K}_j)\subseteq\mathcal{K}_{j+1},因此 \mathbf{r}_j\in\mathcal{K}_{j+1}。此外,因為 \mathbf{r}_j\perp\mathcal{P}_j,即 \mathbf{r}_j\perp\mathcal{K}_j,推論 \mathbf{r}_j\notin\mathcal{K}_j,否則 \mathbf{r}_{j}=\mathbf{0}。所以,\hbox{span}\{\mathbf{r}_0,\mathbf{r}_1,\ldots,\mathbf{r}_j\}=\mathcal{K}_{j+1}

接著證明性質(2),考慮 \mathbf{p}_j=\mathbf{r}_j+\beta_{j-1}\mathbf{p}_{j-1}。性質(1)表明 \mathbf{r}_j\in\mathcal{K}_{j+1},性質(2)假設 \mathbf{p}_{j-1}\in\mathcal{K}_j\subseteq\mathcal{K}_{j+1},故 \mathbf{p}_j\in\mathcal{K}_{j+1},也就證明 \mathcal{P}_{j+1}=\mathcal{K}_{j+1}

欲證明性質(3),考慮 \mathbf{p}_j^TA\mathbf{p}_i=\mathbf{r}_j^TA\mathbf{p}_i+\beta_{j-1}\mathbf{p}_{j-1}^TA\mathbf{p}_i。若 i=j-1,由 \beta_{j-1} 的公式可知 \mathbf{p}_j^TA\mathbf{p}_{j-1}=0。若 i<j-1\mathbf{r}_j^TA\mathbf{p}_{i}=0,因為 \mathbf{r}_j\perp\mathcal{K}_jA\mathbf{p}_i\in\mathcal{K}_{i+2}\subseteq\mathcal{K}_j\mathbf{p}_{j-1}^TA\mathbf{p}_i=0,源自性質(3)的假設。所以,\mathbf{p}_i^TA\mathbf{p}_j=0i<j

 
為減少計算量,利用以上性質化簡 \alpha_j\beta_j。因為 \mathbf{r}_j\perp\mathcal{P}_j,即 \mathbf{r}_j\perp\mathbf{p}_{i}i<j,可得 \mathbf{p}_j^T\mathbf{r}_j=(\mathbf{r}_j+\beta_{j-1}\mathbf{p}_{j-1})^T\mathbf{r}_j=\mathbf{r}_j^T\mathbf{r}_j,故

\displaystyle  \alpha_j=\frac{\mathbf{r}_j^T\mathbf{r}_j}{\mathbf{p}_j^TA\mathbf{p}_j}

因為 \mathbf{p}_0=\mathbf{r}_0,推知 \alpha_0=\eta_0 (見附註[1]),共軛梯度法的第一次調整量與最陡下降法相同。使用 A\mathbf{p}_j=-\frac{1}{\alpha_j}(\mathbf{r}_{j+1}-\mathbf{r}_j)\beta_j 的分子為 \mathbf{r}_{j+1}^TA\mathbf{p}_j=-\frac{1}{\alpha_j}\mathbf{r}_{j+1}^T(\mathbf{r}_{j+1}-\mathbf{r}_j)=-\frac{1}{\alpha_j}\mathbf{r}_{j+1}^T\mathbf{r}_{j+1},分母為 \mathbf{p}_j^TA\mathbf{p}_j=\frac{1}{\alpha_j}\mathbf{r}_j^T\mathbf{r}_j,故

\displaystyle  \beta_j=\frac{\mathbf{r}_{j+1}^T\mathbf{r}_{j+1}}{\mathbf{r}_j^T\mathbf{r}_j}

 
我們將共軛梯度法整理於下:對於一 n\times n 階實對稱正定矩陣 A 和非零向量 \mathbf{b}\in\mathbb{R}^n,設定初始解 \mathbf{x}_0 (一般設為零)。

共軛梯度法


  1. \mathbf{r}_0\leftarrow\mathbf{b}-A\mathbf{x}_0\mathbf{p}_0\leftarrow\mathbf{r}_0
  2. 對於 j=0,1,\ldots,n-1
  3.     \alpha_j\leftarrow {\mathbf{r}_j^T\mathbf{r}_j}/{\mathbf{p}_j^TA\mathbf{p}_j}
  4.     \mathbf{x}_{j+1}\leftarrow\mathbf{x}_j+\alpha_j\mathbf{p}_j
  5.     \mathbf{r}_{j+1}\leftarrow\mathbf{r}_j-\alpha_jA\mathbf{p}_j
  6.     若 \Vert\mathbf{r}_{j+1}\Vert/\Vert\mathbf{b}\Vert<\epsilon,停止計算。
  7.     \beta_j\leftarrow\mathbf{r}_{j+1}^T\mathbf{r}_{j+1}/\mathbf{r}_j^T\mathbf{r}_j
  8.     \mathbf{p}_{j+1}\leftarrow \mathbf{r}_{j+1}+\beta_j\mathbf{p}_j

若不考慮數值計算引進的捨入誤差,共軛梯度法於第 k 次迭代後得到精確解,k\le n。為控制收斂,設相對殘量 \Vert\mathbf{r}_{j+1}\Vert/\Vert\mathbf{b}\Vert 的容許誤差為 \epsilon。一般而言,共軛梯度法可快速逼近至精確解,收斂速度取決於 A 的條件數。如果 A 的條件數過大 (即接近病態),通過前處理 (preconditioning) 可改善收斂性能。具體地說,將 A\mathbf{x}=\mathbf{b} 以等價的 M^{-1}A\mathbf{x}=M^{-1}\mathbf{b} 取代,其中 M 是實對稱正定矩陣使得 M^{-1}A 有較小的條件數[2]。下圖顯示共軛梯度法 (紅色線段) 與最陡下降法 (綠色線段) 的收斂軌跡,兩個算法產生相同的 \mathbf{x}_1,但共軛梯度法迭代二次 (n=2) 便得到解。

Conjugate_gradient_illustration.svg

共軛梯度法與最陡下降法的收斂軌跡 From Wikimedia

 
共軛梯度法具有兩個標誌性質:(1) 朝向精確解的推進方向有共軛性,\mathbf{p}_i^TA\mathbf{p}_j=0i\neq j,(2) 殘差有正交性,\mathbf{r}_i^T\mathbf{r}_j=0i\neq j。最後我們說明何以共軛梯度法被歸於 Krylov 子空間法。共軛梯度法設定近似解 \mathbf{x}_j=\mathbf{x}_0+\alpha_0\mathbf{p}_0+\alpha_1\mathbf{p}_1+\cdots+\alpha_{j-1}\mathbf{p}_{j-1},可知 \mathbf{x}_j-\mathbf{x}_0\in\mathcal{K}_j,且殘差滿足 \mathbf{r}_j\perp\mathcal{K}_j,此即 Petrov-Galerkin 條件 \mathbf{b}-A\mathbf{x}_j\perp\mathcal{K}_j (見“Krylov 子空間法──線性方程的數值解法 (二):GMRES 與 FOM”)。共軛梯度法與 FOM 有相同的定義條件,差異在於前者的 Krylov 子空間基底為共軛向量集並以遞迴式計算近似解 \mathbf{x}_j,後者則以 Lanczos 算法取代 Arnoldi 算法 (因為 A^T=A) 求得 Krylov 子空間的一組單範正交基底 (orthonormal basis) 之後再計算線性方程解 \mathbf{x}_j-\mathbf{x}_0

 
附註:
[1] 寫出

\displaystyle  f(\mathbf{x}_j+\eta\mathbf{r}_j)=\frac{\eta^2}{2}\mathbf{r}_j^TA\mathbf{r}_j+\eta\mathbf{r}_j^TA\mathbf{x}_j+\frac{1}{2}\mathbf{x}_j^T\mathbf{x}_j-\mathbf{x}_j^T\mathbf{b}-\eta\mathbf{r}_j^T\mathbf{b}

\eta 求導可得

\displaystyle  \frac{df}{d\eta}=\eta\mathbf{r}_j^TA\mathbf{r}_j+\mathbf{r}_j^TA\mathbf{x}_j-\mathbf{r}_j^T\mathbf{b}=\eta\mathbf{r}_j^TA\mathbf{r}_j-\mathbf{r}_j^T\mathbf{r}_j

設上式為零,解出

\displaystyle \eta_j=\frac{\mathbf{r}_j^T\mathbf{r}_j}{\mathbf{r}_j^TA\mathbf{r}_j}

[2] 維基百科:Preconditioner

繼續閱讀:
廣告
本篇發表於 線性代數專欄, 數值線性代數 並標籤為 , , , , 。將永久鏈結加入書籤。

One Response to Krylov 子空間法──線性方程的數值解法 (三):共軛梯度法

  1. TZ 說道:

    感謝, 十分有幫助!

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s