2013-05-22 16 views
15

Uso la scomposizione di Cholesky per campionare variabili casuali da gaussiane a più dimensioni e calcolare lo spettro di potenza delle variabili casuali. Il risultato che ottengo da numpy.linalg.cholesky ha sempre una maggiore potenza in alta frequenza rispetto a scipy.linalg.cholesky.Qual è la differenza tra cholesky in numpy e scipy?

Quali sono le differenze tra queste due funzioni che potrebbero causare questo risultato? Quale è più numericamente stabile?

ecco il codice che uso:

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

Dalla scipy faq [Qual è la differenza tra NumPy e SciPy?] (Http://new.scipy.org/faq.html#what-is-the-difference-between-numpy-and-scipy) : "In ogni caso, SciPy co dispone di versioni più complete dei moduli di algebra lineare e di molti altri algoritmi numerici. " Vedi anche [Perché entrambi 'numpy.linalg' e' scipy.linalg'? Qual è la differenza?] (Http://new.scipy.org/faq.html#why-both-numpy-linalg-and-scipy-linalg-what-s-the-difference). –

risposta

19

scipy.linalg.cholesky ti dà la decomposizione triangolare superiore per impostazione predefinita, mentre np.linalg.cholesky ti dà la versione a basso triangolare. Dalla documentazione per 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). 

Ad esempio:

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

Se modifico il codice per utilizzare la stessa matrice casuale entrambe le volte e da usare linalg.cholesky(C,lower=True) invece, poi ho ottenere risposte come:

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

è possibile calcolare il triangolare superiore con la funzione cholesky di numpy? Il documento ufficiale non sembra aiutare con questa domanda (https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.cholesky.html). –

Problemi correlati