最小二乘问题和非线性优化

最小二乘问题和非线性优化

    • 0.引言
    • 1.最小二乘问题
    • 2.迭代下降法
    • 3.最速下降法
    • 4.牛顿法
    • 5.阻尼法
    • 6.高斯牛顿(GN)法
    • 7.莱文贝格马夸特(LM)法
    • 8.鲁棒核函数

0.引言

转载自此处,修正了一点小错误。

1.最小二乘问题

在求解 SLAM 中的最优状态估计问题时,我们一般会得到两个变量,一个是由传感器获得的实际观测值 z \boldsymbol{z} z,一个是根据目前估计的状态量和观测模型计算出来的预测值 h ( x ) h(\boldsymbol{x}) h(x)。求解最优状态估计问题时通常我们会尝试最小化预测值和观测值计算出的残差平方(使用平方是为了统一正负号的影响),即求解以下最小二乘问题:

x ∗ = arg ⁡ min ⁡ x ∣ ∣ z − h ( x ) ∣ ∣ 2 \boldsymbol{x}^* = \arg\min_{\boldsymbol{x}}||\boldsymbol{z} - h(\boldsymbol{x})||^2 x=argxmin∣∣zh(x)2
如果观测模型是线性模型,则上式转化为线性最小二乘问题:

x ∗ = arg ⁡ min ⁡ x ∣ ∣ z − H x ∣ ∣ 2 \boldsymbol{x}^* = \arg\min_{\boldsymbol{x}}||\boldsymbol{z} - \boldsymbol{H}\boldsymbol{x}||^2 x=argxmin∣∣zHx2
对于线性最小二乘问题,我们可以直接求闭式解: x ∗ = − ( H T H ) − 1 H T z \boldsymbol{x}^* = -(\boldsymbol{H}^T\boldsymbol{H})^{-1}\boldsymbol{H}^T\boldsymbol{z} x=(HTH)1HTz 这里不进行赘述。

实际的问题中,我们通常要最小化不止一个残差,不同残差通常会按其重要性(不确定性)分配一个相应的权重系数,且观测模型也通常是非线性的,即求解以下问题:

e i ( x ) = z i − h i ( x ) i = 1 , 2 , . . . , n ∣ ∣ e i ( x ) ∣ ∣ Σ i 2 = e i T Σ i e i x ∗ = arg ⁡ min ⁡ x F ( x ) = arg ⁡ min ⁡ x ∑ i ∣ ∣ e i ( x ) ∣ ∣ Σ i 2 \begin{aligned} \boldsymbol{e}_i(\boldsymbol{x}) &= \boldsymbol{z}_i - h_i(\boldsymbol{x}) \qquad i = 1, 2, ..., n\\ ||\boldsymbol{e}_i(\boldsymbol{x})||^2_{\Sigma_i} &= \boldsymbol{e}_i^T\boldsymbol{\Sigma}_i\boldsymbol{e}_i\\ \boldsymbol{x}^* &= \arg\min_{\boldsymbol{x}}F(\boldsymbol{x})\\ &= \arg\min_{\boldsymbol{x}}\sum_i||\boldsymbol{e}_i(\boldsymbol{x})||^2_{\Sigma_i} \end{aligned} ei(x)∣∣ei(x)Σi2x=zihi(x)i=1,2,...,n=eiTΣiei=argxminF(x)=argxmini∣∣ei(x)Σi2
我们想要获得一个状态量 x ∗ \boldsymbol{x}^* x,使得损失函数 F ( x ) F(\boldsymbol{x}) F(x) 取得局部最小值。

在具体求解之前,先考虑 F ( x ) F(\boldsymbol{x}) F(x) 的性质,对其进行二阶泰勒展开:

F ( x + Δ x ) = F ( x ) + J Δ x + 1 2 Δ x T H Δ x + O ( ∣ ∣ Δ x ∣ ∣ 3 ) F(\boldsymbol{x} + \Delta\boldsymbol{x}) = F(\boldsymbol{x}) + \boldsymbol{J}\Delta\boldsymbol{x} + \frac{1}{2}\Delta\boldsymbol{x}^T\boldsymbol{H}\Delta\boldsymbol{x} + O(||\Delta\boldsymbol{x}||^3) F(x+Δx)=F(x)+JΔx+21ΔxTHΔx+O(∣∣Δx3)
忽略高阶余项,对二次函数,有以下性质:

如果在点 x ∗ \boldsymbol{x}^* x 处,导数为 0 \boldsymbol{0} 0,则这个点为稳定点,根据 Hessian 矩阵的正定性有不同性质:

  • 如果 H \boldsymbol{H} H 为正定矩阵,则 F ( x ∗ ) F(\boldsymbol{x}^*) F(x) 为局部最小值
  • 如果 H \boldsymbol{H} H 为负定矩阵,则 F ( x ∗ ) F(\boldsymbol{x}^*) F(x) 为局部最大值
  • 如果 H \boldsymbol{H} H 为不定矩阵,则 F ( x ∗ ) F(\boldsymbol{x}^*) F(x) 为鞍点

在实际过程中,一般 F ( x ) F(x) F(x) 比较复杂,我们没有办法直接使其导数为 0 进而求出该点,因此常用的是迭代法,即找到一个下降方向使损失函数随 x \boldsymbol{x} x 的迭代逐渐减少,直到 x \boldsymbol{x} x 收敛到 x ∗ \boldsymbol{x}^* x。下面整理以下常用的几种迭代方法。

2.迭代下降法

上面提到,我们需要找到一个 x \boldsymbol{x} x 的迭代量使得 F ( x ) F(\boldsymbol{x}) F(x) 减少。这个过程分两步:

找到 F ( x ) \boldsymbol{F(x)} F(x) 的下降方向,构建该方向的单位向量 d \boldsymbol{d} d
确定该方向的迭代步长 α \alpha α
迭代后的自变量 x + α d \boldsymbol{x} + \alpha\boldsymbol{d} x+αd 对应的函数值可以用一阶泰勒展开近似(当步长足够小的时候):

F ( x + α d ) = F ( x ) + α J d F(\boldsymbol{x} + \alpha\boldsymbol{d}) = F(\boldsymbol{x}) + \alpha\boldsymbol{Jd} F(x+αd)=F(x)+αJd
因此,不难发现,要保持 F ( x ) F(x) F(x) 是下降的,只需要保证 J d < 0 \boldsymbol{Jd} < 0 Jd<0。以下几种方法都是以不同思路在寻找一个合适的方向进行迭代。

3.最速下降法

基于上一部分,我们知道变化量为 α J d \alpha\boldsymbol{Jd} αJd,其中 J d = ∣ ∣ J ∣ ∣ cos ⁡ θ \boldsymbol{Jd} = ||\boldsymbol{J}||\cos{\theta} Jd=∣∣J∣∣cosθ θ \theta θ 为梯度 J \boldsymbol{J} J d \boldsymbol{d} d 的夹角。当 θ = − π \theta = -\pi θ=π 时, J d = − ∣ ∣ J ∣ ∣ \boldsymbol{Jd} = -||\boldsymbol{J}|| Jd=∣∣J∣∣ 取得最小值,此时方向向量为:

d = − J T ∣ ∣ J ∣ ∣ \boldsymbol{d} = \frac{-\boldsymbol{J}^T}{||\boldsymbol{J}||} d=∣∣J∣∣JT
因此,沿梯度 J \boldsymbol{J} J 的负方向可以使 F ( x ) F(\boldsymbol{x}) F(x),但实际过程中,一般只会在迭代刚开始使用这种方法,当接近最优值时这种方法会出现震荡并且收敛较慢。

4.牛顿法

F ( x ) F(\boldsymbol{x}) F(x) 进行二阶泰勒展开有:

F ( x + Δ x ) = F ( x ) + J Δ x + 1 2 Δ x T H Δ x F(\boldsymbol{x} + \Delta\boldsymbol{x}) = F(\boldsymbol{x}) + \boldsymbol{J}\Delta\boldsymbol{x} + \frac{1}{2}\Delta\boldsymbol{x}^T\boldsymbol{H}\Delta\boldsymbol{x} F(x+Δx)=F(x)+JΔx+21ΔxTHΔx
我们关注的是求一个 Δ x \Delta\boldsymbol{x} Δx 使得 J Δ x + 1 2 Δ x T H Δ x \boldsymbol{J}\Delta\boldsymbol{x} + \frac{1}{2}\Delta\boldsymbol{x}^T\boldsymbol{H}\Delta\boldsymbol{x} JΔx+21ΔxTHΔx 最小,因此可以求导得:

∂ ( J Δ x + 1 2 Δ x T H Δ x ) ∂ Δ x = J T + H Δ x = 0 ⇒ Δ x = − H − 1 J T \begin{aligned} \frac{\partial(\boldsymbol{J}\Delta\boldsymbol{x} + \frac{1}{2}\Delta\boldsymbol{x}^T\boldsymbol{H}\Delta\boldsymbol{x})}{\partial\Delta\boldsymbol{x}} &= \boldsymbol{J}^T + \boldsymbol{H}\Delta\boldsymbol{x} = 0\\ \Rightarrow \Delta\boldsymbol{x} &= -\boldsymbol{H}^{-1}\boldsymbol{J}^T \end{aligned} Δx(JΔx+21ΔxTHΔx)Δx=JT+HΔx=0=H1JT
H \boldsymbol{H} H 为正定矩阵且当前 x \boldsymbol{x} x 在最优点附近时,取 Δ x = − H − 1 J T \Delta\boldsymbol{x} = -\boldsymbol{H}^{-1}\boldsymbol{J}^T Δx=H1JT 可以使函数获得局部最小值。但缺点是残差的 Hessian 函数通常比较难求。

5.阻尼法

在牛顿法的基础上,为了控制每次迭代不要太激进,我们可以在损失函数中增加一项惩罚项,如下所示:

arg ⁡ min ⁡ Δ x { F ( x ) + J Δ x + 1 2 Δ x T H Δ x + 1 2 μ ( Δ x ) T ( Δ x ) } \arg\min_{\Delta\boldsymbol{x}}\left\{F(\boldsymbol{x}) + \boldsymbol{J}\Delta\boldsymbol{x} + \frac{1}{2}\Delta\boldsymbol{x}^T\boldsymbol{H}\Delta\boldsymbol{x} + \frac{1}{2}\mu(\Delta\boldsymbol{x})^T(\Delta\boldsymbol{x})\right\} argΔxmin{F(x)+JΔx+21ΔxTHΔx+21μ(Δx)T(Δx)}
当我们选的 Δ x \Delta\boldsymbol{x} Δx 太大时,损失函数也会变大,变大的幅度由 μ \mu μ 决定,因此我们可以控制每次迭代量 Δ x \Delta\boldsymbol{x} Δx 的大小。同样在右侧部分对 Δ x \Delta\boldsymbol{x} Δx 求导有:

J T + H Δ x + μ Δ x = 0 ( H + μ I ) Δ x = − J T \begin{aligned} \boldsymbol{J}^T + \boldsymbol{H}\Delta\boldsymbol{x} + \mu\Delta\boldsymbol{x} &= 0\\ (\boldsymbol{H} + \mu\boldsymbol{I})\Delta\boldsymbol{x} = -\boldsymbol{J}^T \end{aligned} JT+HΔx+μΔx(H+μI)Δx=JT=0
这个思路我们后面在介绍 LM 法时也会用到。

6.高斯牛顿(GN)法

在前面的整理中实际上我们求解的是一系列残差的和,求单个残差的雅可比比较简单,因此在后续的几种方法中我们关注各个残差的变化。将上述非线性最小二乘问题中的各个残差写成向量形式:

F ( x ) = E ( x ) = [ e 1 ( x ) e 2 ( x ) . . . e n ( x ) ] \boldsymbol{F}(\boldsymbol{x}) =\boldsymbol{E}(\boldsymbol{x}) =\begin{bmatrix} \boldsymbol{e}_1(\boldsymbol{x})\\ \boldsymbol{e}_2(\boldsymbol{x})\\ ...\\ \boldsymbol{e}_n(\boldsymbol{x})\\ \end{bmatrix} F(x)=E(x)= e1(x)e2(x)...en(x)

e ( x ) \boldsymbol{e}(\boldsymbol{x}) e(x) 进行泰勒展开,有:

e ( x + Δ x ) = e ( x ) + J Δ x \boldsymbol{e}(\boldsymbol{x} + \Delta\boldsymbol{x}) = \boldsymbol{e}(\boldsymbol{x}) + \boldsymbol{J}\Delta\boldsymbol{x} e(x+Δx)=e(x)+JΔx
上式中, J \boldsymbol{J} J 是残差矩阵 e ( x ) \boldsymbol{e}(\boldsymbol{x}) e(x) 对状态量的雅可比矩阵。

注意到在原来的线性最小二乘问题中,对每个残差还有一个权重矩阵 Σ \boldsymbol{\Sigma} Σ。这种情况下,我们只需要令 e i ( x ) = Σ i e i ( x ) \boldsymbol{e}_i(\boldsymbol{x}) = \sqrt{\boldsymbol{\Sigma}_i}\boldsymbol{e}_i(\boldsymbol{x}) ei(x)=Σi ei(x) 即可。因此下式中暂不考虑 Σ \boldsymbol{\Sigma} Σ 的影响。

在这种形式下对 e ( x ) \boldsymbol{e}(\boldsymbol{x}) e(x) 进行泰勒展开,有:

∣ ∣ e ( x + Δ x ) ∣ ∣ 2 = e ( x + Δ x ) T e ( x + Δ x ) = ( e ( x ) + J Δ x ) T ( e ( x ) + J Δ x ) = e ( x ) T e ( x ) + Δ x T J T e ( x ) + e ( x ) T J Δ x + Δ x T J T J Δ x \begin{aligned} ||e(\boldsymbol{x} + \Delta\boldsymbol{x}) ||^2&= \boldsymbol{e}(\boldsymbol{x} + \Delta\boldsymbol{x})^T\boldsymbol{e}(\boldsymbol{x} + \Delta\boldsymbol{x})\\ &= (\boldsymbol{e}(\boldsymbol{x}) + \boldsymbol{J}\Delta\boldsymbol{x})^T(\boldsymbol{e}(\boldsymbol{x}) + \boldsymbol{J}\Delta\boldsymbol{x})\\ &= \boldsymbol{e}(\boldsymbol{x})^T\boldsymbol{e}(\boldsymbol{x}) + \Delta\boldsymbol{x}^T\boldsymbol{J}^T\boldsymbol{e}(\boldsymbol{x}) + \boldsymbol{e}(\boldsymbol{x})^T\boldsymbol{J}\Delta\boldsymbol{x} + \Delta\boldsymbol{x}^T\boldsymbol{J}^T\boldsymbol{J}\Delta\boldsymbol{x} \end{aligned} ∣∣e(x+Δx)2=e(x+Δx)Te(x+Δx)=(e(x)+JΔx)T(e(x)+JΔx)=e(x)Te(x)+ΔxTJTe(x)+e(x)TJΔx+ΔxTJTJΔx
上式中,注意到: e ( x ) e(\boldsymbol{x}) e(x) 是一维的,因此 Δ x T J T e ( x ) = e ( x ) T J Δ x \Delta\boldsymbol{x}^T\boldsymbol{J}^T\boldsymbol{e}(\boldsymbol{x}) = \boldsymbol{e}(\boldsymbol{x})^T\boldsymbol{J}\Delta\boldsymbol{x} ΔxTJTe(x)=e(x)TJΔx,因此化简得:

F ( x + Δ x ) = e ( x ) T e ( x ) + 2 e ( x ) T J Δ x + Δ x T J T J Δ x = F ( x ) + 2 e ( x ) T J Δ x + Δ x T J T J Δ x \begin{aligned} F(\boldsymbol{x} + \Delta\boldsymbol{x}) &= \boldsymbol{e}(\boldsymbol{x})^T\boldsymbol{e}(\boldsymbol{x}) + 2\boldsymbol{e}(\boldsymbol{x})^T\boldsymbol{J}\Delta\boldsymbol{x} + \Delta\boldsymbol{x}^T\boldsymbol{J}^T\boldsymbol{J}\Delta\boldsymbol{x}\\ &= F(\boldsymbol{x}) + 2\boldsymbol{e}(\boldsymbol{x})^T\boldsymbol{J}\Delta\boldsymbol{x} + \Delta\boldsymbol{x}^T\boldsymbol{J}^T\boldsymbol{J}\Delta\boldsymbol{x} \end{aligned} F(x+Δx)=e(x)Te(x)+2e(x)TJΔx+ΔxTJTJΔx=F(x)+2e(x)TJΔx+ΔxTJTJΔx
通过这种方式,我们同样将其近似一个二次函数,并且和我们之前展开的结果比较,不难发现在这里我们实际上是用 J T e \boldsymbol{J}^T\boldsymbol{e} JTe 来近似 Jacobian 矩阵,用 J T J \boldsymbol{J}^T\boldsymbol{J} JTJ 来近似 Hessian 矩阵。因此,当 J \boldsymbol{J} J 满秩时,我们可以保证在上式导数为 0 的地方可以确保函数取得局部最小值。同样,在上式右侧部分对 Δ x \Delta\boldsymbol{x} Δx 求导并使其为 0 有:

J T e ( x ) + J T J Δ x = 0 ⇒ J T J Δ x = − J T e ( x ) ⇒ H Δ x = b \begin{aligned} \boldsymbol{J}^T\boldsymbol{e}(\boldsymbol{x}) + \boldsymbol{J}^T\boldsymbol{J}\Delta\boldsymbol{x} = 0\\ \Rightarrow \boldsymbol{J}^T\boldsymbol{J}\Delta\boldsymbol{x} &= -\boldsymbol{J}^T\boldsymbol{e}(\boldsymbol{x})\\ \Rightarrow \boldsymbol{H}\Delta\boldsymbol{x} &= \boldsymbol{b} \end{aligned} JTe(x)+JTJΔx=0JTJΔxHΔx=JTe(x)=b
上式中,我们令 H = J T J , b = − J T e \boldsymbol{H} = \boldsymbol{J}^T\boldsymbol{J}, \boldsymbol{b} = -\boldsymbol{J}^T\boldsymbol{e} H=JTJ,b=JTe。这样我们获得了高斯牛顿法的求解过程:

  • 计算残差矩阵关于状态值的雅可比矩阵 J \boldsymbol{J} J
  • 利用雅可比矩阵和残差构建信息矩阵和信息向量 H , b \boldsymbol{H}, \boldsymbol{b} H,b
  • 计算当次迭代量: Δ x = H − 1 b \Delta\boldsymbol{x} = \boldsymbol{H}^{-1}\boldsymbol{b} Δx=H1b
  • 如果迭代量足够小则结束迭代,否则进入下一次迭代

7.莱文贝格马夸特(LM)法

LM 法是在高斯牛顿法的基础上按照阻尼法的思路加入阻尼因子,即求解以下方程:

( H + μ I ) Δ x = b (\boldsymbol{H} + \mu\boldsymbol{I})\Delta\boldsymbol{x} = \boldsymbol{b} (H+μI)Δx=b
上式中,阻尼因子的作用有:

  • 添加进 H \boldsymbol{H} H 保证 H \boldsymbol{H} H 是正定的

  • μ \mu μ 很大时, Δ x = − ( H + μ I ) − 1 b ≈ − 1 μ b = − 1 μ J T E ( x ) \Delta\boldsymbol{x} = -(\boldsymbol{H}+\mu\boldsymbol{I})^{-1}\boldsymbol{b}\approx-\frac{1}{\mu}\boldsymbol{b}=-\frac{1}{\mu}\boldsymbol{J}^T\boldsymbol{E}(\boldsymbol{x}) Δx=(H+μI)1bμ1b=μ1JTE(x),接近最速下降法

  • μ \mu μ 很小时,则接近高斯牛顿法
    因此,合理的设置阻尼因子,能够达到动态对迭代速度进行调节。阻尼因子的设置分为两部分:

  • 初始值的选取

  • 随迭代量变化的更新策略

先来看初始值选取方法,阻尼因子的大小应该根据 J T J \boldsymbol{J}^T\boldsymbol{J} JTJ 的大小来选取,对 J T J \boldsymbol{J}^T\boldsymbol{J} JTJ 进行特征值分解,有: J T J = V Λ V T \boldsymbol{J}^T\boldsymbol{J} = \boldsymbol{V}\boldsymbol{\Lambda}\boldsymbol{V}^T JTJ=VΛVT Λ = diag ( λ 1 , λ 2 , . . . , λ n ) , V = [ v 1 , . . . , v n ] \boldsymbol{\Lambda} = \text{diag}(\lambda_1, \lambda_2,..., \lambda_n), \boldsymbol{V} = [\boldsymbol{v}_1, ..., \boldsymbol{v}_n] Λ=diag(λ1,λ2,...,λn),V=[v1,...,vn],因此,迭代公式化简为:

( V Λ V T + μ I ) Δ x = b Δ x = ( V Λ V T + μ I ) − 1 b = − ∑ i v i T b λ i + μ v i \begin{aligned} (\boldsymbol{V\Lambda}\boldsymbol{V}^T + \mu\boldsymbol{I})\Delta\boldsymbol{x} &= \boldsymbol{b}\\ \Delta\boldsymbol{x} &= (\boldsymbol{V\Lambda}\boldsymbol{V}^T + \mu\boldsymbol{I})^{-1}\boldsymbol{b}\\ &= -\sum_i\frac{\boldsymbol{v}_i^T\boldsymbol{b}}{\lambda_i + \mu}\boldsymbol{v}_i \end{aligned} (VΛVT+μI)ΔxΔx=b=(VΛVT+μI)1b=iλi+μviTbvi
因此,将 μ \mu μ 选为 λ i \lambda_i λi 接近即可,一个简单的思路是设 μ 0 = τ max ⁡ ( J T J ) i i \mu_0 = \tau\max{(\boldsymbol{J}^T\boldsymbol{J})_{ii}} μ0=τmax(JTJ)ii,实际中一般设 τ ≈ [ 1 0 − 8 , 1 ] \tau \approx [10^{-8}, 1] τ[108,1]

接下来看 μ \mu μ 的更新策略,先定性分析应该怎么更新阻尼因子:

  • Δ x \Delta\boldsymbol{x} Δx F ( x ) F(\boldsymbol{x}) F(x) 增加时,应该提高 μ \mu μ 来降低 Δ x \Delta\boldsymbol{x} Δx 即通过减少步长来降低本次迭代带来的影响
  • Δ x \Delta\boldsymbol{x} Δx F ( x ) F(\boldsymbol{x}) F(x) 减少时,应该降低 μ \mu μ 来提高 Δ x \Delta\boldsymbol{x} Δx 即通过增加步长来提高这次迭代的影响

下面进行定量分析,令 L ( Δ x ) = F ( x ) + 2 E ( x ) T J Δ x + 1 2 Δ x T J T J Δ x L(\Delta\boldsymbol{x}) = F(\boldsymbol{x}) +2\boldsymbol{E}(\boldsymbol{x})^T\boldsymbol{J}\Delta\boldsymbol{x} + \frac{1}{2}\Delta\boldsymbol{x}^T\boldsymbol{J}^T\boldsymbol{J}\Delta\boldsymbol{x} L(Δx)=F(x)+2E(x)TJΔx+21ΔxTJTJΔx考虑以下比例因子:

ρ = F ( x ) − F ( x + Δ x ) L ( 0 ) − F ( Δ x ) \rho = \frac{F(\boldsymbol{x}) - F(\boldsymbol{x} + \Delta\boldsymbol{x})}{L(\boldsymbol{0}) - F(\Delta\boldsymbol{x})} ρ=L(0)F(Δx)F(x)F(x+Δx)
Marquardt 提出一个策略:

  • ρ < 0 \rho < 0 ρ<0,表示当前 Δ x \Delta\boldsymbol{x} Δx 使 F ( x ) F(\boldsymbol{x}) F(x) 增大,表示离最优值还很远,应该提高 μ \mu μ 使其接近最速下降法进行较大幅度的更新
  • ρ > 0 \rho > 0 ρ>0 且比较大,表示当前 Δ x \Delta\boldsymbol{x} Δx 使 F ( x ) F(\boldsymbol{x}) F(x) 减小,应该降低 μ \mu μ 使其接近高斯牛顿法,降低速度更新至最优点
  • 如果 ρ > 0 \rho > 0 ρ>0 但比较小,表示已经到最优点附近,则增大阻尼 μ \mu μ,缩小迭代步长

Marquardt 的具体策略如下:

if rho < 0.25: mu = mu * 2
else if rho > 0.75: mu = mu /3

一个使用 Marquardt 策略进行更新的过程如下:
在这里插入图片描述

可以发现,效果并不算很好,随着迭代次数的增加, μ \mu μ 开始进行震荡,表示迭代量周期性的使 F ( x ) F(\boldsymbol{x}) F(x) 增加又下降。

因此,Nielsen 提出了另一个策略,也是 G2O 和 Ceres 中使用的策略:

if rho > 0:mu = mu * max(1/3, 1 - (2 * rho - 1)^3)v = 2
else:mu = mu * vv = 2 * v

使用这种策略进行优化的例子如下:

可以看到, μ \mu μ 随着迭代的进行可以比较平滑的持续下降直至达到收敛。
在这里插入图片描述

8.鲁棒核函数

在进行最小二乘问题中,我们会遇到一些异常观测值使得观测残差特别大,如果不对这些异常点做处理会影响在优化过程中,优化器会尝试最小化异常的残差项,最后影响状态估计的精度,鲁棒核函数就是用来降低这些异常观测值造成的影响。

将鲁棒核函数直接作用在每个残差项上,将最小二乘问题变成如下形式:

F ( x ) = ∑ i ρ ( ∣ ∣ e i ( x ) ∣ ∣ 2 ) F(\boldsymbol{x}) = \sum_i\rho(||e_i(\boldsymbol{x})||^2) F(x)=iρ(∣∣ei(x)2)
使用鲁棒核函数时求解非线性最小二乘的过程
在这个形式下对带有鲁棒核函数的残差进行二阶泰勒展开:

ρ ( s + Δ s ) = ρ ( x ) + ρ ′ ( x ) Δ s + 1 2 ρ ′ ′ ( x ) Δ 2 s \rho(s + \Delta s) = \rho(x) + \rho'(x)\Delta s + \frac{1}{2}\rho''(x)\Delta^2s ρ(s+Δs)=ρ(x)+ρ(x)Δs+21ρ′′(x)Δ2s
上式中,变化量 Δ s \Delta s Δs 的计算方式如下:

Δ s k = ∣ ∣ e i ( x + Δ x ) ∣ ∣ 2 − ∣ ∣ e i ( x ) ∣ ∣ 2 = ∣ ∣ e i ( x ) + J i Δ x ∣ ∣ 2 − ∣ ∣ e i ( x ) ∣ ∣ 2 = 2 e i ( x ) T J i Δ x + Δ x T J i T J i Δ x \begin{aligned} \Delta s_k &= ||e_i(\boldsymbol{x}+\Delta\boldsymbol{x})||^2 - ||e_i(\boldsymbol{x})||^2\\ &= ||e_i(\boldsymbol{x})+\boldsymbol{J}_i\Delta\boldsymbol{x}||^2 - ||e_i(\boldsymbol{x})||^2\\ &= 2e_i(\boldsymbol{x})^T\boldsymbol{J}_i\Delta\boldsymbol{x}+\Delta\boldsymbol{x}^T\boldsymbol{J}_i^T\boldsymbol{J}_i\Delta\boldsymbol{x} \end{aligned} Δsk=∣∣ei(x+Δx)2∣∣ei(x)2=∣∣ei(x)+JiΔx2∣∣ei(x)2=2ei(x)TJiΔx+ΔxTJiTJiΔx
Δ s \Delta s Δs 代入 ρ ( s + Δ s ) \rho(s + \Delta s) ρ(s+Δs) 可得:

ρ ( s + Δ s ) = ρ ( s ) + ρ ′ ( s ) ( 2 e i ( x ) T J i Δ x + Δ x T J i T J i Δ x ) + 1 2 ρ ′ ′ ( s ) ( 2 e i ( x ) T J i Δ x + Δ x T J i T J i Δ x ) 2 ≈ ρ ( s ) + 2 ρ ′ ( s ) e i ( x ) T J i Δ x + ρ ′ ( s ) Δ x T J i T J i Δ x + 2 ρ ′ ′ ( s ) Δ x T J i T e i ( x ) e i ( x ) T J i Δ x \begin{aligned} \rho(s + \Delta s) =& \rho(s) + \rho'(s)(2e_i(\boldsymbol{x})^T\boldsymbol{J}_i\Delta\boldsymbol{x}+\Delta\boldsymbol{x}^T\boldsymbol{J}_i^T\boldsymbol{J}_i\Delta\boldsymbol{x}) \\ &+ \frac{1}{2}\rho''(s)(2e_i(\boldsymbol{x})^T\boldsymbol{J}_i\Delta\boldsymbol{x}+\Delta\boldsymbol{x}^T\boldsymbol{J}_i^T\boldsymbol{J}_i\Delta\boldsymbol{x})^2\\ \approx& \rho(s) + 2\rho'(s)e_i(\boldsymbol{x})^T\boldsymbol{J}_i\Delta\boldsymbol{x}+\rho'(s)\Delta\boldsymbol{x}^T\boldsymbol{J}_i^T\boldsymbol{J}_i\Delta\boldsymbol{x} \\ &+ 2\rho''(s)\Delta\boldsymbol{x}^T\boldsymbol{J}_i^Te_i(\boldsymbol{x})e_i(\boldsymbol{x})^T\boldsymbol{J}_i\Delta\boldsymbol{x} \end{aligned} ρ(s+Δs)=ρ(s)+ρ(s)(2ei(x)TJiΔx+ΔxTJiTJiΔx)+21ρ′′(s)(2ei(x)TJiΔx+ΔxTJiTJiΔx)2ρ(s)+2ρ(s)ei(x)TJiΔx+ρ(s)ΔxTJiTJiΔx+2ρ′′(s)ΔxTJiTei(x)ei(x)TJiΔx
按照之前的思路,对上式求导并使其为 0,可得:

∑ i J i T ( ρ ′ ( s ) + 2 ρ ′ ′ ( s ) e i ( x ) e i ( x ) T ) J Δ x = − ∑ i ρ ′ ( s ) J i T e i ( x ) \sum_i\boldsymbol{J}_i^T(\rho'(s) + 2\rho''(s)e_i(\boldsymbol{\boldsymbol{x}})e_i(\boldsymbol{x})^T)\boldsymbol{J}\Delta\boldsymbol{x} = -\sum_i\rho'(s)\boldsymbol{J}_i^Te_i(\boldsymbol{x}) iJiT(ρ(s)+2ρ′′(s)ei(x)ei(x)T)JΔx=iρ(s)JiTei(x)
对比之前的矩阵 J T J Δ x = − J i T e ( x ) \boldsymbol{J}^T\boldsymbol{J}\Delta\boldsymbol{x} = -\boldsymbol{J}_i^T\boldsymbol{e}(\boldsymbol{x}) JTJΔx=JiTe(x) 可得,当我们使用了鲁棒核函数之后,只需要将各项残差的核函数的一阶二阶导数值计算出,再按照以上形式对信息矩阵和信息向量进行更新即可。

常用的鲁棒核函数
柯西鲁棒核函数:

ρ ( s ) = c 2 log ⁡ ( 1 + s c 2 ) ρ ′ ( s ) = 1 1 + s c 2 ρ ′ ′ ( s ) = − 1 c 2 ( ρ ′ ( s ) ) 2 \begin{aligned} \rho(s) &= c^2\log{(1+\frac{s}{c^2})}\\ \rho'(s) &= \frac{1}{1+\frac{s}{c^2}}\\ \rho''(s) &= -\frac{1}{c^2}(\rho'(s))^2 \end{aligned} ρ(s)ρ(s)ρ′′(s)=c2log(1+c2s)=1+c2s1=c21(ρ(s))2
其中, c c c 为控制参数。当残差是正态分布,Huber c 选为 1.345,Cauchy c 选为 2.3849。不同鲁棒核函数的效果如下图所示:
在这里插入图片描述

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/80972.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

vue结合three.js加载3D模型报404错误

使用vue结合three.js加载3D模型时报404的错误&#xff0c;加载字体库也会报404错误&#xff0c;同样的方法。 vue项目虽然使用npm install three安装了three&#xff0c;但是有些静态资源时读取不到的&#xff0c;当出现异常的404错误时&#xff0c;比如加载3D模型资源时&…

挑战Open AI!!!马斯克宣布成立xAI.

北京时间7月13日凌晨&#xff0c;马斯克在Twitter上宣布&#xff1a;“xAI正式成立&#xff0c;去了解现实。”马斯克表示&#xff0c;推出xAI的原因是想要“了解宇宙的真实本质”。Ghat GPT横空出世已有半年&#xff0c;国内外“百模大战”愈演愈烈&#xff0c;AI大模型的现状…

Oracle 聚合拼接的常用方式

Oracle常用函数&#xff1a;Oracle Database SQL Language Reference, 12c Release 2 (12.2) 1 listagg LISTAGG Syntax Description of the illustration listagg.eps (listagg_overflow_clause::, order_by_clause::, query_partition_clause::) listagg_overflow_claus…

【Docker】Docker的应用场景,Docker 的优点,Ubuntu Docker 安装,使用 Shell 脚本进行安装

作者简介&#xff1a; 辭七七&#xff0c;目前大一&#xff0c;正在学习C/C&#xff0c;Java&#xff0c;Python等 作者主页&#xff1a; 七七的个人主页 文章收录专栏&#xff1a; 七七的闲谈 欢迎大家点赞 &#x1f44d; 收藏 ⭐ 加关注哦&#xff01;&#x1f496;&#x1f…

uniapp根据高度表格合并

没有发现比较友好的能够合并表格单元格插件就自己简单写了一个,暂时格式比较固定 一、效果如下 二、UI视图+逻辑代码 <template><view><uni-card :is-shadow="false" is-full

Java课题笔记~ 使用 Spring 的事务注解管理事务(掌握)

通过Transactional 注解方式&#xff0c;可将事务织入到相应 public 方法中&#xff0c;实现事务管理。 Transactional 的所有可选属性如下所示&#xff1a; propagation&#xff1a;用于设置事务传播属性。该属性类型为 Propagation 枚举&#xff0c; 默认值为 Propagation.R…

Vue2源码分析-day1

初始化数据 vue中最核心的我们都知道那就是响应式数据&#xff0c;数据的变化视图自动更新。那么我们来new一个我们自己的vue 在index.html文件下加入如下代码&#xff0c;这也是vue最常见的基本结构。data已经有了下面我们来获取data的数据 <script src"./vue.js&qu…

IMV7.0

一、背景 经历了多个版本&#xff0c;基础内容在前面&#xff0c;可以使用之前的基础环境&#xff1a; v1&#xff1a; https://blog.csdn.net/wtt234/article/details/132139454 v2&#xff1a; https://blog.csdn.net/wtt234/article/details/132144907 v3&#xff1a; https…

centos7 部署Tomcat和jpress应用

目录 一、静态、动态、伪静态 二、Web 1.0 和 Web 2.0 三、centos7 部署Tomcat 3.1 安装、配置jdk 3.2 安装 Tomcat 3.3 配置服务启动脚本 3.3.1 创建用户和组 3.3.2 创建tomcat.conf文件 3.3.3 创建服务脚本(tomcat.service) 3.3.4 重新加载守护进程并且测试 四、部…

试图将更改推送到 GitHub,但是远程仓库已经包含了您本地没有的工作(可能是其他人提交的修改)

这通常是由于其他人或其他仓库推送到了相同的分支上&#xff0c;导致您的本地仓库和远程仓库之间存在冲突。 错误信息&#xff1a; To github.com:8upersaiyan/CKmuduo.git ! [rejected] main -> main (fetch first) error: failed to push some refs to github.com:8upers…

【干货】商城系统的重要功能特性介绍

电子商务的快速发展&#xff0c;商城系统成为了企业开展线上销售的重要工具。一款功能强大、用户友好的商城系统能够有效提升企业的销售业绩&#xff0c;提供良好的购物体验。下面就商城系统的重要功能特性作一些简单介绍&#xff0c;帮助企业选择合适的系统&#xff0c;打造成…

etcd

文章目录 etcd单机安装设置键值对watch操作读取键过往版本的值压缩修订版本lease租约&#xff08;过期机制&#xff09;授予租约撤销租约keepAlive续约获取租约信息 事务基于etcd实现分布式锁原生实现官方 concurrency 包实现 服务注册与发现Go 操作 Etcd 参考 etcd etcd 是一…

Java课题笔记~ Spring事务的程序举例环境搭建

举例&#xff1a;购买商品 trans_sale 项目 本例要实现购买商品&#xff0c;模拟用户下订单&#xff0c;向订单表添加销售记录&#xff0c;从商品表减少库存。 实现步骤&#xff1a; Step0&#xff1a;创建数据库表 创建两个数据库表 sale , goods sale 销售表&#xff1a;…

访问器模式(C++)

定义 表示一个作用于某对象结构中的各元素的操作。使得可以在不改变(稳定)各元素的类的前提下定义(扩展)作用于这些元素的新操作(变化)。 应用场景 在软件构建过程中&#xff0c;由于需求的改变&#xff0c;某些类层次结构中常常需要增加新的行为(方法)&#xff0c;如果直接…

途乐证券:沪指强势拉升涨0.63%,券商等板块走强,传媒板块活跃

31日早盘&#xff0c;两市股指全线走高&#xff0c;沪指一度涨超1%收复3300点&#xff0c;上证50指数盘中涨逾2%&#xff1b;随后涨幅有所收窄&#xff1b;两市成交额显着放大&#xff0c;北向资金净买入超90亿元。 到午间收盘&#xff0c;沪指涨0.63%报3296.58点&#xff0c;深…

SQL分类及通用语法数据类型(超详细版)

一、SQL分类 SQL是结构化查询语言&#xff08;Structured Query Language&#xff09;的缩写。它是一种用于管理和操作关系型数据库系统的标准化语言。SQL分类如下&#xff1a; DDL: 数据定义语言&#xff0c;用来定义数据库对象&#xff08;数据库、表、字段&#xff09;DML:…

信息安全:认证技术原理与应用.

信息安全&#xff1a;认证技术原理与应用. 认证机制是网络安全的基础性保护措施&#xff0c;是实施访问控制的前提&#xff0c;认证是一个实体向另外一个实体证明其所声称的身份的过程。在认证过程中&#xff0c;需要被证实的实体是声称者&#xff0c;负责检查确认声称者的实体…

如何使用Word转PDF转换器在线工具?在线Word转PDF使用方法

Word转PDF转换器在线&#xff0c;是一种方便快捷的工具&#xff0c;可帮助您在不需要下载任何软件的情况下完成此任务。无论您是需要在工作中共享文档&#xff0c;还是将文件以PDF格式保存以确保格式不变&#xff0c;都可以依靠这款在线工具轻松完成转换。那么如何使用Word转PD…

QGIS二次开发三:显示Shapefile

Shapefile 为 OGR 所支持的最重要的数据格式之一&#xff0c;自然可以被 QGIS 加载。那么该如何显示Shapefile呢&#xff1f; 一、先上代码 #include <qgsapplication.h> #include <qgsproviderregistry.h> #include <qgsmapcanvas.h> #include <qgsvec…

【Spring】Bean的作用域和生命周期

目录 一、引入案例来探讨Bean的作用域 二、Bean的作用域 2.1、Bean的6种作用域 2.2、设置Bean的作用域 三、Spring的执行流程 四、Bean的声明周期 1、生命周期演示 一、引入案例来探讨Bean的作用域 首先我们创建一个User类&#xff0c;定义一个用户信息&#xff0c;在定义…