使用Python计算多元观测概率时遇奇异矩阵错误排查

首先，咱们来拆解一下你遇到的问题：你用10个10维的随机样本计算协方差矩阵，然后用这个矩阵去计算多元正态分布的概率密度，结果触发了奇异矩阵的错误。这可不是代码逻辑写错了，而是样本数量和特征维度的关系导致的必然问题。

为什么会出现奇异矩阵？ #

样本协方差矩阵的秩有个关键限制：它的秩不会超过 min(样本数-1, 特征数)。你这里样本数是10，特征数也是10，所以协方差矩阵的最大秩只能是9——也就是说，这个矩阵天生就是奇异的（不可逆）。而scipy.stats.multivariate_normal默认要求协方差矩阵是非奇异的，直接运行自然会报错。

几个可行的解决办法 #

1. 增加样本数量（最推荐）

让样本数远大于特征维度，这样协方差矩阵的秩大概率能达到特征数，自然就非奇异了。修改后的代码示例：

from scipy.stats import multivariate_normal
import numpy as np

def get_mean(x, axis=0):
    return x.mean(axis=axis)

def get_sigma(x):
    return np.cov(x, rowvar=False)

def get_probabilities(x, mu, sigma):
    return multivariate_normal.pdf(x, mean=mu, cov=sigma)

# 把样本数从10改成100，远大于10个特征
x = np.random.rand(100,10)
t = np.random.rand(1, 10)
mu = get_mean(x)
sigma = get_sigma(x)
p = get_probabilities(t, mu, sigma)
print(p)

2. 允许使用奇异协方差矩阵

如果你必须用当前的样本量，可以在调用pdf时设置allow_singular=True，这样函数会处理奇异矩阵的情况（本质是在样本张成的子空间里计算概率）：

def get_probabilities(x, mu, sigma):
    # 添加allow_singular参数，允许奇异协方差矩阵
    return multivariate_normal.pdf(x, mean=mu, cov=sigma, allow_singular=True)

不过要注意：奇异协方差意味着数据存在线性相关的特征，概率密度的计算结果可能不是你预期的全空间分布。

3. 给协方差矩阵加正则项

给协方差矩阵的对角线加一个极小的正数，让它变成非奇异矩阵，避免数值上的奇异问题：

def get_sigma(x):
    cov = np.cov(x, rowvar=False)
    # 添加正则项，比如1e-6的单位矩阵，既保证可逆又不破坏原有结构
    return cov + 1e-6 * np.eye(x.shape[1])

这种方法适合样本量接近特征数的场景，是数值计算中常用的稳定手段。

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