Conditional Deep Gaussian Processes: Empirical Bayes Hyperdata Learning

Abstract

It is desirable to combine the expressive power of deep learning with Gaussian Process (GP) in one expressive Bayesian learning model. Deep kernel learning showed success as a deep network used for feature extraction. Then, a GP was used as the function model. Recently, it was suggested that, albeit training with marginal likelihood, the deterministic nature of a feature extractor might lead to overfitting, and replacement with a Bayesian network seemed to cure it. Here, we propose the conditional deep Gaussian process (DGP) in which the intermediate GPs in hierarchical composition are supported by the hyperdata and the exposed GP remains zero mean. Motivated by the inducing points in sparse GP, the hyperdata also play the role of function supports, but are hyperparameters rather than random variables. It follows our previous moment matching approach to approximate the marginal prior for conditional DGP with a GP carrying an effective kernel. Thus, as in empirical Bayes, the hyperdata are learned by optimizing the approximate marginal likelihood which implicitly depends on the hyperdata via the kernel. We show the equivalence with the deep kernel learning in the limit of dense hyperdata in latent space. However, the conditional DGP and the corresponding approximate inference enjoy the benefit of being more Bayesian than deep kernel learning. Preliminary extrapolation results demonstrate expressive power from the depth of hierarchy by exploiting the exact covariance and hyperdata learning, in comparison with GP kernel composition, DGP variational inference and deep kernel learning. We also address the non-Gaussian aspect of our model as well as way of upgrading to a full Bayes inference.

1. Introduction

The deep Gaussian process [1] is a Bayesian learning model which combines both the expressive power of deep neural networks [2] and calibrated uncertainty estimation. The hierarchical composition of Gaussian Processes (GPs) [3] is the origin of expressiveness, but also renders inference intractable, as the marginalization of GPs in the stage of computing evidence is not analytically possible. Expectation propagation [4,5] and variational inference [6,7,8,9] are approximate inference schemes for DGP. The latter has issues of posterior collapse, which turns DGP into a GP with transformed input. References [8,9] address this issue and compositional freedom [10] in such hierarchical learning. Nevertheless, inferential challenges continue to slow the adoption of DGP.

Despite challenges, there has been progresses in understanding this seemingly simple yet profound model. In the case where the GPs in the hierarchy are zero-mean, DGP exhibits pathology, becoming a constant function as the depth increases [11]. Using the fact that the exponential covariance function is strictly convex, references [12,13] studied the conditional statistics for squared distance in function space, suggesting region in hyperparameter space to avoid the pathology. Recently, reference [14] showed the connection between DGP and a deep neural network with bottlenecked layers, and reference [15] suggested that a DGP with a large width may collapse back to a GP.

Others have found ways to work around the challenges of DGPs. The deep kernel learning proposed in [16] gained the Bayesian character of GP and the expressive power of a deep neural network without encountering intractability, as the learning of weight parameters, treated as kernel hyperparameters, is an empirical Bayes. Similar ideas also appeared in [17,18]. Hyperparameter learning in [16] was performed through marginal likelihood, which can in principle prevent overfitting due to the built-in competition between data fitting and model complexity [3]. However, [19] suggested that the lack of Bayesian character in the deep feature extracting net might still result in overfitting if the network has too many parameters.

Here, we propose a conditional DGP model in which the intermediate GPs (all but the exposed GP) in the hierarchical composition are conditioned on a set of hyperdata. These hyperdata are inspired by the inducing points in sparse GP [20,21,22], but they are hyperparameters, not random variables. The conditional DGP is motivated by the expressive power and Bayesian character of DGP [1], and the deep kernel learning with an objective in marginal likelihood [16]. Due to the conditioning on the hyperdata, the intermediate GPs can be viewed as collections of random feature functions centered around the deterministic conditional mean. Thus, the intermediate GPs become approximately deterministic functions when the hyperdata are sufficiently dense. Besides, lifting the intermediate GPs from being zero mean might help avoid pathology too. Mathematically, we defined a marginal prior for the conditional DGP; i.e., all intermediate GPs are marginalized, which assures the Bayesian character when dealing with the feature functions. We then use the moment matching method to approximate the non-Gaussian marginal prior as a GP [23], which connects with observed data and allows the marginal likelihood objective. It should be stressed that the effective kernel depends on the conditional mean and conditional covariance in feature function via the hyperdata, which are optimized in the spirit of empirical Bayes [24]. In the implementation, the hyperdata supporting each intermediate GP are represented as a neural network function, $u={\mathrm{nn}}_{\mathbf{w}}\left(z\right)$ with u and z being the output and input of hyperdata, similarly to the trick used in modeling the mean and variance for data in the variational autoencoder [25].

The paper is organized as follows. Section 2 gives a short survey of the current literature on deep probabilistic models; the usage of moment matching in approximate inference; and the inducing points in GP and DGP. Background on mathematical models of GP and DGP, the marginal prior for DGP and the moment matching method are introduced in Section 3. The conditional DGP with SE kernel in the exposed layer, its mathematical connection with deep kernel learning, the parameter learning and the non-Gaussian aspect, are described in Section 4. A preliminary demonstration on extrapolating two time-series data is in Section 5, followed by a discussion in Section 6.

2. Related Work

In the literature on deep probabilistic models, [26] proposed the conditional neural process in which the mean and variance functions are learned from the encoded representation of context data in a regression setup for target data. Deep Gaussian processes (DGPs) constitute one family of models for composition functions by conditioning input to a GP on the output of another GP [1]. A similar idea appeared in the works of warped GP [27,28]. The implicit process in [29] is a stochastic process embedding the Gaussian distribution into a neural network. Solutions of stochastic differential equation driven by GP are also examples of composite processes [30]. Variational DGP casts the inference problem in terms of optimizing ELBO [6] or EP [5]. However, the multi-modalness of DGP posterior [10,23] may arise from the fact that the hidden mappings in intermediate layers are dependent [9]. Inference schemes capable of capturing the multi-modal nature of DGP posterior were recently proposed by [8,9]. Depth of neural network models and the function expressivity were studied in [31,32], and uncertainty estimates were investigated in [33]. DGP in weight space representation and its variational Bayesian approach to DGP inference were introduced in [34], which were based on the notion of random feature expansion of Gaussian [35] and arcsine [36] kernels. Deep hierarchical SVMs and PCAs were introduced in [37].

Moment matching is a way to approximate a complex distribution with, for instance, a Gaussian by capturing the mean and the second moment. Reference [38] considered a GP regression with uncertain input, and replaced the non-Gaussian predictive distribution with a Gaussian carrying the matched mean and variance. Expectation propagation, in [4], computed the vector of mean and variance parameters of non-Gaussian posterior distributions. Reference [21] approximated the distribution over unseen pixels as a multivariate Gaussian with matched mean and covariance. Moment matching is also extensively applicable to comparing two distributions [39] where the embedded means in RKHS are computed. In generative models, the model parameters are learned from comparing the model and data distributions [40].

Inducing points are an important technique in sparse GP [20,22,41,42] and DGP. In addition to being locally defined as a function’s input and output, [43] introduced a transformation to form a global set of inducing features. One popular transformation uses the basis of Gaussian so that one can recover the local inducing points easily [43]. Transformation using the basis of spherical harmonic functions in [44] allows orthogonal inducing features and connects with the arcsine kernels of Bayesian deep neural network [45]. Reference [46] employed the inter-domain features in DGP inference. Recently, [47] proposed a method to express the local inducing points in the weight space representation. All the methods cited here treated the inducing points or features in a full Bayes approach, as they are random variables associated with an approximate distribution [24].

3. Background

Here, we briefly introduce the notions of the Gaussian process as a model for random continuous function $f\left(\mathbf{x}\right):{\mathbb{R}}^d\mapsto \mathbb{R}$. A deep Gaussian process [1] is a hierarchical composition of Gaussian processes for modeling general composite function ${\mathbf{f}}_L\circ {\mathbf{f}}_{L-1}\circ \cdots {\mathbf{f}}_2\circ {\mathbf{f}}_1\left(\mathbf{x}\right)$ where the bold faced function ${\mathbf{f}}_1:{\mathbb{R}}^d\mapsto {\mathbb{R}}^{H_1}$ has an output consisting of $H_1$ independent GPs, and similarly for ${\mathbf{f}}_2:{\mathbb{R}}^{H_1}\mapsto {\mathbb{R}}^{H_2}$ and so on. The depth and width of DGP are thus denoted by L and $H_{1:L}$, respectively.

3.1. Gaussian Process

In machine learning, the attention is often restricted to the finite set of correlated random variables $\mathbf{f}:=\{f\left({\mathbf{x}}_1\right),\cdots ,f\left({\mathbf{x}}_N\right)\}$ corresponding to the design location $\mathbf{X}={({\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N)}^T$. Denoting $f_i:=f\left({\mathbf{x}}_i\right)$, the above set of random variables is a GP if and only if the following relations,

$$\mathbb{E}\left[f_i\right]=\mu \left({\mathbf{x}}_i\right),\mathbb{E}\left[(f_i-{\mu }_i)(f_j-{\mu }_j)\right]=k({\mathbf{x}}_i,{\mathbf{x}}_j),$$

(1)

are satisfied for all indices $i,j$. For convenience, we can use $f\sim \mathcal{GP}(\mu ,k)$ to denote the above. The mean function $\mu (·):{\mathbb{R}}^d\mapsto \mathbb{R}$ and the covariance function $k(·,·):{\mathbb{R}}^d\times {\mathbb{R}}^d\mapsto \mathbb{R}$ then fully specify the GP. One can proceed to write down the multivariate normal distribution as the pdf

$$p\left(\mathbf{f}\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^N\left|K\right|}}exp[-\frac{1}{2}{\left(\mathbf{f}-\mathbf{m}\right)}^t{K}^{-1}\left(\mathbf{f}-\mathbf{m}\right)].$$

(2)

The covariance matrix K has matrix element $K_{ij}=k({\mathbf{x}}_i,{\mathbf{x}}_j)$, characterizing the correlation between the function values. The covariance function k encodes function properties such as smoothness. The vector $\mathbf{m}:=\mu \left(\mathbf{X}\right)$ represents the mean values at corresponding inputs. Popular covariance functions include the squared exponential (SE) $k({\mathbf{x}}_i,{\mathbf{x}}_j)={\sigma }^{2}exp[-||{\mathbf{x}}_i-{\mathbf{x}}_j{\left|\right|}^2/\left(2{\ell }^2\right)]$ and the family of Matern functions. The signal magnitude $\sigma$ and length scale ℓ are hyper-parameters.

The conditional property of Gaussians allows one to place constraint on the model $p\left(\mathbf{f}\right)$. Given a set of function values $\mathbf{u}=f\left(\mathbf{Z}\right)$, the space of random function f now only includes those passing through these fixed points. Then the conditional pdf $p\left(\mathbf{f}\right|\mathbf{u})$ has the conditional mean and covariance:

$$\begin{array}{cc}\hfill \mathbf{m}& \to \mathbf{m}+K_{\mathbf{xz}}{K}_{\mathbf{zz}}^{-1}\left[\mathbf{u}-\mathbf{m}\right]\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\mathbf{K}}_{\mathbf{X}}& \to {\mathbf{K}}_{\mathbf{X}}-{\mathbf{K}}_{\mathbf{XZ}}{\mathbf{K}}_{\mathbf{Z}}^{-1}{\mathbf{K}}_{\mathbf{ZX}}\hfill \end{array}$$

(4)

where the matrix ${\mathbf{K}}_{\mathbf{XZ}}$ represents the covariance matrix evaluated at $\mathbf{X}$ against $\mathbf{Z}$.

3.2. Deep Gaussian Process

We follow the seminal work in [1] to generalize the notion of GP to the composite functions ${\mathbf{f}}_L\circ {\mathbf{f}}_{L-1}\circ \cdots {\mathbf{f}}_2\circ {\mathbf{f}}_1\left(\mathbf{x}\right)$. In most literature, DGP is defined from a generative point of view. Namely, the joint distribution for the simplest zero-mean DGP with $L=2$ and $H_2=H_1=1$ can be expressed as

$$p(f_2,f_1|\mathbf{X})=p\left(f_2\right|f_1)p\left(f_1\right|\mathbf{X}),$$

(5)

with the conditional defined as $f_2|f_1\sim \mathcal{GP}(0,k(f_1,f_1))$ and $f_1\sim \mathcal{GP}(0,k(\mathbf{X},\mathbf{X})$.

3.3. Marginal Prior, Covariance and Marginal Likelihood

In the above DGP model, the exposed GP for $f_2$ is connected with the data output $\mathbf{y}$, and the intermediate GP for $f_1$ with the data input $\mathbf{X}$. In Bayesian learning, both fs shall be marginalized when computing the evidence. Now we define the marginal prior as

$$p\left(\mathbf{f}\right)=\int d{\mathbf{f}}_{1}p\left({\mathbf{f}}_2\right|{\mathbf{f}}_1)p\left({\mathbf{f}}_1\right|\mathbf{X}).$$

(6)

in which the bold faced ${\mathbf{f}}_1$ representing the set of intermediate function values are marginalized, but the exposed ${\mathbf{f}}_2$ is not. Note that the notation $f\left(\mathbf{x}\right)=f_2\left(f_1\left(\mathbf{x}\right)\right)$ is not ambiguous in a generative view, but may cause some confusion in the marginal view as the label $f_1$ has been integrated out. To avoid confusing with the exposed function $f_2(·)$, we still use $f(·)$ to denote the marginalized composite function unless otherwise stated.

Motivated to write down an objective in terms of marginal likelihood, the moment matching method in [23] was proposed, so one can approximate Equation (6) with a multivariate Gaussian $q\left(\mathbf{f}|\mathbf{X}\right)$ such that the mean and the covariance are matched. In the zero-mean DGP considered in [23], the covariance matching refers to

$${\mathbb{E}}_{\mathbf{f}\sim q}\left[f_i{f}_j\right]={\mathbb{E}}_{{\mathbf{f}}_1}\left[{\mathbb{E}}_{{\mathbf{f}}_2|{\mathbf{f}}_1}\left[f_i{f}_j\right]\right]=\int d{\mathbf{f}}_{2}d{\mathbf{f}}_1{f}_2\left(f_1\left({\mathbf{x}}_i\right)\right)f_2\left(f_1\left({\mathbf{x}}_j\right)\right)p\left({\mathbf{f}}_2\right|{\mathbf{f}}_1)p\left({\mathbf{f}}_1\right|\mathbf{X}).$$

(7)

In the case where the squared exponential kernel is used in both GPs, the approximate marginal prior $q\left(\mathbf{f}|\mathbf{X}\right)=\mathcal{N}(0,K_{\mathrm{eff}})$, with the effective kernel being $k_{\mathrm{eff}}={\sigma }_2^2[1+2\frac{{\sigma }_1^2}{{\ell }_2^2}(1-exp(-|{\mathbf{x}}_i-{\mathbf{x}}_j{|}^2/2{\ell }_1^2{\left)\right)]}^{-\frac{1}{2}}$ [23]. The hyperparameters include the length scale ℓ and signal magnitude $\sigma$ with layer indexed at the subscript.

Consequently, the evidence of the data $\mathbf{X},\mathbf{y}$ associated with the 2-layer DGP can be approximately expressed as

$$p\left(\mathbf{y}|\mathbf{X}\right)\approx \int d\mathbf{f}p\left(\mathbf{y}|\mathbf{f}\right)q\left(\mathbf{f}|\mathbf{X}\right).$$

(8)

Thus, the learning of hyperparameters $\sigma$s and ℓs in the zero-mean DGP model is through the gradient descent on $logp\left(\mathbf{y}|\mathbf{X}\right)$, and the gradient components $\frac{\partial K}{\partial {\ell }_{1,2}}$ and $\frac{\partial K}{\partial {\sigma }_{1,2}}$ are needed in the framework of GPy [48].

4. Model

Following the previous discussion, we shall introduce the model of conditional DGP along with the covariance and marginal prior. The mathematical connection with deep kernel learning and the non-Gaussian aspect of marginal prior will be discussed. The difference between the original DGP and the conditional DGP is that the intermediate GPs in the latter are conditioned on the hyperdata. Learning the hyperdata via the approximate marginal likelihood is, loosely speaking, an empirical Bayesian learning of the feature function in the setting of deep kernel learning.

4.1. Conditional Deep Gaussian Process

In the simple two-layer hierarchy with width $H_1=H_2=1$, the hyperdata $\left\{\mathbf{Z},\mathbf{u}\right\}=\{{\mathbf{z}}_{1:M}\in {\mathbb{R}}^d,u_{1:M}\in \mathbb{R}\}$ are introduced as support for the intermediate GP for $f_1$, and the exposed GP for $f_2$ remains zero-mean and does not condition on any point. Thus, $f_1$ can be viewed as a space of random functions constrained with $f_1\left({\mathbf{z}}_{1:M}\right)=u_{1:M}$, and the Gaussian distribution $p\left(f_1\left({\mathbf{x}}_{1:N}\right)\right|\mathbf{Z},\mathbf{u})$ has its conditional mean and covariance in Equation (3) (with $\mathbf{m}$ on RHS set to zero) and (4), respectively. Following Equation (6), the marginal prior for this conditional DGP can be similarly expressed as

$$p\left(\mathbf{f}\right)=\int d{\mathbf{f}}_{1}p\left({\mathbf{f}}_2\right|{\mathbf{f}}_1)p\left({\mathbf{f}}_1\right|\mathbf{X},\mathbf{Z},\mathbf{u}).$$

(9)

With $f_1$ being conditioned on the hyperdata $\left\{\mathbf{Z},\mathbf{u}\right\}$, one can see that the multivariate Gaussian $p\left(f_1\left({\mathbf{x}}_{1:N}\right)\right|\mathbf{Z},\mathbf{u})$ emits samples in the space of random functions passing through the fixed hyperdata so that Equation (9) is a sum of an infinite number of GPs. Namely,

$$f\sim \sum _{f_1}\mathcal{GP}(0,k_2(f_1\left(\mathbf{X}\right),f_1\left(\mathbf{X}\right))),$$

with $f_1$ under the constraints due to the hyperdata and the smoothness implied in kernel $k_1$. Therefore, f are represented by an ensemble of GPs with same kernel but different feature functions. We shall come back to this point more rigorously in Section 4.2.

Now we shall approximate the intractable distribution in Equation (9) with a multivariate Gaussian $q\left(\mathbf{f}|\mathbf{X},\mathbf{Z},\mathbf{u}\right)$ carrying the matched covariance. The following lemma is useful for the case where the exposed GP for $f_2|f_1$ uses the squared exponential (SE) kernel.

Lemma 1.

(Lemma 3 in [49]) The covariance in $p\left(\mathbf{f}\right)$ (Equation (9)) with the SE kernel $k_2(x,y)={\sigma }_2^{2}exp[-{(x-y)}^2/2{\ell }_2^2]$ in the exposed GP for $f_2|f_1$ can be calculated analytically. With the Gaussian conditional distribution, $p\left({\mathbf{f}}_1\right|\mathbf{X},\mathbf{Z},\mathbf{u})$, supported by the hyperdata, the effective kernel reads

$$k_{\mathrm{eff}}({\mathbf{x}}_i,{\mathbf{x}}_j)=\frac{{\sigma }_2^2}{\sqrt{1+{\delta }_{ij}^2/{\ell }_2^2}}exp\left[-\frac{{(m_i-m_j)}^2}{2({\ell }_2^2+{\delta }_{ij}^2)}\right],$$

(10)

where $m_{i,j}:=m\left({\mathbf{x}}_{i,j}\right)$ and $c_{ij}:=\mathrm{cov}(f_1\left({\mathbf{x}}_i\right),f_1\left({\mathbf{x}}_j\right))$ are the conditional mean and covariance, respectively, at the inputs ${\mathbf{x}}_{i,j}$. The positive parameter ${\delta }_{ij}^2:=c_{ii}+c_{jj}-2c_{ij}$ and the the length scale ${\ell }_2$ dictates how the uncertainty about $f_1$ affects the function composition.

Next, in addition to the hyperparameters such as $\sigma$s and ℓs, the function values $u_{1:M}$ are hyperdata that shall be learned from the objective. With approximating the non-Gaussian marginal prior $p\left(\mathbf{f}|\mathbf{X},\mathbf{Z},\mathbf{u}\right)$ with $q\left(\mathbf{f}|\mathbf{X},\mathbf{Z},\mathbf{u}\right)$, we are able to compute the approximate marginal likelihood as the objective

$$\mathcal{L}=-log\int d\mathbf{f}p\left(\mathbf{y}|\mathbf{f}\right)q\left(\mathbf{f}|\mathbf{X},\mathbf{Z},\mathbf{u}\right).$$

(11)

The learning of all hyperparameter data follows the standard gradient descent used in GPy [48], and the gradient components include the usual ones, such as $\partial K_{\mathrm{eff}}/\partial {\ell }_2$ in exposed layer and those related to the intermediate layer $\partial K_{\mathrm{eff}}/\partial {\ell }_1$ and the hyperdata $\partial K_{\mathrm{eff}}/\partial u_{1:M}$ through chaining with $\partial K_{\mathrm{eff}}/\partial (m_i-m_j)$ and $\partial K_{\mathrm{eff}}/\partial {\delta }_{ij}^2$ via Equations (3) and (4). To exploit the expressive power of neural network during optimization, the hyperdata can be further modeled by a neural network; i.e.,

$$u_{1:M}={\mathrm{nn}}_{\mathbf{w}}\left({\mathbf{z}}_{1:M}\right),$$

(12)

with $\mathbf{w}$ denoting the weight parameters. In such case, the weights $\mathbf{w}$ are learned instead of the hyperdata $u_{1:M}$.

4.2. When Conditional DGP Is Almost a GP

In the limiting case where the probabilistic nature of $f_1$ is negligible, then the conditional DGP becomes a GP with the transformed input; i.e., the distribution $p\left({\mathbf{f}}_1\right|\mathbf{X},\mathbf{Z},\mathbf{u})$ becomes highly concentrated around a certain conditional mean ${\overline{f}}_1\left(\mathbf{x}\right)$. To get insight, we reexamine the covariance in the setting where $f_1$ is almost deterministic. We can reparameterize the random function $f_1$ at two distinct inputs ${\mathbf{x}}_{1,2}$ for the purpose of computing covariance:

$$f_1\left({\mathbf{x}}_{i,j}\right)=m\left({\mathbf{x}}_{i,j}\right)+{ϵ}_{i,j},$$

(13)

where $m\left(\mathbf{x}\right)$ is the conditional mean given the fixed $\mathbf{Z}$ and $\mathbf{u}$. The random character lies in the two correlated random variables, ${({ϵ}_i,{ϵ}_j)}^T\sim \mathcal{N}(0,C)$, corresponding to the weak but correlated signal around zero. Under that assumption, we follow the analysis in [9,38] and prove the following lemma.

Lemma 2.

Consider $p\left(f\right)$ defined in Equation (9), with $f_2|f_1$ being a more general $\mathcal{GP}({\mu }_2,k_2)$ and $f_1$ reparametrized as in Equation (13). The covariance, $\mathrm{cov}(f_2\left(f_1\left({\mathbf{x}}_i\right)\right),f_2\left(f_1\left({\mathbf{x}}_j\right)\right))$, has the following form:

$$[1+\frac{c_{ii}}{2}{\partial }_{m_j}^2+\frac{c_{jj}}{2}{\partial }_{m_i}^2+c_{ij}{\partial }_{m_i{m}_j}^2]k_2(m_i,m_j)+c_{ij}{\mu }_2^{\prime }\left(m_i\right){\mu }_2^{\prime }\left(m_j\right),$$

(14)

where the cs are matrix elements of kernel matrix C associated with the weak random variables ${ϵ}_{i,j}$ in Equation (13). The notation $m_{i,j}:=m\left({\mathbf{x}}_{i,j}\right)$ and prime as derivative is used.

Proof. 

The assumption is that $f_2|f_1\sim \mathcal{GP}({\mu }_2,k_2)$ and that $f_1\left(x\right)$ a weak random function $ϵ\left(x\right)$ overlaying a fixed function $m\left(x\right)$. At any two inputs $x_{i,j}$, we expand the target function f to the second order:

$$f\left(x_{i,j}\right)=f_2\left(f_1\left(x_{i,j}\right)\right)\approx f_2\left(m_{i,j}\right)+{ϵ}_{i,j}{f}_2^{\prime }\left(m_{i,j}\right)+\frac{{ϵ}_{i,j}^2}{2}{f}_2^{″}\left(m_{i,j}\right),$$

(15)

where the shorthanded notation $m_i:=m\left(x_i\right)$ and $ϵ\left(x_i\right):={ϵ}_i$ is used. Note that $({ϵ}_i,{ϵ}_j)$ is bivariate Gaussian with zero mean and covariance matrix $\mathbf{C}$. We use the law of total covariance, $\mathrm{cov}[a,b]=\mathrm{cov}\left(\mathbb{E}\right[a\left|d\right],\mathbb{E}\left[b\right|d\left]\right)+\mathbb{E}\left[\mathrm{cov}\right(a,b\left|d\right)]$ with a, b and d being some random variables. To proceed with the first term, we calculate the conditional mean given the $ϵ$s:

$$\mathbb{E}\left[f\left(x_{i,j}\right)\right|{ϵ}_i,{ϵ}_j]={\mu }_2\left(m_{i,j}\right)+{ϵ}_{i,j}{\mu }^{\prime }\left(m_{i,j}\right).$$

(16)

Then one uses the fact that $f_2|f_1$, $f_2^{\prime }|f_1$ and $f_2^{″}|f_1$ are jointly Gaussian to compute the conditional covariance, which can be expressed in a compact form:

$$\mathrm{cov}\left[f\left(x_i\right)f\left(x_j\right)\right|{ϵ}_i,{ϵ}_j]=\widehat{O}({ϵ}_i,{ϵ}_j)k(m_i,m_j).$$

(17)

The operator $\widehat{O}$ accounts for the fact that $\mathrm{cov}[f_2^{\prime }\left(m_i\right),f_2^{\prime }\left(m_j\right)]={\partial }_{m_i{m}_j}^2{k}_2(m_i,m_j)$ and

$\mathrm{cov}[f_2\left(m_i\right),f_2^{″}\left(m_j\right)]={\partial }_{m_j}^2{k}_2(m_i,m_j)$. Thus, the operator reads

$$\widehat{O}=1+{ϵ}_i{\partial }_{m_i}+{ϵ}_j{\partial }_{m_j}+\frac{{ϵ}_i^2}{2}{\partial }_{m_i}^2+\frac{{ϵ}_j^2}{2}{\partial }_{m_j}^2+{ϵ}_i{ϵ}_j{\partial }_{m_i{m}_j}^2.$$

(18)

Now we are ready to deal with the outer expectation with respect to the $ϵ$s. Note that the covariance $c_{ij}:=\mathbb{E}\left[{ϵ}_i{ϵ}_j\right]=c(x_i,x_j)$ and variance $c_{ii}:=\mathbb{E}\left[{ϵ}_i^2\right]=c(x_i,x_i)$ are matrix elements of $\mathbf{C}$. Consequently, we prove the total covariance in Equation (14). □

Remark 1.

Since the second derivatives ${\partial }_{m_i}^2{k}_2(m_i,m_j)={\partial }_{m_j}^2{k}_2(m_i,m_j)=-{\partial }_{m_i{m}_j}^2{k}_2(m_i,m_j)$ hold for the stationary $k_2$, the above covariance (Equation (14)) with ${\mu }_2=0$ is identical to the effective kernel in Equation (10) in the limit ${\ell }_2^2\gg {\delta }_{ij}^2$, which reads

$$\mathrm{cov}\left(f\left(x_i\right)f\left(x_j\right)\right)\propto [1+\frac{{(m_i-m_j)}^2-{\ell }_2^2}{2{\ell }_2^4}{\delta }_{ij}^2]exp[-\frac{{(m_i-m_j)}^2}{2{\ell }_2^2}].$$

(19)

Such a situation occurs when the inputs $\mathbf{Z}$ in hyperdata are dense enough so that $f_1$ becomes almost deterministic.

Consequently, in the limit when the conditional covariance in ${\delta }^2$ is small compared with the length scale ${\ell }_2^2$, Equation (19) indicates that the effective kernel is the SE kernel with a deterministic input $m\left(\mathbf{x}\right)$, which is equivalent to the deep kernel with SE as the base kernel (see Equation (5) in [16]). On the other hand, when ${\delta }^2$ and ${\ell }_2$ are comparable, the terms within the first bracket in the RHS of Equation (19) are a non-stationary function which may attribute multiple frequencies in the function f. The deep kernel with the spectral mixture kernel (Equation (6) in [16]) as the base is similar to the effective kernel.

4.3. Non-Gaussian Aspect

The statistics of the non-Gaussian marginal prior $p\left(\mathbf{f}|\mathbf{X},\mathbf{Z},\mathbf{u}\right)$ are not solely determined by the moments up to the second order. The fourth moment can be derived in a similar manner in [23] with the help of the theorem in [50]. Relevant discussion of the heavy-tailed character in Bayesian deep neural network can be found in [51,52,53]. See Lemma A1 for the details of computing the general fourth moment in the case where SE kernel is used in $f_2|f_1$ in the conditional 2-layer DGP. Here, we briefly discuss the non-Gaussian aspect, focusing on the variance of covariance, i.e., by comparing ${\mathbb{E}}_p\left[{\left(f\left({\mathbf{x}}_i\right)f\left({\mathbf{x}}_j\right)\right)}^2\right]$ and ${\mathbb{E}}_q\left[{\left(f\left({\mathbf{x}}_i\right)f\left({\mathbf{x}}_j\right)\right)}^2\right]$, with p being the true distribution (Equation (9)) and q being the approximating Gaussian.

In the SE case, one can verify the difference in the fourth order expectation value:

$${\mathbb{E}}_p\left[{\left(f\left({\mathbf{x}}_i\right)f\left({\mathbf{x}}_j\right)\right)}^2\right]-{\mathbb{E}}_q\left[{\left(f\left({\mathbf{x}}_i\right)f\left({\mathbf{x}}_j\right)\right)}^2\right]=\frac{e^{-\frac{{(m_i-m_j)}^2}{1+2{\delta }_{ij}^2}}}{{[1+2{\delta }_{ij}^2]}^{1/2}}-\frac{e^{-\frac{{(m_i-m_j)}^2}{1+{\delta }_{ij}^2}}}{[1+{\delta }_{ij}^2]}\ge 0,$$

where we have used the fact that the inequalities ${(1+2{\delta }^2)}^{-1/2}\ge {(1+{\delta }^2)}^{-1}$ and $exp[-{(1+2{\delta }^2)}^{-1}]\ge exp[-{(1+{\delta }^2)}^{-1}]$ hold. Therefore, the inequality suggests the heavy-tailed statistics of the marginal prior $p(f\left({\mathbf{x}}_i\right),f\left({\mathbf{x}}_j\right))$ over any pair of function values.

5. Results

The works in [54,55] demonstrate that GPs can still have superior expressive power and generalization if the kernels are dedicatedly designed. With the belief that deeper models generalize better than the shallower counterparts [56], DGP models are expected to perform better in fitting and generalization than GP models do if the same kernel is used in both. However, such expectation may not be fully realized, as the approximate inference may lose some power in DGP. For instance, diminishing variance in the posterior over the latent function was reported in [9] regarding the variational inference for DGP [6]. Here, with a demonstration of extrapolating real-world time series data with the conditional DGP, we shall show that the depth, along with optimizing the hyperdata, does enhance the expressive power and the generalization due to the multiple length scale and multiple-frequency character of the effective kernel. In addition, the moment matching method as an approximate inference for conditional DGP does not suffer from the posterior collapse. The simulation codes can be found in the github repository.

5.1. Mauna Loa Data

Figure 1a,b shows fitting and extrapolating the classic carbon dioxide data (yellow marks for training, red for test) with GPs using, respectively, the SE kernel and a mixture of SE, periodic SE and rational quadratic kernels [3].

$$k_{mix}(t,t^{\prime })={\theta }_1^{2}exp[-\frac{{(t-t^{\prime })}^2}{{\theta }_2^2}]+{\theta }_3^{2}exp[-\frac{{sin}^2(t-t^{\prime })}{{\theta }_4^2}]+{\theta }_5^2{[1+\frac{{(t-t^{\prime })}^2}{{\theta }_6^2}]}^{-{\theta }_7}.$$

All the $\theta$s are hyperparameters in the mixture kernel. As a result of the multiple time scales appearing in the data, the vanilla GP fails to capture the short time trend, but the GP with mixture of kernels can still present excellent expressivity and generalization. The log marginal likelihoods (logML) were 144 and 459 for the vanilla GP and kernel mixture GP, respectively. The two-layer zero-mean DGP with SE kernel in both layers performed better than the single-layer counterpart. In Figure 1c, the GP with the SE[SE] effective kernel has excellent fitting with the training data but has extrapolated poorly. The good fitting may have resulted from the fact that the SE[SE] kernel does capture the character of multiple length scales in DGP. The logML for the SE[SE] GP is 338.

Next, we shall see whether improved extrapolation can arise in other deep models or other inference schemes. In Figure 2, the results from DKL and from DGP using the variational inference are shown. Both were implemented in GPFlux [57]. We modified the tutorial code for hybrid GP with three-layer neural network as the code for DKL. The result in Figure 2a does not show good fitting nor good extrapolation, which is to some extent consistent with the simulation of a Bayesian neural network with ReLu activation [32]. As for the DGP using variational inference, the deeper models do not show much improvement compared to the vanilla GP, and the obtained ELBO was 135 for the two-layer DGP (Figure 2b), and it was 127 for three-layer (Figure 2c).

Now we continue to show the performance of our model. In the two-layer model, we have 50 points in hyperdata supporting the intermediate GP. A width-5 tanh neural network is used to represent the hyperdata, i.e., $u_{1:50}={\mathrm{nn}}_{\mathbf{w}}\left(z_{1:50}\right)$. Then, the hyperparameters, including ${\sigma }_{1,2}$, ${\ell }_{1,2}$ and weight parameter $\mathbf{w}$, were learned from gradient descent upon the approximate marginal likelihood. The top panel in Figure 3a displays the prediction and confidence from the posterior over $f_1$, obtained from a GP conditioned on the learned hyperparameters and hyperdata. The logML of the two-layer model was also 338, the same as the SE[SE] GP, and in the bottom panel of Figure 3a one can observe a good fit with the training data. More importantly, the extrapolation shows some high frequency signal in the confidence (shaded region). In comparing ot with Figure 3 of [54], the high-frequency signal only appeared after a periodic kernel was inserted. We attribute the high-frequency signal to the propagation of uncertainty in $f_1$ (top panel) to the exposed layer (see discussion in Section 4.2).

Lastly, the three-layer model using 37 and 23 hyperdata in the $f_1$ and $f_2$ layer, respectively, has its results in Figure 3b. Those hyperdata were parameterized by the same neural network used in the two-layer model. The training had a logML of 253, resulting in a good fit with the data. The extrapolation captured the long term trend in its predictive mean, and the test data were mostly covered in the confidence region. In the latent layers, more expressive patterns overlaying the latent mappings seemed to emerge due to the uncertainty and the depth of the model. The learned ${\sigma }_{1,2,3}\approx (0.49,0.86,2.56)$ and ${\ell }_{1,2,3}\approx (0.014,1.2,0.46)$ show that different layers managed to learn different resolutions.

5.2. Airline Data

The models under consideration can be applied to the airline data too. It can be seen in Figure 4 that the vanilla GP was too simple for the complex time-series data, and the GP with the same kernel composition could both fit and extrapolate well. The logML values were -11.7 and 81.9 for the vanilla GP and kernel mixture GP, respectively. Similarly, the SE[SE] kernel captured the multiple length scale character in the data, resulting in a good fit with logML 20.9, but poor extrapolation.

Here, we display the results using the DKL, variational inference DGP, both two-layer and three-layer, in Figure 5. For the airline data, the DKL with a ReLu neural network as feature extractor panel (Figure 5a) had similar performance to its counterpart on CO${}_2$ data, as did the variational DGP(Figure 5b,c).

Our two-layer model, aided by the probabilistic latent layer supported by 13 hyperdata, showed improved extrapolation along with the high-frequency signal in prediction and confidence. The optimal logML was 28.5, along with the learned ${\sigma }_{1,2}=(3.03,0.73)$, ${\ell }_{1,2}=(0.026,2.19)$ and noise level ${\sigma }_n=0.004$. As shown in Figure 6b, the latent function supported by learned hyperdata shows an increasing trend on top of an oscillating pattern, which led to the periodic extrapolation in the predictive distribution, albeit only the vanilla kernels were used. It is interesting to compare Figure 6b with Figure 6a, as the latter model had 23 hyperdata supporting the latent function, and the vanishing uncertainty learned in the latent function produced an extrapolation that collapsed to zero. The logML in Figure 6a is 7.25 with learned ${\sigma }_{1,2}=(2.18,0.59)$, ${\ell }_{1,2}=(0.3,0.07)$, and noise level ${\sigma }_n=0.02$.

6. Discussion

What did we gain and lose while modifying the original DGP defined in Equation (6) by additionally conditioning the intermediate GPs on the hyperdata? On one hand, when the hyperdata are dense, the conditional DGP is mathematically connected with the deep kernel learning, i.e., a GP with warped input. On the other hand, in the situations when less dense hyperdata are present and the latent GPs are representations of random functions passing through the hyperdata, the conditional DGP can be viewed as an ensemble of deep kernels, and the moment matching method allows us to express it in a closed form. What do we lose in such an approximation? Apparently, the approximate q for the true marginal prior p in Equation (9) cannot account for the heavy-tailed statistics.

In the demonstration, the presence of hyperdata constrains the space of the intermediate functions and moves the mass of the function distribution toward the more probable ones in the process of optimization. Comparing the SE[SE] GP, which represents an approximate version of zero-mean 2-layer DGP, against the conditional DGP model, the constrained space of intermediate functions does not affect the learning significantly, and the generalization is improved. Besides, the uncertainty in the latent layers is not collapsed.

One possible criticism of the present model may result from the empirical Bayes learning of the weight parameters. Although the weight parameters are hyperparameters in both our model and in DKL, it is important to distinguish that the weight parameters in our model parameterize $u_{1:M}$, which supports the intermediate GP, representing an ensemble of latent functions. In DKL, however, the weight parameters fully determine the one latent function, which might lead to overfitting even though marginal likelihood is used as an objective [19]. A possible extension is to consider upgrading the hyperdata to random variables, and the associated mean and variance in $q\left(u_{1:M}\right)$ can also be modeled as neural network functions of $\mathbf{Z}$. The moment matching can then be applied to approximate the marginal prior $\int d{\mathbf{f}}_{1}d\mathbf{u}p\left({\mathbf{f}}_2\right|{\mathbf{f}}_1)p\left({\mathbf{f}}_1\right|\mathbf{X},\mathbf{Z},\mathbf{u})q\left(\mathbf{u}\right)$.

7. Conclusions

Deep Gaussian processes (DGPs), based on nested Gaussian processes (GPs), offer the possibility of expressive inference and calibrated uncertainty, but are limited by intractable marginalization. Approximate inference for DGPs via inducing points and variational inference allows scalable inference, but incurs costs by limiting expressiveness and causing an inability to propagate uncertainty. We introduce effectively deep kernels with optimizable hyperdata supporting latent GPs via a moment-matching approximation. The approach allows joint optimization of hyperdata and GP parameters via maximization of marginal likelihood. We show that the approach avoids mode collapse, connects DGPs and deep kernel learning, effectively propagating uncertainty. Future directions on conditional DGP include that consideration of randomness in the hyperdata and the corresponding inference.