Sampling from a specific level set of a Gaussian

Ellipsoids

A $(d-1)$-dimensional ellipsoid in ${\mathbb{R}}^d$ is defined by the following equation $$\begin{array}{r}(x-v{)}^{\mathrm{\top }}\mathsf{A}(x-v)=1,\end{array}$$ where $v\in {\mathbb{R}}^n$ is the center of the ellipsoid and $\mathsf{A}$ is an $n\times n$ positive-definite symmetric matrix.

Gaussian Level-Sets as Ellipsoids

Consider a multivariate Gaussian distribution $\pi$ with covariance matrix $\mathrm{\Sigma }$ and mean $\mu$ with density $$\begin{array}{r}\pi (x)=(2\pi {)}^{-d/2}\text{det}(\mathrm{\Sigma }{)}^{-1/2}\exp (-\frac{1}{2}(x-\mu {)}^{\mathrm{\top }}{\mathrm{\Sigma }}^{-1}(x-\mu )).\end{array}$$ The maximum density value achieved is $\pi (\mu )$ as can be seen by setting the gradient to zero $${\mathrm{\nabla }}_x\pi (x)=-\pi (x){\mathrm{\Sigma }}^{-1}(x-\mu )=0⟹x=\mu .$$ For any positive value $c\in (0,\pi (\mu ))$, ${\pi }^{-1}(\{c\})$ denotes the $c$-level set of $\pi$. We wish to show that this is a $(d-1)$-dimensional ellipsoid in ${\mathbb{R}}^d$. The level set is the following set of points $${\mathcal{E}}_c:=\{x\in {\mathbb{R}}^d:(2\pi {)}^{-d/2}\text{det}(\mathrm{\Sigma }{)}^{-1/2}\exp (-\frac{1}{2}(x-\mu {)}^{\mathrm{\top }}{\mathrm{\Sigma }}^{-1}(x-\mu ))=c\}$$ Since $\log$ is a monotone function, this is equivalent to $${\mathcal{E}}_c:=\{x\in {\mathbb{R}}^d:-\frac{d}{2}\log (2\pi )-\frac{1}{2}\log det\mathrm{\Sigma }-\frac{1}{2}(x-\mu {)}^{\mathrm{\top }}{\mathrm{\Sigma }}^{-1}(x-\mu )=\log c\}$$ Rearranging terms we get $${\mathcal{E}}_c:=\{x\in {\mathbb{R}}^d:(x-\mu {)}^{\mathrm{\top }}{\mathrm{\Sigma }}^{-1}(x-\mu )=b(c,\mathrm{\Sigma })\}$$ where we have defined the constant $$b(c,\mathrm{\Sigma })=-2\log c-d\log (2\pi )-\log det\mathrm{\Sigma }.$$ Since the LHS is a quadratic form, $b(c,\mathrm{\Sigma })$ must be positive so we can divide through to obtain $${\mathcal{E}}_c:=\{x\in {\mathbb{R}}^d:(x-\mu {)}^{\mathrm{\top }}(b\mathrm{\Sigma }{)}^{-1}(x-\mu )=1\}$$ Therefore they are ellipsoids centered at $v=\mu$ with matrix $\mathsf{A}=(b\mathrm{\Sigma }{)}^{-1}$.

Sampling from the level set

To sample exactly from this level set we follow the following algorithm:

1. Sample $x\sim \mathcal{N}(0_d,{\mathrm{I}}_d)$
2. Set $\hat{x}:=x/\Vert x\Vert$ so that $\hat{x}\sim \mathcal{U}({\mathbb{S}}^{d-1})$
3. Obtain sample $y:=L\hat{x}$ where $L$ is the Cholesky decomposition of $b\mathrm{\Sigma }$, i.e. $b\mathrm{\Sigma }=L{L}^{\mathrm{\top }}$