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
    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$ 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
    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
    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
    May 15, 2018 at 14:15
  • $\begingroup$ Could you also provide a derivation? $\endgroup$ 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 :)

  • $\begingroup$ This is a fantastic answer. Thanks especially for pointing out that the matrix inverse can be avoided entirely. $\endgroup$ 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
    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;


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.