{'fraction_non_alphanumeric': 0.04994262216696949, 'fraction_numerical': 0.01697427560485799, 'mean_word_length': 4.0942637924735115, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Sparse Expansion For Deep Gaussian Processes

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 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.


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;. 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 , 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 MNIST data set, which is publicly available on the web http://yann.lecun.com/exdb/mnist/, 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.

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.
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 , Vecchia approximation (Vecchia, 1988;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. (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.

Hensman and Lawrence
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.

Deep Gaussian Processes
A GP G is a random function characterized by its mean function µ and covariance function k. To be more specific, given any input x, the output G(x) has a Gaussian distribution with mean µ(x) and variance k(x, x), and given any pair of inputs x and x , the covariance between G(x) and G(x ) is k(x, x ): (1) A W -variate GP is simply a W -vector of GPs [G 1 , G 2 , · · · , G W ]. Without loss of generality, we assume that any multivariate GP in the following content has mean function µ = 0 and independently distributed entries. A H-layer DGP f (H) is then the composition of H multi-variate GPs: where f • g denotes the composition function f (g(·)), where f (h) denotes the output vector from layer h: with h = 1, · · · , H. Evidently, the integral (3) is intractable, leading to the intractability of inference and predictions of DGP f (H) conditioned on any observed data X, Y .

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(X, X) which are O(n 3 ) and 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) ...... k(· , u 1 ) k(· , u 2 ) k(· , u m ) Figure 1: Induced approximationĜ can be represented as a one-layer neural network with correlated Gaussian distributed weights w = k(U, U) −1 G(U) and bias µ.
and random Fourier feature approximation in Cutajar et al. (2017), focus on approximating the covariance matrix k(X, 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 G defined in (1) be a GP with constant mean µ and covariance function k. Inspired by the kriging method (Matheron, 1963;Sacks et al., 1989), we can approximate GP G defined in (1) by the following finite-rank approximation: where we call U = {u i ∈ R d } m i=1 the inducing points and G(U) are the values of G on 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 U are well-chosen. Hence, (4) can be computed in a highly efficient way.
Finite-rank approximation (4) can be represented as a one-layer neural network with correlated Gaussian distributed weights w = k(U, U) −1 G(U) and bias µ 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(U, U), which leads to the following equation:  where R U is the Cholesky decomposition of matrix k(U, U) and Z = [R T U ] −1 G(U) ∈ 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 φ act as an extra hidden layer. Because the approximation (5) consists of only m Gaussians, the time and space complexity for its inference are then reduced to O(m 2 n) and 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 f (H) , all activations in the same layer are characterized by the same covariance function k (h) . Then, where, for h = 1, · · · , 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 U are a specific design called sparse grid, then induced approximations (5) and (6) can be further written as a sparse approximation.

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.

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): For any Markov GP 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: whence Markov GP. Markov GPs have two advantages in computations. Firstly, given any ordered inducing points {u 1 < · · · < u n }, the distribution of G(x) at any x, where u i < x < u i+1 , depends only on its left and right inducing points G(u i ) and 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.
To generalize one-dimensional Markov GP to higher dimensions, we can take the tensor product of d one-dimensional Markov GPs to construct a GP in d dimension: Definition 1. A GP is called TMGP if and only if it is a zero mean GP with covariance function of the form with Markov kernel k j for j = 1, · · · , d. A kernel of the form (8) is called tensor Markov kernel (TMK).
Similar to its one-dimensional counterpart, we want to construct an accurate approximation of TMGP with the smallest possible number of inducing variables. The challenge is to define properly the "left", and "right" neighbors on a multidimensional space, analogously to the real line. This leads us to the hierarchical expansion, described next.

Hierarchical Expansion
To construct the hierarchical expansion, we let inducing variables U be an experimental design X 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. Hierarchical expansion of a TMGP G with inducing points X SG l is then simply: are called hierarchical features. Hierarchical expansion (9) yields a sparse approximation of TMGP that is easy to compute as shown in the following theorems. x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 9 Figure 4: Support of hierarchical feature φ 1 corresponding to x 1 in level 1 (left); supports is also regarded as the index set of φ, because instead of labeling each hierarchical feature by number i, we can also label it by the corresponding inducing variable x i ∈ X SG l . As shown in Figure 5, for any i, x i is the only point at which φ i is not differentiable along all dimensions. Hence, we call x i the center of φ i . From this perspective, we can claim that, given any x * , the distribution of the approximation at x * only depends on 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 φ 1 (left) and φ 2 , φ 3 (right) generated by the tensor Brownian The following numerical example illustrates the sparsity of hierarchical expansion: be the sorted level-2 SG. Then, we numerically compute the inverse Cholesky decomposition R −1 2 and we get Let x * 1 = (0.8147, 0.9058), and x * 2 = (0.2785, 0.5469) be two randomly chosen points. We numerically compute the hierarchical features at x * 1 and x * 2 and we get As we can see, matrix R −1 2 is sparse and each of k(x * 1 , X SG 2 )R −1 2 and k(x * 2 , X SG 2 )R −1 2 has two zero elements. The sparsity is even more visible as the number of inducing points increases.
: Output Entry : Zero Activation : Non-zero Activation

Deep Tensor Markov GPs
A H-layer DTMGP T (H) is defined as the composition of H hierarchical expansions of multi-variate TMGP: where, for h = 1, · · · , H, We then have the following structure: . This structure is also illustrated in Figure 6.
If we treat DTMGP as a DNN, then, for h = 1, · · · , H, layer h is equipped with m (h)

Variational Inference
The proposed DTMGP is mathematically equivalent to a deep neural network with biases 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 θ in the parameterized family F Θ : where Θ denotes the set of parameters and A B denotes the entry-wise multiplication of matrices A and B, i.e., the Hadamard product. In other words, the family F Θ consists of BNNs with activations {φ (h) } H h=1 , parameterized biases, and normally distributed weights.
Given observations (X, Y), our goal is to search for the parameters θ * that maximizes the marginal likelihood P(Y X, θ). The DTMGP f θ * can best interpret how (X, Y) is generated. Computing P(Y|X, θ) 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θ with somẽ to the underlying stochastic process generating (X, Y) and then maximizes the following evidence lower bound (ELBO): which is the lower bound of log-likelihood log Pr Y X, θ .
In our framework, the term E f θ log Pr(Y X, θ) is called the negative energy and can be explicitly written as 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 D KL f θ fθ is the KL-divergence between f θ and the prior fθ and it can be explicitly written The computation of E is efficient at any θ ∈ Θ, the training process of DTMGP is simply defined as 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 (12) can be solved efficiently by using the automatic differentiation technique (Neidinger, 2010;Baydin et al., 2018).

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 x follows an unknown distribution F (·|x): is the realization of the underlying system on inputs X train . We then train competing models using training set (X train , Y train ).
Let Y (·) denote a trained model. We first choose a test set X test from the input space.
For each input x ∈ X test , we sample 100 independent realizations from the true underlying system and from the trained model to get data sets , respectively. We then use the two-sample Kolmogorov-Smirnov (KS) statistic: to quantify the similarity between Y and Y on input x. The smaller D x is, the closer the two distributions are. We access the overall performance of Y (·) via the following averaged two-sample KS statistic over test set X test : We repeat each experiment R times and call each a macro-replication from which we obtain D r , r = 1, . . . , R. Finally, we report the mean and standard deviation 1. GP: standard Gaussian process with squared exponential covariance function. This model can be interpreted as an infinitely-wide BNN with one hidden layer (Neal, 1996). (Schmidt et al., 1992;Graves, 2011;Blundell et al., 2015): ReLU deep neural networks (Glorot et al., 2011)  In all experiments, we use VI to train all the competing deep models and use the standard Bayes rule to get the posterior of the standard GP.

BNN-ReLU
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 in Matlab 2021 with Deep Learning Toolbox on a computer with macOS, 3.3 GHz Intel Core i5 CPU, and 8 GB of RAM (2133Mhz).

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: where B(x) is a Brownian sheet, which is defined as a zero-mean GP with covariance function In this experiment, X train is randomly collected from the cube [0, 1] 2 and Y train is from independently sample paths of B(x) on X train . We investigate the performance of each method as the sizes of the training set (X train , Y train ) increase. For each training set, we randomly select n test = 100 points from [0, 1] 2 as the test set, denoted as X test . We run R = 20 macro-replications and computeD and the associated standard deviationσ for each method.
All the competing models in this experiment share the following architecture:

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 where X = (X 1 , · · · , X 13 ), T (X) is the (random) duration of the longest path from a to i, and f (X) = 13 i=1 X −1 i . The closed form of C(X) is unknown, however, 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, we run m simulation replications with m = 2, 4, · · · , 20. The simulation output is stochastic and has an unknown variance. 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 our BNNs be normal distributed with mean zero 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 X test of size n test = 100 by random samples from [0.5, 5] 13 . We computeD of each method, its associated standard deviationσ based on R = 20 macro-replications. Similar to section 5.1, all the competing models in this experiment share the following architecture: i,j ∼ N (0, 0.04) for any i, j, h. In general, all BNNs can be treated as three-layer deep GPs where the first hidden layer consists of ten GPs and the second hidden layer consists of one GP. Moreover, GPs in the first hidden layer are linear combinations of 391 activation functions and GP in the second hidden layer is linear combination of 241 activation functions.
Such an architecture is flexible enough to capture the random process C(X). is 5000m with m = 2, 4, · · · , 20. Here we exclude the GP models as a competing method because the standard GP regression package cannot handle such a big data set. For all m = 2, 4, · · · , 20, DTMGP significantly outperforms all its competitors in simulating the target stochastic process and its advantage tends to be more prominent as m increases. to the previous artificial random field experiment. The only difference is that DTMGP is more stable compared with BNN-ReLU and DGP-RFF. We believe that it is the sparsity of DTMGP that improves the stability of its training.

Experiment on Real Data
In this section we examine the performance of the proposed methodology on a real data set. We use the MNIST dataset (Deng, 2012) and consider training DGP models on conditional generative problem. The MNIST dataset consists of n = 60000 data pairs where, for any i = 1, · · · , n, input x i ∈ {0, 1, · · · , 9} is a digit from zero to nine and outputs y i is the 28-by-28 image of a hand-written x i . We aim at using the Bayesian formula to train a generative model conditioned on the input x. Given an input x in a specific digit class, the trained model should generate a random matrix output that is similar to data Y in the same class. To summarize, the generative problem can be treated as a stochastic process Y (x) ∈ R 28×28 with x ∈ {0, 1, · · · , 9} and our goal is to use deep GP model to reconstruct Y (·) from observed data.
We compare DTMGP with BNN-ReLU and DGP-RFF introduced in section 5. We omit GP model because the data size is beyond its capacity. All competing deep GP models share the following architecture. We first apply the word embedding technique on input x, which uses a linear transformation A to map discrete input x ∈ {0, 1, · · · , 9} to a 100-dimensional the architecture can be written as: 256 ), j = 1, · · · , 512, Z i,j ∼ N (0, 1), 512 ), l, j = 1, · · · , 28, Z i,l,j ∼ N (0, 1). Figure 11 illustrates the above architecture.
We use all 60000 samples from the MNIST 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 = 0, · · · , 9. Because the output in this experiment is 28-by-28 random matrix and we have no way to sample data from the real distribution of MNIST, the previous KS statistics does not apply. Given that the output random matrix is, in fact, hand-written digit, we can directly check the outputs of each model to evaluate their performances.
Conditional samples generated by DTMGP, BNN-ReLU, and DGP-RFF at the 10 th , 20 th , 30 th , and 40 th are shown in the left of Figure 12. It is obvious that, DTMGP generates BNN-ReLU successfully reconstructs some, but not all, conditional samples. Compared with DTMGP, its performance is unstable. This is because the number of parameters in BNN-ReLU is much higher than that of DTMGP, which leads to an unstable training process. Instead, as shown in the right of Figure 12, the ELBO of DTMGP converges to a stable level after the 100 th iteration (the 5 th epoch) and, hence, the conditional samples generated by DTMGP converge after the 5 th epoch.

Conclusions and Discussion
We use hierarchical features to expand DGPs which are composition of TMGPs and we call this expansion DTMGP. Our expansion is easy to compute and leads to a sparse representation of DGPs. Because the supports of hierarchical features are either nested or disjoint and exhibit a hierarchical structure, only a poly-logarithmic number of activations throughout a DTMGP have non-zero output conditioned on any input or operation. Also, hierarchical features can accurately and efficiently capture local features from inputs. This property of hierarchical features enables DTMGP to have high performance in both prediction problems and generative problems as shown in our numerical experiments. Compared with existing DGP models, inference and training of DTMGP is efficient thanks to its sparsity while instances generated 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.

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 X l with increasing order is simply defined as X l = {2 −l , 2 · 2 −l , · · · , 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 ∈ N, we define the following set consisting only of odd numbers: ρ(l) = {1, 3, 5, · · · , 2 l − 1}.
Then given a level-l dyadic point set X l and a level-(l − 1) dyadic point set X l−1 , we can define the sorted incremental set with D 0 = ∅. It is straightforward to check that For any incremental set D l , we can label any point in D l by x l,i , i ∈ ρ(l) so that the label is unique for any point: x l,i = x l ,i if and only if (l, i) = (l , i ). Now we can define the following dyadic set with the required order: For example, We also define the Cartesian product of one-dimensional sorted X * l 's for later use in Algorithms and proofs. We call this sorted Cartesian product full grid (FG) X * l labeled by l = (l 1 , · · · , l d ) ∈ N d : For simplicity, we define the Cartesian product ρ(l) = × d j=1 ρ(l j ) for any l ∈ N d so that any points in X * l can be represented as x l,i = (x l 1 ,i 1 , x l 2 ,i 2 , · · · , x l d ,i d ) for any l ∈ N d and i ∈ ρ(l). Then we have a more compact representation of X * l : For example, the sorted FG X * (1,2) is Finally, a level-l SG satisfying our required order is the union of sorted FGs as follows: where |l| denote the l 1 norm d j=1 |l j |. To be more specific, a point x l,i ∈ X SG l must be ahead of any x l with |l | > |l|. For example, the two-dimensional sorted SG X SG 2 is ) .
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 X SG l is given by

B Algorithms
We first introduce how to compute R −1 , with R the Cholesky decomposition of k(X * l , X * l ), k is a one-dimensional Markov kernel and 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 X * l . For example, [R −1 ] x,x represents entry with row index corresponding to point x ∈ X * l and column index corresponding to point x ∈ X * l . Now let k = d j=1 k j be a d-dimensional TMK and X * l = × d j=1 X * l j be a sorted FG. Then, R −1 l , with R l the Cholesky decomposition of k(X * l , X * l ) can be directly calculated as where denotes the Kronecker product between matrices and R −1 l j is the inverse Cholesky decomposition yielded by Algorithm 1 with input k j and X * l j . Let R X * l ,X * l represent a sub-matrix of R consisting of entries R x,y , x, y ∈ X * l ⊂ X SG l . Now we can present the algorithm of constructing R −1 l We will prove the correctness of Algorithm 1 and Algorithm 2 in later sections.

C Discussion on Algorithm 1
Before presenting the essential idea of Algorithm 1 in Appendix B, we first need to prove the following lemma and theorem.
Algorithm 1: Computing R −1 for Markov kernel k and point set X * l Input : Markov kernel k, level-l dyadic points X * l Output : R −1 search the closest left neighbor x lef t and right neighbor x right of x ,i in N 5 Solve c 1 , c 2 , and c 3 for the following system: Lemma 2. Let k be a TMK defined on interval U and G be the GP generated by k. For any x < y < z in U , suppose c 1 , c 2 , c 3 satisfies the following conditions: Then function φ(·) = c 1 k(·, x) + c 2 k(·, y) + c 3 k(·, z) is compactly supported on [x, z] and φ, φ k = 1 where ·, · k is the inner product generated by k.
Proof. From the definition of TMK in Lemma 1 in the main paper, the first two conditions Algorithm 2: Computing R −1 l for TMK k and level-l SG X SG l Input : TMK k, level-l SG X SG l Output : can be written as So for any s < x φ(s) = c 1 p(x ∧ s)q(x ∨ s) + c 2 p(y ∧ s)q(y ∨ s) + c 3 p(z ∧ s)q(z ∨ s) = p(s) c 1 q(x) + c 2 q(y) + c 3 q(z) = 0, and for any s > z This proves that the support of φ is [x, z].
Therefore, we have We can notice that equation (14) in Algorithm 1 is exactly the three conditions listed in Lemma 2. The purpose of Algorithm 1 is, in fact, generate compactly supported orthonormal basis functions in the reproducing kernel Hilbert space (Hastie et al., 2009) generated by k as we show in the following theorem: Theorem 3. The hierarchical features where R −1 generated by Algorithm 1 are orthonormal under inner product ·, · k : Proof. The theorem is a consequence of the way we sort the dyadic point set X * l and Lemma 2. Remind that X * l has the following multi-resolution structure: We first prove the following claim: 1. Supports of functions {φ ,i : i ∈ ρ( )} are mutually disjoint for any fixed .
In Algorithm 1, for any pair of points x ,i and x ,i with > , x ,i must be processed ahead of x ,i because the outer for loop will run through all the points in D before processing any point in D . Moreover, the left and right neighbors of any x ,i in the iteration for processing x ,i must be in D and D for some , < .
In iteration ( , i), Lemma 2 tells that the support of φ ,i is [x lef t , x right ] because where c 1 , c 2 , c 3 solve equation (14) in the iteration for processing x ,i . When the outer loop is in the iteration, all points in X * −1 have already been added to N. Remind that if we sort points in X * −1 in increasing order, then X * −1 = {2 −( −1) , 2 · 2 −( −1) , · · · , 1 − 2 −( −1) }. So we have the following relation regarding the distance between x ,i = i2 − and X * −1 : On the other hand, for any other point in D , we have: So both x lef t and x right are must be from X * −1 such that Therefore, the supports of hierarchical features {φ ,i : i ∈ ρ( )} are [i2 − − 2 − , i2 − + 2 − ] : i = 1, 3, · · · , 2 − 1 and they are mutually disjoint. We finish the proof of claim 1.
Disjoint supports of hierarchical features with the same index indicates that they are also mutually orthogonal under ·, · k .
We now prove φ ,i and φ ,i are also orthogonal for any < . Without loss of generality, we can assume that x ,i < x ,i . As a direct consequence of the fact that support of φ ,i is we have the following claim: 2. For any x ,i < x ,i with < and i ∈ ρ( ), i ∈ ρ( ), there are only two possible cases regarding the supports of φ ,i and φ ,i : For case (a) supports of φ ,i and φ ,i are disjoint, the inner product φ ,i , φ ,i k = 0.
To analyze case (b), we first use the following identity: In the above equation, point set X ,i consists of the left and right neighbors of x ,i , and x ,i itself. C ,i = [c 1 , c 2 , c 3 ] T where c 1 , c 2 , c 3 are the solution of equation (14) in the iteration for process φ ,i . X ,i and C ,i are defined in the same manner.
we have: where the last line can be derived directly from the first two conditions in Lemma 2 imposed on C ,i and C ,i .
Lemma 2 and Theorem 3 show that the matrix R generated by Algorithm 1 is indeed the Cholesky decomposition of the covariance matrix: Proof. From our analysis in Theorem 3, we know that in the iteration for processing any x ,i , x lef t and x right are from D and D with some , < . According to how we sort X * l , entries [R −1 ] x lef t ,x ,i , [R −1 ] x ,i ,x ,i , and [R −1 ] x right ,x ,i are in the upper triangular part of R −1 .
In Algorithm 1, operation in iteration for processing x ,i is assigning values to three entries of R −1 , namely, [R −1 ] x lef t ,x ,i , [R −1 ] x ,i ,x ,i , and [R −1 ] x right ,x ,i and let any other value on column [R −1 ] :,x ,i equal 0. Therefore, the returned matrix R −1 is a upper triangular matrix.

D Discussion on Algorithm 2
Algorithm 2 in our main paper is a special case of Plumlee (2014) where, for j = 1, · · · , d, φ j ,i j is the one-dimensional hierarchical feature associated to point x j ,i j in one-dimensional set D j as we defined in Lemma 2 and Theorem 3. Similarly, we can use Kronecker product technique to derive that where the last equality is exactly equation (13).
We now count the time complexity of Appendix B, Algorithm 2. In each iteration, the time complexity for constructing R −1 where the first term on the left-hand side is the time complexity for constructing each component R −1 l j , j = 1, · · · , d, in the Kronecker product and the second term is the time complexity for computing the Kronecker product of upper triangular matrices.
Similar to our counting for its non-zero entries, the total time complexity is then

F Proof of Theorem 2:
Proof. To prove this theorem, we need two more algorithms for computing hierarchical features. The first algorithm is developed from Appendix B-Algorithm 1: We first count the number of non-zero entries on the output of Algorithm 4. Given any x * , we first fix a in the outer for-loop. From Lemma 2, we can derive that the only non-zero [φ] x ,i is the ( , i) satisfying the requirement x lef t < x * < x right .
Because of the way we sort X * l = [D 1 D 2 · · · D l ], all the points {x ,i , i ∈ ρ( )} in the inner for-loop are in D and none of their closet neighbor is in D . As a result, for any fixed , there is only one ( , i) satisfies x lef t < x * < x right in the inner-loop and, hence, there is only one non-zero entry is added to φ in Algorithm 4 in each otter iteration. This shows that the number of non-zero on the output of Algorithm 4 is l.
We then propose an algorithm to compute φ(x * ) = R −T l k(X SG l , x * ) as shown in Algorithm 5. Algorithm 5 is a simple extension of Appendix B-Algorithm 2. Based on this algorithm, we can estimate the number of non-zero entries on φ(x * ).
We first estimate the number of non-zero entries for φ l (x * ) in each iteration. Because the number of non-zero entries on each [R T ] −1 l j k j (X * l j , x * j ) is l j , the total number of non-zero entries on tensor φ l (x * ) is O( d j=1 l j ).
of G(y) − k(y, X SG l )K −1 G(X SG l ) for any pair of x, y not sharing the same left and right neighboring points in X SG l . Now for any i and x ∈ (i−1)2 −l , i2 −l , the posterior M i (x) := G(x)−k(x, X i )K −1 i G(X i ) can be treated again as a Markov GP with fixed boundary at x i−1 and x i . Because M i is a Markov GP defined on interval of length 2 −l with fixed boundary condition, it is straightforward to derive that max x E |M i (x)| 2 = O(2 −l ). This holds for any i so we can get when dimension D = 1, max x∈(0,1) where I i denote the inverval (i − 1)2 −l , i2 −l . Now suppose our theorem holds for D = d − 1. We now prove the case D = d. Let d j=1 k j denote TMK k(x, x). For any x, the L 2 distance has the following identity where M i is any m-by-m matrix and v i is any m-by-1 vector for i = 1, 2.
Let ∆ (l) j denotek (l) j −k (l−1) j and l −d denote the vector (l 1 , · · · , l d−1 ). We have the following identity: where the last line is from telescoping sum. Hence, the L 2 distance becomes Proof of Theorem 6. Let G (h) denote the TMGP in layer h of f (H) and letG (h) denote the hierarchical expansion of G (h) . We prove the statement by induction in the number of layers H. When H = 1, the statement is obviously true because it is equivalent to Theorem 5.
Then, we have the following estimate of (s.8): (s.11) where C and C are some number independent of {m (h) } H−1 h=1 , the second line is from the Hölder condition (s.10), the third line is from the induction assumption.
Putting equations (s.9) and (s.11) together, we can have the final result.
From theorem (6), we can see that if a DTMGP is not tremendously deep and the number of inducing points dominates the dimension of output from each layer, then the distance between f (H) and T (H) is relatively close. DGP f (H) is, in fact, an infinite-dimensional non-parametric model while DTMGP T (H) is a parametrized model with finitely many parameters. Therefore, theorem 6 partially reflects that DTMGP can also well approximate other complex stochastic systems.