從滿正交法 (FOM) 推導共軛梯度法 (CG)

本文的閱讀等級:高級

考慮線性方程 A\mathbf{x}=\mathbf{b},其中 An\times n 階可逆實矩陣,\mathbf{b}\in\mathbb{R}^n 是非零向量。在線性方程的迭代解法中,滿正交法 (full orthogonalization method,簡稱 FOM) 適用於一般非對稱矩陣,共軛梯度法 (conjugate gradient method,簡稱 CG) 則限定於實對稱正定矩陣。令 \mathbf{x}_0 為初始猜測解,\mathbf{r}_j=\mathbf{b}-A\mathbf{x}_j 為對應近似解 \mathbf{x}_j 的殘差 (residual),以及 Krylov 子空間 \mathcal{K}_j=\hbox{span}\{\mathbf{r}_0,A\mathbf{r}_0,\ldots,A^{j-1}\mathbf{r}_0\}。FOM 與 CG 同屬 Krylov 子空間法,而且有同樣的定義條件:

  1. \mathbf{x}_j-\mathbf{x}_0\in\mathcal{K}_j
  2. \mathbf{r}_j\perp\mathcal{K}_j

本文從 FOM 推導 CG,證明 CG 是 FOM 的一種簡約表達。關於 FOM 與 CG 的介紹,請參閱“Krylov 子空間法──線性方程的數值解法 (二):GMRES 與 FOM”和 “Krylov 子空間法──線性方程的數值解法 (三):共軛梯度法”。

 
為方便參照,底下抄錄 FOM 與 CG 算法。對於實對稱矩陣 A,FOM 使用的 Arnoldi 算法替換為對稱 Lanczos 算法 (見“Krylov 子空間法──線性方程的數值解法 (一):Arnoldi 與 Lanczos 算法”)。

FOM


  1. \mathbf{r}_0\leftarrow\mathbf{b}-A\mathbf{x}_0\mathbf{q}_1\leftarrow \mathbf{r}_0/\Vert\mathbf{r}_0\Vert
  2. 對於 j=1,2,\ldots,n
  3.     執行 Lanczos 算法的第 j 次迭代,得到 \tilde{T}_jQ_j
  4.     \mathbf{y}_j\leftarrow T_j^{-1}(\Vert\mathbf{r}_0\Vert\mathbf{e}_1)
  5.     \mathbf{x}_j\leftarrow \mathbf{x}_0+Q_j\mathbf{y}_j
  6.     若 \Vert\mathbf{b}-A\mathbf{x}_j\Vert/\Vert\mathbf{b}\Vert<\epsilon,停止計算。

Lanczos 算法滿足 Lanczos 關係式

AQ_j=Q_jT_j+\mathbf{q}_{j+1}b_j\mathbf{e}_j^T=Q_{j+1}\tilde{T}_j

其中 n\times j 階矩陣 Q_j=\begin{bmatrix}  \mathbf{q}_1&\cdots&\mathbf{q}_j  \end{bmatrix} 使得 Q_j^TQ_j=I_jQ_j^T\mathbf{q}_{j+1}=\mathbf{0} (即 \{\mathbf{q}_1,\ldots,\mathbf{q}_{j+1}\} 是一個單範正交集),\tilde{T}_j=\begin{bmatrix}  T_j\\  b_j\mathbf{e}_j^T  \end{bmatrix}\mathbf{e}_j=(0,\ldots,0,1)^Tj 維單位向量,T_jj\times j 階三對角對稱矩陣,如下:

\displaystyle    T_j=\begin{bmatrix}    a_1&b_1& & & \\    b_1&a_2&b_2&&\\     & b_2&a_3&\ddots&\\    && \ddots&\ddots&b_{j-1}\\    &&&b_{j-1}&a_j  \end{bmatrix}

Lanczos 關係式左乘 Q_j^T,可得 Q_j^TAQ_j=T_j。因為 A 是正定矩陣,在 Lanczos 算法停止前,\hbox{rank}T_j=\hbox{rank}(Q_j^TAQ_j)=j,故 T_j 是可逆矩陣。

 
CG


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

CG 使用三個遞迴式計算近似解 \mathbf{x}_j、推進方向 \mathbf{p}_j 與殘差 \mathbf{r}_j,並具有兩個標誌性質:朝向精確解的推進方向有共軛性,\mathbf{p}_i^TA\mathbf{p}_j=0i\neq j,且殘差有正交性,\mathbf{r}_i^T\mathbf{r}_j=0i\neq j

 
為方便閱讀,我將證明過程分為幾個推導步驟。

 
(1) T_j=L_jU_j

因為 T_j 是三對角對稱矩陣,設 LU 分解為 T_j=L_jU_j (見“LU 分解”),如下:

L_j=\begin{bmatrix}  1&& & & \\    c_1&1&&&\\     & c_2&1&&\\    && \ddots&\ddots&\\    &&&c_{j-1}&1  \end{bmatrix},~~U_j=\begin{bmatrix}    d_1&b_1& & & \\    &d_2&b_2&&\\     & &d_3&\ddots&\\    && &\ddots&b_{j-1}\\    &&&&d_j  \end{bmatrix}

乘開再比較各元,可得 c_i=b_i/d_i1\le i\le j-1d_1=a_1d_k=a_k-b_{k-1}c_{k-1}2\le k\le j

 
(2) \mathbf{y}_j=T_j^{-1}(\Vert\mathbf{r}_0\Vert\mathbf{e}_1)

使用高斯─約當法計算 L_j^{-1} (見“高斯─約當法”),以 j=4 為例,

\begin{aligned}  \begin{bmatrix}  L_4&I_4  \end{bmatrix}&=\begin{bmatrix}  1&&&&\vline&1&&&\\  c_1&1&&&\vline&&1&&\\  &c_2&1&&\vline&&&1&\\  &&c_3&1&\vline&&&&1  \end{bmatrix}\\  &\to\begin{bmatrix}  1&&&&\vline&1&&&\\  &1&&&\vline&-c_1&1&&\\  &&1&&\vline&c_1c_2&-c_2&1&\\  &&&1&\vline&-c_1c_2c_3&c_2c_3&-c_3&1  \end{bmatrix}\\  &=\begin{bmatrix}  I_4&L_4^{-1}\end{bmatrix}.\end{aligned}

推廣至 j\times j 階矩陣,

L_j^{-1}=\begin{bmatrix}  1& & & &\\  -c_1&1 & & &\\  c_1c_2& -c_2&1&&\\  \vdots&\ddots&\ddots&\ddots&\\  (-1)^{j-1}\prod_{i=1}^{j-1}c_i&\cdots&c_{j-2}c_{j-1}&-c_{j-1}&1  \end{bmatrix}

將上式代入 T_j^{-1}=U_j^{-1}L_j^{-1}

\mathbf{y}_j=T_j^{-1}(\Vert\mathbf{r}_0\Vert\mathbf{e}_1)=U_j^{-1}L_j^{-1}(\Vert\mathbf{r}_0\Vert\mathbf{e}_1)=U_j^{-1}\mathbf{z}_j

其中

\mathbf{z}_j=\begin{bmatrix}  \zeta_1\\  \zeta_2\\  \vdots\\  \zeta_j  \end{bmatrix}=\begin{bmatrix}  \mathbf{z}_{j-1}\\  \zeta_j  \end{bmatrix}=L_j^{-1}(\Vert\mathbf{r}_0\Vert\mathbf{e}_1)

\zeta_j=(-1)^{j-1}c_1c_2\cdots c_{j-1}\Vert\mathbf{r}_0\Vert=-c_{j-1}\zeta_{j-1}

 
(3) 近似解 \mathbf{x}_j 的遞迴式

n\times j 階矩陣 \tilde{P}_j=\begin{bmatrix}  \tilde{\mathbf{p}}_1&\cdots&\tilde{\mathbf{p}}_j  \end{bmatrix}=Q_jU_j^{-1}。使用 (2),可得

\displaystyle  \begin{aligned}  \mathbf{x}_j&=\mathbf{x}_0+Q_j\mathbf{y}_j\\  &=\mathbf{x}_0+Q_jU_j^{-1}\mathbf{z}_j\\  &=\mathbf{x}_0+\tilde{P}_j\mathbf{z}_j\\  &=\mathbf{x}_0+\begin{bmatrix}  \tilde{P}_{j-1}&\tilde{\mathbf{p}}_j  \end{bmatrix}  \begin{bmatrix}  \mathbf{z}_{j-1}\\  \zeta_j  \end{bmatrix}\\  &=\mathbf{x}_0+\tilde{P}_{j-1}\mathbf{z}_{j-1}+\zeta_j\tilde{\mathbf{p}}_j\\  &=\mathbf{x}_{j-1}+\zeta_j\tilde{\mathbf{p}}_j.  \end{aligned}

考慮

\displaystyle  \begin{aligned}  \tilde{P}_j^TA\tilde{P}_j&=(Q_jU_j^{-1})^TA(Q_jU_j^{-1})\\  &=(U_j^{-1})^TQ_j^TAQ_jU_j^{-1}\\  &=(U_j^{-1})^TT_jU_j^{-1}\\  &=(U_j^{-1})^TL_jU_jU_j^{-1}\\  &=(U_j^{-1})^TL_j  \end{aligned}

因為 \tilde{P}_j^TA\tilde{P}_j 是對稱矩陣,(U_j^{-1})^TL_j 皆為下三角矩陣,推論 \tilde{P}_j^TA\tilde{P}_j 是對角矩陣,即 \tilde{\mathbf{p}}_i^TA\tilde{\mathbf{p}}_j=0i\neq j,表明 \{\tilde{\mathbf{p}}_1,\ldots,\tilde{\mathbf{p}}_j\} 是共軛向量集。

 
(4) 殘差 \mathbf{r}_j 的遞迴式

使用 (3),

\begin{aligned}  \mathbf{r}_j&=\mathbf{b}-A\mathbf{x}_j\\  &=\mathbf{b}-A(\mathbf{x}_{j-1}+\zeta_j\tilde{\mathbf{p}}_j)\\  &=\mathbf{b}-A\mathbf{x}_{j-1}-\zeta_jA\tilde{\mathbf{p}}_j\\  &=\mathbf{r}_{j-1}-\zeta_jA\tilde{\mathbf{p}}_j.  \end{aligned}

另一方面,使用 Lanczos 關係式,

\begin{aligned}  \mathbf{r}_j&=\mathbf{b}-A(\mathbf{x}_0+Q_j\mathbf{y}_j)\\  &=\mathbf{b}-A\mathbf{x}_0-AQ_j\mathbf{y}_j\\  &=\mathbf{r}_0-(Q_jT_j+\mathbf{q}_{j+1}b_j\mathbf{e}_j^T)\mathbf{y}_j\\  &=\mathbf{r}_0-Q_jT_j\mathbf{y}_j-b_j\mathbf{q}_{j+1}\mathbf{e}_j^T\mathbf{y}_j\\  &=\mathbf{r}_0-Q_j(\Vert\mathbf{r}_0\Vert\mathbf{e}_1)-b_j(\mathbf{e}_j^T\mathbf{y}_j)\mathbf{q}_{j+1}\\  &=\mathbf{r}_0-\Vert\mathbf{r}_0\Vert\mathbf{q}_1-b_j(\mathbf{e}_j^T\mathbf{y}_j)\mathbf{q}_{j+1}\\  &=-b_j(\mathbf{e}_j^T\mathbf{y}_j)\mathbf{q}_{j+1}.  \end{aligned}

寫出 U_j=\begin{bmatrix}  U_{j-1}&b_{j-1}\mathbf{e}_{j-1}\\  \mathbf{0}&d_j  \end{bmatrix}。套用分塊矩陣的逆矩陣公式 (見“分塊矩陣的解題案例──逆矩陣與矩陣乘積的特徵值”),

U_j^{-1}=\begin{bmatrix}  U_{j-1}&b_{j-1}\mathbf{e}_{j-1}\\  \mathbf{0}&d_j  \end{bmatrix}^{-1}=\begin{bmatrix}  U_{j-1}^{-1}&-(b_{j-1}/d_j)U_{j-1}^{-1}\mathbf{e}_{j-1}\\  0&1/d_j\end{bmatrix}

因此,b_j(\mathbf{e}_j^T\mathbf{y}_j)=b_j(\mathbf{e}_j^TU_j^{-1}\mathbf{z}_j)=b_j(\zeta_j/d_j)=c_j\zeta_j=-\zeta_{j+1},推得

\mathbf{r}_j=\zeta_{j+1}\mathbf{q}_{j+1}

證明 \{\mathbf{r}_0,\mathbf{r}_1,\ldots,\mathbf{r}_j\} 是一個正交集。

 
(5) 共軛向量 \tilde{\mathbf{p}}_j 的遞迴式

使用 U_j^{-1} 的分塊矩陣公式,

\displaystyle  \begin{aligned}  \tilde{P}_j&=\begin{bmatrix}  \tilde{P}_{j-1}&\tilde{\mathbf{p}}_j  \end{bmatrix}\\  &=Q_jU_j^{-1}\\  &=\begin{bmatrix}  Q_{j-1}&\mathbf{q}_j  \end{bmatrix}\begin{bmatrix}  U_{j-1}^{-1}&-(b_{j-1}/d_j)U_{j-1}^{-1}\mathbf{e}_{j-1}\\  0&1/d_j\end{bmatrix}\\  &=\begin{bmatrix}  Q_{j-1}U_{j-1}^{-1}&-(b_{j-1}/d_j)Q_{j-1}U_{j-1}^{-1}\mathbf{e}_{j-1}+1/d_j\mathbf{q}_j\end{bmatrix}\\  &=\begin{bmatrix}  \tilde{P}_{j-1}&-(b_{j-1}/d_j)\tilde{P}_{j-1}\mathbf{e}_{j-1}+1/d_j\mathbf{q}_j  \end{bmatrix}\\  &=\begin{bmatrix}  \tilde{P}_{j-1}&(1/d_j)(-b_{j-1}\tilde{\mathbf{p}}_{j-1}+\mathbf{q}_j)  \end{bmatrix},  \end{aligned}

推得

\displaystyle  \begin{aligned}  \tilde{\mathbf{p}}_j&=\frac{1}{d_j}\mathbf{q}_j-\frac{b_{j-1}}{d_{j}}\tilde{\mathbf{p}}_{j-1}\\  &=\frac  {1}{d_j\zeta_j}\mathbf{r}_{j-1}-\frac{b_{j-1}}{d_{j}}\tilde{\mathbf{p}}_{j-1}.\end{aligned}

定義 b_0=0,因此 \tilde{\mathbf{p}}_1=(d_1\zeta_1)^{-1}\mathbf{r}_0

 
(6) 共軛向量 \mathbf{p}_j 的遞迴式

伸縮 \tilde{\mathbf{p}}_j 的長度不影響共軛性,令 \mathbf{p}_j=d_{j+1}\zeta_{j+1}\tilde{\mathbf{p}}_{j+1}。給定初始猜測解 \mathbf{x}_0\mathbf{p}_0=\mathbf{r}_0=\mathbf{b}-A\mathbf{x}_0,步驟 (3-5) 的遞迴式可簡化如下:對於 j\ge 1

\begin{aligned}  \mathbf{x}_{j}&=\mathbf{x}_{j-1}+\alpha_{j-1}\mathbf{p}_{j-1}\\  \mathbf{r}_{j}&=\mathbf{r}_{j-1}-\alpha_{j-1}A\mathbf{p}_{j-1}\\  \mathbf{p}_{j}&=\mathbf{r}_{j}+\beta_{j-1}\mathbf{p}_{j-1},  \end{aligned}

其中 \mathbf{r}_i^T\mathbf{r}_j=0\mathbf{p}_i^TA\mathbf{p}_j=0i\neq j。因此,\mathbf{r}_{j-1}^T\mathbf{r}_{j}=\mathbf{r}_{j-1}^T(\mathbf{r}_{j-1}-\alpha_{j-1}A\mathbf{p}_{j-1})=0,解出

\displaystyle  \begin{aligned}  \alpha_{j-1}&=\frac{\mathbf{r}_{j-1}^T\mathbf{r}_{j-1}}{\mathbf{r}_{j-1}^TA\mathbf{p}_{j-1}}\\  &=\frac{\mathbf{r}_{j-1}^T\mathbf{r}_{j-1}}{(\mathbf{p}_{j-1}-\beta_{j-2}\mathbf{p}_{j-2})^TA\mathbf{p}_{j-1}}\\  &=\frac{\mathbf{r}_{j-1}^T\mathbf{r}_{j-1}}{\mathbf{p}_{j-1}^TA\mathbf{p}_{j-1}}.  \end{aligned}

另外,\mathbf{p}_{j}^TA\mathbf{p}_{j-1}=(\mathbf{r}_{j}+\beta_{j-1}\mathbf{p}_{j-1})^TA\mathbf{p}_{j-1}=0,使用 \alpha_{j-1} 的公式,

\displaystyle  \begin{aligned}  \beta_{j-1}&=-\frac{\mathbf{r}_{j}^TA\mathbf{p}_{j-1}}{\mathbf{p}_{j-1}^TA\mathbf{p}_{j-1}}\\  &=-\frac{\mathbf{r}_{j}^T(\mathbf{r}_{j-1}-\mathbf{r}_{j})}{\alpha_{j-1}\mathbf{p}_{j-1}^TA\mathbf{p}_{j-1}}\\  &=\frac{\mathbf{r}_{j}^T\mathbf{r}_{j}}{\mathbf{r}_{j-1}^T\mathbf{r}_{j-1}}.  \end{aligned}

This entry was posted in 線性代數專欄, 數值線性代數 and tagged , , , . Bookmark the permalink.

發表迴響

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

WordPress.com Logo

你正使用 WordPress.com 帳號留言。 登出 / 變更 )

Twitter picture

你正使用 Twitter 帳號留言。 登出 / 變更 )

Facebook照片

你正使用 Facebook 帳號留言。 登出 / 變更 )

Google+ photo

你正使用 Google+ 帳號留言。 登出 / 變更 )

連結到 %s