前言
SLAM 方向对 Gaussian Process (GP)的需求不大,但这两年有好几篇 IROS,ICRA 的论文用高斯过程来拟合轨迹,拟合误差模型等,因此这篇笔记对高斯过程概念和原理进行简单梳理,理清楚 GP 是怎么来的,以及怎么用它。如果想更进一步系统学习下,推荐 MIT 出版的 Gaussian Processes for Machine Learning.
高斯过程是什么
高斯过程可以通俗的认为是函数拟合或回归近似,即给出变量 x \boldsymbol{x} x, 能预测出 f ( x ) f(\boldsymbol{x}) f(x)。高斯分布是针对向量 x \boldsymbol{x} x ,高斯过程则是对函数 f ( x ) f(\boldsymbol{x}) f(x)。更学术一点的描述,高斯过程认为是一个服从均值是函数 m ( x ) m(\boldsymbol{x}) m(x),方差也是函数 k ( x , x ′ ) k\left(\boldsymbol{x}, \boldsymbol{x}^{\prime}\right) k(x,x′) 的高斯分布 ,其中对任意变量 x , x ′ \boldsymbol{x}, \boldsymbol{x}^{\prime} x,x′ 有 m ( x ) = E [ f ( x ) ] , k ( x , x ′ ) = Cov ( f ( x ) , f ( x ′ ) ) m(\boldsymbol{x})=\mathbb{E}[f(\boldsymbol{x})] \text { , } k\left(\boldsymbol{x}, \boldsymbol{x}^{\prime}\right)=\operatorname{Cov}\left(f(\boldsymbol{x}), f\left(x^{\prime}\right)\right) m(x)=E[f(x)] , k(x,x′)=Cov(f(x),f(x′))。
举个简单的例子,一个函数 f ( ⋅ ) f(\cdot) f(⋅) 的任意有限子集 f ( x 1 ) , … , f ( x n ) f\left(\boldsymbol{x}_{1}\right), \ldots, f\left(\boldsymbol{x}_{n}\right) f(x1),…,f(xn) 服从多元高斯分布,则称 f ( ⋅ ) f(\cdot) f(⋅) 是一个高斯过程,其中输入 { x n } n = 1 N \left\{\boldsymbol{x}_{n}\right\}_{n=1}^{N} {xn}n=1N 是任意维的向量。比如, { x n } n = 1 365 \left\{\boldsymbol{x}_{n}\right\}_{n=1}^{365} {xn}n=1365 表示 365 天, f ( x n ) f\left(\boldsymbol{x}_{n}\right) f(xn) 表示某天的温度。
对于 N 组 D 维的输入数据
X
∈
R
N
×
D
\boldsymbol{X} \in \mathbb{R}^{N \times D}
X∈RN×D,每一行表示一个输入(如24个小时时间点),协方差矩阵写为:
K
X
X
=
(
k
(
x
1
,
x
1
)
⋯
k
(
x
1
,
x
N
)
⋮
⋱
⋮
k
(
x
N
,
x
1
)
⋯
k
(
x
N
,
x
N
)
)
\boldsymbol{K}_{\boldsymbol{X X}}=\left( \begin{array}{ccc}{k\left(\boldsymbol{x}_{1}, \boldsymbol{x}_{1}\right)} & {\cdots} & {k\left(\boldsymbol{x}_{1}, \boldsymbol{x}_{N}\right)} \\ {\vdots} & {\ddots} & {\vdots} \\ {k\left(\boldsymbol{x}_{N}, \boldsymbol{x}_{1}\right)} & {\cdots} & {k\left(\boldsymbol{x}_{N}, \boldsymbol{x}_{N}\right)}\end{array}\right)
KXX=⎝⎜⎛k(x1,x1)⋮k(xN,x1)⋯⋱⋯k(x1,xN)⋮k(xN,xN)⎠⎟⎞
k
(
x
,
x
′
)
k\left(\boldsymbol{x}, \boldsymbol{x}^{\prime}\right)
k(x,x′) 为核函数,通常只要求这个核函数是正定函数,好让
K
X
X
\boldsymbol{K}_{\boldsymbol{X X}}
KXX 符合协方差矩阵的性质。比如对于服从零均值的 1 维变量
{
x
n
}
\left\{x_{n}\right\}
{xn} ,通常会用如下的核函数 (rational quadratic covariance function, RQ):
k
(
x
,
x
′
)
=
h
2
(
1
+
(
x
−
x
′
)
2
2
α
l
2
)
−
α
k\left(x, x^{\prime}\right)=h^{2}\left(1+\frac{\left(x-x^{\prime}\right)^{2}}{2 \alpha l^{2}}\right)^{-\alpha}
k(x,x′)=h2(1+2αl2(x−x′)2)−α
其中
h
,
α
,
l
h,\alpha, l
h,α,l 为正实数,我们称之为超参数,需要通过数据来训练得到。
**所以,给定一堆已有数据,怎么预测其他数据呢?**对于给定的输入
X
∈
R
N
×
D
\boldsymbol{X} \in \mathbb{R}^{N \times D}
X∈RN×D 以及对应的输出
f
∈
R
n
\boldsymbol{f} \in \mathbb{R}^{n}
f∈Rn, 当给定新的数据
X
∗
\boldsymbol{X}_{*}
X∗ 时,我们要能预测出对应的
f
∗
\boldsymbol{f}_{*}
f∗. 将训练集和测试集联合起来,得到多元高斯分布:
(
f
f
∗
)
∼
N
(
(
m
X
m
X
∗
)
,
(
K
X
X
K
X
X
∗
K
X
∗
X
K
X
∗
X
∗
)
)
\left( \begin{array}{c}{\boldsymbol{f}} \\ {\boldsymbol{f}_{*}}\end{array}\right) \sim \mathcal{N}\left(\left( \begin{array}{c}{\boldsymbol{m}_{\boldsymbol{X}}} \\ {\boldsymbol{m}_{\boldsymbol{X}_{*}}}\end{array}\right), \left( \begin{array}{c}{\boldsymbol{K}_{\boldsymbol{X X}} \quad \boldsymbol{K}_{\boldsymbol{X}\boldsymbol{X}_{*}}} \\ {\boldsymbol{K}_{\boldsymbol{X}_{*}\boldsymbol{X}} \quad \boldsymbol{K}_{\boldsymbol{X}_{*}\boldsymbol{X}_{*}}}\end{array}\right)\right)
(ff∗)∼N((mXmX∗),(KXXKXX∗KX∗XKX∗X∗))
对于熟悉 SLAM 的同学来说,要从这个分布中计算出
f
∗
\boldsymbol{f}_{*}
f∗, 采用舒尔补就行。从联合分布中计算
f
∗
\boldsymbol{f}_{*}
f∗, 也就是算条件概率:
p
(
f
∗
∣
X
∗
,
X
,
f
)
=
N
(
m
X
∗
+
K
X
∗
X
K
X
X
−
1
(
f
−
m
X
)
,
K
X
∗
X
∗
−
K
X
∗
X
K
X
X
−
1
K
X
X
∗
)
p\left(\boldsymbol{f}_{*} | \boldsymbol{X}_{*}, \boldsymbol{X}, \boldsymbol{f}\right)=\mathcal{N}\left(\boldsymbol{m}_{\boldsymbol{X}_{*}}+\boldsymbol{K}_{\boldsymbol{X}_{*} \boldsymbol{X}} \boldsymbol{K}_{\boldsymbol{X} \boldsymbol{X}}^{-1}\left(\boldsymbol{f}-\boldsymbol{m}_{\boldsymbol{X}}\right), \boldsymbol{K}_{\boldsymbol{X}_{*} \boldsymbol{X}_{*}}-\boldsymbol{K}_{\boldsymbol{X}_{*} \boldsymbol{X}} \boldsymbol{K}_{\boldsymbol{X} \boldsymbol{X}}^{-1} \boldsymbol{K}_{\boldsymbol{X}\boldsymbol{X}_{*}}\right)
p(f∗∣X∗,X,f)=N(mX∗+KX∗XKXX−1(f−mX),KX∗X∗−KX∗XKXX−1KXX∗)
对于不熟悉高斯分布这个操作的,可以参考这个回答。假设你已从训练数据中,得到了超参数,那就能用上述公司得到预测函数的分布,从而对于给定的
X
∗
\boldsymbol{X}_{*}
X∗, 采样出
f
∗
\boldsymbol{f}_{*}
f∗.
需要注意的是,训练的时候我们是用有噪声的数据 y n = f ( x n ) + ϵ n y_{n}=f\left(\boldsymbol{x}_{n}\right)+\epsilon_{n} yn=f(xn)+ϵn 来作为观测模型,其中 ϵ n ∼ N ( 0 , σ ϵ 2 ) \epsilon_{n} \sim \mathcal{N}\left(0, \sigma_{\epsilon}^{2}\right) ϵn∼N(0,σϵ2)。而之前的核函数并没有假设这个噪声,因此可以将这个噪声加到协方差 K X X \boldsymbol{K}_{\boldsymbol{X X}} KXX 的对角元素上。
高斯过程怎么训练
有了上述铺垫,训练其实就是最小二乘拟合那些超参数了。如有如下高斯过程,
f
∼
G
P
(
m
,
k
)
m
(
x
)
=
a
x
2
+
b
x
+
c
,
and
k
(
x
,
x
′
)
=
σ
y
2
exp
(
−
(
x
−
x
′
)
2
2
ℓ
2
)
+
σ
ϵ
2
δ
i
i
′
\begin{aligned} f & \sim \mathcal{G P}(m, k) \\ m(x) &=a x^{2}+b x+c, \quad \text { and } \quad k\left(x, x^{\prime}\right)=\sigma_{y}^{2} \exp \left(-\frac{\left(x-x^{\prime}\right)^{2}}{2 \ell^{2}}\right)+\sigma_{\epsilon}^{2} \delta_{i i^{\prime}} \end{aligned}
fm(x)∼GP(m,k)=ax2+bx+c, and k(x,x′)=σy2exp(−2ℓ2(x−x′)2)+σϵ2δii′
其中,超参数为
θ
=
{
a
,
b
,
c
,
σ
y
,
σ
n
,
ℓ
}
\theta=\left\{a, b, c, \sigma_{y}, \sigma_{n}, \ell\right\}
θ={a,b,c,σy,σn,ℓ},对于训练集,可以最大化似然概率:
L
=
log
p
(
y
∣
x
,
θ
)
=
−
1
2
log
∣
K
∣
−
1
2
(
y
−
m
)
⊤
K
−
1
(
y
−
m
)
−
n
2
log
(
2
π
)
L=\log p(\mathbf{y} | \mathbf{x}, \theta)=-\frac{1}{2} \log |K|-\frac{1}{2}(\mathbf{y}-\boldsymbol{m})^{\top} K^{-1}(\mathbf{y}-\boldsymbol{m})-\frac{n}{2} \log (2 \pi)
L=logp(y∣x,θ)=−21log∣K∣−21(y−m)⊤K−1(y−m)−2nlog(2π)
然后使用高斯牛顿或者共而梯度法求解就行,只是其中涉及到求偏导数:
∂
L
∂
θ
m
=
−
(
y
−
m
)
⊤
Σ
−
1
∂
m
∂
θ
m
∂
L
∂
θ
k
=
1
2
trace
(
K
−
1
∂
K
∂
θ
k
)
+
1
2
(
y
−
m
)
⊤
∂
K
∂
θ
k
K
−
1
∂
K
∂
θ
k
(
y
−
m
)
\begin{aligned} \frac{\partial L}{\partial \theta_{m}} &=-(\mathbf{y}-\boldsymbol{m})^{\top} \Sigma^{-1} \frac{\partial m}{\partial \theta_{m}} \\ \frac{\partial L}{\partial \theta_{k}} &=\frac{1}{2} \operatorname{trace}\left(K^{-1} \frac{\partial K}{\partial \theta_{k}}\right)+\frac{1}{2}(\mathbf{y}-\boldsymbol{m})^{\top} \frac{\partial K}{\partial \theta_{k}} K^{-1} \frac{\partial K}{\partial \theta_{k}}(\mathbf{y}-\boldsymbol{m}) \end{aligned}
∂θm∂L∂θk∂L=−(y−m)⊤Σ−1∂θm∂m=21trace(K−1∂θk∂K)+21(y−m)⊤∂θk∂KK−1∂θk∂K(y−m)
通常采用数值求导或者自动微分就可以了。
结语
主要是 17 年, 1GA9 年 ICRA 有两篇关于高斯过程学习里程计误差模型的论文,以及 GTSAM 那边 Jing Dong 也有几篇把高斯过程用在 SLAM 里的,所以才有了兴趣大致梳理下这个知识。
- 2019 ICRA, Learning Wheel Odometry and Imu Errors for Localization.
- 2017 ICRA,Gaussian Process Estimation of Odometry Errors for Localization and Mapping
这篇博客主要参考了: