# Sampling from a product of two Multivariate Gaussians

I have a multivariate Gaussian defined as follows: $$p(x) = \omega(x)\gamma(x)$$ where $\omega$ and $\gamma$ both are multivariate Gaussians and from which I can sample very efficiently given due to special structure in their covariances: one is diagonal, of the other I know the diagonalization, which is very sparse.

I wonder if there is a way to sample from $p$, given samples from $\omega$ and $\gamma$. The standard way of doing a Cholesky on the joint covariance is not efficient enough.

• Presumably there's an implicit normalization constant in the product and $p$, $\omega$, and $\gamma$ refer to PDFs. In that case $p$ also is a multivariate Gaussian. Whether its variance has any special structure depends on the means and eigensystems of $\omega$ and $\gamma$. As far as I can tell, then, absent any specific information about those relationships, the only possible (valid) answers to this question will simply explain how to sample from an arbitrary multivariate Gaussian. – whuber Sep 18 '14 at 17:14

Rejection sampling might work: sample from one distribution and accept based on the other density

1. Sample $x$ from $\omega(x)$.
2. Return $x$ with probability $\gamma(x) / \gamma(\mu_\gamma)$, otherwise goto 1.

The roles of $\gamma$ and $\omega$ may be switched due to symmetry, so one could try both to see which is more efficient.

However, if neither component distribution is 'close' to the joint distribution, this will not be efficient as most samples get rejected.

• Since $p$ is also multivariate Gaussian, it seems unlikely that this rejection sampling method would be superior to any generic method of sampling from multivariate Gaussians. Under what conditions on $\gamma$ and $\omega$ do you suppose rejection sampling would have better performance? – whuber Sep 18 '14 at 17:16

First, note that the product $$p(x) = \omega(x)\gamma(x),$$ is not a density function, unless you normalise it.

What you can do is to take $\omega(x)$ as an instrumental "likelihood" and $\gamma(x)$ as an instrumental "prior" (not in a strict sense, it is only a trick to construct a tractable sampler). This leads you to a conjugate model (since the normal prior is the conjugate prior of the mean of a normal sampling model):

http://en.wikipedia.org/wiki/Conjugate_prior#Continuous_distributions

This product should be normal, with a certain mean and covariance structure, which is tractable due to conjugacy.

• The (normalized) product is indeed normal, but the last sentence of the question mentions that computing the mean and covariance of the 'posterior' and sampling from that in the standard way is not efficient enough. – Juho Kokkala Sep 18 '14 at 10:09
• @JuhoKokkala I wonder why is it that difficult since the variance of the conjugate posterior is kind of simple. – Coco Sep 18 '14 at 10:20
• That's what I meant with the last paragraph. – bayerj Sep 18 '14 at 11:07
• @bayerj What's the dimension of $x$? – Coco Sep 18 '14 at 11:15
• The dimension can be anything, but practical are a few hundred. However, There will be thousands of these problems in an inner loop of an optimisation, which should run as fast as possible. – bayerj Sep 18 '14 at 11:36