I have given two conditional distributions $p$ (assuming a conditional deterministic relationship): $$p(y_1|y_3)=\mathcal{N}(K_{N_1M}K_{MM}^{-1}y_3,0),$$ $$p(y_2|y_3)=\mathcal{N}(K_{N_2M}K_{MM}^{-1}y_3,0),$$ where $y_1\in \mathbb{R}^{N_1}$, $y_2\in \mathbb{R}^{N_2}$, $y_3\in \mathbb{R}^{M}$ and kernel matrices $K_{N_1M}\in \mathbb{R}^{N_1\times M}$, $K_{N_2M}\in \mathbb{R}^{N_2\times M}$, $K_{MM}\in \mathbb{R}^{M\times M}$, where each entry of the kernel matrices is a squared-exponential covariance function. $y$ are function values and e.g. the covariance matrix $K_{MM}$ corresponds to the input values $x_3$ belonging to the output $y_3$ Furthermore, we know $p(y_3)=\mathcal{N}(0,K_{MM})$.
The question now is how to get the joint distribution $p(y_1,y_2)$ under the assumption that we have conditional independence of $y_1,y_2$ given $y_3$.
What I know so far is that I can write: $$p(y_1,y_2)=\int p(y_1,y_2|y_3)p(y_3)dy_3 =\int p(y_1|y_3)p(y_2|_3)p(y_3)dy_3$$ where the second equation holds due to the additional assumption of conditional independence.
The solution (and here is the missing step I don't understand) should be $$p(y_1,y_2)=\mathcal{N}(0,\begin{bmatrix} K_{N_1M}K_{MM}^{-1}K_{MN_1} & K_{N_1M}K_{MM}^{-1}K_{MN_2} \\ K_{N_2M}K_{MM}^{-1}K_{MN_1} & K_{N_2M}K_{MM}^{-1}K_{MN_2} \end{bmatrix}).$$
I tried to solve the integral by plugging in the density functions of the normal distributions, but since the two deterministic probabilities don't have a density I don't know how to solve it.
Thanks a lot in advance!
It's a lot easier to think about the random vectors $y_i$ directly than about their density functions. You have $y_1=Ay_3$, $y_2=By_3$, and $y_3\sim N(0,M)$, for certain matrices $A$, $B$, and $M$. So $y_1$ and $y_2$ are jointly normal. The expectations of $y_1$ and $y_2$ are $AEy_3 = 0$ and $BEy_3= 0$, respectively. The covariance matrix of $y_1$ is $AMA^T$, that of $y_2$ is $BMB^T$ and the cross covariance matrix of $y_1$ and $y_2$ is $AMB^T$. In more detail, it is the expectation of $y_1 y_2^T $, that is, of $Ay_3 (By_3)^T = Ay_3 y_3^T B^T$. Since the covariance matrix of $y_3$ is $M$, $E y_3 y_3^T=M$, and thus the cross covariance matrix is $AMB^T$. So your big covariance matrix is $$\pmatrix{AMA^T&AMB^T\\BMA^T&BMB^T},$$ which is your expression when written out in full.