22
$\begingroup$

Given two multi-variate gaussians distrubtions, given by mean & covariance, G1(m1,sigma1) & G2(m2,sigma2), what are the formulae to find the product i.e G1 * G2 ? 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
    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.html2012-06-12

4 Answers 4

14

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.

13

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

  • 0
    The formulas in the document also premultiply the final PDF (by c_c) . Do your formula's take this into account?2018-05-15
8

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

  • 0
    This is a fantastic answer. Thanks especially for pointing out that the matrix inverse can be avoided entirely.2018-10-21
0

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, Matrix> MultiVariateGaussianProduct(List, Matrix>> 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, Matrix>(combinedMean, combinedCoVariance);   } 

N.B

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