I want to get random samples from a high-dimensional (say, $10^4$ or more) multivariate Gaussian distribution. Assume a case of zero mean for simplicity, the probability density is given as
$$P({\bf x})=\frac{1}{\sqrt{(2\pi)^n|{\bf \Sigma}|}} \exp(-\frac{1}{2}{\bf x}^T {\bf \Sigma}^{-1} {\bf x})~~~~~~~~~(1)$$
where ${\bf \Sigma}$ is a $n\times n$ covariance matrix.
I here have a low-rank approximation for inverse of the covariance matrix ${\bf Q}={\bf \Sigma}^{-1}$ (i.e., precision matrix, NOT the covariance matrix ${\bf \Sigma}$ itself),
$${\bf Q} \simeq {\bf P}^T {\bf P}$$
where ${\bf P}$ is a $n\times r$ matrix where $r < n$. (Such an approximation can be obtained based on truncated singular value decomposition, for instance)
When $r \ll n$, what is the most efficient way to get the random samples from the Gaussian distribution?
At first, I was thinking that simple Markov chain Monte-Carlo methods like the Metropolis-Hastings (MH) algorithm are directly applicable because $P({\bf x})$ is unimodal and we can relatively easily evaluate $P({\bf x})$ for a given ${\bf x}$ because of fast calculation of ${\bf P}^T {\bf P}$.
However, articles like the one below introduces more complicated sampling algorithms for high-dimensional Gaussian distributions instead of such simple ones. I also want to know why they never tell that application of methods like simple MH algorithm is effective? (Note that the article do not discuss low-rank approximation)
Vono, M., Dobigeon, N., & Chainais, P. (2022). High-dimensional Gaussian sampling: a review and a unifying approach based on a stochastic proximal point algorithm. SIAM Review, 64(1), 3-56.
I am one of the authors of the paper you are citing in your post.
Regarding your question "When ≪, what is the most efficient way to get the random samples from the Gaussian distribution?" :
If it is sufficient for you to target an approximation of your initial Gaussian distribution $P$, then you can obtain approximate samples $\{\theta_t\}_{t=1}^T$ via the Cholesky sampler detailed in Algorithm 3.1. in our paper.
If you want to target $P$, then you can indeed use Metropolis-Hastings algorithm using as a proposal the Gaussian distribution with precision matrix $\mathbf{P}^\top\mathbf{P}$. Note that if your low-rank approximation is poor, the mixing time of your MCMC algorithm might be large.
Regarding your question "I also want to know why they never tell that application of methods like simple MH algorithm is effective?" :
You can indeed use any proposal distribution (e.g. based on a low-rank approximation of your initial Gaussian distribution $P$) within a Metropolis-Hastings algorithm. Nevertheless, both the storage and the computational cost of this procedure can be large. Indeed, you need to store your truncated SVD factor $\mathbf{P}$ which might not be possible on standard laptops if $r$ is not sufficiently small. Also the acceptance ratio can be low if the approximation is not sufficiently good, see for instance https://arxiv.org/pdf/1409.0606.pdf.
In our paper, we only focused on the main research lines that have been considered in the existing literature when Cholesky sampling is too costly (in terms of storage and/or computations). Low-rank approximation approaches encompass Lanczos-based approaches which are discussed in our paper, see for instance Section 3.3.2.