数值线性代数: Krylov子空间法

本文旨在总结线性方程组求解的相关算法，特别是Krylov子空间法的原理及流程。

注1：限于研究水平，分析难免不当，欢迎批评指正。

注2：文章内容会不定期更新。

零、预修

0.1 共轭转置

对于、，若矩阵第行第列元素的共轭等于矩阵第行第列元素，即，则称矩阵是矩阵的共轭转置矩阵，记作。

可以看出，若，则。

0.2 Hermite矩阵

对于，若，则称矩阵为Hermite矩阵。

可以看出，若，则，即实数域的Hermite矩阵即是对称矩阵。

0.3 Hessenberg矩阵

0.4 Cholesky分解

对于正定矩阵，则有，其中，矩阵是下三角矩阵，矩阵是矩阵的共轭转置。

若对称正定矩阵，则有，其中，矩阵是下三角矩阵，矩阵是矩阵的转置。

0.5 Arnoldi分解

设，则有，其中，，，，，。

---

下面如无特殊说明，均仅讨论实数域内的线性方程组。

一、直接法求解线性方程组

1.1 LU分解

设，若对于，均有，则存在唯一的单位下三角矩阵和上三角矩阵，使得，并且。

1.2 Cholesky分解

若对称正定，则有，其中，为单位下三角矩阵，为对角元均为正数的对角矩阵。

二、 总论：迭代法求解线性方程组的一般思路

对于非奇异矩阵，，使用迭代法求解线性方程组过程中，一般需要以下流程进行：

1. 给定一个初始向量；
2. 构造一个递推公式；
3. 不断递推，使其近似收敛于；

下表列出了若干迭代算法的迭代公式。

三、Projection Methods

投影法在较小的线性空间内寻找满足精度要求的近似解，也即在某个线性空间内寻找真解的投影。这其实就是投影法得名的原因。

对于非奇异矩阵，考虑线性方程组，其中，，。

令，首先，构造列满秩矩阵与，其中；然后，设真解在线性空间内的投影为，即，令其满足Petrov-Galerkin条件，即，均有。称为搜索空间，称为约束空间。若时，称为正投影算法，否则称为斜投影算法。

若令，则有，其中，，。可以看出，投影法实际上是将阶线性方程组转化为了阶线性方程组()。

讨论：

•  若

以正则化方法为例，即，，由于，另外，，，即是对称正定矩阵。因此，正则化方法实际上将一般矩阵线性方程组转化为对称正定线性方程组。

• 若

由于，，线性方程组的规模由阶降为了阶。因此，可以通过投影法可将问题转换为更为简单的子问题，然后再进行求解。

以求解非对称线性方程组的Arnoldi方法为例，即，其中，，，。则。

综合上述讨论，可归出结投影法的操作流程：

• 使用投影法将问题转化为较易求解的子问题；
• 使用合适的方法求解子问题。

四、Krylov Subspace Methods

Krylov子空间法本质上也是一种投影法，其核心思想是在较小维度的Krylov子空间内寻找满足精度要求的近似解。也就是说，相对于一般投影法，Krylov子空间法选取了Krylov子空间作为搜索空间。

对于线性方程组，Krylov子空间法可以表述为：给定初始解向量，令，选取为维Krylov子空间，即，寻找，使得满足Petrov-Galerkin条件，即。其中，选择不同的，就对应不同的Krylov子空间法，

若设、分别是、的一组基向量，并令，则有。若可逆，则有。

由此可以看出，Krylov子空间法的核心是如何计算与。对于，目前广泛采用的选取方法如下表，

4.1 Steepest Descent Method

4.2 Conjugate Gradient Method

下面考虑对称正定线性方程组，其中，，，。

使用Arnoldi分解定理，可得到矩阵的步Arnoldi分解，即，其中，，，，是Hessenberg矩阵，是阶单位矩阵，。

因为，，可知，由于对称，则矩阵也是对称矩阵，进一步，结合Hessenberg矩阵的定义，可知是三对角矩阵。另外，对于，则，此时，由于正定，则知，即矩阵也是正定矩阵。因此，若矩阵对称正定，矩阵也是对称正定矩阵，更加准确地说，是对称正定三对角矩阵。

令，并令矩阵第行第列元素，则有

，考虑到，因此，可得到Arnoldi分解的递推公式，

由此，可以看出矩阵是Krylov子空间的一组标准正交基。因此，对应于Krylov子空间法，可取。

若给定的初始向量，令近似解，根据Petrov-Galerkin条件，有，记，则，结合Arnoldi分解，有，考虑到与，即有。

根据Cholesky分解定理，因为对称正定，则，其中，矩阵是下三角矩阵，为对角元均为正数的对角矩阵，。

若选取，则是Krylov子空间的一组标准正交基，此时

记，,则有。

由上述分析，矩阵进行步Arnoldi分解，使用Krylov子空间法，可以得到近似解：

，其中，。

类似地，矩阵进行步Arnoldi分解，按照上面地思路，亦可得到近似解：

，其中，。

根据Arnoldi分解过程，很明显，

，

，

，

同时，，并将、、$\boldsymbol{D}_{m+1}$、分块，则有

$\boldsymbol{L}_{m+1}=\begin{pmatrix} \boldsymbol{L}_{m+1}\left ( 1:m,1:m \right ) & 0\\ \boldsymbol{L}_{m+1}\left ( m+1,1:m \right )&\boldsymbol{L}_{m+1}\left ( m+1,m+1 \right ) \end{pmatrix}$

则有，$\boldsymbol{L}_{m+1}\boldsymbol{D}_{m+1}=\begin{pmatrix} \boldsymbol{L}_{m+1}\left ( 1:m,1:m \right )\boldsymbol{D}_{m+1}\left ( 1:m,1:m \right ) & \boldsymbol{0}\\ \boldsymbol{L}_{m+1}\left ( m+1,1:m \right )\boldsymbol{D}_{m+1}\left ( 1:m,1:m \right )&\boldsymbol{L}_{m+1}\left ( m+1,m+1 \right ) \boldsymbol{D}_{m+1}\left ( m+1,m+1 \right )\end{pmatrix}$

$\boldsymbol{H}_{m+1}\left ( 1:m,m+1 \right )= \boldsymbol{L}_{m+1}\left ( 1:m,1:m \right )\boldsymbol{D}_{m+1}\left ( 1:m,1:m \right )\boldsymbol{T}_{m+1}\left ( 1:m,m+1 \right )$

由于、是对称正定三角矩阵，则知$\boldsymbol{H}_{m+1}\left ( 1:m,m+1 \right )=\left ( 0,0,\cdots ,0,\boldsymbol{H}_{m+1}\left ( m,m+1 \right )\right )^{T}$，，$\boldsymbol{H}_{m+1}\left ( m,m+1 \right )=\sum_{k=1}^{k=m}\boldsymbol{L}_{m+1}\left ( m,k \right )\boldsymbol{D}_{m+1}\left ( k,k \right )\boldsymbol{L}_{m+1}\left ( m+1,k \right )$，

$\boldsymbol{H}_{m+1}\left ( m+1,m \right )=\sum_{k=1}^{k=m}\boldsymbol{L}_{m+1}\left ( m+1,k \right )\boldsymbol{D}_{m+1}\left ( k,k \right ) \boldsymbol{L}_{m+1}\left ( m,k \right )$

,

其中，$\tau _{m+1}=\frac{\boldsymbol{H}_{m+1}\left ( m,m+1 \right )}{\boldsymbol{L}_{m+1}\left ( m,m \right )\boldsymbol{D}_{m+1}\left ( m,m \right )}$。

考虑到Cholesky分解的唯一性，则，$\boldsymbol{D}_{m+1}\left ( 1:m,1:m \right )=\boldsymbol{D}_{m}$，，即

$\boldsymbol{L}_{m+1}=\begin{pmatrix} \boldsymbol{L}_{m} &\boldsymbol{0} \\ \boldsymbol{L}_{m+1}\left ( m+1,1:m \right ) & \boldsymbol{L}_{m+1}\left ( m+1,m+1 \right ) \end{pmatrix}$

$\boldsymbol{T}_{m+1}=\begin{pmatrix} \boldsymbol{T}_{m} & \boldsymbol{T}_{m+1}\left (1:m, m+1 \right )\\ \boldsymbol{0} & \boldsymbol{T}_{m+1}\left ( m+1,m+1 \right ) \end{pmatrix}$

$\boldsymbol{D}_{m+1}^{-1}=\begin{pmatrix} \boldsymbol{D}_{m}^{-1} &\boldsymbol{0} \\ \boldsymbol{0} & \frac{1}{\boldsymbol{D}_{m+1}\left ( m+1,m+1 \right )} \end{pmatrix}$

$\boldsymbol{D}_{m+1}^{-1}\boldsymbol{L}_{m+1}^{-1}\boldsymbol{d}_{m+1}=\begin{pmatrix}\boldsymbol{D}_{m}^{-1} \boldsymbol{L}_{m}^{-1}\boldsymbol{d}_{m}\\ \frac{ \boldsymbol{d}_{m+1}\left ( m+1 \right )-\boldsymbol{L}_{m+1}\left ( m+1,1:m \right )\boldsymbol{L}_{m}^{-1}\boldsymbol{d}_{m}}{\boldsymbol{L}_{m+1}\left ( m+1,m+1 \right )\boldsymbol{D}_{m+1}\left ( m+1,m+1 \right )}\end{pmatrix}$

令，

则，$\beta \boldsymbol{V}_{m+1}\boldsymbol{T}_{m+1}^{-1}\boldsymbol{D}_{m+1}^{-1}\boldsymbol{L}_{m+1}^{-1}\boldsymbol{d}_{m+1}=\beta \boldsymbol{V}_{m}\boldsymbol{T}_{m}^{-1}\boldsymbol{D}_{m}^{-1}\boldsymbol{L}_{m}^{-1}\boldsymbol{d}_{m}+\mu _{m+1}\boldsymbol{p}_{m+1}$，

由上述表达式，可得共轭梯度法中近似解的递推公式：

将近似解代入原方程组，可得残差向量的递推公式：$\boldsymbol{r}_{m+1}=\boldsymbol{b}-\boldsymbol{A} \boldsymbol{\tilde{x}}_{m+1}=\boldsymbol{b}-\boldsymbol{A}\left ( \boldsymbol{\tilde{x}}_{m}+\mu _{m+1}\boldsymbol{p}_{m+1} \right )=\boldsymbol{r}_{m}-\mu _{m+1}\boldsymbol{A}\boldsymbol{p}_{m+1}$

由于$\boldsymbol{r}_{m}=\boldsymbol{b}-\boldsymbol{A}\left ( \boldsymbol{x}_{0}+\boldsymbol{V}_{m}\boldsymbol{y}_{m} \right ) =\boldsymbol{r}_{0}-\boldsymbol{A}\boldsymbol{V}_{m}\boldsymbol{y}_{m}$，将步Arnoldi分解代入，则有，其中，。另外，由于，则有，$\boldsymbol{r}_{m}=\boldsymbol{r}_{0}-\boldsymbol{V}_{m}\boldsymbol{V}_{m}^{T}\boldsymbol{r}_{0}-\boldsymbol{v}_{m+1}\left (\boldsymbol{e}_{m}^{T}\boldsymbol{y}_{m} \right )$。如前面分析知，，所以，$\boldsymbol{V}_{m}\boldsymbol{V}_{m}^{T}\boldsymbol{r}_{0}=\left ( \boldsymbol{r}_{0}^{T}\boldsymbol{r}_{0} \right )\boldsymbol{v}_{1}=\boldsymbol{r}_{0}$。因此，，也就是说，$\boldsymbol{r}_{m}$平行于$\boldsymbol{v}_{m+1}$，且。

由于,则有$\boldsymbol{p}_{m+2}=\frac{\boldsymbol{V}_{m+2}\left ( 1:n,m+2 \right )-\tau _{m+2}\boldsymbol{p}_{m+1}}{\boldsymbol{T}_{m+2}\left ( m+2,m+2 \right )}$，也就是说，。实际上，$\boldsymbol{V}_{m}\boldsymbol{T}_{m}^{-1}\boldsymbol{T}_{m+1}\left ( 1:m,m+1 \right )$是的线性组合，是$\boldsymbol{V}_{m+1}=\left ( \boldsymbol{v}_{1},\boldsymbol{v}_{2},\cdots ,\boldsymbol{v}_{m},\boldsymbol{v}_{m+1} \right )$的线性组合。对于，$\boldsymbol{p}_{i}^{T}\boldsymbol{A}\boldsymbol{p}_{m+1}=\boldsymbol{p}_{i}^{T}\left ( \boldsymbol{r}_{m}-\boldsymbol{r}_{m+1} \right )=\boldsymbol{p}_{i}^{T}\left (\boldsymbol{v}_{m+2}\boldsymbol{e}_{m+1}^{T}\boldsymbol{y}_{m+1} -\boldsymbol{v}_{m+1}\boldsymbol{e}_{m}^{T}\boldsymbol{y}_{m} \right )$，根据矩阵正交性，很容易知道，$\boldsymbol{p}_{i}^{T}\boldsymbol{A}\boldsymbol{p}_{m+1}=0$。也就是说，对于，向量与向量关于矩阵共轭。这其实就是共轭梯度法得名的原因。

由上述分析，可知，$\boldsymbol{p}_{m+1}=\frac{\xi _{m+1}}{\mu _{m+1}}\boldsymbol{r}_{m}+\frac{\eta _{m+1}}{\mu _{m+1}}\boldsymbol{p}_{m}$，其中，，$\eta _{m+1}=-\frac{\tau _{m+1}\mu _{m+1}}{\boldsymbol{T}_{m+1}\left ( m+1,m+1 \right )}$，

至此，可以得到简化后的近似解表达式：$\boldsymbol{\tilde{x}}_{m+1}=\boldsymbol{\tilde{x}}_{m}+\xi _{m+1}\boldsymbol{r}_{m}+\eta _{m+1}\boldsymbol{p}_{m}$

相应地，简化后的的残差向量表达式：

由，可知$\boldsymbol{r}_{m}^{T}\boldsymbol{r}_{m}=\mu _{m+1}\boldsymbol{r}_{m}^{T}\boldsymbol{A}\boldsymbol{p}_{m+1}$。即。

由于，可知$\boldsymbol{p}_{m+1}^{T}\boldsymbol{A}\boldsymbol{p}_{m+1}=\frac{\boldsymbol{p}_{m+1}^{T}\boldsymbol{A}\boldsymbol{V}_{m+1}\left ( 1:n,m+1 \right )}{\boldsymbol{T}_{m+1}\left ( m+1,m+1 \right )}=-\frac{\boldsymbol{r}_{m}^{T}\boldsymbol{A}\boldsymbol{p}_{m+1}}{\boldsymbol{e}_{m}^{T}\boldsymbol{y}_{m}\boldsymbol{T}_{m+1}\left ( m+1,m+1 \right )}$。即。

另外，由于，可知$-\frac{\boldsymbol{p}_{m}^{T}\boldsymbol{A}\boldsymbol{r}_{m}}{\boldsymbol{e}_{m}^{T}\boldsymbol{y}_{m}}=\tau _{m+1}\boldsymbol{p}_{m}^{T}\boldsymbol{A}\boldsymbol{p}_{m}$。即。

相应地，$\boldsymbol{p}_{m+1}=\frac{\xi _{m+1}}{\mu _{m+1}}\boldsymbol{r}_{m}+\frac{\eta _{m+1}}{\mu _{m+1}}\boldsymbol{p}_{m}$。

以上即为整个共轭梯度法的推导过程。

4.3 Preconditioned Conjugate Gradient

参考书籍

Golub G H , Loan C F V .Matrix Computations.Johns Hopkins University Press,1996.

Ford W .Numerical Linear Algebra with Applications using MATLAB. 2014.

徐树方. 数值线性代数(第二版).  北京大学出版社, 2010.

参考文献

Hestenes M R , Stiefel E L .Methods of Conjugate Gradients for Solving Linear Systems. Journal of Research of the National Bureau of Standards (United States), 1952. 

Arnoldi W E .The principle of minimized iterations in the solution of the matrix eigenvalue problem.Quarterly of Applied Mathematics, 1951, 9(1).

Saad Y .Krylov Subspace Methods for Solving Large Unsymmetric Linear Systems.Mathematics of Computation, 1981, 37(155):105-105.