数值计算之 Levenberg-Marquardt算法
前言
本篇是牛顿法的最后一篇，Levenberg-Marquardt算法，也就是阻尼牛顿法。
高斯牛顿法
回顾上篇中的高斯牛顿法：
f
(
x
0
+
Δ
x
)
=
f
(
x
0
)
+
J
(
x
0
)
Δ
x
min
?
Δ
x
1
2
∣
∣
f
(
x
)
∣
∣
2
=
min
?
Δ
x
1
2
(
f
(
x
0
)
+
J
(
x
0
)
Δ
x
)
T
(
f
(
x
0
)
+
J
(
x
0
)
Δ
x
)
J
T
J
Δ
x
=
?
J
T
f
f({\bf x_0+\Delta x}) = f({\bf x_0}) + J({\bf x_0}){\bf \Delta x} \\ \quad \\ \min_{\bf \Delta x} \frac{1}{2} ||f({\bf x})||^2 \\ = \min_{\bf \Delta x} \frac{1}{2} (f({\bf x_0})+J({\bf x_0}){\bf \Delta x})^T(f({\bf x_0})+J({\bf x_0}){\bf \Delta x}) \\ \quad \\ J^TJ{\bf \Delta x} = -J^Tf
f(x0?+Δx)=f(x0?)+J(x0?)ΔxΔxmin?21?∣∣f(x)∣∣2=Δxmin?21?(f(x0?)+J(x0?)Δx)T(f(x0?)+J(x0?)Δx)JTJΔx=?JTf
将牛顿法中的海森矩阵
H
H
H，使用一阶梯度
J
T
J
J^TJ
JTJ进行了替换，因而避免了计算二阶梯度。然而，
J
T
J
J^TJ
JTJ并非正定，可能是奇异阵或者病态矩阵，求逆不稳定，使得迭代的稳定性和收敛性差。
LM算法
阻尼牛顿法
为了提升迭代的稳定性，在优化函数上添加阻尼（机器学习里有时称为正则项），来降低迭代步长：
min
?
Δ
x
(
1
2
∣
∣
f
(
x
)
∣
∣
2
+
1
2
λ
Δ
x
T
Δ
x
)
=
min
?
Δ
x
1
2
(
f
(
x
0
)
+
J
(
x
0
)
Δ
x
)
T
(
f
(
x
0
)
+
J
(
x
0
)
Δ
x
)
+
1
2
λ
Δ
x
T
Δ
x
)
\min_{\bf \Delta x} (\frac{1}{2} ||f({\bf x})||^2+\frac{1}{2}\lambda{\bf \Delta x^T\Delta x}) \\ = \min_{\bf \Delta x} \frac{1}{2} (f({\bf x_0})+J({\bf x_0}){\bf \Delta x})^T(f({\bf x_0})+J({\bf x_0}){\bf \Delta x}) + \frac{1}{2} \lambda {\bf \Delta x^T\Delta x}) \\
Δxmin?(21?∣∣f(x)∣∣2+21?λΔxTΔx)=Δxmin?21?(f(x0?)+J(x0?)Δx)T(f(x0?)+J(x0?)Δx)+21?λΔxTΔx)
在
1
2
∣
∣
f
(
x
)
∣
∣
2
\frac{1}{2}||f({\bf x})||^2
21?∣∣f(x)∣∣2后添加
1
2
λ
Δ
x
T
Δ
x
\frac {1}{2}\lambda {\bf \Delta x}^T {\bf \Delta x}
21?λΔxTΔx，当求出的
Δ
x
\bf \Delta x
Δx比较大时，对函数值进行惩罚。
还是按照高斯牛顿法一样求一阶梯度等于零，获得迭代表达式：
J
T
J
Δ
x
+
λ
Δ
x
=
?
J
T
f
(
J
T
J
+
λ
I
)
Δ
x
=
?
J
T
f
Δ
x
=
?
(
J
T
J
+
λ
I
)
?
1
J
T
f
J^TJ{\bf \Delta x} +\lambda {\bf \Delta x} = -J^Tf \\ (J^TJ+\lambda I){\bf \Delta x}=-J^Tf \\ {\bf \Delta x} = -(J^TJ+\lambda I)^{-1}J^Tf \\
JTJΔx+λΔx=?JTf(JTJ+λI)Δx=?JTfΔx=?(JTJ+λI)?1JTf
可以看出，阻尼牛顿法通过
J
T
J
+
λ
I
J^TJ+\lambda I
JTJ+λI替换
H
H
H，一方面避免二阶梯度计算，另一方面保证了矩阵的可逆性。
阻尼系数
λ
\lambda
λ
上面的惩罚项中，
λ
\lambda
λ用于控制
Δ
x
\bf \Delta x
Δx对迭代步长的影响。当
f
(
x
0
+
Δ
x
)
f(\bf x_0 +\Delta x)
f(x0?+Δx)与一阶泰勒展开
f
(
x
0
)
+
J
(
x
0
)
Δ
x
f({\bf x_0}) +J({\bf x_0}){\bf\Delta x}
f(x0?)+J(x0?)Δx足够近似，迭代步长可以大一些，加速收敛；如果函数值与泰勒展开不够近似，迭代步长要小一些，增加稳定。
因此，可以设计以下策略调整阻尼系数：
ρ
=
f
(
x
0
+
Δ
x
)
?
f
(
x
0
)
J
(
x
0
)
Δ
x
0
i
f
ρ
>
3
4
:
λ
=
1
2
λ
i
f
ρ
<
1
2
:
λ
=
2
λ
\rho = \frac {f({\bf x_0 + \Delta x})-f({\bf x_0})} {J({\bf x_0}){\bf \Delta x_0}} \\ \quad \\ if \quad \rho > \frac{3}{4}: \quad \lambda=\frac{1}{2}\lambda \\ \quad \\ if \quad \rho < \frac{1}{2}: \quad \lambda= 2\lambda
ρ=J(x0?)Δx0?f(x0?+Δx)?f(x0?)?ifρ>43?:λ=21?λifρ<21?:λ=2λ
置信域法
LM算法还可以通过置信域思想实现：
min
?
Δ
x
1
2
∣
∣
f
(
x
)
∣
∣
2
s
.
t
.
Δ
x
T
Δ
x
≤
d
\min_{\bf \Delta x} \frac{1}{2} ||f({\bf x})||^2 \\ \quad \\ s.t. \quad {\bf \Delta x}^T{\bf \Delta x} \le d
Δxmin?21?∣∣f(x)∣∣2s.t.ΔxTΔx≤d
通过限制
Δ
x
\bf \Delta x
Δx的长度来约束迭代步长。这时可以通过拉格朗日乘子法把有约束的优化问题转换成无约束的优化问题：
1
2
∣
∣
f
(
x
0
)
+
J
(
x
0
)
Δ
x
∣
∣
2
+
1
2
λ
(
Δ
x
T
Δ
x
?
d
)
\frac{1}{2}||f({\bf x_0})+J({\bf x_0}){\bf \Delta x}||^2+\frac{1}{2}\lambda({\bf \Delta x}^T{\bf \Delta x} - d)
21?∣∣f(x0?)+J(x0?)Δx∣∣2+21?λ(ΔxTΔx?d)
这个优化问题涉及到KKT收敛性等等判定，最后的迭代表示式形式与阻尼牛顿法相同
置信域
d
d
d
置信域法中的
λ
\lambda
λ是求解拉格朗日乘子法计算出来的，而我们需要调节的是
Δ
x
\bf \Delta x
Δx的置信域
d
d
d：
ρ
=
f
(
x
0
+
Δ
x
)
?
f
(
x
0
)
J
(
x
0
)
Δ
x
0
i
f
ρ
>
3
4
:
d
=
2
d
i
f
ρ
<
1
2
:
d
=
1
2
d
\rho = \frac {f({\bf x_0 + \Delta x})-f({\bf x_0})} {J({\bf x_0}){\bf \Delta x_0}} \\ \quad \\ if \quad \rho > \frac{3}{4}: \quad d=2d \\ \quad \\ if \quad \rho < \frac{1}{2}: \quad d=\frac{1}{2}d
ρ=J(x0?)Δx0?f(x0?+Δx)?f(x0?)?ifρ>43?:d=2difρ<21?:d=21?d
LM算法的流程
最后给出LM算法的整体流程：
- 根据迭代表达式计算
Δ
x
\bf \Delta x
Δx
- 计算
ρ
\rho
ρ，判断是否满足终止条件，不满足则计算
x
=
Δ
x
+
x
0
\bf x=\Delta x+ x_0
x=Δx+x0?，否则迭代结束
- 调节阻尼系数
λ
\lambda
λ，或者调节置信域
d
d
d，返回步骤1
|