I am working with singular Gaussian densities, which are formulated as follows.
Given $$x \sim \mathcal{N}(\mu,\Sigma),\text{ }(1)$$ where $\Sigma \in \mathbb{R}^{n \times n}$ has rank $p <n$ and $p>0$, then one version of its probablity density function is $$ \frac{1}{ \sqrt{ (2 \pi)^p |\Sigma|_+}} \exp{[-\frac{1}{2}(x-\mu)^\top \Sigma^\dagger (x-\mu) ]}, \text{ }(2) $$ where $|\bullet|_+$ denotes the product of non-zero eigenvalues of one matrix (pseudo determinant), and $()^\dagger$ denotes the Moore-Penrose inverse of a matrix. See equation $(2.1)$ in Srivastava's paper.
If now I have a conditional density and a marginal density as $$ Y| \theta \sim \mathcal{N}(A^\top\theta,\Sigma) , \text{ } \theta \sim \mathcal{N}(0,K), \text{ } (3) $$ where $\Sigma$ is the same as in $(1)$, while $K$ is non-singular and of dimension $m$. Then the question is what is the joint density of $$ \begin{bmatrix} Y\\ \theta \end{bmatrix}. $$
One possible choice is $$ \begin{bmatrix} Y\\ \theta \end{bmatrix} \sim \mathcal{N}(0, \begin{bmatrix} \Sigma + A^\top K A &A^\top K \\ K A & K \end{bmatrix}). \text{ }(4) $$ Then it can be found that the above joint density leads to the conditional density and marginal density in $(3)$, by using standard equations for Conditional distributions in Wiki.
However, the choice in $(4)$ seems to have one issue. If I use the formula that $p(Y,\theta)$ equals $$ p(Y|\theta)p(\theta), \text{ } (5) $$ and plug the densities into $(5)$ according to $(2)$ and $(3)$, equation $(5)$ supposes to lead to the same density as in $(4)$. However, the above equation leads to $$ \frac{1}{ \sqrt{ (2 \pi)^{p+m} |\Sigma|_+ |K| }} \exp{(\cdots)}, \text{ }(6) $$ while in general $$ |\begin{bmatrix} \Sigma + A^\top K A &A^\top K \\ K A & K \end{bmatrix}|_ + \not= |\Sigma|_+ |K|, $$ according to this question. Therefore, there seems to be a mismatch between the joint density in $(4)$ and the one produced by $(5)$. I am not sure which one leads to the correct joint density for $\begin{bmatrix}Y \\ \theta \end{bmatrix}$, or both of them are correct due to the non-uniqueness of density functions when covariance is singular?