I consider a Gaussian distribution with density
$$p({\bf x}) \propto \exp\left[-\frac{1}{2}({\bf x} - {\boldsymbol\mu})^\top{\boldsymbol\Sigma}^{-1}({\bf x} - {\boldsymbol\mu})\right]$$
where ${\bf x}, {\boldsymbol\mu}\in\mathcal{S}^{K}$ live on a $K$-dimensional hypersphere. Hence, the normalization constant $C({\boldsymbol\mu}, {\boldsymbol\Sigma})$ is given by
$$C({\boldsymbol\mu}, {\boldsymbol\Sigma}) = \int_{\mathcal{S}^{K}} \text{d}{\bf x}\exp\left[-\frac{1}{2}({\bf x} - {\boldsymbol\mu})^\top{\boldsymbol\Sigma}^{-1}({\bf x} - {\boldsymbol\mu})\right].$$
Question: I am trying to find an (efficient) algorithm to sample from this distribution, particularly in high dimensions ($K \approx 500$).
I so far used rejection sampling, but at least a naive version of that doesn't scale to high dimensions. If it helps, $\boldsymbol\Sigma$ can be considered a diagonal matrix. I highly appreciate any hint on how to approach this!