Today I want to introduce a very interesting paper of Bhattacharya et al[1]. It is a short miscellanea about how to improve the speed and alleviate the time complexity of the sampling of a multivariate normal random vector. But of course, it doesn’t apply generally to every existing form of normal distributions but rather we should use a special structure of the covariance matrix and the mean vector that arises very often in the context of Bayesian statistics.
The paper is really succinct and contains everything that’s necessary (obviously) so I don’t think I can add more explanation to it, nor do I think it is possible to further summarize it. So I’ll just reiterate what’s explained in the paper.
Most of the times, Bayesian inference resorts to the Markov chain Monte Carlo (MCMC), especially Gibbs sampler if there is no special difficulty. Because of the product structure of the likelihood and the prior, i.e.,
the full conditional distribution of becomes as follows:
where is a symmetric positive definite matrix, , and . Normally, we use the Cholesky decomposition to sample from such distribution, i.e.
but the paper points out that the computation of the Cholesky factor has complexity so it quickly becomes a huge burden for the computer as grows larger. So they suggest the following algorithm. I will rewrite what’s written in the paper.

Sample and independently.

Set .

Solve to get .

Set .
The proof that this is equivalent to the one that’s obtained by the Cholesky decomposition can be easily shown by using the ShermanMorrisonWoodbury matrix identity.
The paper also offers an empirical study of how much faster the algorithm becomes compared to the Cholesky decomposition. The value of this new algorithm shines when the regression problem is highdimensional so it is natural that the paper also mentions the highdimensional setting where the priors are aimed to single out the signals from the noises, horse shoe prior in particular.
library(mvnfast) Bhattacharya < function(X,D,alpha,diag_flag=TRUE) { p < nrow(D) n < nrow(X) if (diag_flag) { u < rnorm(p)*sqrt(diag(D)) delta < rnorm(n) v < drop(X%*%u)+delta U < tcrossprod(D,X) w < solve(X%*%U+diag(n),alphav) return(u+drop(U%*%w)) } else { u < drop(rmvn(1,numeric(p),D)) delta < rnorm(n) v < drop(X%*%u)+delta U < tcrossprod(D,X) w < solve(X%*%U+diag(n),alphav) return(u+drop(U%*%w)) } }