查看“︁共轭梯度法的推导”︁的源代码

在[[数值线性代数]]中，[[共轭梯度法]]是一种求解[[对称矩阵|对称]][[正定矩阵|正定]][[线性方程组]]

:<math>\boldsymbol{Ax}=\boldsymbol{b}</math>

的[[迭代方法]]。共轭梯度法可以从不同的角度推导而得，包括作为求解[[最优化]]问题的[[共轭方向法]]的特例，以及作为求解[[特征值]]问题的[[Arnoldi迭代|Arnoldi]]/[[Lanczos迭代|Lanczos]]迭代的变种。

本条目记述这些推导方法中的重要步骤。

==从共轭方向法推导==
{{expand-section|date=2010年6月}}
共轭梯度法可以看作是应用于二次函数最小化的共轭方向法的特例

: <math>f(\boldsymbol{x})=\boldsymbol{x}^\mathrm{T}\boldsymbol{A}\boldsymbol{x}-2\boldsymbol{b}^\mathrm{T}\boldsymbol{x}\text{.}</math>

=== 共轭方向法 ===
在共轭方向上求最小化：

: <math>f(\boldsymbol{x})=\boldsymbol{x}^\mathrm{T}\boldsymbol{A}\boldsymbol{x}-2\boldsymbol{b}^\mathrm{T}\boldsymbol{x}\text{.}</math>

从最初的猜测<math>\boldsymbol{x}_0</math>开始，以及相应的[[误差|残差]] <math>\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}_0</math>, 并通过公式计算[[迭代冪次|迭代]]和残差

: <math>\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}</math>

其中，<math>\boldsymbol{p}_0,\boldsymbol{p}_1,\boldsymbol{p}_2,\ldots</math> 为一系列相互共轭的方向，例如：

: <math>\boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_j=0</math>，对于任何<math>i\neq j</math>

由于共轭方向法没有给出用于选择方向的公式，因此该方法具有误差 <math>\boldsymbol{p}_0,\boldsymbol{p}_1,\boldsymbol{p}_2,\ldots</math>

而共轭方向法也有多种方法分支，包括共轭梯度法和[[高斯消去法|高斯消元法]]。

==从Arnoldi/Lanczos迭代推导==
共轭梯度法可以看作Arnoldi/Lanczos迭代应用于求解线性方程组时的一个变种。

===一般Arnoldi方法===
Arnoldi迭代从一个向量<math>\boldsymbol{r}_0</math>开始，通过定义<math>\boldsymbol{v}_i=\boldsymbol{w}_i/\lVert\boldsymbol{w}_i\rVert_2</math>，其中

:<math>\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}</math>

逐步建立[[Krylov子空间]]

:<math>\mathcal{K}(\boldsymbol{A},\boldsymbol{r}_0)=\{\boldsymbol{r}_0,\boldsymbol{Ar}_0,\boldsymbol{A}^2\boldsymbol{r}_0,\ldots\}</math>

的一组[[标准正交基]]<math>\{\boldsymbol{v}_1,\boldsymbol{v}_2,\boldsymbol{v}_3,\ldots\}</math>。

换言之，对于<math>i>1</math>，<math>\boldsymbol{v}_i</math>由将<math>\boldsymbol{Av}_{i-1}</math>相对于<math>\{\boldsymbol{v}_1,\boldsymbol{v}_2,\ldots,\boldsymbol{v}_{i-1}\}</math>进行[[格拉姆-施密特正交化|Gram-Schmidt正交化]]然后归一化得到。

写成矩阵形式，迭代过程可以表示为

:<math>\boldsymbol{AV}_i=\boldsymbol{V}_{i+1}\boldsymbol{\tilde{H}}_i\text{,}</math>

其中

:<math>\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}</math>

当应用于求解线性方程组时，Arnoldi迭代对应于初始解<math>\boldsymbol{x}_0</math>的残量<math>\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}_0</math>开始迭代，在每一步迭代之后计算<math>\boldsymbol{y}_i=\boldsymbol{H}_i^{-1}(\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{e}_1)</math>和新的近似解<math>\boldsymbol{x}_i=\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{y}_i</math>.

===直接Lanzcos方法===
在余下的讨论中，我们假定<math>\boldsymbol{A}</math>是对称正定矩阵。由于<math>\boldsymbol{A}</math>的对称性, [[Hessenberg矩阵|上Hessenberg矩阵]]<math>\boldsymbol{H}_i=\boldsymbol{V}_i^\mathrm{T}\boldsymbol{AV}_i</math>变成对阵三对角矩阵。于是它可以被更明确地表示为

:<math>\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{.}</math>

这使得迭代中的<math>\boldsymbol{v}_i</math>有一个简短的三项递推式。Arnoldi迭代从而简化为Lanczos迭代。

由于<math>\boldsymbol{A}</math>对称正定，<math>\boldsymbol{H}_i</math>同样也对称正定。因此，<math>\boldsymbol{H}_i</math>可以通过不[[选主元]]的[[LU分解]]分解为

:<math>\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{,}</math>

其中<math>c_i</math>和<math>d_i</math>有简单的递推式：

:<math>\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}</math>

改写<math>\boldsymbol{x}_i=\boldsymbol{x}_0+\boldsymbol{V}_i\boldsymbol{y}_i</math>为

:<math>\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}</math>

其中

:<math>\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}</math>

此时需要观察到

:<math>\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}</math>

实际上，<math>\boldsymbol{p}_i</math>和<math>\zeta_i</math>同样有简短的递推式：

:<math>\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}</math>

通过这个形式，我们得到<math>\boldsymbol{x}_i</math>的一个简单的递推式：

:<math>\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}</math>

以上的递推关系立即导出比共轭梯度法稍微更复杂的直接Lanczos方法。

===从正交性和共轭性导出共轭梯度法===
如果我们允许缩放<math>\boldsymbol{p}_i</math>并通过常数因子补偿缩放的系数，我们可能可以得到以下形式的更简单的递推式：

:<math>\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}</math>

作为简化的前提，我们现在推导<math>\boldsymbol{r}_i</math>的正交性和<math>\boldsymbol{p}_i</math>的共轭性，即对于<math>i\neq j</math>,

:<math>\begin{align}
\boldsymbol{r}_i^\mathrm{T}\boldsymbol{r}_j&=0\text{,}\\
\boldsymbol{p}_i^\mathrm{T}\boldsymbol{Ap}_j&=0\text{.}
\end{align}</math>

各个残量之间满足正交性的原因是<math>\boldsymbol{r}_i</math>实际上可由<math>\boldsymbol{v}_{i+1}</math>乘上一个系数而得。这是因为对于<math>i=0</math>，<math>\boldsymbol{r}_0=\lVert\boldsymbol{r}_0\rVert_2\boldsymbol{v}_1</math>，对于<math>i>0</math>，

:<math>\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}</math>

要得到<math>\boldsymbol{p}_i</math>的共轭性，只需证明<math>\boldsymbol{P}_i^\mathrm{T}\boldsymbol{AP}_i</math>是对角矩阵：

:<math>\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}</math>

是对称的下三角矩阵，因而必然是对角矩阵。

现在我们可以单纯由<math>\boldsymbol{r}_i</math>的正交性和<math>\boldsymbol{p}_i</math>的共轭性推导相对于缩放后的<math>\boldsymbol{p}_i</math>的常数因子<math>\alpha_i</math>和<math>\beta_i</math>。

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

:<math>\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}</math>

类似地，由于<math>\boldsymbol{p}_i</math>的共轭性，必然有<math>\boldsymbol{p}_{i+1}^\mathrm{T}\boldsymbol{Ap}_i=(\boldsymbol{r}_{i+1}+\beta_i\boldsymbol{p}_i)^\mathrm{T}\boldsymbol{Ap}_i=0</math>。于是

:<math>\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}</math>

推导至此完成。

==参考文献==
#{{cite journal|last1 = Hestenes|first1 = M. R.|last2 = Stiefel|first2 = E.|title = Methods of conjugate gradients for solving linear systems|journal = Journal of Research of the National Bureau of Standards|volume = 49|issue = 6|url = http://nvl.nist.gov/pub/nistpubs/jres/049/6/V49.N06.A08.pdf|format = PDF|date = 1952年12月|deadurl = yes|archiveurl = https://web.archive.org/web/20100505153048/http://nvl.nist.gov/pub/nistpubs/jres/049/6/V49.N06.A08.pdf|archivedate = 2010-05-05|access-date = 2010-06-20}}
#{{cite book|last = Saad|first = Y.|title = Iterative methods for sparse linear systems|url = https://archive.org/details/iterativemethods0000saad|edition = 2nd|chapter = Chapter 6: Krylov Subspace Methods, Part I|publisher = SIAM|year = 2003|isbn = 978-0898715347}}

[[Category:数值线性代数]]
[[Category:最优化算法]]