机器学习数学推导（十三） EM 算法与 GMM

EM 算法(Expectation-Maximization)是处理隐变量模型的通用框架——当数据含有未观测的隐变量时,直接最大化似然困难,EM 通过迭代"期望"和"最大化"两步,保证似然单调上升直至收敛。从高斯混合模型的参数估计到图像分割、语音识别,EM 算法展现了强大的理论美感与实用价值。Jensen 不等式保证了算法的收敛性，ELBO（证据下界）提供了优化目标，而 GMM 则是 EM 算法最经典的应用场景。

隐变量与不完全数据

问题背景

完全数据：观测数据 = {_1, , _N} $及隐变量$

不完全数据：只观测，隐变量未知

目标：最大化不完全数据对数似然

困难：对数内有求和（或积分），难以直接优化。

示例：高斯混合模型

观测：数据点
隐变量：所属成分
似然：

对数内的求和使直接求导困难。

Jensen 不等式与对数似然下界

Jensen 不等式：对于凹函数$f $和随机变量$ X $：$ $

对于（凹函数）：

应用到对数似然：引入任意分布：

ELBO (Evidence Lower Bound)：是对数似然的下界。

KL 散度与紧界条件

改写 ELBO：

其中 KL 散度：

紧界条件：当时，，下界等于对数似然。

EM 算法推导

算法框架

初始化：参数

循环（第次迭代）：

E 步（ Expectation）：固定，选择使下界紧：

计算 Q 函数：

M 步（Maximization）：固定，优化参数：

终止：收敛（对数似然或参数变化小于阈值）

Q 函数的意义

解释： Q 函数是完全数据对数似然在隐变量后验分布下的期望。

收敛性证明

定理：EM 算法保证对数似然单调递增：

证明：

其中是熵。第一个不等号由 KL 散度非负,第二个由 M 步的定义。

推论:EM 算法收敛到对数似然的(局部)极大值或鞍点。

高斯混合模型(GMM)

模型定义

生成过程:

选择成分，其中是成分的先验概率
从成分生成数据:

联合分布:

边缘分布(观测数据的似然):

参数：

约束：,

EM 算法 for GMM

E 步：计算责任度 (responsibility)：

物理意义:样本属于成分的后验概率(软分配)。

M 步:更新参数。

Q 函数:

更新（用拉格朗日乘数法处理约束）:

其中是成分的有效样本数。

更新：对求导：

加权平均：样本按责任度加权。

更新：

加权协方差。

GMM 的几何解释

K-means vs GMM: - K-means:硬分配( $Extra close brace or missing open bracez_i \in \{1, \dots, K\\}$ ),球形簇 - GMM:软分配（），椭球簇，概率输出

GMM 的优势: - 允许重叠簇 - 不同簇不同形状(协方差) - 提供不确定性估计

初始化策略

EM 算法对初始值敏感。常用策略:

K-means 初始化:用 K-means 结果初始化, ,
随机初始化:从数据中随机选择个点作为
多次重启:运行多次 EM,选择对数似然最大的结果

GMM 的应用

密度估计

目标:学习数据分布G MM 提供灵活的密度模型:任意连续分布都可用足够多高斯成分逼近(万能逼近性)。

应用:异常检测,生成新样本,概率预测。

聚类分析

硬聚类:

软聚类:保留作为隶属度

BIC 模型选择:选择成分数:

其中是参数数量:。

图像分割

问题:将图像像素聚类为个区域。

特征:RGB 颜色,位置坐标,纹理特征。

GMM 模型:每个成分对应一个区域。

分割结果:像素分配到。

完整实现

import numpy as np
from scipy.stats import multivariate_normal

class GMM:
    def __init__(self, n_components=3, max_iter=100, tol=1e-4, init='kmeans'):
        self.n_components = n_components
        self.max_iter = max_iter
        self.tol = tol
        self.init = init
        
        self.weights = None  # pi_k
        self.means = None  # mu_k
        self.covariances = None  # Sigma_k
        self.converged = False
        
    def _initialize(self, X):
        N, d = X.shape
        K = self.n_components
        
        if self.init == 'kmeans':
            # K-means 初始化
            from sklearn.cluster import KMeans
            kmeans = KMeans(n_clusters=K, random_state=0).fit(X)
            self.means = kmeans.cluster_centers_
            labels = kmeans.labels_
            
            self.weights = np.array([np.sum(labels == k) / N for k in range(K)])
            self.covariances = np.array([
                np.cov(X[labels == k].T) + 1e-6 * np.eye(d) for k in range(K)
            ])
        elif self.init == 'random':
            # 随机初始化
            indices = np.random.choice(N, K, replace=False)
            self.means = X[indices]
            self.weights = np.ones(K) / K
            self.covariances = np.array([np.cov(X.T) + 1e-6 * np.eye(d) for _ in range(K)])
        
    def _e_step(self, X):
        """E 步:计算责任度"""
        N = X.shape[0]
        K = self.n_components
        gamma = np.zeros((N, K))
        
        for k in range(K):
            gamma[:, k] = self.weights[k] * multivariate_normal.pdf(
                X, mean=self.means[k], cov=self.covariances[k]
            )
        
        # 归一化
        gamma_sum = gamma.sum(axis=1, keepdims=True)
        gamma /= gamma_sum + 1e-10
        
        return gamma
    
    def _m_step(self, X, gamma):
        """M 步:更新参数"""
        N, d = X.shape
        K = self.n_components
        
        # 有效样本数
        N_k = gamma.sum(axis=0)
        
        # 更新权重
        self.weights = N_k / N
        
        # 更新均值
        self.means = (gamma.T @ X) / N_k[:, np.newaxis]
        
        # 更新协方差
        for k in range(K):
            diff = X - self.means[k]
            self.covariances[k] = (gamma[:, k, np.newaxis] * diff).T @ diff / N_k[k]
            self.covariances[k] += 1e-6 * np.eye(d)  # 正则化
        
    def _compute_log_likelihood(self, X):
        """计算对数似然"""
        N = X.shape[0]
        K = self.n_components
        
        log_likelihood = 0
        for i in range(N):
            prob = 0
            for k in range(K):
                prob += self.weights[k] * multivariate_normal.pdf(
                    X[i], mean=self.means[k], cov=self.covariances[k]
                )
            log_likelihood += np.log(prob + 1e-10)
        
        return log_likelihood
    
    def fit(self, X):
        """训练 GMM"""
        self._initialize(X)
        
        log_likelihood_old = -np.inf
        
        for iteration in range(self.max_iter):
            # E 步
            gamma = self._e_step(X)
            
            # M 步
            self._m_step(X, gamma)
            
            # 检查收敛
            log_likelihood = self._compute_log_likelihood(X)
            
            if abs(log_likelihood - log_likelihood_old) < self.tol:
                self.converged = True
                print(f"收敛于第{iteration + 1}次迭代")
                break
            
            log_likelihood_old = log_likelihood
        
        return self
    
    def predict(self, X):
        """预测类别(硬分配)"""
        gamma = self._e_step(X)
        return np.argmax(gamma, axis=1)
    
    def predict_proba(self, X):
        """预测概率(软分配)"""
        return self._e_step(X)
    
    def score(self, X):
        """计算对数似然"""
        return self._compute_log_likelihood(X)
    
    def sample(self, n_samples=100):
        """生成样本"""
        # 采样成分
        components = np.random.choice(
            self.n_components, size=n_samples, p=self.weights
        )
        
        # 从每个成分采样
        samples = np.zeros((n_samples, self.means.shape[1]))
        for k in range(self.n_components):
            mask = components == k
            n_k = np.sum(mask)
            if n_k > 0:
                samples[mask] = np.random.multivariate_normal(
                    self.means[k], self.covariances[k], size=n_k
                )
        
        return samples, components

# 示例使用
if __name__ == '__main__':
    import matplotlib.pyplot as plt
    from sklearn.datasets import make_blobs
    
    # 生成数据
    X, y_true = make_blobs(n_samples=500, centers=3, n_features=2, 
                           cluster_std=0.6, random_state=42)
    
    # 训练 GMM
    gmm = GMM(n_components=3, max_iter=100)
    gmm.fit(X)
    
    # 预测
    y_pred = gmm.predict(X)
    gamma = gmm.predict_proba(X)
    
    # 可视化
    plt.figure(figsize=(15, 5))
    
    plt.subplot(131)
    plt.scatter(X[:, 0], X[:, 1], c=y_true, cmap='viridis')
    plt.title('真实标签')
    
    plt.subplot(132)
    plt.scatter(X[:, 0], X[:, 1], c=y_pred, cmap='viridis')
    plt.scatter(gmm.means[:, 0], gmm.means[:, 1], 
                c='red', marker='x', s=200, linewidths=3)
    plt.title('GMM 聚类')
    
    plt.subplot(133)
    plt.scatter(X[:, 0], X[:, 1], c=gamma[:, 0], cmap='coolwarm')
    plt.title('成分 1 的责任度')
    plt.colorbar()
    
    plt.tight_layout()
    plt.savefig('gmm_clustering.png', dpi=150)
    plt.show()
    
    print(f"对数似然: {gmm.score(X):.2f}")
    print(f"权重: {gmm.weights}")
    
    # 生成新样本
    X_new, _ = gmm.sample(n_samples=100)
    plt.figure()
    plt.scatter(X_new[:, 0], X_new[:, 1], alpha=0.6)
    plt.scatter(gmm.means[:, 0], gmm.means[:, 1], 
                c='red', marker='x', s=200, linewidths=3)
    plt.title('GMM 生成的样本')
    plt.savefig('gmm_samples.png', dpi=150)
    plt.show()

Q&A 精选

Q1:EM 算法为什么叫期望最大化?

A:E 步计算 Q 函数(完全数据对数似然的期望),M 步最大化Q 函数。名字直接描述算法步骤。

Q2:EM 算法会收敛到全局最优吗?

A:不保证。 EM 保证对数似然单调上升,但可能收敛到局部最优或鞍点。解决方法: - 多次随机初始化 - 更好的初始化(K-means) - 结合其他优化算法(如梯度下降)

Q3:如何选择 GMM 的成分数 K?

A:模型选择准则: - BIC: - AIC: - 交叉验证:最大化 held-out 对数似然 - 轮廓系数:聚类质量评估

通常选择 BIC 最大的。

Q4:GMM 与 K-means 的关系?

A:K-means 是 GMM 的特殊情况: - 所有协方差（球形簇） - 时,GMM 退化为 K-means（硬分配）

GMM 是 K-means 的软化和推广。

Q5:EM 算法能并行化吗?

A:E 步天然并行（每个样本独立计算），M 步部分并行（统计量计算可并行，参数更新需同步）。分布式实现:Map-Reduce 或参数服务器。

Q6:协方差矩阵奇异怎么办?

A:常见原因:数据维度高,某成分样本少。解决方法: - 正则化: - 对角协方差: - 共享协方差:所有成分用相同 - PCA 降维

Q7:EM 算法的变种有哪些?

A: - 增量 EM:在线学习,逐批更新 - 随机 EM:每次迭代随机采样子集 - 变分 EM:E 步用变分推断近似后验 - 广义 EM:M 步不完全优化(如梯度上升几步)

Q8:GMM 能处理高维数据吗?

A:高维困难(维度灾难): - 协方差参数爆炸 - 样本稀疏,估计不准

解决方法: - 对角协方差(假设特征独立) - 因子分析(低秩协方差) - 特征选择/降维

Q9:如何可视化高维 GMM?

A: - PCA 降维:投影到前 2-3 个主成分 - t-SNE:非线性降维 - 责任度热图:样本×成分的矩阵 - 成分权重:柱状图显示 ---

Q10:EM 算法在其他模型中的应用?

A: - 隐马尔可夫模型:Baum-Welch 算法 - 主题模型:LDA 的变分 EM - 缺失数据插补:EM 迭代估计缺失值 - 混合专家模型:门控网络+专家网络

EM 是处理隐变量的通用框架。

✏️ 练习题与解答

练习 1：E步计算

题目：GMM两分量，，，。观测，计算。解答：

练习 2：M步更新

题目：3个样本，值。更新。解答：

练习 3：EM收敛性

题目：EM算法保证收敛到什么？解答：局部最优（或鞍点）。EM保证似然单调不减，但不保证全局最优。多次随机初始化可提高找到好解的概率。

练习 4：GMM vs K-means

题目：GMM和K-means区别？解答：K-means硬分配，GMM软分配（概率）。GMM建模簇形状（协方差），K-means假设球形簇。

练习 5：缺失数据场景

题目：数据，缺失。如何用EM估计？解答：E步：用当前参数估计（填充期望）。M步：用完整+填充数据更新参数。迭代至收敛。

参考文献

Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B, 39(1), 1-22.
Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer. [Chapter 9: Mixture Models and EM]
Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective. MIT Press. [Chapter 11: Mixture Models and EM]
McLachlan, G., & Krishnan, T. (2007). The EM Algorithm and Extensions (2nd ed.). Wiley.
Neal, R. M., & Hinton, G. E. (1998). A view of the EM algorithm that justifies incremental, sparse, and other variants. Learning in Graphical Models, 89, 355-368.

EM 算法以其优雅的数学结构和广泛的实用价值,成为现代机器学习的基石之一。从高斯混合模型的参数估计到隐马尔可夫模型的训练,从缺失数据处理到半监督学习,EM 算法展现了"化繁为简"的智慧——将困难的不完全数据优化问题分解为简单的"期望"和"最大化"两步。理解 EM 算法,不仅是掌握经典算法,更是领悟概率推断与优化的深刻联系——这一思想贯穿变分推断、蒙特卡洛方法等现代贝叶斯机器学习技术。