TL; DR：

你的数据中有高度相关的特征，这会导致数值错误。有几种可能的方法来解决这个问题，每一种都有优点和缺点。下面提出了gaussian_kde的直接替换类。

诊断

您的dataset（即在创建gaussian_kde对象时，* 而不是在使用它时，您所输入的矩阵可能包含高度依赖的特征。这一事实（可能与低数值分辨率（如float32）和“太多”数据点相结合）导致协方差矩阵具有接近零的特征值，由于数值误差可能会低于零。这反过来又会破坏在数据集协方差矩阵上使用Cholesky分解的代码（有关具体细节，请参见解释）。
假设你的数据集的形状是(dims, N)，你可以通过np.linalg.eigh(np.cov(dataset))[0] <= 0测试这是否是你的问题。如果任何输出出来True，让我成为第一个欢迎你加入俱乐部。

疗程

这个想法是让所有的特征值大于零。
1.将数值分辨率提高到实际的最高浮点数可能是一个简单的解决方案，值得尝试，但可能还不够。
1.考虑到这是由相关的 * 特征 * 引起的，删除数据点并没有太大的 * 先验 * 帮助。有一个渺茫的希望，有更少的数字粉碎将传播更少的错误，并保持特征值大于零。它很容易实现，但它会丢弃数据点。
1.一个更复杂的修复方法是识别高度相关的特征并合并它们或忽略“多余”的特征。这是棘手的，尤其是当维度之间的相关性因示例而异时。
1.也许最干净的方法是保持数据不变，并将有问题的特征值提升为正值。这通常通过两种方式完成：

1. SVD直接解决了其核心问题：分解协方差矩阵，并将负特征值替换为小的正epsilon。这将使您的矩阵回到PD域，引入最小的误差。
   1.如果SVD太昂贵而无法计算，则另一种数值方法是将epsilon * np.eye(D)添加到协方差矩阵。这样做的效果是将epsilon加到每个特征值上。同样，如果epsilon足够小，这不会引入太多的错误。
   如果你采用最后一种方法，你需要告诉gaussian_kde修改它的协方差矩阵。这是我发现的一种相对简单的方法：只需将这个类添加到您的代码库中，并将gaussian_kde替换为GaussianKde（在我这边测试过，似乎工作正常）。

class GaussianKde(gaussian_kde):
        """
        Drop-in replacement for gaussian_kde that adds the class attribute EPSILON
        to the covmat eigenvalues, to prevent exceptions due to numerical error.
        """
    
        EPSILON = 1e-10  # adjust this at will
    
        def _compute_covariance(self):
            """Computes the covariance matrix for each Gaussian kernel using
            covariance_factor().
            """
            self.factor = self.covariance_factor()
            # Cache covariance and inverse covariance of the data
            if not hasattr(self, '_data_inv_cov'):
                self._data_covariance = np.atleast_2d(np.cov(self.dataset, rowvar=1,
                                                             bias=False,
                                                             aweights=self.weights))
                # we're going the easy way here
                self._data_covariance += self.EPSILON * np.eye(
                    len(self._data_covariance))
                self._data_inv_cov = np.linalg.inv(self._data_covariance)
    
            self.covariance = self._data_covariance * self.factor**2
            self.inv_cov = self._data_inv_cov / self.factor**2
            L = np.linalg.cholesky(self.covariance * 2 * np.pi)
            self._norm_factor = 2*np.log(np.diag(L)).sum()  # needed for scipy 1.5.2
            self.log_det = 2*np.log(np.diag(L)).sum()  # changed var name on 1.6.2

说明

如果你的错误是类似的，但不完全是，或者有人感到好奇，这里是我遵循的过程，希望它能有所帮助。
1.异常堆栈指定错误是在Cholesky分解过程中触发的。在我的例子中，这是_compute_covariance方法中的这一行。
1.事实上，在异常之后，通过np.eigh检查covariance和inv_cov属性的特征值，发现covariance具有接近于零的负特征值，因此它的逆特征值很大。由于Cholesky期望PD矩阵，因此这触发了错误。
1.在这一点上，我们可以严重怀疑微小的负特征值是一个数值误差，因为协方差矩阵是PSD。实际上，当协方差矩阵最初从传递给构造函数的数据集计算时，误差源就出现了。在我的例子中，协方差矩阵产生了至少2个几乎相同的列，这导致了由于数值误差而产生的剩余负特征值。
1.什么时候你的数据集会导致一个准奇异协方差矩阵？这个问题在另一个SE post中得到了完美的解决。底线是：If some variable is an exact linear combination of the other variables, with constant term allowed, the correlation and covariance matrices of the variables will be singular. The dependency observed in such matrix between its columns is actually that same dependency as the dependency between the variables in the data observed after the variables have been centered (their means brought to 0) or standardized (if we mean correlation rather than covariance matrix)（荣誉和+1到ttnphns的惊人工作）。

Scipy 1.11.2的验证码

# SciPy imports.
from scipy import linalg, special
from numpy import (atleast_2d,cov,pi)

class GaussianKde(stats.gaussian_kde):
    """
    Drop-in replacement for gaussian_kde that adds the class attribute EPSILON
    to the covmat eigenvalues, to prevent exceptions due to numerical error.
    """

    EPSILON = 1e-10 # adjust this at will

    def _compute_covariance(self):
        """Computes the covariance matrix for each Gaussian kernel using
        covariance_factor().
        """
        self.factor = self.covariance_factor()
        # Cache covariance and Cholesky decomp of covariance
        if not hasattr(self, '_data_cho_cov'):
            self._data_covariance = atleast_2d(cov(self.dataset, rowvar=1,
                                               bias=False,
                                               aweights=self.weights))
            
            self._data_covariance += self.EPSILON * np.eye(len(self._data_covariance))
            
            self._data_cho_cov = linalg.cholesky(self._data_covariance,
                                                 lower=True)

        self.covariance = self._data_covariance * self.factor**2
        self.cho_cov = (self._data_cho_cov * self.factor).astype(np.float64)
        self.log_det = 2*np.log(np.diag(self.cho_cov
                                        * np.sqrt(2*pi))).sum()