2013-05-22 13 views
15

Ich benutze Cholesky-Zerlegung, um Zufallsvariablen aus Multi-Dimension-Gaussian zu erfassen und das Leistungsspektrum der Zufallsvariablen zu berechnen. Das Ergebnis, das ich von numpy.linalg.cholesky bekomme, hat immer eine höhere Leistung in hohen Frequenzen als von scipy.linalg.cholesky.Was ist der Unterschied zwischen Cholesky in numpy und scipy?

Welche Unterschiede zwischen diesen beiden Funktionen könnten möglicherweise zu diesem Ergebnis führen? Welcher ist numerisch stabiler? Hier

ist der Code, den ich verwenden:

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) 
+0

Aus der scipy faq [Was ist der Unterschied zwischen NumPy und SciPy?] (Http://new.scipy.org/faq.html#what-ist-the-difference-between-numpy-and-scipy) : "SciPy co Es enthält mehr Funktionen der linearen Algebra-Module sowie viele andere numerische Algorithmen. " Siehe auch [Warum numpy.linalg' und 'scipy.linalg'? Was ist der Unterschied?] (Http://new.scipy.org/faq.html#why-both-numpy-linalg-and-scipy-linalg-what-s-the-difference). –

Antwort

19

scipy.linalg.cholesky gibt Ihnen die obere Dreieckszerlegung standardmäßig während np.linalg.cholesky gibt Ihnen die untere Dreiecks Version. Aus der Dokumentation für 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). 

Zum Beispiel:

>>> 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]]) 

Wenn ich Ihren Code modifizieren, um die gleiche Zufallsmatrix beide Male zu verwenden und linalg.cholesky(C,lower=True) stattdessen zu verwenden, dann bekomme ich Antworten wie:

>>> 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 
+0

ist es möglich, das obere Dreieck mit der cholesky-Funktion von numpy zu berechnen? Das offizielle Dokument scheint bei dieser Frage nicht hilfreich zu sein (https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.cholesky.html). –

Verwandte Themen