非线性最小二乘问题与Levenberg–Marquardt算法详解

1 最小二乘问题

       给定一组连续函数 $\mathbf{f}:{\mathbb{R}}^n\to {\mathbb{R}}^m,m⩾n$，其最小二乘问题的定义为

x ∗ = arg ⁡ min ⁡ x { F ( x ) } ,   F ( x ) = 1 2 ∑ i = 1 m ( f i ( x ) ) 2 = 1 2 ∥ f ( x ) ∥ 2 = 1 2 f ( x ) T f ( x ) . (1.1) {{\mathbf{x}}^ * } = \arg {\min _{\mathbf{x}}}\left\{ {F({\mathbf{x}})} \right\},{\text{ }}F({\mathbf{x}}) = \tfrac{1}{2}\sum\limits_{i = 1}^m {{{\left( {{f_i}({\mathbf{x}})} \right)}^2}} = \tfrac{1}{2}{\left\| {{\mathbf{f}}({\mathbf{x}})} \right\|^2} = \tfrac{1}{2}{\mathbf{f}}{({\mathbf{x}})^T}{\mathbf{f}}({\mathbf{x}}).\tag{1.1} x∗=argxmin​{F(x)}, F(x)=21​i=1∑m​(fi​(x))2=21​∥f(x)∥2=21​f(x)Tf(x).(1.1)

该问题一个重要的应用是数据的拟合，其中 $f_i(\mathbf{x})$ 可认为是样本值与拟合值之间的残差或者距离，$\mathbf{x}$ 是需要拟合优化的参数。例如，假设有一组需要拟合的函数 $\mathbf{g}(\mathbf{t};\mathbf{x}):{\mathbb{R}}^N\to {\mathbb{R}}^M$，这时自变量为 $\mathbf{t}$，待优化参数为 $\mathbf{x}$，在这里被视为常数；同时有一组样本 $({\mathbf{t}}_k,{\mathbf{y}}_k)$，其中 $1⩽k⩽K$，那么 $f_i(\mathbf{x})$ 可表示为

$$f_i(\mathbf{x})=f_i(\mathbf{x};{\mathbf{t}}_k)=g^{(j)}({\mathbf{t}}_k;\mathbf{x})-y_k^{(j)},\text{ 1}⩽j⩽M,\text{ 1}⩽i⩽M\times K.\text{\hspace{1em}(1.2)}$$

从式（1.2）可以看出，各 $f_i(\mathbf{x})$ 的形式并不要求一致，而且可以同时优化不同的拟合函数，一个样本可能对应多个拟合函数。注意，这里实际上并不要求每个拟合函数包含所有的自变量与参数，也就是说实际的 $f_i(\mathbf{x})$ 个数可能少于 $M\times K$ 个，但一般来说应该比需要拟合的参数的数量多，这在后面会有分析。

       为了求解最小二乘问题，最直接的方法是对 $F(\mathbf{x})$ 求梯度，可得

$$\frac{\partial F(\mathbf{x})}{\partial x_j}=\sum _{i=1}^m{f}_i(\mathbf{x})\frac{\partial f_i(\mathbf{x})}{\partial x_j}\Rightarrow \nabla F(\mathbf{x})={\mathbf{F}}^{\prime }(\mathbf{x})={\mathbf{J}}_f(\mathbf{x}{)}^T\mathbf{f}(\mathbf{x}).\text{\hspace{1em}(1.3)}$$

其中 $J_{f(\mathbf{x})}$ 称为 $f(\mathbf{x})$ 的雅可比（Jacobian）矩阵，其定义为

$${\mathbf{J}}_f(\mathbf{x})=\left[\begin{array}{cccc}\frac{\partial f_1(\mathbf{x})}{\partial x_1}& \frac{\partial f_1(\mathbf{x})}{\partial x_2}& \cdots & \frac{\partial f_1(\mathbf{x})}{\partial x_n}\\ & & & \\ \frac{\partial f_2(\mathbf{x})}{\partial x_1}& \frac{\partial f_2(\mathbf{x})}{\partial x_2}& \cdots & \frac{\partial f_2(\mathbf{x})}{\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial f_m(\mathbf{x})}{\partial x_1}& \frac{\partial f_m(\mathbf{x})}{\partial x_2}& \cdots & \frac{\partial f_m(\mathbf{x})}{\partial x_n}\end{array}\right]\in {\mathbb{R}}^{m\times n}.\text{\hspace{1em}(1.4)}$$

假设 $F(\mathbf{x})$ 连续且可导，那么其极值点对应着梯度为零的位置，即需要求解

$${\mathbf{F}}^{\prime }({\mathbf{x}}^{\ast })={\mathbf{J}}_f({\mathbf{x}}^{\ast }{)}^T\mathbf{f}({\mathbf{x}}^{\ast })=0.\text{\hspace{1em}(1.5)}$$

1.1 线性最小二乘问题

       假设 $\mathbf{f}(\mathbf{x})$ 为线性函数，即可表示为 $\mathbf{f}(\mathbf{x})=\mathbf{A}\mathbf{x}-\mathbf{b}$，那么式（1.1）称为线性最小二乘问题，其求解比较简单，即

$${\mathbf{J}}_f(\mathbf{x})=\mathbf{A}\Rightarrow {\mathbf{F}}^{\prime }({\mathbf{x}}^{\ast })={\mathbf{A}}^T\left(\mathbf{A}{\mathbf{x}}^{\ast }-\mathbf{b}\right)=0\Rightarrow {\mathbf{x}}^{\ast }={\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T\mathbf{b}.\text{\hspace{1em}(1.6)}$$

注意这里基于 ${\mathbf{A}}^T\mathbf{A}$ 可逆的前提，这就要求其是满秩的，即 $R({\mathbf{A}}^T\mathbf{A})=n$。又因为

R ( A T A ) = R ( A ) ⩽ min ⁡ { m ,   n } , R\left( {{{\mathbf{A}}^T}{\mathbf{A}}} \right) = R\left( {\mathbf{A}} \right) \leqslant \min \left\{ {m,{\text{ }}n} \right\}, R(ATA)=R(A)⩽min{m, n},

因此要求 $m\ge n$，且 $\mathbf{A}$ 是列满秩的，这就要求对于每一个参数，至少有一个函数 $f_i(\mathbf{x})$ 与其相关。一般来说，我们不会通过直接求解 ${\mathbf{A}}^T\mathbf{A}$ 的逆来求解式（1.6）的线性方程组，因为其效率低且不稳定，特别是 ${\mathbf{A}}^T\mathbf{A}$ 条件数很大时。另外计算机截断精度也是一个不可忽视的问题。所以通常会通过矩阵分解，如 QR 分解或者 Cholesky 分解等方法来求解线性方程组以获得更好的效率以及稳定性。

1.2 非线性最小二乘问题

       假设 $\mathbf{f}(\mathbf{x})$ 为非线性函数，那么式（1.1）称为非线性最小二乘问题，这时通常无法或者很难直接得到式（1.5）的解析解，因此需要通过迭代的方法查找，即给定一组初始的参数值 ${\mathbf{x}}_0$，然后在此基础上查找下一组参数，依此类推，直至满足相关的收敛条件或者到达最大的迭代次数。查找的策略有很多种，这些方法在易用性、鲁棒性和收敛速度等方面各有千秋，所以彼此并不是割裂的，采用混合的方法或许能够获取更好的性能。注意，非线性最小二乘问题只在乎极值点的查找，其并不保证能得到全局最优解，得到的极值点通常与初始位置相关，但也并不一定是最靠近初始位置的那个极值点。不失一般性，后面主要讨论极小值点的搜索。假设该函数为全局凸函数，那么可以保证找到的极值点是全局最小值点。

1.2.1 线性搜索

       线性搜索（Line Search）属于最简单的搜索方法，即给定一个初始点 $\mathbf{x}$ 和一个搜索方向 $\mathbf{h}$，找到最佳的步长，使得更新后的位置的函数值最小，即

$$\begin{array}{c}\underset{\alpha }{\min }\phi \left(\alpha \right)=F\left(\mathbf{x}+\alpha \mathbf{h}\right)\\ s.t.F\left(\mathbf{x}+\alpha \mathbf{h}\right)⩽F\left(\mathbf{x}\right)\end{array}\text{\hspace{1em}(1.7)}$$

可以看到，线性搜索的性能好坏与给定的搜索方向密切相关，所以通常不会独立使用，而是配合其他一些优化算法，确定一个适合的方向再进行搜索。对于式（1.7），令其导数为零，有

$${\phi }^{\prime }\left({\alpha }^{\ast }\right)={\mathbf{h}}^T{\mathbf{F}}^{\prime }\left(\mathbf{x}+{\alpha }^{\ast }\mathbf{h}\right)=0.\text{\hspace{1em}(1.8)}$$

这是一个重要的性质，即如果式（1.8）有解，那么意味着最佳的更新位置的梯度方向与原来的搜索方向是垂直的。

1.2.2 最速下降法

       由于某一点的梯度方向代表函数在该点上升最快的方向，因此我们有理由认为沿着其梯度的反方向搜索能够更容易找到极小值点。具体来说，给定任意一个非零搜索方向 $\mathbf{h}$，根据泰勒公式展开可得

$$F(\mathbf{x}+\alpha \mathbf{h})=F(\mathbf{x})+\alpha {\mathbf{h}}^T{\mathbf{F}}^{\prime }(\mathbf{x})+o(\alpha )\approx F(\mathbf{x})+\alpha {\mathbf{h}}^T{\mathbf{F}}^{\prime }(\mathbf{x})\text{\hspace{1em}(1.9)}$$

当 $\alpha$ 趋向于无穷小时，余项 $o(\alpha )$ 可以忽略，那么其函数值下降速度为

$$\underset{\alpha \to 0}{\lim }\frac{F(\mathbf{x})-F(\mathbf{x}+\alpha \mathbf{h})}{\alpha \Vert \mathbf{h}\Vert }=-\frac{1}{\Vert \mathbf{h}\Vert }{\mathbf{h}}^T{\mathbf{F}}^{\prime }(\mathbf{x})=-\Vert {\mathbf{F}}^{\prime }(\mathbf{x})\Vert \cos \theta .\text{\hspace{1em}(1.10)}$$

其中 $\theta$ 为搜索方向 $\mathbf{h}$ 和梯度方向的夹角。可以看到，当 $\theta$ 为 $\pi$ 即为梯度的反方向时，函数值下降的速度最快。

       因为每次的搜索方向都为该点梯度的反方向，所以这种搜索方法被称为梯度下降法（Gradient Descent）或者最速下降法（Steepest Descent）。这种方法通常配合线性搜索以获取最佳的步长。由于最速下降法只用到了一阶导数，所以比较简单，易于实现。同时，因为线性搜索可以沿着梯度方向一直前进，所以实际的步长可以比较大，特别是当搜索方向上函数梯度方向变化较小时，这时即使初始点离极值点较远也有比较快的收敛速度，鲁棒性也比较强。因此，最速下降法的应用十分广泛，特别是在深度学习领域，这里不细说。

图1 最速下降法搜索路径示例

       然而，尽管被称为最速下降法，由于其只用到一阶导数，所以当搜索方向上函数梯度方向变化较快时，梯度下降法具有比较差的最终收敛性能，即靠近极值点时，其收敛速度有可能大幅减缓甚至停滞。这是在线性搜索中由式（1.8）所导致的后果，因为式（1.8）所求得的最佳步长要求当前搜索方向与下一个点的梯度方向垂直，而该点的梯度方向的反方向正是下一次搜索的方向，所以在最速下降法中，相邻两次搜索方向是相互垂直的，导致其搜索路径是锯齿形状的，在接近极值点时很难到达精确的位置。

       图 1 展示了两个简单的例子，其中 $n=2$，黑色的圈代表函数的等高线，越往内部其值越小，中心处为极小值点，可以证明梯度的方向与等高线的切线垂直。当函数的等高线为一系列的同心圆时，可以知道在同一条半径上的点的梯度方向是相同的，因此，对于任意一个初始点，沿着其梯度的反方向进行线性搜索，总可以一次性地到达圆心处，这时候其收敛速度是最快的。当函数的等高线为一系列的同心且形状相同的椭圆时，如果初始点刚好落在长轴或者短轴上，那么也可以一次性收敛到极小值点。然而，当初始点在其他位置时，其最终的收敛性能将变得很差。因为由式（1.7）所定义的线性搜索是要找到当前搜索方向上函数值最小的点，因此其下一个点必定为当前搜索方向即当前梯度的反方向与某个椭圆等高线的切点。由于椭圆只有跟长轴和短轴相交的4个顶点的切线的垂直线经过中心即函数的极值点，所以搜索路径只会成锯齿状地无限靠近而无法精确地到达极值点。因此，某些梯度下降算法可能会使用自定义的步长而不是通过线性搜索来寻找最佳步长，这样就可以避免相邻两次的搜索方向相互垂直，但这又有可能导致搜索路径在接近极值点时来回震荡。

1.2.3 牛顿迭代法

       从图 1 右侧我们已经可以看到，在接近极值点时，单纯在梯度反方向上搜索并不一定能获得最好的性能，这也是只使用一阶导数的缺点。那么一个自然的想法就是，在接近极值点时，不再只限定梯度这一个方向，而是在一个较小的邻域上面搜索，并且希望使用更高阶的导数信息。这就是牛顿迭代法（Newton’s method）以及后面一些改进方法的基本思想。由于三阶以上的导数求解比较困难，一般只使用到二阶的导数。

       在当前点的较小邻域内，假设有一个增量 $\mathbf{h}$，根据泰勒公式，可得

$$\begin{array}{c}F(\mathbf{x}+\mathbf{h})=F(\mathbf{x})+{\mathbf{h}}^T{\mathbf{F}}^{\prime }(\mathbf{x})+\frac{1}{2}{\mathbf{h}}^T{\mathbf{F}}^{\prime \prime }(\mathbf{x})\mathbf{h}+o\left({\Vert \mathbf{h}\Vert }^2\right),\\ F(\mathbf{x}+\mathbf{h})\approx L_{\text{n}}(\mathbf{h})=F(\mathbf{x})+{\mathbf{h}}^T{\mathbf{F}}^{\prime }(\mathbf{x})+\frac{1}{2}{\mathbf{h}}^T{\mathbf{F}}^{\prime \prime }(\mathbf{x})\mathbf{h}.\end{array}\text{\hspace{1em}(1.11)}$$

$\mathbf{h}$ 较小时，泰勒公式余项部分可忽略。其实质是用一个经过当前点的二次曲面 $L_{\text{n}}(\mathbf{h})$ 来拟合 $F(\mathbf{x})$ 。那么，我们希望找到一个最佳方向和增量 ${\mathbf{h}}_{\text{n}}$，使得 $L_{\text{n}}(\mathbf{h})$ 最小，有

$$L_{\text{n}}^{\prime }({\mathbf{h}}_{\text{n}})={\mathbf{F}}^{\prime }(\mathbf{x})+\frac{1}{2}\left({\mathbf{F}}^{\prime \prime }(\mathbf{x})+{\mathbf{F}}^{\prime \prime }{(\mathbf{x})}^T\right){\mathbf{h}}_{\text{n}}={\mathbf{F}}^{\prime }(\mathbf{x})+{\mathbf{F}}^{\prime \prime }(\mathbf{x}){\mathbf{h}}_{\text{n}}=0.\text{\hspace{1em}(1.12)}$$

其中

$$\begin{array}{c}{\left[{\mathbf{F}}^{\prime \prime }(\mathbf{x})\right]}_{jk}={\left[{\mathbf{F}}^{\prime \prime }(\mathbf{x})\right]}_{kj}=\sum _{i=1}^m\left[\frac{\partial f_i(\mathbf{x})}{\partial x_j}\frac{\partial f_i(\mathbf{x})}{\partial x_k}+f_i(\mathbf{x})\frac{\partial {}^2{f}_i(\mathbf{x})}{\partial x_j\partial x_k}\right],\\ {\mathbf{H}}_F(\mathbf{x})={\mathbf{F}}^{\prime \prime }(\mathbf{x})={\mathbf{J}}_f(\mathbf{x}{)}^T{\mathbf{J}}_f(\mathbf{x})+\sum _{i=1}^m{f}_i(\mathbf{x}){{\mathbf{f}}^{\prime \prime }}_i(\mathbf{x})\in {\mathbb{R}}^{n\times n}.\end{array}\text{\hspace{1em}(1.13)}$$

${\mathbf{H}}_F(\mathbf{x})$ 称为 $F(\mathbf{x})$ 的 Hessian 矩阵，其是一个实对称矩阵。为了简洁，在不影响理解的情况下，后面有时会省略部分符号的参数 $(\mathbf{x})$，即如 ${\mathbf{f},\mathbf{J}}_f,{\mathbf{H}}_F$ 等。那么，如果 ${\mathbf{H}}_F$ 的逆矩阵存在，式（1.12）的解为

$${\mathbf{h}}_{\text{n}}=-{\mathbf{F}}^{\prime \prime }(\mathbf{x}{)}^{-1}{\mathbf{F}}^{\prime }(\mathbf{x})=-{\mathbf{H}}_F^{-1}{\mathbf{J}}_f^T\mathbf{f}.\text{\hspace{1em}(1.14)}$$

若 ${\mathbf{H}}_F$ 为正定矩阵，则可以保证其逆矩阵存在，且 ${\mathbf{h}}_{\text{n}}$ 是 $L_{\text{n}}(\mathbf{h})$ 的最小值点，如式（1.15）所示。同理，若 ${\mathbf{H}}_F$ 为负定矩阵，${\mathbf{h}}_{\text{n}}$ 是 $L_{\text{n}}(\mathbf{h})$ 的最大值点。

$$\begin{array}{rl}{L}_{\text{n}}({\mathbf{h}}_{\text{n}})& =F(\mathbf{x})-{\mathbf{h}}_{\text{n}}^T{\mathbf{F}}^{\prime \prime }(\mathbf{x}){\mathbf{h}}_{\text{n}}+\frac{1}{2}{\mathbf{h}}_{\text{n}}^T{\mathbf{F}}^{\prime \prime }(\mathbf{x}){\mathbf{h}}_{\text{n}}\\ & =F(\mathbf{x})-\frac{1}{2}{\mathbf{h}}_{\text{n}}^T{\mathbf{F}}^{\prime \prime }(\mathbf{x}){\mathbf{h}}_{\text{n}}<F(\mathbf{x})=L_{\text{n}}(0).\end{array}\text{\hspace{1em}(1.15)}$$

       可以看到，相比于最速下降法只在梯度反方向上搜索，牛顿法是在一个邻域内搜索，而且其最优搜索方向实际上相当于梯度方向关于 ${\mathbf{H}}_F$ 的一个空间变换，因此具有很好的局部收敛性能。然而，牛顿法的局限性也很明显，使用二阶导数是其优点，同时也是缺点，因为计算 ${\mathbf{H}}_F$ 矩阵的难度和计算量都比较大，特别是有些时候我们甚至不知道 $\mathbf{f}(\mathbf{x})$ 的具体形式。而且，${\mathbf{H}}_F$ 并不总是正定或者可逆的。同时，因为其使用的是二次曲面近似，虽然对于二次型函数能够精确表示并一次性收敛，然而对于大多数具有更高阶导数的函数，其泰勒公式余项有时并不能忽略，有可能导致求得的 ${\mathbf{h}}_{\text{n}}$ 超出了所能近似的邻域范围。根据式（1.15），如果二阶项大于泰勒公式余项，那么依然可以保证迭代是往下降方向前进的，但在离极值点较远时所需要的迭代的次数较多；而如果余项部分太大，则反而会导致迭代后的函数值上升。因此，牛顿法的收敛性能和初始点的选择有很大的关系。

       为了解决牛顿法的稳定性问题，不少基于混合的方法被提出，例如我们可以结合最速下降法和牛顿法两者的优点，即当 ${\mathbf{H}}_F$ 正定时，优先采用牛顿法，否则使用最速下降法。通常在远离极值点时，最速下降法往往会更有优势；而当搜索进入后期阶段即靠近极值点时，牛顿法的二次曲面近似效果比较好，这时就能更好地发挥作用。或者，类似于最速下降法，我们可以结合牛顿法和线性搜索，即在 ${\mathbf{h}}_{\text{n}}$ 前面再乘以一个步长，并通过线性搜索来寻找最佳的步长，这样就可以避免迭代出现退化，并且，因为 ${\mathbf{h}}_{\text{n}}$ 不是梯度的方向，因此相邻两次的搜索方向不会相互垂直，在一定程度上避免了最速下降法的收敛困难问题。

1.2.4 高斯-牛顿法

       从式（1.13）可以看出，${\mathbf{H}}_F$ 需要计算各 $f_i(\mathbf{x})$ 的二阶偏导，然而这并不是一件易事，特别是有时候我们并不知道 $f_i(\mathbf{x})$ 的具体形式，通过有限差分法近似一来计算量大，二来计算机精度不足容易造成误差。高斯-牛顿法（Gauss-Newton method）就是要解决这样的问题。

       如果在一个较小的邻域内，各 $f_i(\mathbf{x})$ 近似为线性函数，那么根据泰勒公式，

$$\mathbf{f}(\mathbf{x}+\mathbf{h})\approx \mathbf{f}(\mathbf{x})+{\mathbf{J}}_f(\mathbf{x})\mathbf{h}.\text{\hspace{1em}(1.16)}$$

$$F(\mathbf{x}+\mathbf{h})=\frac{1}{2}\mathbf{f}(\mathbf{x}+\mathbf{h}{)}^T\mathbf{f}(\mathbf{x}+\mathbf{h})\approx L_{\text{gn}}(\mathbf{h})=\frac{1}{2}{\mathbf{f}}^T\mathbf{f}+{\mathbf{f}}^T{\mathbf{J}}_f\mathbf{h}+\frac{1}{2}{\mathbf{h}}^T{\mathbf{J}}_f^T{\mathbf{J}}_f\mathbf{h}.\text{\hspace{1em}(1.17)}$$

这时

$$L_{\text{gn}}^{\prime }(\mathbf{h})={\mathbf{J}}_f^T\mathbf{f}+{\mathbf{J}}_f^T{\mathbf{J}}_f\mathbf{h}.\text{\hspace{1em}(1.18)}$$

可以看到，相比于式（1.12）和（1.13），由于各 $f_i(\mathbf{x})$ 被近似为线性函数，所以其二阶偏导部分为零，整个计算过程就只涉及 $\mathbf{f}(\mathbf{x})$ 的雅可比矩阵即一阶导数了。然而，根据式（1.17），$F(\mathbf{x}+\mathbf{h})$ 仍然是二阶的近似，因此虽然高斯-牛顿法不用计算 ${\mathbf{H}}_{\text{F}}$，但在接近极值点时依然具有与牛顿法相似的收敛性能。

       根据式（1.18），可以得到最佳的搜索方向和步长为

$$L_{\text{gn}}^{\prime }({\mathbf{h}}_{\text{gn}})={\mathbf{J}}_f^T\mathbf{f}+{\mathbf{J}}_f^T{\mathbf{J}}_f{\mathbf{h}}_{\text{gn}}=0\Rightarrow {\mathbf{h}}_{\text{gn}}=-{\left({\mathbf{J}}_f^T{\mathbf{J}}_f\right)}^{-1}{\mathbf{J}}_f^T\mathbf{f}.\text{\hspace{1em}(1.19)}$$

这里同样要求 ${\mathbf{J}}_f^T{\mathbf{J}}_f$ 是正定的。相比于 ${\mathbf{H}}_{\text{F}}$ 有可能是正定、负定、不定等，${\mathbf{J}}_f^T{\mathbf{J}}_f$则可以确定是半正定的，即对于任意一个非零向量 $\mathbf{h}$，有

$$\forall \mathbf{h}\ne 0,{\mathbf{h}}^T{\mathbf{J}}_f^T{\mathbf{J}}_f\mathbf{h}={\left({\mathbf{J}}_f\mathbf{h}\right)}^T\left({\mathbf{J}}_f\mathbf{h}\right)⩾0.\text{\hspace{1em}(1.20)}$$

当 ${\mathbf{J}}_f^T{\mathbf{J}}_f$ 满秩时，则可以证明是正定的，那么式(1.19)成立，同时可以保证

$$\begin{array}{rl}{L}_{\text{gn}}({\mathbf{h}}_{\text{gn}})& =\frac{1}{2}{\mathbf{f}}^T\mathbf{f}+{\mathbf{f}}^T{\mathbf{J}}_f{\mathbf{h}}_{\text{gn}}+\frac{1}{2}{\mathbf{h}}_{\text{gn}}^T{\mathbf{J}}_f^T{\mathbf{J}}_f{\mathbf{h}}_{\text{gn}}\\ & =F(\mathbf{x})-\frac{1}{2}{\mathbf{h}}_{\text{gn}}^T{\mathbf{J}}_f^T{\mathbf{J}}_f{\mathbf{h}}_{\text{gn}}<F(\mathbf{x})=L_{\text{gn}}(0).\end{array}\text{\hspace{1em}(1.21)}$$

因为

R ( J f T J f ) = R ( J f ) ⩽ min ⁡ { m ,   n } . R\left( {{\mathbf{J}}_f^T{{\mathbf{J}}_f}} \right) = R\left( {{{\mathbf{J}}_f}} \right) \leqslant \min \left\{ {m,{\text{ }}n} \right\}. R(JfT​Jf​)=R(Jf​)⩽min{m, n}.

所以一般要求 $m⩾n$，即函数 $f_i(\mathbf{x})$ 的个数要比参数多，在数据拟合中意味着约束方程要比需要优化的参数个数多，并且这些约束方程彼此应该是独立的。因为这些条件一般是满足的，所以式（1.19）通常成立，相比于牛顿法就更有通用性。

       然而，必须注意的是，以上的分析都是基于式（1.16）的，即 $\mathbf{f}(\mathbf{x})$ 在较小邻域内近似为线性函数，而这在非线性最小二乘问题中很明显是没法保证成立的，这就要求邻域足够小，因此其同样不适合初始点离极值点较远的情况。另外，除了线性近似，根据式（1.13），当 ${{\mathbf{f}}^{\prime \prime }}_i(\mathbf{x})$ 变化较小，而 { f i ( x ) } \left\{ {{f_i}({\mathbf{x}})} \right\} {fi​(x)} 分布接近于白噪声时，那么我们可以认为其等号右侧第二项近似为零，即

$$\sum _{i=1}^m{f}_i(\mathbf{x}){{\mathbf{f}}^{\prime \prime }}_i(\mathbf{x})\approx \mathbf{A}\cdot \sum _{i=1}^m{f}_i(\mathbf{x})\approx \mathbf{O}.\text{\hspace{1em}(1.22)}$$

这也符合一些数据拟合问题的情况，这时我们则不需要对 $\mathbf{f}(\mathbf{x})$ 进行线性近似。

1.2.5 Levenberg–Marquardt法

       通过前面的分析我们发现，牛顿法和高斯-牛顿法都需要进行泰勒公式展开近似，这就要求邻域足够小，问题是应该多小？如果邻域太小，而初始点离极值点比较远，我们是不是应该考虑使用最速下降法？如果近似程度比较好我们是不是应该增大搜索邻域？对于这些问题，后来也衍生出一些称为基于信赖域的优化算法，例如 Powell 的狗腿法（Dog-leg method）等。Levenberg（1944）和Marquartdt（1963）两人则在高斯-牛顿法的基础上增加了一个阻尼系数提出了 LM 算法，其虽然不是明确地基于信赖域，但其实质则和信赖域法有异曲同工之妙。

       如式（1.23）所示，

$$L_{\text{lm}}(\mathbf{h})=\frac{1}{2}{\mathbf{f}}^T\mathbf{f}+{\mathbf{f}}^T{\mathbf{J}}_f\mathbf{h}+\frac{1}{2}{\mathbf{h}}^T{\mathbf{J}}_f^T{\mathbf{J}}_f\mathbf{h}+\frac{1}{2}\mu {\mathbf{h}}^T\mathbf{h},\mu >0.\text{\hspace{1em}(1.23)}$$

相比式（1.17），LM 法增加了一项惩罚项，或者正则化项，其中 $\mu$ 称为阻尼系数。惩罚项的意义十分明显，因为高斯-牛顿法需要对 $\mathbf{f}(\mathbf{x})$ 进行线性近似，所以我们必须要控制邻域的大小。当近似效果较好时，我们能在较大的邻域内进行搜索，这时 $\mu$ 应该是一个比较小的值，这样惩罚项在 $L_{\text{lm}}(\mathbf{h})$ 中所占比重较小，$\mathbf{h}$ 则可以相对较大；而当近似效果较差时，我们应该把邻域收紧到一个较小的范围，这时 $\mu$ 应该是一个比较大的值，这样惩罚项在 $L_{\text{lm}}(\mathbf{h})$ 中所占比重较大，$\mathbf{h}$ 则只能较小。

       对 $L_{\text{lm}}(\mathbf{h})$ 进行求导，可得

$$L_{\text{lm}}^{\prime }(\mathbf{h})={\mathbf{J}}_f^T\mathbf{f}+\left({\mathbf{J}}_f^T{\mathbf{J}}_f+\mu \mathbf{I}\right)\mathbf{h}={\mathbf{J}}_f^T\mathbf{f}+\mathbf{H}\mathbf{h}.\text{\hspace{1em}(1.24)}$$

注意这时

$$\forall \mathbf{h}\ne 0,{\mathbf{h}}^T\mathbf{H}\mathbf{h}={\mathbf{h}}^T\left({\mathbf{J}}_f^T{\mathbf{J}}_f+\mu \mathbf{I}\right)\mathbf{h}={\left({\mathbf{J}}_f\mathbf{h}\right)}^T\left({\mathbf{J}}_f\mathbf{h}\right)+\mu {\mathbf{h}}^T\mathbf{h}>0.\text{\hspace{1em}(1.25)}$$

这意味着 $\mathbf{H}$ 是严格正定的，因此其逆矩阵也总是存在，那么

$$L_{\text{lm}}^{\prime }({\mathbf{h}}_{\text{lm}})={\mathbf{J}}_f^T\mathbf{f}+\left({\mathbf{J}}_f^T{\mathbf{J}}_f+\mu \mathbf{I}\right){\mathbf{h}}_{\text{lm}}=0\Rightarrow {\mathbf{h}}_{\text{lm}}=-{\left({\mathbf{J}}_f^T{\mathbf{J}}_f+\mu \mathbf{I}\right)}^{-1}{\mathbf{J}}_f^T\mathbf{f}.\text{\hspace{1em}(1.26)}$$

这时，

$$\begin{array}{c}\mu \to 0,{\mathbf{h}}_{\text{lm}}\approx -{\left({\mathbf{J}}_f^T{\mathbf{J}}_f\right)}^{-1}{\mathbf{J}}_f^T\mathbf{f},\\ \mu \to \infty ,{\mathbf{h}}_{\text{lm}}\approx -\frac{1}{\mu }{\mathbf{J}}_f^T\mathbf{f}=-\frac{1}{\mu }{\mathbf{F}}^{\prime }(\mathbf{x}).\end{array}\text{\hspace{1em}(1.27)}$$

那么，当 $\mu$ 趋近于 0 时，${\mathbf{h}}_{\text{lm}}$ 趋近于高斯-牛顿法的解，说明这时候对 $\mathbf{f}(\mathbf{x})$ 的线性近似效果比较好；当 $\mu$ 趋近于无穷大时，${\mathbf{h}}_{\text{lm}}$ 趋近于一个小步长的最速下降法，说明这时对 $\mathbf{f}(\mathbf{x})$ 的线性近似效果较差，为了保证算法的鲁棒性，我们更加趋向于在梯度的反方向上搜索。然而，怎么界定 $\mu$ 是大是小？同时，$\mathbf{f}(\mathbf{x})$ 在不同位置的线性近似程度肯定不同，我们应该以什么标准来根据优化的情况更新 $\mu$ 值？

       首先我们得确定一个适当的初始值 ${\mu }_0$。如果我们可以确定初始点 ${\mathbf{x}}_0$ 是在极值点附近，那么 ${\mu }_0$ 应该相对比较小，那么搜索方式就更加偏向于高斯-牛顿法；如果初始点 ${\mathbf{x}}_0$ 离极值点可能比较远，那么那么 ${\mu }_0$ 应该大一些，那么搜索方式就更加偏向于最速下降法。根据式（1.26），$\mu$ 的大小是与 ${\mathbf{J}}_f^T{\mathbf{J}}_f$ 的元素大小相关的，因为 ${\mathbf{J}}_f^T{\mathbf{J}}_f$ 是一个半正定实对称矩阵，因此可以证明其特征值 { λ i } \{\lambda_i\} {λi​} 是实数且非负，且不同特征值对应的特征向量必定正交，重根特征值对应的特征向量也总可以经过正交化得到相同数量的相互正交的特征向量。所以，我们可将这些相互正交的特征向量 { v i } \{{\bf{v}}_i\} {vi​} 作为 ${\mathbb{R}}^n$ 空间的一个规范正交基。那么，${\mathbf{h}}_{\text{lm}}$ 可表示为 { v i } \{{\bf{v}}_i\} {vi​} 的线性组合，即有

$${\mathbf{h}}_{\text{lm}}=\sum _{i=1}^n{a}_i{\mathbf{v}}_i=\sum _{i=1}^n{\mathbf{v}}_i^T{\mathbf{h}}_{\text{lm}}{\mathbf{v}}_i=\sum _{i=1}^n{\mathbf{v}}_i{\mathbf{v}}_i^T{\mathbf{h}}_{\text{lm}}.\text{\hspace{1em}(1.28)}$$

根据式（1.26），有

$$\begin{array}{c}\left({\mathbf{J}}_f^T{\mathbf{J}}_f+\mu \mathbf{I}\right){\mathbf{h}}_{\text{lm}}=\sum _{i=1}^n\left({\mathbf{J}}_f^T{\mathbf{J}}_f+\mu \mathbf{I}\right){\mathbf{v}}_i{\mathbf{v}}_i^T{\mathbf{h}}_{\text{lm}}\\ =\sum _{i=1}^n\left({\lambda }_i+\mu \right){\mathbf{v}}_i{\mathbf{v}}_i^T{\mathbf{h}}_{\text{lm}}=-{\mathbf{F}}^{\prime }(\mathbf{x}).\end{array}\text{\hspace{1em}(1.29)}$$

那么，

$$\begin{array}{c}{\mathbf{v}}_j^T\left[\sum _{i=1}^n\left({\lambda }_i+\mu \right){\mathbf{v}}_i{\mathbf{v}}_i^T{\mathbf{h}}_{\text{lm}}\right]{\mathbf{v}}_j=\left({\lambda }_j+\mu \right){\mathbf{v}}_j{\mathbf{v}}_j^T{\mathbf{h}}_{\text{lm}}=-{\mathbf{v}}_j^T{\mathbf{F}}^{\prime }(\mathbf{x}){\mathbf{v}}_j,\\ {\mathbf{h}}_{\text{lm}}=\sum _{j=1}^n{\mathbf{v}}_j{\mathbf{v}}_j^T{\mathbf{h}}_{\text{lm}}=-\sum _{j=1}^n\frac{{\mathbf{v}}_j^T{\mathbf{F}}^{\prime }(\mathbf{x})}{{\lambda }_j+\mu }{\mathbf{v}}_j.\end{array}\text{\hspace{1em}(1.30)}$$

因此，$\mu$ 的大小和特征值 { λ i } \{λ_i\} {λi​} 的大小相关。根据瑞利商的性质，实对称矩阵的对角元素的值介于最小和最大特征值之间，因此我们可以定义初始 ${\mu }_0$ 为

$${\mu }_0=\tau \cdot \max \left\{{\left({\mathbf{J}}_f^T{\mathbf{J}}_f\right)}_{ii}\right\}.\text{\hspace{1em}(1.31)}$$

其中 $\tau$ 可根据前面的分析定义，例如倾向高斯-牛顿法则取较小的数如 $10^{-6}$，倾向最速下降法则取较大的数如 $10^{-3}$ 甚至 $1$。

       确定初始 ${\mu }_0$ 后，接下来是要确定 $\mu$ 的更新规则。根据前面的内容，我们知道 $\mu$ 的存在是为了保证 $\mathbf{f}(\mathbf{x})$ 在一定邻域内的线性近似程度，而式（1.17）是 $\mathbf{f}(\mathbf{x})$ 线性近似的结果，那么我们可以定义 $\mathbf{f}(\mathbf{x})$ 线性近似的程度为

$$\rho =\frac{F(\mathbf{x})-F(\mathbf{x}+{\mathbf{h}}_{\text{lm}})}{L_{\text{gn}}(0)-L_{\text{gn}}({\mathbf{h}}_{\text{lm}})}.\text{\hspace{1em}(1.32)}$$

其中

$$\begin{array}{rl}{L}_{\text{gn}}(0)-L_{\text{gn}}({\mathbf{h}}_{\text{lm}})& =-{\mathbf{f}}^T{\mathbf{J}}_f{\mathbf{h}}_{\text{lm}}-\frac{1}{2}{\mathbf{h}}_{\text{lm}}^T{\mathbf{J}}_f^T{\mathbf{J}}_f{\mathbf{h}}_{\text{lm}}\\ & =-\frac{1}{2}{\mathbf{h}}_{\text{lm}}^T\left(2{\mathbf{J}}_f^T\mathbf{f}+{\mathbf{J}}_f^T{\mathbf{J}}_f{\mathbf{h}}_{\text{lm}}\right)\\ & =-\frac{1}{2}{\mathbf{h}}_{\text{lm}}^T\left[-2\left({\mathbf{J}}_f^T{\mathbf{J}}_f+\mu \mathbf{I}\right){\mathbf{h}}_{\text{lm}}+{\mathbf{J}}_f^T{\mathbf{J}}_f{\mathbf{h}}_{\text{lm}}\right]\\ & =\frac{1}{2}{\mathbf{h}}_{\text{lm}}^T\left({\mathbf{J}}_f^T{\mathbf{J}}_f+2\mu \mathbf{I}\right){\mathbf{h}}_{\text{lm}}>0.\end{array}\text{\hspace{1em}(1.33)}$$

       因此，${\mathbf{h}}_{\text{lm}}$ 可以保证式（1.32）的分母一定是正数。如果 $\rho \le 0$，即 $F({\mathbf{x}+\mathbf{h}}_{\text{lm}})$ 不降反升，那我们可以十分肯定 $\mathbf{f}(\mathbf{x})$ 的线性近似是有问题的，或者说搜索的邻域太大，那么这时我们不应该更新 $\mathbf{x}$，同时应该增大 $\mu$ 值，以缩小搜索的邻域；当 $\rho$ 是一个比较大的正数，最好是接近 $1$，这时我们就可以认为在该邻域内 $\mathbf{f}(\mathbf{x})$ 的线性近似是比较恰当的，而且有理由认为更新 $\mathbf{x}$ 后下一次搜索的邻域可以再扩大一些，因此这时候可以适当减小 $\mu$ 值；当 $\rho >0$ 但是又比较小，虽然更新后的 $\mathbf{x}$ 可以使得函数值下降，但下降的幅度远小于我们的预期，这时我们应该适当缩小下一次搜索的邻域，因此这时候可以适当增大 $\mu$ 值。综上所述，一种经典的 LM 算法的优化过程可以总结如下：

$$\begin{array}{rl}& \text{if }\rho <{\rho }_1,\text{ then }{\mu }_{\text{new}}\text{ = }\beta \ast \mu \\ & \text{if }\rho >{\rho }_2,\text{ then }{\mu }_{\text{new}}\text{ = }\mu /\gamma \\ & \text{if }\rho >0,\text{ then }{\mathbf{x}}_{\text{new}}\text{ = }\mathbf{x}+{\mathbf{h}}_{\text{lm}}\end{array}\text{\hspace{1em}(1.34)}$$

其中 $0<{\rho }_1<{\rho }_2<1$，而且 $\beta ,\gamma >1$。实验表明 LM 算法对这些参数的选择并不是很敏感，比较常用的参数如下

$${\rho }_1=0.25,{\rho }_2=0.75,\beta =2,\gamma =3.$$

       在式（1.34）中，${\mu }_{\text{new}}/\mu$ 在经过阈值 ${\rho }_1$ 和 ${\rho }_2$ 的时候会发生跳变，可能会导致优化过程发生震荡，因此，一些论文里面也提出了一些变化更加平滑的策略，而且也不需要给出明确的阈值。图 2 展示了一种完整的 LM 算法流程，其 ${\mu }_{\text{new}}/\mu$ 关于 $\rho$ 的函数图像如图 3 所示。

图2 LM算法流程示例
图3 图2所示算法（实线）和式(1.34)（虚线）的更新策略

       总体来说，LM 算法综合了最速下降法和高斯-牛顿的的特点，对初始点的要求相对不那么严格，在收敛性能和鲁棒性方面都表现不错，因此在很多参数估计问题中都有应用，比如运动参数估计、相机姿态参数估计等等。总之，非线性最小二乘问题是一个相对复杂的问题，怎么在速度、易用性、鲁棒性等方面进行取舍需要根据实际的问题进行分析。除了前面提到的这几种算法，还有如前面提到的基于信赖域的算法、拟牛顿法等等，以及它们的混合算法和一些改进，在这里不一一细说，有兴趣的可以详细阅读下面的参考资料以及在网上进行搜索。

参考资料：

• H.B. Nielsen. Damping Parameter in Maquardt’s Method.
• K. Madsen, H.B. Nielsen, O. Tingleff. Methods for Non-Linear Least Squares Problems.