How to prove that this joint distribution is Gaussian without using probability densities?

Question: I am wondering if there was a way to prove this result without using probability densities:

If $$\bf x \sim \mathcal N (m, P)$$ and $$\bf y \;|\; x \sim \mathcal N (Hx, R)$$, then $$\begin{pmatrix} \bf x \\ \bf y \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} \bf m \\ \bf Hm \end{pmatrix}, \begin{pmatrix} \bf P & \bf P H^{\top} \\ \bf HP & \bf H P H^{\top} + R \end{pmatrix} \right)$$

I came across the result here on slide 14 some time ago. A sketch of a proof in the univariate case can be found here.

The author of the slides uses the term "Gaussian densities" so the covariance matrices are non-singular. The author's choice of words may be harmless but it made me wonder about how a proof of the result would be when one cannot use probability densities.

Because $$(\bf y \;|\; x) \sim \mathcal N (Hx, R)$$, $$\bf y$$ can be written as $$\bf y = \bf Hx + \bf \epsilon$$ where $$\bf \epsilon \sim \mathcal N (0, R)$$ and independent with $$\bf x$$. It means $$\begin{pmatrix} \bf x \\ \bf \epsilon \end{pmatrix} \sim \mathcal N \left (\begin{pmatrix} \bf m \\ \bf 0 \end{pmatrix}, \begin{pmatrix} \bf P & 0 \\ \bf 0 &\bf R \end{pmatrix} \right )$$
Then $$\begin{pmatrix} \bf x \\ \bf y \end{pmatrix} = \begin{pmatrix} \bf I & 0 \\ \bf H &\bf I \end{pmatrix} \begin{pmatrix} \bf x \\ \bf \epsilon \end{pmatrix}$$ Following the fact that $$AY\sim N(A\mu, A\Sigma A') \text { given that } Y\sim N(\mu, \Sigma)$$ The results are there.