EM算法实现异常：方差快速收缩至零的排查求助

嘿，我来帮你拆解下这个问题——你遇到的方差快速归零、收敛到均值的情况，本质是GMM EM算法常见的数值稳定性问题，你的代码逻辑大体是对的，但缺了一些关键的稳定化处理，先给你一步步梳理：

核心问题是什么？ #

你发现改高斯概率的指数系数从-0.5到-0.3就好用，其实是因为你人为降低了马氏距离对概率的惩罚力度：原本离均值远的数据点概率会骤降，现在降得慢了，权重不会过度集中在中心数据点上，协方差自然不会被拉到零。但这是个临时 workaround，因为你相当于修改了高斯分布的定义，拟合出的不是标准GMM，所以会出现方差高估的情况。

真正的问题出在代码的数值稳定性和EM步骤的细节处理上，我帮你从代码里找了几个关键点：

1. 协方差矩阵容易变奇异（趋近于零） #

当某个高斯分量的均值刚好靠近某个数据点时，这个点的权重会变得极大，导致协方差被单个点主导，快速收缩到零。而你的代码里：

• 没有给协方差加正则化项，当协方差对角元趋近于0时，矩阵会变成奇异矩阵，求逆和行列式计算都会出现数值爆炸，进一步让概率计算异常，形成“方差越小→权重越集中→方差更小”的恶性循环。
• 看你的calculateProbabilityFast方法，当协方差逆矩阵的元素极大时，马氏距离会变得非常大，exp(-0.5*exponent)会趋近于无穷大，这个分量的权重直接拉满，其他分量被压制，最后就只剩一个均值等于数据整体均值、方差趋近于0的分量。

2. 初始化的小瑕疵 #

你用random(-100,100)初始化均值，如果你的数据范围远小于这个区间，很可能出现某个分量一开始就远离所有数据，权重快速降到几乎为0，剩下的分量容易过度拟合数据。

3. 概率计算的数值溢出风险 #

直接计算exp(-0.5*exponent)很容易出现数值溢出：当马氏距离太大时，指数部分会变成极小的负数，exp结果趋近于无穷大，导致权重计算异常。

具体修复方案 #

1. 给协方差加正则化（最关键） #

在MStep更新完协方差后，给对角元加一个很小的正数，保证矩阵始终正定，不会奇异：

// 在协方差归一化之后、更新逆矩阵之前添加
float regularization = 1e-3; // 可以根据数据范围调整，比如1e-4到1e-2之间
for (int j = 0; j < data[0].nDims; ++j) {
    covariances[Gauss][j][j] += regularization;
}
// 再更新逆矩阵
invCov[Gauss] = inverse(covariances[Gauss]);

这一步能彻底解决方差收缩到零的问题，因为它强制协方差矩阵不能太小。

2. 改用对数概率计算，避免数值溢出 #

把概率计算改成对数形式，再转换为概率，能大幅提升数值稳定性：

float calculateProbabilityFast(DataPoint x, int Gauss) {
    // 先计算对数形式的归一化项
    float logNormalizer = (dims/2.0f) * log(TWO_PI) + 0.5f * log(det(covariances[Gauss]));
    // 计算马氏距离平方
    float mahalanobisSq = 0;
    for (int i = 0; i < x.nDims; ++i) {
        float temp = 0;
        for (int j = 0; j < x.nDims; ++j) {
            temp += (x.data[j] - means[Gauss][j]) * invCov[Gauss][i][j];
        }
        mahalanobisSq += temp * (x.data[i] - means[Gauss][i]);
    }
    // 对数概率
    float logProb = -0.5f * mahalanobisSq - logNormalizer;
    // 转换为概率，防止溢出（如果logProb太大，exp会溢出，所以加个上限）
    float prob = exp(logProb);
    if (prob > 1e10f) prob = 1e10f;
    return prob;
}

对数计算能避免大数/小数的溢出问题，尤其是协方差很小时的情况。

3. 优化初始化，贴近数据分布 #

• 均值不要随机生成，而是从数据里随机采样，这样初始位置更合理：

// 替换原来的均值初始化代码
int randomDataIdx = (int) random(data.length); // 假设初始化时可访问数据数组
means[i] = data[randomDataIdx].data.clone();

• 权重初始化为均匀分布，比随机初始化更稳定：

weights[i] = 1.0f / gNo;

4. 加入收敛判断，避免过度迭代 #

当参数变化小于某个阈值时停止迭代，防止数值漂移：
比如记录上一轮的均值和协方差，计算它们的平方差之和，当小于1e-6时就停止迭代。

验证效果 #

做完这些修改后，你可以把指数系数改回原来的-0.5，应该就能正常收敛，不会出现方差归零的问题了，而且拟合出来的是标准的GMM模型，不会有方差高估的情况。

内容的提问来源于stack exchange，提问作者Nabushika