非线性最小二乘法优化
高斯-牛顿法
参考文章:[优化] Gauss-Newton非线性最小二乘算法
算法流程如下图(来自参考文章)所示:
接下来本文使用的数学符号意义与上图一样。其中
x
x
x是需要求解的参数,
f
(
x
)
f(x)
f(x)是一个残差向量。比如有一个优化问题,
y
=
a
s
i
n
(
w
t
+
b
)
+
c
y=asin(wt+b)+c
y=asin(wt+b)+c,给出m个数据
(
t
i
,
y
i
)
(
i
=
0
,
1
,
⋯
,
m
−
1
)
(t_i,y_i)(i=0,1,\cdots,m-1)
(ti,yi)(i=0,1,⋯,m−1),则
x
=
[
a
,
w
,
b
,
c
]
T
f
(
x
)
=
[
y
0
−
(
a
s
i
n
(
w
t
0
+
b
)
+
c
)
,
y
1
−
(
a
s
i
n
(
w
t
1
+
b
)
+
c
)
,
⋯
,
y
m
−
1
−
(
a
s
i
n
(
w
t
m
−
1
+
b
)
+
c
)
]
T
x=[a,w,b,c]^T \\ f(x)=[y_0-(asin(wt_0+b)+c),y_1-(asin(wt_1+b)+c),\cdots,y_{m-1}-(asin(wt_{m-1}+b)+c)]^T
x=[a,w,b,c]Tf(x)=[y0−(asin(wt0+b)+c),y1−(asin(wt1+b)+c),⋯,ym−1−(asin(wtm−1+b)+c)]T
则
∣
∣
f
(
x
)
∣
∣
2
||f(x)||^2
∣∣f(x)∣∣2(向量二范数)就是最小二乘法的损失值。
设损失函数
l
(
x
)
=
1
2
∣
∣
f
(
x
)
∣
∣
2
l(x)=\frac{1}{2}||f(x)||^2
l(x)=21∣∣f(x)∣∣2。
另外,
J
(
x
)
J(x)
J(x)为
f
(
x
)
f(x)
f(x)的雅可比矩阵,假设
x
x
x的长度为n,
f
(
x
)
f(x)
f(x)长度为m,则
J
(
X
)
J(X)
J(X)矩阵大小为(m,n)。
H
=
J
T
J
H=J^TJ
H=JTJ为
f
(
x
)
f(x)
f(x)的黑塞矩阵的近似矩阵。
B
=
−
J
T
f
(
x
)
B=-J^Tf(x)
B=−JTf(x)为损失函数
l
(
x
)
l(x)
l(x)(
1
2
\frac{1}{2}
21只是为了求导后约掉
∣
∣
f
(
x
)
∣
∣
2
||f(x)||^2
∣∣f(x)∣∣2的指数2)的负梯度
−
∂
l
(
x
)
∂
x
-\frac{\partial l(x)}{\partial x}
−∂x∂l(x)。
最速下降法
参考文章:
LM算法
在高斯-牛顿法中引入 μ \mu μ得到LM算法
引入 μ \mu μ的意义
- 高斯牛顿法的缺点
- H有可能不可逆
首先, H = J T J H=J^TJ H=JTJ为半正定对称矩阵(注:形如 A T A A^TA ATA(A为任意矩阵)都是半正定对称矩阵,这个定理是奇异值分解的基础),可以分解为 H = Q Λ Q T H=Q\Lambda Q^T H=QΛQT,其中矩阵 Q Q Q的每个列向量为 H H H的特征向量, Λ \Lambda Λ为对角矩阵,对角元素为对应特征向量的特征值。
因为 H H H为半正定对称矩阵,因此特征值有可能为0,因此不可逆。因为若H可逆,则 H − 1 = Q Λ − 1 Q T H^{-1}=Q\Lambda ^{-1}Q^T H−1=QΛ−1QT,其中 Λ − 1 \Lambda ^{-1} Λ−1对角元素为对应特征值 λ \lambda λ的倒数 1 λ \frac{1}{\lambda} λ1,因此若特征值为0,则 H H H不可逆。 - 步长
Δ
x
\Delta x
Δx可能过大,导致发散
由高斯牛顿法的算法流程可知,其核心是在点 x k x_k xk处利用 l ( x ) l(x) l(x)的泰勒展开,用二次多项式 p k ( x ) p_k(x) pk(x)(注:实际上 p k ( x ) p_k(x) pk(x)不是真正泰勒展开的二次多项式,因为矩阵 H H H只是黑塞矩阵的近似矩阵)近似 f ( x ) f(x) f(x)。
l ( x k + Δ x ) ≈ p k ( x k + Δ x ) = l ( x k ) + ( − B T ) Δ x + 1 2 Δ x T H Δ x l(x_k+\Delta x) \approx p_k(x_k+\Delta x)= l(x_k)+(-B^T)\Delta x+\frac{1}{2}{\Delta x}^T H \Delta x l(xk+Δx)≈pk(xk+Δx)=l(xk)+(−BT)Δx+21ΔxTHΔx
然后求二次多项式 p k ( x ) p_k(x) pk(x)的最小值点 x k + 1 = x k + argmin Δ x p k ( x k + Δ x ) x_{k+1}=x_{k}+\underset{\Delta x}{\operatorname{argmin}}\ {p_k(x_k+\Delta x)} xk+1=xk+Δxargmin pk(xk+Δx),然后 x k + 1 x_{k+1} xk+1则是这一次迭代的结果。
因此当 x k x_k xk与 p k ( x ) p_k(x) pk(x)的最小值点相距很远时,步长 Δ x \Delta x Δx会很大。但泰勒展开一般只在 x k x_k xk的局部区域内能很好的近似原始函数 l ( x ) l(x) l(x),因此步长太大算法可能会发散(损失值不降反升)。
- H有可能不可逆
- 引入(非负数)
μ
\mu
μ解决高斯牛顿法的缺点
- 步长
Δ
x
\Delta x
Δx太大的问题
步长可能太大,那么一个自然的想法就是正则化。因此,修改损失函数为:
p k ( x k + Δ x ) = l ( x k ) + ( − B T ) Δ x + 1 2 Δ x T H Δ x + 1 2 μ Δ x T Δ x p_k(x_k+\Delta x)= l(x_k)+(-B^T)\Delta x+\frac{1}{2}{\Delta x}^T H \Delta x+\frac{1}{2}\mu{\Delta x}^T \Delta x pk(xk+Δx)=l(xk)+(−BT)Δx+21ΔxTHΔx+21μΔxTΔx
正则化系数 μ \mu μ越大,则越能限制步长 Δ x \Delta x Δx的大小。
求解 argmin Δ x p k ( x k + Δ x ) \underset{\Delta x}{\operatorname{argmin}}\ {p_k(x_k+\Delta x)} Δxargmin pk(xk+Δx)的过程如下:
(1) 求导: ω ( Δ x ) = ∂ p k ( x k + Δ x ) ∂ Δ x = ( − B ) + H Δ x + μ Δ x = ( − B ) + ( H + μ I ) Δ x \omega (\Delta x)=\frac{\partial p_k(x_k+\Delta x)}{\partial \Delta x}=(-B)+H\Delta x+\mu \Delta x=(-B)+(H+\mu I) \Delta x ω(Δx)=∂Δx∂pk(xk+Δx)=(−B)+HΔx+μΔx=(−B)+(H+μI)Δx
(2) 令 ω ( Δ x ) = 0 \omega (\Delta x)=0 ω(Δx)=0得:
argmin Δ x p k ( x k + Δ x ) = ( H + μ I ) − 1 B \underset{\Delta x}{\operatorname{argmin}}\ {p_k(x_k+\Delta x)=(H+\mu I)^{-1}B } Δxargmin pk(xk+Δx)=(H+μI)−1B - H不可逆的问题
由上面可知现在 H H H变成了 ( H + μ I ) (H+\mu I) (H+μI),只要 μ > 0 \mu >0 μ>0,则 ( H + μ I ) (H+\mu I) (H+μI)一定可逆。因为:
(1) 首先 ( H + μ I ) (H+\mu I) (H+μI)是对称矩阵(保证了 ( H + μ I ) (H+\mu I) (H+μI)有n个正交特征向量,n为 x x x的长度, ( H + μ I ) (H+\mu I) (H+μI)大小为(n,n))。
(2) 其次 ( H + μ I ) (H+\mu I) (H+μI)与 H H H特征向量相同,并且:假设 H x = λ x Hx=\lambda x Hx=λx,则 ( H + μ I ) x = H x + μ x = ( λ + μ ) x (H+\mu I)x=Hx+\mu x=(\lambda +\mu)x (H+μI)x=Hx+μx=(λ+μ)x。所以 ( H + μ I ) (H+\mu I) (H+μI)的特征值为 H H H对应特征值加 μ \mu μ。又因为 λ ≥ 0 \lambda \ge 0 λ≥0,所以若 μ \mu μ大于0,则 ( H + μ I ) (H+\mu I) (H+μI)的特征值大于0。
(3)结合(1)(2)得若 μ > 0 \mu>0 μ>0,则 ( H + μ I ) (H+\mu I) (H+μI)为对称正定矩阵,所以 ( H + μ I ) (H+\mu I) (H+μI)可逆。
- 步长
Δ
x
\Delta x
Δx太大的问题
如何自动调整 μ \mu μ,LM与高斯牛顿法和最速下降法的关系,算法实现流程
- 如何自动调整
μ
\mu
μ,LM与高斯牛顿法和最速下降法的关系
参考文章:Levenberg–Marquardt算法学习- 其实信赖域法的本质就是看近似函数(比如这里就是泰勒展开的二阶形式)的损失值下降量 Δ L k \Delta L_{k} ΔLk和实际损失函数的损失值下降量 Δ F k \Delta F_{k} ΔFk的比值,如果 Δ F k Δ L k \frac{\Delta F_{k}}{\Delta L_{k}} ΔLkΔFk 约等于1说明近似函数在步长 Δ k \Delta_{k} Δk内与实际损失函数很近似,可以保持这个步长或者扩大步长,否则若 Δ F k Δ L k \frac{\Delta F_{k}}{\Delta L_{k}} ΔLkΔFk约等于0甚至是负数,就缩小步长。(需要保证 Δ L k > 0 \Delta L_{k}>0 ΔLk>0)
- 算法实现流程
参考文章:A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar-
注意: 里面的伪代码中有点错误,g应该是负梯度,也就是 g : = − J T ϵ p g:=-J^T \epsilon_{p} g:=−JTϵp。
-
ϵ T Σ y − 1 ϵ \epsilon^T\Sigma_y^{-1}\epsilon ϵTΣy−1ϵ的作用
参考文章 A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar中提到了这样一段话:
注意,这篇文章里的向量 x x x(是本文中的真实值 y = [ y 0 , y 1 , ⋯ , y m − 1 ] T y=[y_0,y_1,\cdots,y_{m-1}]^T y=[y0,y1,⋯,ym−1]T)与本文的 x x x意义不一样。因此下面本文用 y y y代替这篇文章的 x x x。 ϵ T Σ y − 1 ϵ \epsilon^T\Sigma_y^{-1}\epsilon ϵTΣy−1ϵ的作用是消除不同 y i y_i yi有可能有不同量级的影响。
我们假设 Σ y \Sigma_y Σy为对角矩阵,也就是 y i y_{i} yi之间相互独立,则对角值 σ i \sigma_i σi为 y i y_i yi的方差, σ i \sigma_i σi表示了 y i y_i yi的变化范围(可以理解为量级)。量级越大,那么对应误差 ϵ i \epsilon_i ϵi的值变化范围也会大,因此在优化过程中会重点优化 ϵ i \epsilon_i ϵi。因此我们要避免这种由量级导致的误差过大或过小。因此算法以 ϵ T Σ y − 1 ϵ \epsilon^T\Sigma_y^{-1}\epsilon ϵTΣy−1ϵ作为损失值,代替 ϵ T ϵ \epsilon^T \epsilon ϵTϵ。
若 Σ y \Sigma_y Σy不是对角矩阵,但因为协方差矩阵和协方差矩阵的逆都是正定对称矩阵(只要没有互相关变量)。因此 Σ y − 1 \Sigma_y^{-1} Σy−1可分解为 Q Λ − 1 Q T Q\Lambda^{-1} Q^T QΛ−1QT。而 ϵ T Σ y − 1 ϵ \epsilon^T\Sigma_y^{-1}\epsilon ϵTΣy−1ϵ= ( ϵ T Q ) Λ − 1 ( Q T ϵ ) (\epsilon^TQ)\Lambda^{-1} (Q^T\epsilon) (ϵTQ)Λ−1(QTϵ),把 ( Q T ϵ ) (Q^T\epsilon) (QTϵ)当成新的随机变量。而 ( Q T ϵ ) (Q^T\epsilon) (QTϵ)的协方差矩阵为 Λ \Lambda Λ,因此也实现了消除量级影响。 -
μ \mu μ初始值
在参考文章A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar的伪代码里, μ \mu μ的初始值如下图所示。其中 m a x i = 1 , ⋯ , m ( H i i ) max_{i=1,\cdots,m}(H_{ii}) maxi=1,⋯,m(Hii)(参考文章的 A A A等于本文的 H H H)。这其实是为了让 μ \mu μ和 H H H对角线上的值的数量级一致。因为我们有 H + μ I H+\mu I H+μI,因此 μ \mu μ是加到 H H H的对角线上的。
-
参考文章建议的初始值 :
-