最速下降法matlab程序_共轭梯度法（一）：线性共轭梯度

本文开始介绍一种全新的优化方法——共轭梯度法。

最初，共轭梯度法是用来求解线性方程

的，称为线性共轭梯度法。后来，有人把这种方法扩展到了非线性优化问题中，称为非线性共轭梯度法。本文从基本的线性共轭梯度讲起，介绍它的思想和算法流程。

问题描述

共轭梯度法并不适用于任意线性方程，它要求系数矩阵对称且正定。我们会在后面逐渐提到为什么会有这样苛刻的要求。

比较巧的是，在非线性优化问题中，所有的二阶方法的每次迭代都需要求解一个形如

这样的方程。比如，前面提到的牛顿系列方法，需要求解

，SLAM中常用的高斯牛顿法，需要求解

，这些系数矩阵

和

通常都满足对称正定的要求（如果不正定，则想办法让它们正定，比如加上一个单位阵）。

之前，我们都默认直接对系数矩阵求逆，得到

，而没有考虑求逆的性能代价。即便是通过矩阵分解的方式求逆，计算量仍然很大。线性共轭梯度法的提出，让我们能够通过迭代的方式求

，而不必对系数矩阵求逆。

具体如何实现呢？我们先额外构造一个二次函数

如果对该函数求导，并令导数为零，我们得到

你会发现，使得该函数导数为0的

恰好是线性方程

的解。换句话说，

的极小值点（

是正定矩阵，意味着

有极小值）恰好是线性方程

的解。于是，我们可以把求解线性方程的问题转化为求二次函数极值的问题。

寻找思路

如何求二次函数极值？可以用之前讲过的最速下降法、牛顿法之类的方法吗？

显然，牛顿法我们不用考虑了，因为刚才提到，求解线性方程其实是牛顿法的每次迭代都要面对的问题，所以它们不是一个层次上的问题。最速下降法倒是可以考虑，我们只需要求解当前位置的导数

，然后在该方向上确定一个最优的步长，不断迭代即可。回顾一下线搜索（二）那篇文章中给出的最速下降法的收敛路径，如下图所示。

对于一个简单的二次函数，最速下降法竟然走出这么蜿蜒曲折的路径，这要用在求解线性方程上，恐怕比直接矩阵求逆还要慢。

我们仍然以上图所示的二维目标函数为例，最好的结果是走出下面这条路径。

我们这样解读这条新的路径：整个路径只有两段，即两次迭代，对应着目标函数只有两维。第一次迭代从

到

，沿着坐标轴

的方向前进，达到了目标函数在该方向上可能达到的最优值。第二次迭代从

到

，沿着坐标轴

的方向前进，达到了目标函数在该方向上可能达到的最优值。两次迭代后就达到了极小值点。

重新审视这一路线，可以发现，每次迭代的前进方向是互相正交的，因为都是沿着坐标轴前进的。但最速下降法的搜索路径也是互相正交的，为什么就不行呢？如果我们把目标函数稍微旋转一下，就会发现问题的关键。

这是另一个二次函数，它与前面的二次函数很不一样，椭圆的轴不与坐标轴平行。（这里补充一点，之所以我们画出的二次函数的图形都是椭圆，是因为我们只考察系数矩阵为对称阵的情况。）此时，仍然沿着坐标轴搜索，这就是不折不扣的最速下降法了。这说明，能不能在两次迭代内达到极小值，与是否沿着坐标轴搜索没什么关系，关键在于椭圆的轴是否与搜索方向一致。在第一种情况中，椭圆的轴与坐标轴平行，搜索方向也和坐标轴平行，所以两次迭代即可。第二种情况中，椭圆的轴不与坐标轴平行，搜索方向依旧与坐标轴平行，所以两次迭代是不够的。这给了我们一些提示，如果我们先算出椭圆长短轴的方向，然后沿着这些方向去搜索，就能在两个迭代后达到极小值。

从上面的二维函数推广到N维函数，想象N维空间中的椭圆，我们仍然可以计算出每一维的轴向，然后沿着这些轴搜索，就能确保在N次迭代后达到极小值。那如何求椭圆的轴呢？很简单，椭圆的轴其实就是系数矩阵的特征向量，对

做特征值分解就可以了。说到这里，其实有点怪怪的，特征值分解的计算量恐怕不小吧，为了避免求逆却整出来个特征值分解，岂不是舍近求远。

那么，有没有不求特征值，也能在N次迭代达到极小值点的方法呢？比如，能不能走出下面这样的路径。

第一次迭代前进的方向并不与椭圆的轴平行，在这种情况下，第二次迭代却能一次达到最优点。异想天开？当然不是，共轭梯度法就是要走出这样的路线。

我们用数学工具探究一下这种路径需要满足的条件。对于二次函数中的任意一个点来说，该点与函数极值点的连线正是牛顿法中每次迭代的线搜索方向。回顾线搜索与信赖域中的式(2.14)，套用这里的系数矩阵，可以得到图中的向量

。如果假设起始点

的搜索方向是最速下降方向，即

，于是有

其中，最后两个梯度的乘积等于0是因为最速下降法的每次迭代方向都是垂直的。我们从二维例子中推出了一个看起来有些意义的结果，这个结果表明，连续的两次迭代的前进方法关于系数矩阵

正交。也就是说，其中的一个方向经过

的变换后与另一个方向正交。

接下来的问题是如何把这一结论推广到N维目标函数中去，我们先不关注证明，直接看一下结论。

共轭梯度法

对于一组向量

，如果任意两个向量间满足

则称这组向量与对称正定矩阵

共轭。

可见，刚才我们发现的正是共轭的性质，只不过，推广到高维空间后，任意两个向量之间都要关于

正交，不仅仅是相邻的。

共轭向量的好处正如我们前面发现的，依次沿着每个向量优化，最终可以在N个迭代后达到极小值。书中有严格的证明，但我们无意在此着墨，我认为结合上面的例子便可以理解其中的思想。

现在的问题是，给定任意一个对称正定矩阵，如何求取一组共轭向量？

显然，

的特征向量就是一组共轭向量，因为

对于

的任意两个特征向量

和

都成立，满足式(5.4)，这就从公式上验证了我们前面的发现。不过，之前也说了，实际中是不能用特征向量的，因为计算太复杂。（请注意，这里我们用到了

这一性质，只有对称矩阵有此性质，呼应了文章开头提到的对矩阵

的特殊要求。这一性质贯穿了整个共轭梯度法的证明过程。）

实际中，需要一种不依赖于矩阵分解，可以在每次迭代时以较小的代价计算出

的方法。回顾上文的最后一张图，当时我们假设起始点

的搜索方向是最速下降方向。这其实不是随便假设的，这样做会使得后面的计算变得简单。接下来，我们就从

开始，推导出后续每次迭代所需的

。

满足式(5.4)的共轭向量有无数个，而我们只需要取其中的一个。在

的时刻，我们已知

和

，注意到这两个向量一定是线性无关的，因为

是该方向上使得目标函数最小的步长，所以在该方向上目标函数不可能继续下降。那么就可以在这两个向量构成的平面上生成新的

，令

，将

和

带入式(5.4)，即可求得

，从而得到

。

在此基础上，结合数学归纳法，我们可以证明，令

其中

即可使得生成的所有

满足式(5.4)的要求，从而得到一组共轭向量。有了每次迭代的方向

，就可以计算出

在该方向上的最优步长。计算方式很简单，我们用

替换

中的

，得到

，将其对

求导，令导数等于0，即可得到最优步长

于是，整个共轭梯度法的流程整理如下

其中，用

替代了

。从

开始，

，依次迭代下去。

实际中，为了降低计算量，一般会对式(5.13)稍做改动，得到下面的算法流程

这一版本的共轭梯度算法涉及的矩阵乘法和加法会更少一些。

收敛速度

优化算法的收敛速度是个值得关心的问题。前面提到了，共轭梯度法一定会在N次迭代内收敛到极值。但在N次迭代的过程中，收敛速度怎样呢？

我们用线搜索（二）：收敛性和收敛速度中提到过的向量的矩阵范数来衡量共轭梯度法的收敛速度。比如，评估

随着

的增加的变化率。这里省略证明过程，直接给出结论

其中，

是矩阵

的第

个特征值，

是任意一个关于

的

阶多项式。这个等式可以理解为，算法在第

次迭代后，残差的最大值与一个

阶多项式的值有关，这个

阶多项式是所有

阶多项式里面结果最小的那个，但每个多项式的输入必须是使其取得最大值的

。单从这个不等式其实是看不出什么东西的，但却可以推导出一个惊人的结论。

如果矩阵

的N个特征值只能取

个互不相同的值，这

个值分别是

。我们定义一个

阶多项式

可以发现该多项式满足

。我们再定义另一个多项式

此时有，

，

是它的一个根。而且，

是一个

次多项式。如果我们把这个多项式带入到式(5.32)中，可以得到

其中，第二行用特定的

阶多项式

代替了任意的

阶多项式

，相当于放大了表达式的值。第四行，根据

的定义，输入任意的

，

都等于0。因此，我们最终发现，在第

次迭代后，

就已经收敛到了极小值

。

这是个非常好的性质，它说明，共轭梯度法在N次迭代后收敛其实是最坏情况下的结果，当特征值存在重叠时，只需要更少的迭代次数就能收敛。

不过，以上还没有回答我们想要了解的问题，在迭代过程中，算法的收敛速度如何。下面的式子在一定程度上给出了回答（这里仍然省略证明）

在第

次迭代后，残差不大于

。这个结论也是很有价值的，这里的

是从小到大排列的特征值，

越小，残差就越小。对于同样的

次迭代，如果矩阵

的第

大的特征值比矩阵

的第

大的特征值小，矩阵

就会收敛得更快。另一方面，也可以理解为，整个收敛过程中的收敛速度是变化的，与特征值的大小分布有关。如果特征值的取值呈聚类状分布，比如有5个较大的特征值，还有几个接近于1的特征值，那么迭代过程很可能像下图这样

在7次迭代后就达到了比较小的值。这是因为除了5个较大的特征值，其它特征值都集中在一个聚类中。根据前面得出的结论，如果

只有

个不同大小的特征值，

次迭代后就能收敛。现在的情况与之类似，虽然聚类中的特征值并非完全一样，但我们仍然可以将其当做近似相等的特征值，所以优化在6次迭代后理论上就会接近于极值。当然，从图中可以看到，实际情况下，第7次迭代才接近极值，这也是合理的。

如果我们进一步对式(5.34)放缩，可以得到一个更宽松的收敛性质

在该式中，收敛速度受

的条件数制约，条件数越小，收敛速度越快。在矩阵的条件数中我们介绍了条件数的定义和计算方式，即

。这说明，特征值分布越集中，收敛速度越快，分布越分散，收敛速度越慢。极端情况下，当

时，所有特征值都相等，

的条件数为1，此时只需要一次迭代即可收敛到极值点。

既然共轭梯度法的收敛速度如此依赖于特征值分布，我们理应做点什么，以改善那些特征值分布不佳的情况。

预处理策略

最容易想到，也是最有效的方法就是，对

做个变换，让变换后的

具有更为集中的特征值。比如，令

针对

的优化目标函数为

其中，新优化方程中的系数矩阵为

。如果我们找到合适的

，使得

具有比较集中的特征值以及较小的条件数，那么对

的优化就能很快收敛。

但问题的关键在于如何选择合适的

。极端情况下，如果令

，可以发现，

恰好是

的Cholesky分解，即

。但我们不能这样做，因为对

做Cholesky分解无异于直接求

的逆，计算量太大。一种更快的方案是Incomplete Cholesky分解，该分解得到使得

的稀疏矩阵

，该矩阵可保证

，从而得到良好的特征值分布。

考虑预处理之后，式(5.23)转化为如下算法

从

开始，

，依次迭代下去。不过，预处理虽然改善了收敛速度，却多了求解

的步骤，会增加一部分计算量。

总结

本文介绍了线性共轭梯度法，这是一种以迭代的方式求解线性方程的方法。这种方法的妙处在于，不需要矩阵分解和求逆，只需要上次迭代的步长和当前的梯度，即可计算下次的步长。计算量极低，且不占用额外的内存空间，适用于求解大规模的线性方程。

参考资料

Chapter 16 Preconditioning Illinois math course

上一篇：信赖域（二）：确切方法
下一篇：共轭梯度法（二）：非线性共轭梯度