There are several tricks available for sampling from a multivariate Gaussian without matrix inversion--cholesky/LU decomposition among them. Are there any tricks for calculating the likelihood of a multivariate Gaussian without doing the full matrix inversion?

I'm working in python, using numpy arrays. scipy.stats.multivariate_normal is absurdly slow for the task, taking significantly longer than just doing the matrix inversion directly with numpy.linalg.inv.

So at this point I'm trying to understand what is best practice.