4

How can one efficiently sample $x \in \mathbb{R}^N$ from a multivariate normal distribution $x \sim \exp(- \frac{1}{2}x^T \Sigma^{-1} x)$ given a normalization constraint $x^T x = 1$? In my application, $N$ would typically be order 10 or 20. If the covariance matrix $\Sigma^{-1}$ were the identity, then the target distribution would be uniform in the $(N-1)$-dimensional sphere, and there are straightforward algorithms. I get stuck when $\Sigma$ is not the identity.

One can reframe the problem by working in the eigenbasis of $\Sigma^{-1}$. Effectively, we can assume that $\Sigma^{-1}$ is a diagonal matrix. Then $x^T \Sigma^{-1} x$ has the form $\sum_i \lambda_i x_i^2$, where $\lambda_{i=1,\dots N}$ denote the (real) eigenvalues of $\Sigma^{-1}$. The task can be alternatively written as a problem of sampling nonnegative real numbers $\{q\}_{i=1,\dots,N} \sim \exp(-\sum_i \lambda_i q_i / 2)$ subject to the constraint ("given that") $\sum_i q_i = 1$, but note the change of variables $x \rightarrow q$ introduces a Jacobian factor to the integration measure on $d q$ that makes it difficult to "integrate out" (marginalize over?) some of the $q_i$ variables.

This question has a similar flavor, but is not quite the same. Here, my constraint on $x$ is quadratic.

@whuber Made a very interesting connection with "sum of scaled $\Gamma(1/2)$ variables", with approximation methods discussed here: Generic sum of Gamma random variables

Kipton Barros
  • 141
  • What can you say about the distribution of the eigenvalues? There might be one good algorithm if one eigenvalue is much larger than the others, and another good algorithm if the eigenvalues are all similarly sized. – Matt F. Jun 30 '22 at 02:15
  • 1
    To a first approximation, let's say there are 10 eigenvalues, uniformly spaced between 0 and 4 – Kipton Barros Jun 30 '22 at 05:59
  • So a nice test case (for which I have no ideas yet) is $$\Sigma= \begin{pmatrix}1&1&1\0&3&1\0&0&4\end{pmatrix}$$ with eigenvalues $1,3,4$. – Matt F. Jun 30 '22 at 06:50
  • There are two possible distributions and the wording of your question doesn't make it entirely clear which one you are after. Do you want to sample from the conditional distribution of $N(0,\Sigma)$ conditional on $x^Tx=1$? Or do you want to sample values from $x/\sqrt{x^Tx}$ where $x$ is $N(0,\Sigma)$? Both satisfy the constraint but the two distributions are not the same. – Gordon Smyth Jun 30 '22 at 07:53
  • @GordonSmyth, this is about the first of your questions. – Matt F. Jun 30 '22 at 08:11
  • Correct, thank you for clarifying. Also, I would suggest to work in the eigenbasis, this way you can effectively \Sigma as a diagonal matrix, e.g. diagm([1,2,3,4]), then the components of x squared become the q_i variables. – Kipton Barros Jun 30 '22 at 12:34
  • I have used the eigenbasis of the covariance matrix to decompose a multinormal distribution many times in my own research, but I do not think that approach will lead to an accurate simulation algorithm for your problem. The theoretical approximations to sums of gammas with different scales are all quite rough. – Gordon Smyth Jun 30 '22 at 22:08
  • Do you have any choice about the normalization contraint applied? From a statistical point of view, we would usually choose to normalize such that $x^T\Sigma^{-1}x=1$ because that constraint has much better statistical properties. – Gordon Smyth Jul 01 '22 at 07:58
  • Thanks for the correction. It's fine to redefine $\Sigma$ that way. This problem is coming from quantum physics, and I tried to translate it into statistic language. – Kipton Barros Jul 01 '22 at 16:34
  • If you can redefine the constraint as $x^T\Sigma^{-1}x$ then the simulation problem becomes vastly easier because it transforms to the uniform-on-sphere problem. – Gordon Smyth Jul 01 '22 at 20:41
  • It's fine to use the constraint $x^\intercal \Sigma^{-1}x =1$ that @GordonSmyth wrote? But that means all x fulfilling the constraint have the same density so the problem gets really easy. – all feedback welcome Jul 01 '22 at 20:58
  • @GordonSmyth exactly, in that case you can just sample $x$ and compute $x' = \frac x{\left\lVert \Sigma^{-\frac12}x\right\rVert}$. That works since if you were to reparameterize by $x=\Sigma^\frac12\epsilon$, with white noise $\epsilon$, you'd get the desired uniform density over $\epsilon' = \frac\epsilon{\lVert \epsilon\rVert}$ due to rotation-symmetry. – all feedback welcome Jul 01 '22 at 21:23
  • Hi, Gordon and "all feedback" -- you're right! Unfortunately I mispoke: in my application, the constraint is x^T x = 1, and cannot be rescaled. For my application, I ended up doing some form of MCMC sampling of x, where some components of x are exponentially more represented than others, and this works pretty well. Thanks again! – Kipton Barros Jul 04 '22 at 16:36

0 Answers0