共轭梯度

数值优化 线性代数

---

本文介绍和总结共轭梯度(conjugate gradient)算法的如下内容：

共轭梯度算法的背景和由来
梯度算法的几何理解和共轭梯度算法的流程
基于SSOR预条件子的共轭梯度算法
非线性共轭梯度算法

背景

共轭梯度算法的由来，源于求解如下所示的二次型问题：
$f(x ) =\frac{1}{2}x^TAx-b^Tx+c\tag{1}$
如果对于问题的求解，可以被转化为如上所示的二次型表达式(Quadratic Form)的极值，那么求解这个问题的方法也称为二次规划(Quadratic Programming)
二次规划问题，有一个极好的性质，这个性质直接导致了该问题的求解有较好的几何解释，这个性质就是如果矩阵$A$为对称矩阵，那么把这个函数对列向量$x$求导取极值，得到的梯度方向向量$x$满足：
$f^{'}(x)=Ax-b=0 \implies Ax=b \tag{2}$
这个$Ax=b$是经典的线性代数求解线性方程组的问题，求极值即是求解线性方程组，而矩阵$A$的性质，直接决定了梯度方向的情况。
对一个从$n$维到$1$维的函数映射$f(x)$，对列向量求导即：

$$f^{'}(x)=\begin{bmatrix} \frac{\partial f(x)}{\partial x_{1}} \\ \vdots\\ \frac{\partial f(x)}{\partial x_{n}} \end{bmatrix}\tag{3}$$
如果仅仅只看二次型中的二次项，那么
$$g(x)=x^TAx=\sum_{i=0}^{n}\sum_{j=0}^{n}a_{ij}x_ix_j \\ =... + \sum_{i=0}^{n}a_{ik}x_ix_k+\sum_{j=0}^{n}a_{kj}x_kx_j-a_{kk}x_kx_k + ...\\ =...+\sum_{i=0,i\ne k}^{n}a_{ik}x_ix_k + \sum_{j=0,j\ne k}^{n}a_{kj}x_kx_j+a_{kk}x_kx_k+...\\ \implies \frac{\partial g(x)}{\partial x_k}=\sum_{i=0}^{n}a_{ik}x_i+\sum_{j=1}^{n}a_{kj}x_j \tag{4}$$
仔细观察上述表达式，可以发现，求导结果的第$k$个元素等于矩阵$A$的第$k$行和第$k$列分别与向量点积的结果，所以
$$\frac{\partial g(x)}{\partial x}=A^Tx+Ax \implies \frac{\partial f(x)}{\partial x}=\frac{1}{2}(A^T+A)x-b\tag{5}$$
如果A是对称矩阵的形式，那么，通过求导求解极值相当于求解方程组
$$Ax=b \tag{6}$$
二次规划是一些问题的最简单的建模表达形式，比如，在EDA(电子电气自动化)软件开发中，有这样的需求，就是要对各个电子元器件在电路板或者FPGA的电路平面进行自动化布局。每个元器件可以抽象为一些有面积的方型块，然后通过设定一些约束，定义损失函数，最后，通过不断迭代减少损失函数，更新各个电路元件的坐标信息，最终得到各个电子元器件在电路中的较优的全局布局结果。这些图可以有多层，而且损失函数通常也是非线性的，迭代求解通常基于各种梯度下降算法。期间需要考虑各个图层每一块面积的覆盖率，重叠程度等等参数。
如果对问题进行全面的建模，那么自然会比较复杂，如果只考虑最简单的情况，通过把元器件视为一个点，同时通过引入锚点，那么这个问题可以被简化为二次规划问题，通常二次规划问题用于电子元器件的预布局，然后才引入非线性函数进一步处理。
如果优化的目标可以被简化为让各个元器件之间的$x$和$y$坐标之间的距离尽可能分散，那么此时损失函数等于
$$cf=\sum_{i,j}w_{ij}((x_i-x_j)^2+(y_i-y_j)^2)+ \\ \sum_{f,H(f)}w_{fH(f)}((x_f - x_{H(f)})^2+(y_f-y_{H(f)})^2) \tag{7}$$
如果令$x=(x_1,x_2,...,x_n)^T$，$y=(y_1,y_2,...,y_n)^T$，那么，完全可以表达成两个二次型表达式，其中后半部分的$w_{fH(f)}$表示的是锚点到可动点的权重，实际优化该函数，锚点部分只会对$A$的对角线，以及$b^T$有影响，优化的结果就是在锚点的约束下，各个点的距离之间尽可能小。
由于问题的特殊形式，问题的求解完全也可以把$x$和$y$左边分别分成两个部分：
$$Ax=b_x \tag{8}$$ $$Ay=b_y \tag{9}$$
也就是等价于求解线性方程组了

几何理解和共轭梯度算法的流程

共轭梯度算法来源于二次型问题的求解，二次型问题可以在几何上直观的理解，即使是最速梯度下降(Steepest Descent)，也可以这样来理解。

二次型的图形

二次型的几何图解完全取决于矩阵$A$，如果$A$是正定的矩阵，那么图像将会是：

也就是是一个碗状朝上的结构，因为对于正定矩阵来说，任何非零向量，都有：

$$x^TAx>0 \tag{10}$$
这个定义的本质是，所有的特征向量所对应的特征值都是正的，因为：
$$x^TAx=x^T\lambda x>0 \implies \lambda > 0 \tag{11}$$
对于二次型来说，求导的梯度方向由$Ax=b$决定，$A$正定(positive definite)，意味着各个向量$x$变换到$A$的列空间中得到的向量$Ax$都是单调递增的，所以二次型表达式(如果是二维的)的几何图形就是中间凹陷，四周往上的碗状。同理，可以定义半正定，负定等，那么图形就是：

对于既非负定也非正定的情况，就是马鞍型，这个不管使用什么梯度方法，理论上都可能不会收敛。对于更高维度的情况，其实也是类似的，这里边正定或者负定的属性，决定了梯度方向是不是单调递增/递减的。对二次型，后边都是基于对称正定矩阵来分析的，那些梯度方法也源于此来进行基本的几何分析。

Steepest Descent

对于最速梯度下降(Steepest Descent)而言，选定了初始的方向向量$x_0$后，后续每一次迭代所选择的方向向量由当前的梯度方向决定：

$$x_{i+1}=x_i+\alpha r_i \tag{12}$$
其中$r_i$是方向向量，默认就是取梯度的负方向。这个过程的本质，其实就是线搜索(line search)，如下图所示：

上图是二次型在变量维度为二维的时候的共势等高图，线是梯度负方向，每次迭代的本质就是在沿着这条线走，直到到达某点后停止，得到此时的向量$x_{i+1}$，一般情况下，步长$\alpha$是可以固定的，特别是在函数是非线性，且非常复杂的时候，要想每一次都计算最好的步长，可能是非常复杂且不划算的，不过对于二次型问题而言，我们却可以设法找到最优的$\alpha$，通过查看上图的各个结果向量，为了更接近与最小点，$r_{i+1}$应该和$r_i$正交，此时，才是沿着当前的搜索方向的最优解。因为这个时候，向量已经沿着该方向尽可能往这个方向走，且已经不可能贡献再多的效用了，再走就又回去了。
对于二次型，梯度负方向等于：
$$f^{'}_i(x)=Ax_i-b \implies r_i=-f^{'}_i(x)=b-Ax_i \tag{13}$$
利用正交关系进行推导，可以得到：
$$r_{i+1}^Tr_i=(b-Ax_{i+1})^Tr_i=(b-A(x_i+\alpha r_i))^Tr_i=0 \\ (b-Ax_i)^Tr_i-\alpha r_i^TA^Tr_i=r^T_i-\alpha r_i^TAr_i=0$$
所以
$$\alpha=\frac{r_i^Tr_i}{r_i^TAr_i} \tag{14}$$
这里边每次给出$r_i$的时候都相当于直接求导了，如果根据迭代的形式给出下一个方向向量，那么
$$r_{i+1}=b-Ax_{i+1}=b-A(x_i+\alpha r_i)=r_i-\alpha Ar_i \tag{15}$$
用迭代的形式给出每一次的梯度方向向量，虽然也许能够节省计算开销，但是限于实际机器的浮点精度，可能会产生累积误差，所以一般也需要在迭代了几十次后，重新严格计算梯度，以避免传递误差过大。
这样我们就可以得到最速梯度下降算法的流程是：

import numpy as np
def steepest_descent(A, b, x_initial, max_step, thresh=0.00001):
    assert(isinstance(A, np.matrix))
    assert (isinstance(b, np.matrix))
    assert (isinstance(x_initial, np.matrix))
    x = x_initial
    for _ in range(max_step):
        r = b - A * x
        alpha = (r.transpose() * r) / (r.transpose() * A * r)
        x = x + r * alpha
        dist = np.sqrt(np.sum(np.square(b - A * x)))
        if dist < thresh:
            break
    return x
if __name__ == '__main__':
    N = 100
    Ar = np.mat(np.random.rand(N, N))
    As = Ar * Ar.transpose()  # get positive definite matrix
    bn = np.mat(np.random.rand(N, 1))
    xi = np.mat(np.zeros((N, 1)))
    xr = steepest_descent(As, bn, xi, 1000)
    print('1000:', np.sqrt(np.sum(np.square(bn - As * xr))))
    xr = steepest_descent(As, bn, xi, 10000)
    print('10000:', np.sqrt(np.sum(np.square(bn - As * xr))))
    xr = steepest_descent(As, bn, xi, 20000)
    print('20000:', np.sqrt(np.sum(np.square(bn - As * xr))))

上边的Python代码体现了SD的最直接流程，梯度都是直接计算的，实际真的要应用，还要考虑稀疏矩阵，计算效率问题，否则效率太低了，也没有使用价值。当然，这个只用于二次型问题。
对一些非线性函数问题，就需要通过别的方式求导，而且步长$\alpha$也需要确定，或者设定个固定值！如下所示是该代码的一次运行结果：

1000: 1.9453821317304896
10000: 0.9116628408830175
20000: 0.6950563935826073

Conjugate Gradient

运行SD的例子，会发现，随机生成的一个对称矩阵，迭代了几万次仍然不收敛，可以见得，按照SD的理论公式求解$Ax=b$效果有些不理想，收敛过慢。在数值优化里边，一个矩阵迭代求解方法是否收敛，取决于这个矩阵是否严格对角占优，即

$$A=(a_ij), \forall i \in [1, n], |a_{ii}| \gt \sum_{i \ne j}|a_{ij}| \tag{16}$$
越是严格对角占优，越是收敛得越快。
SD这个方法，使用梯度方向，每次选择梯度方向都和上一次的正交，这样迭代的结果取决于初始点的选择和$A$，如果初始点不理想，可能一直在锯齿状的接近极值点，但是迭代多次后仍然不能命中。
对于当前的二次型而言，不是多么复杂的问题，变量的维度也许就是$n$而已，存不存在这样的迭代方法，使得每一次迭代，都消减一个维度方向，后边的方向和前边的方向都完全正交呢？这样的话，我只需要$n$步迭代就可以直接命中极值点了。这个思想，我觉得就是CG背后的想法。如果选择$n$个线性无关的向量$u_1,u_2,...,u_n$，那么依据施密特(Gram-Schmidt Process)正交化的过程，每个方向向量为：
$$d_i=u_i+\sum_{k=0}^{i-1}\beta_{ik}d_k \tag{17}$$

Conjugate Direction

在CG中，引入了共轭方向的概念。这个概念本质上就是对正交的推广，任何向量在欧几里得空间内正交意味着：

$$x^Ty=0 \tag{18}$$
如果引入矩阵$A$，且
$$x^TAy=0 \tag{19}$$
那么称向量$x$和$y$是$A$正交的，这里边$A$都是方阵，看起来就像是把变量线变换到$A$的行列空间中，然后这些向量在$A$的行列空间中正交。这些两两$A$正交的向量被称为共轭方向向量。如果给个图解，那么就会如下图所示：

前后两个方向向量是$A$正交(A-orthogonal)的，意味着$d_i^TAd_j=0$,即使这两个向量在当前空间不正交，在$A$的空间中正交也行，可以想象，当前空间中的所有向量都变换到$A$的空间后，着两个向量在$A$空间中就是正交的。

二次型的标准化

上边的过程容易让人联想起线性代数中的二次型标准化过程，二次型标准化，首先通过平移变换，让二次型变成了如下的形式：

$$f(x)=x^TAx \tag{20}$$
然后找寻$A$的$n$个标准正交的特征向量$C=(c_1,c_2,...,c_n)$，然后做如下变换：
$$Cy=x \implies f(x)=x^TAx=y^TC^TACy=y^T\Lambda y \tag{21}$$
最终$A$变成了对角矩阵的形式，这个时候对应的二次型的几何图形就是很标准的碗状了
个人觉得，$A$正交也是希望在$A$标准化后的空间中正交。其实标准化前后，所代表的线性空间是不变的。在几何理解里边，都以标准化后的空间中的向量来分析理解。

CG的推导过程

二次型的基本迭代流程还是不变的：

$$x_{k+1}=x_{k}+\alpha_{k}d_k \tag{22}$$
如果采用迭代求解梯度，那么下边的等式仍然成立：
$$r_{k+1}=r_k-\alpha_{k}Ar_k \tag{23}$$
对等式22向$\alpha_{k}$求导：
$$\frac{\partial f(x_{k+1})}{\partial \alpha_{k}}=f^{'}(x_{k+1})^T\frac{\partial x_{k+1}}{\partial \alpha_k} \\ =-r_{k+1}^Td_k=-(Ax-Ax_{k+1})^Td_k=-d_k^TAe_{k+1}=0 \\ \implies d_k^TAe_{k+1}=0 \tag{24}$$
可以看到，最优的步长是使得下一个误差方向和当前的方向$A$正交，由于
$$e_{k+1}=e_k+\alpha_kd_k \tag{25}$$
所以
$$d_k^TA(e_k+\alpha_kd_k)=0 \implies \alpha_k=-\frac{d_k^TAe_k}{d_k^TAd_k}=\frac{d_k^Tr_k}{d_k^TAd_k}\\ e_k=x_k-x, Ax=b \tag{26}$$
最优的步长因子和当前的梯度和方向向量均相关。从上边的表达式知道，当前反向和剩余误差方向是$A$正交的，而我们也希望：
$$d_k^TAd_{k+1}=0 \tag{27}$$
如果下一个方向向量和当前的也$A$正交就好了，当然了，如果下一个方向向量$d_{k+1}$刚好选择到了剩余的误差向量$e_{k+1}$，那么下一次迭代就会直接命中极值点了，但是一般这个是不可能的，因为$e_{k+1}$是依赖于$\alpha_k$的，会有循环求解的问题，除非我们已经知道解了，否则我们是得不到$e_{k+1}$的，而且在高维空间中，和$d_k$向量$A$正交的向量是有很多个的，我们希望的是找到一个$d_{k+1}$满足上边表达式即可，于此同时，也满足
$$d_k^TAd_m=0, \forall m \lt k \tag{28}$$
这样每次迭代都减少了一个维度，那么，最多$n$步就可以收敛了。
这个过程，在变量为二维的情况下理解，大概如下图所示：

这些椭圆是三维图形在二维的等势面，其实就是$f(x)$相等的值所构成的一个椭圆。
一开始选择梯度作为初始方向，也就是等势椭圆的切线，下一步，由于取$A$正交关系，$d_1$和$e_1$是同向的，两步就命中结果了。对应到$A$的空间中，这个椭圆就是个圆。变量是四维的情况，也可以用下图表示：

变量为三维的情况，加上因变量实际上是四个维度的，由于二次型也可以描述能量关系，所以相等能量的各个三维自变量构成的图形在三维就是个椭球体。映射到$A$的空间中其实就是个标准的球型。每次迭代搜寻的方向向量，都是和过去所有的方向向量$A$正交的，三维的情况可以看得比较明显，那就是$d_1^Td_0=0, e_1^Td_0=0, e_1 \ne d_1$，也就是方向向量在高维度的情况下，即使和过去的方向向量都$A$正交，但是也不能说这个方向向量就是当前的错误向量。但是从事实出发，错误向量$e$随着每次迭代，都减少了各个方向向量的那部分，所以，错误向量$e_0$即，一开始，初始点到极值点的错误向量可以表示为：
$$e_0=\sum_{j=0}^{n-1}\delta_jd_j \tag{29}$$
而等式28中又限定了各个方向向量的关系，所以可以：
$$d_k^TAe_0=\sum \delta_kd_k^TAd_j=\delta_kd_k^TAd_k \\ \implies \delta_k=\frac{d_k^TAe_0}{d_k^TAd_k}=\frac{d_k^TA(e_0 + \sum_{j=0}^{k-1}\alpha_jd_j)}{d_k^TAd_k}=\frac{d_k^TAe_k}{d_k^TAd_k} \\ \implies \delta_k=\frac{d_k^Tr_k}{d_k^TAd_k}=-\alpha_k \tag{30}$$
所以
$$e_0=-\sum_{j=0}^{n-1}\alpha_jd_j \\ \implies e_k=e_0 + \sum_{j=0}^{k-1}\alpha_jd_j=-\sum_{j=k}^{n-1}\alpha_jd_j \tag{31}$$
上边的推到也证明，只要我们满足条件

条件1：各个方向向量$d_i$和$d_j$之间两两$A$正交
条件2：通过求导让每一次迭代的步长因子$\alpha$取得最优

那么每一次迭代求解的过程，其实就是相当于减少了这个方向向量上的误差，每一步的误差向量可以被表示为方向向量的连加的形式，这样只需要迭代$n$步后就会直接收敛到极值点。这个和我们在几何上边的理解是完全一致的。
只是，我们该如何选择这些方向向量$d_k$呢？
如果仅仅只是随便选择一组线性无关的构造基$(u_1,u_2,...u_n)$，然后如同等式17那样构造方向向量，并满足$A$正交，那么：

$$\forall i > j, d_i^TAd_j=u_i^TAd_j+\sum_{k=0}^{i-1}\beta_{ik}d_k^TAd_j=0 \\ \implies \beta_{ij}=-\frac{u_i^TAd_j}{d_j^TAd_j} \tag{32}$$
通过这一种方式，也可以得到如何去构造所有的方向向量，只是，任何方向向量$d_i$都和过去的$i-1$个向量相关，而且每次迭代还需要求解系数，这样就直接导致了计算所需的内存和时间开销巨大，没有实用价值。所以这一组构造基不应该随便选。而应该让等式32得出的各个系数拥有良好的关系，减少计算所需的开销。
事实上，共轭方法，在一开始提出的时候，确实存在实现上的困难，不过通过仔细查看这些关系，我们还没有把梯度向量引入。
在计算过程中，我们始终要得到当前点的梯度，在上边两个条件满足的情况下，观察错误向量，梯度等之间的关系：
$$\forall i \lt j, -d_i^TAe_j=-\sum_{j=p}^{n-1}\delta_jd_i^TAd_j=-d_i^Tr_j \\ \implies d_i^Tr_j=0 \tag{33}$$
也就是说，任何迭代后期的梯度$r_j$都是和之前的方向$d_i$在当前空间正交的，进一步的
$$d_i^Tr_j=u_i^Tr_j+\sum_{k=0}^{i-1}\beta_{ik}d_k^Tr_j=0 \\ \implies u_i^Tr_j=d_i^Tr_j \ne 0, j = i \\ \implies u_i^Tr_j=d_i^Tr_j=0, \forall j \gt i \tag{34}$$
结合前边的所有关系，可以得到如下图所示的关系图：

图中所示，$u_2$和$d_2$到$r_2$的投影是相等的，所以末端所在的平面和$d_0$以及$d_1$所在的平面共面。通过上边的各个关系可以知道，没有必要取选择别的基${u_k}$，直接让$u_k=r_k$就是最优的选择，这个时候，各个向量的关系图将会变成如下所示：

这样做的直接结果有：
$$\alpha_k=\frac{d_k^Tr_k}{d_k^TAd_k}=\frac{r_k^Tr_k}{d_k^TAd_k} \tag{35}$$
$$\beta_{ij}=-\frac{u_i^TAd_j}{d_j^TAd_j}=-\frac{r_i^TAd_j}{d_j^TAd_j},\forall i\gt j \\ r_i^Tr_{j+1}=r_i^T(r_j-\alpha_jAd_j) \\ \implies r_i^TAd_j=\frac{1}{\alpha_j}(r_i^Tr_j-r_i^Tr_{j+1}) \\ \implies r_i^TAd_j=0,\beta_{ij}=0,\forall i>j+1 \\ \implies r_i^TAd_j=-\frac{1}{\alpha_j}r_i^Tr_i,\\ \beta_{i,j+1}=\frac{r_{i}^Tr_{i}}{r_{i-1}^Tr_{i-1}},i=j+1 \tag{36}$$
综合前边所属，CG的流程可以被表示为：

$d_0=r_0=b-Ax_0$
for loop until reach max_step or cost_function below threshold:
$\qquad \alpha_k=\frac{r_k^Tr_k}{d_k^TAd_k}$
$\qquad x_{k+1}=x_k+\alpha_kd_k$
$\qquad r_{k+1}=r_k-\alpha_kAd_k\space or \space r_{k+1}=b-Ax_{k+1}$
$\qquad \beta_{k+1}=\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}$
$\qquad d_{k+1}=r_{k+1}+\beta_{k+1}d_k$
return $x_n$

使用Python代码实现：

import numpy as np
def conjugate_gradient(A, b, x_initial, max_step, threshold=0.00001):
    assert(isinstance(A, np.matrix))
    assert(isinstance(b, np.matrix))
    assert(isinstance(x_initial, np.matrix))
    r_old = b - A * x_initial
    d = r_old
    x = x_initial
    for i in range(max_step):
        alpha = (r_old.transpose() * r_old) / (d.transpose() * A * d)
        x = x + d * alpha
        # r_new = b - A * x
        r_new = r_old - A * d * alpha
        beta = (r_new.transpose() * r_new) / (r_old.transpose() * r_old)
        d = r_new + d * beta
        r_old = r_new
        cf = np.sqrt(np.sum(np.square(b - A * x)))
        if cf < threshold:
            print("Using step: ", i)
            break
    return x
if __name__ == '__main__':
    N = 200
    Ar = np.mat(np.random.rand(N, N))
    As = Ar * Ar.transpose()
    bn = np.mat(np.random.rand(N, 1))
    xi = np.mat(np.random.rand(N, 1))
    xr = conjugate_gradient(As, bn, xi, 1000)
    print('1000:', np.sqrt(np.sum(np.square(bn - As * xr))))
    xr = conjugate_gradient(As, bn, xi, 10000)
    print('10000:', np.sqrt(np.sum(np.square(bn - As * xr))))

运行上述程序的一个输出：

Using step:  410
1000: 6.43707958798e-06
Using step:  410
10000: 6.43707958798e-06

大约需要两倍于维度的迭代次数才可以收敛到要求的范围，虽然理论上100维的变量值需要100次迭代，但是实际上各种浮点数运算误差导致不可能实现这样理想化的结果。通过对比会发现，迭代效果要比SD要好。当然了，实际使用的时候，判断损失函数其实没必要通过那样的形式。而且可以几十次迭代后才计算一下等等。

基于SSOR预条件子的共轭梯度算法

方程组的迭代求解

上边的迭代算法，本质上其实还是在求解形如$Ax=b$的方程。毕竟最终的结果，都是让$e=Ax-b$尽可能的接近于零。一般情况下，对于$Ax=b$的求解，通常直接通过求逆$x=A^{-1}b$得到，这就存在求矩阵的逆可能不存在等等问题。如果$A$不是方阵，且$m \gt \gt n$，那么就是在求解超定方程组，方程组可能没有解，求解$Ax=b$的过程就会等效于最小二乘，因为，为了让误差最小，其实就是让最终的的解$Ax=b^{'}$得到的$b^{'}$离$b$最近，那么$b^{'}$只能是$b$到$A$的列空间的投影，所以它们之间的差和$A^T$的列空间正交，所以有

$$A^T(Ax-b)=0 \implies x = (A^TA)^{-1}A^Tb \tag{37}$$
即使$A$是非方阵，但是$A^TA$却是对称方阵，所以求解逆是可行的，这个时候其实还是转换成了
$$A^TAx=A^Tb,b^{'}=A^Tb, A^{'}=A^TA \implies A^{'}x=b^{'} \tag{38}$$
所以还是可以使用上边的梯度迭代方法来求解。所以可以见得，线性系统问题，即可以表达为求解$Ax=b$的系统问题，还是比较多的。
求解的方法，除了直接求逆，以及高斯的LU分解，其它的基本就是迭代方法了。最早的迭代有高斯-赛德尔迭代，SOR等，它们都是首先把$A$拆分:
$$A=D+L+U \tag{39}$$
也就是把矩阵拆解成对角，上三角，下三角的形式，对于高斯-赛德尔方法，把$(D+L+U)x=b$变成了：
$$(L+D)x_{k+1}=-Ux_k+b \\ x_{k+1}=D^{-1}(b-Ux_k-Lx_{k+1})\tag{40}$$
提取成对角，上下三角矩阵的原因，恐怕还是因为逆矩阵好算。上边的过程，每次迭代其实都是循环回代的过程，所以等式40右边有$x_{k+1}$，其实是说要用当前解出来的那些新的变量值的意思。
对于SOR来说，则引入了松弛系数$w$，迭代过程是
$$x_{k+1} =(wL+D)^{-1}[(1-w)Dx_k-wUx_k]+w(D+wL)^{-1}b \tag{41}$$
上边的等式其实是用了$w(D+L+U)x=wb$转换而来。
这些方法都有其收敛条件，在数值优化中有分析。
对于稀疏矩阵，或者是对称正定矩阵来说，还是一般使用本篇总结的梯度迭代方法。只不过值得留意的却是，$Ax=b$竟然对应于一个二次型的极值求解过程

Precondition

在数值分析里边，矩阵迭代求解的收敛速度取决于矩阵条件数：

$$cond(A)=||A|| \cdot ||A||^{-1} \tag{42}$$
其中$||A||$是矩阵范数，数值分析里边研究了矩阵迭代的收敛很大程度上取决于条件数，而预条件方法就是为了减少条件数，快速收敛。详细的分析很复杂，详情可看数值分析教材，这里只阐述用于共轭梯度算法的预条件过程。
预条件是希望找到一个$M$，将问题转化为求解：
$$M^{-1}Ax=M^{-1}b \tag{43}$$
其中$M$是$n \times n$的可逆矩阵，称为预条件矩阵。在改进的共轭梯度方法中，这个矩阵$M$是对称正定的。这个矩阵试图对$A$逆转，从而让$A$的条件数降低。此外，还引入了一种广义$M$内积$(v,w)_M=v^TMw$来取代欧几里得内积，所有内积默认都改成$M$内积，原始的共轭梯度方法仍然成立，因为矩阵$M^{-1}A$相对于新的内积仍然是对称正定矩阵：
$$(M^{-1}Av,w)_M=v^TAM^{-1}Mw\\ =v^TAw=v^TMM^{-1}Aw=(v,M^{-1}Aw)_M \tag{44}$$
那么，由于CG的推到基本都是基于$v^TAw$的形式，如果把等式44的内积结果应用到所有推导的过程，就有：
$$z_k=M^{-1}b-M^{-1}Ax_k=M^{-1}r_k \\ \alpha_k=\frac{(z_k, z_k)_M}{(d_k, M^{-1}Ad_k)_M} \\ x_{k+1}=x_k+\alpha_kd_k \\ z_{k+1}=z_{k}-\alpha_kM^{-1}Ad_k \\ \beta_k=\frac{(z_{k+1},z_{k+1})_M}{(z_k,z_k)_M} \\ d_{k+1}=z_{k+1}+\beta_kd_k \tag{45}$$
这里边一个有意思的地方在于，这个方法，把$v^TM^{-1}w$的欧几里得的内积，替换成了$(v, M^{-1}Aw)_M$，定义内积在$M$的逆空间里边，因为$M$内积仍然满足普通内积该有的性质，所以可以完全替代欧几里得内积，公式仍然不发生变化。
进一步转换可以得到：
$$(z_k,z_k)_M=z_k^TMz_k=z_k^Tr_k \\ (d_k,M^{-1}Ad_k)_M=d_k^TMM^{-1}Ad_k=d_k^TAd_k \\ (z_{k+1},z_{k+1})_M=z_{k+1}^Tr_{k+1} \tag{46}$$
对于SSOR(对称连续过松弛)预条件子来说：
$$M=(D+wL)D^{-1}(D+wU),w \in (0,2) \\ \implies M^{-1}=(D+wU)^{-1}D(D+wL)^{-1} \tag{47}$$
总结上边的过程可以得到，SSOR预条件子的共轭梯度算法的流程是：

$r_0=b-ax_0,z_0=M^{-1}r_0,d_0=z_0$
for loop until reach max_step or cost_function below threshold:
$\qquad \alpha_k=\frac{z_k^Tr_k}{d_k^TAd_k}$
$\qquad x_{k+1}=x_k+\alpha_kd_k$
$\qquad r_{k+1}=r_k - \alpha_kAd_k$
$\qquad z_{k+1}=M^{-1}r_{k+1}$
$\qquad \beta_{k+1}=\frac{z_{k+1}^Tr_{k+1}}{z_k^Tr_k}$
$\qquad d_{k+1}=z_{k+1}+\beta_{k+1}d_k$
return $x_n$

对应的Python代码：

import numpy as np
def get_ssor_precondition_matrix(A, w):
    UD = np.triu(A)
    LD = np.tril(A)
    dim = A.shape[0]
    D = np.mat(np.zeros((dim, dim)))
    for i in range(dim):
        D[i, i] = A[i, i]
    for i in range(dim):
        for j in range(i+1, dim):
            UD[i, j] = w * UD[i, j]
    for i in range(dim):
        for j in range(0, i):
            LD[i, j] = w * LD[i, j]
    # 对上下三角矩阵求逆矩阵，其实不必用通用的求逆方法，不停回代即可
    return np.mat(np.linalg.inv(UD)) * D * np.mat(np.linalg.inv(LD))
def precondition_conjugate_gradient(A, b, x_initial, max_step,
                                    threshold=0.00001, w=0.2):
    assert(isinstance(A, np.matrix))
    assert(isinstance(b, np.matrix))
    assert(isinstance(x_initial, np.matrix))
    r_old = b - A * x_initial
    M_inv = get_ssor_precondition_matrix(A, w)
    z_old = M_inv * r_old
    d = z_old
    x = x_initial
    for i in range(max_step):
        alpha = (z_old.transpose() * r_old) / (d.transpose() * A * d)
        x = x + d * alpha
        # r_new = b - A * x
        r_new = r_old - A * d * alpha
        z_new = M_inv * r_new
        beta = (z_new.transpose() * r_new) / (z_old.transpose() * r_old)
        d = z_new + d * beta
        r_old = r_new
        z_old = z_new
        cf = np.sqrt(np.sum(np.square(b - A * x)))
        if cf < threshold:
            print("Using step: ", i)
            break
    return x
if __name__ == '__main__':
    N = 200
    Ar = np.mat(np.random.rand(N, N))
    As = Ar * Ar.transpose()
    bn = np.mat(np.random.rand(N, 1))
    xi = np.mat(np.random.rand(N, 1))
    xr = precondition_conjugate_gradient(As, bn, xi, 1000, 0.00001, 0.05)
    print('w=0.05:', np.sqrt(np.sum(np.square(bn - As * xr))))
    xr = precondition_conjugate_gradient(As, bn, xi, 1000, 0.00001, 0.5)
    print('w=0.5:', np.sqrt(np.sum(np.square(bn - As * xr))))
    xr = precondition_conjugate_gradient(As, bn, xi, 1000, 0.00001, 1)
    print('w=1:', np.sqrt(np.sum(np.square(bn - As * xr))))

运行三次的结果：

runfile('C:/Users/zczx1/.spyder-py3/temp.py', wdir='C:/Users/zczx1/.spyder-py3')
Using step:  378
w=0.05: 6.71578034331e-06
Using step:  405
w=0.5: 9.67223926338e-06
Using step:  573
w=1: 8.09438315554e-06
runfile('C:/Users/zczx1/.spyder-py3/temp.py', wdir='C:/Users/zczx1/.spyder-py3')
Using step:  371
w=0.05: 6.1508910381e-06
Using step:  401
w=0.5: 8.51261715479e-06
Using step:  602
w=1: 7.57385104633e-06
runfile('C:/Users/zczx1/.spyder-py3/temp.py', wdir='C:/Users/zczx1/.spyder-py3')
Using step:  373
w=0.05: 6.44798434081e-06
Using step:  406
w=0.5: 6.55294163482e-06
Using step:  580
w=1: 8.13783429322e-06

虽然随机生成的矩阵不太能说明问题，但是这个只能说，在松弛系数比较小的时候，PCG和CG相比，相对来说减少了迭代的次数。但是当松弛系数较大的时候，还不如CG。

非线性共轭梯度算法

算法流程

上边的总结，都是在二次型问题上进行的分析。可以看到，虽然算法本身的实现比较简单，但是算法背后的数学思路如果进行严谨地推导的话，还是有些复杂的，不过因为有了几何的理解，所以也还是比较直观。
但是如果这个问题的表达式函数$f(x)$是非线性的，共轭梯度的思想仍然可以被推广，只是形式上不那么精确。而事实上，就算严格按照上边的算法流程，通过上边的实践也会发现，由于round off误差，其实不能达到理论的计算效果，所以在某种程度上近似求解也是可以的。回顾共轭梯度算法和核心，其实就两个：
1. 依据步长和方向向量更新解向量，即$x_{k+1}=x_k+\alpha_kd_k$
2. 利用梯度更新方向向量，即$d_{k+1}=r_{k+1}+\beta_{k+1}d_{k}$
其中，$r_k=-f^{'}(x_k),d_0=-f^{'}(x_0)$，我们只需要考虑怎么选择$\alpha_k$和$\beta_{k+1}$即可
对$\alpha_k$，是希望找到这样的步长因子使得$f^{'}(x_k+\alpha_kd_k)^Td_k=0$，本质上还是希望下一个梯度方向和当前的方向向量正交，这个时候应该就是最优的解，和SD中的线搜索是类似的。不过，在实际的实现中，可能只是设置一个固定值，然后乘以一个比例，这个比例随着迭代的进行可能会变得越来越小。所以:

$$\alpha_k=argmin(f(x_k + \alpha_kd_k) \space or \space \alpha_k=C\times radio \tag{48}$$
而对于$\beta_{k+1}$，不同的学者提出了使用不同的公式，比如：
$$Fletcher-Reeves: \beta_{k+1}^{FR}=\frac{r_{i+1}^Tr_{i+1}}{r_i^Tr_i} \tag{49}$$
$$Polak-Ribiere: \beta_{k+1}^{PR}=max(0, \frac{r^T_{k+1}(r_{k+1}-r_{k})}{r_k^Tr_k}) \tag{50}$$
也还有很多其它的计算方法，不过显然，这些只是一种近似求解的方法。因为非线性的损失函数各种各样，精确的求解不现实。这样，非线性CG的流程如下:

$d_0=r_0=-f^{'}(x_0)$
for loop until max step or $f^{'}(x) \lt thresh$:
$\qquad \alpha_k=argmin(f(x_k+\alpha_kd_k))$
$\qquad x_{k+1}=x_k+\alpha_kd_k$
$\qquad r_{k+1}=-f^{'}(x_{k+1})$
$\qquad \beta_{k+1}^{PR}=max(0, \frac{r^T_{k+1}(r_{k+1}-r_{k})}{r_k^Tr_k})$
$\qquad d_{k+1}=r_{k+1}+\beta_{k+1}d_k$

非线性函数中最优$\alpha$的选择

寻求最优的步长因子的过程其实就是广义线搜索，根据泰勒级数展开，有：

$$f(x+\alpha d) \approx f(x)+\alpha [\frac{d}{d\alpha}f(x+\alpha d)]_{\alpha=0}+\frac{\alpha^2}{2}[\frac{d^2}{d\alpha^2}f(x+\alpha d)]_{\alpha=0} \\ = f(x)+\alpha[f^{'}(x)]^Td+\frac{\alpha^2}{2}d^Tf^{''}(x)d \tag{51}$$
$$\frac{d}{d\alpha}f(x+\alpha d) \approx [f^{'}(x)]^Td+\alpha d^Tf^{''}(x)d=0 \\ \implies \alpha=-\frac{f^{'}(x)^Td}{d^Tf^{''}(x)d} \tag{52}$$
从等式52中，可以看到二次型中求解最优$\alpha$的影子，如果是二次型，那么就是$\alpha=\frac{r^Td}{d^TAd}$的形式了，不过对于非线性问题来说，就需要计算二阶导数。二阶导数为：
$$f^{''}(x)= \begin{bmatrix} \\ \frac{\partial^2f}{\partial x_1 \partial x_1} & \frac{\partial^2f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2f}{\partial x_1 \partial x_n}\\ \frac{\partial^2f}{\partial x_2 \partial x_1} & \frac{\partial^2f}{\partial x_2 \partial x_2} & \cdots & \frac{\partial^2f}{\partial x_2 \partial x_n}\\ \vdots & & \ddots & \\ \frac{\partial^2f}{\partial x_n \partial x_1} & \frac{\partial^2f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2f}{\partial x_n \partial x_n}\\ \end{bmatrix} \tag{52}$$
其实就是求解Hessian矩阵。这个计算量可谓真的大，而且随着迭代，这个矩阵的具体数值，可能还是需要不断计算的，这个开销几乎是不可忍受的，所以严格的计算步长因子基本是不可取的。何况有些时候，可能都求不出二阶导，也无怪乎实际应用中，直接给个固定步长乘以比例了事。与其算这个矩阵那么复杂，不如多迭代几次。

非线性函数求导简介

非线性CG需要对函数进行求导，而求导的方法，目前有:
1. 公式法
2. 符号微分法
3. 自动微分法
这些根据具体场景来实现，只是非线性CG迭代流程就是如上所示的过程。
值得注意的是，非线性函数的这个过程，最终收敛到的可能只是局部极值，和初始值的选择相关。如果初始值在一个类二次型的区域，那么就可以收敛到这块区域的极值点。

总结

本文总结了本人学习EDA软件中布局算法中使用到的共轭梯度相关的内容。共轭梯度算法，应用到实际工程当中，可能还是使用非线性的部分内容，只是由于也是近似的求解，所以和随机梯度下降相比也没有太多的优势，毕竟都是和初始点相关，所以引入更多的随机性，也许可以取得更好的效果！

参考资料

1. An Introduction to the Conjugate Gradient Method Without the Agonizing Pain
2. Numerical Analysis

@fsfzp888
2019 年 01月