python - 如何在 Python 中进行 Cholesky 分解

我正在尝试进行就地 Cholesky 分解，但是该功能实际上并未在 scipy 中实现，或者有一些我不理解的地方。我在这里发帖以防万一是后者。这是我正在做的事情的简化示例:

import numpy
import scipy.linalg

numpy.random.seed(0)
X = numpy.random.normal(size=(10,4))
V = numpy.dot(X.transpose(),X)
R = V.copy()
scipy.linalg.cholesky(R,False,overwrite_a=True)
print V
print R

我认为应该发生的是 R 被上三角矩阵覆盖。但是，当我运行此代码时，我的 V 和 R 是相同的(cholesky 没有覆盖 R)。我误解了 overwrite_a 的目的，还是犯了其他错误？如果这只是 scipy 中的一个错误，是否有解决方法或其他方法可以在 Python 中进行就地 Cholesky 分解？

最佳答案

如果你足够勇敢，你可以避免 scipy 并低级调用 linalg.lapack_lite.dpotrf

import numpy as np

# prepare test data
A = np.random.normal(size=(10,10))
A = np.dot(A,A.T)
L = np.tril(A)

# actual in-place cholesky
assert L.dtype is np.dtype(np.float64)
assert L.flags['C_CONTIGUOUS']
n, m = L.shape
assert n == m
result = np.linalg.lapack_lite.dpotrf('U', n, L, n, 0)
assert result['info'] is 0

# check if L is the desired L cholesky factor
assert np.allclose(np.dot(L,L.T), A)
assert np.allclose(L, np.linalg.cholesky(A))

您必须了解 lapack 的 DPOTRF 、fortran 调用约定、内存布局。不要忘记在退出时检查 result['info'] == 0 。尽管如此，您会看到它只是一行代码，通过丢弃 linalg.cholesky 完成的所有健全性检查和复制，这也可能更有效。

关于python - 如何在 Python 中进行 Cholesky 分解，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/14408873/