标签归档：线性方程组

Notes on Conjugate Gradient Method – 2

发表评论
版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

打算再将这个主页做一做，先更一波旧笔记。还是那句话，这不是教程，不适合初学。本文中记号大致与前一篇Notes on Conjugate Gradient Method一致。

当\mat{A}对称正定时，共轭梯度法是一个很好的选择。

此时的极小化问题其实是一个\mat{x}与\mat{x}^{*}间的距离极小化问题

(1)   \begin{eqnarray*} &&\min_{\mat{x}}\frac{1}{2}\mat{x}^{T}\mat{A}\mat{x}-\mat{b}^{T}\mat{x}\\ \Leftrightarrow&&\min_{\mat{x}}\lVert \mat{x}-\mat{x}^{*} \rVert_{\mat{A}} \end{eqnarray*}

以\mat{A}为度量，在一个\mat{A}正交标架下，上述距离极小化问题成为一个简单的欧式距离极小化问题，只需在依次各\mat{A}正交方向上极小化距离即可。
这需要解决以下问题：
\begin{enumerate}
\item 一个搜索方向上如何极小化\lVert \mat{x}-\mat{x}^{*} \rVert_{\mat{A}}。
\item 如何高效的构造出与之前搜索过的子空间\mat{A}正交的新搜索方向；
\end{enumerate}

第一个问题。由简单的几何直观，距离的极小化等价于

(2)   \begin{eqnarray*} &&\langle\mat{e}^{(n)},D_{n}\rangle_{\mat{A}}=0\\ \Leftrightarrow&&\langle\mat{r}^{(n)},D_{n}\rangle=0 \end{eqnarray*}

其中D_{n}=\mathrm{span}\{ \mat{d}^{(0)},\dots,\mat{d}^{(n-1)} \}表示\mat{x}=\mat{x}^{(n)}时已经搜索过的子空间。

本式D_{n}换成\mat{d}^{(n)}同样成立。所以，只需搜索到\mat{r}^{(n)}正交于搜索方向\mat{d}^{(n)}时为止即可达到极小化距离的目的。

第二个问题。我们令D_{n}=\mathcal{K}_{m}(\mat{A},\mat{r}^{(0)}),其中\mat{r}^{(0)}为初始残差，于是有D_{n}=\mathrm{span}\{ \mat{r}^{(0)},\mat{A}D_{n-1} \}.这样式2就意味着

(3)   \begin{eqnarray*} \langle\mat{r}^{(n)},D_{n-1}\rangle_{\mat{A}}=0 \end{eqnarray*}

这带来一个极大的便利，即算法可以通过\mat{r}^{(n)}快速构造出与D_n子空间\mat{A}正交的\mat{d}^{(n)},为此只需要将\mat{r}^{(n)}中\mat{d}^{(n-1)}方向的分量（\mat{A}内积意义上）扣除即可。

这样我们可以得到共轭梯度法的算法流程：每次搜索都搜索到残差方向与搜索方向正交为止，新的搜索方向由当前残差用\mat{A}内积与\{ \mat{d}^{(0)},\dots,\mat{d}^{(n-1)} \}做正交化得到，实际上只需要扣除\mat{d}^{(n-1)}方向的分量即可。
算法描述为

(4)   \begin{eqnarray*} \mat{x}^{(n+1)}&=&\mat{x}^{(n)}+\alpha^{(n)}\mat{d}^{(n)}\\ \mat{d}^{(n)}&=&\mat{r}^{(n)}-\beta_{n,n-1}\mat{d}^{(n-1)} \end{eqnarray*}

其中\alpha由\mat{r}^{(n+1)}和\mat{d}^{(n)}的正交性给出

(5)   \begin{eqnarray*} &&\mat{r}^{(n+1)^T}\mat{d}^{(n)}=0\\ \Leftrightarrow&&(\mat{e}^{(n)}+\alpha_{n}\mat{d}^{(n)})^T\mat{A}\mat{d}^{(n)}=0\\ \Leftrightarrow&&\alpha_{n}=-\frac{\langle \mat{d}^{(n)},\mat{e}^{(n)} \rangle_{\mat{A}}}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}}=\frac{\langle \mat{d}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}}\\ &&=\frac{\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}} \end{eqnarray*}

而\beta则这样计算出

(6)   \begin{eqnarray*} &&\mat{r}^{(n)}=-\mat{A}\mat{e}^{(n)}=\mat{r}^{(n-1)}-\alpha_{n-1}\mat{A}\mat{d}^{(n-1)}\\ \Leftrightarrow&&\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle = -\alpha_{n-1} \langle \mat{r}^{(n)},\mat{d}^{(n-1)} \rangle_{\mat{A}}\\ \Rightarrow&&\beta_{n,n-1}=-\frac{\langle \mat{r}^{(n)},\mat{d}^{(n-1)} \rangle_{\mat{A}}}{\langle \mat{d}^{(n-1)},\mat{d}^{(n-1)} \rangle_{\mat{A}}}\\ &&=\frac{1}{\alpha^{(n-1)}}\frac{\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n-1)},\mat{d}^{(n-1)} \rangle_{\mat{A}}}\\ &&=\frac{\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{r}^{(n-1)},\mat{r}^{(n-1)} \rangle} \end{eqnarray*}

共轭梯度法在加入precondition时会遇到一个问题，即如何保证\mat{M}^{-1}\mat{A}仍是对称矩阵。
对于对称正定矩阵\mat{M},存在一个分解\mat{M}=\mat{E}\mat{E}^T.容易证明\mat{M}^{-1}\mat{A}和\mat{E}^{-1}\mat{A}\mat{E}^{-T}具有相同的本征值，因而具有同样的precondition效果。于是可以将问题转化为

(7)   \begin{eqnarray*} \mat{E}^{-1}\mat{A}\mat{E}^{-T}\hat{\mat{x}}=\mat{E}^{-1}\mat{b},\hat{\mat{x}}=\mat{E}^{T}\mat{x} \end{eqnarray*}

再用CG求解。\mat{E}^{-1}\mat{A}\mat{E}^{-T}的几何意义是将问题转入一个新的坐标系下求解，这个坐标系下二次型的超椭球被变换的更接近超球了一些。这样要求计算出\mat{E},但实际上我们可以不用计算\mat{E}而直接用\mat{M}完成计算。为此只需要将原来的迭代公式做坐标变换，写出新坐标系下的迭代公式即可。要注意迭代公式在该变换下并不完全保持形式一致，具体的对应关系如下所示

Rendered by QuickLaTeX.com
有稍许微分几何或相对论背景的同学不难看出这其实是一组在坐标系变换下具有不同变换性质的量之变换式，包含了协/逆变矢量，标量的变换性质，以及如何将当前坐标系下定义，非形式不变的量\langle \mat{r},\mat{r}^{\prime} \rangle提升为任意坐标系下均保持形式不变的量\langle \mat{r},\mat{M}^{-1}\mat{r}^{\prime} \rangle的操作.

替换后可将\mat{E}全部消除掉，剩下仅用\mat{M}表达的迭代公式

(8)   \begin{eqnarray*} \mat{d}^{(0)}&=&\mat{M}^{-1}\mat{r}^{(0)}\\ \mat{d}^{(n)}&=&\mat{M}^{-1}\mat{r}^{(n)}-\beta_{n,n-1}\mat{d}^{(n-1)}\\ \mat{x}^{(n+1)}&=&\mat{x}^{(n)}+\alpha^{(n)}\mat{d}^{(n)}\\ \alpha_{n}&=&\frac{\langle \mat{r}^{(n)},\mat{M}^{-1}\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}}\\ \beta_{n,n-1}&=&\frac{\langle \mat{r}^{(n)},\mat{M}^{-1}\mat{r}^{(n)} \rangle}{\langle \mat{r}^{(n-1)},\mat{M}^{-1}\mat{r}^{(n-1)} \rangle} \end{eqnarray*}