I am working on using sequential Monte Carlo methods to estimate some kind of financial models. However, this method is really computationally expensive, and the bottleneck boils down to sampling from multivariate Gaussian distribution. Say if the number of samples (particles) that I use is 5000, and the dimension of the problem is 5. It means I will have 5000 different mean vectors with shape (5,), 5000 different covariance arrays with shape (5,5), and I need to sample 5000 samples from those 5000 pairs. How do I do it in a fast way? The fastest solution I have so far is like:
seeds = np.random.multivariate_normal(np.zeros(5), np.eye(5), size = 5000)
Ls = [np.linalg.cholesky(covs[i]) for i in range(5000)]
results = means + np.einsum('nij,njk->nik', Ls, seeds[:,:,np.newaxis])[:,:,0]
where np.einsum is the function in numpy doing einstein sum, and means and covs are the given 5000 pairs that mentioned above.
I think the step of doing cholesky in a for loop is pretty slow, but I couldn't find any way of doing vectorized computation in this case. Any suggestions would be welcome.