Scipy 高斯核密度估计：矩阵不是正定的

Question

我正在尝试使用scipy在特定点估计数据集的密度。

from scipy.stats import gaussian_kde
import numpy as np

我有一个3D点的数据集A（这只是一个最简单的例子。我的实际数据具有更多维度和更多样本）。

A = np.array([[0.078377  , 0.76737392, 0.45038174],
       [0.65990129, 0.13154658, 0.30770917],
       [0.46068406, 0.22751313, 0.28122463]])

我想要估算密度的点

，并且需要润色内容使其更加易懂。

B = np.array([[0.40209377, 0.21063273, 0.75885516],
       [0.91709997, 0.79303252, 0.65156937]])

但是我似乎无法使用gaussian_kde函数，因为

result = gaussian_kde(A.T)(B.T)

返回值

LinAlgError: Matrix is not positive definite

如何修复这个错误？如何获取我的样本的密度？

- usernumber

Answer 1

简介：

数据中存在高度相关的特征导致数值误差，有多种方法可以解决，每种方法都有其优缺点。下面提供了一种用于替换gaussian_kde的类。

诊断

您的数据集（即在创建gaussian_kde对象时输入的矩阵，而不是在使用它时输入的矩阵）可能包含高度依赖的特征。这个事实（可能与低数值分辨率（如float32）和“过多”的数据点相结合）导致协方差矩阵具有接近零的特征值，由于数值误差，这些特征值可能会变成负值。这反过来会破坏使用Cholesky分解矩阵的代码，该矩阵是数据集的协方差矩阵（见详细说明）。

假设您的数据集具有形状(dims, N)，您可以通过np.linalg.eigh(np.cov(dataset))[0] <= 0测试是否存在此问题。如果任何输出结果为True，那么欢迎您加入这个俱乐部。

---

处理方式

目标是将所有特征值提高到大于零。

1. 增加数值分辨率到最高的浮点数可能是一个容易的解决方法，值得一试，但可能不够。

2. 考虑到这是由相关的特征引起的，因此在先验情况下删除数据点并没有太大帮助。有一点希望，因为少一些数据点可以减少传播误差，并保持特征值大于零。这很容易实现，但会丢弃数据点。

3. 更复杂的解决方法是识别高度相关的特征并将它们合并或忽略“多余”的特征。这很棘手，特别是如果维度之间的关联在不同实例之间变化。

4. 可能最干净的方法是保持数据不变，并将问题特征值提高到正值。通常有两种方法:

5. SVD直接解决了问题的核心：分解协方差矩阵并用小正数epsilon替换负特征值。这会使您的矩阵返回到PD域，引入最小误差。

6. 如果计算SVD太昂贵，一种替代方法是将epsilon * np.eye(D)添加到协方差矩阵中。这具有将epsilon添加到每个特征值中的效果。同样，如果epsilon足够小，这不会引入太多误差。

如果您选择最后一种方法，您需要告诉gaussian_kde修改其协方差矩阵。我发现这是一个相对简洁的方法：只需将此类添加到您的代码库中，并用GaussianKde替换gaussian_kde（已在我的端口测试，看起来运行良好）。

    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 方法内的这行代码。

2. 确实，在异常之后，通过 np.eigh 检查 covariance 和 inv_cov 属性的特征值表明： covariance 具有接近零的负特征值，因此其逆矩阵具有一个很大的特征值。由于Cholesky期望PD矩阵，这就触发了一个错误。

3. 此时我们可以强烈怀疑微小的负特征值是数值误差，因为协方差矩阵是PSD的。事实上，错误的来源来自协方差矩阵最初从传递给构造函数的数据集计算出来，这里。在我的情况下，协方差矩阵至少产生了 2 个几乎相同的列，这导致了由于数值误差而产生的剩余负特征值。

4. 何时您的数据集会导致一个近似奇异的协方差矩阵？这个问题在另一个SE帖子中得到了完美的解答。底线是：如果一些变量是其他变量的精确线性组合，允许常数项，则变量的相关和协方差矩阵将是奇异的。在这样的矩阵中观察到的依赖关系实际上与在将变量居中后（将它们的平均值变为0）或标准化后（如果我们意味着相关而不是协方差矩阵），在数据中观察到的变量之间的依赖关系相同 （祝贺并+1 给ttnphns的出色工作）。

- fr_andres

Answer 2

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()

- lmbfrank