LM优化算法,是一种非线性优化算法,其可以看作是梯度下降和高斯牛顿法的结合。综合了梯度下降对初值不敏感和高斯牛顿在最优值附近收敛速度快的特点。
本人非数学专业,且对算法理解可能不到位,详细的算法推导及各个优化算法之间的关系,非常推荐看METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS ,其介绍更详细也更专业。
以下简要介绍LM算法,然后给出opencv中有关的实现。
本文代码也是来自opencv相关代码,位置在。
当然,有关LM算法的实现,在opencv的usac模块中有更为正式的实现,在该模块中,对许多算法进行了集成和优化,后续再进行研究。
对于一个最小二乘问题,有:
F
(
x
)
=
1
2
∑
i
=
1
m
(
f
i
(
x
)
)
2
f
i
(
x
)
为
残
差
F(x) = \frac12\sum_{i=1}^m(f_i(x))^2\\ f_i(x) 为残差
F(x)=21?i=1∑m?(fi?(x))2fi?(x)为残差
我们的目的是求解残差f(x)的一组系数,使目标函数(或代价函数)F(x)达到最小值。通常的做法是,给定一个初值x0,优化算法不断的迭代,产生x1,x2,…,最终收敛于X。
因此,在这个收敛的过程中,我们需要两样东西,即方向h和步长α
x
k
+
1
=
x
k
+
α
h
d
h
d
即
为
收
敛
的
方
向
α
即
为
收
敛
的
步
长
x_{k+1} = x_k + \alpha h_d \\ h_d即为收敛的方向 \\ \alpha即为收敛的步长
xk+1?=xk?+αhd?hd?即为收敛的方向α即为收敛的步长
所以,我们看到的梯度下降、高斯牛顿以及本文说的LM算法,其思想都是一致的。
梯度下降
这里给出梯度下降的公式,有:
x
k
+
1
=
x
k
?
α
F
˙
(
x
)
F
˙
(
x
)
为
F
(
x
)
一
阶
导
数
x_{k+1} = x_k - \alpha \dot{F}(x) \\ \dot{F}(x)为F(x)一阶导数
xk+1?=xk??αF˙(x)F˙(x)为F(x)一阶导数
高斯牛顿
高斯牛顿是在牛顿迭代的基础上,用雅可比矩阵的平方来近似Hessian矩阵,求解起来将会更加高效(大多数情况下,Hessian矩阵太难求)。这样做的缺点就是,要求迭代初值必须在最优解附近,否则以上假设将不成立。这里给出其公式:
x
k
+
1
=
x
k
+
α
h
g
n
J
T
J
h
g
n
=
?
J
T
f
x_{k+1} = x_k + \alpha h_{gn} \\ J^TJh_{gn} = -J^Tf \\
xk+1?=xk?+αhgn?JTJhgn?=?JTf
Levenberg–Marquardt
给出其公式:
x
k
+
1
=
x
k
+
α
h
l
m
(
J
T
J
+
μ
I
)
h
l
m
=
?
J
T
f
x_{k+1} = x_k + \alpha h_{lm} \\ (J^TJ + \mu I)h_{lm} = -J^Tf
xk+1?=xk?+αhlm?(JTJ+μI)hlm?=?JTf
在有些情况下,JtJ不可逆,导致高斯牛顿无法使用。LM通过将JtJ加上一个μI(I为单位阵),保证了正规方程的左侧一定可逆。
在算法实际的应用中,通过调节μ的大小,可以使算法在梯度下降和高斯牛顿两者之间切换,如:
- 使用较大的μ值,算法可以看作梯度下降,适合当前优化位置距离最优解较远的情况
- 使用较小的μ值,算法可以看作高斯牛顿,适合当前优化位置接近最优解。
LM算法通过求解正规方程,得到每次迭代的步长和方向。因此,算法的主要目的是求解正规方程左右侧的项,然后通过SVD分解即可得到参数更新的步长及方向,即:
- JtJ - 雅可比矩阵的转置与其自己相乘
- JErr - 雅可比矩阵的转置与残差矩阵(向量)相乘
算法框架在外部计算雅可比矩阵和残差矩阵,然后在算法内部通过求解正规方程,得到每次参数更新的步长及方向。
然后利用更新的参数,在外部计算残差,然后判断残差是否朝着收敛方向进行。
算法的整体流程
求解器update流程
求解器updata()内部根据不同的状态进行相应的计算,具体流程如下:
说明
- J表示雅可比矩阵
- err表示残差矩阵
- 求解器内外同步表示求解器内部和外部相应的参数指向相同,便于在求解器外部进行雅可比矩阵和残差矩阵的计算
算法的实现主要步骤都体现在上述流程图中。