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.

  • 106
  • 7
  • 491
  • 1
  • 4
  • 6
  • 1
    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. – cardinal Jun 11 '12 at 23:41
  • 2
    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. – Michael Hardy Jun 11 '12 at 23:54
  • I mean the d-dimensional multivariate case of this http://www.tina-vision.net/docs/memos/2003-003.pdf – oracle3001 Jun 11 '12 at 23:55
  • 1
    Essentially the maths being conducted in this matlab function (in the case where there are two d-dimensional gaussian distributions. http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/doc/voicebox/gausprod.html – oracle3001 Jun 12 '12 at 00:00

5 Answers5


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.

  • 321
  • 2
  • 3

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):


  • 113
  • 4
  • 1,466
  • 1
  • 14
  • 22

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 :)

  • 342
  • 2
  • 6
  • This is a fantastic answer. Thanks especially for pointing out that the matrix inverse can be avoided entirely. – David J. Harris Oct 21 '18 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. – Cupitor Jan 21 '22 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;


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

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



  • 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)$

Betterthan Kwora
  • 284
  • 4
  • 11