# Product of Two Multivariate Gaussians Distributions

Given two multivariate gaussians distributions, given by mean and covariance, $$G_1(x; \mu_1,\Sigma_1)$$ and $$G_2(x; \mu_2,\Sigma_2)$$, what are the formulae to find the product i.e. $$p_{G_1}(x) p_{G_2}(x)$$ ?

And if one was looking to implement this in c++, what would an efficient way of doing it?

Go easy, I am primarily a computer scientist and not a pure mathematician.

Any help much appreciated.

• What do you mean by "find the product"? Do you want to do the distribution of the product or something else? Also, what "product" are you interested in? Is $G_1 \cdot G_2$ an inner (i.e., dot) product? An outer product? Something else? Recall that $G_1$ and $G_2$ are vectors, so, in particular, the inner product wouldn't make sense if $G_1$ and $G_2$ are of differing dimensions. Jun 11, 2012 at 23:41
• I suspect what the question was intended to mean is this: What is the distribution of the product of two random variables, whose distributions are those Gaussian distributions? Probably they were intended to be independent---that's an assumption people often forget to mention. Definitely the poster should clarify. Jun 11, 2012 at 23:54
• I mean the d-dimensional multivariate case of this tina-vision.net/docs/memos/2003-003.pdf Jun 11, 2012 at 23:55
• Essentially the maths being conducted in this matlab function (in the case where there are two d-dimensional gaussian distributions. ee.ic.ac.uk/hp/staff/dmb/voicebox/doc/voicebox/gausprod.html Jun 12, 2012 at 0:00

An alternative expression of the PDF proportional to the product is:

$\Sigma_3 = \Sigma_1(\Sigma_1 + \Sigma_2)^{-1}\Sigma_2$

$\mu_3 = \Sigma_2(\Sigma_1 + \Sigma_2)^{-1}\mu_1 + \Sigma_1(\Sigma_1 + \Sigma_2)^{-1}\mu_2$

The advantage of this form for computation is that it requires only one matrix inverse.

Denoting the product by $G_3 = (\mu_3, \Sigma_3)$, the formulas are:

$\Sigma_3 = (\Sigma_1^{-1}+\Sigma_2^{-1})^{-1}$

$\mu_3 = \Sigma_3\Sigma_1^{-1}\mu_1 + \Sigma_3\Sigma_2^{-1}\mu_2$

as found in the Matrix cookbook (Section 8.1.8):

http://compbio.fmph.uniba.sk/vyuka/ml/old/2008/handouts/matrix-cookbook.pdf

• The formulas in the document also premultiply the final PDF (by c_c) . Do your formula's take this into account?
– Ben
May 15, 2018 at 14:15
• Could you also provide a derivation? Nov 26, 2020 at 14:05

I depends on the information you have and the quantities you want to get out.

• If you have the covariance matrices themselves then you should use the formula $$\Sigma_3 = \Sigma_1(\Sigma_1 + \Sigma_2)^{-1}\Sigma_2$$ $$\mu_3 = \Sigma_2(\Sigma_1 + \Sigma_2)^{-1}\mu_1 + \Sigma_1(\Sigma_1 + \Sigma_2)^{-1}\mu_2$$ The computationally efficient and numerically stable way to do this would be to take the Cholesky decomposition of $\Sigma_1 + \Sigma_2$ (the Cholesky decomposition is probably a standard part of whatever matrix library you're using). $$LL^T = \Sigma_1 + \Sigma_2$$ Then compute \begin{align*} \tilde \Sigma_1 &= L^{-1}\Sigma_1 & \tilde \Sigma_2 &= L^{-1}\Sigma_2\\ \tilde \mu_1 &= L^{-1}\mu_1 & \tilde \mu_2 &= L^{-1}\mu_2 \end{align*} Which is efficient because $L$ is lower triangular (make sure to make use of built-in linear solve functions of your matrix library). The full solution is $$\Sigma_3 = \tilde \Sigma_1^T \tilde\Sigma_2\\ \mu_3 = \tilde \Sigma_2^T \tilde \mu_1 + \tilde \Sigma_1^T \tilde \mu_2$$

• If however you have the inverse covariances, because Gaussian distributions are expressed in terms of the inverse covariance, the computation can be even more efficient. In that case you should compute $$\Sigma_3^{-1} = \Sigma_1^{-1} + \Sigma_2^{-1}\\ \mu_3 = \Sigma_3(\Sigma_1^{-1}\mu_1 + \Sigma_2^{-1}\mu_2)$$ When you compute the expression for the mean use a built in linear solve function; it can be more efficient and numerically stable than actually computing the inverse of $\Sigma_3^{-1}$.

The C++ implementation is up to you :)

• This is a fantastic answer. Thanks especially for pointing out that the matrix inverse can be avoided entirely. Oct 21, 2018 at 23:16
• @Patrick, I am not sure your first explicit formula is right. $\Sigma_1(\Sigma_1+\Sigma_2)^{-1}\Sigma_2$ is not necessarily symmetric. Jan 21 at 17:23

Since the poster referred to c++ in their question, here's a code-based answer in a language with a similar syntax, viz. c#:

public Tuple<Vector<double>, Matrix<double>> MultiVariateGaussianProduct(List<Tuple<Vector<double>, Matrix<double>>> vm)
{

//v: Mean Vector
//m: CoVariance Matrix

// m-1
var mSumInv = vm[0].Item2.Inverse();
// v/m
var mInvV = mSumInv*vm[0].Item1;

for (int i = 1; i < vm.Count; i++)
{
// m-1 +
var mInv= vm[i].Item2.Inverse();

mSumInv += mInv;
// v/m +
mInvV += mInv * vm[i].Item1;

}

//(m-1)-1
var combinedCoVariance = mSumInv.Inverse();
// m*(v/m)
var combinedMean =combinedCoVariance* mInvV;

return new Tuple<Vector<double>, Matrix<double>>(combinedMean, combinedCoVariance);

}


N.B

• This method allows for an indetermate number of distributions.
• I used MathNet's implementation of Matrices/Vectors.

Following up on @benno's answer, this can be generalized to more than two Gaussians. The product of $$K$$ Gaussians, indexed by $$k$$, is proportional to a Gaussian with the following covariance $$\Sigma$$ and mean $$\mu$$:

$$\Sigma = \Big(\sum^K_{k=1}\Sigma_k^{-1}\Big)^{-1}$$

$$\mu = \Big(\sum^K_{k=1}\Sigma_k^{-1}\Big)^{-1} \Big(\sum^K_{k=1} \Sigma_k^{-1} \mu_k\Big)$$