The Venerable Gaussian

The Gaussian distribution underpins a broad range of methods in data analysis. This is often due to its remarkable analytical properties, which make working with it relatively easy.

To use it effectively however, one often has to trudge through a fair bit of linear algebra and memorise some identities. This article aims to collect several of these in one place; intended as both an initial exposition and a reference.

We first introduce the univariate Gaussian in Section 1 and highlight some of its properties before generalising to the multivariate case in Section 1.2. We then derive the matrix inversion lemma and matrix determinant lemma in Section 2, which sets us up for a treatment of Gaussian joint distributions in Section 3 in which we derive the equations for the conditional of a joint Gaussian distribution. Finally, we conclude by looking at linear Gaussian systems in the multivariate case in Section 4.

Introducing the Gaussian Distribution

Rather than giving the formula for the Gaussian distribution straight away, we will work up to it, pointing out some of its useful properties as we go along. For clarity, we’ll begin with the univariate case before generalising to the full Gaussian distribution.

At its core, a Gaussian function is the exponential of a negative quadratic:

\[\begin{equation}\label{eq:gaussian_function} f(x) = e^{-\frac{x^{2}}{2}}, \end{equation} \]

where we have divided by \(2\) for mathematical convenience later. Plotting this function gives us the familiar bell curve:

This function is symmetric around the peak, which lies at \(x=0\) in the expression above, and quickly tapers off 0: In fact, the Gaussian is notable for how quickly it tapers off when far from the mode. This distribution is said ot have light tails, meaning that there is not much support in the tails of the distribution far from the mean. ^0[0] for large \(\vert x \vert \gg 1\).

The expression above is not a valid density function however, since it is not normalised (it does not integrate to 1). In order to find the proper normalisation term, we must compute the Gaussian integral:

\[I_{1} \doteq \int_{-\infty}^{\infty} e^{-\frac{x^{2}}{2}}dx, \]

where the \(1\) in the subscript of \(I_{1}\) indicates that the above distribution has unit variance—we will later denote more general Gaussian integrals for distributions with variance \(\sigma^{2}\) as \(I_{\sigma^{2}}\).

This integral is not trivial to compute, but there is a neat trick 1: There are actually several approaches, with the approach using polar coordinates, due to Poisson, arguably being the most straightforward. For another, more geometric take on this integral, see this blog. ^1[1] to do it: by squaring the integral, we can change the variables to polar coordinates \((x, y) = (r \cos \theta, r \sin \theta)\), which will result in two slightly simpler elementary integrals.

We therefore square the integral, and substitute two dummy integration variables \(\alpha\) and \(\beta\) such that \(x = \alpha\beta\), giving:

\[I^{2}_{1} = \left(\int_{-\infty}^{\infty}e^{-\frac{x^{2}}{2}}dx\right)^{2} = \int_{-\infty}^{\infty}e^{-\frac{\alpha^{2}}{2}}d\alpha \int_{-\infty}^{\infty}e^{-\frac{\beta^{2}}{2}}d\beta = \iint_{-\infty}^{\infty}e^{-\frac{1}{2}(\alpha^{2} + \beta^{2})}d\alpha d\beta. \]

We can now interpret the above using polar coordinates, and so substitute the integral measure \(d\alpha d\beta = rdrd\theta\) while also recalling that \(r = \sqrt{\alpha^{2} + \beta^{2}}\), which leaves us with two elementary integrals to compute:

\[\begin{align*} I^{2}_{1} = \iint_{-\infty}^{\infty}e^{-\frac{1}{2}(\alpha^{2} + \beta^{2})} d\alpha d\beta &= \int^{\infty}_{0}r\ dr \int_{0}^{2\pi} e^{-\frac{r^{2}}{2}}d\theta \\ &= 2\pi \int^{\infty}_{0} r e^{-\frac{r^{2}}{2}}dr = 2\pi \left[-e^{-\frac{r^{2}}{2}}\right]_{r=0}^{r=\infty} \\ &= 2\pi. \end{align*} \]

Remembering to take the square root again, the solution to the Gaussian integral comes to:

\[\begin{equation} \label{eq:std_gaussian_integral} I_{1} = \int^{\infty}_{-\infty}e^{-\frac{x^{2}}{2}}dx = \sqrt{2\pi}. \end{equation} \]

This is the normalisation term that we need, in order to turn the Gaussian function of Equation \(\ref{eq:gaussian_function}\) into a Gaussian probability density function (PDF) with zero mean and unit variance. Hence, dividing Equation \(\ref{eq:gaussian_function}\) by \(\sqrt{2\pi}\), we get the (univariate) standard normal distribution:

\[\begin{equation} \label{eq:univariate_standard_gaussian} \gN(\ervx; 0, 1) \doteq \frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}. \end{equation} \]

To generalise this to Gaussian distributions with different positive variances \(\sigma^{2} \gt 0\), we can divide the \(x\) term in the Gaussian function by \(\sigma^{2}\), giving \(f_{\sigma^{2}}(x) = \exp(-\frac{x^{2}}{2\sigma^{2}})\). However, we now need to re-calculate the normalisation term. Fortunately, since the variance does not depend on the integrand \(x\), we can just factor it out of the integral, make the \(u\)-substitution \(u = \frac{x}{\sqrt{\sigma^{2}}}\), and make use of our result from Equation \(\ref{eq:std_gaussian_integral}\):

\[I_{\sigma^{2}} \doteq \int_{-\infty}^{\infty}e^{-\frac{x^{2}}{2\sigma^{2}}}dx = \sqrt{\sigma^{2}}\int_{-\infty}^{\infty}e^{-\frac{u^{2}}{2}}du = \sqrt{2\pi \sigma^{2}}. \]

Dividing the Gaussian function by this new normalisation term, the more general Gaussian density function with variance \(\sigma^{2}\) is therefore 2: It is conventional to denote the variance as \(\sigma^{2}\), which lets us easily recover the standard deviation or ‘scale’ of \(x\) as \(\sigma\). ^2[2]

\[\gN(\ervx; 0, \sigma^{2}) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^{2}}{2\sigma^{2}}}. \]

Plotting this new Gaussian density function, the bell curve is still symmetric around \(x=0\), however \(\sigma^{2}\) now characterises its broadness, and it tapers off for large \(x\); \(\vert x \vert \gg \sqrt{\sigma^{2}}\).

Generalising this further, we can translate the mean of the distribution 3: Intuitively: shift the centre of the bell ^3[3] by \(\mu \in \R\), which gives the full definition of a univariate Gaussian density:

Definition (Univariate Gaussian). The univariate Gaussian distribution is defined as

\[\begin{equation} \label{eq:univariate_gaussian} \gN(\ervx; \mu, \sigma^{2}) \doteq \frac{1}{\sqrt{2\pi}\sigma}\exp \left(-\frac{(x - \mu)^2}{2\sigma^{2}}\right), \end{equation} \]

where \(\mu\in\R\) is the mean of the distribution of \(\ervx\) (and coincidentally also the mode of the distribution), \(\sigma^{2} \in \R\) is the variance of \(\ervx\), and \(\sigma\) is the standard deviation of \(\ervx\).
\(\blacksquare\)

The fact that the centre value \(\mu\) is the mean of the distribution can be seen by taking the expectation of the above:

\[\begin{align*} \E[\ervx] \doteq \int_{-\infty}^{\infty} \gN(x; \mu, \sigma^{2}) x\ dx &= \frac{1}{\sqrt{2\pi}\sigma} \int_{-\infty}^{\infty}e^{-\frac{(x - \mu)^{2}}{2\sigma^{2}}}x\ dx \\ &= \frac{1}{I_{\sigma^{2}}}\int_{-\infty}^{\infty}e^{-\frac{\omega^{2}}{2\sigma^{2}}}(\mu + \omega)d\omega \\ &= \frac{\mu I_{\sigma^{2}}}{I_{\sigma^{2}}} + \frac{1}{I_{\sigma^{2}}}\int_{-\infty}^{\infty}\left(e^{-\frac{\omega^{2}}{2\sigma^{2}}}\omega\right)d\omega \\ &= \mu, \end{align*} \]

where on the second line above we made the substitution \(\omega \doteq x - \mu\). On the last line we noticed that the integrand of the second term is odd with respect to the sign flip of the integration variable \(\omega \leftrightarrow -\omega\), and hence it integrates to 0. In other words, the negative bits, when \(\omega < 0\), cancel out the positive bits \(\omega > 0\):

As some final basic properties of the univariate Gaussian, we note that the variance 4: of any distribution, not just the Gaussian ^4[4] is the expected squared distance of samples drawn from the distribution to its mean, \(\mu\). That is, for samples \(\{\ervx_{i}\}_{i=1}^{N} \sim X = \gN(\mu, \sigma^2)\), we have

\[\begin{align*} \Var(X) &= \E\left[(X - \mu)^2\right] = \E[X^2] - 2\mu\E[X] + \mu^2 \\ &= \E[X^2] - \mu^2 \\ &\approx \frac{1}{N}\sum_{i=1}^{N} \ervx_{i}^2 - \mu^{2}, \end{align*} \]

where the first line follows from the linearity of expectation, and the second from the definition that \(\E[X] = \mu\).

As alluded to above when mentioning how quickly the Gaussian function tapers off for values far from the mean, the Gaussian distribution is particularly unsuitable for modelling populations drawn from heavy tailed distributions, in which some events occur very far from the mean of the distribution. To see this, consider the following diagram, which shows the probability mass placed in each consecutive, 1 standard deviation ‘bin’ away from the mean:

It can finally be worth pointing out that we have symmetry in the random variable \(\ervx\) and the mean \(\mu\):

\[\gN(\ervx; \mu, \sigma^{2}) = \gN(\mu; \ervx, \sigma^{2}). \]

Univariate Linear Gaussian Systems

Here we will take a brief look at a simple case of a univariate linear Gaussian system to gain some intuition about working with Gaussians. For the multivariate case, see Section 4.

We might use a linear Gaussian system when we have a parameter of interest, \(\ervx\) which is possibly hidden or latent. We are interested in modelling the value of this parameter, but we can only make noisy observations of another value \(\ervy\), which we hypothesise to be linearly related to \(\ervx\). We can formalise this as follows:

Definition (Univariate linear Gaussian system). A linear Gaussian system arises when a Gaussian conditional distribution is centred at (some linear function of) another Gaussian distributed random variable. For clarity of exposition, in this section we will use an identity function; see Section 4 for the more general case. Thus, we have

\[\begin{align*} \label{eq:univariate_linear_gaussian_sys} p(\ervx) &= \gN(\ervx; \mu, \sigma^2) \\ p(\ervy \vert \ervx) &= \gN(\ervy; \ervx, \nu^2). \end{align*} \]

Having observed \(\ervy\), we are usually interested in inferring 5: That is, applying Bayes rule to find the posterior. ^5[5] \(p(\ervx \vert \ervy)\). The fact that the Gaussian remains closed under multiplication 6: products of exponentiated quadratics remain exponentiated quadratics. We will derive a more genearl expression for the product of two Gaussian distributions later ^6[6] allows us to deduce the following:

\[\begin{align} \label{eq:gaussian_posterior} p(\ervx \vert \ervy) &= \frac{p(\ervx)p(\ervy \vert \ervx)}{\int p(\ervx) p(\ervy \vert \ervx) d\ervx} = \frac{p(\ervx, \ervy)}{p(\ervy)} \\ &= \gN(\ervx; m, s^2), \end{align} \]

where

\[\begin{align*} s^2 &\doteq (\sigma^{-2} + \nu^{-2})^{-1}, \\ m &\doteq \big(\sigma^{-2}\mu + \nu^{-2}\ervy\big)s^{2}. \end{align*} \]

\(\blacksquare\)

Proof. We can omit the normalising constants in the Gaussians, since these can be factored out of both the integral in the denominator of Eq. \(\ref{eq:gaussian_posterior}\) as well as the product in the numerator, at which point they will cancel out.

Thus, we take

\[p(\ervx) \propto \exp\left(-\frac{1}{2}\frac{(\ervx-\mu)^{2}}{\sigma^{2}}\right), \]

and

\[p(\ervy \vert \ervx) \propto \exp\left(-\frac{1}{2}\frac{(\ervy - \ervx)^{2}}{\nu^{2}}\right). \]

To proceed, note that \(p(\ervx \vert \ervy) \propto p(\ervx, \ervy) = p(\ervy \vert \ervx)p(\ervx)\), where the final product is the exponential of a negative quadratic. By writing this quadratic as a perfect square, we can more easliy match terms and recover the mean and variance of the Gaussian \(p(\ervx \vert \ervy)\); \(m\) and \(s^{2}\).

Multiplying these two exponents gives \(\exp\left(-\frac{1}{2}\left(\frac{(\ervx-\mu)^{2}}{\sigma^{2}} + \frac{(\ervx - \ervy)^{2}}{\nu^{2}}\right)\right)\). Since we want to recover the first two moments required to specify the resulting Gaussian, \(\exp\left(-\frac{1}{2}\frac{(\ervx - m)^{2}}{s^{2}}\right)\), writing the perfect square in the form \(a\left(x + \frac{b}{2a}\right)^{2}\) makes it trivial to find the values of \(m\) and \(s^{2}\). Proceeding,

\[\begin{align*} &\frac{(\ervx - \mu)^{2}}{\sigma^{2}} + \frac{(\ervx - \ervy)^{2}}{\nu^{2}} \\ &= \frac{\ervx^{2} - 2\ervx\mu + \mu^{2}}{\sigma^{2}} + \frac{\ervx^{2} - 2\ervx\ervy + \ervy^{2}}{\nu^{2}} \\ &= \left(\frac{1}{\sigma^{2}} + \frac{1}{\nu^{2}}\right)\ervx^{2} - 2\left(\frac{\mu}{\sigma^{2}} + \frac{\ervy}{\nu^{2}}\right) \ervx + \text{const} \\ &= \underbrace{\left(\frac{1}{\sigma^{2}} + \frac{1}{\nu^{2}}\right)}_{s^{-2}}\left(\ervx - \underbrace{\left(\frac{\mu}{\sigma^{2}} + \frac{\ervy}{\nu^{2}}\right) \left(\frac{1}{\sigma^{2}} + \frac{1}{\nu^{2}}\right)^{-1}}_{m}\right)^{2} + \text{lin.} + \text{const} \end{align*} \]

giving the variance and mean of \(p(\ervx \vert \ervy)\) as

\[\begin{align*} s^{2} &= (\sigma^{-2} + \nu^{-2})^{-1} \\ m &= (\mu \sigma^{-2} + \ervy \nu^{-2})s^{2}. \end{align*} \]

\(\square\)

Introducing the Multivariate Gaussian

Most interesting applications of the Gaussian require us to jointly model two or more random variables. Hence, we will now consider the multivariate Gaussian (which we will simply refer to as the Gaussian).

Definition (Gaussian). The Gaussian distribution is a joint probability density over multiple variables. If the number of variables of interest is \(D\), then we work with the \(D\)-dimensional random vector \(\rvx\), and the Gaussian PDF is defined as

\[\gN(\rvx \vert \vmu, \mSigma) = (2\pi)^{-D/2}\vert \mSigma \vert^{-1/2} \exp\left(-\frac{1}{2}(\rvx - \vmu)^{\top}\mSigma^{-1}(\rvx - \vmu)\right), \]

where we now have a mean vector \(\vmu\), as well as a covariance matrix \(\mSigma\). The notation \(\vert \mSigma \vert\) denotes the determinant of the matrix \(\mSigma\).

\(\blacksquare\)

If \(D = 1\), and \(\mSigma = \mI\) (the identity matrix) then we can see that we recover the univariate Gaussian exactly.

Recall that in the univariate case, we had a quadratic expression in the exponent; in the multivariate case we now have a quadratic form:

\[\Delta^{2} = (\rvx - \vmu)^{\top}\mSigma^{-1}(\rvx - \vmu). \]

For \(D > 1\) and \(\mSigma\) still the identity, the quantity \(\Delta\) is just the Euclidean distance (or \(\ell_2\) norm) between the random variable and the mean.

Interpreting the Covariance Matrix

For more general covariance matrices, the quadratic form corresponds to the Mahalanobis distance between the data vector \(\rvx\) and the mean vector \(\vmu\). Intuitively, the Mahalanobis distance corresponds to the Euclidean distance in a transformed coordinate system, with a shift of \(\vmu\), and a rotation by some rotation matrix \(\rmU\).

The rotation matrix \(\rmU\) comes from the eigendecomposition of the covariance matrix 7: You might notice that this form looks slightly different from the usual eigendecomposition, where we have a matrix inverse rather than a transpose on the last term: this is due to the properties of the covariance matrix, which is a real, symmetrical matrix. ^7[7] : \(\mSigma = \rmU\mLambda\rmU^{\top}\). Here \(\rmU\) is the orthonormal matrix of eigenvectors (i.e. with \(\rmU^{\top}\rmU = \rmI\)) and \(\mLambda\) is the diagonal matrix of eigenvalues.

The covariance matrix’s inverse can then be written 8: Using the fact that \((\rmA\rmB)^{-1} = \rmB^{-1}\rmA^{-1}\) ^8[8] as

\[\mSigma^{-1} = \rmU^{-\top}\mLambda^{-1}\rmU^{-1} = \rmU\mLambda^{-1}\rmU^{\top} = \sum_{i=1}^{D}\frac{1}{\lambda_{i}}\rvu_{i}\rvu_{i}^{\top}, \]

where \(\rvu_{i}\rvu_{i}^{\top}\) is the outer product of the \(i\)th eigenvector. Substituting this into the quadratic form in the exponent of the Gaussian, we get

\[\begin{align*} (\rvx - \vmu)^{\top}\mSigma^{-1}(\rvx - \vmu) &= (\rvx - \vmu)^{\top} \left(\sum_{i=1}^{D}\frac{1}{\lambda_{i}}\rvu_{i}\rvu_{i}^{\top}\right)(\rvx - \vmu) \\ &= \sum_{i=1}^{D}\frac{1}{\lambda_{i}}(\rvx - \vmu)^{\top}\rvu_{i}\rvu_{i}^{\top}(\rvx - \vmu) = \sum_{i=1}^{D}\frac{y_{i}^{2}}{\lambda_{i}}, \end{align*} \]

where we have defined the scalar \(y_{i} \doteq \rvu_{i}^{\top}(\rvx - \vmu)\). We may interpret \(\{y_i\}\) as a new coordinate system, defined by the orthonormal basis vectors \(\rvu_{i}\) that are shifted and rotated wrt. the original \(\ervx_{i}\) coordinates.

Recalling that the equation of an ellipse (corresponding to \(D=2\)) is

\[\frac{y_{1}^{2}}{\lambda_{1}} + \frac{y^{2}_{2}}{\lambda_{2}} = 1, \]

we can see that contours of equal probability density for a Gaussian lie along ellipses. This is often used as a way to visualise high-dimensional Gaussians, where we pick pairs of eigenvectors, usually with large eigenvalues, and plot the resulting ellipse.

Note that in order for the Gaussian to be well defined, we need all the eigenvalues of the covariance matrix to be strictly positive—if they weren’t, we would not be able to normalise the distribution. A suitable covariance matrix with positive eigenvalues is said to be positive definite, and it is common when working with Gaussians to check for the size of these eigenvalues: lots of values close to 0 may predict numerical instabilities.

Natural Parametrisation

We have so far been using the first two moments of a Gaussian distribution, the mean \(\vmu\) and covariance \(\mSigma\), as its parameters. However, the Gaussian is part of the exponential family of distributions, meaning that we can define a Gaussian in an altenative form—information form— using the so-called natural or canonical parameters.

Many procedures and manipulations involving Gaussians simplify vastly when using the information form of the Gaussian, as opposed to the moment form. For instance certain Gaussian identities, such as the product or KL divergence between two Gaussians, are much simpler to derive in information form. The information form is also helpful in other situations such as when updating beliefs over variables in Bayesian networks, or solving linear Gaussian systems.

We can obtain the natural parameters of the Gaussian by rearranging the moment form. Denoting the normalising terms as \(p \doteq (2\pi)^{-D/2}\vert \mSigma \vert^{-1/2}\), and expanding the brackets in the quadratic form, we get:

\[\begin{align*} \gN(\rvx; \vmu, \mSigma) &= p \cdot \exp\left(-\frac{1}{2}(\rvx - \vmu)^{\top}\mSigma^{-1}(\rvx - \vmu)\right) \\ &= \exp\left( \log p - \frac{1}{2}\vmu^{\top}\mSigma^{-1}\vmu + \rvx^{\top}\mSigma^{-1}\vmu - \frac{1}{2}\rvx^{\top}\mSigma^{-1}\rvx \right) \\ &\doteq \exp\left(g + \rvx^{\top}\rvh - \frac{1}{2}\rvx^{\top}\rmK\rvx \right) \end{align*} \]

where we have defined

\[\begin{equation} \label{eq:natural_parameters} \rmK \doteq \mSigma^{-1},\hspace{2em}\rvh \doteq \mSigma^{-1}\vmu,\hspace{2em}g \doteq \log p - \frac{1}{2} \vmu^{\top}\rmK\vmu. \end{equation} \]

\(\rmK\) is referred to as the precision matrix, and is also commonly denoted \(\mLambda\). \(\rvh\) is referred to as the precision-adjusted mean.

Generally for the Gaussian, there will always be a one-to-one mapping between the moments and the natural parameters. We can go back to the moment form by inverting the definitions above:

\[\begin{equation}\label{eq:natural_to_moment} \mSigma = \rmK^{-1}\hspace{2em}\vmu = \rmK^{-1}\rvh\hspace{2em}p=\exp(g+\frac{1}{2}\vmu^{\top}\mSigma{-1}\vmu) \end{equation} \]

KL Divergence Between Two Gaussians

The KL divergence between two multivariate Gaussian distributions is given by

\[\begin{align*} \KL[&\gN(\rvx \vert \vmu_{1}, \mSigma_{1})\Vert \gN(\rvx \vert \vmu_{2}, \mSigma_{2})] =\\ &= \frac{1}{2}\left[\Tr(\mSigma_{2}^{-1}\mSigma_{1}) + (\vmu_{2} - \vmu_{1})^{\top}\mSigma^{-1}_{2}(\vmu_{2} - \vmu_{1}) - D + \log \left(\frac{\det(\mSigma_{2})}{\det(\mSigma_{1})}\right)\right]. \end{align*} \]

In the scalar case, this reduces to

\[\KL[\gN(\ervx \vert \mu_{1}, \sigma_{1}) \Vert \gN(\ervx \vert \mu_{2}, \sigma_{2})] = \log \frac{\sigma_{2}}{\sigma_{1}} + \frac{\sigma^{2}_{1} + (\mu_{1} - \mu_{2})^{2}}{2\sigma_{2}^{2}} - \frac{1}{2}. \]

Inverting a Partitioned Matrix

Evaluating a Gaussian density involves inverting a covariance matrix 9: As well as computing its determinant. For most real-world applications, we compute the Cholesky decomposition of \(\mSigma\) which lets us re-use our computational work in both the determinant and matrix inversion calculation. The partitioned matrix inverse however has many useful corollaries, and is thus useful for understanding Gaussian manipulations. ^9[9] :

\[\gN(\rvx; \vmu, \mSigma) = (2\pi)^{-D/2}\vert \mSigma \vert^{-1/2} \exp \left(-\frac{1}{2}(\rvx - \vmu)^{\top} \mSigma^{-1} (\rvx - \vmu) \right). \]

Often, our problem gives this covariance matrix some structure; in particular, we can often partition it as follows:

\[\mSigma = \begin{bmatrix}\mSigma_{11} & \mSigma_{12} \\ \mSigma_{21} & \mSigma_{22}\end{bmatrix}. \]

We will assume that the sub-matrices \(\mSigma_{11}\) and \(\mSigma_{22}\) are themselves invertible—in many contexts this is not an inconvenience since these sub-matrices are themselves covariance matrices (e.g. for a prior or likelihood).

To find \(\mSigma^{-1}\), we can begin by block-diagonalising it 10: that is; making the off-diagonal terms \(\mSigma_{12}\) and \(\mSigma_{21}\) zero. The resulting block diagonal matrix has two, square main-diagonal blocks (of the same size as \(\mSigma_{11}\) and \(\mSigma_{22}\), respectively), and all off-diagonal blocks are zero matrices ^10[10] . To zero-out the top right block, we left-multiply \(\mSigma\) by a matrix \(\mathbf{A}\) such that \(\mSigma_{12} = \mathbf{A}_{11}\mSigma_{12} + \mathbf{A}_{12}\mSigma_{22} = \mathbf{0}\), leaving other entries unchanged. Setting \(\rmA_{11} = \rmI\) and \(\rmA_{12} = -\mSigma_{12}\mSigma_{22}^{-1}\) gives us what we need:

\[\underbrace{\begin{bmatrix} \mathbf{I} & -\mathbf{\Sigma_{12}}\mSigma_{22}^{-1} \\ \mathbf{0} & \mathbf{I}\end{bmatrix}}_{\mathbf{A}} \begin{bmatrix}\mSigma_{11} & \mSigma_{12} \\ \mSigma_{21} & \mSigma_{22}\end{bmatrix} = \begin{bmatrix} \mSigma_{11} - \mSigma_{12}\mSigma_{22}^{-1}\mSigma_{21} & \mathbf{0} \\ \mSigma_{21} & \mSigma_{22}\end{bmatrix}. \]

Similarly, to zero-out the bottom-left block, we right-multiply \(\mSigma\) by a matrix \(\mathbf{B}\) such that \(\mSigma_{21} = \mSigma_{21}\rmB_{11} + \mSigma_{22}\rmB_{21} = \mathbf{0}\) to obtain our block-diagonal matrix, which we will call \(\mathbf{\Xi}\)

\[\begin{bmatrix} \mSigma_{11} - \mSigma_{12}\mSigma_{22}^{-1}\mSigma_{21} & \mathbf{0} \\ \mSigma_{21} & \mSigma_{22}\end{bmatrix} \underbrace{ \begin{bmatrix} \mathbf{I} & \mathbf{0} \\ -\mSigma_{22}^{-1}\mSigma_{21} & \mathbf{I} \end{bmatrix}}_{\mathbf{B}} = \underbrace{\begin{bmatrix} \mSigma_{11} - \mSigma_{12}\mSigma_{22}^{-1}\mSigma_{21} & \mathbf{0} \\ \mathbf{0} & \mSigma_{22}\end{bmatrix}}_{\mathbf{\Xi}}. \]

In summary, we have \(\mathbf{A}\mSigma\mathbf{B} = \mathbf{\Xi}\). Each of these matrices are individually invertible:

\[\mathbf{B}^{-1}\mSigma^{-1}\mathbf{A}^{-1} = \mathbf{\Xi}^{-1}, \]

and rearranging gives us an expression for \(\mSigma^{-1}\):

\[\mSigma^{-1} = \mathbf{B}\mathbf{\Xi}^{-1}\mathbf{A}. \]

Substituting our expressions for \(\mathbf{A}\), \(\mathbf{B}\) and \(\mathbf{\Xi}\), we get the following expression for the inverse:

\[\begin{align} \mSigma^{-1} &= \begin{bmatrix} \mathbf{I} & \mathbf{0} \\ -\mSigma_{22}^{-1}\mSigma_{21} & \mathbf{I} \end{bmatrix} \begin{bmatrix} (\mSigma/\mSigma_{22})^{-1} & \mathbf{0} \\ \mathbf{0} & \mSigma_{22} \end{bmatrix} \begin{bmatrix} \mathbf{I} & -\mSigma_{12}\mSigma_{22}^{-1} \\ \mathbf{0} & \mathbf{I}\end{bmatrix} \nonumber \\ &= \begin{bmatrix} (\mSigma/\mSigma_{22})^{-1} & \mathbf{0} \\ -\mSigma_{22}^{-1}\mSigma_{21}(\mSigma/\mSigma_{22})^{-1} & \mSigma_{22}^{-1} \end{bmatrix} \begin{bmatrix} \mathbf{I} & -\mSigma_{12}\mSigma_{22}^{-1} \\ \mathbf{0} & \mathbf{I}\end{bmatrix} \nonumber \\ &= \begin{bmatrix} (\mSigma/\mSigma_{22})^{-1} & -(\mSigma/\mSigma_{22})^{-1}\mSigma_{12}\mSigma_{22}^{-1} \\ -\mSigma_{22}^{-1}\mSigma_{21}(\mSigma/\mSigma_{22})^{-1} & \mSigma_{22}^{-1} -\mSigma_{22}^{-1}\mSigma_{21}(\mSigma/\mSigma_{22})^{-1}\mSigma_{12}\mSigma_{22}^{-1} \end{bmatrix} \label{eq:block_inv_1}. \end{align} \]

Note that had we done the block-diagonalisation the other way (zeroing-out the bottom-left block first) we would have obtained an expression in terms of \((\mSigma/\mSigma_{11})\):

\[\begin{equation} \label{eq:block_inv_2} \mSigma^{-1} = \begin{bmatrix} \mSigma_{11}^{-1} + \mSigma^{-1}_{11}\mSigma_{12}(\mSigma/\mSigma_{11})^{-1}\mSigma_{21}\mSigma_{11}^{-1} & -\mSigma_{11}^{-1}\mSigma_{12}(\mSigma/\mSigma_{11})^{-1} \\ -(\mSigma/\mSigma_{11})^{-1}\mSigma_{21}\mSigma_{11}^{-1} & (\mSigma/\mSigma_{11})^{-1} \end{bmatrix}. \end{equation} \]

Matrix Inversion Lemma

This soup of matrices above, giving the equation for the inverse of a partitioned matrix, might not immediately seem very helpful, however there are a few useful corollaries that we can draw from them. One example is the matrix inversion lemma:

We can now consider some examples where the matrix inversion lemma is useful in machine learning applications.

Example (Woodbury identity). Suppose that the matrix \(\rmM\) is partitioned such that:

• \(\rmA = \mSigma\) is an \(N \times N\) diagonal matrix,
• \(\rmB = \rmC^\top = \rmX \in \R^{N\times D}\), where we will assume that \(N \gg D\), and
• \(\rmD^{-1} = -\rmI\).

Substituting these terms into Equation \(\ref{eq:matrix_inversion}\) of the matrix inversion lemma gives us

\[\begin{equation} \label{eq:woodbury} (\mSigma + \rmX\rmX^\top)^{-1} = \mSigma^{-1} - \mSigma^{-1}\rmX(\rmI + \rmX^{\top}\mSigma^{-1}\rmX)^{-1}\rmX^{\top}\mSigma^{-1}. \end{equation} \]

The term \(\mSigma + \rmX\rmX^{\top}\) corresponds to doing a low-rank 11: since we have assumed that \(D\) is very small compared to \(N\); it may even be \(D = 1\) ^11[11] update to the matrix \(\mSigma\). The left hand side shows a naive re-inversion of this new matrix, which will take \(O(N^{3})\) time. However, if we keep a previously inverted \(\mSigma^{-1}\) in memory, we can quickly find the inverse of the updated matrix in \(O(D^{3})\) time by applying the right hand side of Equation \(\ref{eq:woodbury}\).

Crucially this might make the difference between an algorithm which is impractical to run on even the largest clusters, and one which is readily run on commodity hardware.
\(\blacksquare\)

The Sherman-Morrison formula is another useful example.

Example (Sherman-Morrison formula for rank one updates). This is just a special case of the above where \(D = 1\), corresponding to performing a rank-one update to an inverse matrix—this is common when data is coming in incrementally, and we’d like to re-use computation. This time, let our matrix \(\rmM\) be partitioned as follows:

• \(\rmA\) is some large matrix for which we already know the inverse \(\rmA^{-1}\),
• \(\rmB = \rvu\) (a column vector),
• \(\rmC = \rvv^\top\) (a row vector), and
• \(D = -1\) (a single scalar value). We then have

\[\begin{align*} (\rmA + \rvu\rvv^{\top})^{-1} &= \rmA^{-1} + \rmA^{-1}\rvu(-1 - \rvv^{\top}\rmA^{-1}\rvu)^{-1}\rvv^{\top}\rmA^{-1} \\ &= \rmA^{-1} - \frac{\rmA^{-1}\rvu\rvv^{\top}\rmA^{-1}}{1 + \rvv^{\top}\rmA^{-1}\rvu}. \end{align*} \]

\(\blacksquare\)

Another result that we can derived from the partitioned inverse formula is the matrix determinant lemma, which when used in conjunction with the Sherman-Morrison formula allows us to conveniently calculate low-rank updates to both a matrix inverse and determinant.

Similar to our statement of the Sherman-Morrison formula, we can highlight a special case of the matrix determinant lemma for computing a rank-one update to a determinant, where \(\rmB = \rvu\), \(\rmC = \rv^\top\) and \(D = -1\), giving

\[\vert \rmA + \rvu\rvv^\top \vert = (1 + \rvv^\top \rmA^{-1}\rvu)\vert \rmA \vert. \]

Gaussian Joint Distributions

A \(D\)-dimensional Gaussian random vector \(\mathbf{x} \sim \mathcal{N}(\vmu, \mSigma)\) can be partitioned into two sub-vectors \(\mathbf{x} = \begin{bmatrix}\mathbf{x}_1 & \mathbf{x}_{2}\end{bmatrix}^\top\) where \(\mathbf{x}_1 \in \mathbb{R}^N\), \(\mathbf{x}_{2}\in\mathbb{R}^M\) and \(N + M = D\).

In other words, we can interpret \(p(\mathbf{x})\) as the joint distribution \(p(\mathbf{x}_1, \mathbf{x}_2)\).

There are two common manipulations that we often want to perform on such a joint distribution:

1. marginalising out one of the sub-vectors e.g. \(p(\mathbf{x}_1) = \int p(\mathbf{x}_1, \mathbf{x}_2) d\mathbf{x}_2\),
2. computing a conditional e.g. \(p(\mathbf{x}_1 \vert \mathbf{x}_2)\).

For multivariate Gaussians, both of these can done in closed form.

Computing the marginal is particularly easy. Partitioning the parameters \(\vmu, \mSigma\) of our Gaussian random vector just as we partitioned \(\mathbf{x}\), we get:

\[\begin{equation} \label{eq:partitioned_gaussian_params} \vmu = \begin{bmatrix}\vmu_1 \\ \vmu_2 \end{bmatrix}, \hspace{20pt} \mSigma = \begin{bmatrix}\mSigma_{11} & \mSigma_{12} \\ \mSigma_{21} & \mSigma_{22}\end{bmatrix}, \hspace{20pt} \mathbf{\Lambda} \doteq \mSigma^{-1} = \begin{bmatrix} \mathbf{\Lambda}_{11} & \mathbf{\Lambda}_{12} \\ \mathbf{\Lambda}_{21} & \mathbf{\Lambda}_{22} \end{bmatrix}, \end{equation} \]

where \(\vmu_1 \in \mathbb{R}^N, \mSigma_{12} \in \mathbb{R}^{N\times M}\) etc., the marginals can be found by just reading off the corresponding parameters. Do also note that the partitions in \(\mLambda\) do not correspond straightforwardly to the partitions of \(\mSigma\); for instance, \(\mLambda_{11}\) is not merely the inverse of \(\mSigma_{11}\).

In other words, we simply have to select the relevant sub-vectors and sub-matrices of the joint Gaussian’s parameters \(\vmu\) and \(\mSigma\), respectively, to parametrise the marginals.

Finding the conditional is slightly more involved. We state the result here before proving it.

Proof idea: The Gaussian that we’re after has the form

\[\mathcal{N}(\mathbf{x}_1 \vert \vmu_{1\vert 2}, \mSigma_{1\vert 2}) = (2\pi)^{-d_1/2}\vert \mSigma_{1\vert 2}\vert^{-1/2} \exp\left(-\frac{1}{2}(\mathbf{x}_1 - \vmu_{1\vert 2})^\top \mSigma_{1\vert 2}^{-1}(\mathbf{x}_{1} - \vmu_{1\vert 2})\right). \]

In particular, the log-density of the above is a quadratic in both \(\rvx_1\) and \(\rvx_2\) 12: the \(\vmu_{1\vert 2}\) term must in general make reference to \(\rvx_2\); since \(\rvx_1\) is not independent from \(\rvx_2\). ^12[12] and the parameters we need to define \(p(\rvx_1\vert \rvx_2)\) (i.e. \(\vmu_{1\vert 2}\) and \(\mSigma_{1\vert 2}\)) appear in this quadratic.

Factorising the joint, \(\log p(\mathbf{x}_1, \mathbf{x}_2) = {\color{MidnightBlue}{\log p(\mathbf{x}_1 \vert \mathbf{x}_2)}} + \color{RedOrange}{\log p(\mathbf{x}_2)}\), we see that the first term on the right is the sum of a quadratic in both \(\rvx_1\) and \(\rvx_2\), and the second is a quadratic in \(\rvx_2\) only (ignoring constants and linear terms). Hence we can read off the parameters we need (the first two moments of the posterior, namely \(\color{MidnightBlue}{\vmu_{1\vert 2}}\) and \(\color{MidnightBlue}{\mSigma_{1\vert 2}}\)) from this first quadratic, discarding the rest.

Proof. We begin by explicitly writing the log-density of the joint:

\[\log p(\rvx_1, \rvx_2) = -\frac{1}{2}\log \vert \mSigma \vert - \frac{1}{2} \underbrace{ \begin{bmatrix}\mathbf{x}_1 - \vmu_1 \\ \mathbf{x}_2 - \vmu_2 \end{bmatrix}^\top \begin{bmatrix}\mSigma_{11} & \mSigma_{12} \\ \mSigma_{21} & \mSigma_{22}\end{bmatrix}^{-1} \begin{bmatrix}\mathbf{x}_1 - \vmu_1 \\ \mathbf{x}_2 - \vmu_2\end{bmatrix}}_{Q} + \text{const}, \]

and henceforth refer to this quadratic form as \(Q\). Since we want to end up with a term in \(\color{RedOrange}{\mathbf{x}_2}\), \(\color{RedOrange}{\vmu_{2}}\) and \(\color{RedOrange}{\mSigma_{22}}\) only, we block-diagonalise the partitioned matrix \(\mSigma\) (exactly as we did while deriving the partitioned matrix inverse formulae above); making sure to keep \(\mSigma_{22}\) on its own:

\[\begin{align*} Q &= \begin{bmatrix}\mathbf{x}_1 - \vmu_1 \\ \mathbf{x}_2 - \vmu_2 \end{bmatrix}^\top \begin{bmatrix}\mathbf{I} & \mathbf{0} \\ -\mSigma_{22}^{-1}\mSigma_{21} & \mathbf{I}\end{bmatrix} \begin{bmatrix}(\mSigma/\mSigma_{22})^{-1} & \mathbf{0} \\ \mathbf{0} & \mSigma_{22}^{-1} \end{bmatrix} \begin{bmatrix}\mathbf{I} & -\mSigma_{12}\mSigma_{22}^{-1} \\ \mathbf{0} & \mathbf{I}\end{bmatrix} \begin{bmatrix}\mathbf{x}_1 - \vmu_1 \\ \mathbf{x}_2 - \vmu_2\end{bmatrix} \\ &= \begin{bmatrix} \color{MidnightBlue}{\mathbf{x}_1 - \vmu_1 - \mSigma^{-1}_{22}\mSigma_{21}(\mathbf{x}_2 - \vmu_2)} \\ \color{RedOrange}{\mathbf{x}_2 - \vmu_2}\end{bmatrix}^\top \begin{bmatrix}\color{MidnightBlue}{(\mSigma/\mSigma_{22})^{-1}} & \mathbf{0} \\ \mathbf{0} & \color{RedOrange}{\mSigma_{22}^{-1}} \end{bmatrix} \begin{bmatrix}\color{MidnightBlue}{\mathbf{x}_1 - \vmu_1 - \mSigma_{12}\mSigma^{-1}_{22}(\mathbf{x}_2 - \vmu_2)} \\ \color{RedOrange}{\mathbf{x}_2 - \vmu_2}\end{bmatrix} \end{align*} \]

At this point, since the notation is already running over the lines, we will pre-emptively substitute \(\color{MidnightBlue}{\vmu_{1\vert 2} \doteq \vmu_1 + \mSigma_{12}\mSigma^{-1}_{22}(\mathbf{x}_2 - \vmu_2)}\). This should not be too much of a leap of faith, since we arrive at the factorised form with just a few more multiplications:

\[\begin{align*} Q &= \begin{bmatrix} \color{MidnightBlue}{\mathbf{x}_1 - \vmu_{1\vert 2}} \\ \color{RedOrange}{\mathbf{x}_2 - \vmu_2}\end{bmatrix}^\top \begin{bmatrix}\color{MidnightBlue}{(\mSigma/\mSigma_{22})^{-1}} & \mathbf{0} \\ \mathbf{0} & \color{RedOrange}{\mSigma_{22}^{-1}} \end{bmatrix} \begin{bmatrix}\color{MidnightBlue}{\mathbf{x}_1 - \vmu_{1 \vert 2}} \\ \color{RedOrange}{\mathbf{x}_2 - \vmu_2}\end{bmatrix} \\ &= \begin{bmatrix} \color{MidnightBlue}{\big(\mathbf{x}_1 - \vmu_{1\vert 2}\big)^\top(\mSigma/\mSigma_{22})^{-1}} & \color{RedOrange}{(\mathbf{x}_2 - \vmu_2)^\top\mSigma^{-1}_{22}}\end{bmatrix} \begin{bmatrix}\color{MidnightBlue}{\mathbf{x}_1 - \vmu_{1\vert 2}} \\ \color{RedOrange}{\mathbf{x}_2 - \vmu_2}\end{bmatrix} \\ &= {\color{MidnightBlue}{(\mathbf{x}_1 - \vmu_{1\vert 2})^\top(\mSigma/\mSigma_{22})^{-1}(\mathbf{x}_1 - \vmu_{1\vert 2})}} + {\color{RedOrange}{(\mathbf{x}_2 - \vmu_2)^\top\mSigma^{-1}_{22}(\mathbf{x}_2 - \vmu_2)}}. \\ \end{align*} \]

With the expression in the form we need, we can now just read off the covariance term for the posterior \(\mSigma_{1 \vert 2}\), which comes to the Schur complement of \(\mSigma\) wrt. \(\mSigma_{22}\):

\[\mSigma_{1\vert 2} = (\mSigma/\mSigma_{22}) \doteq \mSigma_{11} - \mSigma_{12}\mSigma_{22}^{-1}\mSigma_{21}, \]

and we already have \(\vmu_{1\vert 2}\) from the substitution made above.

While we omitted the linear and constant terms above (we only worked with the exponential terms in the factorisation of \(p(\vx_1 \vert \vx_2)\)), we can check our value for \(\mSigma_{1\vert 2}\) using the terms we omitted. Substituting \(\vert \mSigma \vert\) with the matrix determinant lemma (Equation \(\ref{eq:matrix_determinant_lemma}\)), we get

\[\begin{align*} (2\pi)^{(d_1 + d_2)/2}\vert \mSigma \vert^{1/2} &= (2\pi)^{(d_1 + d_2)/2}(\vert \mSigma / \mSigma_{22}\vert \vert \mSigma_{22}\vert)^{1/2} \\ &= {\color{MidnightBlue}{(2\pi)^{\frac{d_1}{2}}\vert \mSigma / \mSigma_{22}\vert^{\frac{1}{2}}}} {\color{RedOrange}{(2\pi)^{\frac{d_2}{2}}\vert \mSigma_{22}\vert^{\frac{1}{2}}}}, \end{align*} \]

as required.
\(\square\)

Working with Gaussians

As we saw above, when working with Gaussians we can often sidestep the need to do any integration when marginalising out terms, by instead considering the quadratic form in the exponent of the Gaussian that we are after, and matching coefficients to find the parameters of the distribution we’re after.

This turns out to be a useful pattern that is worth dwelling on. Considering the partitioning in Equation \(\ref{eq:partitioned_gaussian_params}\), the quadratic form in question is:

\[\begin{align} -\frac{1}{2}&(\rvx - \vmu)^{\top}\mSigma^{-1}(\rvx - \vmu) = \nonumber \\ &-\frac{1}{2}{\color{MidnightBlue}(\rvx_{1} - \vmu_1)^\top}\mLambda_{11}{\color{MidnightBlue}(\rvx_1 - \vmu_1)} - \frac{1}{2}{\color{MidnightBlue}(\rvx_1 - \vmu_1)^\top} \mLambda_{12}{\color{RedOrange}(\rvx_2 - \vmu_2)} \nonumber \\ &-\frac{1}{2}{\color{RedOrange}(\rvx_2 - \vmu_2)^\top} \mLambda_{21}{\color{MidnightBlue}(\rvx_1 - \vmu_1)} - \frac{1}{2}{\color{RedOrange}(\rvx_2 - \vmu_2)^\top} \mLambda_{22}{\color{RedOrange}(\rvx_2 - \vmu_2)} \label{eq:general_quadratic_form}. \end{align} \]

As an example of how we can use this, consider the above as a function of \(\rvx_1\). This would correspond to the quadratic form in the exponent of \(p(\rvx_1 \vert \rvx_2)\), and we could obtain the parameters of this Gaussian by inspection as we did above 13: This is often referred to as ‘completing the square’ ^13[13] .

The term in the exponent of a general Gaussian distribution, \(\gN(\rvx \vert \vmu, \mSigma)\) is:

\[\begin{equation} \label{eq:gaussian_exponent} -\frac{1}{2}(\rvx - \vmu)^{\top} \mSigma^{-1}(\rvx - \vmu) = -\frac{1}{2} \rvx^{\top} \mSigma^{-1}\rvx + \rvx^{\top}\mSigma^{-1}\vmu + \text{const}, \end{equation} \]

and by merely equating the general quadratic form to the right hand side of the above, we can equate the matrix of coefficients entering the second order term in \(\rvx_1\) to the inverse covariance \(\mSigma^{-1}\), and the coefficient of the linear term in \(\rvx_1\) to \(\mSigma^{-1}\vmu\), from which we can isolate \(\vmu\).

Returning to our example, we can derive the alternative forms for \(\vmu_{1\vert 2}\); that is, the mean of \(p(\rvx_1 \vert \rvx_2)\) that we gave in terms of the precision in Equations \(\ref{eq:conditional_mean_precision_1}\), \(\ref{eq:conditional_mean_precision_2}\). Treating \(\rvx_2\) as a constant conditioning term, we can pick out all the terms that are second order in \(\rvx_1\) in Equation \(\ref{eq:general_quadratic_form}\) to get

\[-\frac{1}{2}\rvx_1^\top \mLambda_{11}\rvx_1, \]

from which we can immediately conclude that the covariance of \(p(\rvx_1 \vert \rvx_2)\) is:

\[\mSigma_{1 \vert 2} = \mLambda_{11}^{-1}. \]

Now for the mean that we were originally after 14: remember, we needed to find \(\mSigma^{-1}\) first to recover the mean \(\vmu\) from the coefficient of the linear term ^14[14] , if we consider all the terms in Equation \(\ref{eq:general_quadratic_form}\) that are linear in \(\rvx_1\), we get 15: Using the fact that \(\mLambda_{21}^{\top} = \mLambda_{12}\) ^15[15]

\[\rvx_1^\top{\color{MidnightBlue}\left(\mLambda_{11}\vmu_1 - \mLambda_{12}(\rvx_2 - \vmu_2)\right)} \doteq \rvx_1^\top {\color{MidnightBlue}\mSigma^{-1}_{1\vert 2}\vmu_{1\vert 2}}. \]

Hence, left-multiplying both sides of the coefficients by \(\mSigma_{1\vert 2}\) gives us the conditional mean in the from of Equations \(\ref{eq:conditional_mean_precision_1}\), \(\ref{eq:conditional_mean_precision_2}\) that we were after:

\[\begin{align*} \vmu_{1\vert 2} &= \mSigma_{1 \vert 2}\big(\mLambda_{11}\vmu_1 - \mLambda_{12}(\rvx_2 - \vmu_2)\big) \\ &= \vmu_1 - \mLambda_{11}^{-1}\mLambda_{12}(\rvx_2 - \vmu_2). \end{align*} \]

\(\square\)

Of course, applying the Schur complement identities for the inverse of a partitioned matrix derived previously lets us write the precision matrices \(\mLambda_{11}\) and \(\mLambda_{12}\) in the above in terms of the covariance matrices, recovering the result from the previous section.

Linear Gaussian Systems

One problem that arises in many contexts is the following: suppose that we have two random variables. The first is a Gaussian random variable \(\rvx \in \R^{m}\), that might explain or cause something of interest. The second, \(\rvy \in \R^{d}\), is a Gaussian-noise corrupted observation of the thing of interest. If we also postulate that \(\rvx\) and the thing of interest are related by a linear function, then we have the following setup:

\[\begin{align*} p(\rvx) &= \gN(\rvx\vert \vmu, \mLambda^{-1}) \\ p(\rvy \vert \rvx) &= \gN(\rvy \vert \rmA\rvx + \rvb, \rmL^{-1}), \end{align*} \]

for some \(\rmA\in\R^{d\times m}\), \(\rvb\in\R^{d}\), and covariance \(\rmL^{-1}\) which is independent of \(\rvx\). This is called a linear Gaussian system.

We’re often interested in reversing this relationship: namely to ask, given a noisy observation \(\rvy\) of the thing I’m interested in, what are some probable values of \(\rvx\) that might have caused this observation?

It turns out that the Gaussian is particularly amenable to this sort of analysis, giving a convenient closed-form solution:

Proof idea: Since the prior and likelihood are both Gaussian, the joint \(p(\rvx, \rvy) = p(\rvx)p(\rvy \vert \rvx)\) is also a Gaussian whose mean and covariance parameters can be partitioned into the elements relevant to \(p(\rvx)\) and those relevant for \(p(\rvy \vert \rvx)\). We therefore first want to find the partitioned parameters of this joint.

From then, obtaining the parameters for \(p(\rvy)\) is a trivial case of reading off the relevant terms in the partitioned mean and covariance 16: see the marginalisation result above ^16[16] , while for \(p(\rvx \vert \rvy)\) we can complete the square and inspect the coefficients, as outlined in the previous section, to find the terms we need.

Proof. We begin by writing out the log of the joint, focusing on terms that depend on \(\rvx\) or \(\rvy\):

\[\begin{align} \log p(\rvx, \rvy) &= \log p(\rvx) + \log p(\rvy \vert \rvx) \nonumber \\ &= -\frac{1}{2}(\rvx - \vmu)^{\top}\mLambda(\rvx - \vmu) \nonumber \\ &\hspace{1em} -\frac{1}{2}(\rvy - \rmA\rvx - \rvb)^{\top}\rmL (\rvy - \rmA\rvx - \rvb) + \mathrm{const} \label{eq:log_joint}. \end{align} \]

We can verify that the joint is indeed Gaussian since this is a quadratic function of the components of \(\begin{bmatrix}\rvx & \rvy\end{bmatrix}^{\top}\). For conciseness later, let us define \(\rvz \doteq \begin{bmatrix}\rvx & \rvy\end{bmatrix}^{\top}\), where \(p(\rvz) = \gN(\rvz\vert \vmu_{\rvz}, \mSigma_{\rvz})\).

We first focus on finding the precision of this joint, which will be the coefficient of the second order term of \(\begin{bmatrix}\rvx & \rvy\end{bmatrix}^{\top}\). Expanding the brackets in the above, while keeping only the second order terms gives us:

\[\begin{align*} -\frac{1}{2}&\rvx^{\top}\big(\mLambda + \rmA^{\top}\rmL\rmA\big)\rvx - \frac{1}{2}\rvy^{\top}\rmL\rvy + \frac{1}{2}\rvy^{\top}\rmL\rmA\rvx + \frac{1}{2}\rvx^{\top}\rmA^{\top}\rmL\rvy \\ &= -\frac{1}{2}\begin{bmatrix}\rvx \\ \rvy \end{bmatrix}^{\top}\underbrace{\begin{bmatrix}\mLambda + \rmA^{\top}\rmL\rmA & -\rmA^{\top}\rmL \\ -\rmL \rmA & \rmL\end{bmatrix}}_{\rmR} \begin{bmatrix}\rvx \\ \rvy \end{bmatrix}. \end{align*} \]

Hence \(\rmR\) is the precision of \(p(\rvx, \rvy)\). To apply our results for the conditional of a joint Gaussian however, we need to invert the precision to find the covariance. Using the formula for the inverse of a partitioned matrix we derived earlier (Equation \(\ref{eq:block_inv_1}\)), we get

\[\begin{align*} \rmR^{-1} &= \begin{bmatrix}\mLambda^{-1} & \mLambda^{-1}\rmA^{\top} \\ \rmA\mLambda^{-1} & \rmL^{-1} + \rmA\mLambda^{-1}\rmA^{\top}\end{bmatrix} \\ &= \mSigma_{\rvz}. \end{align*} \]

To find the partitioned mean of the joint, we identify the linear terms in Equation \(\ref{eq:log_joint}\) giving

\[\rvx^{\top}\mLambda \vmu - \rvx^{\top}\rmA^{\top}\rmL\rvb + \rvy^{\top}\rmL\rvb = \begin{bmatrix}\rvx \\ \rvy\end{bmatrix}^{\top}\begin{bmatrix}\mLambda \vmu - \rmA^{\top}\rmL\rvb \\ \rmL \rvb\end{bmatrix}. \]

Recall that the linear term in the Gaussian exponent of the joint \(p(\rvz)\) is \(\rvz^{\top}\mSigma_{\rvz}^{-1}\vmu_{\rvz}\) (c.f. Equation \(\ref{eq:gaussian_exponent}\)). Matching terms,

\[\begin{align*} \mSigma_{\rvz}^{-1}\vmu_{\rvz} &= \begin{bmatrix}\mLambda \vmu - \rmA^{\top}\rmL\rvb \\ \rmL \rvb\end{bmatrix} \\ \vmu_{\rvz} &= \rmR^{-1} \begin{bmatrix}\mLambda \vmu - \rmA^{\top}\rmL\rvb \\ \rmL \rvb\end{bmatrix} \\ \end{align*} \]

Substituting our value for \(\rmR^{-1}\), we obtain

\[\vmu_{\rvz} = \begin{bmatrix}\vmu \\ \rmA\vmu + \rvb\end{bmatrix}. \]

With the partitioned parameters of the joint in hand, we can now trivially find the parameters of the marginal distribution \(p(\rvy)\) selecting the corresponding partitions:

\[\begin{align*} \E[\rvy] &= \rmA \vmu + \rvb \\ \Cov[\rvy] &= \rmL^{-1} + \rmA \mLambda^{-1}\rmA^{\top}. \end{align*} \]

For the parameters of the conditional \(p(\rvx \vert \rvy)\), we refer to Equations \(\ref{eq:conditional_mean_precision_1}\) and \(\ref{eq:conditional_cov}\), to obtain

\[\begin{align*} \E[\rvx \vert \rvy] &= (\mLambda + \rmA^{\top}\rmL\rmA)^{-1}\big(\rmA^{\top}\rmL(\rvy - \rvb) + \mLambda \vmu\big) \\ \Cov[\rvx \vert \rvy] &= (\mLambda + \rmA^{\top}\rmL\rmA)^{-1}. \end{align*} \]

\(\square\)

@article{robeyns2021gpmath
    author={Maxime Robeyns},
    title={The Venerable Gaussian},
    year={2023},
    url={https://maximerobeyns.com/gaussian},
}