使用Distributions.jl拟合MvNormal报PosDefException如何解决

报错原因说明 #

你遇到的报错本质和输入矩阵非方阵无关，是fit(MvNormal)内部计算得到的样本协方差矩阵非正定导致的：光谱数据的特征维度（波数数量，对应矩阵行数）通常远大于样本数量（对应矩阵列数），此时计算得到的经验协方差矩阵必然不满秩，存在零特征值甚至数值计算引入的负特征值，无法通过Cholesky分解校验。

一、调整数据/参数满足正定要求 #

你可以通过以下3种常用方案解决：

• 协方差加岭修正
  手动计算均值和协方差后，给协方差矩阵对角元素添加极小的正数偏移，强制矩阵正定：

  using LinearAlgebra, Distributions
  
  μ = vec(mean(myMatrix, dims=2)) # 计算各波数的均值
  Σ = cov(myMatrix, dims=2) # 计算样本协方差矩阵
  # 给对角加偏移，1e-6可根据你的数据量级调整，保证大于最小负特征值的绝对值即可
  Σ_reg = Σ + 1e-6 * I 
  mvndist = MvNormal(μ, Σ_reg)
  

• 光谱降维预处理
  光谱数据本身冗余度极高，先用PCA将特征维度降至小于样本数，既可以解决协方差非正定问题，还能过滤噪声，是光谱场景下更推荐的方案：

  using MultivariateStats, Distributions
  
  # 拟合PCA模型，保留99%方差，且输出维度不超过样本数-1
  pca = fit(PCA, myMatrix, maxoutdim=size(myMatrix,2)-1, pratio=0.99)
  reduced_data = predict(pca, myMatrix) # 得到降维后的矩阵
  mvndist = fit(MvNormal, reduced_data)
  
  # 后续如果需要生成原维度的样本，可通过PCA重建
  # gen_sample_reduced = rand(mvndist)
  # gen_sample_original = reconstruct(pca, gen_sample_reduced)
  

• 使用收缩协方差估计
  替换默认的经验协方差为带正则化的收缩协方差估计，输出结果天然满足正定要求：

  using CovarianceEstimation, Distributions
  
  μ = vec(mean(myMatrix, dims=2))
  # 使用带常数收缩的线性估计，自动保证协方差正定
  Σ = cov(LinearShrinkage(ConstantShrinkage()), myMatrix, dims=2)
  mvndist = MvNormal(μ, Σ)
  

二、替代拟合方案 #

如果你不想手动调整数据或协方差，可直接使用内置正则参数的拟合接口：

• 直接调用fit_mle的正定向容参数
  Distributions.jl的fit_mle方法内置了posdef_tol参数，会自动给协方差矩阵添加最小必要偏移满足正定要求，无需手动计算：

  using Distributions
  mvndist = fit_mle(MvNormal, myMatrix, posdef_tol=1e-6)
  

• 贝叶斯估计方案
  如果需要更鲁棒的高维正态分布拟合，可使用Turing、Soss等概率编程包，给协方差矩阵添加LJK等正定先验，自动规避Cholesky分解失败问题，适合样本量极小的场景。

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