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

  • $\begingroup$ 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. $\endgroup$
    – Matt F.
    Jun 30, 2022 at 2:15
  • 1
    $\begingroup$ To a first approximation, let's say there are 10 eigenvalues, uniformly spaced between 0 and 4 $\endgroup$ Jun 30, 2022 at 5:59
  • $\begingroup$ 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$. $\endgroup$
    – Matt F.
    Jun 30, 2022 at 6:50
  • $\begingroup$ 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. $\endgroup$ Jun 30, 2022 at 7:53
  • $\begingroup$ @GordonSmyth, this is about the first of your questions. $\endgroup$
    – Matt F.
    Jun 30, 2022 at 8:11


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.