高斯牛顿法 Guass-Newton

青旅半醒 2022-06-11 05:16 547阅读 0赞

求最优估计  *x* ∗   x∗，使得误差(残差)向量的  *ϵ*=*f*(*x* ∗ )−*z*  ϵ=f(x∗)−z的平方和  *S*(*x*)=*ϵ* *T* *ϵ*  S(x)=ϵTϵ最小，即求

*x* ∗ =argmin *x* *ϵ* *T* *ϵ*=argmin *x* *S*(*x*)=argmin *x* ∥*f*(*x*)−*z*∥ 2 2  (1)  (1)x∗=arg⁡minxϵTϵ=arg⁡minxS(x)=arg⁡minx‖f(x)−z‖22

最理想的情况，误差*ϵ*=0  ϵ=0，此时*f*(*x* ∗ )=*z*  f(x∗)=z，故最优化问题[(1)][1]  (1)等价于“解方程”问题：

*f*(*x*)=*z* (2)  (2)f(x)=z

大部分非线性最优化算法需要迭代求解

*x*:=*x*\+*δ*  x:=x+δ

若已有初值  *x*  x，每次迭代需要求解一个最优增量  *δ*  δ，使得迭代后的“预测输出”  *f*(*x*\+*δ*)  f(x+δ)与实际的观测值  *z*  z的误差  ∥*f*(*x*\+*δ*)−*z*∥  ‖f(x+δ)−z‖最小。

对误差向量线性化，有

*f*(*x*\+*δ*)−*z*≈*f*(*x*)\+*J**δ*−*z*  f(x+δ)−z≈f(x)+Jδ−z

其中，  *J*=∂*f*(*x*)∂*x*   J=∂f(x)∂x为  *f*(*x*)  f(x)的雅克比矩阵。令上式为  0  0，就变成了解线性方程问题

*J**δ*=*ϵ*,*ϵ*=*z*−*f*(*x*)  Jδ=ϵ,ϵ=z−f(x)

实际应用中，上式往往是超定的，故用线性最小二乘法求解，

*J* *T* *J**δ*=*J* *T* *ϵ*⇔*δ* ∗ =argmin *δ* ∥*J**δ*−*ϵ*∥ 2 2  (3)  (3)JTJδ=JTϵ⇔δ∗=arg⁡minδ‖Jδ−ϵ‖22

故非线性最优化问题变成了迭代求解线性方程的问题。上述算法又称为\\emph\{高斯牛顿法\}。

如果已知观测*z*  z的协方差的矩阵Σ  Σ，应该对指标函数按方差Σ  Σ加权，方差大的观测分量权重小，对结果的影响小。引入方差加权后，[(1)][1]  (1)中的优化问题变成

*x* ∗ =argmin *x* *S*(*x*)=argmin *x* *ϵ* *T* Σ −1 *ϵ*=argmin *x* ∥*f*(*x*)−*z*∥ 2 Σ −1   (4)  (4)x∗=arg⁡minxS(x)=arg⁡minxϵTΣ−1ϵ=arg⁡minx‖f(x)−z‖Σ−12

要解决上述问题，则在每一次迭代过程中求解

*δ* ∗ =argmin *δ* ∥*J**δ*−*ϵ*∥ 2 Σ −1   (5)  (5)δ∗=arg⁡minδ‖Jδ−ϵ‖Σ−12

设信息矩阵Σ −1   Σ−1有Cholesky分解

Σ −1 =*A* *T* *A* (6)  (6)Σ−1=ATA

则  [(5)][5]  (5)的指标函数可以改写成

∥*J**δ*−*ϵ*∥ 2 Σ −1  =(*J**δ*−*ϵ*) *T* *A* *T* *A*(*J**δ*−*ϵ*)=∥(*A**J*)*δ*−*A**ϵ*∥ 2 2  (7)  (7)‖Jδ−ϵ‖Σ−12=(Jδ−ϵ)TATA(Jδ−ϵ)=‖(AJ)δ−Aϵ‖22

只需令  *J* ~ =*A**J*  J~=AJ，  *ϵ* ~ =*A**ϵ*  ϵ~=Aϵ，  [(7)][7]  (7)变成  ∥*J* ~ *δ*−*ϵ* ~ ∥ 2 2   ‖J~δ−ϵ~‖22，故加权问题  [(5)][5]  (5)可以转化为非加权问题  [(3)][3]  (3)。利用  [(3)][3]  (3)中“  ⇔  ⇔”左右两侧的等价关系，问题  [(5)][5]  (5)中最优增量  *δ*  δ由下式确定：

*J* *T* Σ −1 *J**δ*=*J* *T* Σ −1 *ϵ* (8)  (8)JTΣ−1Jδ=JTΣ−1ϵ

在实际操作中，不需要进行Cholesky分解  [(6)][6]  (6)，只需在每一次迭代过程中求解  [(8)][8]  (8)即可。

可以证明，[(8)][8]  (8)中*J* *T* Σ −1 *J*  JTΣ−1J为能量函数*S*(*x*)  S(x)的Hessian矩阵；*J* *T* Σ −1   JTΣ−1为能量函数的梯度。

[1]: http://www.cnblogs.com/rongyilin/archive/2012/12/21/2827898.html#mjx-eqn-eqgen_opt
[5]: http://www.cnblogs.com/rongyilin/archive/2012/12/21/2827898.html#mjx-eqn-eqgen_opt_lin3
[7]: http://www.cnblogs.com/rongyilin/archive/2012/12/21/2827898.html#mjx-eqn-eqgen_opt_lin4
[3]: http://www.cnblogs.com/rongyilin/archive/2012/12/21/2827898.html#mjx-eqn-eqgen_opt_lin2
[6]: http://www.cnblogs.com/rongyilin/archive/2012/12/21/2827898.html#mjx-eqn-eqgen_opt_chol
[8]: http://www.cnblogs.com/rongyilin/archive/2012/12/21/2827898.html#mjx-eqn-eqgen_opt_lin_weighted