2013-05-22 16 views
15

Używam dekompozycji Cholesky'ego do próbkowania zmiennych losowych z wielowymiarowego Gaussa i obliczania widma mocy zmiennych losowych. Wynik jaki uzyskałem z numpy.linalg.cholesky zawsze ma wyższą moc w wysokich częstotliwościach niż od scipy.linalg.cholesky.Jaka jest różnica między cholesky in numpy a scipy?

Jakie są różnice między tymi dwiema funkcjami, które mogą powodować ten wynik? Który z nich jest bardziej stabilny liczbowo?

Oto kod używam:

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

Z scipy faq [Jaka jest różnica między NumPy i SciPy?] (Http://new.scipy.org/faq.html#what-jest-difting-between-numpy-and-scipy) : "W każdym razie, SciPy co zawiera więcej w pełni funkcjonalnych wersji modułów algebry liniowej, a także wiele innych algorytmów numerycznych. " Zobacz także [Dlaczego oba 'numpy.linalg' i' scipy.linalg'? Jaka jest różnica?] (Http://new.scipy.org/faq.html#why-both-numpy-linalg-and-scipy-linalg- what-s-the-difference). –

Odpowiedz

19

scipy.linalg.cholesky daje Ci przewagę trójkątne rozkładu domyślnie, natomiast np.linalg.cholesky daje Ci wersji niższej trójkątny. Od docs dla 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). 

Na przykład:

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

Gdybym zmodyfikować kod, aby korzystać z tej samej macierzy losowych oba razy i używać linalg.cholesky(C,lower=True) zamiast, potem uzyskać odpowiedzi jak:

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

Czy można obliczyć górną trójkątną funkcję chropowatą? Oficjalny dokument nie wydaje się pomóc w tym pytaniu (https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.cholesky.html). –

Powiązane problemy