前言

经常在论文中看到一些熟悉又陌生的算法名，例如EM，例如GMM，HMM等等，这些算法好像都是知道个大概，但是要我自己去实现，感觉总是这里缺先验那里缺先验，不能完整地复述出来，就是不是自己的东西，有必要像矩阵求导那样，推导一个完整的式子来。

期望最大化算法

英文原名：Expectation-Maximization，缩写名：EM

问题描述

给定某个概率分布X\mathcal{X}的系列观测值X(θ)X(\theta)，其中θ\theta是某个未知的，需要被优化的参数，现假设待观测值Z XZ~\mathcal{X}，参数估计要求算法使得似然函数L(θ;X)L(\theta;X)取得最大值。即使得在出现观测值XX的概率最大：

L(θ;X)=p(X∣θ)=∑z∈Zp(X,Z∣θ)zL(\theta;X)=p(X|\theta)=\sum\limits_{z\in\mathcal{Z}}p(X,Z|\theta)z

算法核心

由于不能直接得到隐变量θ\theta的值，且随机变量ZZ的分布情况也未知，不能直接使用求导等方式求最值，可通过迭代求解。EM只是提供了一个计算的框架，而没有给出最优化的具体方法，这需要使用最优化的理论来解决之。

期望最大化算法就分为两步：期望+最大化

求解期望

直接使用期望的定义式求出在当前似然函数L(X,Z∣θ)L(X,Z|\theta)情况下的数学期望

EZ∼X(θ)[p(X,Z∣θ)]\displaystyle E_{Z\sim\mathcal{X}(\theta)}[p(X,Z|\theta)]

使用中，通常遇到高斯分布的随机变量，因此为了简化上述表示，且不影响最值的位置，可求对数期望：

EZ∼X(θ)[ln⁡p(X,Z∣θ)]\displaystyle E_{Z\sim\mathcal{X}(\theta)}[\ln p(X,Z|\theta)]

最大化期望

在当前（对数）期望表达式已知的情况下，在θ\theta的定义域内，使用搜索算法找到一个最优解，使得：

θ˜=arg max⁡θEZ∼X(θ)[p(X,Z∣θ)]\displaystyle \~\theta=\argmax_{\theta}E_{Z\sim\mathcal{X}(\theta)}[p(X,Z|\theta)]

然后用这个搜索得到的最优θ˜\~\theta代替之前分布参数θ\theta，获得新的参数分布：

θ:=θ˜\theta:=\~\theta

然后重新计算期望，并进行下一步迭代。

使用小结

本算法有一点强化学习的味道，表演家得到结果，评论家对表演结果作出评价，然后从最好的表演结果中选择一个作为下一次表演的基准，使模型能够作出最佳的表演行为（最大化期望）。

应用实例推导

高斯混合模型GMM是期望最大化算法的直接应用，其问题描述为：

现有一线性组合的高斯叠加模型G(X)=∑k=1KλkN(μk,σk)G(X)=\sum\limits_{k=1}^{K}\lambda_kN(\mu_k,\sigma_k)，通过先验知识和前期观测操作得到了一系列观测数据集{xi}\{x_i\}，常见的一种情况是，在高斯分布类别数目KK已知的情况下，求解各项高斯分布函数的参数。

RANSAC方法

随机抽样一致性方法Random Sample Consensus，是用于从带有较多噪声的数据样本中提取准确模型分布、并完成模型拟合的稳健方法。该方法最早于1981年由Fischler等人[1]提出，用于解决一类高度噪声化的位置确定问题Locatioon Determination Problem。传统的最小二乘方法将使用全部数据样本点，而不管该样本点是否符合模型的拟合，因此高度偏离模型的点将会对拟合结果产生显著的不利影响；即使使用加权的最小二乘方法，也需要人工确定点的权重，一定程度上不是自动化实现的。RANSAC的提出解决了这一难题。

算法表述

设某个模型需要的最少标定点为nn，现有一个带噪声的标注集PP，其中标注点数多于nn。按照以下逻辑进行模型的拟合：

  1. 从标注集PP中随机抽取nn个点的子集S1S_1，用其完成模型参数拟合，得到结果θ1\theta_1；
  2. 对子集S1S_1的补集，使用一次模型θ1\theta_1计算残差，将其中残差小于某个容许值ϵ\epsilon的点加入至当前集合S1S_1中，并排除其中残差大于容许值的点，得到一致性集合S1∗S_1^*；
  3. 判断一致性集合的基数，当∣S1∗∣|S_1^*|大于预设的阈值tt后，用集合S1∗S_1^*计算新的模型θ2\theta_2，并重复执行步骤2；
  4. 否则∣S1∗∣<t|S_1^*| < t时，重新随机抽取nn个点的子集S2S_2，重复执行步骤1。
  5. 如果在某个预设的步骤数中，取得的一致性集合基数都小于阈值tt，宣告模型拟合失败，强制终止。否则，以最后一步取得的一致性集合为拟合数据集S∗S^*，完成输出模型θ∗\theta^*的拟合。

说明

关于RANSAC算法部分超参数的设置，论文原文[1]中提到了以下几点注意事项：

  1. 最大迭代次数不宜过少，以避免在随机抽样时最不利情况下错过模型的拟合，应当保证抽样次数k≥log⁡(1−z)log⁡(1−b)k\geq \frac{\log(1-z)}{\log(1-b)}，其中zz是需要的置信度，b=wnb=w^n是nn个抽样点计算得到的残差都恰好位于容许值以内的概率。
  2. 关于基数阈值tt，作者认为应当足够大，以保证一致性集合包含足够多的点，便于解出正确的模型与后期的模型平滑。

双目立体视觉

这玩意在运动物体成像中也有应用，被称为homography estimation，即同质性估计，是一种计算机视觉中的基础问题，用于计算两个平面之间的映射关系，常用于计算相机的运动轨迹，或者计算两个平面之间的变换关系。

双目系统极线几何关系

如下图，对极几何关系是双目系统中最为重要的几何关系之一。
对极几何关系
首先介绍几个基本的概念：

基本概念

  1. 基线：两个相机的光心之间的距离，图中的O1O2O_1O_2；
  2. 极点：两个相机的光心在另一个相机的成像平面上的投影点，图中的e1,e2e_1,e_2；
  3. 极平面：两个相机的光心和空间点P所在的平面，图中的O1O2PO_1O_2P；
  4. 极线：极平面与图像平面的交线，图中的l1,l2l_1,l_2；
  5. 极平面簇：基线随空间点P变化所形成的平面簇，所有的极平面相交于基线。

极线约束

称图中p2p_2极点与极线l1l_1满足对极关系，同理可称图中的p1p_1极点和极线l2l_2满足对极关系。在双目立体视觉中，极点不能与极点形成对应，而只能与极线对应，由此引出了极线约束的概念。极线约束是指，极点所对应的极点必定在对极线上。该约束关系极大了减少了搜索的范围，从而提高了立体视觉算法效率。约束方程如下：

{slpl=HlPsrpr=HrP\displaystyle \left\{ \begin{aligned} s_l\mathbf{p_l=H_lP}\\ s_r\mathbf{p_r=H_rP} \end{aligned} \right.

式中，Hl,Hr\mathbf{H_l,H_r}是左右两个相机的内外参矩阵之乘积，即变换矩阵，可得：

Hl=Kl[Rl∣tl]≜[Ml∣ml]Hr=Kr[Rr∣tr]≜[Mr∣mr]\displaystyle \mathbf{H_l=K_l[R_l|t_l]}\triangleq[\mathbf{M}_l|\mathbf{m}_l]\\ \mathbf{H_r=K_r[R_r|t_r]}\triangleq[\mathbf{M}_r|\mathbf{m}_r]

利用中间量P=[X,Y,Z,1]T=[PXYZT,1]\mathbf{P}=[X,Y,Z,1]^\mathrm{T}=[\mathbf{P}_{XYZ}^\mathrm{T},1]消元，可得：

(slpl−ml)Ml−1=(srpr−mr)Mr−1  ⟹  srpr−slMrMl−1pl=mr−MrMl−1ml\begin{aligned} &(s_l\mathbf{p_l}-\mathbf{m}_l)\mathbf{M}_l^{-1}&=(s_r\mathbf{p_r}-\mathbf{m}_r)\mathbf{M}_r^{-1}\\ &\implies s_r\mathbf{p}_r-s_l\mathbf{M_rM_l^{-1}p_l}&=\mathbf{m_r-M_rM_l^{-1}m_l} \end{aligned}

上式即为极线约束的第一种形式，可根据图像坐标p1,p2\mathbf{p_1,p_2}的齐次化消去尺度因子sl,srs_l,s_r。经过齐次化后，可以得到第二种形式，以下为推导过程：

记上式右边的向量为m\mathbf{m}，则根据齐次向量相等关系（即向量与自己的向量积为0），对左式向量积m\mathbf{m}，可得：

m×(srpr−slMrMl−1pl)=0  ⟹  m×pr=m×sMrMl−1pl\begin{aligned} &\mathbf{m}\times(s_r\mathbf{p_r}-s_l\mathbf{M_rM_l^{-1}p_l})&=0\\ &\implies \mathbf{m}\times\mathbf{p_r}&=\mathbf{m}\times s\mathbf{M_rM_l^{-1}p_l} \end{aligned}

可见左式得到的结果是一个向量，且这个向量必定与pr\mathbf{p}_r正交，因此，其与pr\mathbf{p}_r的内积为0，即：

prTm×pr=prTm×sMrMl−1pl  ⟹  prTm×sMrMl−1pl=0  ⟹  prTm×MrMl−1pl=0\begin{aligned} &\mathbf{p_r^\mathrm{T}m}\times\mathbf{p_r}&=\mathbf{p_r^\mathrm{T}m}\times s\mathbf{M_rM_l^{-1}p_l}\\ &\implies\mathbf{p_r^\mathrm{T}m}\times s\mathbf{M_rM_l^{-1}p_l}&=0\\ &\implies\mathbf{p_r^\mathrm{T}m}\times\mathbf{M_rM_l^{-1}p_l}&=0 \end{aligned}

可见消元完成的结果，可以写为二次型的形式，本质上该式是极平面方程的表达式，即：

prTFpl=0\mathbf{p_r^\mathrm{T}}\mathbf{Fp_l}=0

两者的对极线关系可以非常明显地写出：

l2=Fpl⟷pll1=FTpr⟷prl_2=\mathbf{Fp_l}\longleftrightarrow p_l\\ l_1=\mathbf{F^\mathrm{T}p_r}\longleftrightarrow p_r

约束方程

本部分参考文献 Determining the Epipolar Geometry and its Uncertainty: A Review。
下面推导矩阵F\mathbf{F}的具体形式，从上式中可以看出，该矩阵只与左右两个摄像机的内参数和两者摆放的相对位置（结构参数），进一步地，选取左相机图像坐标系为参考坐标系，则可以将两个相机的成像方程写为：

Hl=Al[I∣0]Hr=Ar[R∣t]\begin{aligned} \mathbf{H_l}=\mathbf{A_l[I|0]}\\ \mathbf{H_r}=\mathbf{A_r[R|t]} \end{aligned}

则成像结果为：

{srpr=ArRPXYZ+Artslpl=AlPXYZ\left\{\begin{aligned} s_r\mathbf{p}_r&=\mathbf{A_rRP_{XYZ}+A_rt}\\ s_l\mathbf{p}_l&=\mathbf{A_lP_{XYZ}} \end{aligned}\right.

注意到这里的R\mathbf{R}是一个4×34\times3的矩阵，其中每列是正交的，因此其不是方阵，对上述方程作变形可得：

{srAr−1pr−t=RPXYZslRAl−1pl=RPXYZ  ⟹  srAr−1pr−slRAl−1pl=t\left\{\begin{aligned} s_r\mathbf{A_r^{-1}p}_r-\mathbf{t}&=\mathbf{RP_{XYZ}}\\ s_l\mathbf{RA_l^{-1}p}_l&=\mathbf{RP_{XYZ}} \end{aligned}\right.\\ \implies s_r\mathbf{A_r^{-1}p}_r-s_l\mathbf{RA_l^{-1}p}_l=\mathbf{t}

同理，对左式向量积右式使得上式归零，可得：

[t]×(Ar−1pr−sRAl−1pl)=0\begin{aligned} [\mathbf{t}]_\times (\mathbf{A_r^{-1}p}_r-s\mathbf{RA_l^{-1}p}_l)=0 \end{aligned}

然后，内积括号里面的第一项，使得前一项归零，从而消去尺度因子ss，可得：

prTAr−T[t]×RAl−1pl=0\mathbf{p_r^\mathrm{T}A_r^\mathrm{-T}[t]_\times RA_l^\mathrm{-1}p_l}=0

从而得到《机器视觉》（张广军著）书上111页，或者论文中第163页公式（1）的结果式，称该结果为共平面约束，所得矩阵F=Ar−T[t]×RAl−1\mathbf{F=A_r^\mathrm{-T}[t]_\times RA_l^\mathrm{-1}}称为基本矩阵。此式即为三角测量、双目视觉中最重要的约束公式。

参考文献


[1] Fischler M. A. , Bolles R. C. .Random sample Consensus[J/OL].Commun. ACM,1981,24(6):381-395