Title: DEEP RECURRENT GAUSSIAN PROCESS WITH VARIA-TIONAL SPARSE SPECTRUM APPROXIMATION

Abstract: Modeling sequential data has become more and more important in practice. Some applications are autonomous driving, virtual sensors and weather forecasting. To model such systems, so called recurrent models are frequently used. In this paper we introduce several new Deep Recurrent Gaussian Process (DRGP) models based on the Sparse Spectrum Gaussian Process (SSGP) and the improved version, called Variational Sparse Spectrum Gaussian Process (VSSGP). We follow the recurrent structure given by an existing DRGP based on a specific variational sparse Nyström approximation, the recurrent Gaussian Process (RGP). Similar to previous work, we also variationally integrate out the input-space and hence can propagate uncertainty through the Gaussian Process (GP) layers. Our approach can deal with a larger class of covariance functions than the RGP, because its spectral nature allows variational integration in all stationary cases. Furthermore, we combine the (Variational) Sparse Spectrum ((V)SS) approximations with a well known inducing-input regularization framework. For the DRGP extension of these combined approximations and the simple (V)SS approximations an optimal variational distribution exists. We improve over current state of the art methods in prediction accuracy for experimental data-sets used for their evaluation and introduce a new data-set for engine control, named Emission.

1. INTRODUCTION:

Modeling sequential data for simulation tasks in the context of machine learning is hard for several reasons. Their internal structure poses the problem of modeling short term behavior and long term behavior together for different types of data variables, where the data variables themselves might differ in the information gain in the chosen time frequency. Recurrent models (Hochreiter & Schmidhuber, 1997; Nelles, 2013; Pascanu et al., 2013) have proven to perform well on these tasks. They consist of output-data and input-data structured sequentially for shifted discrete time steps. The general form of a recurrent model is given by

EQUATION

where x i is an external input, y i is an output observation, h i is a latent hidden representation or state at time i = H + 1, . . . , N , where N ∈ N is the number of data samples, H ∈ N is the chosen time horizon, ζ, ψ are non-linear functions modeling transition and observation and h i , y i are transition and observation noise, which are adjusted for the specific problem (details on dimensions and ranges will be specified in upcoming sections). In control and dynamical system identification previous work on Bayesian recurrent approaches for modeling sequential data usually makes use of (non-)linear auto-regressive with exogenous inputs models ((N)ARX) and state-space models (SSM), for both see Nelles (2013) . The general recurrent model given in Equation ( 1) and (2) represents both cases. This can be recognized by its general recurrent and hierarchical structure. This work deals with deep learning in a recurrent fashion for modeling sequential data in a Bayesian non-parametric approach by using GPs. To make a connection to the general recurrent model, the deep structure arises by defining ζ in Equation (1) in a deep manner (Pascanu et al., 2013, Section 3) . To achieve scalability, GPs normally make use of sparse approximations for the covariance function.

This paper proposes DRGP models based on (V)SS approximations, denoted by DRGP-(V)SS. For reproducibility of the experimental results, we provide the code online1 . Therefore, we follow the same deep recurrent structure as introduced in Mattos et al. (2016) . To summarize, the contributions of this paper are the following:

• Extension of the sparse GP based on the SS approximation introduced by Lázaro-Gredilla et al. (2010) and the improved VSS approximation by Gal & Turner (2015) to DRGPs; • Improvement of regularization properties of the variational bounds through the combination of the (V)SS approximations with the inducing-point (IP) regularization of Titsias & Lawrence (2010); • Propagation of uncertainty through the hidden layers of our DGPs by variationally integrating out the input-space; • The existence of an optimal variational distribution in the sense of a functional local optimum of the variational bounds of our DRGPs models is established.

The DRGP of Mattos et al. ( 2016) is limited to a small class of deterministic covariance functions, because the covariance functions variational expectation has to be analytically tractable. Using the (V)SS approximations instead, we can derive a valid approximation for every stationary covariance function, because the basis functions expectation is always tractable. We show that this approach improves over state of the art approaches in prediction accuracy on several cases of the experimental data-sets used in Mattos et al. (2016) ; Svensson et al. (2016) ; Al-Shedivat et al. (2017) ; Salimbeni & Deisenroth (2017); Dörr et al. (2018) in a simulation setting. For scalability, Distributed Variational Inference (DVI) (Gal et al., 2014) is recommended and can lower the complexity from O(N M2 Q max (L + 1)) down to O(M 3 ) for N ≤ M , M the sparsity parameter, (L + 1) the amount of GPs and Q max is the maximum over all input dimensions used in our defined deep structure for ζ and ψ. Therefore, the number of cores must scale suitably with the number of training-data.

2. RELATED WORK TO GAUSSIAN PROCESSES WITH SSM AND DGPS:

An SSM with GPs (GP-SSM) for transition and observation functions is used by Wang et al. (2005) , where the uncertainty in the latent states is not accounted for, which can lead to overconfidence. 

y i = f (L+1) (f (L) (f (L-1) (. . . (f (1) (x i ) + h (1) i ) . . . ) + h (L-1) i ) + h (L) i ) + y i ,

where the index i = 1, . . . , N is not necessarily the time and where we define h

(1)

i def = f (1) (x i ) + h (1) i , h (l+1) i def = f (l) (h (l) i ) + h (l) i

, for l = 2 . . . , L -1, where L ∈ N is the number of hidden layers. The noise h (l) i , y

i is assumed Gaussian and the functions f (l) are modeled with GPs for l = 1, . . . , L + 1. To obtain computational tractability, in most cases variational approximation and inference is used. Damianou & Lawrence (2013) introduced these kind of DGPs based on the sparse variational approximation following Titsias (2009); Titsias & Lawrence (2010). Based on this, Dai et al. (2016) introduced a DGP with a variationally auto-encoded inference mechanism and which scales on larger data-sets. Cutajar et al. (2016) introduced a DGP for the so called Random Fourier Features (RFF) approach (Rahimi & Recht, 2008) , where the variational weights for each GPs are optimized along with the hyperparameters. This approach does not variationally integrate out the latent inputs to carry through the uncertainty and no existence of an optimal variational distribution for the weights is proven to reduce the amount of parameters to optimize in training. Furthermore, Salimbeni & Deisenroth (2017) introduced an approximation framework for DGPs, which is similar to the single GP of Hensman & Lawrence (2014) , but does not force independence between the GP layers and which scales to billions of data. Two state of the art approaches for DRGPs have been introduced by Mattos et al. (2016) , the RGP, which we call DRGP-Nyström, based on the variational sparse Nyström/inducing-point approximation introduced by Titsias (2009); Titsias & Lawrence (2010), as well as Al-Shedivat et al. ( 2017), which we call GP-LSTM, based on deep kernels via a Long-short term memory (LSTM) network (Hochreiter & Schmidhuber, 1997) , a special type of Recurrent Neural Network (RNN). DRGP-Nyström uses a recurrent construction, where the auto-regressive structure is not realized directly with the observed output-data, but with the GPs latent output-data and uses a Variational Inference (VI) framework, named Recurrent Variational Bayes (REVARB). The structure acts like a standard RNN, where every parametric layer is a GP. So additionally uncertainty information can be carried through the hidden layers. GP-LSTM is a combination of GPs and LSTMs. LSTMs have proven to perform well on modeling sequential data. LSTMs try to overcome vanishing gradients by placing a memory cell into each hidden unit. GP-LSTM uses special update rules for the hidden representations and the hyperparameters through a semi-stochastic optimization scheme. It combines a GP with the advantages of LSTMs by defining structured recurrent deep covariance functions, also called deep kernels, which reduces the time and memory complexities of the linear algebraic operations (Wilson et al., 2016) .

3. GAUSSIAN PROCESSES AND VARIATIONAL SPARSE SPECTRUM GP:

Loosely speaking, a GP can be seen as a Gaussian distribution over functions. We will first introduce GPs and GP regression and then recall the SSGP by (Lázaro-Gredilla et al., 2010) and its improved version VSSGP by (Gal & Turner, 2015) . Based on these, we derive new variational approximations. We use the notation a, f x , y, (italic) for random variables, a, f(x), y (upright) for realizations/samples and data.

3.1. GAUSSIAN PROCESSES:

A stochastic process f is a GP if and only if any finite collection of random variables f

X def = [f x1 , . . . , f x N ]

T forms a Gaussian random vector (Rasmussen, 2006) . A GP is completely defined by its mean function m : R Q → R, x → m(x), Q ∈ N the input-dimension, and covariance function (Kallenberg, 2006, Lemma 11.1) , where

k : R Q × R Q → R, (x, x ) → k(x, x )

m(x) def = E [f x ] , k(x, x ) def = cov(f x , f x ) = E [(f x -m(x))(f x -m(x ))] ,

and the GP will be written as f ∼ GP(m, k). Be aware of that a valid covariance function must produce a positive definite matrix

K N N def = (k(x i , x j )) N i,j=1 ∈ R N ×N , when filling in combinations of data-input points x i , i = 1, . . . , N . Let y def = [y 1 , . . . , y N ] T ∈ R N , X def = [x 1 , . . . , x N ] T ∈ R N ×Q

and we assume y i = f(x i ) + y i , where y i ∼ N (0, σ 2 noise ), for i = 1, . . . , N , and our aim is to model any set of function values f def = [f(x 1 ), . . . , f(x N )] T ∈ R N at X as samples from a random vector f X . Moreover, we assume the prior f X |X ∼ N (0, K N N ), meaning that any set of function values f given X are jointly Gaussian distributed with mean 0 ∈ R N and a covariance matrix K N N . The predictive distribution p f x * |x * ,X,y for a test point x * ∈ R Q , where

K N * def = (k(x i , x * )) N i=1 ∈ R N ,

and analogously K * * , K * N , can be derived through the joint probability model and conditioning as

f x * |x * , X, y ∼ N (K * N (K N N + σ 2 noise I N ) -1 y, K * * -K * N (K N N + σ 2 noise I N ) -1 K N * ).

In preview of our experiments in Section 5 and the following sections, we choose a specific covariance function, the spectral mixture (SM) covariance function (Wilson & Adams, 2013 )

EQUATION

with an amplitude σ 2 power and length scales l q , p q ∈ R, q = 1, . . . , Q. As p q → ∞, this corresponds to the squared exponential (SE) covariance function in the limit (Gal & Turner, 2015) .

3.2. VARIATIONAL SPARSE SPECTRUM GP:

We introduce the SSGP following Gal & Turner (2015) . For a stationary covariance function k on

R Q × R Q there exists a function ρ : R Q → R, τ → ρ(τ ), such that k(x, x ) = ρ(x -x ) for x, x ∈ R Q .

Bochner's theorem states that any stationary covariance function k can be represented as the Fourier transform of a positive finite measure µ (Stein, 2012). Then ρ(τ ), using τ = xx , can be expressed via Monte Carlo Approximation (MCA) following Gal & Turner (2015) , Section 2, as

EQUATION

We refer to z m as the spectral points, b m as the spectral phases and u m as the pseudo-input/inducing points for m = 1, . . . , M . Choosing the probability density like in Gal & Turner (2015) , Proposition 2, we approximate the SM covariance function with a scaling matrix

L def = diag([2πl q ] Q q=1 ) ∈ R Q×Q , a scaling vector p def = [p -1 1 , . . . , p -1 Q ] T ∈ R Q and Z def = [z 1 , . . . , z M ] T ∈ R Q×M via φ(x, Z) def = 2σ 2 power M -1 cos(2π(L -1 z 1 + p) T (x -u 1 ) + b 1 ), . . . , cos(2π(L -1 z M + p) T (x -u M ) + b M ) T ∈ R M .

(5)

Here we sample b ∼ Unif [0, 2π], z ∼ N (0, I Q ) and we set K(SM)

N N def = ΦΦ T with Φ def = [φ(x 1 , Z), . . . , φ(x N , Z)]

T ∈ R N ×M . In Gal & Turner (2015) the SSGP was improved to VSSGP by variationally integrating out the spectral points and instead of optimizing the spectral points, additionally optimizing the variational parameters. We follow the scheme of Gal & Turner (2015) , Section 4, for the 1-dimensional output case y ∈ R N ×1 . By replacing the covariance function with the sparse covariance function k(SM) and setting the priors to

EQUATION

for m = 1, . . . , M , where we have y|a, Z, U , X ∼ N (Φa,

EQUATION

highlighting U just in the integral, to be notationally conform to Gal & Turner (2015) , Section 3. Now, to improve the SSGP to VSSGP, variational distributions are introduced in terms of

q Z def = M m=1 q zm , where z m ∼ N (α m , β m ), q a , where a ∼ N (m, s),

with β m ∈ R Q×Q diagonal, for m = 1, . . . , M , and s ∈ R M ×M diagonal. From here on we use variational mean-field approximation to derive the approximate models with different lower bounds to the log ML introduced by Gal & Turner (2015) :

log(p(y|X)) ≥ E qaq Z [log(p(y|a, Z, U, X))] -KL(q a ||p a ) -KL(q Z ||p Z )).

(8) As usual, KL defines the Kullback-Leibler (KL) divergence. By proving the existence of an optimal distribution q opt a for a Gal & Turner (2015) , Proposition 3, in the sense of a functional local optimum of the right hand side of ( 8), where a ∼

N (A -1 Ψ T 1 y, σ 2 noise A -1 ), with A = Ψ 2 + σ 2 noise I M , Ψ 1 = E q Z [Φ] ∈ R N ×M , Ψ 2 = E q Z Φ T Φ ∈ R M ×M ,

we can derive the optimal bound case.

3.3. (V)SSGP WITH REGULARIZATION PROPERTIES VIA INDUCING POINTS (IP):

As a first contribution of this paper we combine two approximation schemes to two new methods (V)SSGP-IP. We want to point out that the (V)SSGP does not have the same regularization properties as the GP of Titsias & Lawrence (2010), when optimizing the parameters U, because the prior in (6) of the weights a is defined generically via Bishop (2006) , Equations (2.113) -(2.115 ). These parameters U, following Gal & Turner (2015) , are similar to the sparse pseudo-input approximation Titsias ( 2009), but in the lower bound in (8) they are simply used without being linked to the weights a. We now define them as

EQUATION

where K N N , K N M , K M M are defined through the given covariance function in Equation ( 3). We can show that for these definitions the integral in Equation ( 7) can be marginalized straightforward for the weights a. We then obtain that our data-samples y are coming from a Gaussian distribution with mean 0 and the true covariance matrix K N N , plus the discrepancy of the two wellknown sparse covariance matrices for the Sparse Spectrum (not exactly, as we have ΦK -1 M M Φ T ) and the Nyström case, plus the noise assumption. This expression can not be calculated efficiently, but shows that we obtain a GP approximation, which can be seen as a trade-off between these two sparse approximations. Following Titsias & Lawrence (2010), Section 3.1, the optimal variational distribution for a collapses by reversing Jensen's inequality and is similar to the one obtained in Titsias & Lawrence (2010). The resulting bounds (SSGP-IP without KL(q Z ||p Z ) and Ψ 1 = Φ, Ψ 2 = Φ T Φ) can be calculated in the same way as Titsias & Lawrence (2010) until Equation ( 14) in closed form:

EQUATION

Consequently, the resulting bound in (11) has an extra regularization property compared to the right hand side of (8), which is reflected in the different form of A = Ψ 2 + σ 2 noise K M M , which involves K M M , the chosen covariance matrix filled in with the pseudo-input points U, and three extra terms

log(|K M M |) 2 , tr(K N N ) 2σ 2 noise , tr(K -1 M M K M N K N M ) 2σ 2 noise .

3.4. VARIATIONAL APPROXIMATION OF THE INPUT-SPACE FOR (V)SSGP(-IP):

As a second contribution of this paper we marginalize also the input-space. This is not straightforward, as it is not clear whether for the (V)SS covariance function in Equation ( 4) and ( 5) these expressions even exist. To prevent misunderstanding, we will write from now on

H = [h 1 , . . . , h N ] T ∈ R N ×Q instead of X,

to highlight that H is now a set of latent variables (later on the hidden states of the DGP). Therefore, we introduce priors and variational distributions

p H def = N i=1 p hi , where h i ∼ N (0, I Q ), q H def = N i=1 q hi , where h i ∼ N (µ i , λ i ), λ i ∈ R Q×Q diagonal, for i = 1, . . . , N .

As a consequence, for VSSGP-IP we overall derive statistics

Ψ 0 = tr (E q H [K N N ]) = N σ 2 power , Ψ 1 = E q Z q H [Φ] ∈ R N ×M , Ψ 2 = E q Z q H Φ T Φ ∈ R M ×M and Ψ reg = E q H [K M N K N M ]

as defined in the Appendix 6.3. These statistics are essentially the given matrices Φ, Φ T Φ, K M N K N M from the beginning, but every input h i and every spectral point z m is now replaced by a mean µ i , α m and a variance λ i , β m resulting in matrices of the same size. The property of positive definiteness is preserved. The SSGP-IP model derives by being not variational over the spectral points. This extra step allows to propagate uncertainty between the hidden layers of a DGP, as we gain an extra variance parameter for the inputs. For the IP case we get the lower bound:

EQUATION

For the extension of the SSGP and VSSGP in the optimal bound case in (8) (we only focus on this case in the following), we again have

A = Ψ 2 + σ 2 noise I M and eliminate log(|K M M |) 2 , Ψ0 2σ 2 noise , tr(K -1 M M Ψreg) 2σ 2

noise in the lower bound (12).

4. DRGP WITH VARIATIONAL SPARSE SPECTRUM APPROXIMATION:

In this section we want to combine our newly derived (V)SS-(IP) approximations in the Sections 3. 

EQUATION

N N ), l = 1, . . . , L

EQUATION

N N ), l = L + 1, with ζ, ψ in Equation ( 1) and (2), h (l) i ∼ N (0, (σ

(l) noise ) 2 ), y i ∼ N (0, (σ (L+1) noise ) 2 ) and N = N -H x , for i = H x + 1, . . . , N . The matrix K (l)

N N represents a covariance matrix coming from our chosen k in Equation ( 4) and (5). A set of input-data Ĥ(l) = [ ĥ(l) Hx+1 , . . . , ĥ(l) N ] T is specified as ĥ(l)

EQUATION

i-H h , [x i-1 , . . . ,

EQUATION

i-H h , h

EQUATION

where ĥ(1)

i ∈ R H h +HxQ , ĥ(l) i ∈ R 2H h for l = 2, . . . , L, ĥ(L+1) i ∈ R H h , for i = H x + 1, . . . , N .

For simplification we set H def = H x = H h in our experiments. Now we use the new approximations in the Sections 3.3, 3.4, to derive first, for the setting (V)SSGP, the new joint probability density

EQUATION

H h ] T ∈ R 2H h -Hx . Here the priors are

p a (l) , where a (l) ∼ N (0, I M ), p Z (l) def = M m=1 p z (l) m , where z (l) m ∼ N (0, I Q ), p h(l) , where h(l) ∼ N (0, I 2H h -Hx ),

the product of them is defined as P REVARB , and the variational distributions are q a (l) , where a (l) ∼ N (m (l) , s (l) ), q Z (l)

def = M m=1 q z (l) m , where z (l) m ∼ N (α (l) m , β (l) m ), q h (l) def = N i=1+H h -Hx q h (l) i

, where h

(l) i ∼ N (µ (l) i , λ (l) i ),

the product of them is defined as Q REVARB , and where

β (l) m ∈ R Q×Q is diagonal, for i = 1 + H x -H h , . . . , N , m = 1 . . . , M , l = 1, . . . , L + 1.

For the setting (V)SSGP-IP we choose no variational distribution for a (l) , but, similar to the assumptions in ( 9), (10) for l = 1, . . . , L + 1, the prior assumptions

p a (l) |U (l) , where a (l) |U (l) ∼ N (0, K (l) M M ), p h (l)

Hx+1: |a (l) ,Z (l) , Ĥ(l) ,U (l) , where h (l)

Hx+1: |a (l) , Z (l) , Ĥ(l) , U (l) ∼ N (Φ (l) (K (l) M M ) -1 a (l) , K (l) N N -K (l) N M (K (l) M M ) -1 K (l) M N + (σ (l)

noise ) 2 I N ). This defines our models for the cases DRGP-VSS, DRGP-VSS-IP. In the case, where we are not variational over the spectral-points, we derive the simplified versions DRGP-SS, DRGP-SS-IP.

4.2. DRGP-(V)SS(-IP) VARIATIONAL EVIDENCE LOWER BOUND (ELBO):

Using standard variational approximation techniques (Blei et al., 2017) , the Recurrent Variational Bayes lower bound for the (V)SS approximations, denoted as REVARB-(V)SS, is given by

EQUATION

and for the (V)SS-IP approximations, denoted as REVARB-(V)SS-IP, is given by

EQUATION

where •, • means, taking the expectation under the integral.

Additionally, for the simple (V)SS approximations in 14, the optimal bound L REVARB (V)SS-opt can be obtained immediately, analogously to Gal & Turner (2015) , Proposition 3, and by the fact that the bound decomposes into a sum of independent terms for a (l) . Maximizing the lower bounds is equivalent to minimizing the KL-divergence of Q REVARB and the true posterior. Therefore, this is a way to optimize the approximated model parameter distribution with respect to the intractable, true model parameter posterior. Calculating L REVARB (V)SS-opt/IP requires O(N M 2 Q max (L + 1)), where

Q max = max l=1...,L+1 Q (l) , Q (l) def

= dim( ĥ(l) i ) and ĥ(l) i from Equation (13). DVI can reduce the complexity to O(M 3 ) if the number of cores scales suitably with the number of training-data, see Appendix 6.3.6 for a detailed description. A detailed derivation of the REVARB-(V)SS(-IP) lower bounds can be found in the Appendix 6.3.4.

4.3. MODEL PREDICTIONS:

After model optimization based on the lower bounds in the Equation ( 14) and ( 15), model predictions for new ĥ(l) * in the REVARB-(V)SS(-IP) framework can be obtained based on the approximate, variational posterior distribution Q REVARB . They are performed iteratively with approximate uncertainty propagation between each GP layer. We derive q h (l) * from previous time-steps and it is per definition Gaussian with mean and variance derived from previous predictions for l = 1, . . . , L. These kind of models propagate the uncertainty of the hidden GP layers' outputs, not of the observed output-data and are relevant for good model predictions. The detailed expressions for the mean and variance of the predictive distribution involved during the predictions can be found in the Appendix 6.3.5. 2010) is named GP-SS, the VSSGP of Gal & Turner (2015) is named GP-VSS. The setting in this system identification task is simulation. This means that, together with past exogenous inputs, no past measured output observations (but perhaps predicted output observations) are used to predict next output observations. To enable a fair comparison, all methods are given access to a predefined amount of data. Details about the methods, their configuration, as well as the benchmark data-sets can be found in the Appendix 6.1 and 6.2.

5.1. IMPLEMENTATION:

Our methods DRGP-(V)SS(-IP) were implemented in Python, using the library Theano, and in Matlab R2016b. For the optimization/training we used Python, Theano. Theano allows us to take full advantage of the automatic differentiation to calculate the gradients. For simulation and visualization we used Matlab R2016b. We further implemented in Matlab R2016b the methods DRGP-Nyström, GP-SS, GP-DTC, GP-FITC, GP-full and used these implementations for the experiments. For GP-VSS, GP-LSTM, LSTM, RNN, DGP-RFF and DGP-DS we used the published code 3456 . For GP-SSM the published code is only applicable for small range of state dimension and a small time horizon, so we just show the results from their paper. For PR-SSM we also show the results from their paper. 

5.2. BENCHMARK DATA-SETS MODEL LEARNING AND COMPARISON:

In Figure 1 we show a comparison of the latent model states before and after training, the simulation results, as well as the simulated latent states for two data-sets, Drive and Actuator, for the model DRGP-VSS. We initialize the states with the output training-data for all layers with minor noise (first column) and after training we obtain a trained state (second column). Unlike Mattos et al. ( 2016) Figure 2 , (l), we get good variance predictions for all data-sets. We used our own implementation for Mattos et al. ( 2016), which gave the same good results as for our methods. Therefore, we think it is an implementation issue. The test RMSE results for all methods are summarized in Table 1 .

The results show, that on most data-sets DRGP-(V)SS(-IP) improve slightly in comparison to other methods. In order to evaluate the reproducing quality of our results, we provide a robustness test for our methods and DRGP-Nyström on the data-sets Drive and Damper in Figure 2 , 3. We run the optimization for different time horizons. For every method we visualized a boxplot with whiskers from minimum to maximum with 10 independent runs. For our models we obtain good results compared with DRGP-Nyström on these data-sets, in particular for the setting of time horizons of Table 1 with H = 10. In Figure 4 

5.3. CONCLUSION:

In this paper we introduced four new DRGPs based on the SS approximation introduced by Lázaro-Gredilla et al. ( 2010) and the improved VSS approximation by Gal & Turner (2015) . We combined the (V)SS approximations with the variational inducing-point approximation from Titsias & Lawrence ( 2010), also integrated variationally over the input-space and established the existence of an optimal variational distribution for a (l) in the sense of a functional local optimum of the variational bounds REVARB-(V)SS-(IP). We could show that our methods slightly improve on the data-sets used in this paper compared to the RGP from Mattos et al. ( 2016) and other state-of-the-art methods, where moreover our sparse approximations are also practical for dimensionality reduction as shown in Titsias & Lawrence (2010) and can be further expanded to a deep version in this application (Damianou & Lawrence, 2013). Furthermore, Hoang et al. ( 2017) introduced a generalized version of the (V)SS approximation, which should be adaptable for our case.

6. APPENDIX:

This additional material provides details about the derivations and configurations of the proposed DRGP-(V)SS(-IP) in the Sections 3-4. and elaborates on the methods and the employed data-sets in the experiments in Section 5.

6.1 DATA-SET DESCRIPTION In this section we introduce the data-sets we used in our experiments. We chose a large number of data-sets in training size going from 250 to 12500 data-points in order to show the performance for a wide range. We will begin with the smallest, the Drive data-set, which was first introduced by Wigren (2010). It is based on a system which has two electric motors that drive a pulley using a flexible belt. The input is the sum of voltages applied to the motors and the output is the speed of the belt. The data-set Dryer7 describes a system where air is fanned through a tube and heated at an inlet. The input is the voltage over the heating device (a mesh of resistor wires). The output is the air temperature measured by a thermocouple. The third data-set Ballbeam 78 describes a system where the input is the angle of a beam and the output the position of a ball. Actuator is the name of the fourth data-set, which was described by Sjöberg et al. (1995) and which stems from an hydraulic actuator that controls a robot arm, where the input is the size of the actuator's valve opening and the output is its oil pressure. The Damper data-set, introduced by Wang et al. (2009) , poses the problem of modeling the input-output behavior of a magneto-rheological fluid damper and is also used as a case study in the System Identification Toolbox of Mathworks Matlab. The data-set Power Load9 , used in Al-Shedivat et al. (2017) , consists of data where the power load should be predicted from the historical temperature data. This data-set was used for 1-step ahead prediction, where past measured output observations are used to predict current or next output observations, but we will use it here for free simulation. We down-sampled by starting with the first sample and choosing every 4th datapoint, because the original data-set with a size of 38064 samples and a chosen time-horizon of 48 is too large for our implementation, which is not parallelized so far. The newly provided data-set Emission10 contains an emission-level of nitrogen oxide from a driving car as output and as inputs the indicated torque, boost pressure, EGR (exhaust gas recirculation) rate, injection, rail pressure and speed. The numerical characteristics of all data-sets are summarized in The methods for the system identification tasks and their references are summarized in Table 2 . For LSTM, RNN we chose the same setting for the amount of hidden layers and time horizon as for our DRGPs. We used Adam optimizer with a learning rate of 0.01. As activation function we chose tanh and NARX structure H x = H y for the first layer. We tested with 8, 16, 32, 64, 128 hidden units (every hidden layer of a RNN is specified by a hidden unit parameter) for all training data-sets.

For GP-LSTM we chose the same setting for the amount of hidden layers, pseudo-input points and time horizon as for our DRGPs. We used Adam optimizer with a learning rate of 0.01. As activation function we chose tanh and NARX structure H x = H y for the first layer. We tested with 8, 16, 32, 64, 128 hidden units (every hidden layer of a RNN is specified by a hidden unit parameter) for all training data-sets. For the data-sets with training size smaller or equal to 2000 we chose the version GP-LSTM in Al-Shedivat et al. ( 2017) and for the ones larger than 2000 the scalable version MSGP-LSTM. We did not pre-train the weights.

For DGP-RFF we tested for L = 1, . . . , 3 number of GPs. We used Adam optimizer with a learning rate of 0.01. The batch size was chosen to be the training data-size and the dimensions of the GP layers to be 5. We set the flags q Omega fixed=1000 and learn Omega fixed=var fixed, so fixing the spectral-point parameters, and iterated until 2000 just for this GP.

For DGP-DS we tested for L = 1, . . . , 3 number of GPs. We used natural gradients for the last GP with gamma 0.1 and Adam optimizer for the others with a learning rate of 0.01. The batch size was chosen to be the training data-size and the dimensions of the GP layers to be 5. All GPs which use the Nyström approximation were initialized for the pseudo-inputs points with a random subset of size M from the input-data and trained with SE covariance function. For the ones which use the (V)SS approximations, which includes our methods, we trained with a spectral-point initialization sampled from N (0, I Q (l) ), an initialization for the pseudo-input points with a random subset of size M from the input-data (or setting them all to zero; not in the IP case). For our methods DRGP-(V)SS(-IP) and GP-VSS we fixed the length scales p (l) q = ∞, for all q, l. So all GPs with

6.3. VARIATIONAL APPROXIMATION AND INFERENCE FOR DRGP-(V)SS(-IP):

In the Sections 6.3.1-6.3.3 we derive the statistics Ψ 0 , Ψ 1 , Ψ 2 and Ψ reg for the four model cases DRGP-(V)SS(-IP). In Section 6.3.4 we derive the resulting variational lower bounds. In Section 6.3.5 we show the mean and variance expressions of the predictive distributions of our DRGP models.

In the following we use the abbreviations and formulas

• B.1, (Gal & Turner, 2015 , see Section 4.1),

• B.2, (Rasmussen, 2006, see A.7.),

• 1 2 (cos(a -b + x -y) + cos(a + b + x + y)) = cos(x + a) cos(y + b).

• JI, for Jensen's inequality.

6.3.1 DRGP-VSS(-IP), THE STATISTICS Ψ 1 , Ψ 2

For the versions DRGP-VSS(-IP) the statistics are In this section we derive the different variational lower bounds for our models DRGP-(V)SS(-IP). We first show the bound L REVARB (V)SS without optimal variational distribution for a (l) . Then the bounds L REVARB (V)SS-opt with optimal variational distribution for a (l) follows, as well as the L REVARB (V)SS-IP case. We use the simplified notation dAdZdH = da (1) . . . da (L+1) dZ (1) . . . dZ (L+1) dh (1) . . . dh (L) . q(a (l) )q(Z (l) )q(h (l) )p(h (l) Hx+1: |a (l) , Z (l) , Ĥ(l) , U (l) )da (l) dZ (l) dh (l)

EQUATION

EQUATION

(Ψ 2 ) n mm = E qz m q hn Φ T nm Φ nm B.1 = E q hn 2σ 2 power (M ) -1 ( 1 2 + 1 2 e -2 hT nm β m hnm cos(2( αT m (h n -u m ) + b m )) B.2 = σ 2 power (M ) -1 1 + Σ2 m Znm e -

- L l=1     N 2 - N i=1+Hx-H h log(2πλ (l) i ) 2 + H h i=1+Hx-H h log(2π) 2 + λ (l) i + µ (l) i 2 2     - L+1 l=1

KL(q a (l) ||p a (l) ) -KL(q Z (l) ||p Z (l) . Using the optimal distribution for a (l) ∼ N ((A (l) ) -1 (Ψ where for the optimal distribution case for a (l) we have with A (l) = Ψ Hx+1: with y Hx+1: . In the DRGP-(V)SS-IP case we make predictions for each layer l and new ĥ(l) * through Hx+1: with y Hx+1: . Prediction for a new input in all cases has time complexity O((L + 1)M 3 ), which comes from the iterative prediction through all GP layers and the calculation of the statistics, see Appendix 6.3.6.

EQUATION

:

DRGP-(V)SS code available from http://github.com/RomanFoell/DRGP-VSS.

Available on the website https://arxiv.org/.

GP-VSS code available from https://github.com/yaringal/VSSGP.

GP-LSTM, LSTM, RNN code available from https://github.com/alshedivat/keras-gp.

DGP-DS code available from https://github.com/ICL-SML/Doubly-Stochastic-DGP.

DGP-RFF code available from https://github.com/mauriziofilippone/deep_gp_ random_features.

Received from http://homes.esat.kuleuven.be/ ˜tokka/daisydata.html.

Description can be found under http://forums.ni.com/t5/NI-myRIO/ myBall-Beam-Classic-Control-Experiment/ta-p/3498079.

Originally received from Global Energy Forecasting Kaggle competitions organized in 2012.

Available from http://github.com/RomanFoell/DRGP-VSS.

Figures and Tables:

3, 3.4, overall resulting in four GP cases: SSGP, VSSGP, SSGP-IP, VSSGP-IP, with the framework introduced in Mattos et al. (2016), to derive our DRGP models: DRGP-SS, DRGP-VSS, DRGP-SS-IP, DRGP-VSS-IP.4.1 DRGP-(V)SS(-IP) MODEL DEFINITIONChoosing the same recurrent structure as inMattos et al. (2016), where now i represents the time, y Hx+1: = [y Hx+1 , . . . , y N ] T ∈ R N , we have

Figure1: Simulation results visualized for the data-sets Drive (first row) and Actuator (second row) for the method DRGP-VSS. First column represents the initial hidden states, blue: first layer, and red: second layer. The second column represents the learned hidden states. The third column shows the simulation results, blue: real data, black: simulation, grey: ±2 times Standard Deviation (SD) and the fourth column the predicted hidden states.

Figure 2: Data-set Drive, Boxplot, 10 runs.

Figure 4: RMSE values for different number of layers L for the data-sets Drive, Actuator, Damper.

For the data-sets Drive and Actuator we chose for our methods DRGP-(V)SS(-IP) the setting L = 2 hidden layers, M = 100 spectral-points and time-horizon H h = H x = 10, which was also used in Mattos et al. (2016), Al-Shedivat et al. (2017) and Dörr et al. (2018) for free simulation (using M pseudo-input points instead of spectral-points). For these two data-sets we filled the results from Mattos et al. (2016); Al-Shedivat et al. (2017) into Table 1. Further, for our methods DRGP-(V)SS(-IP) and DRGP-Nyström we chose on the data-sets Ballbeam and Dryer L = 1, M = 100 and H h = H x = 10. For the data-set Damper we chose L = 2, M = 125 and H h = H x = 10. For the data-set Power Load we chose L = 1, M = 125 and H h = H x = 12. For the data-set Emission we chose L = 1, M = 125 and H h = H x = 10. The other sparse GP, the full GP, DGP-RFF and DGP-DS were trained with NARX structure H x = H y (for the first layer) with the same time horizon as for our DRGPs and with the same amount of pseudo-input points or spectral points.

2 αT m Cnm αm cos(2( αT m (c nm -u m ) + b m )) , for m, m = 1, . . . , M , m = m , with -SS(-IP), THE STATISTICS Ψ 1 , Ψ 2For the versions DRGP-SS(-IP) the statistics are(Ψ 1 ) nm = E q hn [Φ nm ] = E q hn 2σ 2 power (M ) -1 cos(ẑ T m (h n -u m ) + b m ) λnẑm cos(ẑ T m (µ n -u m ) + b m ),for m = 1, . . . , M , n = 1, . . . , N with ẑm = 2π(L -1 z m + p), and 6.3.4 DRGP-(V)SS(-IP), LOWER BOUNDS

log(p(y Hx+1: |X))= logp y Hx+1: , a (l) , Z (l) , h (l) , U (l) Hx+1: , a (L+1) , Z (L+1) , U (L+1) , a (l) , Z (l) , h(l) , U (l) Hx+1: , a (L+1) , Z (L+1) , U (L+1) , a (l) , Z (l) , h (l) , U(l)

a (l) ||p a (l) ) -KL(q Z (l) ||p Z (l) ) def = L REVARB VSS ,where we derive L REVARB SS by being not variational over Z

noise ) 2 (A (l) ) -1 ), withA (l) = Ψ (l) 2 + (σ (l)noise ) 2 I M , respective y Hx+1:6.3.5 PREDICTIONSPredictions for each layer l and new ĥ(l) * are performed in the simple DRGP-(V)SS case withE q f * (l) f * (l) = Ψ (l) 1 * m (l) , V q f * (l) f * (l) = m (l) T Ψ

. . . , L, and fully analog for l = L + 1 by replacing µ (l)

. . . , L, and fully analog for l = L + 1 by replacing µ (l)

Summary of RMSE values for the free simulation results on test data. Best values per dataset are bold. All values are calculated on the original data, unless the data-set Power Load, where the RMSE is shown for the normalized data. Here we have full recurrence for our methods, DRGP-Nyström and GP-LSTM, LSTM, RNN and with auto-regressive part (first layer) for all other GPs. For the column non-rec we turned off the auto-regressive part in the first layer for our methods, DRGP-Nyström and GP-LSTM, LSTM, RNN and also turned off the auto-regressive part (first layer) for all other GPs.

Table content: <table><tr><td colspan="11">methods-data Emission non-rec Power Load Damper Actuator non-rec Ballbeam Dryer Drive non-rec</td></tr><tr><td>DRGP-VSS</td><td>0.104</td><td>0.062</td><td>0.457</td><td>5.825</td><td>0.357</td><td>0.388</td><td colspan="3">0.084 0.109 0.229</td><td>0.268</td></tr><tr><td>DRGP-VSS-IP</td><td>0.119</td><td>0.064</td><td>0.544</td><td>6.112</td><td>0.441</td><td>0.546</td><td colspan="3">0.071 0.107 0.302</td><td>0.293</td></tr><tr><td>DRGP-SS</td><td>0.108</td><td>0.062</td><td>0.497</td><td>5.277</td><td>0.329</td><td>0.563</td><td colspan="3">0.081 0.108 0.226</td><td>0.253</td></tr><tr><td>DRGP-SS-IP</td><td>0.118</td><td>0.065</td><td>0.631</td><td>5.129</td><td>0.534</td><td>0.547</td><td colspan="3">0.076 0.107 0.297</td><td>0.261</td></tr><tr><td>DRGP-Nyström</td><td>0.109</td><td>0.059</td><td>0.493</td><td>6.344</td><td>0.368</td><td>0.415</td><td colspan="3">0.082 0.109 0.249</td><td>0.289</td></tr><tr><td>GP-LSTM</td><td>0.096</td><td>0.091</td><td>0.529</td><td>9.083</td><td>0.430</td><td>0.730</td><td colspan="3">0.062 0.108 0.320</td><td>0.530</td></tr><tr><td>LSTM</td><td>0.098</td><td>0.061</td><td>0.530</td><td>9.370</td><td>0.440</td><td>0.640</td><td colspan="3">0.062 0.090 0.400</td><td>0.570</td></tr><tr><td>RNN</td><td>0.098</td><td>0.066</td><td>0.548</td><td>9.012</td><td>0.680</td><td>0.690</td><td colspan="3">0.063 0.121 0.560</td><td>0.590</td></tr><tr><td>DGP-DS</td><td>0.106</td><td>0.062</td><td>0.543</td><td>6.267</td><td>0.590</td><td>0.576</td><td colspan="3">0.066 0.085 0.422</td><td>0.571</td></tr><tr><td>DGP-RFF</td><td>0.092</td><td>0.069</td><td>0.550</td><td>5.415</td><td>0.520</td><td>0.750</td><td colspan="3">0.074 0.093 0.446</td><td>0.732</td></tr><tr><td>PR-SSM</td><td>N/A</td><td>N/A</td><td>N/A</td><td>N/A</td><td>0.502</td><td>N/A</td><td colspan="3">0.073 0.140 0.492</td><td>N/A</td></tr><tr><td>GP-SSM</td><td>N/A</td><td>N/A</td><td>N/A</td><td>8.170</td><td>N/A</td><td>N/A</td><td>N/A</td><td>N/A</td><td>N/A</td><td>N/A</td></tr><tr><td>GP-VSS</td><td>0.130</td><td>0.058</td><td>0.514</td><td>6.554</td><td>0.449</td><td>0.767</td><td colspan="3">0.120 0.112 0.401</td><td>0.549</td></tr><tr><td>GP-SS</td><td>0.128</td><td>0.060</td><td>0.539</td><td>6.730</td><td>0.439</td><td>0.777</td><td colspan="3">0.077 0.106 0.358</td><td>0.556</td></tr><tr><td>GP-DTC</td><td>0.137</td><td>0.061</td><td>0.566</td><td>7.474</td><td>0.458</td><td>0.864</td><td colspan="3">0.122 0.105 0.408</td><td>0.540</td></tr><tr><td>GP-FITC</td><td>0.126</td><td>0.057</td><td>0.536</td><td>6.754</td><td>0.433</td><td>0.860</td><td colspan="3">0.084 0.108 0.403</td><td>0.539</td></tr><tr><td>GP-full</td><td>0.122</td><td>0.066</td><td>0.696</td><td>9.890</td><td>0.449</td><td>1.037</td><td colspan="3">0.128 0.106 0.444</td><td>0.542</td></tr></table>

Summary of the data-sets for the system identification tasks

Table content: <table><tr><td>parameters \ data-sets</td><td colspan="3">N N train N test</td><td>input dimension</td><td>output dimension</td></tr><tr><td>Drive</td><td>500</td><td>250</td><td>250</td><td>1</td><td>1</td></tr><tr><td>Dryer</td><td>1000</td><td>500</td><td>500</td><td>1</td><td>1</td></tr><tr><td>Ballbeam</td><td>1000</td><td>500</td><td>500</td><td>1</td><td>1</td></tr><tr><td>Actuator</td><td>1024</td><td>512</td><td>512</td><td>2</td><td>1</td></tr><tr><td>Damper</td><td>3499</td><td colspan="2">2000 1499</td><td>1</td><td>1</td></tr><tr><td>Power Load</td><td>9518</td><td colspan="2">7139 2379</td><td>11</td><td>1</td></tr><tr><td colspan="4">Emission 12500 10000 2500</td><td>6</td><td>1</td></tr></table>

The separation of the data-sets Drive, Dryer, Ballbeam, Actuator, Damper, Power Load in training-and test-data was given by the papers we use for comparison. The separation of the Emission data-set was chosen by ourself.

Table content: <table><tr><td>Under review as a conference paper at ICLR 2019</td></tr><tr><td>6.2 NONLINEAR SYSTEM IDENTIFICATION</td></tr></table>

Summary of the methods for the system identification tasks

Table content: <table><tr><td>method</td><td>references</td></tr><tr><td>DRGP-VSS</td><td>introduced in this paper</td></tr><tr><td>DRGP-VSS-IP</td><td>introduced in this paper</td></tr><tr><td>DRGP-SS</td><td>introduced in this paper</td></tr><tr><td>DRGP-SS-IP</td><td>introduced in this paper</td></tr><tr><td>DRGP-Nyström</td><td>(Mattos et al., 2016)</td></tr><tr><td>GP-LSTM</td><td>(Al-Shedivat et al., 2017)</td></tr><tr><td colspan="2">LSTM (Hochreiter &amp; Schmidhuber, 1997)</td></tr><tr><td>RNN</td><td>see e.g. (Al-Shedivat et al., 2017)</td></tr><tr><td>DGP-DS</td><td>(Salimbeni &amp; Deisenroth, 2017)</td></tr><tr><td>DGP-RFF</td><td>Cutajar et al. (2016)</td></tr><tr><td>PR-SSM</td><td>(Dörr et al., 2018)</td></tr><tr><td>GP-SSM</td><td>(Svensson et al., 2016)</td></tr><tr><td>GP-VSS</td><td>(Gal &amp; Turner, 2015)</td></tr><tr><td>GP-SS</td><td>(Lázaro-Gredilla et al., 2010)</td></tr><tr><td>GP-DTC</td><td>(Williams &amp; Seeger, 2000)</td></tr><tr><td>GP-FITC</td><td>(Snelson &amp; Ghahramani, 2006)</td></tr><tr><td>GP-full</td><td>(Rasmussen, 2006)</td></tr></table>

Peer review (iclr2019):
 Overall Score: 7/10. Confidence Score: 7/10. Detailed Comments: This paper introduces various Deep Recurrent Gaussian Process (DRGP) models based on the Sparse Spectrum Gaussian Process (SSGP) models and the Variational Sparse Spectrum Gaussian Process (VSSGP) models. This is a good paper and proposed models are very sound so I recommend for acceptance although as main weakness I can say that is very technical so it can be difficult to follow. Adding more intuitive ideas, motivation and maybe a figure for each step would be a solution. Apart from that it is a really good paper, congratulations. Related to: RNN models and Sparse Nystrom approximation. Strengths: Models are very sound, solutions are solid, the proposed methodology is correct and the empirical results and experiments are valid and properly done. Weaknesses: It is too difficult to follow and it is written in an extreme technical way. More intuitions and a proper motivation both in the abstract and introduction may be put in order to make the paper easier to read and, hence, more used by researchers and data scientists. Does this submission add value to the ICLR community? : Yes it does, the experiments show the efficiency of the proposed methods in some scenarios and are valid methodologies. Quality: Is this submission technically sound?: Yes it is. Are claims well supported by theoretical analysis or experimental results?: Experimental results prove empirically the methods and appendixes show the analysis performed in a clear and elegant way. Is this a complete piece of work or work in progress?: Complete piece of work. Are the authors careful and honest about evaluating both the strengths and weaknesses of their work?: Yes, and I would enfatize that I have liked that some experiments are won by other methods such as GP-LSTM, they are very honest. Clarity: Is the submission clearly written?: Yes, but it is difficult for newcomers due to the reasons that I have stated before. Is it well organized?: Yes it is. Does it adequately inform the reader?: Yes it is. Originality: Are the tasks or methods new?: Yes, they are sound. Is the work a novel combination of well-known techniques?: Yes it is. Is it clear how this work differs from previous contributions?: Yes. Is related work adequately cited?: Yes, being a strength of the paper. Significance: Are the results important?: I would argue that they are and are a clear alternative to consider in order to solve these problems. Are others likely to use the ideas or build on them?: If the paper is written in a more friendly way, yes. Does the submission address a difficult task in a better way than previous work?: Yes I think. Does it advance the state of the art in a demonstrable way?: Yes, empirically. Arguments for acceptance: Models are very sound, solutions are solid, the proposed methodology is correct and the empirical results and experiments are valid and properly done Arguments against acceptance: Clarity of the paper. Minor issues and typos: -> (V)SS not defined before being used. -> Abstract should be rewritten adding a motivation and focusing more on the problems being solved and less in the details of the solutions. -> Recurrent indexes that go backwards (i) of Eq. 1. should be explained why are going backwards before being used like that. Newcomers may be confused. -> Section 2 writing style lacks a bit of cohesion, relating the paragraphs may be a solution. -> Q is not defined in section 3.1 paragraph 1. -> A valid covariance function must produce a PSD matrix, put that in section 3.1. -> I do not see how U marginalizes in Eq. 7, kind of confused about that, I think that it should be p(y|X,U). -> Section 3.4 statistics should be explained. Reading thread and authors response rebuttal decision: ================================================= I consider that the authors have perfomed a good rebuttal and reading the other messages and the authors response I also consider that my issue with clarity is solved. Hence, I upgrade my score to 7 and recommend the paper for publication.