深度解读：最小二乘问题的一系列优化算法及代码实现（一）

2023年伊始，带领团队攻关公司一个核心项目：高精度人体3D骨架及数字人，涉及到大量关于数值最优化方面的知识，其中最典型的是关于最小二乘问题的求解方法，以及如何将一个问题转换为最小二乘问题。这里专门写几篇博文进行总结，同时也会结合著名的开源非线性优化算法库Ceres进行对应说明，有助于大家在以后的工作中遇到类似问题后，对于Ceres中的一些选项有更深入的理解，如有不足或错误之处，还请多多指教。【下面链接是在线实时与微软3D骨架算法的对比效果，如果觉得我们的3D骨架算法效果还可以，还请您高抬贵手，点个大大的赞，非常感谢！】

首先，我们来看看最小二乘的定义：通过最小化代价函数 来估计观测数据的数学模型的最优参数。举个例子，假设我们获取 个观测样本 ，通过构建适当的数学模型(函数) ，并估计最佳的 参数，使得尽量满足这批样本的分布形式，整个过程可以描述成如下形式：

根据构建的 形式，可以分为线性最小二乘问题(比如： 和非线性最小二乘问题。我们身边很多问题都可以转化为最小二乘问题，比如各种形式的曲线拟合，视觉SLAM中前端视觉里程计关于机器人位姿( 、 )估计。这里顺带讲一个问题，为什么我们最小化残差的平方和呢？为什么不是绝对值之和？或者四次方之和？直接原因是平方操作不显著增加多项式次数，且处处可导。而绝对值操作在原点不可导，用四阶也是可以的，只不过好像没那个必要，因为你对一个数进行平方后再平方，给人的感觉是画蛇添足。

如何求解上述问题呢？一个很自然的方式是求 关于 的一阶导数 ，并令 ，可以得到多个极值点(注意：也有可能是鞍点)，然后逐个比较其值大小即可。但是，很多时候求解 并不是意见容易的事(如果 的函数形式比较简单，比如线性函数，则可以直接求解)。换个角度，既然很难从全局维度求解方程，那如果能从某一个初始值 开始，不断更新当前的待优化的变量 ，让 在局部区域内一直下降，直到最终 下降幅度小于一定值，则认为到达了极小值，也就得到了优化变量的值。所以总体思路如下：

【step1】：给定某个初始值 ;
【step2】：对应第 次迭代，计算增量 ，并使得 达到局部极小值；
【step3】：如果 小于设定阈值，则停止迭代，否则令 ，并返回step2。

关于上述思路，我们就把一个求解方程的问题转化为不断求解一个 使得 下降的问题，从更深层次的角度可以认为是把一个全局问题转化为局部问题。因为，求解 并不关注函数局部分布形式(也就是说我们并不需要一点点追溯函数局部值)，但是迭代求解，每次都是在函数的局部区间内寻找一个能使 下降的增量，使得我们只需要关心 在迭代值处的局部性质而非全局性质，这种转换思路是一种非常常见的最优化思想。

为了方便后续代码演示，这里我们举一个实际的曲线拟合例子，考虑如下方程：

其中， 、和 为待估计的曲线参数， 为故意引入高斯噪声。现在假设我们有 个观测样本 ，想依据这组观测数据估计出最佳的 、 和 参数，构建出的最小二乘问题如下所示：

在讲解各种优化方法之前，首先对优化方法的类型做一个大概介绍，参考《Numerical Optimization》所述和上述介绍，可以发现一个共同点，均需要给出一个初始值，从这个初始值开始，寻找一个变化量(增量 )，使得代价函数逐渐下降，然后达到最优。可以发现，在整个过程中，如何寻找增量是其最重要的步骤。增量是一个向量，涉及到大小和方向，也就是说，我们需要将大小和方向都确定后，才能进行迭代计算。纵观优化算法领域发展历程，均是在围绕如何求解出合理的增量，因而衍生出两种策略，线性搜索方法和信赖域方法，Ceres中的枚举类型MinimizerType也对应给出LINE_SEARCH和TRUST_REGION。

线性搜索方法的总体策略是：先确定优化变量的更新方向，然后在该方向上确定一个能使代价函数下降最大的步长。那如何确定其方向呢？这里先不做具体展开讲解，在后面的梯度下降法、牛顿法、共轭梯度法以及高斯牛顿法等等详细说明，这些方法都充分利用到了函数的一阶梯度和二阶梯度信息。在Ceres库中枚举类型LineSearchDirectionType列出了诸如：STEEPEST_DESCENT、NONLINEAR_CONJUGATE_GRADIENT、LBFGS、BFGS四种线性搜索方向优化方法。假设现在已经确定好优化方向 ，接下来就是要确定步长 ，在Ceres库中枚举类型LineSearchType列出了两种方法：ARMIJO、WOLFE。接下来具体讲解一下这两种方法，相当于求解下式：

下图是随着 变化，代价函数的可能变化趋势(注意：我这里说的是可能，实际上可能会让你崩溃的不知所措)：

最暴力的确定步长的方法就是从0开始，逐渐增大 值，然后计算代价函数值，直到找出能使代价函数最小的 值。如果你的计算机性能足够高或者不要求优化速度，这么操作问题不大。但是，数学家们是不能忍受的，毕竟这么干好像不需要高大上的数学操作，但是这种暴力搜索试的次数太多，效率肯定低。随着 的增大，实际上代价函数的变化趋势是起伏不定的，如下图所示：

上图中的 是略微向下倾斜的直线，斜率与初始位置的导数成正比。代价函数高于虚线的部分是不可能存在极小值的，以此将步长分为“可接受（acceptable）”的区域和“不可接受”的区域。公式表达如下：

其中 为Jacobian Matrix。这里容我延伸一下关于Jacobian Matrix的一些知识，一句话概括：Jacobian Matrix 是函数在其任意一点邻域内的线性近似，其本质是利用一阶泰勒展开公式。关于雅克比矩阵更多的精彩讲解可以参看文献【1】和【2】。为了图文对应，这里说明一下：

上式中的 是虚线斜率的缩放比例，一般取 左右，这意味着虚线l的下降速度远慢于代价函数初始阶段的下降速度，因而该约束限制的区域足以保证代价函数获得充分的下降，因此称为充分下降条件（sufficient decrease condition），又称Armijo Condition。可以发现，Armijo condition这个约束过于松弛，也就是说acceptable区域太宽泛，任意小的步长 都可以满足该条件，最终可能导致优化算法的效率并不高。

想提高优化算法效率，就需要进一步压缩acceptable区域，且需要压缩到能使代价函数下降比较大的acceptable区域之内。基于此想法，可以想到，极小值点只可能出现在导数比较小的区域，对于导数比较大的acceptable区域其实可以直接不考虑，因为代价函数正在处于快速下降阶段，如果在此处取 值，就相当于使下降早停，我们更希望找一个能合适的 值，使代价函数的降低到导数平缓区域，也即极小值附近。如下图所示：

针对上述分析，我们可以加入一个曲率条件(Curvature Condition，其实就是用导数大小来衡量)，如下：

是虚线desired slope斜率的放缩比例，针对不同算法取值有差异，比如牛顿法和拟牛顿法一般取0.9，对于共轭梯度法取0.1。将Armijo condition和curvature condition相结合，取acceptable区域的交集，即为Wolfe condition。这里顺便提一下，在Ceres官方文档中清楚说明：为了保证满足BFGS和LBFGS线搜索方向算法的假设，应该使用WOLFE线搜索算法。Wolfe conditions具有尺度不变性，当代价函数存在尺度缩放时，这些约束条件也会跟随函数值和梯度值的变化而变化，可接受区域的计算结果保持不变。

Wolfe Condition需要实时计算导数值，可能会降低算法效率，如果只依赖于代价函数本身则会简单很多。因此，提出了Goldstein Condition，其核心思想就是：将acceptable区域限制在两条直线之间，公式如下：

其中： 。这里解释一下，上式右边部分等价于等价于Armijo Condition，保证能充分减小，而上式左边部分的引入主要是为了从下面控制步长范围，也就是下图中的acceptable steplengths。

如《Numerical Optimization》中3.1章节所说，Goldstein条件相对于Wolfe条件的一个缺点是(7)式左半部分可以排除代价函数的所有极小值。然而，Goldstein和Wolfe条件有很多共同之处，它们的收敛理论也非常相似。Goldstein条件常用于牛顿型方法(诸如牛顿法、高斯牛顿法)，但不适用于维持正定Hessian近似的拟牛顿方法中(Quasi-Newton：BFGS、LBFGS等)。

在满足Armijo、Wolfe还是Goldstein Condition条件下，我们该如何进行步长搜索呢。步长搜索最简单的办法是回溯法(实质就是暴力搜索法)，回溯法既不采用Goldstein conditions，也不采用Wolfe conditions，而是只考虑充分下降条件。从步长 为1开始，测试是否满足充分下降条件，如果满足则停止，否则将步长缩小，再次测试，直到满足充分下降条件为止。这种方法之所以可行，是因为充分下降条件在步长足够小时一定能够成立。实际使用中，大部分情况都可以在步长为1或者少量的迭代后停止，从而获得较大的函数值下降。仔细看看整个过程，采用暴力搜索法有一个非常大的问题就是初始值选择问题，比如说 从初始值为1合适吗？会一定使得代价函数减小吗？从后面即将讲到的梯度下降法可以看出， 的初始值选择对于整个优化过程至关重要，有时候甚至需要试错很多次才能选择一个比较合适的值。

结合之前所述诸多约束条件(Armijo、Wolfe、Goldstein)，无论什么样的代价函数，可接受的步长一定落在某个区间内，我们可以先确定一个足够大的区间，然后逐渐缩小它，直到找到可接受的步长。如何确定一个足够大的区间，其实并不那么简单。实际上，如果我们知道充分下降条件代表的那条直线 与代价函数 的交点，该点就可以作为区间的右端点，直接求解这个交点并不是一件容易的事，因为很多时候代价函数形式非常复杂，比如 次多项式。既然直接求解不行，那就采用迭代办法，可以从一个较小的 开始，如果满足Armijo条件，则令 ，其中 是一个大于1的常数，直到找到一个不满足armijo条件的 。然后针对区间 ，我们逐渐去缩小它，直至最终找到可接受的步长。

我们如何去逐渐缩小呢？并最终找到一个合适的步长呢？针对此问题，Ceres中的枚举类型LineSearchInterpolationType给出了三种方法：BISECTION、QUADRATIC和CUBIC。最简单的方法就是二分法(BISECTION)，直接取中值。但是接着问题又来了，二分之后，我们保留哪半边？换句话说，我们如何保证 内部存在一个满足Wolfe conditions的子区间？让我们思考一下Wolfe conditions，它实际上规定了两件事情，一是函数值必须下降，二是只接受导数较平缓的区域。既然和导数大小有关，那么我们可以试着考虑二分点的导数值。如果二分点的导数值很小，直接满足Wolfe conditions，那就找到了可接受的步长。如果二分点的导数值很大，那么我们就在二分点和其中一个端点间继续二分。至于选择哪个端点，是整个二分法的最核心步骤。理论上，由于函数的形状千变万化，假设某次迭代的二分点为 ，其两侧 都可能存在可接受的步长。但只有一侧可以保证有，那一定就是函数值下降的那一侧。如果代价函数在 处的导数为正，那么向左侧走函数值下降；反之，如果 的导数为负，那么向右侧走函数值下降。于是，我们用 来判断 是否是右端点，通过 来判断 是不是函数下降的那一侧端点，如果成立， 不在函数下降的那一侧，则此时 变为 ，搜索区间变为 。按照这种策略，逐步缩小区间范围，直到二分点满足Wolfe conditions为止。二分法实现简单，但没有利用到更多的信息。一种更好的方法是对于区间 插值或拟合。当我们知道区间的两个端点和它们各自的函数值时，如果再知道其中一个端点的导数值(可以通过数值微分求解，在Ceres的NumericDiffMethodType中给出三种求解方式：中心差分CENTRAL、向前差分FORWARD、自适应微分RIDDERS，这里顺带提一下：尽量使用中心差分法，如果要更高准确度或无法给出目标函数确定一个合适的静态相对步长，建议Ridders方法)，就可以拟合一条二次曲线，取该二次曲线的极小值点，作为新的步长，可能会比二分查找更好一些。如果再考虑上另一个端点的导数值，就可以拟合一条三次曲线，取该三次曲线的极小值点，作为新的步长。利用到的函数和导数信息越多，拟合的曲线就越接近原函数，介于优化效率，一般最多只会拟合到三次曲线，也即拟合二次(QUADRATIC)或三次曲线(CUBIC)即可。

信赖域方法是先用一个简单的模型 近似在当前 值下的代价函数，然后初步确定一个信赖域 ，在该半径范围内寻找一个能使代价函数下降最多的增量 ，相当于求解下式：

可以认为，信赖域方法是先确定最大步长，再确定方向和实际步长。数值分析中常用泰勒公式进行展开近似(二阶泰勒近似)，如下式：

然后再选择一个信赖域半径，并在该半径范围内确定出一个最佳方向和步长。(这里容我延伸一下关于Hessian Matrix的一些知识，Hessian Matrix主要用来判断函数极值点的情况，判断梯度向量为零的点究竟是最高点、最低点还是鞍点。更多的精彩解释请参看文献【3】，搞图像处理的大佬们肯定知道Harris角点检测是就需要根据Hessian Matrix来分析)。如下图所示：

上图中，右侧的虚线是实际代价函数的等高线，实线是近似的代价函数的等高线，且最上面一条代表 =12(表示第12次迭代时的近似函数)，最下面左右对称的两条代表 =1，中间的虚线同心圆表示信赖域半径。根据近似的代价函数 ，可以在外侧的圆上找到一个最优的方向，即图中较长的 (对应本文的是 )，然后判断该 是否足以使原代价函数 产生足够的下降。如果不行，则在半径更小的圆上再次尝试，即图中内侧较短的 。可以发现，每次缩小信赖域半径，找到的优化方向都不一样。读到这里，可以发现，如何确定合适的半径 是信赖域类型优化算法的重点，我们可以思考下，什么样的信赖域半径才是合适的？一个很自然的想法就是如果我们的近似代价函数 在半径 内是接近真实的代价函数 ，那我们用 来估计 的状态，通过求解 的极小值，从而确定出最佳的 。因此，我们可以这样定义一个衡量值：

针对(10)式，假设我们在某个信赖域半径 下求出一个 ，然后我们可以利用上式计算出 的值，如果 接近于1，说明二阶近似很接近代价函数，这一步长 是可以接受的。如果 接近于0甚至小于0，说明二阶近似与真实的代价函数差距较大，我们需要减小信赖域半径 ，并重新计算步长 。反之，我们需要增大信赖域半径 ，关于如何在确定出 之后求出 ，诸如Levenberg–Marquardt、DogLeg等，在后面会详细解说。

介绍完两种优化算法类型，大家有没有思考过一个问题？为什么会出现信赖域这种方法类型，那肯定是线性搜索类型的优化算法存在一些先天性不足，然后总得有人要来填坑吧，所以信赖域类型优化算法就出现了。我们来回顾下，在信赖域方法的每次迭代中，先确定一个信赖域半径，然后在该半径内计算目标函数的二阶近似的极小值。如果该极小值使得目标函数取得了充分的下降，则进入下一个迭代，并扩大信赖域半径，如果该极小值不能令目标函数取得充分的下降，则说明当前信赖域区域内的二阶近似不够可靠，需要缩小信赖域半径，重新计算极小值。如此迭代下去，直到满足收敛所需的条件。

上图中，contours of 即为代价函数的轮廓，为了说明信赖域算法相比线性搜索的优越性，这里把代价函数形式设置的复杂一些，表现就会非常扭曲。左上方的黑点是第 次迭代的起始位置，也即 。以 为圆心，根据信赖域半径 画圆，得到左上方的圆形信赖域，我们在 处构建代价函数的近似模型：

其中 为Jacobian Matrix， 为Hessian Matrix或其Hessian的近似矩阵。(11)式中 作为变量，并使得其最小化。在二维情况下， 函数的等高线是图中的椭圆虚线，可以看到，在 附近， 的轮廓与代价函数 比较接近，但是在信赖域之外，就差得比较远了。 的极小值(上图中右边的黑点)与 的极小值相差甚远，如果此时我们采用线性搜索方法，则可以看到其搜索方向直接指向 函数的极小值点，也就是上图中的“Line Search Direction”，但是从上图中明显可以看出，如果沿着线搜索方向，几乎与代价函数等高线平齐，并不能使代价函数充分下降。如果采用信赖域方法，且任然用 作为近似函数，但是有信赖域半径的限制，最终找到图中标记的“Trust region step”向量，此向量即为使得在该信赖域范围限制下 充分下降的步长。可以发现，对比线性搜素方法，此时信赖域算法展示出他的优越性。在Ceres中的TrustRegionStrategyType中列举了LEVENBERG_MARQUARDT和DOGLEG两种信赖域算法。

从下一小节开始，正式开始讲解上面遗留的一个重要问题，如何确定增量 的更新方向，首先讲解的是：梯度下降法和共轭梯度法。

[1]. 雅可比矩阵几何意义的直观解释及应用 - mathinside的个人空间 - OSCHINA - 中文开源技术交流社区
[2]. https://www.bilibili.com/video/BV1NJ411r7ja/?vd_source=5f11e9196e87fd56a0b5843d68ab6b65
[3]. hessian矩阵的特征向量有什么含义？