Product of Two Multivariate Gaussians Distributions 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.
 A: 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.
A: 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
A: 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 :)
A: 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;

    }

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


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

}

N.B 


*

*This method allows for an indetermate number of distributions.

*I used  MathNet's implementation of Matrices/Vectors.

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