高斯牛顿法 Guass-Newton



求最优估计 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)等价于“解方程”问题：

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)中的优化问题变成

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)的指标函数可以改写成

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

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

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

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

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