python - cholesky在numpy和scipy之间有什么区别？

利用Cholesky分解对多维高斯随机变量进行采样，计算出随机变量的功率谱。我从numpy.linalg.cholesky得到的结果在高频下总是比从scipy.linalg.cholesky得到的更高的功率。
这两个函数之间可能导致这个结果的区别是什么？哪一个在数值上更稳定？
这是我使用的代码：

n = 2000

m = 10000

c0 = np.exp(-.05*np.arange(n))

C = linalg.toeplitz(c0)

Xn = np.dot(np.random.randn(m,n),np.linalg.cholesky(C))

Xs = np.dot(np.random.randn(m,n),linalg.cholesky(C))

Xnf = np.fft.fft(Xn)

Xsf = np.fft.fft(Xs)

Xnp = np.mean(Xnf*Xnf.conj(),axis=0)

Xsp = np.mean(Xsf*Xsf.conj(),axis=0)

最佳答案

默认情况下，scipy.linalg.cholesky为您提供上三角分解，而np.linalg.cholesky为您提供下三角分解。从docs forscipy.linalg.cholesky：

cholesky(a, lower=False, overwrite_a=False)
    Compute the Cholesky decomposition of a matrix.

    Returns the Cholesky decomposition, :math:`A = L L^*` or
    :math:`A = U^* U` of a Hermitian positive-definite matrix A.

    Parameters
    ----------
    a : ndarray, shape (M, M)
        Matrix to be decomposed
    lower : bool
        Whether to compute the upper or lower triangular Cholesky
        factorization.  Default is upper-triangular.
    overwrite_a : bool
        Whether to overwrite data in `a` (may improve performance).

例如：

>>> scipy.linalg.cholesky([[1,2], [1,9]])
array([[ 1.        ,  2.        ],
       [ 0.        ,  2.23606798]])
>>> scipy.linalg.cholesky([[1,2], [1,9]], lower=True)
array([[ 1.        ,  0.        ],
       [ 1.        ,  2.82842712]])
>>> np.linalg.cholesky([[1,2], [1,9]])
array([[ 1.        ,  0.        ],
       [ 1.        ,  2.82842712]])

如果我修改了您的代码以同时使用相同的随机矩阵并改为使用linalg.cholesky(C,lower=True)，那么我会得到如下答案：

>>> Xnp
array([ 79621.02629287+0.j,  78060.96077912+0.j,  77110.92428806+0.j, ...,
        75526.55192199+0.j,  77110.92428806+0.j,  78060.96077912+0.j])
>>> Xsp
array([ 79621.02629287+0.j,  78060.96077912+0.j,  77110.92428806+0.j, ...,
        75526.55192199+0.j,  77110.92428806+0.j,  78060.96077912+0.j])
>>> np.allclose(Xnp, Xsp)
True

关于python - cholesky在numpy和scipy之间有什么区别？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/16699163/