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.
 A: Rejection sampling  might work: sample from one distribution and accept based on the other density


*

*Sample $x$ from $\omega(x)$.

*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.
A: 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.
