Year: 2018

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.

heading: 1 INTRODUCTION
text: 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
hi = ζ(hi−1, . . . ,hi−H ,xi−1, . . . ,xi−H) + h i , yi = ψ(hi, . . . ,hi−H) + y i ,
(1) (2)
where xi is an external input, yi is an output observation, hi 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 hi , 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 Lá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); Dö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(NM2Qmax(L+ 1)) down toO(M3) for N ≤M , M the sparsity parameter, (L+ 1) the amount of GPs and Qmax 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.

heading: 2 RELATED WORK TO GAUSSIAN PROCESSES WITH SSM AND DGPS
text: 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. Turner et al. (2010) solved this problem, but they have complicated approximate training and inference stages and the model is hard to scale. Frigola et al. (2014) used a GP for transition, while the observation is parametric. Svensson et al. (2016) used an approximation of the spectral representation by Bochner’s theorem in a particular form and with a reduced rank structure for the transition function. They realize inference in a fully Bayesian approach over the amplitudes and the noise parameters. The construction of Eleftheriadis et al. (2017) involves a variational posterior that follows the same Markov properties as the true state with rich representational capacity and which has a simple, linear, time-varying structure. Dörr et al. (2018) introduced a GP-SSM with scalable training based on doubly stochastic variational inference for robust training. Our models extend the GP-SSM framework by defining the transition as a DRGP based on our newly derived (V)SS approximations in the Sections 3.3, 3.4, where the latent (non-observed) output-data is learned as a hidden state. We refer to the report Föll et al. (2017)2 for a detailed but preliminary formulation of our models and experiments. Following Damianou & Lawrence (2013), a Deep Gaussian Process (DGP) is a model assuming
yi = f (L+1)(f(L)(f(L−1)(. . . (f(1)(xi) + 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)(xi) +
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,
1DRGP-(V)SS code available from http://github.com/RomanFoell/DRGP-VSS. 2Available on the website https://arxiv.org/.
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-Nyström, based on the variational sparse Nyströ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-Nyströ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).

heading: 3 GAUSSIAN PROCESSES AND VARIATIONAL SPARSE SPECTRUM GP
text: 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 (Lá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, fx, y, (italic) for random variables, a, f(x), y (upright) for realizations/samples and data.

heading: 3.1 GAUSSIAN PROCESSES
text: A stochastic process f is a GP if and only if any finite collection of random variables fX def = [fx1 , . . . , fxN ] T forms a Gaussian random vector (Rasmussen, 2006). A GP is completely defined by its mean function m : RQ → R,x 7→ m(x), Q ∈ N the input-dimension, and covariance function k : RQ × RQ → R, (x,x′) 7→ k(x,x′) (Kallenberg, 2006, Lemma 11.1), where
m(x) def = E [fx] , k(x,x ′) def = cov(fx, fx′) = E [(fx −m(x))(fx′ −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 matrixKNN def = (k(xi,xj)) N i,j=1 ∈ RN×N , when filling in combinations of data-input points xi, i = 1, . . . , N . Let y def= [y1, . . . , yN ]T ∈ RN , X def = [x1, . . . ,xN ]
T ∈ RN×Q and we assume yi = f(xi) + yi , where yi ∼ N (0,σ2noise), for i = 1, . . . , N , and our aim is to model any set of function values f def = [f(x1), . . . , f(xN )]
T ∈ RN at X as samples from a random vector fX. Moreover, we assume the prior fX|X ∼ N (0,KNN ), meaning that any set of function values f given X are jointly Gaussian distributed with mean 0 ∈ RN and a covariance matrix KNN . The predictive distribution pfx∗ |x∗,X,y for a test point x ∗ ∈ RQ, where KN∗ def = (k(xi,x
∗))Ni=1 ∈ RN , and analogously K∗∗, K∗N , can be derived through the joint probability model and conditioning as
fx∗ |x∗,X,y ∼ N (K∗N (KNN + σ2noiseIN )−1y,K∗∗ −K∗N (KNN + σ2noiseIN )−1KN∗).
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)
k(x,x′) = σ2power (∏Q q=1 exp ( −1 2 l−2q (xq − x′q)2 )) cos ( 2π ∑Q q=1 p−1q (xq − x′q) ) , (3)
with an amplitude σ2power and length scales lq , pq ∈ R, q = 1, . . . , Q. As pq →∞, this corresponds to the squared exponential (SE) covariance function in the limit (Gal & Turner, 2015).

heading: 3.2 VARIATIONAL SPARSE SPECTRUM GP
text: We introduce the SSGP following Gal & Turner (2015). For a stationary covariance function k on RQ × RQ there exists a function ρ : RQ → R, τ 7→ ρ(τ ), such that k(x,x′) = ρ(x − x′) for x,x′ ∈ RQ. 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 τ = x− x′, can be expressed via Monte Carlo Approximation (MCA) following Gal & Turner (2015), Section 2, as
ρ(τ ) = ∫ RQ e2πiz T τdµ(z)
≈ σ2power
M ∑M m=1 2 cos(2πzTm(x− um) + bm) cos(2πzTm(x′ − um) + bm) def = k̃(x,x′).
(4)
We refer to zm as the spectral points, bm as the spectral phases and um as the pseudo-input/inducing points form = 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πlq] Q q=1) ∈ RQ×Q, a scaling vector p def= [p−11 , . . . ,p −1 Q ] T ∈ RQ and Z def= [z1, . . . , zM ]T ∈ RQ×M via
φ(x,Z) def = √ 2σ2powerM −1 [ cos(2π(L−1z1 + p) T (x− u1) + b1), . . . ,
cos(2π(L−1zM + p) T (x− uM ) + bM ) ]T ∈ RM . (5)
Here we sample b ∼ Unif [0, 2π], z ∼ N (0, IQ) and we set K̃(SM)NN def = ΦΦT with Φ def= [φ(x1,Z), . . . , φ(xN ,Z)] T ∈ RN×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 ∈ RN×1. By replacing the covariance function with the sparse covariance function k̃(SM) and setting the priors to
pZ def = ∏M
m=1 pzm , where zm ∼ N (0, IQ), pa, where a ∼ N (0, IM ), (6)
for m = 1, . . . ,M , where we have y|a,Z,U ,X ∼ N (Φa,σ2noiseIN ) with U def = [u1, . . . ,uM ] T ∈ RQ×M (we do not define priors on p, L−1, U, b def= [b1, . . . ,bM ]T ∈ RM ), we can expand the marginal likelihood (ML) to
p(y|X) = ∫ p(y|a,Z,U,X)p(a)p(Z)dadZ, (7)
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
qZ def = ∏M
m=1 qzm , where zm ∼ N (αm,βm), qa, where a ∼ N (m, s),
with βm ∈ RQ×Q diagonal, for m = 1, . . . ,M , and s ∈ RM×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)) ≥ EqaqZ [log(p(y|a,Z,U,X))]−KL(qa||pa)−KL(qZ ||pZ)). (8) As usual, KL defines the Kullback-Leibler (KL) divergence. By proving the existence of an optimal distribution qopta 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ΨT1 y,σ2noiseA−1), with A = Ψ2 + σ2noiseIM , Ψ1 = EqZ [Φ] ∈ RN×M , Ψ2 = EqZ [ ΦTΦ ] ∈ RM×M , we can derive the optimal bound case.

heading: 3.3 (V)SSGP WITH REGULARIZATION PROPERTIES VIA INDUCING POINTS (IP)
text: 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
pa|U , where a|U ∼ N (0,KMM ), py|a,Z,U ,X , where y|a,Z,U ,X ∼ N (ΦK−1MMa,KNN −KNMK −1 MMKMN + σ 2 noiseIN ), (9) (10)
where KNN , KNM , KMM 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 KNN , plus the discrepancy of the two wellknown sparse covariance matrices for the Sparse Spectrum (not exactly, as we have ΦK−1MMΦ
T ) and the Nyströ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(qZ ||pZ) and Ψ1 = Φ, Ψ2 = ΦTΦ) can be calculated in the same way as Titsias & Lawrence (2010) until Equation (14) in closed form:
log(p(y|X)) ≥ − (N −M) 2 log(σ2noise)− N 2 log(2π)− y Ty 2σ2noise +
log(|KMM | ∣∣A−1∣∣)
2
+ yTΨ1A −1ΨT1 y 2σ2noise − tr(KNN ) 2σ2noise + tr(K−1MMKMNKNM ) 2σ2noise −KL(qZ ||pZ).
(11)
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 + σ2noiseKMM , which involves KMM , the chosen covariance matrix filled in with the pseudo-input points U, and three extra terms
log(|KMM |) 2 , tr(KNN ) 2σ2noise , tr(K−1MMKMNKNM ) 2σ2noise .

heading: 3.4 VARIATIONAL APPROXIMATION OF THE INPUT-SPACE FOR (V)SSGP(-IP)
text: 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 = [h1, . . . ,hN ]
T ∈ RN×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
pH def = ∏N
i=1 phi , where hi ∼ N (0, IQ), qH
def = ∏N
i=1 qhi , where hi ∼ N (µi,λi),
λi ∈ RQ×Q diagonal, for i = 1, . . . , N . As a consequence, for VSSGP-IP we overall derive statistics Ψ0 = tr (EqH [KNN ]) = Nσ 2 power,
Ψ1 = EqZqH [Φ] ∈ RN×M , Ψ2 = EqZqH [ ΦTΦ ] ∈ RM×M and Ψreg = EqH [KMNKNM ] as defined in the Appendix 6.3. These statistics are essentially the given matrices Φ, ΦTΦ,KMNKNM from the beginning, but every input hi and every spectral point zm 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:
log(p(y|H)) ≥ − (N −M) 2 log(σ2noise)− N 2 log(2π)− y Ty 2σ2noise +
log(|KMM | ∣∣A−1∣∣)
2
+ yTΨ1A −1ΨT1 y
2σ2noise − Ψ0 2σ2noise + tr(K−1MMΨreg) 2σ2noise −KL(qZ ||pZ)−KL(qH ||pH).
(12)
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 + σ2noiseIM and eliminate log(|KMM |) 2 , Ψ0 2σ2noise
, tr(K−1MMΨreg)
2σ2noise in the lower bound (12).

heading: 4 DRGP WITH VARIATIONAL SPARSE SPECTRUM APPROXIMATION
text: In this section we want to combine our newly derived (V)SS-(IP) approximations in the Sections 3.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.

heading: 4.1 DRGP-(V)SS(-IP) MODEL DEFINITION
text: Choosing the same recurrent structure as in Mattos et al. (2016), where now i represents the time, yHx+1: = [yHx+1, . . . , yN ] T ∈ RN̂ , we have
ζ : h (l) i = f (l)(ĥ (l) i ) + h(l) i , with prior f (l) Ĥ(l) |Ĥ(l) ∼ N (0,K(l) N̂N̂ ), l = 1, . . . , L
ψ : yi = f (l)(ĥ (l) i ) + y i , with prior f (l) Ĥ(l) |Ĥ(l) ∼ N (0,K(l) N̂N̂ ), l = L+ 1,
with ζ, ψ in Equation (1) and (2), h (l) i ∼ N (0, (σ (l) noise) 2), yi ∼ N (0, (σ (L+1) noise ) 2) and N̂ = N−Hx, for i = Hx + 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 Ĥ(l) = [ĥ(l)Hx+1, . . . , ĥ (l) N ] T is specified as
ĥ (l) i def = 
[ h (1)
i−1,xi−1
]T def = [[
h (1) i−1, . . . ,h (1) i−Hh ] , [xi−1, . . . ,xi−Hx ] ]T , l = 1[
h (l) i−1,h (l−1) i
]T def = [[
h (l) i−1, . . . ,h (l) i−Hh
] , [ h
(l−1) i , . . . ,h (l−1) i−Hh+1
]]T , l = 2, . . . , L
h (L) i def = [ h (L) i , . . . ,h (L) i−Hh+1 ]T , l = L+ 1,
(13)
where ĥ(1)i ∈ RHh+HxQ, ĥ (l) i ∈ R2Hh for l = 2, . . . , L, ĥ (L+1) i ∈ RHh , for i = Hx + 1, . . . , N . For simplification we set H def= Hx = Hh 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
p yHx+1:[a(l),Z(l),h(l),U(l)] L+1 l=1 |X = ∏L+1 l=1 p h (l) Hx+1: |a(l),Z(l),Ĥ(l),U(l)pa(l)pZ(l)ph̃(l) ,
with ph̃(L+1) = 1, h (L+1) Hx+1: def = yHx+1:, h (L+1) = {} and h̃(l) = [h(l)1+Hx−Hh , . . . ,h (l) Hh ]T ∈ R2Hh−Hx . Here the priors are
pa(l) , where a (l) ∼ N (0, IM ), pZ(l)
def = ∏M
m=1 p z (l) m , where z(l)m ∼ N (0, IQ),
ph̃(l) , where h̃ (l) ∼ N (0, I2Hh−Hx),
the product of them is defined as PREVARB, and the variational distributions are
qa(l) , where a (l) ∼ N (m(l), s(l)), qZ(l)
def = ∏M
m=1 q z (l) m , where z(l)m ∼ N (α(l)m ,β (l) m ), qh(l) def = ∏N
i=1+Hh−Hx q h (l) i , where h(l)i ∼ N (µ (l) i , λ (l) i ),
the product of them is defined as QREVARB, and where β (l) m ∈ RQ×Q is diagonal, for i = 1 +Hx −Hh, . . . , 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
pa(l)|U(l) , where a (l)|U (l) ∼ N (0,K(l)MM ),
p h
(l) Hx+1:
|a(l),Z(l),Ĥ(l),U(l) , where h (l) Hx+1: |a(l),Z(l), Ĥ(l),U (l) ∼
N (Φ(l)(K(l)MM ) −1a(l),K (l) NN −K (l) NM (K (l) MM ) −1K (l) MN + (σ (l) noise) 2IN ).
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.

heading: 4.2 DRGP-(V)SS(-IP) VARIATIONAL EVIDENCE LOWER BOUND (ELBO)
text: 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
log(p(yHx+1:|X)) ≥ EQREVARB [ ∑L+1
l=1 log(p(h (l) Hx+1: |a(l),Z(l), Ĥ(l),U(l)))]−KL(QREVARB||PREVARB) def = LREVARB(V)SS , (14)
and for the (V)SS-IP approximations, denoted as REVARB-(V)SS-IP, is given by
log(p(yHx+1:|X)) ≥
Ep a(l)|U(l)
[ ∑L+1
l=1 exp(〈log(N (h(l)Hx+1:|Φ (l)(K (l) MM ) −1a(l), (σ (l) noise) 2IN ))〉QREVARB )]
− ∑L+1
l=1
Ψ (l) 0
2(σ(l)noise) 2 −
tr((K(l)MM ) −1Ψ (l) reg)
2(σ(l)noise) 2
−KL(QREVARB||PREVARB) def = LREVARB(V)SS(-IP), (15)
where 〈·, ·〉 means, taking the expectation under the integral. Additionally, for the simple (V)SS approximations in 14, the optimal bound LREVARB(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 QREVARB 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 LREVARB(V)SS-opt/IP requires O(NM2Qmax(L+ 1)), where Qmax = maxl=1...,L+1Q(l), Q(l) def = dim(ĥ
(l) i ) and ĥ (l) i from Equation (13). DVI can reduce the complexity to O(M3) 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.

heading: 4.3 MODEL PREDICTIONS
text: After model optimization based on the lower bounds in the Equation (14) and (15), model predictions for new ĥ(l)∗ in the REVARB-(V)SS(-IP) framework can be obtained based on the approximate, variational posterior distribution QREVARB. 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.

heading: 5 EXPERIMENTS
text: In this section we want to compare our methods DRGP-SS, DRGP-VSS (optimal bound case) and DRGP-SS-IP, DRGP-VSS-IP, against other well known sparse GPs and the full GP with NARX structure, the DRGP-Nyström of Mattos et al. (2016), the GP-LSTM of Al-Shedivat et al. (2017), the LSTM of Hochreiter & Schmidhuber (1997), a simple RNN, the DGP-DS of Salimbeni & Deisenroth (2017) and the DGP-RFF of Cutajar et al. (2016), both with NARX structure for the first layer, the GP-SSM of Svensson et al. (2016) and the PR-SSM Dörr et al. (2018). The full GP is named GP-full, the FITC approximation of Snelson & Ghahramani (2006) is named GP-FITC, the DTC approximation of Williams & Seeger (2000) is named GP-DTC, the SSGP of Lázaro-Gredilla et al. (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.

heading: 5.1 IMPLEMENTATION
text: 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-Nyström, GP-SS, GP-DTC, GPFITC, 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 code3456. 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.
3GP-VSS code available from https://github.com/yaringal/VSSGP. 4GP-LSTM, LSTM, RNN code available from https://github.com/alshedivat/keras-gp. 5DGP-DS code available from https://github.com/ICL-SML/Doubly-Stochastic-DGP. 6DGP-RFF code available from https://github.com/mauriziofilippone/deep_gp_
random_features.

heading: 5.2 BENCHMARK DATA-SETS MODEL LEARNING AND COMPARISON
text: 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-Nyströ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-Nyström on these data-sets, in particular for the setting of time horizons of Table 1 with H = 10. In Figure 4 the RMSE results for different layers L on the data-sets Drive, Actuator and Damper are shown. We can observe that different layers L are favored.

heading: 5.3 CONCLUSION
text: In this paper we introduced four new DRGPs based on the SS approximation introduced by LázaroGredilla 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.

heading: ACKNOWLEDGMENTS
text: We would like to acknowledge support for this project from the ETAS GmbH.

heading: 6 APPENDIX
text: 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.

heading: 6.1 DATA-SET DESCRIPTION
text: 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 Ballbeam78 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 Sjö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 Table 2. 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.
7Received from http://homes.esat.kuleuven.be/˜tokka/daisydata.html. 8Description can be found under http://forums.ni.com/t5/NI-myRIO/
myBall-Beam-Classic-Control-Experiment/ta-p/3498079. 9Originally received from Global Energy Forecasting Kaggle competitions organized in 2012.
10Available from http://github.com/RomanFoell/DRGP-VSS.

heading: 6.2 NONLINEAR SYSTEM IDENTIFICATION
text: The methods for the system identification tasks and their references are summarized in Table 2. 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 Hh = Hx = 10, which was also used in Mattos et al. (2016), Al-Shedivat et al. (2017) and Dö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-Nyström we chose on the data-sets Ballbeam and Dryer L = 1, M = 100 and Hh = Hx = 10. For the data-set Damper we chose L = 2, M = 125 and Hh = Hx = 10. For the data-set Power Load we chose L = 1, M = 125 and Hh = Hx = 12. For the data-set Emission we chose L = 1, M = 125 and Hh = Hx = 10. The other sparse GP, the full GP, DGP-RFF and DGP-DS were trained with NARX structure Hx = Hy (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. 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 Hx = Hy 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 Hx = Hy 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 Nyströ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, IQ(l)), an initialization for the pseudo-input points with a random subset of sizeM 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
(V)SS approximations were also initialized as SE covariance function, see Equation 5. For all methods we used automatic relevance determination, so each input dimension has its own length scale. For our methods DRGP-(V)SS(-IP) and DRGP-Nyström the noise parameters and the hyperparameters were initialized by σ(l)noise = 0.1, σ (l) power = 1 and the length scales by either
l (l) q = √ max(Ĥ (l) q )−min(Ĥ(l)q ) or l(l)q = max(Ĥ(l)q ) −min(Ĥ(l)q ), for all q, l, where Ĥ(l)q is the
data-vector containing the q-th input-dimension values of every input-data point ĥ(l)i , for all i. Furthermore, we initialized the latent hidden GP states with the output-observation data y. The other standard GPs were also initialized with σnoise = 0.1, σpower = 1 and the same initialization for length scales with respect to the NARX structure input data as before. For LSTM, RNN we used the initialization for the weights provided by Keras, a Deep Learning library for Theano and TensorFlow. For GP-LSTM we used the initialization for the weights provided by Keras and σnoise = 0.1, σpower = 1 and for the length scale initialization we chose lq = 1 for all input-dimensions. For DGP-RFF we used the initialization coming from the implementation with σ(l)power = 1 and for the length scale initialization with l(l)q = 0.5 log(Q(l)) for all l, q. To our knowledge DGP-RFF has no parameter σ(l)noise for all l. For DGP-DS we used the initialization σ(l)noise = 0.1 and σ (l) power = 1 and for the length scale initialization we chose l(l)q = 1 for all l, q. For all our implementations, see Section 5.1, we used the positive transformation x′ = log(1 + exp(x))2 for the calculation of the gradients in order for the parameters to be constrained positive and with L-BFGS optimizer, either from Matlab R2016b with fmincon, or Python 2.7.12 with scipy optimize. All methods were trained on the normalized data x 7→ x−µσ2 (for every dimension independently), several times (same amount per data-set: the initializations are still not deterministic, e.g. for pseudo-inputs points and spectral points) with about 50 to 100 iterations and the best results in RMSE value on the test-data are shown in Table 1. For our methods DRGP-(V)SS(-IP) and DRGP-Nyström we fixed σ(l)noise, σ (l) power for all l (optional the spectral points/pseudo-input points for DRGP-(V)SS; for the IP case we never excluded the pseudo-input points because we would fall back to the DRGP-(V)SS case; for DRGP-Nyström we always included the pseudo-input points) during the first iterations to independently train the latent states. For all other GPs we also tested with fixed and not fixed σnoise = 0.1, except GP-LSTM and DGP-DS. For DRGP-VSS(-IP) we fixed β(l)m for all m, l to small value around 0.001, as well as the spectral phases bm for all m, l sampling from Unif [0, 2π] (this seems to work better in practice). The limitations for β(l)m also holds for GP-VSS as well. We want to remark at this point that setting um = 0 for all m = 1, . . . ,M worked sometimes better than choosing a subset from the input-data (not in the IP case). This seems to be different to Gal & Turner (2015), who pointed out: ’These are necessary to the approximation. Without these points (or equivalently, setting these to 0), the features would decay quickly for data points far from the origin (the fixed point 0).’ For GP-SSM we show the result of the data-set Damper from their paper in Table 1. For PR-SSM we show the results from their paper in Table 1. We show additional results for the data-sets Drive, Actuator and Emission with missing autoregressive part for the first layer for our methods DRGP-(V)SS(-IP) and DRGP-Nyström in Table 1, named non-rec. For the sparse GP, the full GP, GP-LSTM, DGP-RFF and DGP-DS and the data-sets Drive, Actuator and Emission we show the results with missing auto-regressive part in Table 1, just modeling the data with exogenous inputs. Here we want to examine the effect of the auto-regressive part of the first layer for the DRGP models on the RMSE. GP-SSM and PR-SSM are not listed for this setting of recurrence.

heading: 6.3 VARIATIONAL APPROXIMATION AND INFERENCE FOR DRGP-(V)SS(-IP)
text: 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.),
• 12 (cos(a− b+ x− y) + cos(a+ b+ x+ y)) = cos(x+ a) cos(y + b).
• JI, for Jensen’s inequality.

heading: 6.3.1 DRGP-VSS(-IP), THE STATISTICS Ψ1, Ψ2
text: For the versions DRGP-VSS(-IP) the statistics are
(Ψ1)nm = Eqzmqhn [Φnm]
B.1 = Eqhn [√ 2σ2power(M) −1e− 1 2 h̄ T nmβmh̄nm cos(α̂ T m(hn − um) + bm) ] B.2 = Σ1mZnme − 12 α̂ T mCnmα̂m cos(α̂ T m(cnm − um) + bm),
for m,= 1, . . . ,M , n = 1, . . . , N with
h̄nm = 2πL −1(hn − um),
α̂m = 2π(L −1αm + p),
cnm = Cnm(βm(2π) 2L−2um + λ −1 n µn),
Cnm = (βm(2π) 2L−2 + λ−1n ) −1, vnm = um − µn, Vnm = (2π) −2L2β−1m + λn, Znm = 1√ |Vnm| e− 1 2v T nmV −1 nmvnm ,
Σ1m = √ 2σ2power ∏Q q=1 ( l2q βmq ) (M)−1,
and
Ψ2 = ∑N n=1(Ψ2) n, where
(Ψ2) n mm′ = Eqzmqzm′ qhn
[ ΦTnmΦnm′ ] B.1 = Eqhn [ 2σ2power(M) −1e− 1 2 (h̄ T nmβmh̄nm+h̄ T nm′βm′ h̄nm′ )
cos(α̂ T m(hn − um) + bm) cos(α̂ T m′(hn − um′) + bm′) ] B.2 = Σ2mm′Z n mm′ ( e− 1 2 ᾱ T mm′D n mm′ ᾱmm′ cos(ᾱTmm′d n mm′ − τ̄mm′ + b̄mm′)
+e− 1 2
+ α
T mm′D n mm′ + αmm′ cos( + α T
mm′d n mm′ − + τmm′ +
+ bmm′) ) .
for m,m′ = 1, . . . ,M , m 6= m′, with
bmm′ = bm − bm′ , + bmm′ = bm + bm′ ,
τ̄mm′ = α̂ T mum − α̂ T m′um′ ,
+ τmm′ = α̂ T mum + α̂ T m′um′ , ᾱmm′ = α̂m − α̂m′ , + αmm′ = α̂m + α̂m′ , bmm′ = Bmm′(2π) 2L−2(βmum + βm′um′),
βmm′ = βm + βm′ , Bmm′ = (2π) −2L2β−1mm′ ,
dnmm′ = D n mm′(B −1 mm′bmm′ + λ −1 n µn), Dnmm′ = (B −1 mm′ + λ −1 n ) −1,
wnmm′ = bmm′ − µn, Wnmm′ = Bmm′ + λn,
umm′ = um − um′ , Umm′ = (2π) −2L2(β−1m + β −1 m′ ),
Znmm′ = 1√
|Wnmm′ | |Umm′ | e− 1 2 (w n mm′ T Wn mm′ −1wn mm′+u T mm′U −1 mm′umm′ ),
Σ2mm′ = σ 2 power ∏Q q=1 ( l2q√
βmqβm′q
) (M)−1,
and
(Ψ2) n mm = Eqzmqhn [ ΦTnmΦnm ] B.1 = Eqhn [ 2σ2power(M) −1( 1
2 +
1 2 e−2h̄ T nmβmh̄nm cos(2(α̂ T m(hn − um) + bm)) ] B.2 = σ2power(M) −1 ( 1 + Σ̃2mZ̃nme −2α̂TmC̃nmα̂m cos(2(α̂ T m(cnm − um) + bm)) ) ,
for m,m′ = 1, . . . ,M , m = m′, with
Σ̃2m = √∏Q q=1 ( l2q βmq ) 2−Q.

heading: 6.3.2 DRGP-SS(-IP), THE STATISTICS Ψ1, Ψ2
text: For the versions DRGP-SS(-IP) the statistics are
(Ψ1)nm = Eqhn [Φnm]
= Eqhn [√ 2σ2power(M) −1 cos(ẑTm(hn − um) + bm) ]
B.2 = √
2σ2power(M) −1e− 1 2 ẑ T mλnẑm cos(ẑTm(µn − um) + bm),
for m = 1, . . . ,M , n = 1, . . . , N with
ẑm = 2π(L −1zm + p),
and
Ψ2 = ∑N n=1(Ψ2) n, where
(Ψ2) n mm′ = Eqhn [ ΦTnmΦnm′ ] B.1 = Eqhn [ 2σ2power(M) −1 cos(ẑTm(hn − um) + bm) cos(ẑTm′(hn − um′) + bm′) ]
= σ2power(M) −1 ( e− 1 2 z̄ T mm′λnz̄mm′ cos(z̄Tmm′µn − ρ̄mm′ + b̄mm′)
+e− 1 2
+ z T
mm′λn + zmm′ cos( + z T mm′µn − + ρmm′ + + bmm′)) ) ,
for m,m′ = 1, . . . ,M with ρ̄mm′ = ẑ T mum − ẑTm′um′ ,
+ ρmm′ = ẑ T mum + ẑ T m′um′ , z̄mm′ = ẑm − ẑm′ , + zmm′ = ẑm + ẑm′ ,
for the other variables, see the defined variables in the DRGP-VSS case.
6.3.3 DRGP-(V)SS-IP, THE STATISTICS Ψreg AND Ψ0
For Ψreg = EqH [KMNKNM ] see Ψ2 in Section 6.3.1 but setting bm = 0, αm = 0 and βm = 1 for all m = 1, . . . ,M .
Ψ0 naturally is given by Ψ0 = tr (EqH [KNN ]) = Nσ 2 power because of the chosen SE covariance function.

heading: 6.3.4 DRGP-(V)SS(-IP), LOWER BOUNDS
text: In this section we derive the different variational lower bounds for our models DRGP-(V)SS(-IP). We first show the bound LREVARB(V)SS without optimal variational distribution for a(l). Then the bounds LREVARB(V)SS-opt with optimal variational distribution for a(l) follows, as well as the LREVARB(V)SS-IP case. We use the simplified notation dAdZdH = da(1) . . . da(L+1)dZ(1) . . . dZ(L+1)dh(1) . . . dh(L).
log(p(yHx+1:|X))
= log (∫ p ( yHx+1:, [ a(l),Z(l),h(l),U (l) ]L+1 l=1 |X ) dAdZdH ) = log
(∫ QREVARB QREVARB p ( yHx+1:,a (L+1),Z(L+1),U (L+1), [ a(l),Z(l),h(l),U (l) ]L l=1 |X ) dAdZdH ) JI ≥ ∫ QREVARB log p ( yHx+1:,a (L+1),Z(L+1),U (L+1), [ a(l),Z(l),h(l),U (l) ]L l=1 |X )
QREVARB  dAdZdH =
L+1∑ l=1 ∫ q(a(l))q(Z(l))q(h(l))p(h (l) Hx+1: |a(l),Z(l), Ĥ(l),U (l))da(l)dZ(l)dh(l)
− L∑ l=1 N̂2 − N∑ i=1+Hx−Hh log(2πλ (l) i ) 2 + Hh∑ i=1+Hx−Hh log(2π) 2 + ( λ (l) i + ( µ (l) i )2) 2  − L+1∑ l=1 KL(qa(l) ||pa(l))−KL(qZ(l) ||pZ(l))
= −N̂ 2 L+1∑ l=1 log(2π(σ (l) noise) 2)− yTHx+1:yHx+1: 2 ( σ(L+1)noise
)2 + yTHx+1:Ψ (L+1) 1 m
(L+1)( σ(L+1)noise
)2 − tr ( Ψ(L+1)2 (s (L+1) + m(L+1)(m(L+1))T ) ) 2 ( σ(L+1)noise
)2 +
L∑ l=1 − 1 2 ( σ
(l) noise
)2 (( N∑ i=Hx+1 λ (l) i ) + (µ(l))THx+1:µ (l) Hx+1: )
+
( (µ
(l) Hx+1: )TΨ (l) 1 m
(l) )
( σ
(l) noise
)2 − tr ( Ψ (l) 2 (s (l) + m(l)(m(l))T ) )
2 ( σ
(l) noise
)2
−N̂ 2 + N∑ i=1+Hx−Hh log(2πλ (l) i ) 2 − Hh∑ i=1+Hx−Hh log(2π) 2 −
( λ
(l) i +
( µ
(l) i )2) 2  − L+1∑ l=1 KL(qa(l) ||pa(l))−KL(qZ(l) ||pZ(l))
def = LREVARBVSS ,
where we derive LREVARBSS by being not variational overZ(l). Using the optimal distribution for a(l) ∼ N ((A(l))−1(Ψ(l)1 )Tµ (l) Hx+1: , (σ (l) noise) 2(A(l))−1), with A(l) = Ψ(l)2 + (σ (l) noise) 2IM , respective yHx+1:
for L+ 1 we obtain
≥ −N̂ −M 2 L+1∑ l=1 log (( σ (l) noise )2) − yTHx+1:yHx+1: 2 ( σ(L+1)noise
)2 + yTHx+1:Ψ (L+1) 1 (A (L+1))−1(Ψ(L+1)1 ) TyHx+1:
2 ( σ(L+1)noise )2 + log(|(A(L+1))−1|)2 +
L∑ l=1 − 1 2 ( σ
(l) noise
)2 (( N∑ i=Hx+1 λ (l) i ) + (µ (l) Hx+1: )Tµ (l) Hx+1: )
+ (µ
(l) Hx+1: )TΨ (l) 1 (A (l))−1(Ψ (l) 1 ) Tµ (l) Hx+1: 2 ( σ
(l) noise )2 + log(|(A(l))−1|)2 −N̂
2 + N∑ i=1+Hx−Hh log(2πλ (l) i ) 2 − Hh∑ i=1+Hx−Hh log(2π) 2 −
( λ
(l) i +
( µ
(l) i )2) 2  − (L+ 1)N̂
2 log(2π)− L+1∑ l=1 KL(qZ(l) ||pZ(l))
def = LREVARBVSS-opt ,
where we derive LREVARBSS-opt by being not variational over Z(l).
The IP regularization case with A(l) = Ψ(l)2 + (σ (l) noise) 2K (l) MM and QREVARB, PREVARB just defined with the priors and variational distributions for Z(l) and H(l), is given by:
log(p(yHx+1:|X)) = . . . see Titsias & Lawrence (2010) until Equation 14, JI ≥ Ep a(l)|U(l) [ ∑L+1 l=1 exp(〈log(N (h(l)Hx+1:|Φ (l)(K (l) MM ) −1a(l), (σ (l) noise) 2IN ))〉QREVARB )]
− ( 2 ( σ(l)noise )2)−1∑L+1 l=1 Ψ (l) 0 − tr((K (l) MM ) −1Ψ(l)reg)−KL(QREVARB||PREVARB)
= −N̂ −M 2 L+1∑ l=1 log (( σ (l) noise )2) − yTHx+1:yHx+1: 2 ( σ(L+1)noise
)2 + log(|K(L+1)MM |)2 + yTHx+1:Ψ (L+1) 1 (A (L+1))−1(Ψ(L+1)1 ) TyHx+1:
2 ( σ(L+1)noise )2 + log(|(A(L+1))−1|)2 − Ψ (L+1) 0 + tr((K (L+1) MM ) −1Ψ (L+1) reg )
2 ( σ(L+1)noise )2 −
L∑ l=1 − 1 2 ( σ
(l) noise
)2 (( N∑ i=Hx+1 λ (l) i ) + (µ (l) Hx+1: )Tµ (l) Hx+1: ) + log(|K(l)MM |) 2
+ (µ
(l) Hx+1: )TΨ (l) 1 (A (l))−1(Ψ (l) 1 ) Tµ (l) Hx+1: 2 ( σ
(l) noise )2 + log(|(A(l))−1|)2 − Ψ (l) 0 + tr((K (l) MM ) −1Ψ (l) reg)
2 ( σ(l)noise
)2 − N̂2 + N∑
i=1+Hx−Hh
log(2πλ (l) i )
2
− Hh∑
i=1+Hx−Hh
log(2π)
2 −
( λ
(l) i +
( µ
(l) i )2) 2 − (L+ 1)N̂2 log(2π)− L+1∑ l=1 KL(qZ(l) ||pZ(l))
def = LREVARBVSS-IP ,
where again we derive LREVARBSS-IP by being not variational over Z(l).

heading: 6.3.5 PREDICTIONS
text: Predictions for each layer l and new ĥ(l)∗ are performed in the simple DRGP-(V)SS case with
Eq f∗(l)
[ f∗(l) ] = Ψ
(l) 1∗m (l),
Vq f∗(l)
[ f∗(l) ] = ( m(l) )T ( Ψ
(l) 2∗ −
( Ψ
(l) 1∗
)T Ψ
(l) 1∗
) m(l) + tr ( s(l)Ψ
(l) 2∗
) ,
where for the optimal distribution case for a(l) we have with A(l) = Ψ(l)2 + (σ (l) noise) 2IM
m(l) opt = ( A(l) )−1 ( Ψ
(l) 1
)T µ
(l) Hx+1: , s(l) opt = (σ (l) noise)
2 ( A(l) )−1 ,
for 1, . . . , L, and fully analog for l = L+ 1 by replacing µ(l)Hx+1: with yHx+1:. In the DRGP-(V)SS-IP case we make predictions for each layer l and new ĥ(l)∗ through
Eq f∗(l)
[ f∗(l) ] = Ψ
(l) 1∗Λ (l),
Vq f∗(l)
[ f∗(l) ] = ( Λ(l) )T ( Ψ
(l) 2∗ −
( Ψ
(l) 1∗
)T Ψ
(l) 1∗
) Λ(l)
+ Ψ0 − tr ( (K (l) MM ) −1Ψ (l) reg∗ − (σ(l)noise) 2 ( A(l) )−1 Ψ (l) 2∗ ) ,
where A(l) = Ψ(l)2 + (σ (l) noise) 2K (l) MM and
Λ(l) opt = ( A(l) )−1 ( Ψ
(l) 1
)T µ
(l) Hx+1: ,
for 1, . . . , L, and fully analog for l = L+ 1 by replacing µ(l)Hx+1: with yHx+1:.
Prediction for a new input in all cases has time complexity O((L + 1)M3), which comes from the iterative prediction through all GP layers and the calculation of the statistics, see Appendix 6.3.6.

heading: 6.3.6 DISTRIBUTED VARIATIONAL INFERENCE FOR DRGP-(V)SS(-IP)
text: We refer to (Gal et al., 2014), Equation (4.3), for a comparison. Calculating the optimal REVARB-(V)SS and REVARB-(V)SS-IP requires O(NM2Qmax(L + 1)), where Qmax = max
l=1...,L+1 Q(l), Q(l) def= dim(ĥi (l) ) and ĥ(l)i is coming from the Equations in (13) for
a fixed chosen i and l = 1, . . . , L+ 1. In this section we show how we can reduce the complexity of inference in the REVARB-(V)SS(-IP) setting with distributed inference to O(M3) if the number of cores scales suitably with the number of training-data. We show this for LREVARB(V)SS-IP , because this is the most complex bound, and all other bounds reduce to special cases of this. LREVARB(V)SS-IP in Appendix 6.3.4, separated for each hidden layer and the output layer (µHx+1: and λHx+1:
replaced by yHx+1:), can be written as LREVARB(V)SS-IP = L+1∑ l=1 Bl, where we have the KL terms
KL(qZ(l) ||pZ(l)) = 1
2 ∑M m=1 tr ( β(l)m + α (l) m (α (l) m ) T ) − 1 2 ∑M m=1 log(|β(l)m |)− MQ(l) 2 ,
for l = 1, . . . , L+ 1, the terms
−N̂ 2
+ ∑N
i=1+Hx−Hh
log(2πλ (l) i ) 2 − ∑Hh i=1+Hx−Hh log(2π) 2 −
( λ
(l) i +
( µ
(l) i )2) 2 ,
for l = 1, . . . , L and
Bl = − N̂ −M
2 log
(( σ
(l) noise
)2) − (∑N i=Hx+1 λ (l) i ) + µTHx+1:µHx+1:
2 ( σ
(l) noise )2 + log(|K(L+1)MM |)2 + µTHx+1:Ψ (l) 1 ( A(l) )−1 (Ψ (l) 1 ) TµHx+1:
2 ( σ
(l) noise
)2 + log(|(A(l))−1|)2 − Ψ (l) 0 + tr((K (l) MM ) −1Ψ (l) reg) 2 ( σ(l)noise )2 − N̂
2 log(2π), for all l = 1, . . . , L + 1 and which can be separated further into Bl = ∑N i=Hx+1 Bli, a sum of N̂ = N −Hx independent terms, extracting log(|(A(l)) −1|) 2 , log(|K(l)MM |) 2 and tr((K(l)MM ) −1Ψ(l)reg)
2(σ(l)noise) 2 , by
Bli = − 1− M N̂
2 log
(( σ
(l) noise
)2) −
( σ (l) power )2 + λ (l) i + µ (l) i µ (l) i
2 ( σ
(l) noise )2 + µ (l) i (Ψ (l) 1 )·i(A(l))−1(Ψ(l)1 )T·iµ (l) i
2 ( σ
(l) noise )2 + log(2πλ(l)i )2 − log(2π)2 , where (Ψ(l)1 )·i means, taking the i-th column of Ψ(l)1 for l = 1, . . . , L+ 1, i = Hx + 1 . . . , N . We further inspect
A(l) = Ψ (l) 2 + (σ (l) noise) 2K (l) MM = ∑N i=Hx+1 ( Ψi2 )(l) + (σ (l) noise) 2K (l) MM ,
and Ψ(l)reg = ∑N i=Hx+1 ( Ψireg )(l) .
These terms and the sums of Ψ(l)2 and Ψ (l) reg can be computed on different cores in parallel without communication. Only the 3(L+ 1) inversions and determinants of A(l) and K(l)MM now are responsible for the complexity, which can also be computed on 3(L+ 1) cores. Summing this bound over i and l, we obtain the total complexity of O(M3) per single iteration with N̂(L+ 1) cores.

Review summaries:
This paper addresses the problem of modeling sequential data based on one of the deep recurrent Gaussian process (DRGP) structures proposed by Mattos et al (2016). This structure acts like a recurrent neural net where every layer is defined as a GP. One of the main limitations of the original method proposed by Mattos et al (2016) is that it is limited to a small set of covariance functions, as the variational expectations over these have to be analytically tractable.