下载掌阅APP，畅读海量书库
立即打开
现实工程问题很多表现为非线性优化，（无约束）非线性优化问题的一般形式为
即在可行域
D
中求解全局最优点
，使对于任意
x
有
。
当可行域 D = R n 时，则为无约束非线性优化问题。
无约束非线性优化求解得到的一般是局部最优点，常用求解方法包括梯度下降法、最速下降法、牛顿法、高斯牛顿法、列文伯格-马夸尔特算法等，其共同的求解思路是从起始点 x 0 出发通过迭代公式 x k +1 = x k + t k p k 产生序列{ x i }逼近最优解，其中 p k 为搜索方向， t k 为搜索步长。
特别地，存在一类典型的小残差无约束非线性优化，称为 非线性最小二乘 （Non-linear Least Square，NLS）问题。例如无人驾驶SLAM中的定位，相机标定中利用最小化重投影误差优化相机内外参数等。
此类非线性最小二乘问题，观测样本残差（Residual） r i 与决策向量 x ∈ R n 的映射关系 r i = r i （ x ）： R n → R 通常是非线性的，例如三角函数、指数函数等。
r i = r i （ x ）： R n → R ， i =1，2，…， m
全部 m 个观测样本的残差 r i 构成残差向量 r ：
r = r （ x ）： R n → R m ， x ∈ R n
目标函数 f （ x ）采用（加权）残差平方和：
则非线性最小二乘问题就是寻优化参数
x
的某个解
，使得目标函数
f
（
x
）最小。
非线性最小二乘问题的迭代解法为：
（1）给定起始点 x 0 ，开始迭代；
（2）寻找增量Δ x k ，使得 f （ x k +Δ x k ）极小；
（3）若满足停止准则，则迭代停止；否则迭代更新， x k +1 = x k +Δ x k 。
可见，非线性最小二乘问题的关键在于确定增量Δ x k ，具体求解方法包括一阶的梯度下降法、最速下降法以及二阶的牛顿法、高斯牛顿法、列文伯格-马夸尔特算法等。
需要注意，上述一阶和二阶算法都需要假设目标函数 f （ x ）光滑可微，即函数 f （ x ）在 x 邻域可以进行泰勒展开：
将泰勒展开式（1-17）忽略二次及更高的项，只保留至一阶，则
可见只要▽ f （ x ）Δ x 为负，则目标函数 f （ x ）将是下降的。因此，选定负梯度-▽ f （ x ）作为搜索方向，有
其中，标量 t 为（最优）步长，该方法称为梯度下降法（Gradient Descent，GD）。
梯度下降法伪码如下。
算法1-1中步骤4确定了迭代方向，步骤8确定了步长。步长选择很重要，因为虽然迭代方向 p k 是局部的下降最快方向，但是如果步长 t 选择不合适也会存在问题，例如图1-1。
图1-1 迭代步长
图1-1（a）的步长较小，图1-1（c）的步长较大，步长过大和过小都效果不佳。可见步长 t 直接决定收敛速度，为此需要求解最优步长。
步骤8搜索最优步长时的目标函数值为
f （ x k +1 ）= f （ x k + t k p k ）
其中，
x
k
和
p
k
均为固定的值，因此目标函数值是单变量
t
k
的函数，所以步骤8就是寻找最优步长
t
k
的最优值
。
这是一维搜索问题，求解方法很多，例如成功/失败法、0.618黄金分割法以及插值法、斐波那契法等，也可以通过解析方法
来确定最优步长
。
算法1-1中步骤5以梯度的模
作为停止准则，实际数值计算中以梯度为基准作为停止准则未必恰当，可用的准则包括：
（1）基于梯度的，
。
（2）基于目标函数的，| f （ x k +1 ）- f （ x k ）|≤ ε 。
（3）基于自变量的，
。
（4）基于上述组合而来的各类变种。
以下演示梯度下降法步骤。
具体地，设目标函数为
，超参数
ε
=0.5。
迭代第一次：
（1）设置起始点， x 0 =（0，0） T 。
（2）计算梯度，根据▽ f =（1+4 x 1 +2 x 2 ，-1+2 x 1 +2 x 2 ） T ，可得▽ f （ x 0 ）=（1，-1） T 。
进行停止准则判断，根据
，需要继续。
（3）取负梯度为迭代方向， p 0 =-▽ f （ x 0 ）=（-1，1） T 。
（4）从 x 0 =（0，0） T 出发，沿 p 0 =（-1，1） T 方向，对 f （ x k + t p k ）搜索最优步长。此时自变量为
代入目标函数
，有
f （ x ）=- t -（ t ）+2 t 2 +2（- t ） t + t 2 = t 2 -2 t
选择解析方法求解最优步长
t
0
，根据
解得
t 0 =1
（5）更新迭代。
迭代第二次：
（1）此时 x 1 =（-1，1） T 。
（2）由梯度▽ f =（1+4 x 1 +2 x 2 ，-1+2 x 1 +2 x 2 ）T，得▽ f （ x 1 ）=（-1，-1） T 。
进行停止准则判断，
，需要继续。
（3）迭代方向 p 1 =-▽ f （ x 1 ）=（1，1） T 。
（4）从 x 1 =（-1，1） T 出发，沿 p 1 =（1，1） T 方向，对 f （ x k + t p k ）搜索最优步长。此时自变量为
代入目标函数
，有
f （ x ）=5 t 2 -2 t -1
选择解析方法求解 t 1 。
（5）更新迭代。
迭代第三次：
（1）此时 x 2 =（-0.8，1.2） T 。
（2）由梯度▽ f =（1+4 x 1 +2 x 2 ，-1+2 x 1 +2 x 2 ）T，得▽ f （ x 2 ）=（0.2，-0.2） T 。
进行停止准则判断，
，满足条件，可以停止。
输出 x 2 ：
x 2 =（-0.8，1.2） T
即为
极小值的数值近似解。
以上三次迭代过程如图1-2所示。
图1-2 梯度下降法迭代过程
从图1-2可以发现梯度下降法的缺点，收敛路径在极值点附近呈现锯齿形状：开始时目标函数下降较快，接近极值点时收敛速度变慢。特别地，当目标函数等值线接近较为扁平的椭圆时收敛更慢，如图1-3所示。因此，实用中通常联合其他方法并用：前期使用梯度下降法，接近极值时改用收敛更快的其他方法。
图1-3 目标函数下降速度
梯度下降法的另一个原生缺点是，某位置的负梯度方向通常只在该点邻域附近有效；从某一轮次迭代的整个过程来看梯度可能是实时变化的。为了解决这个问题，产生了最速下降法。
梯度下降法以负梯度方向-▽ f （ x ）作为目标函数值 f （ x ）下降最快方向，每次迭代时将变量 x 沿着搜索方向移动单位步长Δ x =-▽ f （ x ）或者最优步长Δ x =- t ▽ f （ x ），使得目标函数 f （ x ）逐渐收敛。但是梯度下降法假设迭代过程中梯度固定不变，而实际上梯度可能随时变化，因此需要考虑梯度的梯度，即二阶导数，也即Hessian矩阵。迭代的收敛速度取决于函数 f 的Hessian矩阵的条件数。
如果决策向量各个分量相互独立，即∀ i ≠ j ⇒∂ 2 f /∂ x i ∂ x j =0，则Hessian矩阵为对角阵，此时使用梯度下降法逼近的效果不错。但是，如果迭代过程中梯度发生了变化，梯度下降方向不再是最佳增量方向，最佳方向与基于梯度不变假设解得的增量方向存在差异，差异大小可以使用Hessian矩阵的条件数即Hessian矩阵的最大特征值与最小特征值之比来衡量（此时Hessian矩阵不是对角阵）：如果条件数并不大，则梯度下降法也能取得较好的效果，但是对于较大的条件数，由于梯度不变假设不成立，此时可以使用最速下降法或者其他方法。
最速下降法（Steepest Descent，SD）也属于一阶算法，同样基于一阶泰勒展开式：
其中， v 为增量搜索方向， v 不再要求一定是负梯度方向。最速下降法将搜索方向定义为单位范数步长内使目标函数下降最多的方向：
其中，
为某种范数（Norm），最速下降法的搜索方向受到该范数性质的限制。如果取欧氏范数
，则最速下降法基本类似于梯度下降法。
由于最速下降法需要满足下降方向上移动单位步长处于范数集合内这个条件，因此根据所用范数不同将会得到不同结果。以向量1-范数和矩阵2-范数为例，其单位球分别为正菱形体和椭圆体：
（1）向量1-范数搜索方向为-sign{▽ f （ x ）} e i ，其中 e i 为标准基向量。所以对于向量1-范数来说，最快方向是沿着坐标轴下降，如图1-4所示。
图1-4 向量1-范数
（2）矩阵2-范数
的物理意义为最大拉伸系数，即在最大奇异值对应方向上取得最大值，此时最速下降方向为关于矩阵
A
的，
p
=-
A
-1
▽
f
（
x
）。可以看出，最速下降方向Δ
x
在满足椭圆体内的条件下尽可能地在-▽
f
（
x
）方向上延伸，如图1-5所示。
图1-5 矩阵2-范数
采用矩阵2-范数时合理选取能矩阵 A 能够使计算更有效率，Hessian矩阵就是很好的选择。使用Hessian矩阵2-范数的无约束极值求解方法称为牛顿法，它是最速下降法的一种。
前述各种梯度算法，只使用到一阶微分，以下引入二阶微分。
将泰勒展开式（1-17）忽略三阶及更高阶的项，保留至二阶，则
将式（1-22）对Δ x 求导，并令其为0，可得关于Δ x 的增量方程：
从而得到迭代公式：
其中， J f 为目标函数 f 的雅可比矩阵， H f 为目标函数 f 的Hessian矩阵，该方法称为牛顿法（Newton Method）。牛顿法用到二阶导数，需要计算Hessian矩阵。
以下演示牛顿法的算法步骤。
具体地，设目标函数为
，超参数
ε
=0.5。
第一次迭代：
（1）设置起始点， x 0 =（0，0） T 。
（2）此时▽
f
（
x
0
）=（-2，0），
，需要继续迭代。
函数 f 的雅可比矩阵和Hessian矩阵为
计算可得
（3）更新迭代。
第二次迭代：
（1）此时 x 1 =（1，1） T 。
（2）函数 f 的雅可比矩阵和Hessian矩阵为
此时的▽
f
（
x
1
）=（0，0），
，满足条件。
输出 x 1 ：
x 1 =（1，1） T
即为
极小值的近似解。
本例中，因为 f 是二次函数，二阶导数 H 为常数矩阵，所以从任意点出发仅需一步即可获得精确解。对于非二次的目标函数 f ，得到的是近似解。
牛顿法典型缺点在于：
（1）如果初始点选取不当则可能不收敛于极小点，因此有必要添加控制步长的逻辑，即增加步长项♣，使得
x k +1 ← x k -♣ H -1 （ x k ） J （ x k ）
（2）Hessian矩阵未必可逆，对此可以添加正则项 μ k ，使得
x k +1 ← x k -（ H -1 （ x k ）+ μ k I ） T J （ x k ）
牛顿法收敛较快，但是计算Hessian矩阵需要求解二阶导数，如果问题的规模较大，则需要大量运算。为了避免计算Hessian矩阵，产生了高斯牛顿法和列文伯格-马夸尔特算法。
高斯牛顿法（Gauss-Newton Method，GN）也是二阶方法，但是并不直接计算Hessian矩阵，而是通过雅可比矩阵来近似Hessian矩阵，从而只需要计算一阶导数。
高斯牛顿法步骤为：
可见高斯牛顿法最重要的是通过增量方程 H r Δ x = g r 求解Δ x ，以下从三个角度分析。
不考虑整体的目标函数 f （ x ），而是单独将第 i 项残差 r i = r i （ x ）： R n → R 进行泰勒展开，假设二阶及更高阶的项可以忽略，保留至一阶：
则目标函数 f （ x ）为
其中， J r 为残差向量 r = r （ x ）： R n → R m 的雅可比矩阵。
将式（1-26）对Δ x 求导，并令其为0，得到增量方程：
引入两个助记符号：
即增量方程为
即得高斯牛顿法的迭代公式：
高斯牛顿法也可以直接从整体的目标函数 f （ x ）得到。
考查目标函数 f （ x ）的一阶导数▽ f ：
其中， J r （ x ）为残差向量 r = r （ x ）： R n → R m 的雅可比矩阵，从而
再由目标函数 f （ x ）的二阶导数
得到
假设其中二阶项
很小，可以忽略，有
从而 f （ x ）： R n → R 的二阶 H f 可以通过 r （ x ）： R n → R m 的一阶 J r 近似得到。
至此也得到高斯牛顿法的迭代公式：
高斯牛顿法还可以从分项角度叙述，这样可能更有利于理解。
对于函数
f
=
f
（·），以
和
分别表示雅可比矩阵
J
和Hessian矩阵
H
中的分量元素，则对于目标函数
考查标准牛顿法时的迭代公式：
x k +1 ← x k - H -1 （ x k ） J （ x k ）
其中，雅可比矩阵 J 的元素为
即
其中，Hessian矩阵 H 的元素为
依然假设其中的二阶以及更高阶成分可以舍去，则
至此同样可得高斯牛顿法的迭代公式：
前述目标函数为向量点乘的欧氏距离 r T r ，自然地会考虑马氏距离 r T Wr ，不过马氏距离中的 W ∈ R m × m 是协方差矩阵，平差应用不需要如此复杂，可以将矩阵 W 稍微加强，假设为对角阵，这样成为决策变量 x 的权值矩阵，从而目标函数为
由于 W 为常值对称矩阵，以下不再从细节上考虑 W 和 W T 的区别。
此时，其一阶导数为
其中， J r （ x ）为多元向量值函数 r （·）： R n → R m 的雅可比矩阵，即
此时二阶导数为
即
假设其中二阶项可以忽略，则
至此可得迭代公式：
高斯牛顿法之所以舍弃Hessian矩阵中的二阶导数，主要原因在于牛顿法中Hessian矩阵的二阶项难以计算，或者计算的耗费较大，而计算梯度时已经获得现成的 J ，因此为了简化，采用一阶逼近二阶。
注意： 这种近似存在前提条件即残差项 r i （ x ）接近于零，或者接近于线性函数，从而▽ 2 r （ x ）接近零。也就是说，高斯牛顿法主要解决所谓的小残差问题。
一旦Hessian矩阵 H 中被忽略的二阶项并非很小，则高斯牛顿法收敛将会变慢，这种收敛速度变慢的情况在当前步骤存在较大残差项时更是经常出现。例如，待优化的参数向量接近鞍点位置，高斯牛顿法近似就不会十分准确。为了避免这种情况，需要进行合理的参数化，以及去除错误约束，或者使用鲁棒核函数（损失函数），以减小极端残差对优化过程的影响。
总结以上，高斯牛顿法的优点是：
（1）易于实现，形式简单；
（2）对于小残差问题，即 r 较小或接近线性，具有较快的局部收敛速度；
（3）对于线性最小二乘问题，可以一步到达极值点。
当然高斯牛顿法也存在缺点，主要有：
（1）对于并非特别严重的大残差问题，局部收敛速度较慢；
（2）对于残差很大或者 r 的非线性程度很大的问题，算法可能不收敛；
（3）每次迭代时要求雅可比矩阵满秩，否则算法是未定义的；
（4）缺乏控制步长的逻辑，对病态雅可比矩阵也比较敏感。
对于上述缺点部分，可以考虑增加线性搜索策略，保证目标函数每一步下降，从而对于几乎所有非线性最小二乘问题，都具有局部收敛性及总体收敛，此即所谓的阻尼高斯牛顿法。具体为
其中， μ k 为一维搜索因子。
另外需要注意高斯牛顿法的停止规则是 f （ x k +1 ）- f （ x k ）＜ ε 。
前述高斯牛顿法通过一阶雅可比矩阵 J r 合成得到二阶Hessian矩阵 H f ，考虑以下问题：
（1）即使Hessian矩阵 H f 正定，但是泰勒展开仅在 x 邻域成立，如果解得Δ x 较大，则已经超出模型假设，这个增量结果是不足信的。
（2）Hessian矩阵未必正定，也就不一定有逆矩阵。工程应用中计算所得的
H
f
常常只是半正定的，因为可以证明
至少是半正定的。此时如果用高斯牛顿法，不可逆的严重性稍微缓解，但是对于其他几种算法影响比较大。
（3）即使Hessian矩阵 H f 存在逆矩阵，依然难免出现数值计算上的病态，此时Δ x 对数值敏感，不稳定，导致逼近方向出现偏差，严重影响优化方向，算法未必收敛。
针对上述问题，逐步出现了各种改进方法。
式（1-17）的泰勒展开仅在较小邻域成立，自然地，会考虑为增量Δ x 限制一个范围，这个范围称为信赖域（Trust Region，TR），此类方法统称为信赖域法。
将高斯牛顿法的目标函数
通过系数矩阵 D ∈ R n × n 和信赖半径 π 约束Δ x 范围，有
则成为约束型优化问题：
即得到列文伯格-马夸尔特（Levenberg-Marquardt，L-M）法。
算法1-3的关键在于步骤2，这是一个带约束的平差问题，通过拉格朗日乘子 λ ，构造拉格朗日方程：
可以转换为无约束优化问题：
min f （ x ）
考查该问题的目标函数 f （ x ）：
右端对Δ x 求导，并令其为0，得到关于增量线性方程：
即
则得迭代公式
所得结果比高斯牛顿法结果多出了正则化项 λ D T D ，从而一定程度上避免二乘问题线性方程组系数矩阵的非奇异问题和病态问题，提供更稳定、准确的增量Δ x 。
以下对其中的 λ 和 D 两个参数做以下说明。
迭代公式中乘子 λ 的作用为：
（1）如果乘子 λ 较大，对角矩阵 D T D 占据主导，说明二阶近似不佳，此时算法接近梯度下降法；
（2）如果乘子 λ 较小，Hessian矩阵 H f 占据主导，说明二阶近似模型在该范围内工作较好，此时接近高斯牛顿法。
信赖域法中的系数矩阵 D 主要基于经验选择：
（1）列文伯格的方法，将 D 取为单位矩阵，对三维空间就意味着增量Δ x 被限制在一个圆球内，只有圆球中的增量才被认为有效；
（2）马夸尔特的方法，将 D 取为非负对角阵；
（3）一般工程应用中，
D
可以取方阵
对角元素的平方根，如此增量Δ
x
被限制在一个椭圆球内。
列文伯格的方法将 D 强化为单位矩阵，自然地 D T D 也是单位矩阵，因此可以用 μ I 代替 λ D T D ，其中标量 μ 称为阻尼系数，对于大的Δ x 起到惩罚作用，这样从信赖域法得到阻尼法（Damped Method）。
此时，无约束问题为
目标函数为
增量方程为
即
迭代公式为
显然，这就是前面1.3.4节中提到的带有阻尼系数的高斯牛顿法。
迭代公式中阻尼系数 μ 作用为：
（1）对于 μ ＞0，使得 H f + μ I 正定，保证Δ x 为梯度下降的方向；
（2）对于很大的阻尼 μ ，Δ x ≈- μ -1 g f ，此时接近梯度下降法，此种情况一般出现在距离最终结果较远时；
（3）对于较小的阻尼
μ
，
，Hessian矩阵占主导，二阶近似模型在此范围效果较好，此时接近高斯牛顿法。
列文伯格控制增量的重点在于阻尼系数 μ ，因此合理选择阻尼系数 μ 很重要。
初始时，如果起始点 x 0 能够较好地近似结果， μ 0 应该选取一个较小的值，例如1e-6；否则可以选择1e-3或者更大，例如1。此外，还有一种策略将初始值 μ 0 取为Hessian矩阵的最大特征值：
至于迭代过程中阻尼系数 μ 的调节，关键在于控制信赖域。最直观的就是依据实际函数和近似函数之间的差异：差异较小，说明泰勒展开式近似效果比较好，信赖域可以适当放大；如果差异较大，则需要缩小信赖域。因此定义增益率（Gain Ratio）为
其中，
f
为真实函数，分子即真实目标函数
f
经历了增量Δ
x
后的变化情况；
为泰勒展开近似的函数，分母即近似函数经历了增量Δ
x
后的变化情况。分母部分可以通过一阶泰勒展开得到：
即
分子越大，
f
下降很多，
ρ
就越大，说明此时泰勒展开近似是比较准确的，应该减小阻尼系数。分子越小，
f
下降较少，
ρ
就越小，说明增量Δ
x
太大，超出邻域很多，此时泰勒展开
不能成为一个好的近似，需要增大阻尼系数。
以下为两种调节阻尼的推荐方案。
（1）马夸尔特的调节方案。
马夸尔特的方案中：
①如果增益率 ρ 接近于1，近似效果较好；
② ρ 越小说明实际函数比泰勒展开下降得少，需要增大阻尼系数 μ ，缩小信赖域；
③ ρ 越大，说明实际函数比泰勒展开下降得多，需要减小阻尼系数 μ ，加大信赖域。
（2）尼尔森的调节方案。
尼尔森的调节方案中，对于错误的增量即 ρ ＜0时，也就是函数值 f 不仅不降反而增大，因此迅速增大阻尼系数 μ ，使得下一步Δ x 变得更小。
两种调节方案对比：尼尔森的调节方案比马夸尔特的调节方案增加一个参数 ν 。其中常数的设置并不敏感，主要是选择合适的值避免振荡。
以下为使用尼尔森阻尼方案的列文伯格-马夸尔特法。
算法1-4的终止迭代有三个条件：
（1）梯度判据，
是一个很小的正数；
（2）增量判据，
，即在
x
数值较大时使用前者（相对量）作为停止条件，在
x
数值较小时使用后者（绝对量）作为停止条件；
（3）最大迭代次数， k ≥ k max 。
L-M算法比高斯牛顿法更健壮，即使初始点距离局部最优点较远也能够到达最优点，不过，收敛速度方面要比高斯牛顿法稍慢一些。