17

我使用 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)
4

1 回答 1

24

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
于 2013-05-22T18:51:42.667 回答