# example: 2d normal distribution
mu = np.array([0., 0.]) # mean
cov = np.array([[1, 0.8], [0.8 , 1]]) # cov. matrix
# Cholesky decomposition, A = lower triangular matrix
A = np.linalg.cholesky(cov)
# number of random vectors
n = 1000
# random numbers following standard normal distr.
# Ndim rows, N columns (here Ndim = 2)
z = np.random.normal(size=(2, n))
# matrix with Ndim rows (here Ndim = 2):
# [[mu1, mu1, ...], [mu2, mu2, ...]]
mean = np.repeat(mu, n).reshape(2,n)
# matrix of random vectors
# vector i: (rand[1][i], rand[2][i])
rand = mean + np.dot(A, z)