I am currently reading the paper entitled "Variational Auto-Encoded Deep Gaussian Processes" by Dai et al, a copy of which may be found here.
The paper proposes a stacking of Gaussian Process Latent Variable Models in an ANN like fashion and introducing a variational distribution paramaterised by a Multilayer Perceptron to aid tractability.
Before asking my question, I shall introduce some preliminaries for your convenience(paraphrased from the paper). If this is too long, I am happy to shorten it and rely more on the paper as a reference.
The marginal likelihood of the DGP model is given in Equation 4 in the paper as
$p(\displaystyle \mathbf{Y}) = \int p(\mathbf{Y} \mid \mathbf{X}_{1}) \prod_{l=2}^{L} p(\mathbf{X}_{l-1} \mid \mathbf{X}_{l})p(\mathbf{X}_{L})d\mathbf{X}_{1} \ldots d\mathbf{X}_{L}$
where $\mathbf{Y} \in \mathbb{R}^{N \times D}$ is a matrix of observed data, $L$ is the number of layers of latent variables, for layer $l$, $\mathbf{X}_{l} \in \mathbb{R}^{N \times Q_{l}}$ is a latent space representation of feature dimensionality $Q_{l} < Q_{l-1}$.
For layer 1, the input is $\mathbf{Y}$ so for $l=2$, $\mathbf{X}_{l-1} = \mathbf{Y}$
A variational lower bound of the above log marginal distribution is given as per Equation 5 in the paper
$\displaystyle \mathcal{L} = \sum_{l=1}^{L} \mathop{\mathbb{E}}\big[\log p(\mathbf{X}_{l-1}\mid\mathbf{X}_{l})\big]_{q(\mathbf{X}_{l-1})q(\mathbf{X}_{l})} + \sum_{l=1}^{L-1} H(q(\mathbf{X}_{l})) - KL(q(\mathbf{X}_{L}) \mid\mid p(\mathbf{X}_{L}))$
Where $H(.)$ is the Shannon Entropy and $KL(. \mid\mid .)$ is the Kullback Leibler divergence. To make inference tractable, the variational distribution $q(.)$ is introduced above, as is defined as follows
$\displaystyle q(\mathbf{X}_{l}) = \prod_{n=1}^{N} \mathcal{N}(\mathbf{x}^{(n)}_{l} \mid \boldsymbol\mu^{(n)}_{l}, \boldsymbol\Sigma^{(n)}_{l})$
where $\boldsymbol\mu_{l}(.)$ is the output of a Multilayer Perceptron taking $\mathbf{X}_{l}$ as input and $\boldsymbol\Sigma_{l}$ is the posterior variance, assumed to be diagonal and the same over all datapoints.
However, at this point the variational lower bound is still intractable, due to the expectation in the first term of $\mathcal{L}$. As such, auxiliary variables $\mathbf{U}_{l} \in \mathbb{R}^{M \times Q_{l}}$ are introduced, as are noise free observations $\mathbf{F}_{l} \in \mathbb{R}^{N \times Q_{l-1}}$ and the first term of $\mathcal{L}$ reformulated as follows(though at this point, it it not clear to me how $\mathbf{F}_{l}$ differs from $\mathbf{X}_{l-1}$), see Equation 10 in the paper
$\displaystyle\mathop{\mathbb{E}}\big[\log p(\mathbf{X}_{l-1}\mid\mathbf{X}_{l})\big]_{q(\mathbf{X}_{l-1})q(\mathbf{X}_{l})} \geq \mathop{\mathbb{E}}\big[ \log p(\mathbf{X}_{l-1} \mid \mathbf{F}_{l}) - KL(q(\mathbf{U}_{l} \mid \mathbf{X}_{l-1}) \mid\mid p(\mathbf{U}_{l})) \big]_{p(\mathbf{F}_{l} \mid \mathbf{U}_{l}, \mathbf{X}_{l})q(\mathbf{U}_{l} \mid \mathbf{X}_{l-1})q(\mathbf{X}_{l-1})q(\mathbf{X}_{l})}$
Finally, the authors give a distributed(in the parallel computation sense) form of the variational lower bound $\mathcal{L}$ with the first term taking the following form
$\displaystyle Tr\big(\mathop{\mathbb{E}}\big[ \mathbf{X}^{T}_{l-1}\mathbf{X}_{l-1} \big]_{q(\mathbf{X}_{t-1})}\big) = \sum_{n=1}^{N} (\boldsymbol\mu^{(n)}_{l-1})^{T}\boldsymbol\mu_{l-1}^{(n)} + Tr\big( \boldsymbol\Sigma_{l-1}^{(n)} \big)$
and the second term taking the following form
$\displaystyle Tr\big( \boldsymbol\Lambda_{l}^{-1}\boldsymbol\Psi_{l}^{T} \mathop{\mathbb{E}}\big[\mathbf{X}_{l-1}\mathbf{X}_{l-1}\big]_{q(\mathbf{X}_{l-1})} \boldsymbol\Psi_{l} \big) = Tr\big(\boldsymbol\Lambda_{l}^{-1}(\boldsymbol\Psi_{l}^{T}\mathbf{R}_{l-1}^{T})(\boldsymbol\Psi_{l}\mathbf{R}_{l-1})\big) + Tr\bigg(\boldsymbol\Lambda_{l}^{-1} \bigg(\sum_{n=1}^{N}\boldsymbol\Psi_{l}^{(n)}\alpha_{l-1}^{(n)}\bigg)\bigg(\sum_{n=1}^{N}\boldsymbol\Psi_{l}^{(n)}\alpha_{l-1}^{(n)}\bigg)^{T}\bigg)$
where $\boldsymbol\Lambda_{l} = \mathbf{K}_{\mathbf{U}_{l}\mathbf{U}_{l}} + \mathop{\mathbb{E}}\big[\mathbf{K}^{T}_{\mathbf{F}_{l}\mathbf{U}_{l}}\mathbf{K}_{\mathbf{F}_{l}\mathbf{U}_{l}}\big]_{q(\mathbf{X}_{l})}$ for covariance matrix $\mathbf{K}_{<...>}$ generated by a covariance kernel such as the exponentiated quadratic. Similarly, $\boldsymbol\Psi_{l} = \mathop{\mathbb{E}}\big[\mathbf{K}_{\mathbf{F}_{l}\mathbf{U}_{l}}\big]_{q(\mathbf{X}_{l})}$.
Additionally, $\mathbf{R}_{l-1} = \big[(\boldsymbol\mu^{(1)}_{l-1})^{T} \dots (\boldsymbol\mu^{(N)}_{l-1})^{T} \big]$ and $\alpha_{l-1}^{(n)} = \sqrt{Tr\big(\boldsymbol\Sigma_{l-1}^{(n)}\big)}$ and $\mathbf{A}_{l=1} = diag\big(\alpha_{l-1}^{(1)} \dots \alpha_{l-1}^{N}\big)$
Finally, the second term of the above distributed form is obtained by making the following observation
$\mathop{\mathbb{E}}\big[ \mathbf{X}_{l-1}\mathbf{X}_{l-1}^{T} \big] = \mathbf{R}_{l-1}^{T}\mathbf{R}_{l-1} + \mathbf{A}_{l-1}\mathbf{A}_{l-1}$
Points of Confusion - Question(s)
Given the above formulation, there are a few computational issues that I am at present not entirely clear on.
Firstly, how in general does one take the expectation of a covariance matrix generated by a kernel such as the exponentiated quadratic?
For example, in the above formulation, the computation of the quantity $\boldsymbol\Psi_{l}$ is not clear to me. By the definitions given in the paper, to evaluate this quantity one takes the expectation w.r.t. the variational distribution $q(\mathbf{X}_{l})$, which itself is defined to be the product of Normal PDF evaluations of the GP latent variables given their MLP encodings - a scalar.
It is not clear to me how to use this to take the expectation of an arbitrary covariance/cross-covariance matrix $\mathbf{K} \in \mathbb{R}^{N \times M}$. As such, the evaluation of the quantity defined by the expectation in the second term of $\boldsymbol\Lambda_{l}$ is also unclear to me as we are again attemnpting to take the expectation of a cross-covariance matrix w.r.t. the variational distribution $q(\mathbf{X}_{l})$.
My central question is, given a formulation like the one set out above, how does one handle these expectation terms?
In addition, if anybody is familiar with this work and/or area of research, it would also be of great assistance to clarify the form of $\mathbf{F}_{l}$ and $\mathbf{U}_{l}$ for an arbitrary layer. To my understanding, for $l=1$, the data layer, $\mathbf{F}$ is a subset of the data points(removed from the dataset) and $\mathbf{U}$ the corresponding latent variables. However, what about arbitrary $l$?
Any assistance would be greatly appreciated.
