Conditional Multivariate Gaussian Identity I'm trying to verify the form of a multivariate Gaussian provided in a paper I'm reading. It should be pretty elementary.
Let $Y=X+\varepsilon$ where $X\sim N(0,C)$ and $\varepsilon\sim N(0,\sigma^2\mathbf{I})$. The authors then claim that
$$
X|Y,C,\sigma^2 \sim N(\mu,\Sigma),
$$
where
$$
\mu := C(C+\sigma^2\mathbf I)^{-1}Y\\
\Sigma:=\sigma^2C(C+\sigma^2\mathbf I)^{-1}.
$$
My first thought was to consider the joint distribution
$$
\begin{pmatrix}
X\\
Y
\end{pmatrix}\sim N\Big(\begin{pmatrix}
0\\
0
\end{pmatrix},\begin{pmatrix}
C & C\\
C^\top & \sigma^2\mathbf I+C
\end{pmatrix}\Big)
$$
and apply the conditional Gaussian identities. Unfortunately this approach gives me the right $\mu$, but I can't see how their form of $\Sigma$ comes about. Any thoughts?
 A: This is a correct representation of the conditional variance.
Since
$$\begin{pmatrix}
X\\
\epsilon
\end{pmatrix}\sim N\Big(\begin{pmatrix}
0\\
0
\end{pmatrix},\begin{pmatrix}
C & \mathbf O\\
\mathbf O & \sigma^2\mathbf I
\end{pmatrix}\Big)$$
and
$$\begin{pmatrix}
X\\
Y
\end{pmatrix} =  \begin{pmatrix}
\mathbf 1^\text{T} &\mathbf 0^\text{T} \\
\mathbf 1^\text{T} &\mathbf 1^\text{T}
\end{pmatrix}\begin{pmatrix}
X\\
\epsilon
\end{pmatrix}$$
the distribution of $\begin{pmatrix}
X\\
Y
\end{pmatrix}$ is
$$\begin{pmatrix}
X\\
Y
\end{pmatrix}\sim N\Big(\begin{pmatrix}
0\\
0
\end{pmatrix},\underbrace{\begin{pmatrix}
\mathbf 1^\text{T} &\mathbf 0^\text{T} \\
\mathbf 1^\text{T} &\mathbf 1^\text{T}
\end{pmatrix}\begin{pmatrix}
C & \mathbf O\\
\mathbf O & \sigma^2\mathbf I
\end{pmatrix}\begin{pmatrix}
\mathbf 1 &\mathbf 1 \\
\mathbf 0 &\mathbf 1
\end{pmatrix}}_{\begin{pmatrix}
C & C\\
C & \sigma^2\mathbf I
\end{pmatrix}
}\Big)$$
indeed. With
$$\mathbb E[X|Y] = 0 + C (C+\sigma^2\mathbf I)^{-1}
Y $$
and
$$\text{var}(X|Y) = C - C (C+\sigma^2I)^{-1} C
$$
Applying the Woodbury matrix inversion lemma
$$(A+B)^{-1}=A^{-1}-A^{-1}(B^{-1}+A^{-1})^{-1}A^{-1}$$
one gets that
\begin{align*} C - C (C+\sigma^2I)^{-1} C &= C - C (C^{-1}-
C^{-1}(C^{-1}+\sigma^{-2}\mathbf I)^{-1}C^{-1})C\\ &= C - C
+(C^{-1}+\sigma^{-2}\mathbf I)^{-1}\\ &= 
(C^{-1}\mathbf I+\sigma^{-2}C^{-1}C)^{-1}\\ &= \sigma^2 C (\sigma^2\mathbf I+C)^{-1}
\end{align*}
The apparent lack of symmetry in the expression may sound suspicious but actually$$C (\sigma^2\mathbf I+C)^{-1} = (\sigma^2\mathbf I+C)^{-1} C$$
