# A Sparse Expansion For Deep Gaussian Processes

Liang Dinga, Rui Tuob and Shahin Shahrampourc
a School of Data Science, Fudan University, Shanghai, China
b Industrial & Systems Engineering, Texas A&M University, College Station, TX
c Mechanical & Industrial Engineering, Northeastern University, Boston, MA

###### Abstract

In this work, we use Deep Gaussian Processes (DGPs) as statistical surrogates
for stochastic processes with complex distributions. Conventional inferential
methods for DGP models can suffer from high computational complexity as they
require large-scale operations with kernel matrices for training and
inference. In this work, we propose an efficient scheme for accurate inference
and efficient training based on a range of Gaussian Processes, called the
_Tensor Markov Gaussian Processes_ (TMGP). We construct an induced
approximation of TMGP referred to as the _hierarchical expansion_. Next, we
develop a deep TMGP (DTMGP) model as the composition of multiple hierarchical
expansion of TMGPs. The proposed DTMGP model has the following properties: (1)
the outputs of each activation function are deterministic while the weights
are chosen independently from standard Gaussian distribution; (2) in training
or prediction, only $\mathcal{O}({\text{polylog}}(M))$ (out of $M$) activation
functions have non-zero outputs, which significantly boosts the computational
efficiency. Our numerical experiments on synthetic models and real datasets
show the superior computational efficiency of DTMGP over existing DGP models.

Keywords: deep Gaussian processes; Markov Gaussian processes; inducing
variables; sparse expansion.

##  1 Introduction

This work is partially motivated by surrogate modeling of stochastic computer
simulations (Ankenman et al., 2010; Plumlee and Tuo, 2014). The ultimate goal
in this field is to approximate an underlying stochastic process, defined by a
complex stochastic computer simulation code, using a surrogate stochastic
process that is computationally efficient. Among existing statistical
surrogates, Deep Gaussian Processes (DGP) proposed by Damianou and Lawrence
(2013) has recently become a popular choice owing to its exceptional
performance and flexibility in fitting the complex distributions from the
stochastic processes in real applications (Radaideh and Kozlowski, 2020; Sauer
et al., 2022).

DGPs are multi-layer compositions of multi-variate Gaussian Processes (GP).
The deep probabilistic structures of DGPs allows a Bayesian formulation to
model complex stochastic systems, such as computer vision, natural language
processing, etc.. However, the cost of DPGs’ expressiveness and flexibility is
the extreme difficulties in training and inference, because the propagation of
randomness throughout the layers of DGPs is nonlinear and correlated among
layers. Existing approaches for overcoming these difficulties can be separated
into two classes – one class focuses on approximating the Bayesian
formulations of DGPs, such as variational inference (VI) (Blei et al., 2017),
expectation propagation (EP) (Bui et al., 2016), or Vecchia approximation
(Sauer et al., 2022a), while the other class considers simplified
representations of DGPs’ structures, such as inducing-variable approximation
(Hensman and Lawrence, 2014) and random Fourier feature expansion (Cutajar et
al., 2017).

Simplifications based on inducing variables and random Fourier features
(Rahimi and Recht, 2008) partially address the computational issue of DGPs.
However, they are still hard to implement. This is because these
simplifications of DGPs are deep Neural Networks (DNN) with kernel based
activations, correlated random weights and bias parameters, which is
equivalent to Bayesian Neural Networks (BNN) (MacKay, 1992; Neal, 1996) with
dense and highly correlated structures. Training and inference of these
simplifications are at least as hard as existing BNN models.

In this work, we focus on introducing an accurate and efficient simplification
of DGPs. We propose a sparse reduced-rank approximation, referred as
_Hierarchical Expansion_ , for DGPs. Our expansion is specialized for DPGs
which are the compositions of multi-variate tensor Markov GPs (TMGP) (Ding and
Zhang, 2020), so we call this expanded DGP deep TMGP (DTMGP). One of our main
contributions lies in constructing a sparse representation of DGPs. There are
only a poly-logarithmic number of activations in our model with non-zero
outputs each time we run the model. Because of this sparse property, training
of DTMGPs are much faster and easier compared with existing DGP models,
avoiding commonly seen numerical issues in training deep models. Another
contribution is to show that hierarchical expansion is also highly accurate.
Unlike many other GPs/DGPs sparse approximations, the difference between
sample paths generated by any DGP and those generated by its hierarchical
expansion is small. This property allows DTMGPs to model complex stochastic
systems with satisfactory performance as shown in our numerical experiments.

We use VI to train DTMGP on data generated from simulations and real data and
compare it with other existing DGP models. We conduct simulation studies by
running experiments on an artificial random field and simulator of a
stochastic activity network (Pasupathy and Henderson, 2011). We also use the
RC-49 data set, which is publicly available on the web
https://paperswithcode.com/dataset/rc-49, to demonstrate the performance of
the proposed approach. In these numerical studies, we find that the training
process of DTMGP is more stable, and the training loss converges faster
compared with its competitors. Besides, we compare the similarities between
instances from the underlying systems and instances generated from competing
statistical surrogate models. The results show that DTMGP outperforms the
alternative DGPs in all experiments.

The remainder of this article is organized as follows. We will review the
related literature in Section 2. We introduce general DGP models and an
approximation of DGPs called induced approximation in Section 3. The
methodology and detailed implementation of DTMGP are introduced in Section 4.
Simulation studies are given in Section 5 and experiments on real data are
given in 6 . Concluding remarks are made in Section 7. The Appendix includes
the required mathematical tools and technical proofs.

##  2 Literature Review

The state-of-the-art deep learning techniques have brought probabilistic
modeling with deep neural network structure in popularity. Based on the
concept of deep belief network (DBN) (Hinton et al., 2006), Damianou and
Lawrence (2013) generalized the Restricted Boltzmann Machine (Hinton, 2010),
which is a DBN with binary output, to deep Gaussian Processes based on
Gaussian process mapping. DGPs have been a popular tool in several
applications. Their prowess has been demonstrated on many classification tasks
(Damianou and Lawrence, 2013; Fei et al., 2018; Yang and Klabjan, 2021).
Compared with traditional DNNs, the flexibility in uncertainty quantification
of DGPs makes them ideal candidate for surrogate modeling (Radaideh and
Kozlowski, 2020; Sauer et al., 2022). DGPs are commonly used statistical
surrogates in many applications such as Bayesian optimization (Hebbal et al.,
2021), calibration (Marmin and Filippone, 2022), multi-fidelity analysis (Ko
and Kim, 2021), healthcare (Li et al., 2021), and etc.

However, training and inference of DGPs are difficult. Current attempts can be
separated into two classes. One focuses on designing more efficient and
accurate training algorithm while the other one on constructing DGP
architectures with sparsity. Efficient training and inference algorithms
include expectation propagation (Bui et al., 2016), doubly stochastic
variational inference (Salimbeni and Deisenroth, 2017), stochastic gradient
Hamiltonian Monte Carlo (Havasi et al., 2018), elliptical slice sampling
(Sauer et al., 2022), Vecchia approximation (Vecchia, 1988; Sauer et al.,
2022a; Katzfuss and Guinness, 2021), and etc. Reformulation of DGPs to
simplified models is an alternative approach. DGPs are reformulated as an
variational model in Tran et al. (2016). Low-rank approximations for GPs in
Banerjee et al. (2008); Cressie and Johannesson (2008) are also extended to
DGPs. Hensman and Lawrence (2014); Dai et al. (2016) extend the inducing-
variable method to DGPs. Cutajar et al. (2017) uses random Fourier feature
expansion to show that DGPs are equivalent to infinitely wide BNN with
corresponding activations.

In this work, we consider DGPs with Markov structure and, based on the Markov
structure, we design an accurate and efficient sparse expansion for DGPs.
Sidén and Lindsten (2020) proposes compositions of _discrete_ Gaussian Markov
random field on graph and called deep graphical models of this structure _Deep
Gaussian Markov Random Fields_ (DGMRF). We must point out that our DTMGPs is
essentially different from DGMRF because DTMGPs are _not_ deep graphical
models. More specifically, DGMRF treats every activation as a random variable
and imposes a graphical Markov structure on these activations, i.e., any two
non-adjacent randomly distributed activations are conditionally independent
given all other activations. On the other hand, every activation in DTMGP is
one of the orthonormal basis functions of a _continuous_ Gaussian Markov
random field and these orthonormal basis functions constitute a so-called
_hierarchical expansion_ of the field.

##  3  General DGP Models

###  3.1 Deep Gaussian Processes

A GP $\mathcal{G}$ is a random function characterized by its mean function
$\mu$ and covariance function $k$. To be more specific, given any input
$\mathbb{x}$, the output $\mathcal{G}(\mathbb{x})$ has a Gaussian distribution
with mean $\mu(\mathbb{x})$ and variance $k(\mathbb{x},\mathbb{x})$, and given
any pair of inputs $\mathbb{x}$ and $\mathbb{x}^{\prime}$, the covariance
between $\mathcal{G}(\mathbb{x})$ and $\mathcal{G}(\mathbb{x}^{\prime})$ is
$k(\mathbb{x},\mathbb{x}^{\prime})$:

$\displaystyle\mathcal{G}(\mathbb{x})\sim\mathcal{N}\big{(}\mu(\mathbb{x}),k(\mathbb{x},\mathbb{x})\big{)},$
(1)
$\displaystyle\textbf{Cov}\big{(}\mathcal{G}(\mathbb{x}),\mathcal{G}(\mathbb{x}^{\prime})\big{)}=k(\mathbb{x},\mathbb{x}^{\prime}).$

A $W$-variate GP is simply a $W$-vector of GPs
$[\mathcal{G}_{1},\mathcal{G}_{2},\cdots,\mathcal{G}_{W}]$. Without loss of
generality, we assume that any multivariate GP in the following content has
mean function $\mu=0$ and independently distributed entries. A $H$-layer DGP
$\mathbb{f}^{(H)}$ is then the composition of $H$ multi-variate GPs:

$\mathbb{f}^{(H)}(\mathbb{x}^{*})=\mathcal{G}^{(H)}\circ\cdots\circ\mathcal{G}^{(2)}\circ\mathcal{G}^{(1)}(\mathbb{x}^{*}),$
(2)

where $f\circ g$ denotes the composition function $f(g(\cdot))$,
$\mathcal{G}^{(h)}$ denotes a $W^{(h)}$-variate GP
$[\mathcal{G}_{1}^{(h)},\mathcal{G}_{2}^{(h)},\cdots,\mathcal{G}^{(h)}_{W^{(h)}}]$,
and $W^{(h)}$ is called the width of layer $h$. The formulation (2) of DGP can
be viewed as a DNN with random activations $\\{\mathcal{G}_{i}^{(h)}\\}$.

While DGP yields a non-parametric and Bayesian formulation of DNNs,
computations for inference/prediction of DGPs are cumbersome. For example,
given a set of input-output pairs
$(\mathbb{X},\mathbb{Y})\in\mathbb{R}^{n\times d}\times\mathbb{R}^{n}$, the
conditional probability distribution $\mathbb{P}(\mathbb{Y}\big{|}\mathbb{X})$
induced by (2) is

$\mathbb{P}\bigl{(}\mathbb{Y}\big{|}\mathbb{X}\bigr{)}=\int\mathbb{P}\bigl{(}\mathbb{Y}\big{|}\mathbb{f}^{(H)}\bigr{)}\mathbb{P}\bigl{(}\mathbb{f}^{(H-1)}\big{|}\mathbb{f}^{(H-2)}\bigr{)}\cdots\mathbb{P}\bigl{(}\mathbb{f}^{(1)}\big{|}\mathbb{X}\bigr{)}d{\mathbb{f}^{(1)}}\cdots
d{\mathbb{f}^{(H)}},$ (3)

where $\mathbb{f}^{(h)}$ denotes the output vector from layer $h$:

$\mathbb{f}^{(h)}=\mathcal{G}^{(h)}\circ\mathcal{G}^{(h-1)}\circ\cdots\circ\mathcal{G}^{(1)},$

with $h=1,\cdots,H$. Evidently, the integral (3) is intractable, leading to
the intractability of inference and predictions of DGP $\mathbb{f}^{(H)}$
conditioned on any observed data $\bigl{(}\mathbb{X},\mathbb{Y}\bigr{)}$.

###  3.2 Induced Approximation of DGPs

The inference of GP, which acts as a single activation in DGP, is also time
and space consuming. A challenge for GPs lies in their computational
complexity and storage for the inverse of covariance matrix
$k(\mathbb{X},\mathbb{X})$ which are $\mathcal{O}(n^{3})$ and
$\mathcal{O}(n^{2})$, respectively, when the covariance matrix is of size $n$.
To alleviate this computational problem, current methods, such as inducing-
variable approximation in Hensman and Lawrence (2014); Dai et al. (2016) and
random Fourier feature approximation in Cutajar et al. (2017), focus on
approximating the covariance matrix $k(\mathbb{X},\mathbb{X})$ by some low
rank matrices.

Instead of approximating the covariance matrix, we approximate the GP
activations in a DGP directly. In the following content, we will construct a
reduced-rank approximation called induced approximation. Induced approximation
can bring a more flexible and clear deep neural network representation of
DGPs.

Without loss of generality, let $\mathcal{G}$ defined in (1) be a GP with
constant mean $\mu$ and covariance function $k$. Inspired by the kriging
method (Matheron, 1963; Sacks et al., 1989), we can approximate GP
$\mathcal{G}$ defined in (1) by the following finite-rank approximation:

$\hat{\mathcal{G}}(\mathbb{x}^{*})\coloneqq\mu+k(\mathbb{x}^{*},\mathbb{U})\big{[}k(\mathbb{U},\mathbb{U})\big{]}^{-1}\mathcal{G}(\mathbb{U})$
(4)

where we call $\mathbb{U}=\\{\mathbb{u}_{i}\in\mathbb{R}^{d}\\}_{i=1}^{m}$ the
inducing points and $\mathcal{G}(\mathbb{U})$ are the values of $\mathcal{G}$
on $\mathbb{U}$. According to theory in (Ritter, 2000; Wang et al., 2020; Ding
and Zhang, 2020), the required number of inducing points for achieving the
optimal statistical error is much lower than the sample size, provided that
the inducing points $\mathbb{U}$ are well-chosen. Hence, (4) can be computed
in a highly efficient way.

Figure 1: Induced approximation $\hat{\mathcal{G}}$ can be represented as a
one-layer neural network with correlated Gaussian distributed weights
$\mathbb{w}=\big{[}k(\mathbb{U},\mathbb{U})\big{]}^{-1}\mathcal{G}(\mathbb{U})$
and bias $\mu$.

Finite-rank approximation (4) can be represented as a one-layer neural network
with correlated Gaussian distributed weights
$\mathbb{w}=\big{[}k(\mathbb{U},\mathbb{U})\big{]}^{-1}\mathcal{G}(\mathbb{U})$
and bias $\mu$ as shown in Figure 1. However, the correlated Gaussian
distributed weights still make training and inference difficult. In order to
write the finite-rank approximation (4) in the form of a neural network with
independently distributed weights, we can simply apply Cholesky decomposition
on kernel matrix $k(\mathbb{U},\mathbb{U})$, which leads to the following
equation:

$\hat{\mathcal{G}}(\mathbb{x}^{*})=\mu+k(\mathbb{x}^{*},\mathbb{U})R^{-1}_{\mathbb{U}}\mathbb{Z}\coloneqq\mu+\bm{\phi}^{T}(\mathbb{x}^{*})\mathbb{Z},$
(5)

where $R_{\mathbb{U}}$ is the Cholesky decomposition of matrix
$k(\mathbb{U},\mathbb{U})$ and
$\mathbb{Z}=[R_{\mathbb{U}}^{T}]^{-1}\mathcal{G}(\mathbb{U})\in\mathbb{R}^{m}$
are i.i.d. distributed standard Gaussian random variables. We call
approximation (5) induced approximation and Figure 2 shows that induced
approximation can be represented as a two-layer neural network where functions
$\bm{\phi}$ act as an extra hidden layer.

Figure 2:  A two-layer neural network representation of $\hat{\mathcal{G}}$
with deterministic weights $\mathbb{W}=\big{[}R_{\mathbb{U}}^{T}\big{]}^{-1}$,
i.i.d. standard Gaussian weights $\mathbb{Z}$, and bias $\mu$.

Because the approximation (5) consists of only $m$ Gaussians, the time and
space complexity for its inference are then reduced to $\mathcal{O}(m^{2}n)$
and $\mathcal{O}(mn)$, respectively.

The induced approximation of a DGP is to approximate each random activation by
corresponding induced approximation. For example, suppose that in a DGP
$\mathbb{f}^{(H)}$, all activations in the same layer are characterized by the
same covariance function $k^{(h)}$. Then, the induced approximation
$\tilde{\mathbb{f}}^{(H)}$ of $\mathbb{f}^{(H)}$ is

$\displaystyle\tilde{\mathbb{f}}^{(H)}(\mathbb{x}^{*})=\big{[}\mathbb{Z}^{(H)}\bm{\phi}^{(H)}(\cdot)+\bm{\mu}^{(H)}\big{]}\circ\cdots\circ\big{[}\mathbb{Z}^{(1)}\bm{\phi}^{(1)}(\mathbb{x}^{*})+\bm{\mu}^{(1)}\big{]},$
(6)
$\displaystyle[\bm{\phi}^{(h)}(\cdot)]^{T}=k^{(h)}(\cdot,\mathbb{U}^{(h)})R^{-1}_{{h}}\quad\text{for}\
h=1,\cdots,H$

where, for $h=1,\cdots,H$, $\mathbb{Z}^{(h)}\in\mathbb{R}^{W^{(h)}\times
m^{(h)}}$ is a matrix with i.i.d. standard Gaussian entries, $U^{(h)}$ are the
total $m^{(h)}$ inducing variables for kernel $k^{(h)}$,
$R_{h}\in\mathbb{R}^{m^{(h)}\times m^{(h)}}$ is the Cholesky decomposition of
the covariance matrix $k^{(h)}(\mathbb{U}^{(h)},\mathbb{U}^{(h)})$, and
$\bm{\mu}^{(h)}\in\mathbb{R}^{W^{(h)}}$ is the mean vector and can be treated
as the bias of layer $h$.

In the next section, we will show that if kernel $k$ is in a class of kernel
function called tensor Markov kernel and inducing points $\mathbb{U}$ are a
specific design called sparse grid, then induced approximations (5) and (6)
can be further written as a sparse approximation.

##  4 Methodology of DTMGP

In this section, we first introduce the concept of TMGPs. We then introduce
hierarchical expansion of TMGPs and how it leads to mutually orthogonal
feature functions with hierarchical supports. Based on hierarchical expansion,
we can introduce the implementation and training of DTMGPs.

###  4.1 Tensor Markov GPs

Hierarchical expansion is applied to a class of GPs called tensor Markov GP.
In one dimension, a Markov GP is characterized by the following Lemma in
Marcus and Rosen (2006):

######  Lemma 1 (Marcus and Rosen (2006) Lemma 5.1.8).

Let $I\subseteq\mathbb{R}$. Let $\mathcal{G}$ be a zero mean GP defined on $I$
with continuous positive definite kernel $k$. Then, $\mathcal{G}$ is a Markov
GP and $k$ is Markov kernel if and only if there exist positive functions $p$
and $q$ on $I$ with $p/q$ strictly increasing such that

$k(x,x^{\prime})=p(\min(x,x^{\prime}))q(\max(x,x^{\prime})),\quad
x,x^{\prime}\in I.$ (7)

For any Markov GP $\mathcal{G}$ we have that, conditioned on its distribution
at a point $x$, its distributions at any point $x_{l}<x$ and any point
$x_{u}>x$ are independent:
$\mathcal{G}(x_{l})\perp\mathcal{G}(x_{u})\big{|}\mathcal{G}(x)$, whence
Markov GP. Markov GPs have two advantages in computations. Firstly, given any
ordered inducing points $\\{u_{1}<\cdots<u_{n}\\}$, the distribution of
$\mathcal{G}(x)$ at any $x$, where $u_{i}<x<u_{i+1}$, depends only on its left
and right inducing points $\mathcal{G}(u_{i})$ and $\mathcal{G}(u_{i+1})$.
Secondly, for any induced approximation, accuracy of the approximation at $x$
only depends on the distances between $x$ and its left and right neighbors,
respectively. Therefore, evenly distributed inducing variables can achieve an
accurate approximation that is also cheap to compute.

The following purely additive and purely tensor structure of one-dimensional
Markov kernel can both extend Markov GPs to multi-dimension:

$k(\mathbb{x},\mathbb{x}^{\prime})=\sum_{j=1}^{d}k_{j}(x_{j},x_{j}^{\prime}),\quad
k(\mathbb{x},\mathbb{x}^{\prime})=\prod_{j=1}^{d}k_{j}(x_{j},x_{j}^{\prime}).$
(8)

We generalize the additive and tensor form in (8) to the following GP, which
has additive and tensor structure simultaneously:

######  Definition 1.

A GP is called TMGP if and only if it is a zero mean GP with covariance
function of the form

$k(\mathbb{x},\mathbb{x}^{\prime})=\sum_{l=1}^{s}\prod_{j\in\bm{U}_{l}}k_{l,j}(x_{j},x_{j}^{\prime}),$
(9)

where $\bm{U}_{l}\subseteq\\{1,2,\cdots,d\\}$ and $k_{l,j}$ is 1-D Markov
kernel for any $l,j$. A kernel of the form (10) is called tensor Markov kernel
(TMK).

A TMGP $\mathcal{G}$ with kernel (9) can be factorized as sum of $s$ mutually
independent GPs: $\mathcal{G}=\sum_{l=1}^{s}\mathcal{G}_{l}$ where each
$\mathcal{G}_{l}$ is with kernel function
$\prod_{j\in\bm{U}_{l}}^{d}k_{l,j}(x_{j},x_{j}^{\prime})$. To achieve an
accurate approximation of TMGP using the smallest number of inducing
variables, we can only focus on TMGP with a purely tensor structure. By
combining mutually independent approximations, we can arrive at a solution.
Therefore, in the following content, we only study approximation for TMGP with
purely tensor structure:

$k(\mathbb{x},\mathbb{x}^{\prime})=\prod_{j=1}^{d}k_{j}(x_{j},x_{j}^{\prime}),$
(10)

and the general TMGP with kernel (9) is a straightforward extension that
involves combining approximations of GPs with purely tensor kernels.

Commonly used TMKs with purely tensor structure include Laplace kernel
$k(\mathbb{x},\mathbb{x}^{\prime})=\exp(\sum_{j=1}^{d}\theta_{j}|x_{j}-x_{j}^{\prime}|)$
and Brownian sheet kernel
$k(\mathbb{x},\mathbb{x}^{\prime})=\prod_{j=1}^{d}(1+\theta_{j}\min\\{x_{j},x_{j}^{\prime}\\})$.
The challenge in extending the inducing approximation of 1-D Markov GP to
multiple dimensions lies in properly defining the ”left” and ”right” neighbors
in a multidimensional space, similar to how it’s done on the real line. To
address this, we introduce the hierarchical expansion method, which is
explained in the following section.

###  4.2 Hierarchical Expansion

To construct the hierarchical expansion, we let inducing variables
$\mathbb{U}$ be an experimental design $\mathbb{X}^{\rm SG}_{l}$ called
level-l _Sparse Grid_ (SG) (Bungartz and Griebel, 2004; Plumlee, 2014), where
level $l$ determines the number of inducing variables. Furthermore, the SG
inducing variables must be sorted in a specific order. Detailed introduction
of SGs and the required order are provided in Appendix A, and MATLAB codes for
generating SGs satisfying the requirement can be found in the _Sparse Grid
Designs_ package (Plumlee, 2021). Examples of two-dimensional SGs are shown in
Figure 3. From the examples, we can see that the incremental points from the
next level of a SG exhibit a hierarchical structure – higher level SG consists
of local SGs of smaller levels.

Figure 3: Two-dimensional level-2 SG (left), level-3 SG (middle) and level-4
SG (right) on cube $(0,1)^{2}$. The incremental points from the second level,
third level and fourth level are labeled by *, o and +, respectively.

Hierarchical expansion of a TMGP $\mathcal{G}$ with inducing points
$\mathbb{X}^{\rm SG}_{l}$ is then simply:

$\mathcal{G}(\mathbb{x}^{*})\approx k(\mathbb{x}^{*},\mathbb{X}^{\rm
SG}_{l})R_{l}^{-1}\mathbb{Z}\coloneqq\bm{\phi}^{T}(\mathbb{x}^{*})\mathbb{Z},$
(11)

where $k(\mathbb{x}^{*},\mathbb{X}^{\rm
SG}_{l})=[k(\mathbb{x}^{*},\mathbb{x}_{1}),\cdots,k(\mathbb{x}^{*},\mathbb{x}_{m_{l}})]$
is the covariance vector, $k(\mathbb{X}^{\rm SG}_{l},\mathbb{X}^{\rm
SG}_{l})=[k(\mathbb{x}_{i},\mathbb{x}_{j})]_{i,j\leq m_{l}}$ is the covariance
matrix, $R_{l}$ is the Cholesky decomposition of $k(\mathbb{X}^{\rm
SG}_{l},\mathbb{X}^{\rm SG}_{l})$, $\mathbb{Z}=[Z_{1},\cdots,Z_{m_{l}}]$ are
i.i.d. standard Gaussian random variables and
$\bm{\phi}=[\phi_{i}]_{i=1}^{m_{l}}$ are called hierarchical features.
Hierarchical expansion (11) yields a sparse approximation of TMGP that is easy
to compute as shown in the following theorems.

######  Theorem 1.

The number of non-zero entries on $R_{l}^{-1}$ is $\mathcal{O}(m_{l})$ and
$R^{-1}_{l}$ can be computed in $\mathcal{O}(m_{l})$ operations.

######  Theorem 2.

Given any input point $\mathbb{x}^{*}$, the number of non-zero entries on the
vector of hierarchical features $\bm{\phi}(\mathbb{x}^{*})$ is
$\mathcal{O}([\log m_{l}]^{2d-1})$.

The algorithm for computing $R^{-1}_{l}$ in $\mathcal{O}(m_{l})$ operations is
given in Appendix B, and proofs for Theorems 1 and 2 are provided in
supplementary material. Since we can compute $R^{-1}_{l}$ in
$\mathcal{O}(m_{l})$ operations and the numbers of non-zero entries on
$R^{-1}_{l}$ is $\mathcal{O}(m_{l})$, hierarchical expansion can be computed
in $\mathcal{O}(m_{l})$ operations. Moreover, if matrix $R^{-1}_{l}$ is given,
because the numbers of non-zero entries on $\bm{\phi}(\mathbb{x}^{*})$ is only
$\mathcal{O}([\log m_{l}]^{2d-1})$, the computational time of the hierarchical
expansion can be further reduced to $\mathcal{O}([\log m_{l}]^{2d-1})$. The
sparsity of hierarchical feature vector $\bm{\phi}$ relies on the fact that
the supports of $\\{\phi_{i}\\}_{i=1}^{m_{l}}$ are either nested or disjoint
and, hence, form a hierarchical structure as shown in Figure 4.

Figure 4:  Support of hierarchical feature $\phi_{1}$ corresponding to
$\mathbb{x}_{1}$ in level 1 (left); supports of $\\{\phi_{i}\\}_{i=2}^{5}$
corresponding to $\\{\mathbb{x}_{i}\\}_{i=2}^{5}$ (middle figures); supports
of $\\{\phi_{i}\\}_{i=6}^{9}$ corresponding to
$\\{\mathbb{x}_{i}\\}_{i=6}^{9}$ in level 3 (right).

$\mathbb{X}^{\rm SG}_{l}$ is also regarded as the index set of $\bm{\phi}$,
because instead of labeling each hierarchical feature by number $i$, we can
also label it by the corresponding inducing variable
$\mathbb{x}_{i}\in\mathbb{X}^{\rm SG}_{l}$. As shown in Figure 5, for any $i$,
$\mathbb{x}_{i}$ is the only point at which $\phi_{i}$ is not differentiable
along all dimensions. Hence, we call $\mathbb{x}_{i}$ the center of
$\phi_{i}$. From this perspective, we can claim that, given any
$\mathbb{x}^{*}$, the distribution of the approximation at $\mathbb{x}^{*}$
only depends on $\mathcal{O}([\log m_{l}]^{2d-1})$ inducing variables. As a
result, the hierarchical expansion generalizes the one-dimensional Markov
property while still remaining a high accuracy approximate. We discuss the
approximation error of hierarchical expansion in supplementary material.

Figure 5:  Hierarchical features $\phi_{1}$ (left) and $\phi_{2},\phi_{3}$
(right) generated by the tensor Brownian Bridge kernel
$k(\mathbb{x},\mathbb{x}^{\prime})=\prod_{j=1}^{2}\min\\{x_{j},x_{j}^{\prime}\\}(1-\max\\{x_{j},x_{j}^{\prime}\\})$.

The following numerical example illustrates the sparsity of hierarchical
expansion:

######  Example 1.

Let two-dimensional TMK
$k(\mathbb{x},\mathbb{x}^{\prime})=\exp\\{-\sum_{j=1}^{2}|x_{j}-x_{j}^{\prime}|\\}$.
Let

$\mathbb{X}^{\rm
SG}_{2}=\left[\left(\frac{1}{2},\frac{1}{2}\right),\left(\frac{1}{2},\frac{1}{4}\right),\left(\frac{1}{2},\frac{3}{4}\right),\left(\frac{1}{4},\frac{1}{2}\right),\left(\frac{3}{4},\frac{1}{2}\right)\right],$

be the sorted level-2 SG. Then, we numerically compute the inverse Cholesky
decomposition $R_{2}^{-1}$ and we get

$R^{-1}_{2}={\tiny\begin{bmatrix}1.0000&-1.2416&-1.2416&-1.2416&-1.2416\\\
0&1.5942&0.0000&0.0000&0.0000\\\ 0&0&1.5942&0.0000&0.0000\\\
0&0&0&1.5942&0.0000\\\ 0&0&0&0&1.5942\end{bmatrix}}.$

Let $\mathbb{x}^{*}_{1}=(0.8147,0.9058)$, and
$\mathbb{x}^{*}_{2}=(0.2785,0.5469)$ be two randomly chosen points. We
numerically compute the hierarchical features at $\mathbb{x}^{*}_{1}$ and
$\mathbb{x}^{*}_{2}$ and we get

$\displaystyle k(\mathbb{x}^{*}_{1},\mathbb{X}^{\rm
SG}_{2})R^{-1}_{2}=\big{[}0.4865,0,0.3918,0.0000,0.3918\big{]},$
$\displaystyle k(\mathbb{x}^{*}_{2},\mathbb{X}^{\rm
SG}_{2})R^{-1}_{2}=\big{[}0.7646,0.5291,0,0,0.0933\big{]}.$

As we can see, matrix $R^{-1}_{2}$ is sparse and each of
$k(\mathbb{x}^{*}_{1},\mathbb{X}^{\rm SG}_{2})R^{-1}_{2}$ and
$k(\mathbb{x}^{*}_{2},\mathbb{X}^{\rm SG}_{2})R^{-1}_{2}$ has two zero
elements. The sparsity is even more visible as the number of inducing points
increases.

###  4.3 Deep Tensor Markov GPs

A $H$-layer DTMGP $\mathcal{T}^{(H)}$ is defined as the composition of $H$
hierarchical expansions of multi-variate TMGP:

$\mathcal{T}^{(H)}(\mathbb{x}^{*})=[\mathbb{Z}^{(H)}\bm{\phi}^{(H)}(\cdot)+\bm{\mu}^{(H)}]\circ\cdots\circ[\mathbb{Z}^{(1)}\bm{\phi}^{(1)}(\mathbb{x}^{*})+\bm{\mu}^{(1)}],$
(12)

where, for $h=1,\cdots,H$, $\mathbb{Z}^{(h)}$ is a $W^{(h)}$-by-$m^{(h)}$
matrix with i.i.d. standard Gaussian entries, $W^{(h)}$ is the dimension of
the output (width of the neural network) of layer $h$, $m^{(h)}$ is the number
of hierarchical features for approximating a $W^{(h)}$-variate TMGP,
$\bm{\phi}^{(h)}$ consists of the $m^{(h)}$ hierarchical features for
approximating the TMGP, acting as the activation functions in layer $h$, and
$\bm{\mu}^{(h)}\in\mathbb{R}^{W^{(h)}}$ is the parameterized mean for the TMGP
in layer $h$.

Let us now focus on layer $h$ of $\mathcal{T}^{(H)}$ to compare DTMGP with
DNN. Let
$\mathcal{T}^{(h)}=(\mathcal{T}^{(h)}_{1},\cdots,\mathcal{T}^{(h)}_{W^{(h)}})$
denote the $W^{(h)}$-dimensional output of layer $h$. We then have the
following structure:

$\mathcal{T}^{(h-1)}\rightarrow\underbrace{\big{(}\phi^{(h)}_{1}(\mathcal{T}^{(h-1)}),\cdots,\phi^{(h)}_{m^{(h)}}(\mathcal{T}^{(h-1)})\big{)}}_{\text{\\#
non-zero entries:}\ \mathcal{O}\big{(}[\log
m^{(h)}]^{2W^{(h-1)}-1}\big{)}}\rightarrow\bigg{[}\mathcal{T}_{j}^{(h)}=\mu^{(h)}_{j}+\sum_{i=1}^{m^{(h)}}Z_{j,i}\phi^{(h)}_{i}(\mathcal{T}^{(h-1)})\bigg{]}_{j=1}^{W^{(h)}},$

where $\phi^{(h)}_{i}(\mathcal{T}^{(h-1)})$ is the $i^{\rm th}$ entry of
vector $k^{(h)}(\mathcal{T}^{(h-1)},\mathbb{X}^{\rm
SG}_{l^{(h)}})[R^{(h)}]^{-1}$, $k^{(h)}$ is the TMK for layer $h$,
$\mathbb{X}^{\rm SG}_{l^{(h)}}$ is a level-$l^{(h)}$ SG of size $m^{(h)}$ for
induced approximation, and $R^{(h)}$ is the Cholesky decomposition of
covariance matrix $k^{(h)}(\mathbb{X}^{\rm SG}_{l^{(h)}},\mathbb{X}^{\rm
SG}_{l^{(h)}})$. This structure is also illustrated in Figure 6.

Figure 6: Only $\mathcal{O}\big{(}[\log m^{(h)}]^{2W^{(h-1)}-1}\big{)}$
activations are non-zero in layer $h$ of DTMGP $\mathcal{T}^{(H)}$.

If we treat DTMGP as a DNN, then, for $h=1,\cdots,H$, layer $h$ is equipped
with $m^{(h)}$ activations
$\\{k^{(h)}(\cdot,\mathbb{x}):\mathbb{x}\in\mathbb{X}^{\rm SG}_{l^{(h)}}\\}$,
a sparse linear filter $[R^{(h)}]^{-1}$ which converts the dense input to
sparse local hierarchical features, bias $\bm{\mu}^{(h)}$, and i.i.d. standard
Gaussian weights. Moreover, sparse linear filters
$\\{[R^{(h)}]^{-1}\\}_{h=1}^{H}$ are fixed for a specific hierarchical
expansion. Therefore, $\\{[R^{(h)}]^{-1}\\}_{h=1}^{H}$ can be computed and
stored in advance. Due to Theorems 1 and 2, the space complexity for storing
$\\{[R^{(h)}]^{-1}\\}_{h=1}^{H}$ is linear in the number of inducing points
and the time complexity of any operation in layer $h$ of DTMGP is
$\mathcal{O}([\log m^{(h)}]^{2W^{(h-1)}-1})$.

###  4.4 Variational Inference

The proposed DTMGP is mathematically equivalent to a deep neural network with
biases $\\{\bm{\mu}^{(h)}\\}_{h=1}^{H}$ and independently and normally
distributed weights $\\{\mathbb{Z}^{(h)}\\}_{h=1}^{H}$. So we can parameterize
DTMGP as a BNN with independently and normally distributed weights.
Specifically, we regard the variances and means of all Gaussian weights, and
the biases as parameters of the BNN. A DTMGP can be represented as a
stochastic process $f_{\bm{\theta}}$ in the parameterized family
$\mathcal{F}_{\bm{\Theta}}$:

$\displaystyle\mathcal{F}_{\bm{\Theta}}=\bigg{\\{}\left[\big{(}\bm{\Sigma}^{(H)}\odot\mathbb{Z}^{(H)}+\mathbb{m}^{(H)}\big{)}\bm{\phi}^{(H)}+\bm{\mu}^{(H)}\right]\circ\cdots\circ\left[\big{(}\bm{\Sigma}^{(1)}\odot\mathbb{Z}^{(1)}+\mathbb{m}^{(1)}\big{)}\bm{\phi}^{(1)}+\bm{\mu}^{(1)}\right]:$
$\displaystyle\quad\quad\quad\big{\\{}\bm{\Sigma}^{(h)},\mathbb{m}^{(h)},\bm{\mu}^{(h)}\big{\\}}_{h=1}^{H}\in\bm{\Theta}\bigg{\\}},$

where $\bm{\Theta}$ denotes the set of parameters

$\bm{\Theta}\coloneqq\bigg{\\{}\bm{\Sigma}^{(h)}\in\mathbb{R}_{+}^{W^{(h)}\times
m^{(h)}},\mathbb{m}^{(h)}\in\mathbb{R}^{W^{(h)}\times
m^{(h)}},\bm{\mu}^{(h)}\in\mathbb{R}^{W^{(h)}},h=1,\cdots,H\bigg{\\}},$

and $A\odot B$ denotes the entry-wise multiplication of matrices $A$ and $B$,
i.e., the Hadamard product. In other words, the family
$\mathcal{F}_{\bm{\Theta}}$ consists of BNNs with activations
$\\{\bm{\phi}^{(h)}\\}_{h=1}^{H}$, parameterized biases, and normally
distributed weights.

Given observations $(\mathbb{X},\mathbb{Y})$, our goal is to search for the
parameters $\bm{\theta}^{*}$ that maximizes the marginal likelihood
$\mathbb{P}(\mathbb{Y}\big{|}\mathbb{X},\bm{\theta})$. The DTMGP
$f_{\bm{\theta}^{*}}$ can best interpret how $(\mathbb{X},\mathbb{Y})$ is
generated. Computing $\mathbb{P}(\mathbb{Y}|\mathbb{X},\bm{\theta})$ of any
DGP involves intractable integral as shown in (3), but we can obtain a
tractable lower bound using VI. VI first assigns a prior
$f_{\tilde{\bm{\theta}}}$ with some
$\tilde{\bm{\theta}}=\\{\tilde{\bm{\Sigma}}^{(h)},\tilde{\mathbb{m}}^{(h)},\tilde{\bm{\mu}}^{(h)}\\}_{h=1}^{H}$
to the underlying stochastic process generating $(\mathbb{X},\mathbb{Y})$ and
then maximizes the following _evidence lower bound_ (ELBO):

$\mathcal{E}(\bm{\theta})=\mathbb{E}_{f_{\bm{\theta}}}\bigl{[}\log\Pr(\mathbb{Y}\big{|}\mathbb{X},\bm{\theta})\bigr{]}-\mathcal{D}_{KL}\bigl{(}f_{\bm{\theta}}\big{\|}f_{\tilde{\bm{\theta}}}\bigr{)},$
(13)

which is the lower bound of log-likelihood
$\log\Pr\big{(}\mathbb{Y}\big{|}\mathbb{X},\bm{\theta}\big{)}$.

In our framework, the term
$\mathbb{E}_{f_{\bm{\theta}}}\bigl{[}\log\Pr(\mathbb{Y}\big{|}\mathbb{X},\bm{\theta})\bigr{]}$
is called the _negative energy_ and can be explicitly written as

$\mathbb{E}_{f_{\bm{\theta}}}\bigl{[}\log\Pr(\mathbb{Y}\big{|}\mathbb{X},\bm{\theta})\bigr{]}=\mathbb{E}_{\mathbb{Z}^{(1)},\cdots,\mathbb{Z}^{(H)}}\bigl{[}\log\Pr(\mathbb{Y}\big{|}\mathbb{X},\mathbb{Z}^{(1)},\cdots,\mathbb{Z}^{(H)},\bm{\theta})\bigr{]}.$

This term now can be efficiently estimated by stochastic optimization based on
Monte Carlo methods, such as doubly-stochastic approximation (Hensman and
Lawrence, 2014; Dai et al., 2014; Salimbeni and Deisenroth, 2017). The term
$\mathcal{D}_{KL}\bigl{(}f_{\bm{\theta}}\big{\|}f_{\tilde{\bm{\theta}}}\bigr{)}$
is the _KL-divergence_ between $f_{\bm{\theta}}$ and the prior
$f_{\tilde{\bm{\theta}}}$ and it can be explicitly written

$\mathcal{D}_{KL}\bigl{(}f_{\bm{\theta}}\big{\|}f_{\tilde{\bm{\theta}}}\bigr{)}=\frac{1}{2}\sum_{h=1}^{H}\sum_{i=1}^{W^{(h)}}\sum_{j=1}^{m^{(h)}}\frac{\big{|}[\mathbb{m}^{(h)}]_{i,j}-[\tilde{\mathbb{m}}^{(h)}]_{i,j}\big{|}^{2}}{[\tilde{\bm{\Sigma}}^{(h)}]^{2}_{i,j}}+\frac{[\bm{\Sigma}^{(h)}]^{2}_{i,j}}{[\tilde{\bm{\Sigma}}^{(h)}]^{2}_{i,j}}-1-\log\frac{[\bm{\Sigma}^{(h)}]_{i,j}^{2}}{[\tilde{\bm{\Sigma}}^{(h)}]_{i,j}^{2}}.$

The computation of $\mathcal{E}$ is efficient at any
$\bm{\theta}\in\bm{\Theta}$, the training process of DTMGP is simply defined
as

$\bm{\theta}^{*}=\arg\max_{\bm{\theta}\in\bm{\Theta}}\left\\{\mathbb{E}_{f_{\bm{\theta}}}\bigl{[}\log\Pr(\mathbb{Y}\big{|}\mathbb{X})\bigr{]}-\mathcal{D}_{KL}\bigl{(}f_{\bm{\theta}}\big{\|}f_{\tilde{\bm{\theta}}}\bigr{)}\right\\}.$
(14)

This optimization problem can also be treated as a penalized regression
problem where the negative energy term is the model fitness level and the KL-
divergence term is the penalty of being far from the DTMGP prior. Equation
(14) can be solved efficiently by using the automatic differentiation
technique (Neidinger, 2010; Baydin et al., 2018).

##  5 Simulation Studies

In this section, we run simulations on systems with stochastic outputs to
access DTMGP’s capacity in modeling stochastic processes. More specifically,
in each experiment, output $Y$ of the underlying system at each input
$\mathbb{x}$ follows an _unknown_ distribution $F(\cdot|\mathbb{x})$:

$Y(\mathbb{x})\sim F(y|\mathbb{x}).$

A training set $(\mathbb{X}_{\rm train},\mathbb{Y}_{\rm train})$ is
independently collected where $\mathbb{Y}_{\rm train}=Y(\mathbb{X}_{\rm
train})$ is the realization of the underlying system on inputs
$\mathbb{X}_{\rm train}$. We then train competing models using training set
$(\mathbb{X}_{\rm train},\mathbb{Y}_{\rm train})$.

Let $\widehat{Y}(\cdot)$ denote a trained model. We first choose a test set
$\mathbb{X}_{\rm test}$ from the input space. For each input
$\mathbb{x}\in\mathbb{X}_{\rm test}$, we sample $100$ independent realizations
from the true underlying system and from the trained model to get data sets
$\\{{Y}_{i}(\mathbb{x})\\}_{i=1}^{100}$ and
$\\{\widehat{Y}_{i}(\mathbb{x})\\}_{i=1}^{100}$, respectively. We then use the
two-sample _Kolmogorov–Smirnov_ (KS) statistic:

${D}_{\mathbb{x}}=\sup_{y}\left|\frac{1}{100}\sum_{i=1}^{100}\mathbb{1}_{\\{Y_{i}(\mathbb{x})\leq
y\\}}-\frac{1}{100}\sum_{i=1}^{100}\mathbb{1}_{\\{\widehat{Y}_{i}(\mathbb{x})\leq
y\\}}\right|$

to quantify the similarity between $Y$ and $\widehat{Y}$ on input
$\mathbb{x}$. The smaller $D_{\mathbb{x}}$ is, the closer the two
distributions are. We access the overall performance of $\widehat{Y}(\cdot)$
via the following averaged two-sample KS statistic over test set
$\mathbb{X}_{\rm test}$:

${D}=\frac{1}{|\mathbb{X}_{\rm test}|}\sum_{\mathbb{x}\in\mathbb{X}_{\rm
test}}D_{\mathbb{x}}.$

We repeat each experiment $R$ times and call each a macro-replication from
which we obtain ${D}_{r}$, $r=1,\ldots,R$. Finally, we report the mean and
standard deviation

$\displaystyle\bar{D}=\frac{1}{R}\sum_{r=1}^{D}D_{r},$
$\displaystyle\hat{\sigma}=\biggl{[}\frac{1}{R}\sum_{r=1}^{R}\biggl{(}{D}_{r}-\bar{D}\biggr{)}^{2}\biggr{]}^{1/2}.$

In all experiments, activations of DTMGP are hierarchical features generated
from TMGPs with the Laplace covariance function
$k(\mathbb{x},\mathbb{x}^{\prime})=\exp\\{-\frac{\|\mathbb{x}-\mathbb{x}^{\prime}\|_{1}}{d}\\}$.
We compare DTMGP with the following existing models:

  1. 1.

GP: standard Gaussian process with squared exponential covariance function.

  2. 2.

BNN-ReLU (Schmidt et al., 1992; Graves, 2011; Blundell et al., 2015): ReLU
deep neural networks (Glorot et al., 2011) with parameterized Gaussian
distributed weights and biases. VI is used for training BNN-ReLU.

  3. 3.

DGP-RFF (Cutajar et al., 2017): A Random Fourier feature expansion of DGP with
each layer represented by Gaussian processes with the Gaussian kernel
$\exp\\{-\sum_{j=1}^{d}w_{j}|x_{j}-x^{\prime}_{j}|^{2}\\}$. DGP-RFF is trained
with VI approximation provided in Cutajar et al. (2017).

  4. 4.

DGP-VEC (Sauer et al., 2022b) DGP where each activation is GP with
Matérn-$1/2$ kernel. Posterior of the model is approximated through the
incorporation of the Vecchia approximation (Katzfuss and Guinness, 2021).

Note that all deep models can be represented by the same architecture and the
difference among them lies in their activations. For DTMGP, BNN-ReLU, and DGP-
RFF, they are BNN with different deterministic activations and random weights.
For DGP-VEC, it is composition of GPs. Therefore, for each data set, we adopt
the same architectures for fair comparisons. In other words, except for the
specific activations, all deep models are with equal number of layers, widths,
weights, and biases.

We conduct experiments with two simulation models: a two-dimensional non-
Gaussian random field in Section 5.1, and the expected revenue of a
13-dimensional Stochastic Activity Network problem in Section 5.2. All the
experiments are implemented on a computer with macOS, 3.3 GHz Intel Core i5
CPU, and 8 GB of RAM (2133Mhz).

###  5.1 A Non-Gaussian Random Field

In this section, we use the following non-Gaussian and non-stationary random
field in two dimensions to access the sample and computational efficiency of
the proposed methodology:

$\displaystyle
Y(\mathbb{x})=\frac{1}{1+\exp\\{B(\mathbb{x})\\}},\quad\mathbb{x}\in[0,1]^{2}$

where $B(\mathbb{x})$ is a _Brownian sheet_ , which is defined as a zero-mean
GP with covariance function
$\prod_{d=1}^{2}\big{(}1+\min\\{x_{d},x_{d}^{\prime}\\}\big{)}$.

Figure 7:  A sample path of non-Gaussian random field
$\frac{1}{1+\exp\\{B(\mathbb{x})\\}}$. Right: architecture shared by all deep
models to capture the random field.

In this experiment, $\mathbb{X}_{\rm train}$ is randomly collected from
$[0,1]^{2}$ and $\mathbb{Y}_{\rm train}$ is from independently sample paths of
$B(\mathbb{x})$ on $\mathbb{X}_{\rm train}$. We investigate the performance of
each method as the sizes of the training set $(\mathbb{X}_{\rm
train},\mathbb{Y}_{\rm train})$ increase. For each training set, we randomly
select $n_{\rm test}=100$ points from $[0,1]^{2}$ as the test set, denoted as
$\mathbb{X}_{\rm test}$. We run $R=20$ macro-replications and compute
$\hat{D}$ and the associated standard deviation $\hat{\sigma}$ for each
method.

All deep models are with two hidden layers where the first layers consists of
181 activations and 1 output, and the second layer consist of 127 activations
and 1 output. Figure 7 shows the realization of the random field and
architecture shared by all deep models. We let the i.i.d. prior of each
coefficient in DTMGP be normal distributed with mean 0 and variance 1.

Our experiments reveal that DGP-VEC has a significantly longer training time
compared to all other models. Specifically, Table 1 demonstrates that while
DTMGP has the fastest convergence rate to stable training error, DGP-VEC
exhibits the slowest. Consequently, to ensure a fair comparison, we report our
results at different training times. First, we report the KS statistics of all
models when DTMGP’s training error begins to converge, and then we report the
KS statistics of all models when DGP-VEC’s training error begins to converge.

model |  n=40 |  n=80 |  n=120 |  n=160 |  n=200
---|---|---|---|---|---
DTMGP |  28.7 sec |  53.5 sec |  89.8 sec |  138.1 sec |  166.4 sec
DGP-VEC |  44.2 sec |  97.9 sec |  255.3 sec |  384.2 sec |  516.7 sec
Table 1: Time Required to Convergence of Training Error

Figure 8:  Left: $\bar{D}$ of each model, the shaded areas and error bars
represent 1 std $\hat{\sigma}$ of all models. Left: results reported when KS
statistics of DTMGP begins to converge. Right: results reported when KS
statistics of DGP-VEC begins to converge or traing time reaches 4 hours.

As shown in Figure 7, when given enough training time, DGP-VEC outperforms all
other models. However, it’s worth noting that the performance of DTMGP is
still comparable to that of DGP-VEC. If we limit the training time for all
models to end once the training error of DTMGP begins to converge, we find
that DTMGP - which has the fastest convergence rate - outperforms all other
models, including DGP-VEC. This is despite the fact that DGP-VEC requires a
training time that is three to four times longer than DTMGP, as it relies on
MCMC for training. As it is widely acknowledged that VI is more efficient than
MCMC, this serves as an important reminder of the tradeoffs we must consider
when selecting a training method.

###  5.2 Stochastic Activity Network

In this subsection, we consider the stochastic activity network (SAN) where
the arcs are labeled from 1 through 13. The detailed explaition of SAN is
available in Avramidis and Wilson (1996). As shown in the left of Figure 9 ,
each arc $i$ in the SAN is associated with a task with random duration $D_{i}$
and task durations are mutually independent. Suppose that $D_{i}$ is
exponentially distributed with mean $X_{i}$ for each $i$. Suppose that we can
control $X_{i}>0$ for each $i$, but there is an associated cost. In
particular, the overall cost is defined as

$C(\mathbb{X})=T(\bm{\mathbb{X}})+f(\mathbb{X})$

where $\mathbb{X}=(X_{1},\cdots,X_{13})$, $T(\mathbb{X})$ is the (random)
duration of the longest path from a to i, and
$f(\mathbb{X})=\sum_{i=1}^{13}X_{i}^{-1}$. Closed form of $C(\mathbb{X})$ is
unknown but MATLAB simulator of this problem is available in the SimOpt
library (Pasupathy and Henderson, 2011).

In this experiment, data are collected from a maximin Latin hypercube design
(LHD) that maximizes the minimum distance between points (van Dam et al.,
2007). The LHD consists of 5000 sample points from the cube $[0.5,5]^{13}$. At
each sample point, $m$ simulation replications are run with $m=2,4,\cdots,20$.
In order to select reasonable prior for the target stochastic process, we
first normalize all the output so that all output data in the training set and
testing set are distributed on $[0,1]$. We then run $m=10$ replications at
each design point during the training process to estimate the variance at each
point. Based on the sample estimates, we let the i.i.d. prior of each
coefficient in DTMGP be normal distributed with mean 0 and variance 0.04.

We examine the performance of each method as the number of replications $m$
increase. Given $m$, we construct the test set $\mathbb{X}_{\rm test}$ of size
$n_{\rm test}=100$ by random samples from $[0.5,5]^{13}$. We compute $\hat{D}$
of each method, its associated standard deviation $\hat{\sigma}$ based on
$R=20$ macro-replications.

Figure 9:  Left: structure of SAN. Middle: i.i.d. samples of $C(\mathbb{X})$
with $X_{i}=2.5$, $i\neq 7$ and $X_{7}\in[0.5,5]$. Right: architecture shared
by all competing deep GPs to capture $C(\mathbb{X})$.

Similar to section 5.1, all the competing models in this experiment share the
same architecture. All deep models have two hiden layers. The first layer
consists of 391 activations and its output is 10-dimensional. The second layer
consists of 241 activations and its output is 1-dimensional. Such an
architecture is flexible enough to capture the random process.

model |  m=4 |  m=8 |  m=12 |  m=16 |  m=20
---|---|---|---|---|---
DTMGP |  9.5 min |  17.8 min |  24.3 min |  36.7 min |  46.2 min
DGP-VEC |  1.8 hr |  3.2 hr |  $>4$ hr |  $>4$ hr |  $>4$ hr
Table 2: Time Required to Convergence of Training Error

Figure 10:  $\bar{D}$ of each model, the shaded areas and error bars represent
1 std $\hat{\sigma}$ of all models. Left: results reported when KS statistics
of DTMGP begins to converge. Right: results reported when KS statistics of
DGP-VEC begins to converge.

As before, the training time of DGP-VEC is much longer than other models. When
$m>10$, the training time required for DGP-VEC is longer than 4 hours, which
is unaccetable compared with DTMGP. Therefore, to ensure a fair comparison, we
also report our results at different training times. First, we report the KS
statistics of all models when DTMGP’s training error begins to converge, and
then we report the KS statistics of all models when DGP-VEC’s training error
begins to converge or the training time reaches 4 hours.

As demonstrated in Figure 10, DGP-VEC outperforms all other models when given
sufficient training time. However, DGP-VEC’s training time exceeds 4 hours
when $m>10$ or when the total training data size is greater than 50000.
Therefore, we only present the results of DGP-VEC for $m\leq 10$. On the other
hand, DTMGP demonstrates comparable performance, with a training time of
roughly one-tenth that required by DGP-VEC. Similarly to Section 5.1, when we
restrict the training time for all models to stop once the training error of
DTMGP starts to stabilize, we observe that DTMGP outperforms all other models,
and DGP-VEC performs even worse than BNN-ReLU. This is because the training
time is not sufficient for MCMC.

##  6 Experiment on Real Data

This section aims to evaluate the performance of the proposed methodology on a
real dataset. To this end, we utilize the RC-49 dataset (Ding et al., 2021),
which is a synthetic collection of 49 3-D chair models rendered at yaw angles
ranging from 0°to 90°. The RC-49 dataset consists of $n=176400$ data pairs
$(\mathbb{X},\mathbb{Y})=\\{(x_{i},\mathbb{y}_{i})\\}_{i=1}^{n}$ where, for
any $i=1,\ldots,n$, input $x_{i}\in[0,90)$ is a specific yaw angles and
outputs $\mathbb{y}_{i}$ is a $64\times 64$ image of a chair rendered. The
objective is to use the Bayesian formula to train a generative model, which is
conditioned on the input variable $x$. When provided with input yaw angle $x$,
the trained model should generate a random matrix output that is similar to
the data $\mathbb{Y}$ at the same angle. In essence, the generative problem is
a stochastic process $Y(x)\in\mathbb{R}^{64\times 64}$ where $x\in[0,90)$, and
our aim is to employ a deep GP model to reconstruct $Y(\cdot)$ from observed
data.

We compare DTMGP with BNN-ReLU and DGP-RFF introduced in Section 5. Other deep
GP models are omitted due to their limited capacity for handling the large
data size and complex model architecture. All competing deep GP models share
the following architecture. We first apply the embedding technique on input
$x$, which uses a linear transformation $A$ to map the single variable input
$x\in[0,90)$ to a 100-dimensional vector
$\mathbb{x}^{\prime}\in\mathbb{R}^{100}$. We then adopt a four-layer BNN
architecture where $\mathbb{x}^{\prime}$ is the input layer, the first hidden
layer consists of 201 activations and a 256-dimensional output, the second
hidden layer consists of 513 activations and a 512-dimensional output, and the
third hidden layer consists of 1025 activations and the final output is a
$64\times 64$ random matrix. Specifically, the architecture can be written as
follows:

$\displaystyle
f^{(1)}_{j}=\sum_{i=1}^{201}Z_{i,j}^{(1)}\phi_{i}^{(1)}(A\mathbb{x}),\quad
j=1,\cdots,256,\ Z_{i,j}^{(1)}\sim\mathcal{N}(0,1),A\in\mathbb{R}^{100\times
1},$ $\displaystyle
f^{(2)}_{j}=\sum_{i=1}^{513}Z_{i,j}^{(2)}\phi^{(2)}_{i}(f^{(1)}_{1},\cdots,f^{(1)}_{256}),\quad
j=1,\cdots,512,\ Z_{i,j}^{(2)}\sim\mathcal{N}(0,1),$
$\displaystyle\hat{Y}_{l,j}=\sum_{i=1}^{1025}Z_{i,l,j}^{(3)}\phi_{i}^{(3)}(f^{(2)}_{1},\cdots,f^{(2)}_{512}),\quad
l,j=1,\cdots,28,\ Z_{i,l,j}^{(3)}\sim\mathcal{N}(0,1).$

Figure 11 provides an illustration of the above architecture.

Figure 11:  All competing deep GPs share the same architecture.

We use 80000 samples from the RC-49 dataset to train DTMGP and its competitors
BNN-RELU and DGP-RFF. We train the models by mini-batch SGD and, in each
epoch, we sample their random outputs at input $x$, for $x\in[0,90)$. Because
the output in this experiment is a $64\times 64$ random matrix and we are
unable to sample data from the true distribution of RC-49, the previous KS
statistics cannot be used to access the accuracy. Instead, we directly examine
the output images of each model to evaluate their performances.

Figure 12:  Results at the three epoch (left), six epoch (middle), and nine epoch (right) during training with input $x=10,20,\cdots,90$.  Figure 13:  Generated images of DTMGP with x=60 Model |  1 Epoch |  3 Epoch |  6 Epoch |  9 Epoch
---|---|---|---|---
DTMGP |  1099.3 sec |  3115.2 sec |  6117.6 sec |  8953.1 sec
BNN-ReLU |  971.5 sec |  2511.9 sec |  4997.2 sec |  7301.1 sec
DGP-RFF |  1081.4 sec |  2987.7 sec sec |  5823.4 sec |  8682.8 sec
Table 3: Training Time

Conditional samples generated by DTMGP, BNN-ReLU, and DGP-RFF at the $3^{\rm
rd}$, $6^{\rm th}$ and $9^{\rm th}$ epochs are shown in Figure 12. The results
demonstrate that DTMGP generates samples that are closer to the actual data
and performs the best, while other models fail to reconstruct the conditional
samples. This is because, as observed in previous experiments, the
hierarchical structure of DTMGP can effectively capture the local features
while the ReLU or cosine activations of BNN-ReLU or DGP-RFF lack the ability
to do so.

As depicted in Figure 13, DTMGP exhibits a faster convergence rate towards
stability with more consistent generated results. The rapid convergence
suggests that the model can be trained in few epochs, while the consistent
performance implies that the model is reliable. This may be attributed to the
sparsity introduced by our expansion approach, which facilitates feature
detection and speeds up parameter determination.

Table 3 demonstrates that the training time per epoch of both BNN-ReLU and
DGP-RFF is faster than that of DTMGP. This difference could be primarily
attributed to the fact that BNN-ReLU and DGP-RFF are built using built-in
functions and packages in the MATLAB, while our code for DTMGP does not take
this advantage. Nonetheless, the training time per epoch of DTMGP is still
comparable to that of BNN-ReLU and DGP-RFF.

##  7 Conclusions and Discussion

We utilize hierarchical features to expand DGPs, which comprise TMGPs. Our
expansion is known as DTMGP and is computationally efficient, resulting in a
sparse representation of DGPs. The hierarchical nature of the features means
that their supports are either nested or disjoint, enabling a poly-logarithmic
number of activations throughout a DTMGP to have non-zero output conditioned
on any input or operation. Furthermore, hierarchical features can effectively
capture local features from inputs, resulting in high performance in
prediction and generative problems, as demonstrated in our experiments. In
comparison to existing DGP models, DTMGP’s sparsity leads to efficient
inference and training, while generated instances from DTMGP are also
relatively accurate under the KS test.

The current paper can potentially be extended in several ways. Firstly, how to
apply hierarchical expansion on DGPs with more general structures, such as
DGPs with Mátern covariance functions, can be studied in future research.
Secondly, it was pointed out by Stein (2014) that, approximating the
likelihood function of GP by low-rank methods, such as inducing point and
random features, may have an adverse effect on the performance. How large this
effect is for the proposed DGP models should be studied in further
investigation.

## Acknowledgements

The authors are grateful to the editors and referees for very helpful
comments.

## Funding

Tuo’s research is supported by NSF DMS-1914636 and CCF-1934904 and 2022 Texas
A&M Institute of Data Science Career Initiation Fellow Program. Part of Ding’s
work was conducted when he was a postdoctoral researcher at Texas A&M
University and was partially supported by the Texas A&M Institute of Data
Science (TAMIDS) Postdoctoral Fellowships.

Appendix

##  Appendix A Sparse Grid Designs with Required Order

The class of SG adopted in this work is called _Hyperbolic-cross_ SG, which is
the union of full grids (FG) with dyadic structures. Without loss of
generality, we can assume that all design points are collected from the
hypercube $(0,1)^{d}$, and we start constructing SGs satisfying the required
order with the one-dimensional case, which is also called dyadic point set.

A one-dimensional level-$l$ dyadic point set $\mathbb{X}_{l}$ with increasing
order is simply defined as $\mathbb{X}_{l}=\\{2^{-l},2\cdot
2^{-l},\cdots,1-2^{-l}\\}$. However, this order cannot result in a sparse
representation of hierarchical expansion. We need to sort the points according
to the their levels. Firstly, given a level $l\in\mathbb{N}$, we define the
following set consisting only of odd numbers:

$\bm{\rho}(l)=\\{1,3,5,\cdots,2^{l}-1\\}.$

Then given a level-$l$ dyadic point set $\mathbb{X}_{l}$ and a level-$(l-1)$
dyadic point set $\mathbb{X}_{l-1}$ , we can define the sorted incremental set

$\mathbb{D}_{l}=\mathbb{X}_{l}-\mathbb{X}_{l-1}=\big{\\{}i2^{-l}:i\in\bm{\rho}(l)\big{\\}},$

with $\mathbb{D}_{0}=\emptyset$. It is straightforward to check that

$\displaystyle\mathbb{X}_{l}=\bigcup_{\ell=1}^{l}\mathbb{D}_{\ell},\quad\mathbb{D}_{l}\bigcap\mathbb{D}_{k}=\emptyset,\
\text{if}\ l\neq k.$

For any incremental set $\mathbb{D}_{l}$, we can label any point in
$\mathbb{D}_{l}$ by $x_{l,i}$, $i\in\rho(l)$ so that the label is unique for
any point: $x_{l,i}\neq x_{l^{\prime},i^{\prime}}$ if and only if
$(l,i)\neq(l^{\prime},i^{\prime})$. Now we can define the following dyadic set
with the required order:

$\mathbb{X}_{l}^{*}=\big{[}\mathbb{D}_{1},\mathbb{D}_{2},\cdots,\mathbb{D}_{l}\big{]}=[x_{\ell,i}:\ell\leq
l,i\in\bm{\rho}(\ell)].$

For example,
$\mathbb{X}^{*}_{3}=[\frac{1}{2},\frac{1}{4},\frac{3}{4},\frac{1}{8},\frac{3}{8},\frac{5}{8},\frac{7}{8}]$
with $\mathbb{D}_{1}=[\frac{1}{2}]$,
$\mathbb{D}_{2}=[\frac{1}{4},\frac{3}{4}]$, and
$\mathbb{D}_{3}=[\frac{1}{8},\frac{3}{8},\frac{5}{8},\frac{7}{8}]$,
respectively.

We also define the Cartesian product of one-dimensional sorted
$\mathbb{X}^{*}_{l}$’s for later use in Algorithms and proofs. We call this
sorted Cartesian product full grid (FG) $\mathbb{X}^{*}_{\mathbb{l}}$ labeled
by $\mathbb{l}=(l_{1},\cdots,l_{d})\in\mathbb{N}^{d}$:

$\mathbb{X}_{\mathbb{l}}^{*}=\bigtimes_{j=1}^{d}\mathbb{X}^{*}_{l_{j}}=\big{[}(x_{\ell_{1},i_{1}},x_{\ell_{2},i_{2}},\cdots,x_{\ell_{d},i_{d}}):{\ell_{j}\leq
l_{j}},i_{j}\in\bm{\rho}(\ell_{j}),j=1,\cdots,d\big{]}.$

For simplicity, we define the Cartesian product
$\bm{\rho}(\mathbb{l})=\times_{j=1}^{d}\bm{\rho}(l_{j})$ for any
$\mathbb{l}\in\mathbb{N}^{d}$ so that any points in
$\mathbb{X}^{*}_{\mathbb{l}}$ can be represented as
$\mathbb{x}_{\mathbb{l},\mathbb{i}}=(x_{l_{1},i_{1}},x_{l_{2},i_{2}},\cdots,x_{l_{d},i_{d}})$
for any $\mathbb{l}\in\mathbb{N}^{d}$ and
$\mathbb{i}\in\bm{\rho}(\mathbb{l})$. Then we have a more compact
representation of $\mathbb{X}^{*}_{\mathbb{l}}$:

$\mathbb{X}_{\mathbb{l}}^{*}=\big{[}\mathbb{x}_{\bm{\ell},\mathbb{i}}:\bm{\ell}\leq\mathbb{l},\mathbb{i}\in\bm{\rho}(\mathbb{\ell})\big{]}.$

For example, the sorted FG $\mathbb{X}^{*}_{(1,2)}$ is

$\mathbb{X}^{*}_{(1,2)}=[(\frac{1}{2},\frac{1}{2}),(\frac{1}{2},\frac{1}{4}),(\frac{1}{2},\frac{3}{4})].$

Finally, a level-$l$ SG satisfying our required order is the union of sorted
FGs as follows:

$\mathbb{X}^{\rm
SG}_{l}=\bigcup_{|\mathbb{l}|=l+d-1}\mathbb{X}^{*}_{\mathbb{l}}=\big{[}\mathbb{x}_{\mathbb{l},\mathbb{i}}:|\mathbb{l}|\leq
l,\mathbb{i}\in\bm{\rho}(\mathbb{l})\big{]},$

where $|\mathbb{l}|$ denote the $l_{1}$ norm $\sum_{j=1}^{d}|l_{j}|$. To be
more specific, a point $\mathbb{x}_{\mathbb{l},\mathbb{i}}\in\mathbb{X}^{\rm
SG}_{l}$ must be ahead of any $\mathbb{x}_{\mathbb{l}^{\prime}}$ with
$|\mathbb{l}^{\prime}|>|\mathbb{l}|$. For example, the two-dimensional sorted
SG $\mathbb{X}^{\rm SG}_{2}$ is

$\mathbb{X}^{\rm
SG}_{2}=\big{[}\underbrace{(\frac{1}{2},\frac{1}{2})}_{\text{level}\
|(1,1)|=2},\underbrace{(\frac{1}{2},\frac{1}{4}),(\frac{1}{2},\frac{3}{4})}_{\text{level}\
|(1,2)|=3},\underbrace{(\frac{1}{4},\frac{1}{2}),(\frac{3}{4},\frac{1}{2}}_{\text{level}\
|(2,1)|=3})\big{]}.$

Unlike FG, whose size increases exponentially in dimension $d$, the size of
SGs increases relatively mildly in $d$. Lemma 3.6 in Bungartz and Griebel
(2004) stipulates that the sample size of a $d$-dimensional level-$l$
$\mathbb{X}_{l}^{\rm SG}$ is given by

$|\mathbb{X}_{l}^{\rm
SG}|=\sum_{\ell=0}^{l-1}2^{\ell}\binom{\ell+d-1}{d-1}=2^{l}\cdot\left(\frac{l^{d-1}}{(d-1)!}+\mathcal{O}(l^{d-2})\right)=\mathcal{O}(2^{l}l^{d-1}).$
(15)

##  Appendix B Algorithms

We first introduce how to compute $R^{-1}$, with $R$ the Cholesky
decomposition of $k(\mathbb{X}^{*}_{l},\mathbb{X}^{*}_{l})$, $k$ is a one-
dimensional Markov kernel and $\mathbb{X}^{*}_{l}$ is a one-dimensional sorted
level-$l$ dyadic point set. To present the algorithm for computing the inverse
of Cholesky decomposition, we also label the entries on matrices by points in
$\mathbb{X}^{*}_{l}$. For example, $[R^{-1}]_{x,x^{\prime}}$ represents entry
with row index corresponding to point $x\in\mathbb{X}^{*}_{l}$ and column
index corresponding to point $x^{\prime}\in\mathbb{X}^{*}_{l}$.

Input : Markov kernel $k$, level-$l$ dyadic points $\mathbb{X}_{l}^{*}$

Output : $R^{-1}$

1 Initialize
$R^{-1}\leftarrow\mathbb{0}\in\mathbb{R}^{(2^{l}-1)\times(2^{l}-1)}$,
$\mathbb{N}=\\{-\infty,\infty\\}$, define $k(\pm\infty,x)=0,\forall x$

2 for  _$\ell\leftarrow 1$ to $l$_ do

3 for  _$i\in\bm{\rho}(\ell)$_ do

4 search the closest left neighbor $x_{left}$ and right neighbor $x_{right}$
of $x_{\ell,i}$ in $\mathbb{N}$

5 Solve $c_{1},c_{2}$, and $c_{3}$ for the following system:  $\displaystyle
c_{1}k(x_{left},x_{left})+c_{2}k(x_{\ell,i},x_{left})+c_{3}k(x_{right},x_{left})=0,$
(16) $\displaystyle
c_{1}k(x_{left},x_{right})+c_{2}k(x_{\ell,i},x_{right})+c_{3}k(x_{right},x_{right})=0,$
$\displaystyle\mathbb{E}\big{[}\big{(}c_{1}\mathcal{G}(x_{left})+c_{2}\mathcal{G}(x_{\ell,i})+c_{3}\mathcal{G}(x_{right})\big{)}^{2}\big{]}=1$

6 If $x_{left}\neq-\infty$, $[R^{-1}]_{x_{left},x_{\ell,i}}=c_{1}$

7 If $x_{right}\neq\infty$, $[R^{-1}]_{x_{right},x_{\ell,i}}=c_{3}$

8 Let $[R^{-1}]_{x_{\ell,i},x_{\ell,i}}=c_{2}$,
$x_{\ell,i}\rightarrow\mathbb{N}$

9 end for

10

11 end for

12Return $R^{-1}$

Algorithm 1 Computing $R^{-1}$ for Markov kernel $k$ and point set
$\mathbb{X}_{l}^{*}$

Now let $k=\prod_{j=1}^{d}k_{j}$ be a $d$-dimensional TMK and
$\mathbb{X}^{*}_{\mathbb{l}}=\bigtimes_{j=1}^{d}\mathbb{X}^{*}_{l_{j}}$ be a
sorted FG. Then, $R^{-1}_{\mathbb{l}}$, with $R_{\mathbb{l}}$ the Cholesky
decomposition of $k(\mathbb{X}^{*}_{\mathbb{l}},\mathbb{X}^{*}_{\mathbb{l}})$
can be directly calculated as

$R_{\mathbb{l}}^{-1}=\bigotimes_{j=1}^{d}R^{-1}_{l_{j}},$ (17)

where $\bigotimes$ denotes the Kronecker product between matrices and
$R_{l_{j}}^{-1}$ is the inverse Cholesky decomposition yielded by Algorithm 1
with input $k_{j}$ and $\mathbb{X}^{*}_{l_{j}}$.

Let $R_{\mathbb{X}^{*}_{\mathbb{l}},\mathbb{X}^{*}_{\mathbb{l}}}$ represent a
sub-matrix of $R$ consisting of entries $R_{\mathbb{x},\mathbb{y}}$,
$\mathbb{x},\mathbb{y}\in\mathbb{X}^{*}_{\mathbb{l}}\subset\mathbb{X}^{\rm
SG}_{l}$. Now we can present the algorithm of constructing $R_{l}^{-1}$

Input : TMK $k$, level-$l$ SG $\mathbb{X}_{l}^{\rm{SG}}$

Output : $R_{l}^{-1}$

1 Initialize $R^{-1}\leftarrow\mathbb{0}\in\mathbb{R}^{m_{l}\times m_{l}}$

2 for  _all $\mathbb{l}\in\mathbb{N}^{d}$ with $l\leq|\mathbb{l}|\leq l+d-1$ _
do

3 Compute $R_{\mathbb{l}}^{-1}$ associated to
$(k,\mathbb{X}^{*}_{\mathbb{l}})$ via Algorithm 1 and (17)

4 Update $R_{l}^{-1}$ via
$\displaystyle[R_{l}^{-1}]_{\mathbb{X}^{*}_{\mathbb{l}},\mathbb{X}^{*}_{\mathbb{l}}}\leftarrow{}$
$\displaystyle[R_{l}^{-1}]_{\mathbb{X}^{*}_{\mathbb{l}},\mathbb{X}^{*}_{\mathbb{l}}}+(-1)^{l+d-1-|{\mathbb{l}}|}\binom{d-1}{l+d-1-|{\mathbb{l}}|}R_{\mathbb{l}}^{-1}$
(18)

5 end for

6Return $R_{l}^{-1}$

Algorithm 2 Computing $R_{l}^{-1}$ for TMK $k$ and level-$l$ SG
$\mathbb{X}^{\rm SG}_{l}$

We will prove the correctness of Algorithm 1 and Algorithm 2 in the
supplementary material.

## References

  * Ankenman et al. (2010) Ankenman, B., B. L. Nelson, and J. Staum (2010).  Stochastic kriging for simulation metamodeling.  Oper. Res. 58(2), 371–382.
  * Avramidis and Wilson (1996) Avramidis, A. N. and J. R. Wilson (1996).  Integrated variance reduction strategies for simulation.  Operations Research 44(2), 327–346.
  * Banerjee et al. (2008) Banerjee, S., A. E. Gelfand, A. O. Finley, and H. Sang (2008).  Gaussian predictive process models for large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70(4), 825–848.
  * Baydin et al. (2018) Baydin, A. G., B. A. Pearlmutter, A. A. Radul, and J. M. Siskind (2018).  Automatic differentiation in machine learning: A survey.  Journal of Machine Learning Research 18(153), 1–43.
  * Blei et al. (2017) Blei, D. M., A. Kucukelbir, and J. D. McAuliffe (2017, Apr).  Variational inference: A review for statisticians.  Journal of the American Statistical Association 112(518), 859–877.
  * Blundell et al. (2015) Blundell, C., J. Cornebise, K. Kavukcuoglu, and D. Wierstra (2015).  Weight uncertainty in neural network.  In International Conference on Machine Learning, pp. 1613–1622. PMLR.
  * Bui et al. (2016) Bui, T., D. Hernandez-Lobato, J. Hernandez-Lobato, Y. Li, and R. Turner (2016, 20–22 Jun).  Deep Gaussian processes for regression using approximate expectation propagation.  Volume 48 of Proceedings of Machine Learning Research, New York, New York, USA, pp. 1472–1481. PMLR.
  * Bungartz and Griebel (2004) Bungartz, H.-J. and M. Griebel (2004).  Sparse grids.  Acta Numerica 13, 147–269.
  * Cressie and Johannesson (2008) Cressie, N. and G. Johannesson (2008).  Fixed rank kriging for very large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70(1), 209–226.
  * Cutajar et al. (2017) Cutajar, K., E. V. Bonilla, P. Michiardi, and M. Filippone (2017, 06–11 Aug).  Random feature expansions for deep Gaussian processes.  Volume 70 of Proceedings of Machine Learning Research, International Convention Centre, Sydney, Australia, pp. 884–893. PMLR.
  * Dai et al. (2014) Dai, B., B. Xie, N. He, Y. Liang, A. Raj, M.-F. F. Balcan, and L. Song (2014).  Scalable kernel methods via doubly stochastic gradients.  In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger (Eds.), Advances in Neural Information Processing Systems 27, pp. 3041–3049. Curran Associates, Inc.
  * Dai et al. (2016) Dai, Z., A. Damianou, J. González, and N. Lawrence (2016).  Variational auto-encoded deep Gaussian processes.  Proceedings of the International Conference on Learning Representations (ICLR) 2016.
  * Damianou and Lawrence (2013) Damianou, A. and N. Lawrence (2013).  Deep Gaussian processes.  Volume 31 of Proceedings of Machine Learning Research, Scottsdale, Arizona, USA, pp. 207–215. PMLR.
  * Ding and Zhang (2020) Ding, L. and X. Zhang (2020).  Sample and computationally efficient simulation metamodeling in high dimensions.  arXiv preprint arXiv:2010.06802.
  * Ding et al. (2021) Ding, X., Y. Wang, Z. Xu, W. J. Welch, and Z. J. Wang (2021).  Ccgan: Continuous conditional generative adversarial networks for image generation.  In International conference on learning representations.
  * Fei et al. (2018) Fei, J., J. Zhao, S. Sun, and Y. Liu (2018).  Active learning methods with deep Gaussian processes.  In International Conference on Neural Information Processing, pp. 473–483. Springer.
  * Glorot et al. (2011) Glorot, X., A. Bordes, and Y. Bengio (2011, 11–13 Apr).  Deep sparse rectifier neural networks.  Volume 15 of Proceedings of Machine Learning Research, Fort Lauderdale, FL, USA, pp. 315–323. JMLR Workshop and Conference Proceedings.
  * Graves (2011) Graves, A. (2011).  Practical variational inference for neural networks.  In J. Shawe-Taylor, R. S. Zemel, P. L. Bartlett, F. Pereira, and K. Q. Weinberger (Eds.), Advances in Neural Information Processing Systems 24, pp. 2348–2356. Curran Associates, Inc.
  * Havasi et al. (2018) Havasi, M., J. M. Hernández-Lobato, and J. J. Murillo-Fuentes (2018).  Inference in deep Gaussian processes using Stochastic Gradient Hamiltonian Monte Carlo.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (Eds.), Advances in Neural Information Processing Systems, Volume 31. Curran Associates, Inc.
  * Hebbal et al. (2021) Hebbal, A., L. Brevault, M. Balesdent, E.-G. Talbi, and N. Melab (2021).  Bayesian optimization using deep Gaussian processes with applications to aerospace system design.  Optimization and Engineering 22(1), 321–361.
  * Hensman and Lawrence (2014) Hensman, J. and N. D. Lawrence (2014).  Nested variational compression in deep Gaussian processes.
  * Hinton (2010) Hinton, G. (2010).  A practical guide to training restricted Boltzmann machines (version 1).  Technical Report UTML TR 2010-003, University of Toronto 9.
  * Hinton et al. (2006) Hinton, G., S. Osindero, and Y.-W. Teh (2006, 08).  A fast learning algorithm for deep belief nets.  Neural Computation 18, 1527–54.
  * Katzfuss and Guinness (2021) Katzfuss, M. and J. Guinness (2021).  A general framework for vecchia approximations of gaussian processes.  Statistical Science 36(1), 124–141.
  * Ko and Kim (2021) Ko, J. and H. Kim (2021).  Deep Gaussian process models for integrating multifidelity experiments with non-stationary relationships.  IISE Transactions (just-accepted), 1–28.
  * Li et al. (2021) Li, Y., S. Rao, A. Hassaine, R. Ramakrishnan, D. Canoy, G. Salimi-Khorshidi, M. Mamouei, T. Lukasiewicz, and K. Rahimi (2021).  Deep Bayesian Gaussian processes for uncertainty estimation in electronic health records.  Scientific Reports 11(1), 1–13.
  * MacKay (1992) MacKay, D. J. (1992).  A practical Bayesian framework for backpropagation networks.  Neural Computation 4(3), 448–472.
  * Marcus and Rosen (2006) Marcus, M. B. and J. Rosen (2006).  Markov Processes, Gaussian Processes, and Local Times.  Cambridge University Press.
  * Marmin and Filippone (2022) Marmin, S. and M. Filippone (2022).  Deep gaussian processes for calibration of computer models.  Bayesian Analysis 1(1), 1–30.
  * Matheron (1963) Matheron, G. (1963).  Principles of geostatistics.  Econ. Geol. 58(8), 1246–1266.
  * Neal (1996) Neal, R. M. (1996).  Bayesian learning for neural networks, Volume 118.  Springer Science & Business Media.
  * Neidinger (2010) Neidinger, R. D. (2010).  Introduction to automatic differentiation and matlab object-oriented programming.  SIAM review 52(3), 545–563.
  * Pasupathy and Henderson (2011) Pasupathy, R. and S. G. Henderson (2011).  Simopt: A library of simulation optimization problems.  In Proceedings of the 2011 Winter Simulation Conference (WSC), pp. 4075–4085.
  * Plumlee (2014) Plumlee, M. (2014).  Fast prediction of deterministic functions using sparse grid experimental designs.  J. Amer. Statist. Assoc. 109(508), 1581–1591.
  * Plumlee (2021) Plumlee, M. (2021).  Sparse grid designs.  https://www.mathworks.com/matlabcentral/fileexchange/45668-sparse-grid-designs.  MATLAB Central File Exchange.
  * Plumlee and Tuo (2014) Plumlee, M. and R. Tuo (2014).  Building accurate emulators for stochastic simulations via quantile kriging.  Technometrics 56(4), 466–473.
  * Radaideh and Kozlowski (2020) Radaideh, M. I. and T. Kozlowski (2020).  Surrogate modeling of advanced computer simulations using deep Gaussian processes.  Reliability Engineering & System Safety 195, 106731.
  * Rahimi and Recht (2008) Rahimi, A. and B. Recht (2008).  Random features for large-scale kernel machines.  In Advances in Neural Information Processing Systems, pp. 1177–1184.
  * Ritter (2000) Ritter, K. (2000).  Average-case analysis of numerical problems.  Number 1733. Springer Science & Business Media.
  * Sacks et al. (1989) Sacks, J., S. B. Schiller, and W. J. Welch (1989).  Designs for computer experiments.  Technometrics 31(1), 41–47.
  * Salimbeni and Deisenroth (2017) Salimbeni, H. and M. Deisenroth (2017).  Doubly stochastic variational inference for deep Gaussian processes.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Eds.), Advances in Neural Information Processing Systems, Volume 30. Curran Associates, Inc.
  * Sauer et al. (2022a) Sauer, A., A. Cooper, and R. B. Gramacy (2022a).  Vecchia-approximated deep gaussian processes for computer experiments.
  * Sauer et al. (2022b) Sauer, A., A. Cooper, and R. B. Gramacy (2022b).  Vecchia-approximated deep gaussian processes for computer experiments.  Journal of Computational and Graphical Statistics, 1–14.
  * Sauer et al. (2022) Sauer, A., R. B. Gramacy, and D. Higdon (2022).  Active learning for deep gaussian process surrogates.  Technometrics, 1–15.
  * Schmidt et al. (1992) Schmidt, W. F., M. A. Kraaijveld, and R. P. Duin (1992).  Feedforward neural networks with random weights.  In Pattern Recognition, 1992. Vol. II. Conference B: Pattern Recognition Methodology and Systems, Proceedings., 11th IAPR International Conference on, pp. 1–4. IEEE.
  * Sidén and Lindsten (2020) Sidén, P. and F. Lindsten (2020).  Deep Gaussian Markov random fields.  In International Conference on Machine Learning, pp. 8916–8926. PMLR.
  * Stein (2014) Stein, M. L. (2014).  Limitations on low rank approximations for covariance matrices of spatial data.  Spatial Statistics 8, 1–19.
  * Tran et al. (2016) Tran, D., R. Ranganath, and D. M. Blei (2016).  The variational Gaussian process.  Proceedings of the International Conference on Learning Representations (ICLR) 2016.
  * van Dam et al. (2007) van Dam, E. R., B. Husslage, D. den Hertog, and H. Melissen (2007).  Maximin Latin hypercube designs in two dimensions.  Oper. Res. 55(1), 158–169.
  * Vecchia (1988) Vecchia, A. V. (1988).  Estimation and model identification for continuous spatial processes.  Journal of the Royal Statistical Society: Series B (Methodological) 50(2), 297–312.
  * Wang et al. (2020) Wang, W., R. Tuo, and C. F. J. Wu (2020).  On prediction properties of kriging: Uniform error bounds and robustness.  J. Amer. Statist. Assoc. 115(530), 920–930.
  * Yang and Klabjan (2021) Yang, J. and D. Klabjan (2021).  Bayesian active learning for choice models with deep Gaussian processes.  IEEE Transactions on Intelligent Transportation Systems 22(2), 1080–1092.