Graph Optimization
1、非线性最小二乘问题
SLAM问题通常可以表述为非线性最小二乘问题
F
(
x
)
=
1
2
∑
(
x
i
,
x
j
)
∈
C
e
i
j
(
x
i
,
x
j
,
z
i
j
)
T
Ω
i
j
e
i
j
(
x
i
,
x
j
,
z
i
j
)
x
∗
=
arg
min
x
F
(
x
)
\begin{align*} F(\mathbf{x}) &=\frac{1}{2}\sum_{(\mathbf{x}_{i},\mathbf{x}_{j})\in{C}}\mathbf{e}_{ij}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{z}_{ij})^{T}\boldsymbol\Omega_{ij}\mathbf{e}_{ij}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{z}_{ij})\tag{1}\\ \mathbf{x}^{\ast} &=\underset{\mathbf{x}}{\arg\min}\,F(\mathbf{x})\tag{2} \end{align*}
F(x)x∗=21(xi,xj)∈C∑eij(xi,xj,zij)TΩijeij(xi,xj,zij)=xargminF(x)(1)(2)
其中
x
=
[
x
1
T
,
x
2
T
,
⋯
,
x
n
T
]
T
∈
R
N
\mathbf{x}=[\mathbf{x}_{1}^{T},\mathbf{x}_{2}^{T},\cdots,\mathbf{x}_{n}^{T}]^{T}\in\mathbb{R}^{N}
x=[x1T,x2T,⋯,xnT]T∈RN为全部参数组成的向量,
x
k
,
k
=
1
,
⋯
,
n
\mathbf{x}_{k},k=1,\cdots,n
xk,k=1,⋯,n表示一个参数块。
C
C
C是全部参与求和的参数块组合,
e
i
j
(
x
i
,
x
j
,
z
i
j
)
∈
R
M
i
j
\mathbf{e}_{ij}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{z}_{ij})\in\mathbb{R}^{M_{ij}}
eij(xi,xj,zij)∈RMij称为误差函数,
Ω
i
j
∈
R
M
i
j
×
M
i
j
\boldsymbol\Omega_{ij}\in\mathbb{R}^{M_{ij}\times{M_{ij}}}
Ωij∈RMij×Mij称为信息矩阵,为对称正定矩阵。注意,这里为方便推导,假定误差函数的自变量是两个参数块,但在实际情况下,误差函数自变量可以仅包含一个参数块,也可以包含更多的参数块,此时该问题的求解方式可以很方便地通过拓展下面的推导得到。
非线性最小二乘问题可以通过图来构建,其中顶点表示参数快,边表示误差函数,此时称之为图优化问题
2、求解非线性最小二乘问题
略写误差函数中的观测量,并将误差函数视作全部参数的函数,即采用下面的记法
e
i
j
(
x
i
,
x
j
,
z
i
j
)
=
e
i
j
(
x
i
,
x
j
)
=
e
i
j
(
x
)
\mathbf{e}_{ij}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{z}_{ij})=\mathbf{e}_{ij}(\mathbf{x}_{i},\mathbf{x}_{j})=\mathbf{e}_{ij}(\mathbf{x})
eij(xi,xj,zij)=eij(xi,xj)=eij(x)
通过泰勒展开得到误差函数的一阶近似
e
i
j
(
x
i
+
Δ
x
i
,
x
j
+
Δ
x
j
)
=
e
i
j
(
x
+
Δ
x
)
≈
e
i
j
+
J
i
j
Δ
x
(3)
\mathbf{e}_{ij}(\mathbf{x}_{i}+\Delta\mathbf{x}_{i},\mathbf{x}_{j}+\Delta\mathbf{x}_{j}) =\mathbf{e}_{ij}(\mathbf{x}+\Delta\mathbf{x}) \approx\mathbf{e}_{ij}+\mathbf{J}_{ij}\Delta\mathbf{x}\tag{3}
eij(xi+Δxi,xj+Δxj)=eij(x+Δx)≈eij+JijΔx(3)
其中
e
i
j
=
e
i
j
(
x
)
\mathbf{e}_{ij}=\mathbf{e}_{ij}(\mathbf{x})
eij=eij(x),
J
i
j
=
lim
Δ
x
→
0
e
i
j
(
x
+
Δ
x
)
−
e
i
j
(
x
)
Δ
x
∈
R
M
i
j
×
N
\mathbf{J}_{ij}=\lim_{\Delta{x}\rightarrow\mathbf{0}}\frac{\mathbf{e}_{ij}(\mathbf{x}+\Delta\mathbf{x})-\mathbf{e}_{ij}(\mathbf{x})}{\Delta\mathbf{x}}\in\mathbb{R}^{M_{ij}\times{N}}
Jij=limΔx→0Δxeij(x+Δx)−eij(x)∈RMij×N
考虑到参数
x
\mathbf{x}
x可能位于非欧式空间
D
o
m
(
x
)
\mathrm{Dom}(\mathbf{x})
Dom(x),而摄动量
Δ
x
\Delta\mathbf{x}
Δx位于欧式空间
R
N
\mathbb{R}^{N}
RN,故
(
7
)
(7)
(7)中的向量加法很可能导致
x
+
Δ
x
∉
D
o
m
(
x
)
\mathbf{x}+\Delta\mathbf{x}\notin\mathrm{Dom}(\mathbf{x})
x+Δx∈/Dom(x)。为解决该问题,可采用广义加法
⊕
:
D
o
m
(
x
)
×
R
N
→
D
o
m
(
x
)
\oplus:\mathrm{Dom}(\mathbf{x})\times\mathbb{R}^{N}\rightarrow\mathrm{Dom}(\mathbf{x})
⊕:Dom(x)×RN→Dom(x)将
(
7
)
(7)
(7)改写为
e
i
j
(
x
i
⊕
Δ
x
i
,
x
j
⊕
Δ
x
j
)
=
e
i
j
(
x
⊕
Δ
x
)
≈
e
i
j
+
J
i
j
Δ
x
\mathbf{e}_{ij}(\mathbf{x}_{i}\oplus\Delta\mathbf{x}_{i},\mathbf{x}_{j}\oplus\Delta\mathbf{x}_{j}) =\mathbf{e}_{ij}(\mathbf{x}\oplus\Delta\mathbf{x}) \approx\mathbf{e}_{ij}+\mathbf{J}_{ij}\Delta\mathbf{x}
eij(xi⊕Δxi,xj⊕Δxj)=eij(x⊕Δx)≈eij+JijΔx
其中
J
i
j
=
lim
Δ
x
→
0
e
i
j
(
x
⊕
Δ
x
)
−
e
i
j
(
x
)
Δ
x
∈
R
M
i
j
×
N
\mathbf{J}_{ij}=\lim_{\Delta{x}\rightarrow\mathbf{0}}\frac{\mathbf{e}_{ij}(\mathbf{x}\oplus\Delta\mathbf{x})-\mathbf{e}_{ij}(\mathbf{x})}{\Delta\mathbf{x}}\in\mathbb{R}^{M_{ij}\times{N}}
Jij=limΔx→0Δxeij(x⊕Δx)−eij(x)∈RMij×N
记
F
i
j
(
x
i
,
x
j
,
z
i
j
)
=
1
2
e
i
j
(
x
i
,
x
j
,
z
i
j
)
T
Ω
i
j
e
i
j
(
x
i
,
x
j
,
z
i
j
)
F_{ij}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{z}_{ij})=\frac{1}{2}\mathbf{e}_{ij}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{z}_{ij})^{T}\boldsymbol\Omega_{ij}\mathbf{e}_{ij}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{z}_{ij})
Fij(xi,xj,zij)=21eij(xi,xj,zij)TΩijeij(xi,xj,zij),同样略去误差观测量,并将
F
i
j
F_{ij}
Fij视作全部参数的函数,即采用下面的记法
F
i
j
(
x
i
,
x
j
,
z
i
j
)
=
F
i
j
(
x
i
,
x
j
)
=
F
i
j
(
x
)
F_{ij}(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{z}_{ij})=F_{ij}(\mathbf{x}_{i},\mathbf{x}_{j})=F_{ij}(\mathbf{x})
Fij(xi,xj,zij)=Fij(xi,xj)=Fij(x)
计算
F
i
j
(
x
+
Δ
x
)
F_{ij}(\mathbf{x}+\Delta\mathbf{x})
Fij(x+Δx)
F
i
j
(
x
⊕
Δ
x
)
=
1
2
e
i
j
(
x
⊕
Δ
x
)
T
Ω
i
j
e
i
j
(
x
⊕
Δ
x
)
=
1
2
(
e
i
j
+
J
i
j
Δ
x
)
T
Ω
i
j
(
e
i
j
+
J
i
j
Δ
x
)
=
1
2
(
e
i
j
T
Ω
i
j
e
i
j
⏟
2
c
i
j
+
2
e
i
j
T
Ω
i
j
J
i
j
⏟
b
i
j
T
Δ
x
+
Δ
x
T
J
i
j
T
Ω
i
j
J
i
j
⏟
H
i
j
Δ
x
)
=
c
i
j
+
b
i
j
T
Δ
x
+
1
2
Δ
x
T
H
i
j
Δ
x
\begin{align*} F_{ij}(\mathbf{x}\oplus\Delta\mathbf{x}) &=\frac{1}{2}\mathbf{e}_{ij}(\mathbf{x}\oplus\Delta\mathbf{x})^{T}\boldsymbol\Omega_{ij}\mathbf{e}_{ij}(\mathbf{x}\oplus\Delta\mathbf{x})\\ &=\frac{1}{2}(\mathbf{e}_{ij}+\mathbf{J}_{ij}\Delta\mathbf{x})^{T}\boldsymbol\Omega_{ij}(\mathbf{e}_{ij}+\mathbf{J}_{ij}\Delta\mathbf{x})\\ &=\frac{1}{2}\left(\underbrace{\mathbf{e}_{ij}^{T}\boldsymbol\Omega_{ij}\mathbf{e}_{ij}}_{2c_{ij}}+2\underbrace{\mathbf{e}_{ij}^{T}\boldsymbol\Omega_{ij}\mathbf{J}_{ij}}_{\mathbf{b}_{ij}^{T}}\Delta\mathbf{x}+\Delta\mathbf{x}^{T}\underbrace{\mathbf{J}_{ij}^{T}\boldsymbol\Omega_{ij}\mathbf{J}_{ij}}_{\mathbf{H}_{ij}}\Delta\mathbf{x}\right)\\ &=c_{ij}+\mathbf{b}_{ij}^{T}\Delta\mathbf{x}+\frac{1}{2}\Delta\mathbf{x}^{T}\mathbf{H}_{ij}\Delta\mathbf{x}\tag{4} \end{align*}
Fij(x⊕Δx)=21eij(x⊕Δx)TΩijeij(x⊕Δx)=21(eij+JijΔx)TΩij(eij+JijΔx)=21
2cij
eijTΩijeij+2bijT
eijTΩijJijΔx+ΔxTHij
JijTΩijJijΔx
=cij+bijTΔx+21ΔxTHijΔx(4)
将
(
4
)
(4)
(4)代入
(
1
)
(1)
(1)计算
F
(
x
+
Δ
x
)
F(\mathbf{x}+\Delta\mathbf{x})
F(x+Δx)
F
(
x
⊕
Δ
x
)
=
∑
(
x
i
,
x
j
)
∈
C
F
i
j
(
x
⊕
Δ
x
)
=
∑
(
x
i
,
x
j
)
∈
C
c
i
j
+
b
i
j
Δ
x
+
1
2
Δ
x
T
H
i
j
Δ
x
=
c
+
b
T
Δ
x
+
1
2
Δ
x
T
H
Δ
x
\begin{align*} F(\mathbf{x}\oplus\Delta\mathbf{x}) &=\sum_{(\mathbf{x}_{i},\mathbf{x}_{j})\in\mathcal{C}}F_{ij}(\mathbf{x}\oplus\Delta\mathbf{x})\\ &=\sum_{(\mathbf{x}_{i},\mathbf{x}_{j})\in\mathcal{C}}c_{ij}+\mathbf{b}_{ij}\Delta\mathbf{x}+\frac{1}{2}\Delta\mathbf{x}^{T}\mathbf{H}_{ij}\Delta\mathbf{x}\\ &=c+\mathbf{b}^{T}\Delta\mathbf{x}+\frac{1}{2}\Delta\mathbf{x}^{T}\mathbf{H}\Delta\mathbf{x}\tag{5} \end{align*}
F(x⊕Δx)=(xi,xj)∈C∑Fij(x⊕Δx)=(xi,xj)∈C∑cij+bijΔx+21ΔxTHijΔx=c+bTΔx+21ΔxTHΔx(5)
其中 b = ∑ ( x i , x j ) ∈ C b i j ∈ R n \mathbf{b}=\sum_{(\mathbf{x}_{i},\mathbf{x}_{j})\in\mathcal{C}}\mathbf{b}_{ij}\in\mathbb{R}^{n} b=∑(xi,xj)∈Cbij∈Rn, H = ∑ ( x i , x j ) ∈ C H i j ∈ R N × N \mathbf{H}=\sum_{(\mathbf{x}_{i},\mathbf{x}_{j})\in\mathcal{C}}\mathbf{H}_{ij}\in\mathbb{R}^{N\times{N}} H=∑(xi,xj)∈CHij∈RN×N
实际上,
b
\mathbf{b}
b等于
F
F
F在
x
\mathbf{x}
x点处的梯度,而采用
H
\mathbf{H}
H近似
F
F
F在
x
\mathbf{x}
x点处的海森矩阵
b
=
∑
(
x
i
,
x
j
)
∈
C
[
∇
e
i
j
(
x
)
Ω
i
j
e
i
j
(
x
)
]
=
∇
F
(
x
)
H
=
∑
(
x
i
,
x
j
)
∈
C
[
∇
e
i
j
(
x
)
Ω
i
j
∇
e
i
j
(
x
)
T
]
→
∇
2
F
(
x
)
=
∑
(
x
i
,
x
j
)
∈
C
[
⋯
+
∇
e
i
j
(
x
)
Ω
i
j
∇
e
i
j
(
x
)
T
]
\begin{align*} \mathbf{b}&=\sum_{(\mathbf{x}_{i},\mathbf{x}_{j})\in\mathcal{C}}\left[\nabla\mathbf{e}_{ij}(\mathbf{x})\boldsymbol{\Omega}_{ij}\mathbf{e}_{ij}(\mathbf{x})\right]=\nabla{F}(\mathbf{x})\\ \mathbf{H}&=\sum_{(\mathbf{x}_{i},\mathbf{x}_{j})\in\mathcal{C}}\left[\nabla\mathbf{e}_{ij}(\mathbf{x})\boldsymbol{\Omega}_{ij}\nabla\mathbf{e}_{ij}(\mathbf{x})^{T}\right]\rightarrow\nabla^{2}F(\mathbf{x})=\sum_{(\mathbf{x}_{i},\mathbf{x}_{j})\in\mathcal{C}}\left[\cdots+\nabla\mathbf{e}_{ij}(\mathbf{x})\boldsymbol{\Omega}_{ij}\nabla\mathbf{e}_{ij}(\mathbf{x})^{T}\right] \end{align*}
bH=(xi,xj)∈C∑[∇eij(x)Ωijeij(x)]=∇F(x)=(xi,xj)∈C∑[∇eij(x)Ωij∇eij(x)T]→∇2F(x)=(xi,xj)∈C∑[⋯+∇eij(x)Ωij∇eij(x)T]
2.1、梯度下降法
梯度下降法通过取增量方向为目标函数的负梯度方向来保证目标函数的一阶近似下降
保留
(
5
)
(5)
(5)至一阶项
F
(
x
⊕
Δ
x
)
=
c
+
b
T
Δ
x
F(\mathbf{x}\oplus\Delta\mathbf{x})=c+\mathbf{b}^{T}\Delta\mathbf{x}
F(x⊕Δx)=c+bTΔx
此时取
Δ
x
=
−
λ
b
(6)
\Delta\mathbf{x}=-\lambda\mathbf{b}\tag{6}
Δx=−λb(6)
即可保证目标函数的一阶近似下降,其中
λ
>
0
\lambda>0
λ>0称为步长
梯度下降法求解非线性最小二乘问题的具体步骤如下:
- 令 k = 0 k=0 k=0,给定初始值 x 0 \mathbf{x}_{0} x0, λ \lambda λ
- 若 k k k达到最大迭代次数,则停止迭代;否则根据 x k \mathbf{x}_{k} xk求出当前的和 b k \mathbf{b}_{k} bk
- 令 Δ x k = − λ b k \Delta\mathbf{x}_{k}=-\lambda\mathbf{b}_{k} Δxk=−λbk
- 如果 Δ x k \Delta\mathbf{x}_{k} Δxk足够小,则停止迭代;否则令 x k + 1 = x k ⊕ Δ x k \mathbf{x}_{k+1}=\mathbf{x}_{k}\oplus\Delta\mathbf{x}_{k} xk+1=xk⊕Δxk, k = k + 1 k=k+1 k=k+1,返回第2步
2.2、高斯牛顿法
高斯牛顿法采用 H \mathbf{H} H近似 F ( x ) F(\mathbf{x}) F(x)在 x \mathbf{x} x点处的海森矩阵,避免了海森矩阵的计算,又可以通过求解二阶近似的最小值实现更快的收敛
因为
H
\mathbf{H}
H矩阵半正定,故根据
(
8
)
(8)
(8),
Δ
x
\Delta\mathbf{x}
Δx的函数
F
(
x
⊕
Δ
x
)
F(\mathbf{x}\oplus\Delta\mathbf{x})
F(x⊕Δx)为凸函数,因此其最小值在
∂
F
(
x
⊕
Δ
x
)
∂
Δ
x
=
0
\frac{\partial{F(\mathbf{x}\oplus\Delta\mathbf{x})}}{\partial\Delta\mathbf{x}}=\mathbf{0}
∂Δx∂F(x⊕Δx)=0处取得,即要求
Δ
x
\Delta\mathbf{x}
Δx满足下式时取得
H
Δ
x
=
−
b
(7)
\mathbf{H}\Delta\mathbf{x}=-\mathbf{b}\tag{7}
HΔx=−b(7)
(
11
)
(11)
(11)称为增量方程,求解增量方程是整个优化问题的核心。
高斯牛顿法求解非线性最小二乘问题的具体步骤如下:
- 令 k = 0 k=0 k=0,给定初始值 x 0 \mathbf{x}_{0} x0,
- 若 k k k达到最大迭代次数,则停止迭代;否则根据 x k \mathbf{x}_{k} xk求出当前的 H k \mathbf{H}_{k} Hk和 b k \mathbf{b}_{k} bk
- 求解增量方程: H k Δ x k = − b k \mathbf{H}_{k}\Delta\mathbf{x}_{k}=-\mathbf{b}_{k} HkΔxk=−bk:当 H \mathbf{H} H正定时,增量方程可以通过Cholesky分解高效求解;当 H \mathbf{H} H非正定时,可取 Δ x k = Δ x k − 1 \Delta\mathbf{x}_{k}=\Delta\mathbf{x}_{k-1} Δxk=Δxk−1
- 如果 Δ x k \Delta\mathbf{x}_{k} Δxk足够小,则停止迭代;否则令 x k + 1 = x k ⊕ Δ x k \mathbf{x}_{k+1}=\mathbf{x}_{k}\oplus\Delta\mathbf{x}_{k} xk+1=xk⊕Δxk, k = k + 1 k=k+1 k=k+1,返回第2步
2.3、Levenberg-Marquardt算法
Levenberg-Marquardt算法,简称LM算法。考虑到高斯牛顿法中对海森矩阵的近似不一定准确,二阶泰勒展开式也只能才展开点附近有较好的效果,为结合梯度下降法和高斯牛顿法的优点,可以将增量方程改写为
(
H
+
λ
I
)
Δ
x
=
−
b
(8)
(\mathbf{H}+\lambda\mathbf{I})\Delta\mathbf{x}=-\mathbf{b}\tag{8}
(H+λI)Δx=−b(8)
此方程称为LM算法的增量方程
此外,通过 ρ = F ( x ⊕ Δ x ) − F ( x ) ∇ F ( x ) T Δ x \rho=\frac{F(\mathbf{x}\oplus\Delta\mathbf{x})-F(\mathbf{x})}{\nabla{F(\mathbf{x})^{T}\Delta\mathbf{x}}} ρ=∇F(x)TΔxF(x⊕Δx)−F(x)反映一阶近似的准确度,如果一阶近似准确度高,则 λ \lambda λ值增大,更接近梯度下降法,如果一阶近似准确度低,即函数非线性的性质更加突出,则 λ \lambda λ值减小,更接近高斯牛顿法。
LM算法求解非线性最小二乘问题的具体步骤如下:
- 令 k = 0 k=0 k=0,给定初始值 x 0 \mathbf{x}_{0} x0, λ 0 \lambda_{0} λ0
- 若 k k k达到最大迭代次数,则停止迭代;否则根据 x k \mathbf{x}_{k} xk求出当前的 H k , b k , ρ k \mathbf{H}_{k},\mathbf{b}_{k},\rho_{k} Hk,bk,ρk
- 若 ρ k > 3 4 \rho_{k}>\frac{3}{4} ρk>43,则 λ k + 1 = 2 λ k \lambda_{k+1}=2\lambda_{k} λk+1=2λk;若 ρ k < 1 4 \rho_{k}<\frac{1}{4} ρk<41,则 λ k + 1 = 1 2 λ k \lambda_{k+1}=\frac{1}{2}\lambda_{k} λk+1=21λk;若 1 4 ⩽ ρ k ⩽ 3 4 \frac{1}{4}\leqslant\rho_{k}\leqslant\frac{3}{4} 41⩽ρk⩽43,则 λ k + 1 = λ k \lambda_{k+1}=\lambda_{k} λk+1=λk
- 求解 ( H k + λ k + 1 I ) Δ x k = − b k (\mathbf{H}_{k}+\lambda_{k+1}\mathbf{I})\Delta\mathbf{x}_{k}=-\mathbf{b}_{k} (Hk+λk+1I)Δxk=−bk:因为 ( H k + λ k + 1 I ) (\mathbf{H}_{k}+\lambda_{k+1}\mathbf{I}) (Hk+λk+1I)恒为正定,故可以通过Cholesky分解高效求解
- 如果 Δ x k \Delta\mathbf{x}_{k} Δxk足够小,则停止迭代; x k + 1 = x k ⊕ Δ x k \mathbf{x}_{k+1}=\mathbf{x}_{k}\oplus\Delta\mathbf{x}_{k} xk+1=xk⊕Δxk, k = k + 1 k=k+1 k=k+1,返回第2步
3、核函数
注意到,当某条边的误差 e i j \mathbf{e}_{ij} eij很大时, F i j ( x ) = e i j T Ω i j e i j F_{ij}(\mathbf{x})=\mathbf{e}_{ij}^{T}\boldsymbol{\Omega}_{ij}\mathbf{e}_{ij} Fij(x)=eijTΩijeij会很大,其梯度 ∇ F i j ( x ) = J i j T Ω i j e i j \nabla{F}_{ij}(\mathbf{x})=\mathbf{J}_{ij}^{T}\boldsymbol{\Omega}_{ij}\mathbf{e}_{ij} ∇Fij(x)=JijTΩijeij也会很大,而算法会根据梯度更大幅度地调整这条边所连接的节点的估计值,而掩盖这些节点与其他边的关系。
可以定义核函数 ρ ( ⋅ ) \rho(\cdot) ρ(⋅),在非线性最小二乘问题中用 ρ ( F i j ( x ) ) \rho(F_{ij}(\mathbf{x})) ρ(Fij(x))代替 F i j ( x ) F_{ij}(\mathbf{x}) Fij(x)
最常用的核函数是Huber核函数,其定义为
ρ
(
x
)
=
{
x
0
⩽
x
⩽
δ
2
2
δ
(
x
−
δ
2
)
x
>
δ
2
\rho(x)=\left\{ \begin{array}{ll} x&0\leqslant{x}\leqslant\delta^{2}\\ 2\delta(\sqrt{x}-\frac{\delta}{2})&{x}>\delta^{2} \end{array} \right.
ρ(x)={x2δ(x−2δ)0⩽x⩽δ2x>δ2
Huber函数图像如下
![Huber](https://img-blog.csdnimg.cn/direct/532b19246b434def8859841b2503803b.jpeg#pic_center)
Huber核函数是二阶可导的,其一阶导数和二阶导数为
ρ
′
(
x
)
=
{
1
0
⩽
x
⩽
δ
2
δ
x
x
>
δ
2
ρ
′
′
(
x
)
=
{
0
0
⩽
x
⩽
δ
2
−
1
2
δ
(
x
)
3
x
>
δ
2
\begin{align*} \rho^{\prime}(x)=\left\{ \begin{array}{ll} 1&0\leqslant{x}\leqslant\delta^{2}\\ \frac{\delta}{\sqrt{x}}&{x}>\delta^{2} \end{array} \right.\quad\quad\quad \rho^{\prime\prime}(x)=\left\{ \begin{array}{ll} 0&0\leqslant{x}\leqslant\delta^{2}\\ -\frac{1}{2}\frac{\delta}{(\sqrt{x})^{3}}&{x}>\delta^{2} \end{array} \right.\\ \end{align*}
ρ′(x)={1xδ0⩽x⩽δ2x>δ2ρ′′(x)={0−21(x)3δ0⩽x⩽δ2x>δ2
对
ρ
(
F
i
j
(
x
⊕
Δ
x
)
)
\rho(F_{ij}(\mathbf{x}\oplus\Delta\mathbf{x}))
ρ(Fij(x⊕Δx))做如下展开
ρ
(
F
i
j
(
x
⊕
Δ
x
)
)
=
ρ
(
F
i
j
(
x
)
)
+
∇
ρ
(
F
i
j
(
x
)
)
T
Δ
x
+
1
2
Δ
x
T
∇
2
ρ
(
F
i
j
(
x
)
)
Δ
x
\rho(F_{ij}(\mathbf{x}\oplus\Delta\mathbf{x}))=\rho(F_{ij}(\mathbf{x}))+\nabla\rho(F_{ij}(\mathbf{x}))^{T}\Delta\mathbf{x}+\frac{1}{2}\Delta\mathbf{x}^{T}\nabla^{2}\rho(F_{ij}(\mathbf{x}))\Delta\mathbf{x}
ρ(Fij(x⊕Δx))=ρ(Fij(x))+∇ρ(Fij(x))TΔx+21ΔxT∇2ρ(Fij(x))Δx
其中
∇
ρ
(
F
i
j
(
x
)
)
=
∂
ρ
∂
F
i
j
∣
F
i
j
(
x
)
∇
F
i
j
(
x
)
=
ρ
′
(
F
i
j
(
x
)
)
b
i
j
=
J
i
j
T
(
ρ
′
(
F
i
j
(
x
)
)
Ω
i
j
)
e
i
j
∇
2
ρ
(
F
i
j
(
x
)
)
=
∂
2
ρ
∂
F
i
j
2
∣
F
i
j
(
x
)
∇
F
i
j
(
x
)
∇
F
i
j
(
x
)
T
+
∂
ρ
∂
F
i
j
∣
F
i
j
(
x
)
∇
2
F
i
j
(
x
)
=
ρ
′
′
(
F
i
j
(
x
)
)
b
i
j
b
i
j
T
+
ρ
′
(
F
i
j
(
x
)
)
H
i
j
→
ρ
′
(
F
i
j
(
x
)
)
H
i
j
=
J
i
j
T
(
ρ
′
(
F
i
j
(
x
)
)
Ω
i
j
)
J
i
j
\begin{align*} \nabla\rho(F_{ij}(\mathbf{x}))&=\left.\frac{\partial\rho}{\partial{F}_{ij}}\right|_{F_{ij}(\mathbf{x})}\nabla{F}_{ij}(\mathbf{x})=\rho^{\prime}(F_{ij}(\mathbf{x}))\mathbf{b}_{ij}=\mathbf{J}_{ij}^{T}\left(\rho^{\prime}(F_{ij}(\mathbf{x}))\boldsymbol{\Omega}_{ij}\right)\mathbf{e}_{ij}\\ \nabla^{2}\rho(F_{ij}(\mathbf{x}))&=\left.\frac{\partial^{2}\rho}{\partial{F}_{ij}^{2}}\right|_{F_{ij}(\mathbf{x})}\nabla{F}_{ij}(\mathbf{x})\nabla{F}_{ij}(\mathbf{x})^{T}+\left.\frac{\partial\rho}{\partial{F}_{ij}}\right|_{F_{ij}(\mathbf{x})}\nabla^{2}{F}_{ij}(\mathbf{x})=\rho^{\prime\prime}(F_{ij}(\mathbf{x}))\mathbf{b}_{ij}\mathbf{b}_{ij}^{T}+\rho^{\prime}(F_{ij}(\mathbf{x}))\mathbf{H}_{ij}\rightarrow\rho^{\prime}(F_{ij}(\mathbf{x}))\mathbf{H}_{ij}=\mathbf{J}_{ij}^{T}\left(\rho^{\prime}(F_{ij}(\mathbf{x}))\boldsymbol{\Omega}_{ij}\right)\mathbf{J}_{ij} \end{align*}
∇ρ(Fij(x))∇2ρ(Fij(x))=∂Fij∂ρ
Fij(x)∇Fij(x)=ρ′(Fij(x))bij=JijT(ρ′(Fij(x))Ωij)eij=∂Fij2∂2ρ
Fij(x)∇Fij(x)∇Fij(x)T+∂Fij∂ρ
Fij(x)∇2Fij(x)=ρ′′(Fij(x))bijbijT+ρ′(Fij(x))Hij→ρ′(Fij(x))Hij=JijT(ρ′(Fij(x))Ωij)Jij
上式第二行中,因为考虑到可能存在
ρ
′
′
(
F
i
j
(
x
)
)
<
0
\rho^{\prime\prime}(F_{ij}(\mathbf{x}))<0
ρ′′(Fij(x))<0,故为确保
∇
2
ρ
(
F
i
j
(
x
)
)
\nabla^{2}\rho(F_{ij}(\mathbf{x}))
∇2ρ(Fij(x))正定性与
H
i
j
\mathbf{H}_{ij}
Hij一致,舍去含
ρ
′
′
(
F
i
j
(
x
)
)
\rho^{\prime\prime}(F_{ij}(\mathbf{x}))
ρ′′(Fij(x))的项。由上面的分析知,对某条边应用核函数后,在优化过程中只需修改其对应的信息矩阵即可。
令 b i j ← ∇ ρ ( F i j ( x ) ) \mathbf{b}_{ij}\leftarrow\nabla\rho(F_{ij}(\mathbf{x})) bij←∇ρ(Fij(x)), H i j ← ∇ 2 ρ ( F i j ( x ) ) \mathbf{H}_{ij}\leftarrow\nabla^{2}\rho(F_{ij}(\mathbf{x})) Hij←∇2ρ(Fij(x)),从而可以构建新的增量方程。
通过上述过程,可以解决因误差 e i j \mathbf{e}_{ij} eij过大而导致的问题。
附录
§1、标量函数的二阶泰勒展开
若一个多元标量函数
f
:
R
N
→
R
f:\mathbb{R}^{N}\rightarrow\mathbb{R}
f:RN→R在
a
\mathbf{a}
a点处二阶可微,则存在
a
\mathbf{a}
a的邻域
U
(
x
,
δ
)
=
{
x
+
Δ
x
∣
∥
Δ
x
∥
<
δ
}
U(\mathbf{x},\delta)=\{\mathbf{x}+\Delta\mathbf{x}|\|\Delta\mathbf{x}\|<\delta\}
U(x,δ)={x+Δx∣∥Δx∥<δ},使得在该邻域中
f
(
x
+
Δ
x
)
=
f
(
x
)
+
J
f
(
x
)
Δ
x
+
1
2
Δ
x
T
H
f
(
x
)
Δ
x
+
o
(
∥
Δ
x
∥
2
)
(A1)
f(\mathbf{x}+\Delta\mathbf{x})=f(\mathbf{x})+\mathbf{J}f(\mathbf{x})\Delta\mathbf{x}+\frac{1}{2}\Delta\mathbf{x}^{T}\mathbf{H}f(\mathbf{x})\Delta\mathbf{x}+\mathrm{o}(\|\Delta\mathbf{x}\|^{2})\tag{A1}
f(x+Δx)=f(x)+Jf(x)Δx+21ΔxTHf(x)Δx+o(∥Δx∥2)(A1)
其中
J f ( x ) = ∇ f ( x ) T = lim Δ x → 0 f ( x + Δ x ) − f ( x ) Δ x T ∈ R 1 × N \mathbf{J}f(\mathbf{x})=\nabla{f}(\mathbf{x})^{T}=\lim_{\Delta\mathbf{x}\rightarrow\mathbf{0}}\frac{f(\mathbf{x}+\Delta\mathbf{x})-f(\mathbf{x})}{\Delta\mathbf{x}^{T}}\in\mathbb{R}^{1\times N} Jf(x)=∇f(x)T=limΔx→0ΔxTf(x+Δx)−f(x)∈R1×N称为雅可比矩阵;
H f ( x ) = ∇ 2 f ( x ) = lim Δ x → 0 ∇ f ( x + Δ x ) − ∇ f ( x ) Δ x T ∈ R N × N \mathbf{H}f(\mathbf{x})=\nabla^{2}{f}(\mathbf{x})=\lim_{\Delta\mathbf{x}\rightarrow\mathbf{0}}\frac{\nabla{f}(\mathbf{x}+\Delta\mathbf{x})-\nabla{f}(\mathbf{x})}{\Delta\mathbf{x}^{T}}\in\mathbb{R}^{N\times{N}} Hf(x)=∇2f(x)=limΔx→0ΔxT∇f(x+Δx)−∇f(x)∈RN×N称为海森矩阵。
( A 1 ) \mathrm{(A1)} (A1)称为标量函数 f f f的二阶泰勒展开。
下面根据上述理论推导当 f \mathbf{f} f的定义域为非欧式空间 D o m ( x ) \mathrm{Dom}(\mathbf{x}) Dom(x)时的二阶泰勒展开
记
h
(
Δ
x
)
=
f
(
x
⊕
Δ
x
)
:
R
N
→
R
h(\Delta\mathbf{x})=f(\mathbf{x}\oplus\Delta\mathbf{x}):\mathbb{R}^{N}\rightarrow\mathbb{R}
h(Δx)=f(x⊕Δx):RN→R,
⊕
:
D
o
m
(
x
)
×
R
N
→
D
o
m
(
x
)
\oplus:\mathrm{Dom}(\mathbf{x})\times\mathbb{R}^{N}\rightarrow\mathrm{Dom}(\mathbf{x})
⊕:Dom(x)×RN→Dom(x),若
h
h
h在
0
\mathbf{0}
0点处二阶可微,则存在
0
\mathbf{0}
0的邻域
U
(
0
,
δ
)
=
{
Δ
x
∣
∥
Δ
x
∥
<
δ
}
U(\mathbf{0},\delta)=\{\Delta\mathbf{x}|\|\Delta\mathbf{x}\|<\delta\}
U(0,δ)={Δx∣∥Δx∥<δ},使得在该邻域中
h
(
Δ
x
)
=
h
(
0
)
+
J
h
(
0
)
T
Δ
x
+
1
2
Δ
x
T
H
h
(
0
)
Δ
x
+
o
(
∥
Δ
x
∥
2
)
(A2)
h(\Delta\mathbf{x})=h(\mathbf{0})+\mathbf{J}h(\mathbf{0})^{T}\Delta\mathbf{x}+\frac{1}{2}\Delta\mathbf{x}^{T}\mathbf{H}h(\mathbf{0})\Delta\mathbf{x}+\mathrm{o}(\|\Delta\mathbf{x}\|^{2})\tag{A2}
h(Δx)=h(0)+Jh(0)TΔx+21ΔxTHh(0)Δx+o(∥Δx∥2)(A2)
其中
J h ( 0 ) = ∇ h ( 0 ) T = lim Δ x → 0 h ( Δ x ) − h ( 0 ) Δ x T = lim Δ x → 0 f ( x ⊕ Δ x ) − f ( x ) Δ x T ∈ R 1 × N \mathbf{J}h(\mathbf{0})=\nabla{h}(\mathbf{0})^{T}=\lim_{\Delta\mathbf{x}\rightarrow\mathbf{0}}\frac{h(\Delta\mathbf{x})-h(\mathbf{0})}{\Delta\mathbf{x}^{T}}=\lim_{\Delta\mathbf{x}\rightarrow\mathbf{0}}\frac{f(\mathbf{x}\oplus\Delta\mathbf{x})-f(\mathbf{x})}{\Delta\mathbf{x}^{T}}\in\mathbb{R}^{1\times N} Jh(0)=∇h(0)T=limΔx→0ΔxTh(Δx)−h(0)=limΔx→0ΔxTf(x⊕Δx)−f(x)∈R1×N,可以发现 J h ( 0 ) \mathbf{J}h(\mathbf{0}) Jh(0)和 ( A 1 ) \mathrm{(A1)} (A1)中 J f ( x ) \mathbf{J}f(\mathbf{x}) Jf(x)具有类似的结构,故记 J f ( x ) = ∇ f ( x ) T = J h ( 0 ) \mathbf{J}f(\mathbf{x})=\nabla{f}(\mathbf{x})^{T}=\mathbf{J}h(\mathbf{0}) Jf(x)=∇f(x)T=Jh(0),称 J f ( x ) \mathbf{J}f(\mathbf{x}) Jf(x)为对增量 Δ x \Delta\mathbf{x} Δx的雅可比矩阵;
H h ( 0 ) = ∇ 2 h ( 0 ) = lim Δ x → 0 ∇ h ( Δ x ) − ∇ h ( 0 ) Δ x T = lim Δ x → 0 ∇ f ( x ⊕ Δ x ) − ∇ f ( x ) Δ x T ∈ R N × N \mathbf{H}h(\mathbf{0})=\nabla^{2}h(\mathbf{0})=\lim_{\Delta\mathbf{x}\rightarrow\mathbf{0}}\frac{\nabla{h}(\Delta\mathbf{x})-\nabla{h}(\mathbf{0})}{\Delta\mathbf{x}^{T}}=\lim_{\Delta\mathbf{x}\rightarrow\mathbf{0}}\frac{\nabla{f}(\mathbf{x}\oplus\Delta\mathbf{x})-\nabla{f}(\mathbf{x})}{\Delta\mathbf{x}^{T}}\in\mathbb{R}^{N\times{N}} Hh(0)=∇2h(0)=limΔx→0ΔxT∇h(Δx)−∇h(0)=limΔx→0ΔxT∇f(x⊕Δx)−∇f(x)∈RN×N,可以发现 H h ( 0 ) \mathbf{H}h(\mathbf{0}) Hh(0)和 ( A 1 ) \mathrm{(A1)} (A1)中 H f ( x ) \mathbf{H}f(\mathbf{x}) Hf(x)具有类似的结构,故记 H f ( x ) = ∇ 2 f ( x ) = H h ( 0 ) \mathbf{H}f(\mathbf{x})=\nabla^{2}f(\mathbf{x})=\mathbf{H}h(\mathbf{0}) Hf(x)=∇2f(x)=Hh(0),称 H f ( x ) \mathbf{H}f(\mathbf{x}) Hf(x)为对增量 Δ x \Delta\mathbf{x} Δx的海森矩阵。
综上,
(
A
2
)
\mathrm{(A2)}
(A2)也可以写成
f
(
x
⊕
Δ
x
)
=
f
(
x
)
+
J
f
(
x
)
Δ
x
+
1
2
Δ
x
T
H
f
(
x
)
Δ
x
+
o
(
∥
Δ
x
∥
2
)
(A3)
f(\mathbf{x}\oplus\Delta\mathbf{x})=f(\mathbf{x})+\mathbf{J}f(\mathbf{x})\Delta\mathbf{x}+\frac{1}{2}\Delta\mathbf{x}^{T}\mathbf{H}f(\mathbf{x})\Delta\mathbf{x}+\mathrm{o}(\|\Delta\mathbf{x}\|^{2})\tag{A3}
f(x⊕Δx)=f(x)+Jf(x)Δx+21ΔxTHf(x)Δx+o(∥Δx∥2)(A3)
(
A
3
)
\mathrm{(A3)}
(A3)称为
D
o
m
(
x
)
\mathrm{Dom}(\mathbf{x})
Dom(x)上的标量函数
f
\mathbf{f}
f的二阶泰勒展开。
为简便起见,常将上面的雅可比矩阵和海森矩阵分别用 J ( ⋅ ) \mathbf{J}(\cdot) J(⋅), H ( ⋅ ) \mathbf{H}(\cdot) H(⋅)代替