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.

  • 1
    $\begingroup$ 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. $\endgroup$
    – cardinal
    Commented Jun 11, 2012 at 23:41
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Jun 11, 2012 at 23:54
  • $\begingroup$ I mean the d-dimensional multivariate case of this tina-vision.net/docs/memos/2003-003.pdf $\endgroup$
    – oracle3001
    Commented Jun 11, 2012 at 23:55
  • 1
    $\begingroup$ 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 $\endgroup$
    – oracle3001
    Commented Jun 12, 2012 at 0:00

5 Answers 5


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


  • 1
    $\begingroup$ The formulas in the document also premultiply the final PDF (by c_c) . Do your formula's take this into account? $\endgroup$
    – Ben
    Commented May 15, 2018 at 14:15
  • $\begingroup$ Could you also provide a derivation? $\endgroup$ Commented 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 :)

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

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

  • 1
    $\begingroup$ This dovetails beautifully with the last bit of @Patrick's answer. If you store the inverses instead of the covariance matrices directly, the product of K PDFs can be computed with a single matrix inversion. $\endgroup$
    – hosford42
    Commented Jun 28, 2023 at 21:36

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.

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .