4

我正在尝试进行适当的 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 分解？

4

3 回答 3

4

再试一次，用

>>> scipy.__version__
'0.11.0'
>>> np.__version__
'1.6.2'

它完美地工作：

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

输出是：

[[ 11.22274635   5.10611692   0.70263037   3.14603088] # V before
 [  5.10611692   8.94518939  -3.17865941   1.64689675]
 [  0.70263037  -3.17865941   7.35385131  -2.23948391]
 [  3.14603088   1.64689675  -2.23948391   8.25112653]]
[[ 11.22274635   5.10611692   0.70263037   3.14603088] # R before
 [  5.10611692   8.94518939  -3.17865941   1.64689675]
 [  0.70263037  -3.17865941   7.35385131  -2.23948391]
 [  3.14603088   1.64689675  -2.23948391   8.25112653]]
[[ 11.22274635   5.10611692   0.70263037   3.14603088] # V after
 [  5.10611692   8.94518939  -3.17865941   1.64689675]
 [  0.70263037  -3.17865941   7.35385131  -2.23948391]
 [  3.14603088   1.64689675  -2.23948391   8.25112653]]
[[ 3.35003677  1.52419728  0.20973811  0.93910339] # R after
 [ 0.          2.57332704 -1.35946252  0.08375069]
 [ 0.          0.          2.33703292 -0.99382158]
 [ 0.          0.          0.          2.52478036]]
于 2013-01-18T22:56:29.760 回答
2

为了完整起见，以下基于 Warren 注释的代码解决了该问题：

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('F')
print V.flags['C_CONTIGUOUS']
print R.flags['F_CONTIGUOUS']
scipy.linalg.cholesky(R,False,overwrite_a=True)
print V
print R

我所做的只是将 R 的底层存储格式从 C 顺序更改为 Fortran 顺序，这会导致overwrite_a按预期运行。

如果您的矩阵是真正有价值的，您可以通过简单地将转置传递给 cholesky 来执行此操作而无需复制。如果您的矩阵是复值矩阵，我认为转置修复将不起作用，因为复数厄米矩阵的转置通常不等于其共轭转置。但是，如果您确保您的矩阵首先使用 fortran order ( R.flags['F_CONTIGUOUS'] == True)，那么您应该没有问题。但是，请参阅在 Numpy 中设置 *default* 数据​​顺序（C 与 Fortran）以了解该策略的困难。

于 2013-01-22T22:03:51.800 回答
2

如果您足够勇敢，您可以避免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由此完成的复制也可能更有效。

于 2013-01-19T12:02:18.530 回答