【问题标题】:Numpy Cholesky decomposition LinAlgErrorNumpy Cholesky 分解 LinAlgError
【发布时间】:2014-03-03 11:49:08
【问题描述】:

在我尝试对周期性边界条件的二维数组的方差-协方差矩阵执行 cholesky 分解时,在某些参数组合下,我总是得到LinAlgError: Matrix is not positive definite - Cholesky decomposition cannot be computed。不确定是numpy.linalg 还是实现问题,因为脚本很简单:

sigma = 3.
U = 4

def FromListToGrid(l_):
    i = np.floor(l_/U)
    j = l_ - i*U
    return np.array((i,j))

Ulist = range(U**2)

Cov = []
for l in Ulist:
    di = np.array([np.abs(FromListToGrid(l)[0]-FromListToGrid(i)[0]) for i, x in enumerate(Ulist)])
    di = np.minimum(di, U-di)

    dj = np.array([np.abs(FromListToGrid(l)[1]-FromListToGrid(i)[1]) for i, x in enumerate(Ulist)])
    dj = np.minimum(dj, U-dj)

    d = np.sqrt(di**2+dj**2)
    Cov.append(np.exp(-d/sigma))
Cov = np.vstack(Cov)

W = np.linalg.cholesky(Cov)

尝试移除潜在奇点也未能解决问题。非常感谢任何帮助。

【问题讨论】:

  • 你做了什么来消除奇点?像 Cov = Cov + numpy.diag(numpy.repeat(delta, k)) 这样的东西有用吗? (基本上是给Cov加了一个小的对角矩阵。这里delta是一个小浮点数,k是Cov的维度)
  • 我只是有 Cov = Cov + d*np.identity(k)。但是查看原始矩阵,似乎没有任何值接近于零..

标签: python numpy scipy linear-algebra covariance


【解决方案1】:

深入挖掘问题,我尝试打印 Cov 矩阵的特征值。

print np.linalg.eigvalsh(Cov)

答案竟然是这样的

[-0.0801339  -0.0801339   0.12653595  0.12653595  0.12653595  0.12653595 0.14847999  0.36269785  0.36269785  0.36269785  0.36269785  1.09439988 1.09439988  1.09439988  1.09439988  9.6772531 ]

啊哈!注意到前两个负特征值了吗?现在,一个矩阵是正定的当且仅当它的所有特征值都是正的。因此,矩阵的问题不在于它接近“零”,而在于它是“负数”。为了扩展 @duffymo 的类比,这是线性代数,相当于尝试取负数的平方根。

现在,让我们尝试执行相同的操作,但这次使用 scipy。

scipy.linalg.cholesky(Cov, lower=True)

这说明不了更多的事情

numpy.linalg.linalg.LinAlgError: 12-th leading minor not positive definite

这说明了更多的事情,(虽然我真的不明白为什么它在抱怨 12 次未成年人)。

底线,矩阵不是很接近“零”,而是更像“负”

【讨论】:

【解决方案2】:

问题在于您提供给它的数据。根据求解器,矩阵是奇异的。这意味着对角元素为零或接近零,因此不可能进行反转。

如果您能提供一个小版本的矩阵,诊断会更容易。

零对角线并不是创建奇点的唯一方法。如果两行彼此成比例,则解决方案中不需要两者;他们是多余的。它比只在对角线上寻找零更复杂。

如果你的矩阵是正确的,你有一个非空的空空间。您需要将算法更改为 SVD 之类的东西。

请参阅下面的评论。

【讨论】:

  • 嗯。当调用上述矩阵的np.diagonal(Cov) 时,它会输出一个1 的数组。另外,由上述脚本生成的 16x16 矩阵是我找到的返回错误消息的最小矩阵,但可能还是太大而无法粘贴到这里?
  • 我不知道。您了解错误告诉您的数学意义,对吗?它是除以零的线性代数。它与您的矩阵有关,与 NumPy 或您的编码无关。我建议您需要仔细检查以确保正确填充该矩阵。如果您确定这一点,但仍然出现错误,我会说您应该将算法更改为类似 SVD 的东西,如果您告诉它如何处理,它将处理奇异矩阵。
猜你喜欢
  • 2013-02-12
  • 2020-05-29
  • 2017-11-04
  • 2020-11-05
  • 2015-06-20
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2018-11-27
相关资源
最近更新 更多