Mixture Models and EM

Book Notes on Pattern Recognition and Machine Learning (PRML)

K-means Clustering

K-means 聚类是一种无监督学习，从数据本身的特征出发，不依靠标签将数据划分为不同类别。K-means 的算法流程很直观，因此我们先介绍 K-means 聚类的实现流程，再从数学角度刻画分析

K-means 算法流程如下：

首先选择 K 值，即最终将数据聚为 K 类。
初始化，在每个类中选择一个数据点将其划分进去，作为该类的中心
对每个数据点，将其划分到距类中心最近的类中
所有点划分完成后，根据划分好的数据更新这 K 个类中心点的位置
重新分配每一个数据点到最近的类中
重复步骤 4 和 5直到收敛

在介绍完算法流程后，我们用数学形式去刻画该方法。我们假设数据集为 \(\{x_1,...,x_n\}\) ，需要将其聚为 k 类。每个类别 k 的中心点记为 \(\mu_k\) ，由此我们可以写出距离公式作为损失函数，即

\[L = \sum_{i=1}^n\sum_{j=1}^k\delta_{ij}||x_i-\mu_j||^2\]

这里的 \(\delta_{ij}\) 为 0-1 函数，如果数据 i 被划分到类别 j，那么 \(\delta_{ij}=1\) ；否则为零。我们的目标就是找到每个数据的划分方式 \(\delta\) ，以及各聚类的中心点 \(\mu\) ，使得loss最小。我们可以通过迭代的方式来实现目标。我们最开始先确定初始的中心点 \(\mu\)

接下来第一步是确定划分方式 \(\delta\) 。此时我们已经有上一步得到的各聚类中心点 \(\mu\) ，那么对于某个数据 \(x_i\) 而言，我们写出划分方式对应的更新公式

\[\delta_{ik} = 1 \quad if\ k=argmin_{j}||x_i-\mu_j||^2 \qquad (else\ \ 0)\]

在迭代的第二步，我们要修改中心点的位置。我们对损失函数求导，可以得到其关于 \(\mu\) 的导数

\[\frac{\partial L}{\partial \mu_j} = 2\sum_{i=1}^n\delta_{ij}(\mu_j-x_i)=0 \Rightarrow \mu_j = \frac{\sum_{i}\delta_{ij}x_i}{\sum_i\delta_{ij}}\]

此外，我们也可以利用在线的随机算法来实现聚类中心 \(\mu\) 的更新。对某个数据 \(x_i\) ，我们更新距其最近的聚类中心 \(\mu_j\) 的值，即

\[\mu_j^{new} = \mu_{j}^{old} + \eta(x_i-\mu_{j}^{old})\]

最后是关于 K-means 算法的一些细节

收敛的标准是数据点不再移动到其它类，同时各类的中心点位置不变。此时对应的Loss可能只是极小值，不一定是全局最小（这与初始聚类中心的选取有关）
算法的复杂度：每一次迭代是 O(KN) ，K 是类别数，N 是数据总数
聚类数 K 的选取：尝试不同的 K ，计算数据点到类中心的平均距离。该距离会急剧下降直到合适的 K 值，然后缓慢变化

Expectation Maximization Algorithm

EM 算法是用来寻找含隐变量的概率模型极大似然解的一种方法。和上面一样，我们还是先给出 EM 算法的流程，再从数学角度刻画分析。

我们记观察变量为 X，隐变量为 Z，模型参数统记为 \(\theta\) 。EM 算法用来求解后验概率 \(p(X\vert\theta)\) 的极大似然，包含 E 和 M 两步。

E-step: 计算 \(Q(\theta\vert\theta^{(t)})\) ，这里

\[Q(\theta\vert\theta^{(t)}) = E_{Z\sim P(Z\vert X,\theta^{(t)})}[\ln P(X,Z\vert \theta)]\]

M-step: 选择使 \(Q(\theta\vert\theta^{(t)})\) 最大的参数，用来更新当前模型参数，即

\[\theta^{(t+1)} = argmax_{\theta}Q(\theta\vert\theta^{(t)})\]

可以看到，这是一个迭代的过程，因此整体的 EM 算法流程可写为

初始化，给出初始的模型参数 \(\theta^{(0)}\)
E-step，在给定模型 \(\theta\) 下，计算不同隐变量值 Z 的概率，进而计算 \(Q(\theta\vert\theta^{(t)})\)
M-step，利用计算的 Q 值更好地修改模型
重复 2 和 3 步骤直至收敛

在介绍完 EM 算法的流程后，很自然我们有两个问题。E-step 中的优化目标 \(Q(\theta\vert\theta^{(t)})\) 是从何而来的呢？EM算法为什么是正确的呢？

我们首先来看第一个问题，E-step中的优化目标从何而来？对此我们有两种解释方式

首先我们想用极大似然估计，用观测数据去确定模型参数，即

\[max_{\theta}\ln P(X\vert \theta) = max_{\theta} \ln \sum_{Z}P(X,Z\vert \theta)\]

这里 Z 作为未知的隐变量，表示 X 中各个数据所属的类别。由于 Z 未知，我们无法直接得到原优化目标里的 \(P(X,Z\vert \theta)\) 。但如果给出 X 和 \(\theta\) ，我们可以求出 Z 对应的后验分布，即 \(P(Z\vert X,\theta)\) 。所以我们先根据现有模型去猜测 Z 的分布，再用这个猜测的 Z 去计算极大似然，从而改善模型参数。这样我们就得到了 E-step中的优化目标，即 \(max_{\theta} E_{Z\sim P(Z\vert X,\theta^{(t)})}[\ln P(X,Z\vert \theta)]\)

第二种解释方式是梯度更新的角度。我们的原目标为极大似然 \(max_{\theta} \ln P(X\vert\theta)\) ，我们对优化目标求导可以得到

\[\Delta\theta = \eta\frac{1}{P(X\vert\theta)}\frac{\partial}{\partial\theta}P(X\vert\theta) = \eta\frac{1}{P(X\vert\theta)}\sum_{Z}\frac{\partial}{\partial\theta}P(X,Z\vert\theta)\] \[\Delta\theta = \eta\frac{1}{P(X\vert\theta)}\sum_{Z}P(X,Z\vert\theta)\frac{\partial}{\partial\theta}\ln P(X,Z\vert\theta)\] \[\Delta\theta = \eta\sum_{Z}P(Z\vert X,\theta)\frac{\partial}{\partial\theta}\ln P(X,Z\vert\theta)\]

在这里我们进行近似，用 \(\theta^{(t)}\) 代替上式中的 \(\theta\) ，就可以得到

\[\Delta\theta = \eta\sum_{Z}P(Z\vert X,\theta^{(t)})\frac{\partial}{\partial\theta}\ln P(X,Z\vert\theta)\]

也就是E-step中的优化目标

\[\Delta\theta = \eta\frac{\partial}{\partial\theta} E_{Z\sim P(Z\vert X,\theta^{(t)})}[\ln P(X,Z\vert \theta)]\]

在解释了 E-step 中优化目标的来源后，下一个问题是为什么 EM 算法是正确的。对此，我们同样有两种证明方式

第一种证明用到了 Jensen 不等式，即对于凸函数 \(f(\cdot)\) 而言，我们有 \(E[f(x)] \ge f(E[x])\) ，而

\[\ln P(X\vert \theta) = \ln \sum_{Z}P(X,Z\vert \theta) = \ln \sum_{Z}P(Z\vert X,\theta^{(t-1)})\frac{P(X,Z\vert \theta)}{P(Z\vert X,\theta^{(t-1)})}\]

其中 \(\sum_{Z}P(Z\vert X,\theta^{(t-1)}) = 1\) ，再结合 Jensen 不等式，我们就可以得到

\[\ln P(X\vert \theta) \ge \sum_{Z}P(Z\vert X,\theta^{(t-1)}) [\ln P(X,Z\vert \theta) - \ln P(Z\vert X,\theta^{(t-1)})]\]

即

\[\ln P(X\vert \theta) \ge Q(\theta\vert\theta^{(t)}) + constant\]

这意味着我们在优化 \(Q(\theta\vert\theta^{(t)})\) 的过程中，实际上是在不断提升原目标的下界，所以相当于也在提升原目标。

第二种证明利用了条件概率，即

\[\forall Z, \qquad \log P(X\vert \theta) = \log P(X,Z\vert\theta) - \log P(Z\vert X,\theta)\]

我们改写原目标可以得到

\[\ln P(X\vert \theta) = [\sum_{Z'}P(Z'\vert X,\theta^{(t)})]\ln P(X\vert \theta)\]

既然条件概率中的 Z 可以是任意值，那么我们就令 Z 的值为 Z’，代入上式可以得到

\[\ln P(X\vert \theta) = [\sum_{Z'}P(Z'\vert X,\theta^{(t)})][\log P(X,Z'\vert\theta) - \log P(Z'\vert X,\theta)] = Q(\theta\vert\theta^{(t)})+H(\theta\vert\theta^{(t)})\]

这里的 \(H(\cdot)\) 类似与信息熵的定义，我们有 \(H(\theta\vert\theta^{(t)}) \ge H(\theta^{(t)}\vert\theta^{(t)})\) ，因此我们就可以得到

\[\ln P(X\vert\theta)-\ln P(X\vert\theta^{(t)}) \geq Q(\theta\vert\theta^{(t)}) - Q(\theta^{(t)}\vert\theta^{(t)})\]

同样说明我们是在不断提升原目标的下界。

EM for Gaussian Mixture

高斯混合分布可以写成不同高斯分布的线性叠加，即

\[p(x) = \sum_{i=1}^K \pi_iN(x\vert \mu_i,\Sigma_i)\]

接下来我们引入隐变量 z，z 是 K 维的0-1向量，向量中只有一个元素为1，其余元素为0，即 \(\sum_{i=1}^Kz_i=1\) 。为了得到和上式相似的形式，我们用混合系数 \(\pi_i\) 来刻画隐变量 z 的分布，可以写出

\[p(z_i=1)=\pi_i, \quad 0\le \pi_i\le 1, \quad \sum_{i=1}^K\pi_i=1\]

现在我们就可以写出隐变量 z 的边缘分布 \(p(z)\)

\[p(z) = \prod_{i=1}^K \pi_i^{z_i}\]

当我们给出一个确定的 z 时，也就是说我们知道了 z 的具体某个元素 \(z_i\) 为一，此时 x 的条件分布同样服从高斯，\(x\vert z_i=1 \sim N(\mu_i,\Sigma_i)\)

我们写出整体的条件分布

\[p(x\vert z) = \prod_{i=1}^KN(x\vert \mu_i,\Sigma_i)^{z_i}\]

在同时得到边缘分布和条件分布后，我们可以写出 x 和 z 的联合分布，进而求出 x 的分布

\[p(x) = \sum_zp(z)p(x\vert z) = \sum_z\prod_{i=1}^K (\pi_iN(x|\mu_i,\Sigma_i))^{z_i} = \sum_{i=1}^K \pi_iN(x\vert \mu_i,\Sigma_i)\]

在引入隐变量 z 后，我们得到了和原高斯分布等价的形式。我们引入隐变量的目的就是将目标从 x 的边缘分布 \(p(x)\) 转换到 x 和 z 的联合分布 \(p(x,z)\) 。

在正式介绍 EM 在高斯混合模型的应用之前，我们还需要引入一个概念，隐变量 z 的后验分布 \(p(z\vert x)\) 。我们用 \(\gamma(z_k)\) 来表示 \(p(z_k=1\vert x)\) ，根据贝叶斯定理可以写出

\[\gamma(z_k) = \frac{p(z_k=1)p(x\vert z_k=1)}{\sum_{i=1}^Kp(z_i=1)p(x\vert z_i=1)} = \frac{\pi_kN(x\vert \mu_k,\Sigma_k)}{\sum_{i=1}^K\pi_iN(x\vert \mu_i,\Sigma_i)}\]

假设我们有数据集 \(\{x_1,...,x_n\}\) ，我们想用高斯混合模型进行建模，那么此时的模型包含三个参数 \(\pi, \mu, \Sigma\) 。考虑用极大似然来优化模型，我们有

\[\ln p(x\vert \pi, \mu, \Sigma) = \sum_{i=1}^n\ln\{\sum_{j=1}^K \pi_jN(x_i\vert \mu_j,\Sigma_j)\}\]

首先将其对 \(\mu_k\) 求导并令其为零，可以得到

\[\sum_{i=1}^n\gamma(z_{ik})\Sigma^{-1}_k(x_i-\mu_k) = 0 \rightarrow \mu_k = \frac{\sum_{i=1}^n\gamma(z_{ik})x_i}{\sum_{i=1}^n\gamma(z_{ik})}\]

其次将其对 \(\Sigma_k\) 求导并令其为零，可以得到

\[\Sigma_k = \frac{\sum_{i=1}^n\gamma(z_{ik})(x_i-\mu_k)(x_i-\mu_k)^T}{\sum_{i=1}^n\gamma(z_{ik})}\]

最后我们需要考虑 \(\sum_{i=1}^K\pi_i=1\) 的限制条件，利用拉格朗日乘子法将问题转化为

\[\ln p(x\vert \pi, \mu, \Sigma)+ \lambda(\sum_{i=1}^K\pi_i-1)\]

对 \(\pi_k\) 求导并令其为零，得到 \(\pi_k\) 的计算公式

\[\lambda + \sum_{i=1}^n\frac{N(x_i\vert \mu_k,\Sigma_k)}{\sum_{j=1}^K\pi_jN(x_i\vert \mu_j,\Sigma_j)}=0 \rightarrow \pi_k = \frac{\sum_{i=1}^n\gamma(z_{ik})}{n}\]

最后我们总结下 EM 算法在高斯混合模型的应用

初始化，确定模型初始值 \(\pi_k, \mu_k, \Sigma_k\)
E-step: 根据现有模型值，计算 \(\gamma\)

\[\gamma(z_{nk}) = \frac{\pi_kN(x_n\vert \mu_k,\Sigma_k)}{\sum_{i=1}^K\pi_iN(x_n\vert \mu_i,\Sigma_i)}\]

M-step: 根据计算的 \(\gamma\) 值重新估计更新模型参数

\[\mu_k = \frac{\sum_{i=1}^n\gamma(z_{ik})x_i}{\sum_{i=1}^n\gamma(z_{ik})}\] \[\Sigma_k = \frac{\sum_{i=1}^n\gamma(z_{ik})(x_i-\mu_k)(x_i-\mu_k)^T}{\sum_{i=1}^n\gamma(z_{ik})}\] \[\pi_k = \frac{\sum_{i=1}^n\gamma(z_{ik})}{n}\]

收敛之后，我们就可以去估计 log-likelihood

\[\ln p(x\vert \pi, \mu, \Sigma) = \sum_{i=1}^n\ln\{\sum_{j=1}^K \pi_jN(x_i\vert \mu_j,\Sigma_j)\}\]

可以发现，在上述介绍中，并没有明显体现出隐变量 z 的作用。在这里我们从隐变量的角度重新看待 EM 算法。在引入隐变量后，我们可以写出联合分布的极大似然

\[p(x,z\vert\pi, \mu,\Sigma) = \prod_{i=1}^n\prod_{j=1}^K\{\pi_jN(x_i\vert \mu_j,\Sigma_j)\}^{z_{ij}}\]

其对应的 log-likelihood 为

\[\ln p(x,z\vert\pi, \mu,\Sigma) = \sum_{i=1}^n\sum_{j=1}^K z_{ij}\{\ln \pi_j + \ln N(x_i\vert \mu_j,\Sigma_j)\}\]

和原问题中的极大似然 \(\ln p(x\vert \pi, \mu, \Sigma)\) 相比，最显著的特点是我们将求和操作从 log 里面拿到了外面，这可以大大简化问题的求解。因为我们引入了隐变量，最后还需要通过期望消除，所以最后的优化目标为

\[E_z[\ln p(x,z\vert\pi, \mu,\Sigma)] = \sum_{i=1}^n\sum_{j=1}^K \gamma(z_{ij})\{\ln \pi_j + \ln N(x_i\vert \mu_j,\Sigma_j)\}\]

该优化目标下参数 \(\pi, \mu, \Sigma\) 的更新方式和上文介绍的一模一样。

Connection between K-means and EM

与 K-means 算法相比，EM 算法达到收敛需要迭代更多次数，并且每次迭代也需要更多的计算量。所以通常的做法是先用 K-means 来确定模型的初始化参数，再用 EM 算法改进提升。

实际上，K-means 算法可以看成是一种特殊极限情况下的 EM 算法结合高斯混合模型。

我们假设所有数据的协方差矩阵由共享的常数 \(\epsilon\) 控制，即 \(\Sigma = \epsilon I\) ，此时 x 在第 i 个高斯模型下的分布可以写出来

\[p(x\vert \mu_i,\Sigma_i) = \frac{1}{(2\pi\epsilon)^{1/2}}exp\{-\frac{1}{2\epsilon}||x-\mu_i||^2\}\]

进一步写出后验概率 \(\gamma(z_{nk})\)

\[\gamma(z_{nk}) = \frac{\pi_kexp\{-||x_n-\mu_k||^2/2\epsilon\}}{\sum_{i=1}^K\pi_iexp\{-||x_n-\mu_i||^2/2\epsilon\}}\]

当 \(\epsilon\) 逐渐减小直至趋近于零时，\(\vert\vert x_n-\mu_i\vert\vert^2\) 最小的那一项最晚降为零，因此只有这一项为一，其余都很快地变为零。即此时 \(\gamma(z_{nk}) \rightarrow \delta_{nk}\) ，每个数据点都被唯一地划分到均值最小的对应聚类中。此时 M-step 中均值的更新退化为

\[\mu_k = \frac{\sum_{i=1}^n\gamma(z_{ik})x_i}{\sum_{i=1}^n\gamma(z_{ik})} \rightarrow \mu_k = \frac{\sum_{i}\delta_{ik}x_i}{\sum_i\delta_{ik}}\]

再写出 \(\epsilon \rightarrow 0\) 下的log-likelihood 的期望

\[E_z[\ln p(x,z\vert\pi, \mu,\Sigma)] = -\frac{1}{2}\sum_{i=1}^n\sum_{j=1}^K \delta_{ij}||x_i-\mu_j||^2 + constant\]

此时最大化上式中的log-likelihood等价于最小化k-means算法中的loss