本文是关于 coursera 上 《Robotics: Estimation and Learning》 课程的笔记。
前面讲了一维和多维高斯分布的相关知识。
但是在某些情况下,使用 单高斯模型(Gaussian single model, GSM) 会有一些局限。
在现实世界中我们需要学习的目标可能符合这样的分布 :
如上图所示,当你用单高斯模型去拟合它时,得到这样的曲线。
显然它不能很好地表征目标。这样的目标有多种模式,或者缺乏对称性。
你将看到混合高斯模型的表现力则很好,好到可以建模任意的分布。
简单地说,混合高斯模型是一些高斯模型的和。
下面这个图中,那些彩色的线是十个随机的高斯模型曲线,黑线则是它们的和。
一个混合后的曲线可以有非常怪异的形状,这是无法用一个简单的方程描述出来的。
所以,如果我们选择了合适的高斯模型,我们就可以表征任意一个分布 。
看回小球颜色的例子。
右边的 2 2 2 维图展示了红色和绿色通道的颜色分布,我们可以用一个 2 D 2\text{D} 2D 高斯模型去建模这个分布:
或者我们可以用两个高斯模型的混合去建模:
用混合的方式,看起来更好滴描述了这两个颜色的分布情况。
现在用我们最喜爱的工具 ---- 数学😂,来更加正规地描述混合高斯模型 。
高斯混合模型(Gaussian Mixture Model,GMM)
让
g
g
g 代表单高斯模型的概率密度函数,则混合高斯模型可以被写成:
p
(
x
)
=
∑
k
=
1
K
w
k
g
k
(
x
∣
μ
k
,
Σ
k
)
(1)
p(\textbf{\textit{x}}) = \sum^{K}_{k=1} w_k g_k \left( \textbf{\textit{x}} | \boldsymbol{\mu}_k , \Sigma_k \right)\tag{1}
p(x)=k=1∑Kwkgk(x∣μk,Σk)(1)上面表示
k
k
k 个不同参数的高斯模型的加权和,这些权重都是正的,其中:
g
k
g_k
gk 表示均值为
μ
k
\boldsymbol{\mu}_k
μk,协方差为
Σ
k
\Sigma_k
Σk 的单高斯模型。
w
k
w_k
wk 表示混合系数(权值)。
权值的和为 1 1 1,保证了混合高斯模型的密度函数积分仍然是 1 1 1。
如果我们可以使用任意多的单模型,并且这些模型有任意小的方差,我们在理论上就可以表示任意形状的分布,这正是混合高斯模型强大的地方 。
虽然混合高斯模型很强,但是它的代价是需要很多参数: μ = { μ 1 , μ 2 , … , μ K } Σ = { Σ 1 , Σ 2 , … , Σ K } w = { w 1 , w 2 , … , w K } K : Number of Components \boldsymbol{\mu} = \left\{ \boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \dots , \boldsymbol{\mu}_K \right\} \\ \Sigma = \left\{ \Sigma_1, \Sigma_2, \dots, \Sigma_K \right\} \\ \textbf{\textit{w}} = \left\{ w_1, w_2, \dots, w_K \right\} \\ K: \text{Number of Components} μ={μ1,μ2,…,μK}Σ={Σ1,Σ2,…,ΣK}w={w1,w2,…,wK}K:Number of Components
在同样的维数下,混合高斯模型比单高斯模型有更多的参数。
当混合的模型数量增加,均值和协方差矩阵的数量也随之增加 。
而且还多了一个新的参数:权值。
还有高斯模型的数量
K
K
K 本身也是一个参数,你需要确定它。
这么多的参数带来一些副作用:
首先,参数估计变得更困难,后面会看到,我们没有混合高斯模型的解析解。
第二,犯错的概率变大了,特别是我们可能遇到过拟合的问题。谨记,我们需要对模型的复杂程度特别小心。
这里为了简单起见,我们只考虑给定的
K
K
K ,以及恒定权值
w
w
w 的情况:
μ
=
{
μ
1
,
μ
2
,
…
,
μ
K
}
Σ
=
{
Σ
1
,
Σ
2
,
…
,
Σ
K
}
w
=
w
=
1
/
K
K
:
指定的数字
\boldsymbol{\mu} = \left\{ \boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \dots , \boldsymbol{\mu}_K \right\} \\ \Sigma = \left\{ \Sigma_1, \Sigma_2, \dots, \Sigma_K \right\} \\ \textbf{\textit{w}} = w = 1/K \\ K: \text{指定的数字}
μ={μ1,μ2,…,μK}Σ={Σ1,Σ2,…,ΣK}w=w=1/KK:指定的数字在此基础上探讨估计模型的均值和协方差矩阵。
将参数量简化后,我们来尝试做似然估计,求解混合高斯模型的参数。
简单回顾一下,和前面一样,我们的似然函数为:
μ
^
,
Σ
^
=
arg
max
μ
,
Σ
p
(
{
x
i
}
∣
μ
,
Σ
)
μ
=
{
μ
k
}
Σ
=
{
Σ
k
}
k
=
1
,
2
,
…
,
K
\hat{\boldsymbol{\mu}},\hat{\Sigma} = \arg \max_{\boldsymbol{\mu},\Sigma} p(\{\textbf{\textit{x}}_i\} |\boldsymbol{\mu},\Sigma ) \\[1em] \boldsymbol{\mu} = \left\{ \boldsymbol{\mu}_k\right\} \\ \Sigma = \left\{ \Sigma_k \right\} \\ k = 1,2,\dots,K
μ^,Σ^=argμ,Σmaxp({xi}∣μ,Σ)μ={μk}Σ={Σk}k=1,2,…,K
最大似然估计意味着,我们希望找到模型的一组参数,使得产生观测数据的可能性最大。
我们不最大化联合概率,而是假设观测都是独立的,然后最大化这些概率的乘积:
μ
^
,
Σ
^
=
arg
max
μ
,
Σ
∏
i
=
1
N
p
(
x
i
∣
μ
,
Σ
)
\hat{\boldsymbol{\mu}},\hat{\Sigma} = \arg \max_{\boldsymbol{\mu},\Sigma} \prod^N_{i=1} p(\textbf{\textit{x}}_i |\boldsymbol{\mu},\Sigma )
μ^,Σ^=argμ,Σmaxi=1∏Np(xi∣μ,Σ)
然后取对数,并且转为最大化对数概率的和: μ ^ , Σ ^ = arg max μ , Σ ∑ i = 1 N ln p ( x i ∣ μ , Σ ) ‾ (2) \hat{\boldsymbol{\mu}},\hat{\Sigma} = \arg \max_{\boldsymbol{\mu},\Sigma} \sum^N_{i=1} \ln \underline{p(\textbf{\textit{x}}_i |\boldsymbol{\mu},\Sigma )}\tag{2} μ^,Σ^=argμ,Σmaxi=1∑Nlnp(xi∣μ,Σ)(2)
到目前为止,结果看起来和单高斯模型一样。
但是混合高斯模型的概率密度函数是多个单高斯分布的加权和。
对于式 (2) 下划线的部分,混合高斯模型的情况为:
p
(
x
)
‾
=
∑
k
=
1
K
w
k
g
k
(
x
∣
μ
k
,
Σ
k
)
(3)
\underline{p(\textbf{\textit{x}})} = \sum^{K}_{k=1} w_k g_k \left( \textbf{\textit{x}} | \boldsymbol{\mu}_k , \Sigma_k \right)\tag{3}
p(x)=k=1∑Kwkgk(x∣μk,Σk)(3)
K
K
K 表示单高斯模型的数量;
g
k
g_k
gk 是第
k
k
k 个单高斯模型的概率密度函数,其均值为
μ
k
\boldsymbol{\mu}_k
μk,协方差为
Σ
k
\Sigma_k
Σk。
w
k
w_k
wk 是第
k
k
k 个单高斯模型的权值。(这里
w
k
w_k
wk 都是
1
/
K
1/K
1/K)
把式 (3) 代入式 (2) 后,得到这个样子: μ ^ , Σ ^ = arg max μ , Σ ∑ i = 1 N ln { 1 K ∑ k = 1 K g k ( x i ∣ μ k , Σ k ) } (4) \hat{\boldsymbol{\mu}},\hat{\Sigma} = \arg \max_{\boldsymbol{\mu},\Sigma} \sum^N_{i=1} \ln \left\{ \frac{1}{K} \sum^{K}_{k=1} g_k \left( \textbf{\textit{x}}_i | \boldsymbol{\mu}_k , \Sigma_k \right) \right\}\tag{4} μ^,Σ^=argμ,Σmaxi=1∑Nln{K1k=1∑Kgk(xi∣μk,Σk)}(4)其中 g k ( x ) = 1 ( 2 π ) D / 2 ∣ Σ k ∣ 1 / 2 exp { − 1 2 ( x − μ k ) T Σ k − 1 ( x − μ k ) } g_k(\textbf{\textit{x}})= \frac{1}{(2\pi)^{D/2} |\Sigma_k|^{1/2}} \exp\left\{ -\frac{1}{2}(\textbf{\textit{x}} - \boldsymbol{\mu}_k)^T \Sigma_k^{-1}(\textbf{\textit{x}} - \boldsymbol{\mu}_k)\right\} gk(x)=(2π)D/2∣Σk∣1/21exp{−21(x−μk)TΣk−1(x−μk)}
到这里再也无法解析地化简这个方程了,因为对数方程里出现了高斯模型的和。这意味着我们只能通过迭代计算,估计参数。并且得到的结果可能不是全局最优的。
Don’t worry,下面来简要地讲讲EM算法,这种算法能在特定条件下给出合理的解。
期望最大化算法(Expectation-Maximization,EM)
这里先用一个特例来看一下如何计算混合高斯模型的参数,然后扩展到 EM 的一般情况。
首先来看看用 EM 算法计算 GMM 的参数,我们需要两个新的东西:
∙
\bullet
∙ 初始的
μ
\boldsymbol{\mu}
μ 和
Σ
\Sigma
Σ。(可以通过猜测来设定)
∙
\bullet
∙ 潜变量
z
z
z(Latent variable)。
利用一个猜测值去求解我们要估计的参数,听起来可能有点奇怪。
在这类称为非凸的优化问题里,有许多次优解,称为局部最小值,而初始的猜测将影响我们找到的解。猜测初始值也是一门技术活,本教程暂时不研究这个。
现在来看看潜变量:
z
k
i
=
g
k
(
x
i
∣
μ
k
,
Σ
k
)
∑
k
=
1
K
(
x
i
∣
μ
k
,
Σ
k
)
z^{i}_k = \frac{g_k \left( \textbf{\textit{x}}_i | \boldsymbol{\mu}_k,\Sigma_k \right)} {\sum^K_{k=1} \left( \textbf{\textit{x}}_i | \boldsymbol{\mu}_k,\Sigma_k \right)}
zki=∑k=1K(xi∣μk,Σk)gk(xi∣μk,Σk)
z k i z^{i}_k zki 是关于第 k k k 个高斯模型的第 i i i 个数据点的潜变量。它定义为第 k k k 个高斯模型在该数据点处的相对比,本质上表征了第 i i i 个数据点由第 k k k 个高斯模型生成的概率。
(每个数据点上都有 k k k 个潜变量,所以共有 k × i k\times i k×i 个潜变量)
来看一个 K K K 等于 2 2 2 时的 1 1 1 维高斯模型的例子:
我们有数据点
x
i
\textbf{\textit{x}}_i
xi,还有两个高斯模型
g
1
g_1
g1 和
g
2
g_2
g2。
如果:
g
1
g_1
g1 产生
x
i
\textbf{\textit{x}}_i
xi 的概率是
p
1
p_1
p1;
g
2
g_2
g2 产生
x
i
\textbf{\textit{x}}_i
xi 的概率是
p
2
p_2
p2;
则在数据点 x i \textbf{\textit{x}}_i xi 处,潜变量 z 1 i = p 1 p 1 + p 2 z^i_1 = \dfrac{p_1}{p_1+p_2} z1i=p1+p2p1, z 2 i = p 2 p 1 + p 2 z^i_2 = \dfrac{p_2}{p_1+p_2} z2i=p1+p2p2
这个例子中 p 1 p_1 p1 比 p 2 p_2 p2 大,所以 x i \textbf{\textit{x}}_i xi 更可能由 g 1 g_1 g1 产生。
现在给定所有数据点关于所有模型的这个潜变量,我们重新定义带权重的均值和协方差矩阵。以
z
z
z 为权值,如果一个数据点属于第
k
k
k 个高斯模型的概率很小,则这个点对计算第
k
k
k 个高斯模型的参数的贡献就很小。
μ
^
k
=
1
z
k
∑
i
=
1
N
z
k
i
x
i
Σ
^
k
=
1
z
k
∑
i
=
1
N
z
k
i
(
x
i
−
μ
^
)
(
x
i
−
μ
^
)
T
z
k
=
∑
i
=
1
N
z
k
i
\begin{aligned} \hat{\boldsymbol{\mu}}_k &= \frac{1}{z_k} \sum^N_{i=1}z^i_{k} \textbf{\textit{x}}_i \\ \hat{\Sigma}_k &= \frac{1}{z_k} \sum^N_{i=1}z^i_{k} ( \textbf{\textit{x}}_i - \hat{\boldsymbol{\mu}} ) ( \textbf{\textit{x}}_i - \hat{\boldsymbol{\mu}} )^T \\ z_k &= \sum^N_{i=1}z^i_{k} \end{aligned}
μ^kΣ^kzk=zk1i=1∑Nzkixi=zk1i=1∑Nzki(xi−μ^)(xi−μ^)T=i=1∑Nzki
这个求解的形式和多维高斯分布很像,只是多了一个权重项
z
z
z,可以回顾一下前面作为比较。
综合之前讨论的东西,我们可以迭代计算模型参数和潜变量。
首先,我们初始化参数
μ
\boldsymbol{\mu}
μ 和
Σ
\Sigma
Σ:
然后用初始的
μ
\boldsymbol{\mu}
μ 和
Σ
\Sigma
Σ 来计算潜变量
z
z
z:
当
z
z
z 更新后,我们就固定
z
z
z,再用它来更新
μ
\boldsymbol{\mu}
μ 和
Σ
\Sigma
Σ:
往复地进行这两个步骤,直到参数的变化量变得非常小,得到估计值
μ
^
k
\hat{\boldsymbol{\mu}}_k
μ^k 和
Σ
^
k
\hat{\Sigma}_k
Σ^k:
我们将不会在这里证明这个算法,但是迭代最终会收敛到一个局部最优值上。
至此已经介绍了如何利用迭代的 EM 算法来估计 GMM 的参数,通过一个初始设定的均值和协方差矩阵,及引入潜变量,我们神奇地将问题变成了在局部最优意义下可解的问题。
下面介绍更一般化的 EM 算法。
回顾上面公式 (4),这个极大似然估计不好求解,求和里面又有对数,对数里面又有求和。
如果我们能找到似然函数的下边界,并且使下边界最大化,这就相当于将似然函数最大化。
现在要做的是,把 EM 算法看作是对目标函数下界的最大化过程。我们需要找到一个表达式,来描述目标函数的下界。
为了找到这个表达式,这里用到了 Jensen不等式 。它是这个算法的基础,EM算法的主要思想可以用Jensen不等式解释。
Jensen不等式
考虑图所示的一个凸函数,在定义域中取两个点 x 1 x_1 x1 和 x 2 x_2 x2,以及他们的中点 x 1 + x 2 2 \dfrac{x_1 + x_2}{2} 2x1+x2。
现在来比较两个值,首先是中点的函数值 f ( x 1 + x 2 2 ) \textcolor{e05d89}{f(\dfrac{x_1 + x_2}{2})} f(2x1+x2)(红色的点),它是曲线上的一个点。
第二个是两点函数值的平均值 f ( x 1 ) + f ( x 2 ) 2 \textcolor{02a0da}{\dfrac{f(x_1) + f(x_2)}{2}} 2f(x1)+f(x2)(蓝色的点),它落在 f ( x 1 ) f(x_1) f(x1) 和 f ( x 2 ) f(x_2) f(x2) 连线的直线上。
容易看到在这个区间里,直线上的点永远大于曲线上的点。于是有图中那条不等式: f ( x 1 + x 2 2 ) ≤ f ( x 1 ) + f ( x 2 ) 2 f(\dfrac{x_1 + x_2}{2}) \leq \dfrac{f(x_1) + f(x_2)}{2} f(2x1+x2)≤2f(x1)+f(x2)
我们把这个想法一般化,下面这个不等式在任意大于 0 0 0 的权值下都成立: f ( a 1 x 1 + a 2 x 2 ) ≤ a 1 f ( x 1 ) + a 2 f ( x 2 ) a 1 + a 2 = 1 a 1 ≥ 0 a 2 ≥ 0 f(a_1 x_1 + a_2 x_2) \leq a_1 f(x_1) + a_2 f(x_2) \\[1.2em] \begin{aligned} a_1+a_2 & = 1 \\ a_1 & \geq 0 \\ a_2 & \geq 0 \end{aligned} f(a1x1+a2x2)≤a1f(x1)+a2f(x2)a1+a2a1a2=1≥0≥0
这是
2
2
2 个点的情况,事实上对于多个点,权值都为正的情况下同样适用。
甚至是在多维的情况下,只要我们的函数是凸的:
f
(
∑
a
i
x
i
)
≤
∑
a
i
f
(
x
i
)
∑
a
i
=
1
a
i
≥
0
f\left( \sum a_i x_i \right) \leq \sum a_i f(x_i)\\[1.2em] \begin{aligned} \sum a_i &= 1 \\ a_i & \geq 0 \end{aligned}
f(∑aixi)≤∑aif(xi)∑aiai=1≥0
如果函数是凹的,我们只需要翻转这个不等式。
我们将要处理的函数是对数函数,它总是凹的:
ln
(
∑
a
i
p
i
)
≥
∑
a
i
ln
p
i
∑
a
i
=
1
a
i
≥
0
\textcolor{blue}{\ln \left( \sum a_i p_i \right)} \textcolor{red}{\geq} \textcolor{green}{\sum a_i \ln p_i} \\[1.2em] \begin{aligned} \sum a_i &= 1 \\ a_i & \geq 0 \end{aligned}
ln(∑aipi)≥∑ailnpi∑aiai=1≥0
(这里要正确区分凹凸函数的概念,请查询相关资料。为了方便解释仅从形状上描述凹凸性)
利用Jensen不等式,我们能获得目标函数的下界。例如这里 ln ( ∑ a i p i ) \textcolor{blue}{\ln \left( \sum a_i p_i \right)} ln(∑aipi) 的下界就是 ∑ a i ln p i \textcolor{green}{\sum a_i \ln p_i} ∑ailnpi。
我们引入我们永远也无法确切知道的隐变量,隐变量的含义在一些特定的应用中是清楚的,例如前面提到的混合高斯模型的参数估计。
我们将对隐变量取边缘概率,根据边缘概率(marginal probability)的定义: p ( X ∣ θ ) = ∑ Z p ( X , Z ∣ θ ) p(X|\theta) = \sum_Z p(X,Z|\theta) p(X∣θ)=Z∑p(X,Z∣θ) (可以理解为 X X X 是样本集合, θ \theta θ 是所有参数集合?)
像之前一样,我们考虑概率的对数: ln p ( X ∣ θ ) = ln ∑ Z p ( X , Z ∣ θ ) \ln p(X|\theta) = \ln \sum_Z p(X,Z|\theta) lnp(X∣θ)=lnZ∑p(X,Z∣θ)
利用Jensen不等式,我们能获得这样的下界: ln ∑ Z p ( X , Z ∣ θ ) = ln ∑ Z q ( Z ) p ( X , Z ∣ θ ) q ( Z ) ≥ ∑ Z q ( Z ) ln p ( X , Z ∣ θ ) q ( Z ) \ln \sum_Z p(X,Z|\theta) = \textcolor{blue}{\ln \sum_Z q(Z) \frac{ p(X,Z|\theta)}{q(Z)}} \geq \textcolor{green}{\sum_Z q(Z)\ln \frac{ p(X,Z|\theta)}{q(Z)}} lnZ∑p(X,Z∣θ)=lnZ∑q(Z)q(Z)p(X,Z∣θ)≥Z∑q(Z)lnq(Z)p(X,Z∣θ) 注: q ( Z ) q(Z) q(Z) 是 Z Z Z 上的有效概率分布。
我们无法直接最大化概率和的对数,但是我们可以利用函数的下界,对于它我们是可以计算的。
获得下界意味着找到了隐变量 z z z 的一个分布 q ( Z ) q(Z) q(Z),这个下界又恰好是一个关于旧的参数猜测的函数。
现在来看一下怎样具体实施这套算法:
首先有一个初始的猜测 θ 0 \theta_0 θ0,并且通过利用猜测,得到的 q q q,进而得到了下界 G G G。
如我说过的,对数形式的概率可以被看作是参数的函数。
假设 F F F 和 G G G 是如图所示的 θ \theta θ 的函数, F F F 和 G G G 并不完全一致 ,但是他们在 θ 0 \theta_0 θ0 附近有一些局部的相似性。
记得我们的任务是,希望求一个参数 θ \theta θ 的最大似然估计。
给定当前的下界 G G G,我们可以找到一个更好的参数,记作 θ ∗ \theta^{\ast} θ∗:
这仍然不是 F F F 的最大值,但它比之前的估计更优。
由于我们有了一个更好的参数的猜测,我们可以计算改进后的下界 G G G:
同样的,新的 G G G 也不那么像 F F F,但它在 θ 1 \theta_1 θ1 附近和 F F F 很接近。
我们又可以根据改进后的 G G G 再来计算一个更好的 θ \theta θ。
如果我们重复这个过程,参数最终会收敛到一个局部最优值:
总结一下,在 E 步骤,我们通过给定的参数计算对数概率的下界:
再次说明,这意味着我们获得了在给定参数 ( θ ^ \hat{\theta} θ^) 上的隐变量的分布 ( q ( Z ) q(Z) q(Z))
在 M 步骤,我们在当前的下界上最大化参数:
重复这个过程,直到产生的变化小于我们的阈值。
我们已经学习了如何利用 EM 算法寻找复杂问题的最大似然解 。
前面课程中提到的混合高斯模型的参数估计,是应用这个算法的一个例子。
在后面探讨定位(localization)时你将看到EM算法的另一个应用。