共軛梯度法的推導

在數值線性代數中，共軛梯度法是一種求解對稱正定線性方程組

$$\boldsymbol{Ax}=\boldsymbol{b}$$

迭代方法。共軛梯度法可以從不同的角度推導而得，包括作為求解最優化問題的共軛方向法的特例，以及作為求解特徵值問題的Arnoldi/Lanczos迭代的變種。

本條目記述這些推導方法中的重要步驟。

## 從共軛方向法推導

共軛梯度法可以看作是應用於二次函數最小化的共軛方向法的特例

$$f(\boldsymbol{x})=\boldsymbol{x}^\mathrm{T}\boldsymbol{A}\boldsymbol{x}-2\boldsymbol{b}^\mathrm{T}\boldsymbol{x}\text{.}$$

### 共軛方向法

在共軛方向上求最小化：

$$f(\boldsymbol{x})=\boldsymbol{x}^\mathrm{T}\boldsymbol{A}\boldsymbol{x}-2\boldsymbol{b}^\mathrm{T}\boldsymbol{x}\text{.}$$

從最初的猜測$\boldsymbol{x}_0$開始，以及相應的殘差 $\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}_0$, 並通過公式計算迭代和殘差

$$\begin{align} \alpha_i&=\frac{\boldsymbol{p}_i^\mathrm{T}\boldsymbol{r}_i}{\boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_i}\text{,}\\ \boldsymbol{x}_{i+1}&=\boldsymbol{x}_i+\alpha_i\boldsymbol{p}_i\text{,}\\ \boldsymbol{r}_{i+1}&=\boldsymbol{r}_i-\alpha_i\boldsymbol{Ap}_i \end{align}$$

其中，$\boldsymbol{p}_0,\boldsymbol{p}_1,\boldsymbol{p}_2,\ldots$ 為一系列相互共軛的方向，例如：

$\boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_j=0$，對於任何$i\neq j$

由於共軛方向法沒有給出用於選擇方向的公式，因此該方法具有誤差 $\boldsymbol{p}_0,\boldsymbol{p}_1,\boldsymbol{p}_2,\ldots$

而共軛方向法也有多種方法分支，包括共軛梯度法和高斯消元法。

## 從Arnoldi/Lanczos迭代推導

共軛梯度法可以看作Arnoldi/Lanczos迭代應用於求解線性方程組時的一個變種。

### 一般Arnoldi方法

Arnoldi迭代從一個向量$\boldsymbol{r}_0$開始，通過定義$\boldsymbol{v}_i=\boldsymbol{w}_i/\lVert\boldsymbol{w}_i\rVert_2$，其中

$$\boldsymbol{w}_i=\begin{cases} \boldsymbol{r}_0 & \text{if }i=1\text{,}\\ \boldsymbol{Av}_{i-1}-\sum_{j=1}^{i-1}(\boldsymbol{v}_j^\mathrm{T}\boldsymbol{Av}_{i-1})\boldsymbol{v}_j & \text{if }i>1\text{,} \end{cases}$$

逐步建立Krylov子空間

$$\mathcal{K}(\boldsymbol{A},\boldsymbol{r}_0)=\{\boldsymbol{r}_0,\boldsymbol{Ar}_0,\boldsymbol{A}^2\boldsymbol{r}_0,\ldots\}$$

一組標準正交基$\{\boldsymbol{v}_1,\boldsymbol{v}_2,\boldsymbol{v}_3,\ldots\}$。

換言之，對於$i>1$，$\boldsymbol{v}_i$由將$\boldsymbol{Av}_{i-1}$相對於$\{\boldsymbol{v}_1,\boldsymbol{v}_2,\ldots,\boldsymbol{v}_{i-1}\}$進行Gram-Schmidt正交化然後歸一化得到。

寫成矩陣形式，迭代過程可以表示為

$$\boldsymbol{AV}_i=\boldsymbol{V}_{i+1}\boldsymbol{\tilde{H}}_i\text{,}$$

其中

$$\begin{align} \boldsymbol{V}_i&=\begin{bmatrix} \boldsymbol{v}_1 & \boldsymbol{v}_2 & \cdots & \boldsymbol{v}_i \end{bmatrix}\text{,}\\ \boldsymbol{\tilde{H}}_i&=\begin{bmatrix} h_{11} & h_{12} & h_{13} & \cdots & h_{1,i}\\ h_{21} & h_{22} & h_{23} & \cdots & h_{2,i}\\ & h_{32} & h_{33} & \cdots & h_{3,i}\\ & & \ddots & \ddots & \vdots\\ & & & h_{i,i-1} & h_{i,i}\\ & & & & h_{i+1,i} \end{bmatrix}=\begin{bmatrix} \boldsymbol{H}_i\\ h_{i+1,i}\boldsymbol{e}_i^\mathrm{T} \end{bmatrix}\text{,}\\ h_{ji}&=\begin{cases} \boldsymbol{v}_j^\mathrm{T}\boldsymbol{Av}_i & \text{if }j\leq i\text{,}\\ \lVert\boldsymbol{w}_{i+1}\rVert_2 & \text{if }j=i+1\text{,}\\ 0 & \text{if }j>i+1\text{.} \end{cases} \end{align}$$

當應用於求解線性方程組時，Arnoldi迭代對應於初始解$\boldsymbol{x}_0$的殘量$\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}_0$開始迭代，在每一步迭代之後計算$\boldsymbol{y}_i=\boldsymbol{H}_i^{-1}(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)$和新的近似解$\boldsymbol{x}_i=\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{y}_i$.

### 直接Lanzcos方法

在餘下的討論中，我們假定$\boldsymbol{A}$是對稱正定矩陣。由於$\boldsymbol{A}$的對稱性, 上Hessenberg矩陣$\boldsymbol{H}_i=\boldsymbol{V}_i^\mathrm{T}\boldsymbol{AV}_i$變成對陣三對角矩陣。於是它可以被更明確地表示為

$$\boldsymbol{H}_i=\begin{bmatrix} a_1 & b_2\\ b_2 & a_2 & b_3\\ & \ddots & \ddots & \ddots\\ & & b_{i-1} & a_{i-1} & b_i\\ & & & b_i & a_i \end{bmatrix}\text{.}$$

這使得迭代中的$\boldsymbol{v}_i$有一個簡短的三項遞推式。Arnoldi迭代從而簡化為Lanczos迭代。

由於$\boldsymbol{A}$對稱正定，$\boldsymbol{H}_i$同樣也對稱正定。因此，$\boldsymbol{H}_i$可以通過不選主元的LU分解分解為

$$\boldsymbol{H}_i=\boldsymbol{L}_i\boldsymbol{U}_i=\begin{bmatrix} 1\\ c_2 & 1\\ & \ddots & \ddots\\ & & c_{i-1} & 1\\ & & & c_i & 1 \end{bmatrix}\begin{bmatrix} d_1 & b_2\\ & d_2 & b_3\\ & & \ddots & \ddots\\ & & & d_{i-1} & b_i\\ & & & & d_i \end{bmatrix}\text{,}$$

其中$c_i$和$d_i$有簡單的遞推式：

$$\begin{align} c_i&=b_i/d_{i-1}\text{,}\\ d_i&=\begin{cases} a_1 & \text{if }i=1\text{,}\\ a_i-c_ib_i & \text{if }i>1\text{.} \end{cases} \end{align}$$

改寫$\boldsymbol{x}_i=\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{y}_i$為

$$\begin{align} \boldsymbol{x}_i&=\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{H}_i^{-1}(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)\\ &=\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{U}_i^{-1}\boldsymbol{L}_i^{-1}(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)\\ &=\boldsymbol{x}_0+\boldsymbol{P}_i\boldsymbol{z}_i\text{,} \end{align}$$

其中

$$\begin{align} \boldsymbol{P}_i&=\boldsymbol{V}_{i}\boldsymbol{U}_i^{-1}\text{,}\\ \boldsymbol{z}_i&=\boldsymbol{L}_i^{-1}(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)\text{.} \end{align}$$

此時需要觀察到

$$\begin{align} \boldsymbol{P}_i&=\begin{bmatrix} \boldsymbol{P}_{i-1} & \boldsymbol{p}_i \end{bmatrix}\text{,}\\ \boldsymbol{z}_i&=\begin{bmatrix} \boldsymbol{z}_{i-1}\\ \zeta_i \end{bmatrix}\text{.} \end{align}$$

實際上，$\boldsymbol{p}_i$和$\zeta_i$同樣有簡短的遞推式：

$$\begin{align} \boldsymbol{p}_i&=\frac{1}{d_i}(\boldsymbol{v}_i-b_i\boldsymbol{p}_{i-1})\text{,}\\ \alpha_i&=-c_i\zeta_{i-1}\text{.} \end{align}$$

通過這個形式，我們得到$\boldsymbol{x}_i$的一個簡單的遞推式：

$$\begin{align} \boldsymbol{x}_i&=\boldsymbol{x}_0+\boldsymbol{P}_i\boldsymbol{z}_i\\ &=\boldsymbol{x}_0+\boldsymbol{P}_{i-1}\boldsymbol{z}_{i-1}+\zeta_i\boldsymbol{p}_i\\ &=\boldsymbol{x}_{i-1}+\zeta_i\boldsymbol{p}_i\text{.} \end{align}$$

以上的遞推關係立即導出比共軛梯度法稍微更複雜的直接Lanczos方法。

### 從正交性和共軛性導出共軛梯度法

如果我們允許縮放$\boldsymbol{p}_i$並通過常數因子補償縮放的係數，我們可能可以得到以下形式的更簡單的遞推式：

$$\begin{align} \boldsymbol{x}_i&=\boldsymbol{x}_{i-1}+\alpha_{i-1}\boldsymbol{p}_{i-1}\text{,}\\ \boldsymbol{r}_i&=\boldsymbol{r}_{i-1}-\alpha_{i-1}\boldsymbol{Ap}_{i-1}\text{,}\\ \boldsymbol{p}_i&=\boldsymbol{r}_i+\beta_{i-1}\boldsymbol{p}_{i-1}\text{.} \end{align}$$

作為簡化的前提，我們現在推導$\boldsymbol{r}_i$的正交性和$\boldsymbol{p}_i$的共軛性，即對於$i\neq j$,

$$\begin{align} \boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_j&=0\text{,}\\ \boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_j&=0\text{.} \end{align}$$

各個殘量之間滿足正交性的原因是$\boldsymbol{r}_i$實際上可由$\boldsymbol{v}_{i+1}$乘上一個係數而得。這是因為對於$i=0$，$\boldsymbol{r}_0=\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{v}_1$，對於$i>0$。

$$\begin{align} \boldsymbol{r}_i&=\boldsymbol{b}-\boldsymbol{Ax}_i\\ &=\boldsymbol{b}-\boldsymbol{A}(\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{y}_i)\\ &=\boldsymbol{r}_0-\boldsymbol{AV}_i\boldsymbol{y}_i\\ &=\boldsymbol{r}_0-\boldsymbol{V}_{i+1}\boldsymbol{\tilde{H}}_i\boldsymbol{y}_i\\ &=\boldsymbol{r}_0-\boldsymbol{V}_i\boldsymbol{H}_i\boldsymbol{y}_i-h_{i+1,i}(\boldsymbol{e}_i^\mathrm{T}\boldsymbol{y}_i)\boldsymbol{v}_{i+1}\\ &=\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{v}_1-\boldsymbol{V}_i(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)-h_{i+1,i}(\boldsymbol{e}_i^\mathrm{T}\boldsymbol{y}_i)\boldsymbol{v}_{i+1}\\ &=-h_{i+1,i}(\boldsymbol{e}_i^\mathrm{T}\boldsymbol{y}_i)\boldsymbol{v}_{i+1}\text{.}\end{align}$$

要得到$\boldsymbol{p}_i$的共軛性，只需證明$\boldsymbol{P}_i^\mathrm{T}\boldsymbol{AP}_i$是對角矩陣：

$$\begin{align} \boldsymbol{P}_i^\mathrm{T}\boldsymbol{AP}_i&=\boldsymbol{U}_i^{-\mathrm{T}}\boldsymbol{V}_i^\mathrm{T}\boldsymbol{AV}_i\boldsymbol{U}_i^{-1}\\ &=\boldsymbol{U}_i^{-\mathrm{T}}\boldsymbol{H}_i\boldsymbol{U}_i^{-1}\\ &=\boldsymbol{U}_i^{-\mathrm{T}}\boldsymbol{L}_i\boldsymbol{U}_i\boldsymbol{U}_i^{-1}\\ &=\boldsymbol{U}_i^{-\mathrm{T}}\boldsymbol{L}_i \end{align}$$

是對稱的下三角矩陣，因而必然是對角矩陣。

現在我們可以單純由$\boldsymbol{r}_i$的正交性和$\boldsymbol{p}_i$的共軛性推導相對於縮放後的$\boldsymbol{p}_i$的常數因子$\alpha_i$和$\beta_i$。

由於$\boldsymbol{r}_i$的正交性，必然有$\boldsymbol{r}_{i+1}^\mathrm{T}\boldsymbol{r}_i=(\boldsymbol{r}_i-\alpha_i\boldsymbol{Ap}_i)^\mathrm{T}\boldsymbol{r}_i=0$。於是

$$\begin{align} \alpha_i&=\frac{\boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_i}{\boldsymbol{r}_i^\mathrm{T}\boldsymbol{Ap}_i}\\ &=\frac{\boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_i}{(\boldsymbol{p}_i-\beta_{i-1}\boldsymbol{p}_{i-1})^\mathrm{T}\boldsymbol{Ap}_i}\\ &=\frac{\boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_i}{\boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_i}\text{.} \end{align}$$

類似地，由於$\boldsymbol{p}_i$的共軛性，必然有$\boldsymbol{p}_{i+1}^\mathrm{T}\boldsymbol{Ap}_i=(\boldsymbol{r}_{i+1}+\beta_i\boldsymbol{p}_i)^\mathrm{T}\boldsymbol{Ap}_i=0$。於是

$$\begin{align} \beta_i&=-\frac{\boldsymbol{r}_{i+1}^\mathrm{T}\boldsymbol{Ap}_i}{\boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_i}\\ &=-\frac{\boldsymbol{r}_{i+1}^\mathrm{T}(\boldsymbol{r}_i-\boldsymbol{r}_{i+1})}{\alpha_i\boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_i}\\ &=\frac{\boldsymbol{r}_{i+1}^\mathrm{T}\boldsymbol{r}_{i+1}}{\boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_i}\text{.} \end{align}$$

推導至此完成。