# Proof of the Derivation of the marginal and conditional Gaussian

Given a marginal Gaussian distribution for x and a conditional Gaussian distribution for y given x in the form $$p(x) = N(x|\mu, \Lambda^{-1})$$ $$p(y|x) = N(y|Ax + b, L^{-1})$$

the marginal distribution of y and the conditional distribution of x given y are given by $$p(y) = N(y|A\mu +b, L^{-1} + A\Lambda^{-1}A^T)$$ $$p(x|y) = N(x|\Sigma{A^TL(y-b) + \Lambda \mu}, \Sigma)$$

where $\Sigma = (\Lambda + A^TLA)^{-1}$.

This is taken from the textbook Pattern Recognition and Machine Learning.

The textbook has a proof but it is very brief and I could not fully understand the proof. Does anyone have a simple and detailed proof that can help me understand this fact? Thanks.

Let's write the RV $X$, $Y$ as $$X = \mu + \varepsilon _x \\ Y = AX + b + \varepsilon _y$$ with $\varepsilon _x \sim \mathcal N (0, \Lambda ^{-1})$, $\varepsilon _y \sim \mathcal N (0, L ^{-1})$. Now pluging in X in the second equation above gives

$$Y = A\mu + b + A\varepsilon _x + \varepsilon _y.$$

This is a linear combination of normal distributed random variables and as such itself normal distributed with expectation $A\mu +b$ and covariance matrix $A\Lambda ^{-1}A^T + L^{-1}$ ($var(AX)=Avar(X)A^T$). From this you get $p(y)$.

The second fact is a bit more complicated and involves some tedious calculation. From Bayes Theorem it follows that $$p(x|y) \propto p(y|x)p(x) \\ \propto \exp ((y-Ax-b)^TL(y-Ax-b) + (x-\mu)^T\Lambda(x-\mu)).$$

If you multiply everything out and then factor out $x$ you get to something proportional to $$\exp( (x-(\Lambda + A^TLA)^{-1}A^TL(y-b))^T(\Lambda + A^TLA)(x-(\Lambda + A^TLA)^{-1}AL(y-b)))$$ which is proportional to your given normal.

I would like to also give a algebraic response to this question showing that by multiplying the Gaussians together and integrating you can get this result. I did this using the covariance matrix defined as $L, \Lambda$ instead of $L^{-1}$ and $\Lambda^{-1}$ because that is my preference.

\begin{align} y|x \sim & N(Ax + b, L) \\ x \sim & N(\mu, \Lambda) \end{align} \begin{multline} y^T L^{-1} y - 2y^T L^{-1}(Ax + b) + (Ax+b)^T L^{-1} (Ax + b) + \\ x^T\Lambda^{-1} x - 2x^T\Lambda^{-1} \mu + \mu^T \Lambda^{-1} \mu = \\ y^T L^{-1} y- 2y^T L^{-1} Ax - 2y^T L^{-1} b + x^T A^T L^{-1} A x + 2x^T A^T L^{-1} b + b^T L^{-1} b \\ x^T\Lambda^{-1} x - 2x^T\Lambda^{-1} \mu + \mu^T \Lambda^{-1} \mu \end{multline}

As before make quadratic forms out of the $x$ terms to integrate them out: \begin{align} x^T ( \Lambda^{-1} + A^T L^{-1} A) x - 2x^T ( \Lambda^{-1} \mu + A^T L^{-1} y - A^T L^{-1} b) + \mu^T \Lambda^{-1} \mu + b^T L^{-1} b \end{align} \begin{align} V &= [A^TL^{-1} A + \Lambda^{-1} ]^{-1} \\ u &= [ \Lambda^{-1} \mu + A^T L^{-1} y - A^T L^{-1} b]\\ h & = Vu \\ h^TV^{-1}h &= u^T V u \\ c & = \mu^T \Lambda^{-1} \mu + b^T L^{-1} b \\ \end{align}

$x$ is Gaussian and integrates to its normalizing constant. The term $c- K$, $K = u^TVu$ factors into the $y$ exponential giving, \begin{align} y^T L^{-1} y -2y^TL^{-1} b + c - K \end{align}

And,

\begin{multline} u^T V u = \\ [ \Lambda^{-1} \mu + A^T L^{-1} y - A^T L^{-1} b]^T V [ \Lambda^{-1} \mu + A^T L^{-1} y - A^T L^{-1} b] = \\ \mu^T \Lambda^{-1} V \Lambda^{-1} \mu + y^T L^{-1} A V A^T L^{-1} y + b^T L^{-1} AV A^T L^{-1} b \\ + 2y^T L^{-1} A V \Lambda^{-1} \mu - 2\mu^T \Lambda^{-1} V A^T L^{-1} b - 2y^T L^{-1} A V A^T L^{-1} b \end{multline}

Together this becomes:

\begin{multline} y^T [ L^{-1} - L^{-1} A V A^T L^{-1} ] y - 2y^T [L^{-1} b + L^{-1} A V \Lambda^{-1} \mu - L^{-1} A V A^T L^{-1} b] + \\ \mu^T[ \Lambda^{-1} - \Lambda^{-1} V \Lambda^{-1} ]\mu + b^T[L^{-1} - L^{-1} A V A^T L^{-1}] b - 2 b^T L^{-1} A V \Lambda^{-1} \mu \end{multline}

Which simplifies slightly,

$$y^T R^{-1} y - 2y^T g + b^T R^{-1} b + \mu^T[ \Lambda^{-1} - \Lambda^{-1} V\Lambda^{-1}]\mu - 2 b^T L^{-1} A V \Lambda^{-1} \mu$$ \begin{align} R^{-1} = L^{-1} - L^{-1} A V A^T L^{-1} = [L + A \Lambda A^T]^{-1} \\ g = [L^{-1} b + L^{-1} A V \Lambda^{-1} \mu - L^{-1} A V A^T L^{-1} b] \end{align} Furthermore, $Rg = A\mu +b$:

\begin{align} R g &= R [L^{-1} b + L^{-1} A V \Lambda^{-1} \mu - L^{-1} A V A^T L^{-1} b] \\ &= R[ (L^{-1} - L^{-1} A VA^T L^{-1}) b + L^{-1} A V \Lambda^{-1} \mu ] \label{amu1} \\ & = b + R L^{-1} A V \Lambda^{-1} \mu \label{amu2} \end{align}

After canceling out the $b$ term on both sides I get I solve for the other term, equating coefficients:

\begin{align} R L^{-1} A V \Lambda \mu = & [L + A \Lambda A^T] L^{-1} A V \Lambda^{-1} \mu = A \mu \label{amu3} \\ = & A + A\Lambda A^T L^{-1} A = A [\Lambda^{-1} + A^T L^{-1} A] \Lambda \label{amu4} \\ & \text{(using the definition of $V$)} \notag \\ \end{align}

\begin{align} I + \Lambda A^T L^{-1} A = & I + A^T L^{-1} A \Lambda \label{amu5}\\ \Lambda A^T L^{-1} A = & \Lambda A^T L^{-1} A \label{amu6} \end{align} And therefore $Rg = A\mu + b$ and $R^{-1} A \mu = L^{-1} AV \Lambda^{-1} \mu.$

Simplifying equation the equation for $y$ I get, $$y^T R^{-1} y - 2y^TR^{-1} (A\mu + b)+ b^TR^{-1} b + \mu^T[ \Lambda^{-1} - \Lambda^{-1} V \Lambda^{-1} ]\mu - 2 b^T L^{-1} A V \Lambda^{-1} \mu.$$

Notice, $$\Lambda - \Lambda^{-1} V \Lambda^{-1} = A^T R^{-1} A$$

To show this: \begin{align} A^T R^{-1} A & = \Lambda^{-1} - \Lambda^{-1} V \Lambda^{-1} \\ & = \Lambda^{-1} - \Lambda^{-1} I ( I \Lambda I + F^{-1}) I \Lambda^{-1} \end{align}

$F^{-1} = A^TL^{-1} A$,

\begin{align} \Lambda^{-1} - \Lambda^{-1} I ( I \Lambda^{-1} I + F^{-1}) I \Lambda^{-1} & = [\Lambda + F] ^{-1} \end{align}

By the Woodbury identity.

\begin{align} A^T R^{-1} A & = [\Lambda + F] ^{-1}\\ A^{-1} R A^{-1} & = \Lambda + F \\ R & = A \Lambda A^T + L \end{align} Since both sides are equal, the `sandwich' term is equal to $A^T R^{-1} A$ I get the following:

$$y^T R^{-1} y - 2y^TR^{-1} (A\mu + b)+ b^TR^{-1} b + \mu^T A^T R^{-1} A \mu - 2 b^T R^{-1} A\mu.$$

Which gives the results that $y$ is distributed normally.

$$y \sim \mathcal{N} (A\mu + b, R)$$

The other part is much easier, and the algebra is not very tedious. The first the long equation given for $y$ can be simplified when conditioning on $x$ since everything not involving $x$ can be subsumed into the constant of proportionality. This leaves only,

$$-2x^T A^T L^{-1} y + x^T A^T L^{-1} x + 2 x^T A^T L^{-1} b + X^T \Lambda^{-1} x - 2x^T \Lambda ^{-1} \mu$$

Then refactoring this term,

$$x^T [ A^T L^{-1} A + \Lambda^{-1} ] x - 2x^T [A^T L^{-1} (y - b) + \Lambda^{-1} \mu] + K$$ Where $K$ is some constant. Call $[ A^T L^{-1} A + \Lambda^{-1} ] = \Sigma^{-1}$. Therefore, $$x \sim \mathcal{N} ( \Sigma [ A^T L^{-1}(y - b) + \Lambda^{-1} \mu], \Sigma)$$

Remember I defined my covariance matrix to the inverse of what was given in the question.