Levenberg–Marquardt（LM）详解

1、基础概念
- 1.1、信赖域法
- 1.2、泰勒展开
- 1.2、正定矩阵(positive definite matrix)
- 1.3、雅克比矩阵(Jacobian matrix)
- 1.4、黑塞矩阵(Hessian matrix)
- 1.5、范数(norm)
- 1.6、非线性最小二乘问题(Non-linear least squares problems)
2、LM算法
- 2.1、高斯-牛顿法（Gauss-Newton Method）
- 2.2、Levenberg-Marquardt(LM) 算法

1、基础概念

1.1、信赖域法

在最优化算法中，都是要求一个函数的极小值，每一步迭代中，都要求目标函数值是下降的，而信赖域法，顾名思义，就是从初始点开始，先假设一个可以信赖的最大位移 s，然后在以当前点为中心，以 s 为半径的区域内，通过寻找目标函数的一个近似函数（二次的）的最优点，来求解得到真正的位移。在得到了位移之后，再计算目标函数值，如果其使目标函数值的下降满足了一定条件，那么就说明这个位移是可靠的，则继续按此规则迭代计算下去；如果其不能使目标函数值的下降满足一定的条件，则应减小信赖域的范围，再重新求解。其数学模型如下所示：

{minmk(s)=fk+gkTs+12sTGkss.t.∣∣s∣∣≤hk\left\{\begin{matrix} min &m_k(s)=f_k+{g_k}^Ts +\frac{1}{2}s^TG_ks & \\ s.t. & ||s|| \leq h_k& \end{matrix}\right.{mins.t.mk(s)=fk+gkTs+21sTGks∣∣s∣∣≤hk

其中，第一个式子就是我们用于模拟目标函数的二次模型，其自变量为 s，也就是我们要求的位移。gkg_kgk 为梯度， GkG_kGk 为Hesse矩阵，袁亚湘的书上说，如果Hesse矩阵不好计算，可以利用“有限差分”来近似 GkG_kGk, 或者用拟牛顿方法来构造Hesse矩阵的近似矩阵。
第二个式子中的 hkh_khk 是第 k 次迭代的信赖域上界（或称为信赖域半径），因此第二个式子表示的就是位移要在信赖域上界范围内。
通过衡量二次模型与目标函数的近似程度，可以作出判定是否需要扩大信了：

第 k 次迭代的实际下降量为： Δfk=fk−f(xk+sk)Δf_k=f_k−f(x_k+s_k)Δfk=fk−f(xk+sk)
第 k 次迭代的预测下降量为： Δmk=fk−m(sk)Δm_k=f_k−m(s_k)Δmk=fk−m(sk)

定义比值： rk=ΔfkΔmkr_k=\frac{Δf_k}{Δm_k}rk=ΔmkΔfk

这个比值可以用于衡量二次模型与目标函数的近似程度，显然 r 值越接近1越好。

1.2、泰勒展开

泰勒公式是一个用函数在某点的信息描述其附近取值的公式。如果函数满足一定的条件，泰勒公式可以用函数在某一点的各阶导数值做系数构建一个多项式来近似表达这个函数。

f(x)=∑i=0nf(i)(x0)i！(x−x0)if(x)=\sum_{i=0}^{n}\frac{f^{(i)}(x_0)}{i！}(x-x_0)^if(x)=∑i=0ni！f(i)(x0)(x−x0)i

其中，x∈(x0−Δx,x0+Δx)x∈(x_0-\Delta x,x_0+\Delta x)x∈(x0−Δx,x0+Δx)。

1.2、正定矩阵(positive definite matrix)

（1）广义定义
设M是n阶方阵，如果对任何非零向量z，都有 zTMz>0z^TMz> 0zTMz>0，其中 zTz^TzT 表示z的转置，就称M为正定矩阵。
例如：B为n阶矩阵，E为单位矩阵，a为正实数。在a充分大时，aE+B为正定矩阵。（B必须为对称阵）。

（2）狭义定义
一个n阶的实对称矩阵M是正定的的条件是当且仅当对于所有的非零实系数向量z，都有 zTMz>0z^TMz> 0zTMz>0。其中zT表示z的转置。

（3）正定矩阵的性质

正定矩阵的行列式恒为正；
实对称矩阵A正定当且仅当A与单位矩阵合同；
若A是正定矩阵，则A的逆矩阵也是正定矩阵；
两个正定矩阵的和是正定矩阵；
正实数与正定矩阵的乘积是正定矩阵。

（4）正定矩阵的特征

对称阵A为正定的充分必要条件是：A的特征值全为正。
对称阵A为正定的充分必要条件是：A的各阶顺序主子式都为正。
任意阵A为正定的充分必要条件是：A合同于单位阵。

1.3、雅克比矩阵(Jacobian matrix)

在向量分析中，雅可比矩阵是函数的一阶偏导数以一定方式排列成的矩阵，其行列式称为雅可比行列式。在代数几何中，代数曲线的雅可比行列式表示雅可比簇：伴随该曲线的一个代数群，曲线可以嵌入其中。它们全部都以数学家卡尔·雅可比命名；
假设某函数从 RnR^nRn 映到 RmR^mRm，其雅可比矩阵是从 RnR^nRn 到 RmR^mRm的线性映射，其重要意义在于它表现了一个多变数向量函数的最佳线性逼近。因此，雅可比矩阵类似于单变数函数的导数。
假设 Rn→RmR^n\rightarrow R^mRn→Rm 是一个从n维欧氏空间映射到到m维欧氏空间的函数。这个函数由m个实函数组成：y1(x1,...,xn),...,ym(x1,...,xn)y_1(x_1,...,x_n),...,y_m(x_1,...,x_n)y1(x1,...,xn),...,ym(x1,...,xn)。这些函数的偏导数(如果存在)可以组成一个m行n列的矩阵，这个矩阵就是所谓的雅可比矩阵：
(∂y1∂x1...∂y1∂xn...∂ym∂x1...∂ym∂xn)\begin{pmatrix} \frac{\partial y_1}{\partial x_1} & ... & \frac{\partial y_1}{\partial x_n}\\ . & . & .\\ \frac{\partial y_m}{\partial x_1} & ... & \frac{\partial y_m}{\partial x_n}\end{pmatrix}⎝⎛∂x1∂y1.∂x1∂ym.......∂xn∂y1.∂xn∂ym⎠⎞。

记作，JF(x1,...,xn)J_F(x_1,...,x_n)JF(x1,...,xn) 或 ∂(y1,...,ym)∂(x1,...,xn)\frac{\partial (y_1,...,y_m)}{\partial (x_1,...,x_n)}∂(x1,...,xn)∂(y1,...,ym)。
如果p是中的一点，F在 p点可微分，根据高等微积分，JF(p)J_F(p)JF(p) 是在这点的导数。在此情况下，这个线性映射即F在点p附近的最优线性逼近，也就是说当x足够靠近点p时，我们有
F(x)≈F(p)+JF(p)⋅(x−p)F(x)\approx F(p)+J_F(p)\cdot (x-p)F(x)≈F(p)+JF(p)⋅(x−p)

1.4、黑塞矩阵(Hessian matrix)

1.5、范数(norm)

（1）lp−范数l_p-范数lp−范数的定义

（2）l0−范数l_0-范数l0−范数的定义

（3）l1−范数l_1-范数l1−范数的定义

（4）l2−范数l_2-范数l2−范数的定义

1.6、非线性最小二乘问题(Non-linear least squares problems)

设误差函数为 f(x)f(x)f(x)，则损失函数 F(x)F(x)F(x) 为：

F(x)=12∑i=1m(f(xi))2=12∣∣f(x)∣∣2=12fT(x)f(x)F(x)=\frac{1}{2}\sum_{i=1}^{m}(f(x_i))^2=\frac{1}{2}||f(x)||^2=\frac{1}{2}f^T(x)f(x)F(x)=21∑i=1m(f(xi))2=21∣∣f(x)∣∣2=21fT(x)f(x)

目标为对于向量函数 f(x):Rn−−>Rmf(x): R_n --> R_mf(x):Rn−−>Rm 其中 m>=n . 我们希望最小化 ∣∣f(x)∣∣||f(x)||∣∣f(x)∣∣，即

x∗=argminxF(x)x^*=argmin_x{F(x)}x∗=argminxF(x)。

首先对于 f(x)f(x)f(x) 进行泰勒展开，有

f(x+h)=f(x)+J(x)∗h+o(∣∣h∣∣2)f(x+h)=f(x)+J(x)*h+o(||h||^2)f(x+h)=f(x)+J(x)∗h+o(∣∣h∣∣2)

其中， (J(x))ij=∂fi∂xj(x)(J(x))_{ij}=\frac{\partial f_i}{\partial x_j}(x)(J(x))ij=∂xj∂fi(x) ，J∈Rm×nJ∈R^{m×n}J∈Rm×n 是雅克比矩阵(Jacobian matrix)，它包含函数 f(x)f(x)f(x) 分量一阶导数的矩阵。
又

∂F∂xj(x)=∑i=1mfi(x)∂fi∂xj(x)\frac{\partial F}{\partial x_j}(x)=\sum_{i=1}^{m}f_i(x)\frac{\partial f_i}{\partial x_j}(x)∂xj∂F(x)=∑i=1mfi(x)∂xj∂fi(x)

故得到 F(x)F(x)F(x) 的一阶偏导为

F′(x)=J(x)Tf(x){F}'(x)=J(x)^Tf(x)F′(x)=J(x)Tf(x)

同理，由对于 F(x)F(x)F(x) 的黑塞矩阵(Hessian matrix) 中第j行k列的元素有，

∂2F∂xj∂xk(x)=∑i=1m(∂fi∂xj(x)∂fi∂xk(x)+fi(x)∂2fi∂xj∂xk(x))\frac{\partial^2 F}{\partial x_j\partial x_k}(x)=\sum_{i=1}^{m}(\frac{\partial f_i}{\partial x_j}(x)\frac{\partial f_i}{\partial x_k}(x)+f_i(x)\frac{\partial ^2f_i}{\partial x_j\partial x_k}(x))∂xj∂xk∂2F(x)=∑i=1m(∂xj∂fi(x)∂xk∂fi(x)+fi(x)∂xj∂xk∂2fi(x))

得到 F(x)F(x)F(x) 的二阶偏导为

F′′(x)=J(x)TJ(x)+∑i=1mfi(x)fi′′(x){F}''(x)=J(x)^TJ(x)+\sum_{i=1}^{m}f_i(x){f}''_i(x)F′′(x)=J(x)TJ(x)+∑i=1mfi(x)fi′′(x)

2、LM算法

2.1、高斯-牛顿法（Gauss-Newton Method）

最快梯度下降法方法忽略了二阶导数项，最终阶段为线性收敛，且速度较慢，因此，更多时候是将其作为优化初始阶段所采用方法。而牛顿法使用了二阶导数，最终阶段收敛速度快，收敛性好，但不适合初始阶段。因此常将其与最快梯度下降法结合使用。由前面的内容，我们已经知道对于 f(x)f(x)f(x) 进行泰勒展开有

f(x+h)≃l(h)≡f(x)+J(x)hf(x+h)\simeq l(h)\equiv f(x)+J(x)hf(x+h)≃l(h)≡f(x)+J(x)h

则有

F(x+h)≃L(h)≡12f(x+h)Tf(x+h)=12l(h)Tl(h)F(x+h)\simeq L(h)\equiv \frac{1}{2}f(x+h)^Tf(x+h)=\frac{1}{2}l(h)^Tl(h)F(x+h)≃L(h)≡21f(x+h)Tf(x+h)=21l(h)Tl(h)

即

F(x+h)≡12f(x)Tf(x)+hTJ(x)Tf(x)+12hTJ(x)TJ(x)h=F(x)+hTJ(x)Tf(x)+12hTJ(x)TJ(x)hF(x+h)\equiv \frac{1}{2}f(x)^Tf(x)+h^TJ(x)^Tf(x)+\frac{1}{2}h^TJ(x)^TJ(x)h=F(x)+h^TJ(x)^Tf(x)+\frac{1}{2}h^TJ(x)^TJ(x)hF(x+h)≡21f(x)Tf(x)+hTJ(x)Tf(x)+21hTJ(x)TJ(x)h=F(x)+hTJ(x)Tf(x)+21hTJ(x)TJ(x)h

故有

L′(h)=J(x)Tf(x)+J(x)TJ(x)h{L}'(h)=J(x)^Tf(x)+J(x)^TJ(x)hL′(h)=J(x)Tf(x)+J(x)TJ(x)h

L′′(h)=J(x)TJ(x){L}''(h)=J(x)^TJ(x)L′′(h)=J(x)TJ(x)

当 h=0h =0h=0 时有，

L′(0)=J(x)Tf(x)=F′(x){L}'(0)=J(x)^Tf(x)={F}'(x)L′(0)=J(x)Tf(x)=F′(x)

当J(x)J(x)J(x)满秩时，有 L′(hgn)=0L'(h_{gn})=0L′(hgn)=0 时有 J(x)TJ(x)hgn=−J(x)Tf(x)J(x)^TJ(x)h_{gn}=-J(x)^Tf(x)J(x)TJ(x)hgn=−J(x)Tf(x)，则对于 F(x)F(x)F(x) 有，

hgnTF′(x)=hgnT(JTf(x))=−hgnT(J(x)TJ(x)hgn)＜0h_{gn}^TF'(x) =h_{gn}^T(J^Tf(x))=-h_{gn}^T(J(x)^TJ(x)h_{gn})＜0hgnTF′(x)=hgnT(JTf(x))=−hgnT(J(x)TJ(x)hgn)＜0

故，我们可以得到梯度下降方向 hd=hgnh_d=h_{gn}hd=hgn ，迭代 x:=x+αhgnx:=x+\alpha h_{gn}x:=x+αhgn。
其中，步径 α\alphaα 通过直线搜索得到，对于经典的 Gauss-Newton method 采用 α=1\alpha=1α=1 。直线搜索的方法保证收敛性，需满足以下两个条件:

{x∣F(x)≤F(x0))}\begin{Bmatrix} x|F(x)\leq F(x_0)) \end{Bmatrix}{x∣F(x)≤F(x0))} 是有界的；
Jacobian J(x) 在所有迭代步骤中都是满秩的。
Newton’ method 和 Gauss-Newton method 方法的搜索方向分别为：
F′′(x)hn=−F′(x)F''(x)h_{n}=-F'(x)F′′(x)hn=−F′(x)
L′′(h)hgn=−L′(0)L''(h)h_{gn}=-L'(0)L′′(h)hgn=−L′(0)
其中，
F′′(x)=L′′(h)+∑i=1mfi(x)fi′′(x)F''(x)=L''(h)+\sum_{i=1}^{m}f_i(x)f''_i(x)F′′(x)=L′′(h)+∑i=1mfi(x)fi′′(x)

故，当 f(x∗)=0f(x^*)=0f(x∗)=0 时，对离最优解 x∗x^*x∗ 很近的 x 来说 L′′(h)约等于F′′(x)L''(h) 约等于 F''(x)L′′(h)约等于F′′(x)。但是牛顿法需要时刻计算 HHH 矩阵，即二阶导数信息，是计算量很大的一件事情，而LM算法则提出用雅可比矩阵（易计算）近似代替H矩阵的计算，使得优化效率得到提升。

2.2、Levenberg-Marquardt(LM) 算法

LM算法关键是用模型函数 fm对待估参数向量 p 在其邻域内做线性近似，忽略掉二阶以上的导数项，从而转化为线性最小二乘问题，它具有收敛速度快等优点。LM算法的通俗描述为，“如果目标函数值增大，则调整某系数再继续求解；如果目标函数值减小，则调整某系数再继续求解”的迭代过程，这种过程与上面所说的信赖域法是非常相似的，所以LM算法属于一种“信赖域法”。但LM算法需要对每一个待估参数求偏导，所以，如果你的目标函数 f 非常复杂，或者待估参数相当地多，那么可能不适合使用LM算法，而可以选择Powell算法——Powell算法不需要求导。LM算法的数学模型为：

{minmk(s)=fk+gkTs+12sTGkss.t.∣∣s∣∣2≤hk\left\{\begin{matrix} min &m_k(s)=f_k+{g_k}^Ts +\frac{1}{2}s^TG_ks & \\ s.t. & ||s||_2 \leq h_k& \end{matrix}\right.{mins.t.mk(s)=fk+gkTs+21sTGks∣∣s∣∣2≤hk

LM算法是介于牛顿法与梯度下降法之间的一种非线性优化方法，对于过参数化问题不敏感，能有效处理冗余参数问题，使代价函数陷入局部极小值的机会大大减小，这些特性使得LM算法在计算机视觉等领域得到广泛应用。LM算法中用雅可比矩阵近似代替H矩阵的计算来提升效率。其对应的搜索方向为

(J(x)TJ(x)+μI)hlm=−g=−J(x)Tf(x)(J(x)^TJ(x)+\mu I)h_{lm}=-g=-J(x)^Tf(x)(J(x)TJ(x)+μI)hlm=−g=−J(x)Tf(x)

其中， μ\muμ 需大于等于0，对应不同的取值范围下亦有不同的效果：

当 μ>0\mu>0μ>0 ，参数矩阵是正定，这保证了 hlmh_{lm}hlm 是一个下降方向。
当 μ\muμ 较大时，可以得到

hlm≃−1μg=−1μF′(x)h_{lm}\simeq -\frac{1}{\mu}g=-\frac{1}{\mu}F'(x)hlm≃−μ1g=−μ1F′(x)

此时，退化为步长较小的梯度下降法。

当 μ\muμ 较小时，hlm==hgnh_{lm}==h_{gn}hlm==hgn, 在迭代后期，当x离x∗x^*x∗ 很近时，这是一个很好的步长，如果F(x*)=0(或非常小)，那么我们可以得到(几乎)二次最终收敛。
当 μ=0\mu=0μ=0 时，退化为高斯牛顿法。

故引入了一中 μ\muμ 的评价标准

ρ=F(x)−F(x+hlm)L(0)−L(hlm)\rho =\frac{F(x)-F(x+h_{lm})}{L(0)-L(h_{lm})}ρ=L(0)−L(hlm)F(x)−F(x+hlm)

当 0.25<ρ<0.750.25< \rho< 0.750.25<ρ<0.75 时，μ\muμ 趋近于1;
当 0.75<ρ0.75< \rho0.75<ρ 时，μ\muμ 小于 1；
当 ρ<0.25\rho< 0.25ρ<0.25 时，μ\muμ 大于1；

参考：
数学知识–Methods for Non-Linear Least Squares Problems（第三章）
信赖域(Trust Region)算法是怎么一回事
选主元的高斯-约当（Gauss-Jordan）消元法解线性方程组/求逆矩阵
LM算法
如何通俗易懂地解释「范数」

Levenberg–Marquardt（LM）相关推荐

如何计算给定一个unigram语言模型_CS224n笔记[5]:语言模型（LM）和循环神经网络（RNNs）...
CS224n笔记[5]:语言模型(LM)和循环神经网络(RNNs) 作者:郭必扬许久没更新了,十分惭愧,翻了翻之前的笔记,才之前上一期我们讲的是"依存分析".本期,我们介绍一下语 ...
NLP基础——语言模型（LM）
文章目录 NLP基础:语言模型(LM) 1. 模型评估(概率估计) 2. 平滑方法 3. LM在拼写纠正(Spell Correction)中的应用 NLP基础:语言模型(LM) 语言模型(LM,La ...
异方差的拉格朗日乘数（LM）检验python
异方差检验是用于判断数据是否存在异方差性的检验方法.在实际数据分析中,数据的方差有可能会随着自变量的变化而发生变化,这就导致了数据点之间的离散程度不同,使得数据的预测能力降低. 常见的异方差检验方法有 ...
计量经济分析：计量经济学中的三大检验（LR， Wald， LM）
前面用Python底层编写进行计量经济分析(一):多元线性回归(参数估计.T检验.拟合优度.F检验)写过在多元线性回归时的参数检验方法t检验和方程整体的F检验.在分析中和实际情况中,我们可能会假定因素 ...
利用函数wavread对语音信号进行采样_AI大语音（一）——语音识别基础（深度解析）...
1 声音特性声音(sound)是由物体振动产生的声波.是通过介质传播并能被人或动物听觉器官所感知的波动现象.最初发出振动的物体叫声源.声音以波的形式振动传播.声音是声波通过任何介质传播形成的运动. ...
对于谷歌应用传统的自动语音识别（ASR）系统的解析
目前,谷歌的各种语音搜索应用还在使用传统的自动语音识别(ASR)系统,它包括一个包括声学模型(AM ).一个发音模型(PM)和一个语言模型(LM),它们都是彼此独立训练的,而且需要研究人员在不同数据集 ...
DeepFlow高效的光流匹配算法（上）
本周主要介绍一篇基于传统光流法而改进的实现快速的稠密光流算法.该算法已经集成到OpenCV中,算法介绍网址:http://lear.inrialpes.fr/src/deepmatching/ 在介绍 ...
图形学中的光和辐射学（Radiometry）
最近重读Realtime Rendering 一书,对书中光照模型等内容有些感受:准备写些笔记记录一下自己的认识.本篇先主要讲讲辐射学(Radiometry)和简单光照模型:后面会讲讲BRDF及其实现 ...
语音识别（ASR）--语音转文字
语音识别(Automatic Speech Recognition)是以语音为研究对象,通过语音信号处理和模式识别让机器自动识别和理解人类口述的语.语音识别技术就是让机器通过识别和理解过程把语音信号转 ...

Levenberg–Marquardt（LM）