努比的cholesky和scipy有什么区别？

Question

利用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)

Answer 1

默认情况下，scipy.linalg.cholesky给出的是上三角分解，而np.linalg.cholesky给出的是下三角分解。从scipy.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