在看DynamicFusion这篇文章时，同学们有可能会主要到文章运到了一个很让人百思不解的solver叫做Pre-conditioned Conjugate Gradient Solver。在****上转了一圈发现，对于conjugate gradient method深入讲解的文章少之又少，所以决定写一篇能让有一定线性代数基础的同学理解得懂的文章。Ok, let’s get started.

迭代优化问题

假设我们有一个线性问题（用$A$A，$b$b来去求$x$x）：
$Ax = b$Ax=b
如果用直接法（direct method）来求得话，我们可以直接得到：
$x = (A^TA)^{-1}A^Tb$x=(ATA)−1ATb
这种方法好但是如果$A$A是一个$R^{n\times n}$Rn×n，$n > 10^8$n>108的矩阵的话，对$A^TA$ATA直接做Inverse显然是不可取的。于是非直接的迭代优化方式就此诞生。简单来说间接法是用迭代的方式更新$x$x来缩小$Ax$Ax和$b$b的距离。
$\delta x = \argmin_{\delta x}(E(x + \delta x)) = \argmin_{\delta x} (\frac{1}{2}||A(x + \delta x)-b||^2)$δx=δxargmin​(E(x+δx))=δxargmin​(21​∣∣A(x+δx)−b∣∣2)
$x' = x + \delta x$x′=x+δx
不断迭代知道更新后的$x'$x′对距离影响与之前相比不大为止（convergence）。

学视觉或者做神经网的朋友可能都会经常接触到一些常用的迭代优化的算法，如：

1. 梯度下降：$\delta x = \alpha \frac{\partial E}{\partial x}$δx=α∂x∂E​，$\alpha$α是更新大小（step size）
2. 高斯牛顿：$\delta x = (J^TJ)^{-1}J^Te$δx=(JTJ)−1JTe，$J$J是$E$E对于$x$x的雅可比
3. Levenberg-Marquardt：$\delta x = (J^TJ + \lambda I)^{-1}J^Te$δx=(JTJ+λI)−1JTe，$\lambda$λ是一个梯度下降和高斯牛顿的过渡值。
4. 共轭梯度下降：$\delta x =x_k + \alpha p_k$δx=xk​+αpk​，$p_k$pk​是Krylov子空间的一个基向量

Krylov子空间

要想理解共轭梯度下降的动机和原理，我们上来要重温一下Krylov子空间和Krylov序列（Krylov sequence）这些经典的线性代数概念。概念很深奥，我尽量去仔细剖析，记得第一次听学校教授讲解Krylov子空间的时候有种要撞墙的冲动（翻出自己笔记。。。）

我们先来看一个1维的公式，假设我们有一个$x$x的迭代过程：
$x_{n+1} = ax_{n}+b$xn+1​=axn​+b
那么我们如果展开来看的话就会发现：

$x_{n+1} = a^nx_0 + \frac{n!}{n-1!}a^{n-1}bx_0+ ... + \frac{n!}{n-1!}ab^{n-1}x_0 + b^n$xn+1​=anx0​+n−1!n!​an−1bx0​+...+n−1!n!​abn−1x0​+bn

如果我们将$b$b看为常数，而$a$a为一个变量的话，我们会发现$x_{n+1}$xn+1​所在的维度空间由$n$n个基组成：

$x_{n+1}\in Span\{1, a, a^2,...,a^{n-1}, a^n\}$xn+1​∈Span{1,a,a2,...,an−1,an}

而这些基所组成的数列就叫做Krylov数列（Krylov Sequence），而如果只用从$a^0$a0到$a^i$ai，$i < n$i<n基所组成的空间则叫做Krylov子空间（Krylov subspace)

现在我们就可以把$a$a这个常量延伸到一个多维矩阵中，得到（这里是不是很像迭代过程）：
$x_{n+1} = Ax_n + b$xn+1​=Axn​+b
那么$x_{n+1}$xn+1​所在的空间则为：
$x_{n+1} \in Span \{1, A, A^2, ..., A^{n-1}, A^n\}$xn+1​∈Span{1,A,A2,...,An−1,An}
这里我们可以看出，从$x_0$x0​到$x_{n+1}$xn+1​的过程是一个维度增长的过程，而$x_i$xi​所在的空间又成了$x_{i+1}$xi+1​的子空间。既然这样，那么我们能否去重复利用前一个子空间的基去更好地诠释下个延伸后的下一个空间呢？这就是Krylov数列（Krylov Sequence）的重要所在。

如图所示，$S(x_0)$S(x0​)的维度空间是$S(x_1)$S(x1​)空间的子空间，而$S(x_1)$S(x1​)是$S(x_2)$S(x2​)的子空间，由此类推。由于$S(x_0)$S(x0​)是$S(x_1)$S(x1​)的子空间，所以$S(x_0)$S(x0​)的所有基也是$S(x_1)$S(x1​)的基， 由此类推。

所以我们可以利用迭代过程这个一点点累加的更新性质，以Krylov的形式去更新我们高维度梯度下降的方向从而加速其优化速度。

我们发现从$x_{n} \rightarrow x_{n+1}$xn​→xn+1​的过程中，$S(x_{n+1})$S(xn+1​)的基只比$S(x_n)$S(xn​)的基多了一个，而共轭梯度会根据residual的方向在每步生成一个基然后再这个基上做梯度下降

共轭梯度下降

Conjugate Gradient算法其实就是在更新过程中根据residual的方向来生成一个Krylov子空间的基（$p_k$pk​），然后再这个基（$p_k$pk​）的方向上做梯度下降。
$\delta x = \alpha_k p_k$δx=αk​pk​

这个新的基$p_k$pk​其实是先将residual空间进行对之前基的降维，然后在没有被其他基覆盖过得空间中去一个基梯度下降的方向。

$p_k = r_k - \sum_{i=0}^{k-1}\frac{<p_i, r_k>}{<p_i,p_i>_A}p_i$pk​=rk​−i=0∑k−1​<pi​,pi​>A​<pi​,rk​>​pi​

这个部分很像Gram-Schmidt垂直法降维的过程（可以参考QR Factorization: Gram-Schmidt Orthogonalization）

如图所示，Gram-Schmidt是如何利用projection来从$v_2$v2​提取$e_2$e2​这个基的。

（更新中。。。）