使用scipy计算多元正态CDF时遭遇正定矩阵错误的问题

问题重现 #

当使用scipy.stats.multivariate_normal计算多元正态分布CDF时，即使协方差是对角元全正的对角矩阵（理论上属于正定矩阵），但当其中某个对角元极小时，会触发「input matrix must be symmetric positive definite」错误。比如这段代码：

import numpy as np
from scipy.stats import multivariate_normal
std = np.array([0.00001, 2])
mean = np.array([1.23, 3])
multivariate_normal(mean=mean, cov=np.diag(std**2)).cdf([2,1])

会直接报错，但把std改成[0.001, 2]就能正常返回结果。高维矩阵下这个问题会更突出，只要存在极小的对角元就容易触发错误。

原因分析 #

理论上对角元全正的对角矩阵确实是正定矩阵，但scipy在处理多元正态分布时，内部会对协方差矩阵做Cholesky分解来完成后续计算。当对角元极小（比如1e-10量级）时，浮点计算的精度误差会让分解过程认为矩阵接近奇异，或者出现数值上的非正定情况，进而抛出错误。维度越高，数值误差的累积效应越明显，可接受的最小对角元阈值也就越高。

解决方案 #

这里有几个实用的解决办法：

1. 给协方差矩阵加微小正则化项 #

给对角矩阵的每个对角元加一个极小的正数（比如1e-12），确保数值上的正定，避免浮点误差导致的误判：

import numpy as np
from scipy.stats import multivariate_normal
std = np.array([0.00001, 2])
mean = np.array([1.23, 3])
cov = np.diag(std**2) + 1e-12 * np.eye(len(std))  # 加正则化项
multivariate_normal(mean=mean, cov=cov).cdf([2,1])

2. 利用独立变量的CDF可分离性 #

因为对角协方差矩阵意味着各个变量是独立的，多元正态的联合CDF等于每个单变量正态CDF的乘积。这种方法完全绕开协方差矩阵的数值问题，计算效率也更高：

import numpy as np
from scipy.stats import norm
std = np.array([0.00001, 2])
mean = np.array([1.23, 3])
# 逐个计算单变量CDF再相乘
cdf_value = np.prod(norm.cdf([2, 1], loc=mean, scale=std))

3. 使用更宽松的CDF计算函数 #

可以尝试scipy.stats.mvn.mvnun函数，它对协方差矩阵的正定检查更灵活，适合处理数值上接近奇异的正定矩阵：

import numpy as np
from scipy.stats import mvn
std = np.array([0.00001, 2])
mean = np.array([1.23, 3])
cov = np.diag(std**2)
# mvnun需要传入下界、上界、均值、协方差
lower = np.array([-np.inf, -np.inf])
upper = np.array([2, 1])
cdf_value, _ = mvn.mvnun(lower, upper, mean, cov)

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