共轭梯度法（Conjugate Gradient）及其相关算法
1. 引言
一般来说，解$Ax=b$这种方程，我们的数值计算书[4]中讲得方法有：直接法（高斯消元，行、列主元，LU分解，LDLT分解），迭代法（Jacobi，Gauss-Seidel，松弛法）。分解法固然是最直接的方法，但是对于大型的稀疏矩阵，用分解法分解出来的矩阵，非零项很多，很吃内存。书中所讲的几种迭代法，迭代次数会比较多，而且迭代次数越多，迭代的效率越低。共轭梯度法是很好地解决了这一问题的算法。
2. 梯度下降法
首先考虑这样一个函数的极值问题，$f(\mathbf{x})=\frac{1}{2}\mathbf{x}^TA\mathbf{x}-\mathbf{x}^Tb+c$。与极值密切相关的是函数的导数，对$f(\mathbf{x})$关于$\mathbf{x}$求导，得到，$\frac{df(\mathbf{x})}{d\mathbf{x}}=\frac{1}{2}(A+A^T)\mathbf{x}-b$，当$A=A^T$时，即$A$为对称阵时，$\frac{df(\mathbf{x})}{d\mathbf{x}}=A\mathbf{x}-b$。这样，求解方程$A\mathbf{x}=b$的问题被转化为一个函数求极值的问题。而且形式上要求$A$为对称阵。也就是说共轭梯度法并不是对所有形式的矩阵A都适用，需要A满足对称正定的条件。
梯度下降法的思想是，给某个解$\mathbf{x_n}$，残差$\mathbf{r_n}=b-A\mathbf{x_n}=-f’(\mathbf{x_n})$。按照牛顿法的思想，让$\mathbf{x_n}$沿着负梯度$\mathbf{r_n}$的方向移动到$\mathbf{x_{n+1}}$。这里要求$\mathbf{x_{n+1}}$是$f(\mathbf{x})$在$\mathbf{x_n}+k\mathbf{r_n}$上的最小值。因此，有$\frac{df(\mathbf{x_n}+k\mathbf{r_n})}{dk}=\mathbf{r_n}^T\mathbf{r_{n+1}}=0$。所以，$\mathbf{r_n}^T(b-A(\mathbf{x_n}+k\mathbf{r_n}))=\mathbf{r_n}^T\mathbf{r_n}-k\mathbf{r_n}^TA\mathbf{r_n}=0$，所以，$k=\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{r_n}^TA\mathbf{r_n}}$。这样，算法的流程就很清楚了。
$\alpha_n=\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{r_n}^TA\mathbf{r_n}}$
$\mathbf{x_{n+1}}=\mathbf{x_n}+\alpha_n\mathbf{r_n}$
$\mathbf{r_{n+1}}=\mathbf{r_n}-\alpha_nA\mathbf{r_n}$
但是，存在一个问题，即便$\mathbf{r_j}^T\mathbf{r_n}=0,j<n$，$\mathbf{r_j}^T\mathbf{r_{n+1}}=\mathbf{r_j}^T\mathbf{r_n}-k\mathbf{r_j}^TA\mathbf{r_n}=0$也不一定满足，$\mathbf{r_{n+1}}$也不能保证与前面$n$个$\mathbf{r_i}$都垂直。即后面前进的方向会有部分与前面的某步重合，这使得梯度下降法的计算效率比较低。注意到，式子中有一项$\mathbf{r_j}^TA\mathbf{r_n}$干扰了正交性，所以需要在这上面进行一些调整。（这可能是共轭梯度法设计出来的早先思路）
3. 共轭梯度法
基于此，定义了一个$\mathbf{d_j}$与$\mathbf{d_k}$关于$A$正交，即$\mathbf{d_j}^TA\mathbf{d_k}=0$。
现在，我们尝试构建下CG的计算方法，从最简单的情形考虑起。设$\mathbf{r_0}=b-A\mathbf{x_0}$，$\mathbf{d_0}=\mathbf{r_0}$。此时，$\alpha_0=\frac{\mathbf{r_0}^T\mathbf{r_0}}{\mathbf{d_0}^TA\mathbf{d_0}}$，$\mathbf{x_1}=\mathbf{x_0}+\alpha_0\mathbf{d_0}$，$\mathbf{r_1}=\mathbf{r_0}-\alpha_0A\mathbf{d_0}$，此时，可以看出，$\mathbf{d_0}^T\mathbf{r_0}=\mathbf{r_0}^T\mathbf{r_0}$，$\mathbf{r_0}^T\mathbf{r_1}=\mathbf{r_0}^T\mathbf{r_0}-\alpha_0\mathbf{r_0}^TA\mathbf{d_0}=0$。可以知道$\mathbf{r_1}$与$\mathbf{r_0}$正交，因此由$\mathbf{r_1}$构建出一个新的更新方向就是(Schmit)，$\mathbf{d_1}=\mathbf{r_1}+\beta_0\mathbf{d_0}$，因为要满足关于$A$正交，所以，$\mathbf{d_0}^TA\mathbf{d_1}=\mathbf{d_0}^TA(\mathbf{r_1}+\beta_0\mathbf{d_0})=0$。此时，利用$A\mathbf{d_0}=\frac{\mathbf{r_0}-\mathbf{r_1}}{\alpha_0}$，得$\frac{\mathbf{r_0}^T-\mathbf{r_1}^T}{\alpha_0}\mathbf{r_1}=-\beta_0\mathbf{d_0}^TA\mathbf{d_0}$，所以，$\beta_0=\frac{\mathbf{r_1}^T\mathbf{r_1}}{\alpha_0\mathbf{d_0}^TA\mathbf{d_0}}=\frac{\mathbf{r_1}^T\mathbf{r_1}}{\mathbf{r_0}^T\mathbf{r_0}}$，$\mathbf{d_1}=\mathbf{r_1}+\frac{\mathbf{r_1}^T\mathbf{r_1}}{\mathbf{r_0}^T\mathbf{r_0}}\mathbf{d_0}$。此时，注意到，$\mathbf{r_0}^TA\mathbf{d_1}=\mathbf{d_0}^TA\mathbf{d_1}=0$，此式备用。
下面来到第$n$步，根据简单情形，我们利用归纳法构建行进的方向。首先，假设，当$i,j\leqslant n \quad (n\geqslant1)$时，有这些关系成立。$\mathbf{d_i}^T\mathbf{r_i}=\mathbf{r_i}^T\mathbf{r_i}$，$\mathbf{r_i}^T\mathbf{r_j}=0\quad(i\not=j)$，$\mathbf{d_i}^TA\mathbf{d_j}=0\quad(i\not=j)$，$\mathbf{r_i}^TA\mathbf{d_j}=0\quad(i<j)$，$\mathbf{d_i}=\mathbf{r_i}+\beta_{i-1}\mathbf{d_{i-1}}$。
现在来到第$n+1$步，首先可以确定，$\alpha_n=\frac{\mathbf{d_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}=\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}$，这样就得到，$\mathbf{r_{n+1}}=\mathbf{r_n}-\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}A\mathbf{d_n}$。此时，证明$\mathbf{r_{n+1}}$与$\mathbf{r_i}(i\leqslant n)$的正交性。$\mathbf{r_n}^T\mathbf{r_{n+1}}=\mathbf{r_n}^T\mathbf{r_n}-\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}\mathbf{r_n}^TA\mathbf{d_n}=\mathbf{r_n}^T\mathbf{r_n}-\frac{\mathbf{r_n}^T\mathbf{r_n}}{(\mathbf{r_n}^T+\beta_{n-1}\mathbf{d_{n-1}}^T)A\mathbf{d_n}}\mathbf{r_n}^TA\mathbf{d_n}=\mathbf{r_n}^T\mathbf{r_n}-\mathbf{r_n}^T\mathbf{r_n}=0$，而$\mathbf{r_i}^T\mathbf{r_{n+1}}=\mathbf{r_i}^T\mathbf{r_n}-\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}\mathbf{r_i}^TA\mathbf{d_n}=0+0,(i<n)$，有了$\mathbf{r_{n+1}}$，就可以构造$\mathbf{d_{n+1}}$，设$\mathbf{d_{n+1}}=\mathbf{r_{n+1}}+\beta_n\mathbf{d_n}$，则$\mathbf{d_n}^TA\mathbf{d_{n+1}}=\mathbf{d_n}^TA(\mathbf{r_{n+1}}+\beta_n\mathbf{d_n})=0$，所以利用$A\mathbf{d_n}=\frac{\mathbf{r_n}-\mathbf{r_{n+1}}}{\alpha_n}$，$\beta_n=\frac{\mathbf{d_n}^TA\mathbf{r_{n+1}}}{\mathbf{d_n}^TA\mathbf{d_n}}=\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\mathbf{r_n}^T\mathbf{r_n}}$，$\mathbf{d_{n+1}}=\mathbf{r_{n+1}}+\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\mathbf{r_n}^T\mathbf{r_n}}\mathbf{d_n}$，现在，证明$\mathbf{d_{n+1}}$关于$A$的正交性。$\mathbf{d_j}^TA\mathbf{d_{n+1}}=\mathbf{d_j}^TA\mathbf{r_{n+1}}+\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\mathbf{r_n}^T\mathbf{r_n}}\mathbf{d_j}^TA\mathbf{d_n}=\frac{\mathbf{r_j}^T-\mathbf{r_{j+1}}^T}{\alpha_j}\mathbf{r_{n+1}}+\beta_n\mathbf{d_j}^TA\mathbf{d_n}=0(j<n)$，对于$j=n$的情形，$\mathbf{d_n}^TA\mathbf{d_{n+1}}=\frac{\mathbf{r_n}^T-\mathbf{r_{n+1}}^T}{\alpha_n}\mathbf{r_{n+1}}+\beta_n\mathbf{d_n}^TA\mathbf{d_n}=-\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\alpha_n}+\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\mathbf{r_n}^T\mathbf{r_n}}\mathbf{d_n}^TA\mathbf{d_n}=0$。然后，还有一些辅助的式子的证明，$\mathbf{d_{n+1}}^T\mathbf{r_{n+1}}=\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}+\beta_n\mathbf{d_n}^T\mathbf{r_{n+1}}=\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}+\beta_n\mathbf{d_n}^T(\mathbf{r_n}-\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}A\mathbf{d_n})=\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}$。$\mathbf{r_i}^TA\mathbf{d_{n+1}}=(\sum_{j}^{i-1}a_j\mathbf{d_j})A\mathbf{d_{n+1}}=0$。至此，我们证明了所有的结论在第$n+1$步后也成立。这就是CG的形式的由来，关键在于两点，一个是关于$A$正交，另一个是 $\mathbf{d_i}=\mathbf{r_i}+\beta_{i-1}\mathbf{d_{i-1}}$的形式的确立。
$\alpha_n=\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}$
$\mathbf{x_{n+1}}=\mathbf{x_n}+\alpha_n\mathbf{d_n}$
$\mathbf{r_{n+1}}=\mathbf{r_n}-\alpha_nA\mathbf{d_n}$
$\beta_n=\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\mathbf{r_n}^T\mathbf{r_n}}$
$\mathbf{d_{n+1}}=\mathbf{r_{n+1}}+\beta _n\mathbf{d_n}$
4. 预条件共轭梯度法(Preconditioner Conjugate Gradient)
按照算法写了代码之后，发现迭代的速度并不是非常快。这与谱半径有关，即最大的特征值绝对值。所以，又有PCG(Preconditioner conjugate gradient)来加速收敛。
PCG的思想是，$A$是对称正定矩阵，所以，可以将$A$写成$A=LL^T$的形式。但是$L$有很多非零项，所以可以设计出某些Preconditioner来近似$A$，这样使得$L^{-1}AL^{-T}\approx I$，而求解$I\mathbf{x}=\mathbf{b}$是非常容易的，所以，这样就使收敛变得更容易了。
这里先推导下引入preconditioner后的CG流程形式是怎么来的。
$A\mathbf{x}=b\Leftrightarrow E^{-1}AE^{-T}(E^T\mathbf{x})=E^{-1}\mathbf{b}\Leftrightarrow \hat{A}\hat{\mathbf{x}}=\hat{b}$
，令$M=EE^T$，CG的流程变为：
$\hat{\mathbf{r}}_n=\hat{\mathbf{b}}-\hat{A}\hat{\mathbf{x}}_n=E^{-1}\mathbf{r}_n$
$\mathbf{\hat{d}_n}=E^{-1}\mathbf{d_n}$
$\hat{\alpha}_n=\frac{\hat{\mathbf{r}}_n^T\hat{\mathbf{r}}_n}{\hat{\mathbf{d}}_n^T\hat{A}\hat{\mathbf{d}}_n}=\frac{\mathbf{r}_n^TE^{-T}E^T\mathbf{r}_n}{\mathbf{d}_n^TE^{-T}E^{-1}AE^{-T}E^{-1}\mathbf{d}_n}$
$\hat{\mathbf{x}}_{n+1}=\hat{\mathbf{x}}_n+\hat{\alpha}_n\hat{\mathbf{d}}_n\Rightarrow \mathbf{x}_{n+1}=\mathbf{x}_n+\hat{\alpha}_nE^{-T}E^{-1}\mathbf{d}_n$
$\hat{\beta}_n=\frac{\hat{\mathbf{r}}_{n+1}^T\hat{\mathbf{r}}_{n+1}}{\hat{\mathbf{r}}_n^T\hat{\mathbf{r}}_n}=\frac{\mathbf{r}_{n+1}^TE^{-T}E^{-1}\mathbf{r}_{n+1}}{\mathbf{r}_n^TE^{-T}E^{-1}\mathbf{r}_n}$
$\hat{\mathbf{d}}_{n+1}=\hat{\mathbf{r}}_{n+1}+\hat{\beta}_n\hat{\mathbf{d}}_n\Rightarrow E^{-T}E^{-1}\mathbf{d}_{n+1}=E^{-T}E^{-1}\mathbf{r}_n+\hat{\beta}_nE^{-T}E^{-1}\mathbf{d}_n$
令$E^{-T}E^{-1}\mathbf{r_n}=M^{-1}\mathbf{r_n}=\mathbf{Z_n}$，$\mathbf{d’_n}=E^{-T}E^{-1}\mathbf{d_n}$，得到，
$\mathbf{Z_n}=M^{-1}\mathbf{r_n}$
$\hat{\alpha}_n=\frac{\mathbf{r}_n^T\mathbf{Z_n}}{\mathbf{d’}_n^TA\mathbf{d’}_n}$
$\mathbf{x}_{n+1}=\mathbf{x}_n+\hat{\alpha}_n\mathbf{Z_{n+1}}$
$\hat{\beta}_n=\frac{\mathbf{r}_{n+1}^T\mathbf{Z_{n+1}}}{\mathbf{r}_n^T\mathbf{Z_n}}$
$\mathbf{d’_{n+1}}=\mathbf{Z_n}+\hat{\beta}_n\mathbf{d’_n}$
常用的Preconditioner，
Jacobi preconditioner：
$A=F+D+F^T$，$M=D$
Incomplete cholesky preconditioner[3]：
$A=F+D+F^T$，$L=FE^{-1}+E$，$LL^T=(FE^{-1}+E)(E^{-1}F^T+E)=FE^{-1}E^{-1}F^T+F^T+F+EE$，为了确定$E$的数值，使$FE^{-1}E^{-1}F^T+F^T+F$的对角项等于$D$。注意$LL^T$并不是完全等于$A$的，有应该是0项的成了非0项。
SSOR preconditioner：
$A=F+D+F^T$，$M=(D+L)D^{-1}(D+L)^T$，
带上参数之后，$M=\frac{1}{2-\omega}(\frac{1}{\omega}D+L)(\frac{1}{\omega}D)^{-1}(\frac{1}{\omega}D+L)^T$。
Preconditioner的关键是需要$\mathbf{Z_n}=M^{-1}\mathbf{r_n}$这一步算起来很容易。上面的几个构建的都是可以直接求解（LDLT，LU，追赶法）的矩阵。
从实践的结果看，加preconditioner的确要快。Incomplete cholesky preconditioner要快于SSOR preconditioner。同一个问题，SSOR要132步收敛，ICP要52步收敛。
5. 总结
接下来进行MultiGrid以及FFT的学习，这两个都可以快速地求解线性方程组而且广泛的应用在很多地方，属于高级技巧，必须掌握。印象中，CUDA里有一个FFT解NS方程的example，可以参考一下。
6. 参考资料
[1] An Introduction to the Conjugate Gradient Method Without the Agonizing Pain
[2] 资料[1]的中文翻译
[3] Fluid simulation for computer graphics
[4] 数值计算方法与算法（第二版）
[5] Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods