多变量高斯分布、高斯混合模型和EM算法

多变量高斯分布

先总结一些基本结论。

设有随机变量组成的向量\(X=[X_1,\cdots,X_n]^T\)，均值为\(\mu\in\mathbb{R}^n\)，协方差矩阵\(\Sigma\)为对称正定\(n\)阶矩阵。在此基础上，如果还满足概率密度函数

\[p(x;\mu,\Sigma)=\frac{1}{(2\pi)^{\frac{n}{2}}|\Sigma|^{\frac{1}{2}}}\exp\left(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)\right) \]

则称其满足多变量高斯分布，记为\(x\sim\mathcal{N}(\mu,\Sigma)\)。其中\(\exp()\)前的分式可视为归一化系数。

关于协方差矩阵，在前面的PCA和白化里面有总结，一个随机向量的协方差矩阵就是这样一个矩阵：\(\Sigma_{i,j}=\text{Cov}(X_i,X_j)\)。实际上，随机向量的均值和协方差矩阵有下面关系：

\[\Sigma=E[(X-\mu)(X-\mu)^T]=E(XX^T)-\mu\mu^T \]

上面说到，要求\(\Sigma\)是一个对称的正定矩阵，实际上任何随机向量的协方差矩阵都是半正定的，只不过概率密度函数里面进一步要求\(\Sigma^{-1}\)的存在，所以\(\Sigma\)必然是满秩矩阵，所以肯定是正定的。

若\(\Sigma\)是对角形式的（比如随机向量是二维的，高维也类似）

\[\Sigma=\begin{pmatrix}\sigma_1^2&0\\0&\sigma_2^2\end{pmatrix} \]

则化简概率密度函数，可以发现它就是两个标准差分别为\(\sigma_1,\sigma_2\)的单变量高斯分布的乘积。
如果\(\Sigma\)不是对角的，则不能化简成两个单变量高斯分布的乘积。

直觉上来说，应该可以使用正交变换\(R\)使得协方差矩阵\(\Sigma\)对角化（结合协方差矩阵的意义，这会使得随机向量\(X\)每个元素互不相关）。

为此，试探性地作代换\(X\leftarrow R^T(X-\mu)\)，则根据\(\Sigma=E[(X-\mu)(X-\mu)^T]=E(XX^T)-\mu\mu^T\)可得新的协方差矩阵

\[\Sigma'=E[R^T(X-\mu)(X-\mu^T)R]-0=R^TE[(X-\mu)(X-\mu)^T]R=R^T\Sigma R \]

所以取合适的正交变换矩阵\(R\)就可以使得协方差矩阵对角化。一旦协方差矩阵对角化了，则随机向量总的概率密度函数一定可以写成各个变量自己的高斯分布的乘积。

实际上上述对角化的过程，就是对一组随机变量进行线性变换的过程，因为是线性变换，一定是严格单调的，因此可以使用概率论中的雅可比变换。具体详见《随机信号分析》（高新波）第2.6节。

高斯混合模型和EM算法

关于高斯混合模型的背景，一篇博客觉得讲的比较明白（地址），这里就不抄写一些背景概念了，只总结最主要的部分。

样本集（\(x\)不是随机变量，是样本）为\(\{x^{(1)},\cdots,x^{(m)}\}\)，共\(m\)个样本，每个样本是一个向量。为估计出概率密度函数\(p(x;\phi,\mu,\Sigma)\)中的三个参数\(\phi,\mu,\Sigma\)，引入中间变量\(z\)，\(z\)等于几就意味着\(x\)来自于第几个高斯分布。其中，设有\(k\)个类，向量\(\phi\)的每个元素满足\(\phi_j\geqslant0,\sum\limits_{j=1}^k\phi_j=1\)，且\(\phi_j=p(z^{(i)}=j)\)。

这样，给定了样本集\(\{x^{(1)},\cdots,x^{(m)}\}\)，三个待估计参数为\(\phi,\mu,\Sigma\)，为求其估计值，写出对数似然函数：

\[L(\phi,\mu,\Sigma)=\sum^m_{i=1}\log p(x^{(i)};\phi,\mu,\Sigma)=\sum^m_{i=1}\log\sum^k_{z^{(i)}=1}p(x^{(i)}|Z^{(i)};\mu,\Sigma)p(z^{(i)};\phi) \]

直接令其各偏导为零求解，发现并无闭式解。

但是另一方面，注意到，如果对应于\(x^{(i)}\)的中间变量\(z^{(i)}\)的取值已知，则上面对数似然函数就是

\[L(\phi,\mu,\Sigma)=\sum^m_{i=1}\log p(x^{(i)}|Z^{(i)};\mu,\Sigma)+\log p(z^{(i)};\phi) \]

令其各偏导为零，可以得到闭式解（\(I\)为示性函数）

\[\begin{aligned} \phi_j=&\frac{1}{m}\sum^m_{i=1}I(z^{(i)}=j)\\ \mu_j=&\frac{\sum^m\limits_{i=1} I(z^{(i)}=j)x^{(i)}}{\sum^m\limits_{i=1}I(z^{(i)}=j)}\\ \Sigma_j=&\frac{\sum^m\limits_{i=1}I(z^{(i)}=j)(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum^m\limits_{i=1}I(z^{(i)}=j)}\end{aligned} \]

注意，上面的\(\phi\)是一个向量，长度等于混合的高斯模型的个数。而\(\mu,\Sigma\)是矩阵，其每一列\(\mu_j,\Sigma_j\)是属于第\(j\)个高斯模型的均值、标准差向量（因为每个高斯分布都是多维高斯分布）。

而所谓EM算法，直观上就是使用了一个后验概率代替了上式中的示性函数：对每个样本每个高斯分布进行遍历（每个\(i,j\)），计算后验概率\(w_j^{(i)}\)如下（称之为E-step）

\[w^{(i)}_j=p(z^{(i)}=j|x^{(i)};\phi,\mu,\sigma) \]

根据贝叶斯法则，有

\[p(z^{(i)}=j|x^{(i)};\phi,\mu,\sigma)=\frac{p(x^{(i)}|z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{\sum_l p(x^{(i)}|z^{(i)}=l;\mu,\Sigma)p(z^{(i)}=l;\phi)} \]

然后使用该后验概率更新各待估计的参数（称之为M-step）：

\[\begin{aligned} \phi_j=&\frac{1}{m}\sum^m_{i=1}w_j^{(i)}\\\mu_j=&\frac{\sum^m\limits_{i=1}w_j^{(i)}x^{(i)}}{\sum^m\limits_{i=1}w^{(i)}_j}\\ \Sigma_j=&\frac{\sum^m\limits_{i=1}w_j^{(i)}(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum^m\limits_{i=1}w^{(i)}_j}\end{aligned} \]

算法先选取随机的\(\phi,\mu,\Sigma\)代入E-step，然后执行M-step，反复迭代直至收敛。考虑到可能的局部最优，初始值可以多选取几组。

EM算法的数学解释

上面的EM算法确实有了具体步骤，但是它的得到基本上可以说是比较“直观的”，并且也并没有证明它一定收敛。EM意思是“期望-最大化”，是一类算法。它的数学解释（严密推导及其收敛性的证明）需要用到琴生不等式。

琴生（Jensen）不等式说的是，对于随机变量\(X\)和凸函数\(f(\cdot)\)，有\(E[f(X)]\geqslant f[E(X)]\)，且若\(f(\cdot)\)为严格凸函数，则取等号当且仅当随机变量\(X\)是一个常数，即\(E(X)=X\)的概率为\(1\)。对于简单的一元函数\(f(\cdot)\)，若在其定义域内二阶导数\(f''(x)\geqslant0\)恒成立，则它是凸函数，反之亦然；如果严格大于0，则严格凸。如果\(f(\cdot)\)为凹函数，则琴声不等式中的不等号改变方向。

上面高斯混合模型等问题的对数似然函数为

\[L(\theta)=\sum^m_{i=1}\log p(x^{(i)};\theta)=\sum^m_{i=1}\log\sum^m_{z^{(i)}}p(x^{(i)},z^{(i)};\theta) \]

直接将其最大化，往往没有闭式解。

但是可以对每一个样本\(x^{(i)}\)引入一个任意分布\(Q_i(z)\geqslant0\)，满足\(\sum\limits_zQ_i(z)=1,i=1,\cdots,m\)，把对数似然函数变为

\[L(\theta)=\sum^m_{i=1}\log\sum^m_{z^{(i)}}Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \]

由于函数\(\log(x)\)二阶导数为\(-\frac{1}{x^2}<0\)，所以有\(f[E(x)]\geqslant E[f(x)]\)恒成立，所以对任意的一组分布\(Q_i,i=1,\cdots,m\)按照琴生不等式有

\[L(\theta)=\sum^m_{i=1}\log\sum^m_{z^{(i)}}Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\geqslant\sum^m_{i=1} \sum^m_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \]

对于每个猜测的\(\theta\)，上式给出其似然函数的下界，但是也只给出下界；为了离目标更进一步，我们需要选取精心构造的一组分布\(Q_i\)以使得上式取等号，如果能做到的话，那么每给一个猜测的\(\theta\)值，上式就给出其对数似然函数的值。

由于\(\log(x)\)是严格凹函数，所以取等号当且仅当关于\(z^{(i)}\)的随机变量取常数，即：

\[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c \]

\(c\)为不依赖\(z^{(i)}\)的常数。上式即表明\(Q_i(z^{(i)})\)正比于\(p(x^{(i)},z^{(i)};\theta)\)，但是\(Q_i(z^{(i)})\)是归一的，所以必然有

\[Q_i(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta)}{\sum_z p(x^{(i)},z;\theta)}=\frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)}=p(z^{(i)}|x^{(i)};\theta)\tag{$*$} \]

这就给出了满足取等号条件的一组分布\(Q_i\)。也就是说，对于初始猜测的参数\(\theta\)，先计算出一组分布\(Q_i\)，再根据\(\theta\)和\(Q_i\)计算出该\(\theta\)取值对应的似然函数的值，这称为E-step。剩下来的就是令该\(Q_i\)不变，而将对数似然函数关于\(\theta\)最大化（求关于\(\theta\)的偏导，令其为零），而这一步是有闭式解的（称为M-step），得到更新后的\(\theta\)后，再利用它来计算新的一组分布\(Q_i\)，这样反复迭代。

需要注意的是，在将对数似然函数关于\(\theta\)最大化的过程中，把分布\(Q_i\)当做固定不随\(\theta\)变化的，所以EM算法整体上类似于把似然函数看做关于\(Q_i\)和\(\theta\)的坐标上升法。EM算法做的事情，就是把一个难以直接优化的问题，通过引入中间变量(\(z\))及其分布(\(Q\))来转化为一个迭代问题，每次迭代都有解析解。

关于EM算法的收敛性，可以利用琴生不等式证明对数似然函数每次EM迭代后，其值都会有一个正的增量，若对数似然函数有上界，则必然收敛。

在具体的高斯混合模型这里，参数\(\theta\)是\(\phi,\mu,\Sigma\)，\(w^{(i)}_j\)就是\(Q_i(z^{(i)}=j)\)，因此在高斯混合模型那里说在E-step中令（见\(*\)式）

\[w^{(i)}_j=p(z^{(i)}=j|x^{(i)};\phi,\mu,\theta) \]

而在高斯混合模型M-step中的式子，其实就是对固定\(w_j^{(i)}\)（即分布\(Q_i\)）的对数似然函数求偏导令其为零得到的（可验证）。