1. 高斯过程(边缘化与条件作用)
1.1. 高斯过程的直观解释
1.1.1. 多元高斯分布
高斯分布（也叫做正态分布）是高斯过程的基础构件。而我们最感兴趣的是多元高斯分布，其每个随机变量都呈正态分布，联合分布也是高斯的。一般来说，多元高斯分布由均值向量\(\mu\)和协方差矩阵\(\sum\)定义。
_%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-09-08-44.png)
- 均值向量\(\mu\)描述了该分布的期望值，它的每个组件描述了对应维度的均值
- 协方差矩阵\(\sum\)对每个维度的方差进行建模，并确定不同随机变量之间的关联(协方差矩阵总是对称且半正定的)
- Σ 的对角线由第\(i\)个随机变量的标准差\(\sigma_i\) 组成，而非对角线的元素则描述了每个元素\(\sigma_{ij}\)之间的相关性
- \(X\)符合正态分布。协方差矩阵协方差矩阵\(\sum\)描述了该分布的形状
_%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-09-11-15.png)
1.1.1.1. 二元高斯分布
从图形上来看，该分布以均值为中心，由协方差矩阵决定其形状。下图展示了这些参数对于一个二维高斯分布的影响。每个随机变量的标准差在协方差矩阵的对角线上，而其它的值则显示了它们之间的协方差。
_%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-09-12-27.png)
高斯分布被广泛应用于为真实世界建模，有时在原分布未知的情况下作为替代品，有时用于中心极限定理。接下来我们会进一步讲解如何操纵高斯分布，以及如何从中获得有用的信息。
1.1.2. 边缘化与条件作用
高斯分布有一个很赞的代数性质：它在条件作用和边缘化情况下是封闭的。意思是，经过这些运算后，在结果中得到的分布依旧是高斯分布，这就使得很多统计学和机器学习中的问题变得易解。
边缘化: 对高斯分布进行边缘化和条件作用进行分解, 即\(P(X,Y)=P(X)P(Y|X)=P(Y)P(X|Y)\):
_%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-09-16-36.png)
其中 X 和 Y 代表原始随机变量的子集。给定随机变量 X 和 Y 组成的向量的正态概率分布\(P(X,Y)\), 我们可以用以下方法确定每个变量各自的边缘概率分布:
_%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-09-18-47.png)
这个公式所表达的意思很直接了当：X 和 Y 这两个子集各自只依赖于它们 μ 和 Σ 中对应的值。因此，要从高斯分布中边缘化一个随机变量，我们只需把μ 和Σ里那些对应的变量丢掉就行。
通过边缘化，我们可以获取多元概率分布的一部分信息(这里对y进行边缘化,去掉变量y,得到关于x的边缘分布):
_%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-09-24-21.png)
\[
\begin{aligned}
p_X(x)=\int_y p_{X,Y}(x,y)dy&=
\int_y p_{Y|X}(y|x)p_X(x)dy \\
&=
p_X(x)\int_y p_{Y|X}(y|x)dy \\
&=p_X(x)
\end{aligned}
\] 上式子的意思是:- 如果我们只想考虑X=x的情况, 即只对X=x的概率感兴趣, 那么我们要考虑变量Y所有可能的值, 需要对变量Y进行积分.(下面有图)
条件作用 高斯过程的另一个重要运算是条件作用，它可以用于得到一个变量在另一个变量条件下的概率分布。和边缘化类似，这个运算也是封闭的，会得到一个不同的高斯分布。条件运算是高斯过程的基石，它使贝叶斯推断成为可能。条件作用如下定义:
_%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-09-27-31.png)
新的均值只依赖于作为条件的变量，而协方差矩阵则和这个变量无关, 即:
- \(X \sim \mathcal{N}(\mu_X,\sum_{XX})\)
- 而 \(P(Y|X) \sim \mathcal{N}(\mu_Y+\sum_{YX}\sum_{XX}^{-1}(X-\mu_{X}),\sum_{YY}-\sum_{YX}\sum_{XX}\sum_{XY})\)
边缘化与条件作用的直观感受
- 边缘化可以理解为在高斯分布的一个维度上做累加，这也符合边缘分布的一般定义。
- 条件作用也有个很好的几何表达——我们可以把它想象成在多元分布上切下一刀，从而获得一个维数更少的高斯分布。
下面的例子是关于Y的边缘分布[MARGINALIZATION(Y)], 实际上就是去掉变量x, 把变量x边缘化掉 (注意哦,与前面的描述是相反的, 这个是去掉变量X)
\[ \begin{aligned} p_Y(y)=\int_x p_{X,Y}(x,y)dx&= \int_x p_{X|Y}(x|y)p_Y(y)dx \\ &= p_Y(y)\int_x p_{X|Y}(x|y)dx \\ &=p_Y(y) \end{aligned} \]
上式的意思就是说, 求出来的边缘分布\(P_Y\)就是给定一个y值, 在原来的联合概率分布\(P(X,Y)\)上, 给定y值之后把变量\(x\)在给定的y值方向上做累加.
- 图左是 Y 的边缘分布,是边缘化掉变量 X 的结果，类似于沿着 Y 轴把所有的X值做累加
- 图右是以给定的 X 为条件的关于变量Y的条件分布，类似于在原始分布上切下一刀
_%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-09-48-37.png)
_%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-09-38-08.png)
_%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-10-15-22.png)
图片来自:看得见的高斯过程(原文)
1.2. 相关理论推导
1.2.1. 高斯分布的表示
- 协方差矩阵+均值
- 信息矩阵+信息矢量
常见的是"协方差矩阵+均值"形式, \[p(x)=\eta \exp\{-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)\}\]
其中对称正定矩阵Σ为随机变量x的协方差矩阵，μ为x的均值，简记为 \[p(x) = N(\mu, \Sigma)\]
信息矩阵+信息矢量的形式可以由上式推导而来 \[ \begin{aligned} p(x)&=\eta \exp\{-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)\} \\ &=\eta\exp\{-\frac{1}{2}x^T\Sigma^{-1}x+x^T\Sigma^{-1}\mu\} \end{aligned} \] > 其中,运算中产生的常数项都全部吸收到了 η 中.
现在定义信息矩阵\(\Lambda=\Sigma^{-1}\), 信息矢量\(\xi=\Sigma^{-1}\mu=\Lambda\mu\), 则有 \[ p(x)=\eta\exp\{-\frac{1}{2}x^T\Lambda x+x^T\xi\} \] 即 \[p(x) = N^{-1}(\xi, \Lambda)\]
1.2.2. 联合高斯分布的分解
设随机变量\(X, Y\)满足联合高斯分布\(p(X,Y)\)
于是有 \[ p(X,Y)=p(X)p(Y|X)\]
所谓边缘化操作就是求出上面的分布p(X) * 也就是把变量Y边缘化掉(去掉)
用"协方差矩阵+均值"形式给出联合分布的表示: \[ p(X,Y) = N\Bigg(\begin{pmatrix} \mu_X \\ \mu_Y \end{pmatrix}, \begin{pmatrix} \Sigma_{XX} & \Sigma_{XY} \\ \Sigma_{YX} & \Sigma_{YY} \end{pmatrix}\Bigg) \]
补充:舒尔补公式
- \[ \begin{aligned} \begin{bmatrix} I & -BD^{-1} \\ 0 & I \end{bmatrix} \begin{bmatrix} A & B \\ C & D \end{bmatrix} = \begin{bmatrix} \Delta_D & 0 \\ C & D \end{bmatrix} \end{aligned} \] \[ \Delta_D=A-BD^{-1}C \]
- \[ \begin{aligned} \begin{bmatrix} I & 0 \\ -CA^{-1} & I \end{bmatrix} \begin{bmatrix} A & B \\ C & D \end{bmatrix} = \begin{bmatrix} A & B \\ 0 & \Delta_A \end{bmatrix} \end{aligned} \] \[ \begin{aligned} \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} I & -A^{-1}B \\ 0 & I \end{bmatrix} = \begin{bmatrix} A & 0 \\ C & \Delta_A \end{bmatrix} \end{aligned} \] \[\Delta_A=D-CA^{-1}B\]
将上面两式联合起来,就可以把矩阵变成对角型 \[ \begin{aligned} \begin{bmatrix} I & 0 \\ -CA^{-1} & I \end{bmatrix} \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} I & -A^{-1}B \\ 0 & I \end{bmatrix} = \begin{bmatrix} A & 0 \\ 0 & \Delta_A \end{bmatrix} \end{aligned} \] 也可以通过对角型恢复出原来的矩阵 \[ \begin{aligned} \begin{bmatrix} I & 0 \\ CA^{-1} & I \end{bmatrix} \begin{bmatrix} A & 0 \\ 0 & \Delta_A \end{bmatrix} \begin{bmatrix} I & A^{-1}B \\ 0 & I \end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \end{aligned} \] 就可以快速写出矩阵的逆 \[ \begin{aligned} \begin{bmatrix} A & B \\ C & D \end{bmatrix}^{-1} &= \begin{bmatrix} I & A^{-1}B \\ 0 & I \end{bmatrix}^{-1} \begin{bmatrix} A & 0 \\ 0 & \Delta_A \end{bmatrix}^{-1} \begin{bmatrix} I & 0 \\ CA^{-1} & I \end{bmatrix}^{-1} \\ &= \begin{bmatrix} I & -A^{-1}B \\ 0 & I \end{bmatrix} \begin{bmatrix} A^{-1} & 0 \\ 0 & \Delta_A^{-1} \end{bmatrix} \begin{bmatrix} I & 0 \\ -CA^{-1} & I \end{bmatrix} \end{aligned} \]
联合高斯分布p(X,Y)p(X,Y)的概率密度函数("信息矩阵+信息矢量形式")的表示为: \[
\begin{aligned}
p(X, Y)&=
\eta \exp\Bigg\{-\frac{1}{2}\begin{pmatrix}
X-\mu_X \\
Y-\mu_Y
\end{pmatrix}^T\begin{pmatrix}
\Sigma_{XX} & \Sigma_{XY} \\
\Sigma_{YX} & \Sigma_{YY}
\end{pmatrix}^{-1}\begin{pmatrix}
X-\mu_X \\
Y-\mu_Y
\end{pmatrix}\Bigg\} \\
&\propto
\exp\Bigg\{-\frac{1}{2}\begin{pmatrix}
a \\
b
\end{pmatrix}^T\begin{pmatrix}
A & C^{T} \\
C & D
\end{pmatrix}^{-1}\begin{pmatrix}
a \\
b
\end{pmatrix}\Bigg\}
\end{aligned}
\] _%E9%AB%98%E6%96%AF%E8%BF%87%E7%A8%8B(%E8%BE%B9%E7%BC%98%E5%8C%96%E4%B8%8E%E6%9D%A1%E4%BB%B6%E4%BD%9C%E7%94%A8)/imgs/2020-02-13-12-29-56.png)
于是, 就得到了关于X的边缘分布\(p(X)\)和关于Y的条件概率分布\(p(Y|X)\)
关于X的边缘分布\(p(X)\) \[ \begin{aligned} p(X) &=\eta_1 \exp \Big\{-\frac{1}{2} (X-\mu_X)^{T} \Sigma_{XX}^{-1}(X-\mu_X) \Big\} \sim \mathcal{N}(0,\Sigma_{XX}) \end{aligned} \] > 这说明了,边缘分布\(p(X)\)的协方差就是联合分布\(P(X,Y)\)的协方差矩阵中对应的矩阵块\(\Sigma_{XX}\)
关于Y的条件概率分布 \[ \begin{aligned} p(Y|X)&= \eta_2\exp\{-\frac{1}{2}[Y-(\mu_Y+\Sigma_{YX}\Sigma_{XX}^{-1}(X-\mu_X))]^T\Theta_{A}^{-1}[Y-(\mu_Y+\Sigma_{YX}\Sigma_{XX}^{-1}(X-\mu_X))]\} \\ & \propto \eta_2\exp\{-\frac{1}{2}[b-CA^{-1}a]^T\Theta_{A}^{-1}[b-CA^{-1}a]\} \\ \end{aligned} \]
其中,\(\Theta_A=\Delta_A=D-CA^{-1}B=D-CA^{-1}C^{T}\) 表示矩阵块A在联合概率分布\(P(X|Y)\)的协方差矩阵中的舒尔补
所以, \[ \Theta_A=\Delta_A=\Sigma_{YY}-\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY} \]
这说明了,关于Y的条件分布\(p(Y|X)\)的协方差就是联合分布\(P(X,Y)\)的协方差矩阵关于子块\(A=\Sigma_{XX}\)的舒尔补\(\Delta_A=\Sigma_{YY}-\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY}\)
假定现在我们只有联合概率分布\(p(X,Y)\)的信息矩阵 现在来讨论上面两个分布的信息矩阵
首先利用舒尔补写出联合概率分布\(p(X,Y)\)的信息矩阵 (也就是协方差矩阵的逆) \[ \begin{aligned} \begin{bmatrix} A & B \\ C & D \end{bmatrix}^{-1} &= \begin{bmatrix} I & -A^{-1}B \\ 0 & I \end{bmatrix} \begin{bmatrix} A^{-1} & 0 \\ 0 & \Delta_A^{-1} \end{bmatrix} \begin{bmatrix} I & 0 \\ -CA^{-1} & I \end{bmatrix} \end{aligned} \]
\[ \begin{aligned} \begin{bmatrix} A & C^{T} \\ C & D \end{bmatrix}^{-1} &= \begin{bmatrix} A^{-1}+A^{-1}C^{T}\Delta_{A}^{-1}CA^{-1} & -A^{-1}C^{T}\Delta_{A}^{-1} \\ -\Delta_{A}^{-1}CA^{-1} & \Delta_{A}^{-1} \end{bmatrix} \\ &= \begin{bmatrix} \Sigma_{XX}^{-1}+\Sigma_{XX}^{-1}\Sigma_{XY}\Theta_{A}^{-1}\Sigma_{YX}\Sigma_{XX}^{-1} & -\Sigma_{XX}^{-1}\Sigma_{XY}\Theta_{A}^{-1} \\ -\Theta_{A}^{-1}\Sigma_{YX}\Sigma_{XX}^{-1} & \Theta_{A}^{-1} \end{bmatrix} \\ &= \begin{bmatrix} \Lambda_{XX} & \Lambda_{XY} \\ \Lambda_{YX} & \Lambda_{YY} \end{bmatrix} \end{aligned} \]
- 关于X的边缘分布\(p(X)\)的信息矩阵 > 从前面的讨论, 已经知道了关于X的边缘分布\(p(X)\)的信息矩阵就是联合分布\(P(X,Y)\)的协方差矩阵中对应的矩阵块\(\Sigma_{XX}\), 那么对应的信息矩阵就是\(\Sigma_{XX}^{-1}\)
那么现在的目的就是: 利用原有的联合概率分布信息矩阵, 求出关于X的边缘分布\(p(X)\)的信息矩阵, 也就是边缘化掉变量Y之后的信息矩阵
易得 \[ \begin{aligned} \Sigma_{XX}^{-1}=\Lambda_{XX}-\Lambda_{XY}\Lambda_{YY}^{-1}\Lambda_{YX} \end{aligned} \]
- 关于Y的条件概率分布的信息矩阵 > 从前面的讨论, 已经知道了关于Y的条件概率分布的信息矩阵就是联合分布\(P(X,Y)\)的协方差矩阵关于子块\(A=\Sigma_{XX}\)的舒尔补\(\Theta_A=\Delta_A\), 那么对应的信息矩阵就是\(\Theta_A^{-1}\)
易得 \[ \Theta_A^{-1}=\Lambda_{YY} \]
1.3. 边缘化操作的作用
所谓边缘化，就是求某个联合概率分布的边缘分布。比如对于联合概率\(p(X,Y)\)，对\(Y\)进行边缘化，就是对\(Y\)在整个空间中积分，即 \[ \begin{aligned} p_X(x)=\int_y p_{X,Y}(x,y)dy&= \int_y p_{Y|X}(y|x)p_X(x)dy \\ &= p_X(x)\int_y p_{Y|X}(y|x)dy \\ &=p_X(x) \end{aligned} \]
因此，对\(Y\)边缘化的结果就是X的边缘化分布的\(p(X)\)，仍然是一个高斯函数。伴随着边缘化，\(p(Y|X)\)就是\(p(X,Y)\)对\(X\)的条件化。
最小二乘
在信息矩阵+信息矢量的表示方式下，边缘化和条件化与最小二乘法有密切关系。在许多基于最小二乘的优化问题中，常有如下形式的优化目标： \[ \min_x\Vert e(x)\Vert_W^2=\min_x~e(x)^TW^{-1}e(x) \] 其中\(W\)是\(e(x)\)的协方差矩阵。
高斯牛顿
为了寻找上式的最小值，常使用迭代优化的方法，每一次迭代都会寻找一个增量\(\Delta x\)使目标函数减小。为了求增量，往往会将\(e(x)\)在当前\(x\)处展开为一阶近似（这种处理方式即Gauss-Newton Method），即: \[ e(x+\Delta x)\simeq e(x)+J(x)\Delta x \] 其中, \[ J(x)=\frac{\partial e}{\partial x} \] 于是, 可以生成新的优化目标: \[ \min_{\Delta x}~[e(x)+J(x)\Delta x]^TW^{-1}[e(x)+J(x)\Delta x] \] 这是关于\(\Delta x\)的二次函数，对\(\Delta x\)求导，并令导数等于0，有: \[ \begin{aligned} J(x)^TW^{-1}J(x)\Delta x+J(x)^TW^{-1}e(x)=0 \\ \Longrightarrow J(x)^TW^{-1}J(x)\Delta x = -J(x)^TW^{-1}e(x) \end{aligned} \] 令\(J(x)^TW^{-1}J(x) = H\), 是变量\(\Delta x\)的信息矩阵
令\(-J(x)^TW^{-1}e(x)=b\), 则有 \[ H\Delta x=b \] 这就是非线性优化中的增量方程, 通过不断迭代求解\(\Delta x\), 更新变量\(x+=\Delta x\), 使得目标函数不断下降, 直到达到要求.
SLAM问题求解
在很多优化问题中，待优化的变量有明确意义，比如在SLAM或者SfM问题中，要优化的是所有相机的位姿\(c\)以及地图中所有三维点(路标点)的坐标\(p\)，设\(\Delta x\)由这两个分量的增量构成，即 \[ \Delta x = \begin{pmatrix} \Delta c \\ \Delta p \end{pmatrix} \] 并且求出了对应的信息矩阵\(\Lambda\) \[ \Lambda=\begin{pmatrix} \Lambda{cc} & \Lambda{cp} \\ \Lambda{pc} & \Lambda{pp} \end{pmatrix}, ~~~ b = \begin{pmatrix} b_c \\ b_p \end{pmatrix} \] 于是有 \[ \begin{pmatrix} \Lambda{cc} & \Lambda{cp} \\ \Lambda{pc} & \Lambda{pp} \end{pmatrix} \begin{pmatrix} \Delta c \\ \Delta p \end{pmatrix}= \begin{pmatrix} b_c \\ b_p \end{pmatrix} \] 直接对信息矩阵求逆可以对上述方程求解, 但是当矩阵维度过大的时候, 直接求逆很难, 通常采用各种分解
可以采用上面的舒尔补公式进行消元, 也就是所说的高斯消元法 >(1) 舒尔补公式 >\[ \begin{aligned} \begin{bmatrix} I & -BD^{-1} \\ 0 & I \end{bmatrix} \begin{bmatrix} A & B \\ C & D \end{bmatrix} = \begin{bmatrix} \Delta_D & 0 \\ C & D \end{bmatrix} \end{aligned} >\] >\[ \Delta_D=A-BD^{-1}C >\] 即对方程两边同时左乘 \[ \begin{bmatrix} I & -\Lambda{cp}\Lambda{pp}^{-1} \\ 0 & I \end{bmatrix} \] 得到一个下三角型的形式 \[ \begin{pmatrix} \Lambda{cc}-\Lambda{cp}\Lambda{pp}^{-1}\Lambda{pc} & 0 \\ \Lambda{pc} & \Lambda{pp} \end{pmatrix} \begin{pmatrix} \Delta c \\ \Delta p \end{pmatrix}= \begin{pmatrix} b_c-\Lambda{cp}\Lambda{pp}^{-1}b_c \\ b_p \end{pmatrix} \] 于是原方程可以转换为两个独立方程 \[ \begin{aligned} (\Lambda{cc}-\Lambda{cp}\Lambda{pp}^{-1}\Lambda{pc})\Delta c &= b_c-\Lambda{cp}\Lambda{pp}^{-1}b_c \\ \Lambda_{pp}\Delta p &= b_p -\Lambda_{pc}\Delta c \end{aligned} \]
可以发现，\(\Delta c\)的系数矩阵，与上文中边缘分布\(p(X)\)的信息矩阵和信息矢量有相同的形式。而\(\Delta p\)的系数矩阵，则与关于Y的条件概率分布\(p(Y|X)\)的信息矩阵有相同形式。
也就是说，这里的高斯消元法，等价于对变量\(\Delta c\)做了边缘化，先将\(\Delta p\)边缘化掉,就是消掉路标点p，单独求相机位姿增量\(\Delta c\)，然后再在\(\Delta c\)已知的情况下求\(\Delta p\)。
参考: 1. 高斯分布与边缘化 2. 看得见的高斯过程(译文) 3. 看得见的高斯过程(原文)