I use Cholesky decomposition to sample random variables from multi-dimension Gaussian, and calculate the power spectrum of the random variables. The result I get from numpy.linalg.cholesky always has higher power in high frequencies than from scipy.linalg.cholesky.

What are the differences between these two functions that could possibly cause this result? Which one is more numerically stable?

Here is the code I use:

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)
python numpy scipy
Edit
Report