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 '12 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$ – Michael Hardy Jun 11 '12 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 '12 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 '12 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.

| cite | improve this answer | |

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


| cite | improve this answer | |
  • 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 '18 at 14:15

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

| cite | improve this answer | |
  • $\begingroup$ This is a fantastic answer. Thanks especially for pointing out that the matrix inverse can be avoided entirely. $\endgroup$ – David J. Harris Oct 21 '18 at 23:16

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.
| cite | improve this answer | |

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.