The following problem comes from studying the conditional expectation of a multivariate normal distribution. Let $n\ge2$ be an integer, and let $\Sigma$ be a positive semidefinite, symmetric $n\times n$ matrix of real numbers partitioned as $$\Sigma=\begin{pmatrix}\Sigma_{a,a}&\Sigma_{a,b}\\\Sigma_{b,a}&\Sigma_{b,b}\end{pmatrix},$$ where $\Sigma_{a,a}$ is $1\times1,$ $\Sigma_{a,b}$ is $1\times(n-1)$ and $\Sigma_{b,b}$ is $(n-1)\times(n-1).$ Is it true that $\Sigma_{a,b}=\Sigma_{a,b}\Sigma_{b,b}^+\Sigma_{b,b}?$ (Here, $\Sigma_{b,b}^+$ is the Moore-Penrose pseudoinverse.)
In the CrossValidated post "Conceptual proof that conditional of a multivariate Gaussian is multivariate Gaussian", someone claims this result. The result is clearly true if $\Sigma_{b,b}$ is invertible, in which case $\Sigma_{b,b}^+=\Sigma_{b,b}^{-1}.$ In addition, I have tried two examples, $\Sigma=0$ and $$\Sigma=\begin{pmatrix}\!\!\begin{array}{c|cc}1&1&0\\\hline1&1&0\\0&0&0\end{array}\!\!\end{pmatrix},$$ and in both cases we have $\Sigma_{a,b}-\Sigma_{a,b}\Sigma_{b,b}^+\Sigma_{b,b}=0,$ as desired. However, I could not prove the result in generality using the definition or the properties listed in the Wikipedia page.
In the context of this problem, we partition a multivariate Gaussian vector into two subvectors $X_a$ and $X_b$. The goal is to find a matrix $C$ such that $Z:=X_a- C X_b$ is uncorrelated with $X_b$, so that the equality $$\Sigma_{a,b}=C\,\Sigma_{b,b}\tag1$$ holds. The general solution is to take $$C:=\Sigma_{a,b}\Sigma_{b,b}{}^+,$$ where $\Sigma_{b,b}{}^+$ is the Moore-Penrose inverse of $\Sigma_{b,b}$.
The reason why this works, even when $\Sigma_{b,b}$ is not of full rank, depends crucially on the fact that we are dealing with multivariate Gaussians, in which case $\Sigma_{b,b}$ and $\Sigma_{a,b}$ have a special form. Recall that every multivariate Gaussian vector is an affine transformation of some vector $Z$ of independent standard Gaussians. We can then write the subvectors $X_a$ and $X_b$ in the form $X_a = AZ + \mu_a$, $X_b = BZ + \mu_b$, with $A$ and $B$ matrices of constants. Since the covariance matrix of $Z$ is the identity, we calculate $\Sigma_{a,b} = AB^T$ and $\Sigma_{b,b}=BB^T$.
Using properties of Moore-Penrose inverses, we find $\Sigma_{b,b}{}^+=(BB^T)^+=(B^T)^+B^+$ and verify (1): $$ C\,\Sigma_{b,b}=\Sigma_{a,b}\Sigma_{b,b}{}^+\Sigma_{b,b} =AB^T(B^T)^+\underbrace{B^+BB^T}_{B^T} =A\underbrace{B^T(B^T)^+B^T}_{B^T}=AB^T=\Sigma_{a,b}. $$