Deep recurrent Gaussian process with variational Sparse Spectrum approximation
Roman Föll Bernard Haasdonk Markus Hanselmann Holger Ulmer


Efficient charge modulation in ultrathin -field-effect transistors
    W.G. van der Wiel
    December 30, 2023
==================================================================




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




§ 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 <cit.> have proven to perform well on these tasks. They consist of observations, output-data and input-data, structured sequentially for shifted discrete time steps.

The general form of a recurrent model is given by

𝐡_i     = (𝐡_i-1,…,𝐡_i-H,𝐱_i-1,…,𝐱_i-H) + ϵ_i^𝐡,

𝐲_i     = (𝐡_i,…,𝐡_i-H)+ϵ_i^𝐲,

where 𝐱_i is an external input, 𝐲_i is an output observation, 𝐡_i is a latent hidden representation or state at time i=H+1,…,N, where N∈ℕ is the number of data samples, H∈ℕ is the chosen time horizon, , are non-linear functions modeling transition and observation and ϵ_i^𝐡, ϵ_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 <cit.>. The general recurrent model given in Equation (<ref>) and (<ref>) 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 (<ref>) in a deep manner <cit.>.

To achieve scalability, GPs normally make use of sparse approximations for the covariance function. This paper proposes DRGP models based on (V)SS approximations <cit.>, denoted by DRGP-(V)SS. Therefore, we follow the same deep recurrent structure as introduced in <cit.>. For reproducibility of the experimental results, we provide the code online[DRGP-(V)SS code available at <http://github.com/RomanFoell/DRGP-VSS>.]. To summarize, the contributions of this paper are the following:

	
  * Extension of the sparse GP based on the SS approximation and the improved VSS approximation 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 <cit.>;
	
  * 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 <cit.> 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 <cit.> in a simulation setting. For scalability, distributed variational inference (DVI) <cit.> is recommended and can lower the complexity from 𝒪(NM^2Q_max(L+1)) down to 𝒪(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.



§ RELATED WORK TO GAUSSIAN PROCESSES WITH SSM AND DGPS


An SSM with GPs (GP-SSM) for transition and observation functions is used by <cit.>, where the uncertainty in the latent states is not accounted for, which can lead to overconfidence. <cit.> solved this problem, but they have complicated approximate training and inference stages and the model is hard to scale. <cit.> used a GP for transition, while the observation is parametric. <cit.> 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 <cit.> 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. <cit.> 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 <ref>, <ref>, where the latent (non-observed) output-data is learned as a hidden state. We refer to the reports <cit.>[Available on the websites <https://arxiv.org/>, <https://openreview.net>.] for a detailed but preliminary formulation of the models and experiments presented in this paper.

Following <cit.>, a Deep Gaussian process (DGP) is a model assuming

    𝐲_i    = f^(L+1)(…(f^(1)(𝐱_i) + ϵ_i^𝐡^(1))…) + ϵ_i^𝐲,

where the index i=1,…,N is not necessarily the time and where we define 𝐡_i^(1)def= f^(1)(𝐱_i) + ϵ_i^𝐡^(1), 𝐡_i^(l+1)def= f^(l)(𝐡_i^(l)) + ϵ_i^𝐡^(l), for l=2…,L-1, where L∈ℕ is the number of hidden layers. The noise ϵ_i^𝐡^(l), ϵ_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. <cit.> introduced these kind of DGPs based on the sparse variational approximation following <cit.>. Based on this, <cit.> introduced a DGP with a variationally auto-encoded inference mechanism and which scales on larger data-sets. <cit.> introduced a DGP for the so called random Fourier features (RFF) approach <cit.>, 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,  <cit.> introduced an approximation framework for DGPs, which is similar to the single GP of <cit.>, 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 <cit.>, the RGP, which we call DRGP-Nyström, based on the variational sparse Nyström/inducing-point approximation introduced by <cit.>, as well as <cit.>, which we call GP-LSTM, based on deep kernels via a long-short term memory (LSTM) network <cit.>, 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 <cit.>.









§ GAUSSIAN PROCESSES (GPS) 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 <cit.> and its improved version VSSGP by <cit.>. Based on these, we derive new variational approximations. We use the notation a, f_𝐱, y, (italic) for random variables, 𝐚, f(𝐱), 𝐲 (upright) for realizations and data.



 §.§ Gaussian processes


A stochastic process f is a GP if and only if any finite collection of random variables f_𝐗def=[f_𝐱_1,…,f_𝐱_N]^T forms a Gaussian random vector <cit.>. A GP is completely defined by its mean function m: ℝ^Q→ℝ,𝐱↦ m(𝐱), Q∈ℕ the input-dimension, and covariance function k: ℝ^Q×ℝ^Q→ℝ,(𝐱,𝐱')↦ k(𝐱,𝐱') <cit.>, where

    m(𝐱)    def= 𝐄[f_𝐱],
    
    k(𝐱,𝐱')    def=𝐄[(f_𝐱-m(𝐱))(f_𝐱'-m(𝐱'))],

and the GP will be written as f∼𝒢𝒫(m,k). Be aware of that a valid covariance function must produce a positive definite matrix K_NNdef=(k(𝐱_i,𝐱_j))_i,j=1^N∈ℝ^N× N, when filling in combinations of data-input points 𝐱_i, i=1,…,N.

Let 𝐲def=[y_1,…,y_N]^T∈ℝ^N, 𝐗def=[𝐱_1,…,𝐱_N]^T∈ℝ^N× Q be our obeservations and we assume
y_i = f(𝐱_i) + ϵ_i^y, where  ϵ_i^y∼𝒩(0,_noise^2), for i=1,…,N, and our aim is to model any set of function values 𝐟def=[f(𝐱_1),…,f(𝐱_N)]^T∈ℝ^N at 𝐗 as samples from a random vector f_𝐗. Moreover, we assume the prior f_𝐗|X∼𝒩(0,K_NN), meaning that any set of function values 𝐟 given 𝐗 are jointly Gaussian distributed with mean 0∈ℝ^N  and covariance matrix K_NN.

The predictive distribution p_f_𝐱_∗|x_∗,X,y for a test point 𝐱_∗∈ℝ^Q, where K_N∗def= (k(𝐱_i,𝐱_∗))_i=1^N∈ℝ^N, and analogously K_∗∗, K_∗ N, can be derived through the joint probability model
and conditioning as

    f_𝐱_∗|x_∗,X,y∼𝒩(K_∗ N(K_NN + _noise^2 I_N)^-1𝐲,
    
    K_∗∗- K_∗ N(K_NN + _noise^2 I_N)^-1K_N∗).

In preview of our experiments in Section <ref> and the following ones, we choose a specific covariance function, the spectral mixture (SM) covariance function <cit.>

    k(𝐱,𝐱')def=
       _power^2 e^∑_q=1^Q(x_q-x_q')^2/-2l_q^2cos(∑_q=1^Q2π(x_q-x_q')/p_q),

with an amplitude _power^2 and length scales l_q, p_q∈ℝ, q=1,…, Q. As p_q →∞, this corresponds to the squared exponential (SE) covariance function in the limit <cit.>.



 §.§ Variational Sparse Spectrum GP


We introduce the SSGP following <cit.>. For a stationary covariance function k on ℝ^Q×ℝ^Q there exists a function :ℝ^Q→ℝ,τ↦(τ), such that k(𝐱,𝐱') =(𝐱-𝐱') for 𝐱,𝐱'∈ℝ^Q. Bochner’s theorem states that any stationary covariance function k can be represented as the Fourier transform of a positive finite measure  <cit.>. Then (τ), using τ=𝐱-𝐱', can be expressed via Monte Carlo approximation (MCA) following <cit.>, Section 2, as

    (τ)    = ∫_ℝ^Q e^2π i 𝐳^Tτd(𝐳)
       ≈_power^2/M∑_m=1^M2cos(2π𝐳_m^T(𝐱-𝐮_m)+b_m)
       cos(2π𝐳_m^T(𝐱'-𝐮_m)+b_m)def=k̃(𝐱,𝐱').

We refer to 𝐳_m as the spectral points, b_m as the spectral phases and 𝐮_m as the pseudo-input/inducing points for m=1,…, M. Choosing the probability density like in <cit.>, Proposition 2, we approximate the SM covariance function with a scaling matrix 𝔏def= diag([2πl_q]_q=1^Q)∈ℝ^Q× Q, a scaling vector 𝐩def= [p_1^-1,…,p_Q^-1]^T∈ℝ^Q, 𝐙def=[𝐳_1,…,𝐳_M]^T∈ℝ^Q× M, 𝐔def=[𝐮_1,…,𝐮_M]^T∈ℝ^Q× M via

    ϕ(𝐱_n) def =√(2_power^2)/M[cos(𝐱_n,m)]_m=1^M∈ℝ^M.

Here it is 𝐱_n,m =2π(𝔏^-1𝐳_m+𝐩)^T(𝐱_n-𝐮_m) + b_m, we sample b∼Unif[0,2π], z∼𝒩(0,I_Q), choose 𝐮_m as subset of the input-data and we set K̃_NN^(SM)def= ΦΦ^T with Φdef= [ϕ(𝐱_n)]_n=1^N∈ℝ^N× M. 

In <cit.> 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 <cit.>, Section 4, for the 1-dimensional output case 𝐲∈ℝ^N× 1. By replacing the covariance function with the sparse covariance function k̃^(SM) and setting the priors to

    p_Zdef= ∏_m=1^M p_z_m, where z_m ∼𝒩(0,I_Q),
       p_a, where a∼𝒩(0,I_M),

for m=1,…,M, where we have y|a,Z,U,X∼𝒩(Φ𝐚,_noise^2 I_N) (we do not define priors on 𝐩, 𝔏^-1, 𝐔, 𝐛def=[b_1,…,b_M]^T∈ℝ^M), we can expand the marginal likelihood (ML) to

    p(𝐲|𝐗) = ∫ p(𝐲|𝐚,𝐙,𝐔,𝐗)p(𝐚)p(𝐙)d𝐚 d𝐙,

highlighting 𝐔 just in the integral, to be notationally conform to <cit.>, Section 3.
 
Now, to improve the SSGP to VSSGP, variational distributions are introduced in terms of

    q_Zdef= ∏_m=1^M q_z_m, where z_m∼𝒩(_m,_m),
       q_a, where a∼𝒩(𝐦,𝐬),

with _m∈ℝ^Q× Q diagonal, for m=1,…,M, and 𝐬∈ℝ^M× M diagonal.
From here on we use variational mean-field approximation to derive the approximate models with different lower bounds to the ln-ML introduced by <cit.>:

    ln(p(𝐲|𝐗)) ≥    𝐄_q_aq_Z[ln(p(𝐲|𝐚,𝐙,𝐔,𝐗))]
        - 𝐊𝐋(q_aq_Z||p_ap_Z).

As usual, 𝐊𝐋 defines the Kullback-Leibler (KL) divergence. By proving the existence of an optimal distribution q^opt_a for a <cit.>, Proposition 3, in the sense of a functional local optimum of the right hand side of (<ref>), where  a∼𝒩(A^-1Ψ_1^T𝐲,_noise^2 A^-1), with A = Ψ_2 + _noise^2 I_M, Ψ_1 = 𝐄_q_Z[Φ]∈ℝ^N× M, Ψ_2 = 𝐄_q_Z[Φ^TΦ]∈ℝ^M× M, we can derive the optimal bound case. 



 §.§ (V)SSGP with regularization properties via inducing points (IP)


As a first contribution of this paper we combine two approximation schemes to four new methods (V)SSGP-IP-1/2. We want to point out that the (V)SSGP does not have the same regularization properties as the GP of <cit.>, when optimizing the parameters 𝐔, because the priors in (<ref>), (<ref>) of the weights 𝐚 is defined generically via <cit.>, equations (2.113) - (2.115). These parameters 𝐔, following <cit.>, are similar to the sparse pseudo-input approximation <cit.>, but in the lower bound in (<ref>) they are simply used without being linked to the weights 𝐚.

We now define them as

    p_a|U, where a|U∼𝒩(0,K_MM),

    p_y|a, Z,U,X, where y|a, Z,U,X∼

   𝒩(𝐦_1,2,K_NN-K_NMK_MM^-1K_MN + _noise^2 I_N),
where K_NN, K_NM, K_MM are defined through the given covariance function in Equation (<ref>) and 𝐦_1=Φ𝐚, 𝐦_2=Φ K_MM^-1𝐚.

We can show that for these definitions the integral in Equation (<ref>) can be marginalized straightforward for the weights 𝐚. We then obtain that our data-samples 𝐲 are coming from a Gaussian distribution with mean 0 and the true covariance matrix K_NN, plus the discrepancy of the two well-known sparse covariance matrices for the Sparse Spectrum and the Nyström case, plus the noise assumption:

    K_NN+(ΦK_MM^-1Φ^T-K_NMK_MM^-1K_MN) + _noise^2 I_N.

For 𝐦_1 we have ΦΦ^T and we highlight this case by … in the following throughout the paper, for 𝐦_2 we have Φ K_MM^-1Φ^T, and K_NMK_MM^-1K_MN is the Nyström case.

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 <cit.>, Section 3.1, the optimal variational distribution for a collapses by reversing Jensen's inequality and is similar to the one obtained in <cit.>. The resulting bounds, here for VSSGP-IP-1/2, can be calculated in the same way as <cit.> until Equation (14) in closed form:  

   ln(p(𝐲|𝐗)) ≥-(N-M)/2ln(_noise^2)- N/2ln(2π)

   -𝐲^T𝐲/2_noise^2+ln(|K_MM||A^-1|)/2+ 𝐲^TΨ_1 A^-1Ψ_1^T𝐲/2_noise^2

   -tr(K_NN-K_MM^-1K_MNK_NM)/2_noise^2- 𝐊𝐋(q_Z||p_Z).
SSGP-IP-1/2 is derived by deleting 𝐊𝐋(q_Z||p_Z) and setting Ψ_1= Φ, Ψ_2 = Φ^TΦ.
Consequently, the resulting bound in (<ref>) has an extra regularization property compared to the right hand side of (<ref>), which is reflected in the different form of A = Ψ_2 + _noise^2 K_MM for 𝐦_2, which involves K_MM, the chosen covariance matrix filled in with the pseudo-input points 𝐔, two extra terms for 𝐦_1 and three extra terms for 𝐦_2: 

    ln(|K_MM|)/2,  tr(K_NN)/2_noise^2,  tr(K_MM^-1K_MNK_NM)/2_noise^2.




 §.§ Variational approximation of the input-space for (V)SSGP(-IP-1/2)


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 (<ref>) and (<ref>) these expressions even exist. To prevent misunderstanding, we will write from now on 𝐇 = [𝐡_1,…,𝐡_N]^T∈ℝ^N× Q instead of 𝐗, to highlight that 𝐇 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= ∏_i=1^N p_h_i, where h_i∼𝒩(0,I_Q),

q_H    def= ∏_i=1^N q_h_i, where h_i ∼𝒩(_i,_i),

_i∈ℝ^Q× Q diagonal, for i=1,…,N.

As a consequence, for VSSGP-IP-1/2 we overall derive statistics Ψ_0 = tr(𝐄_q_H[K_NN])= N _power^2, Ψ_1 = 𝐄_q_Zq_H[Φ]∈ℝ^N× M, Ψ_2 = 𝐄_q_Zq_H[Φ^TΦ]∈ℝ^M× M and Ψ_reg = 𝐄_q_H[K_MNK_NM] as defined in the Appendix <ref>. These statistics are essentially the given matrices Φ, Φ^TΦ, K_MNK_NM from the beginning, but every input 𝐡_i and every spectral point 𝐳_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-1/2 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 (V)SSGP-IP-1/2 cases we get the lower bound:

   ln(p(𝐲|𝐇)) ≥-(N-M)/2ln(_noise^2)- N/2ln(2π)

   - 𝐲^T𝐲/2_noise^2 +ln(|K_MM||A^-1|)/2+ 𝐲^TΨ_1 A^-1Ψ_1^T𝐲/2_noise^2

   -Ψ_0/2_noise^2+tr(K_MM^-1Ψ_reg)/2_noise^2- 𝐊𝐋(q_Zq_H||p_Zp_H).
For the extension of the SSGP and VSSGP in the optimal bound case Section <ref>, Equation (<ref>), we again have A = Ψ_2 + _noise^2 I_M and eliminate ln(|K_MM|)/2, Ψ_0/2_noise^2, tr(K_MM^-1Ψ_reg)/2_noise^2 in the lower bound (<ref>). We only focus on the optimal bound cases in the following.


§ DRGP WITH VARIATIONAL SPARSE SPECTRUM APPROXIMATION


In this section, we want to combine our newly derived (V)SS(-IP-1/2) approximations in the sections <ref>, <ref>, overall resulting in six GP cases: SSGP, VSSGP, SSGP-IP-1, VSSGP-IP-1, SSGP-IP-2, VSSGP-IP-2, with the framework introduced in <cit.>, to derive our DRGP models: DRGP-SS, DRGP-VSS, DRGP-SS-IP-1, DRGP-VSS-IP-1, DRGP-SS-IP-2, DRGP-VSS-IP-2.



 §.§ DRGP-(V)SS(-IP-1/2) model definition


Choosing the same recurrent structure as in <cit.>, where now i represents the time, we have

    :      h_i^(l) = f^(l)(𝐡̂_i^(l))+ϵ_i^h^(l),  with prior
       f^(l)_𝐇̂^(l)|Ĥ^(l)∼𝒩(0,K_N̂N̂^(l)),   l = 1,…,L
    :      y_i = f^(l)(𝐡̂_i^(l))+ϵ_i^y,  with prior
       f^(l)_𝐇̂^(l)|Ĥ^(l)∼𝒩(0,K_N̂N̂^(l)),   l = L+1,

with ,  in Equation (<ref>) and (<ref>), ϵ_i^h^(l)∼𝒩(0,(_noise^(l))^2), ϵ_i^y∼𝒩(0,(_noise^(L+1))^2) and N̂=N-H_𝐱, for i = H_𝐱 + 1,…,N. The matrix K_N̂N̂^(l) represents a covariance matrix coming from our chosen k̃ in Equation (<ref>) and (<ref>). A set of input-data 𝐇̂^(l)=[𝐡̂_H_𝐱+1^(l),…,𝐡̂_N^(l)]^T is specified as
 

    𝐡̂_i^(1)def=[[h_i-1^(1),…,h_i-H_h^(1)],[𝐱_i-1,…,𝐱_i-H_𝐱]]^T,
       𝐡̂_i^(l)def=[[h_i-1^(l),…,h_i-H_h^(l)],[h_i^(l-1),…,h_i-H_h + 1^(l-1)]]^T,



 𝐡̂_i^(L+1) def=[h_i^(L),…,h_i-H_h + 1^(L)]^T,      

where 𝐡̂_i^(1)∈ℝ^H_h+H_𝐱Q, 𝐡̂_i^(l)∈ℝ^2H_h for l = 2,…,L, 𝐡̂_i^(L+1)∈ℝ^H_h, for i = H_𝐱 + 1,…,N. For simplification we set Hdef=H_𝐱=H_h in our experiments. We further introduce the notation 𝐲_H_𝐱+1:def= [y_H_𝐱 + 1,…,y_N]^T∈ℝ^N̂, starting the output-data vector 𝐲 from index H_𝐱+1.

Now we use the new approximations in the sections <ref>,  <ref>, to derive first, for the setting (V)SSGP, the new joint probability density

    p_y_H_𝐱+1:[a^(l),Z^(l),h^(l),U^(l)]_l=1^L+1|X=
       ∏_l=1^L+1p_h^(l)_H_𝐱+1:|a^(l),Z^(l),Ĥ^(l),U^(l)p_a^(l)p_Z^(l)p_h̃^(l),

with 𝐡^(L+1)_H_𝐱+1:def= 𝐲_H_𝐱+1:, p_h̃^(L+1)def=1, h^(L+1)def={} (these definitions are just for simplification of the notation for the last GP layer) and 𝐡̃^(l)def= [h_1+H_𝐱-H_h^(l),…,h_H_h^(l)]^T∈ℝ^2H_h-H_𝐱 for l=1,…,L. Here the priors are similar to (<ref>), (<ref>), (<ref>) with

    p_a^(l), where a^(l)∼𝒩(0,I_M),
       p_Z^(l)def= ∏_m=1^M p_z^(l)_m, where z_m^(l)∼𝒩(0,I_Q),
       p_h̃^(l), where h̃^(l)∼𝒩(0,I_2H_h-H_𝐱),

for m=1…,M, l=1,…,L+1, excluding p_h̃^(L+1), and the product of them is defined as P_.8REVARB. The variational distributions are similar to (<ref>), (<ref>), (<ref>)

    q_a^(l), where a^(l)∼𝒩(𝐦^(l),𝐬^(l)),
       q_Z^(l)def= ∏_m=1^M q_z^(l)_m, where z_m^(l)∼𝒩(_m^(l),_m^(l)),
        q_h^(l)def=∏_i=1+H_h-H_𝐱^N q_h^(l)_i, where  h_i^(l)∼𝒩(^(l)_i,^(l)_i),

where _m^(l)∈ ℝ^Q× Q is diagonal, for i=1+H_𝐱-H_h,…,N, m=1…,M, l=1,…,L+1, excluding q_h^(L+1), and the product of them is defined as Q_.8REVARB.

For the setting (V)SSGP-IP-1/2 we choose no variational distribution for a^(l), but, similar to the assumptions in (<ref>), (<ref>) for l = 1,…,L+1, the prior assumptions

    p_a^(l)|U^(l), where a^(l)|U^(l)∼𝒩(0,K_MM^(l)),
       p_h_H_𝐱+1:^(l)|a^(l),Z^(l),Ĥ^(l),U^(l), where 
       h_H_𝐱+1:^(l)|a^(l),Z^(l),Ĥ^(l),U^(l)∼𝒩(𝐦_1,2^(l),
          K_N̂N̂^(l)-K_N̂M^(l)(K_MM^(l))^-1K_MN̂^(l)+(_noise^(l))^2 I_N̂),

where 𝐦_1,2^(l)=Φ^(l)(K_MM^(l))^-1𝐚^(l).

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


 §.§ DRGP-(V)SS(-IP-1/2) evidence lower bound (ELBO)


Using standard variational approximation techniques <cit.>, the recurrent variational Bayes lower bound for the (V)SS approximations, denoted as REVARB-(V)SS, is given by

   ln(p(𝐲_H_𝐱+1:|𝐗))≥𝐄_Q_.8REVARB[∑_l=1^L+1.

   .ln(𝒩(𝐡^(l)_H_𝐱+1:|Φ^(l)𝐚^(l),(_noise^(l))^2 I_N̂))]

    - 𝐊𝐋(Q_.8REVARB||P_.8REVARB)def=L_.8(V)SS^.8REVARB,
and for the (V)SS-IP-1, (V)SS-IP-2 approximations, denoted as REVARB-(V)SS-IP-1, REVARB-(V)SS-IP-2 is given by

   ln(p(𝐲_H_𝐱+1:|𝐗))≥ln(∏_l=1^L+1𝐄_p_a^(l)|U^(l)[ .

   .exp(⟨ln(𝒩(𝐡^(l)_H_𝐱+1:|𝐦_1,2^(l),(_noise^(l))^2 I_N̂))⟩_Q_.8REVARB)])

   -∑_l=1^L+1Ψ_0^(l)/2(_noise^(l))^2-tr((K_MM^(l))^-1Ψ_reg^(l))/2(_noise^(l))^2

   - 𝐊𝐋(Q_.8REVARB||P_.8REVARB)def=L_.8(V)SS-IP-1/2^.8REVARB,
where ⟨·,·⟩ means the expectation under the integral.

Additionally, for the approximations in (<ref>), (<ref>) the optimal bound L^.8REVARB_.8(V)SS(-IP-1/2)-opt can be obtained immediately, analogously to <cit.>, Proposition 3, <cit.>, Section 3.1, and by the fact that the bound decomposes into a sum of independent terms for 𝐚^(l). Maximizing the lower bounds is equivalent to minimizing the KL-divergence of Q_.8REVARB 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^.8REVARB_.8(V)SS(-IP-1/2)-opt requires 𝒪(N̂M^2Q_max(L+1)), where Q_max = max_l=1…,L+1Q^(l), Q^(l)def=(𝐡̂_i^(l)) and 𝐡̂_i^(l) from the equations in (<ref>). DVI can reduce the complexity to 𝒪(M^3) if the number of cores scales suitably with the number of training-data, see Appendix <ref> for a detailed description. A detailed derivation of the REVARB-(V)SS(-IP-1/2) lower bounds can be found in the Appendix <ref>.



 §.§ Model predictions


After model optimization based on the lower bounds in the Equation (<ref>) and (<ref>), model predictions for new 𝐡̂_*^(l) in the REVARB-(V)SS(-IP-1/2) framework can be obtained based on the approximate, variational
posterior distribution Q_.8REVARB. 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 <ref>.



§ EXPERIMENTS


In this section we want to compare our methods DRGP-SS, DRGP-VSS, DRGP-SS-IP-1, DRGP-VSS-IP-1, DRGP-SS-IP-2, DRGP-VSS-IP-2, (optimal bound cases) against other well known sparse GPs with NARX structure and the full GP with NARX structure, the DRGP-Nyström <cit.>, the GP-LSTM <cit.>, the LSTM <cit.>, a simple RNN, the DGP-DS <cit.> and the DGP-RFF <cit.>, both with NARX structure for the first layer, the GP-SSM <cit.> and the PR-SSM <cit.>. The full GP is named GP-full, the FITC approximation <cit.> is named GP-FITC, the DTC approximation <cit.> is named GP-DTC, the SSGP <cit.> is named GP-SS, the VSSGP <cit.> 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 ones) are used to predict next output observations. To enable a fair comparison, all methods are trained on the same amount of data.




 §.§ Model implementation, configuration, data-sets


Our methods DRGP-(V)SS(-IP-1/2) were implemented in Python, using the lib 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 used the published code[GP-VSS code available at <https://github.com/yaringal/VSSGP>.][GP-LSTM, LSTM, RNN code available at <https://github.com/alshedivat/keras-gp>.][DGP-DS code available at <https://github.com/ICL-SML/Doubly-Stochastic-DGP>.][DGP-RFF code available at <https://github.com/mauriziofilippone/deep_gp_random_features>.] for GP-VSS, GP-LSTM, LSTM, RNN, DGP-RFF and DGP-DS. 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-SSM, PR-SSM we show the results from their papers. Details about the methods, their configuration, as well as the benchmark data-sets can be found in the Appendix <ref>, <ref>.



 §.§ Model learning and comparison


In Figure <ref> 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, Actuator, for the model DRGP-VSS. We initialize the states with the output training-data for all layers with minor noise (first row) and after training we obtain a trained state (second row). Unlike <cit.> Figure 2, (l), we get good variance predictions for all data-sets. We used our own implementation for <cit.>, 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 <ref>. 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 <ref>, <ref>. 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 <ref> with H=10. We see, that throughout the different models the same time horizon is favored in terms of robustness. In Figure <ref> the RMSE results for different layers L on the data-sets Drive, Actuator, Damper are shown. We can observe that different layers L are favored. 














 §.§ Conclusion


In this paper we introduced six new DRGPs based on the SS approximation introduced by <cit.> and the improved VSS approximation by <cit.>. We combined the (V)SS approximations with the variational inducing-point (IP) approximation from <cit.>, 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 L^.8REVARB_.8(V)SS(-IP-1/2)-opt, called REVARB-(V)SS(-IP-1/2). We could show that our methods slightly improve on the data-sets used in this paper compared to the RGP from <cit.> and other state-of-the-art methods, where moreover our sparse approximations are also practical for dimensionality reduction as shown in <cit.> and can be further expanded to a deep version in this application <cit.>. Furthermore, <cit.> introduced a generalized version of the (V)SS approximation, which should be adaptable for our case.



  §.§.§ Acknowledgments


We would like to acknowledge support for this project
from the ETAS GmbH.
We would like to acknowledge support for this project from the ETAS GmbH. We are financed for a PhD in cooperation with the University of Stuttgart Germany. We thank Gabriele Santin for reviewing our manuscript.


icml2019










§ APPENDIX


This additional material provides details about the derivations and configurations of the proposed DRGP-(V)SS(-IP-1/2) in the sections <ref>-<ref>. and elaborates on the methods and the
employed data-sets in the experiments in Section <ref>.



 §.§ 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 <cit.>. 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 Dryer[Received from <http://homes.esat.kuleuven.be/ tokka/daisydata.html>.] 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 Ballbeamdryer[Description can be found under <http://forums.ni.com/t5/NI-myRIO/myBall-Beam-Classic-Control-Experiment/ta-p/3498079>.] 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 <cit.> 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 <cit.>, 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 Load[Originally received from Global Energy Forecasting Kaggle competitions organized in 2012.], used in <cit.>, 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 data-point, 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 Emission[Available from <http://github.com/RomanFoell/DRGP-VSS>.] 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 <ref>. 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.


 §.§ Nonlinear system identification




The methods for the system identification tasks and their references are summarized in Table <ref>.

For the data-sets Drive and Actuator we chose for our methods DRGP-(V)SS(-IP-1/2) the setting L=2 hidden layers, M=100 spectral-points and time-horizon H_h = H_𝐱 = 10, which was also used in <cit.>, <cit.> and <cit.> for free simulation (using M pseudo-input points instead of spectral-points). For these two data-sets we filled the results from <cit.> into Table <ref>. 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_𝐱 = 10. For the data-set Damper we chose L=2, M=125 and H_h = H_𝐱 = 10. For the data-set Power Load we chose L=1, M=125 and H_h = H_𝐱 = 12. For the data-set Emission we chose L=1, M=125 and H_h = H_𝐱 = 10.

The other sparse GP, the full GP, DGP-RFF and DGP-DS were trained with NARX structure H_𝐱 = 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.

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_𝐱 = 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_𝐱 = 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 <cit.> 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 𝒩(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_q^(l) = ∞, for all q, l. So all GPs with (V)SS approximations were also initialized as SE covariance function, see Equation <ref>.

For all methods we used automatic relevance determination, so each input dimension has its own length scale. For our methods and DRGP-Nyström the noise parameters and the hyperparameters were initialized by _noise^(l)=0.1, _power^(l)=1 and the length scales by either l_q^(l)= √(max(𝐇̂_q^(l)) - min(𝐇̂_q^(l))) or l_q^(l)= max(𝐇̂_q^(l)) - min(𝐇̂_q^(l)), for all q, l, where 𝐇̂^(l)_q is the data-vector containing the q-th input-dimension values of every input-data point 𝐡̂_i^(l), for all i. Furthermore, we initialized the latent hidden GP states with the output-observation data 𝐲.

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 l_q=1 for all input-dimensions.

For DGP-RFF we used the initialization coming from the implementation with _power^(l)=1 and for the length scale initialization with l_q^(l)=0.5ln(Q^(l)) for all l, q. To our knowledge DGP-RFF has no parameter _noise^(l) for all l.

For DGP-DS we used the initialization _noise^(l)=0.1 and _power^(l)=1 and for the length scale initialization we chose l_q^(l)=1 for all l, q.

For all our implementations, see Section <ref>, we used the positive transformation x'=ln(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↦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 <ref>.

For our methods and DRGP-Nyström we fixed _noise^(l), _power^(l) for all l (optional the spectral points/pseudo-input points for DRGP-(V)SS; for the IP cases we never excluded the pseudo-input points because we would fall back to the DRGP-(V)SS case; for DRGP-Nyströ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-1/2) we fixed _m^(l) for all m, l to small value around 0.001, as well as the spectral phases b_m for all m, l sampling from Unif[0,2π] (this seems to work better in practice). The limitations for _m^(l) also holds for GP-VSS as well.

We want to remark at this point that setting 𝐮_m=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 <cit.>, 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 <ref>.

For PR-SSM we show the results from their paper in Table <ref>.

We show additional results for the data-sets Drive, Actuator and Emission with missing auto-regressive part for the first layer for our methods DRGP-(V)SS(-IP-1/2) and DRGP-Nyström in Table <ref>, 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  <ref>, 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.



 §.§ Variational approximation and inference for DRGP-(V)SS(-IP-1/2)


In the sections <ref>-<ref> we derive the statistics Ψ_0, Ψ_1, Ψ_2 and Ψ_reg for the six model cases DRGP-(V)SS(-IP-1/2). In Section <ref> we derive the resulting variational lower bounds. In Section <ref> 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, <cit.>,
	
  * B.2, <cit.>,
	
  * 1/2(cos(a-b+x-y)+cos(a+b+x+y))=cos(x+a)cos(y+b).
	
  * JI, for Jensen’s inequality.




  §.§.§ DRGP-VSS(-IP-1/2), the statistics , 


For the versions DRGP-VSS(-IP-1/2) the statistics are

    (Ψ_1)_nm    = 𝐄_q_z_mq_h_n[Φ_nm]
       B.1= 𝐄_q_h_n[√(2_power^2(M)^-1)e^-1/2𝐡̅_nm^T_m𝐡̅_nmcos(_m^T(𝐡_n-𝐮_m)+b_m)]
       B.2= Σ_m^1 Z_nme^-1/2_m^T C_nm_mcos(_m^T(𝐜_nm-𝐮_m)+b_m),

for m, = 1,…, M, n = 1,…, N with

    𝐡̅_nm   = 2π𝔏^-1(𝐡_n-𝐮_m),
    _m    = 2π(𝔏^-1_m+𝐩),
    𝐜_nm   = C_nm(_m (2π)^2𝔏^-2𝐮_m + _n^-1_n),
    C_nm   = (_m(2π)^2𝔏^-2 + _n^-1)^-1,
    𝐯_nm   = 𝐮_m-_n,
    V_nm   = (2π)^-2𝔏^2_m^-1  + _n,
    Z_nm   = 1/√(|V_nm|)e^-1/2𝐯_nm^T V_nm^-1𝐯_nm,
    Σ_m^1     = √(2_power^2∏_q=1^Q(l_q^2/_m_q)(M)^-1),

and



Ψ_2 = ∑_n=1^N (Ψ_2)^n, where

    (Ψ_2)_mm'^n    = 𝐄_q_z_mq_z_m'q_h_n[Φ^T_nmΦ_nm']
       B.1= 𝐄_q_h_n[2_power^2(M)^-1e^-1/2(𝐡̅_nm^T_m𝐡̅_nm+𝐡̅_nm'^T_m'𝐡̅_nm').
       .cos(_m^T(𝐡_n-𝐮_m)+b_m)cos(_m'^T(𝐡_n-𝐮_m')+b_m')]
       B.2= Σ_mm'^2Z_mm'^n (e^-1/2_mm'^T D_mm'^n_mm'cos(_mm'^T𝐝_mm'^n - _mm' + b̅_mm').
       . + e^-1/2+_mm'^T D_mm'^n+_mm'cos(+_mm'^T𝐝_mm'^n - +_mm' +  +b_mm')).

for m,m' = 1,…, M, m≠ m', with 

    b_mm'   =b_m-b_m',
    +b_mm'   =b_m+b_m',
    _mm'   =_m^T 𝐮_m- _m'^T 𝐮_m',
    +_mm'   =_m^T 𝐮_m+ _m'^T 𝐮_m',
    _mm'   = _m - _m',
    +_mm'   = _m + _m',
    𝐛_mm'   = B_mm'(2π)^2𝔏^-2(_m𝐮_m+_m'𝐮_m'),
    _mm'   = _m + _m',
    B_mm'   = (2π)^-2𝔏^2_mm'^-1,
    𝐝_mm'^n    = D_mm'^n(B_mm'^-1𝐛_mm'+_n^-1_n),
    D_mm'^n    = (B_mm'^-1 + _n^-1)^-1,
    𝐰_mm'^n    = 𝐛_mm'-_n,
    W_mm'^n    = B_mm'  + _n,
    𝐮_mm'   = 𝐮_m-𝐮_m',
    𝐔_mm'   = (2π)^-2𝔏^2(_m^-1  + _m'^-1),
    Z_mm'^n    = 1/√(|W_mm'^n||𝐔_mm'|)e^-1/2(𝐰_mm'^n^TW_mm'^n^-1𝐰_mm'^n+𝐮_mm'^T𝐔_mm'^-1𝐮_mm'),
    Σ_mm'^2    = _power^2∏_q=1^Q(l_q^2/√(_m_q_m'_q))(M)^-1,

and

    (Ψ_2)_mm^n    = 𝐄_q_z_mq_h_n[Φ^T_nmΦ_nm]
       B.1= 𝐄_q_h_n[2_power^2(M)^-1(1/2+1/2e^-2𝐡̅_nm^T_m𝐡̅_nmcos(2(_m^T(𝐡_n-𝐮_m)+b_m))]
       B.2= _power^2(M)^-1(1+Σ̃_m^2 Z̃_nme^-2_m^TC̃_nm_mcos(2(_m^T(𝐜_nm-𝐮_m)+b_m))),

for m,m' = 1,…, M, m=m', with

    Σ̃_m^2   = √(∏_q=1^Q(l_q^2/_m_q))2^-Q.




  §.§.§ DRGP-SS(-IP-1/2), the statistics , 

For the versions DRGP-SS(-IP-1/2) the statistics are

    (Ψ_1)_nm    = 𝐄_q_h_n[Φ_nm]
        =  𝐄_q_h_n[√(2_power^2(M)^-1)cos(𝐳̂_m^T(𝐡_n-𝐮_m)+b_m)]
       B.2= √(2_power^2(M)^-1) e^-1/2𝐳̂_m^T_n𝐳̂_mcos(𝐳̂_m^T(_n-𝐮_m)+b_m),

for m = 1,…, M, n = 1,…, N with 

    𝐳̂_m = 2π(𝔏^-1𝐳_m+𝐩),

and



Ψ_2 = ∑_n=1^N (Ψ_2)^n, where

    (Ψ_2)_mm'^n    = 𝐄_q_h_n[Φ^T_nmΦ_nm']
       B.1= 𝐄_q_h_n[2_power^2(M)^-1cos(𝐳̂_m^T(𝐡_n-𝐮_m)+b_m)cos(𝐳̂_m'^T(𝐡_n-𝐮_m')+b_m')]
        = _power^2(M)^-1(e^-1/2𝐳̅_mm'^T_n𝐳̅_mm'cos(𝐳̅_mm'^T_n - _mm' + b̅_mm').
       .+e^-1/2+𝐳_mm'^T_n+𝐳_mm'cos(+𝐳_mm'^T_n - +_mm' + +b_mm'))),

for m,m' = 1,…, M with 

    _mm'   =𝐳̂_m^T 𝐮_m- 𝐳̂_m'^T 𝐮_m',
    +_mm'   =𝐳̂_m^T 𝐮_m+ 𝐳̂_m'^T 𝐮_m',
    𝐳̅_mm'   = 𝐳̂_m - 𝐳̂_m',
    +𝐳_mm'   = 𝐳̂_m + 𝐳̂_m',

for the other variables, see the defined variables in the DRGP-VSS case.

Further we have for the covariance matrix K̃_NN=ΦΦ^T the following

    𝐄_q_h_nq_h_n'[Φ_mnΦ_mn'^T]
       B.1= 𝐄_q_h_nq_h_n'[2_power^2(M)^-1cos(𝐳̂_m^T(𝐡_n-𝐮_m)+b_m)cos(𝐳̂_m^T(𝐡_n'-𝐮_m)+b_m)]
       =Ψ_1Ψ_1^T

for n, n'=1,…,N, n ≠ n' and

    =2_power^2(M)^-1(1/2+1/2e^-2𝐳̂_m^T_n𝐳̂_mcos(2(𝐳̂_m^T(_n-𝐮_m)+b_m))),

n, n'=1,…,N, n= n'.


  §.§.§ DRGP-(V)SS-IP-1/2, the statistics  and 



For Ψ_reg = 𝐄_q_H[K_MNK_NM] see Ψ_2 in Section <ref> but setting b_m=0, _m=0 and _m=1 for all m =1,…,M.


Ψ_0 naturally is given by Ψ_0 = tr(𝐄_q_H[K_NN])= N _power^2 because of the chosen SE covariance function.


  §.§.§ DRGP-(V)SS(-IP-1/2), lower bounds


In this section we derive the different variational lower bounds for our models DRGP-(V)SS(-IP-1/2). We first show the bound L^.8REVARB_.8(V)SS without optimal variational distribution for a^(l). Then the bounds L^.8REVARB_.8(V)SS-opt with optimal variational distribution for a^(l) follows, as well as the L^.8REVARB_.8(V)SS-IP-1/2-opt case.

We use the simplified notation d𝐀 d𝐙 d𝐇 = d𝐚^(1)… d𝐚^(L+1) d𝐙^(1)… d𝐙^(L+1) d𝐡^(1)…  d𝐡^(L).

    ln(p(𝐲_H_𝐱+1:|𝐗))
        = ln(∫ p(𝐲_H_𝐱+1:,[𝐚^(l),Z^(l),𝐡^(l),U^(l)]_l=1^L+1|𝐗)d𝐀 d𝐙 d𝐇)
        = ln(∫Q_.8REVARB/Q_.8REVARBp(𝐲_H_𝐱+1:,𝐚^(L+1),Z^(L+1),U^(L+1),[𝐚^(l),Z^(l),𝐡^(l),U^(l)]_l=1^L|𝐗)d𝐀 d𝐙 d𝐇)
       JI≥∫ Q_.8REVARBln(p(𝐲_H_𝐱+1:,𝐚^(L+1),Z^(L+1),U^(L+1),[𝐚^(l),Z^(l),𝐡^(l),U^(l)]_l=1^L|𝐗)/Q_.8REVARB)d𝐀 d𝐙 d𝐇
       = ∑_l=1^L+1∫ q(𝐚^(l))q(Z^(l))q(𝐡^(l))p(𝐡^(l)_H_𝐱+1:|𝐚^(l),Z^(l),𝐇̂^(l),U^(l))d𝐚^(l)d𝐙^(l)d𝐡^(l)
       -∑_l=1^L(N̂/2-∑_i=1+H_𝐱-H_h^Nln(2π_i^(l))/2+∑_i=1+H_𝐱-H_h^H_hln(2π)/2 + (_i^(l) + (_i^(l))^2)/2) 
       - ∑_l=1^L+1𝐊𝐋(q_a^(l)||p_a^(l)) - 𝐊𝐋(q_Z^(l)||p_Z^(l))
        = -N̂/2∑_l=1^L+1ln(2π(_noise^(l))^2)-𝐲_H_𝐱+1:^T𝐲_H_𝐱+1:/2(_noise^(L+1))^2
        + 𝐲^T_H_𝐱+1:Ψ_1^(L+1)𝐦^(L+1)/(_noise^(L+1))^2 - tr(Ψ_2^(L+1)(𝐬^(L+1) + 𝐦^(L+1)(𝐦^(L+1))^T))/2(_noise^(L+1))^2
        +∑_l=1^L(-1/2(_noise^(l))^2((∑_i=H_𝐱+1^N_i^(l))+ (^(l))^T_H_𝐱+1:^(l)_H_𝐱+1:).
       . + ((^(l)_H_𝐱+1:)^TΨ_1^(l)𝐦^(l))/(_noise^(l))^2 - tr(Ψ_2^(l)(𝐬^(l) + 𝐦^(l)(𝐦^(l))^T))/2(_noise^(l))^2.
       .-N/2+∑_i=1+H_𝐱-H_h^Nln(2π_i^(l))/2-∑_i=1+H_𝐱-H_h^H_hln(2π)/2 - (_i^(l) + (_i^(l))^2)/2) 
        - ∑_l=1^L+1𝐊𝐋(q_a^(l)||p_a^(l)) - 𝐊𝐋(q_Z^(l)||p_Z^(l))
       def=L^.8REVARB_.8VSS,

where we derive L^.8REVARB_.8SS by being not variational over Z^(l). Using the optimal distribution for a^(l)∼𝒩((A^(l))^-1(Ψ_1^(l))^T_H_𝐱+1:^(l),(_noise^(l))^2 (A^(l))^-1), with A^(l) = Ψ_2^(l) + (_noise^(l))^2 I_M, respective 𝐲_H_𝐱+1: for L+1 we obtain

    ≥ -N̂-M/2∑_l=1^L+1ln((_noise^(l))^2)-𝐲_H_𝐱+1:^T𝐲_H_𝐱+1:/2(_noise^(L+1))^2
        + 𝐲_H_𝐱+1:^TΨ_1^(L+1) (A^(L+1))^-1(Ψ_1^(L+1))^T𝐲_H_𝐱+1:/2(_noise^(L+1))^2 +ln(|(A^(L+1))^-1|)/2
       +∑_l=1^L(-1/2(_noise^(l))^2((∑_i=H_𝐱+1^N_i^(l)) + (^(l)_H_𝐱+1:)^T_H_𝐱+1:^(l)).
       .+ (_H_𝐱+1:^(l))^TΨ_1^(l) (A^(l))^-1(Ψ_1^(l))^T_H_𝐱+1:^(l)/2(_noise^(l))^2 +ln(|(A^(l))^-1|)/2.
       .-N/2+∑_i=1+H_𝐱-H_h^Nln(2π_i^(l))/2-∑_i=1+H_𝐱-H_h^H_hln(2π)/2 - (_i^(l) + (_i^(l))^2)/2) 
        -(L+1)N̂/2ln(2π)- ∑_l=1^L+1𝐊𝐋(q_Z^(l)||p_Z^(l))
       def=L^.8REVARB_.8VSS-opt,

where we derive L^.8REVARB_.8SS-opt by being not variational over Z^(l).
The IP-1/2-opt regularization case with A^(l) = Ψ_2^(l) + (_noise^(l))^2 I_M for 𝐦_1, and  A^(l) = Ψ_2^(l) + (_noise^(l))^2 K_MM^(l) for 𝐦_2, and Q_.8REVARB, P_.8REVARB defined with the priors and variational distributions for 𝐙^(l) and 𝐇^(l), is given by:

    ln(p(𝐲_H_𝐱+1:|𝐗))
       =…see <cit.> until Equation 14,
       JI≥ln(∏_l=1^L+1𝐄_p_a^(l)|U^(l)[exp(⟨ln(𝒩(𝐡^(l)_H_𝐱+1:|Φ^(l)(K_MM^(l))^-1𝐚^(l),(_noise^(l))^2 I_N))⟩_Q_.8REVARB)])
       -(2(_noise^(l))^2)^-1∑_l=1^L+1Ψ_0^(l)-tr((K_MM^(l))^-1Ψ_reg^(l)) - 𝐊𝐋(Q_.8REVARB||P_.8REVARB)
       = -N̂-M/2∑_l=1^L+1ln((_noise^(l))^2)-𝐲_H_𝐱+1:^T𝐲_H_𝐱+1:/2(_noise^(L+1))^2+ln(|K_MM^(L+1)|)/2
        + 𝐲_H_𝐱+1:^TΨ_1^(L+1) (A^(L+1))^-1(Ψ_1^(L+1))^T𝐲_H_𝐱+1:/2(_noise^(L+1))^2 +ln(|(A^(L+1))^-1|)/2
        -Ψ_0^(L+1)+tr((K_MM^(L+1))^-1Ψ_reg^(L+1))/2(_noise^(L+1))^2
       -∑_l=1^L(-1/2(_noise^(l))^2((∑_i=H_𝐱+1^N_i^(l)) + (^(l)_H_𝐱+1:)^T_H_𝐱+1:^(l))+ln(|K_MM^(l)|)/2.
       .+ (_H_𝐱+1:^(l))^TΨ_1^(l) (A^(l))^-1(Ψ_1^(l))^T_H_𝐱+1:^(l)/2(_noise^(l))^2 +ln(|(A^(l))^-1|)/2.
       .-Ψ_0^(l)+tr((K_MM^(l))^-1Ψ_reg^(l))/2(_noise^(l))^2-N/2+∑_i=1+H_𝐱-H_h^Nln(2π_i^(l))/2.
       .-∑_i=1+H_𝐱-H_h^H_hln(2π)/2 - (_i^(l) + (_i^(l))^2)/2)-(L+1)N̂/2ln(2π)- ∑_l=1^L+1𝐊𝐋(q_Z^(l)||p_Z^(l))
       def=L^.8REVARB_.8VSS-IP-1/2-opt,

where again we derive L^.8REVARB_.8SS-IP-1/2-opt by being not variational over Z^(l).


  §.§.§ Predictions


Predictions for each layer l and new 𝐡̂_*^(l) with Ψ_1*^(l) = 𝐄_q_Z^(l)q_𝐡̂_*^(l)[ϕ(𝐡̂_*^(l))]∈ℝ^N× M, Ψ_2* = 𝐄_q_Z^(l)q_𝐡̂_*^(l)[ϕ(𝐡̂_*^(l))ϕ(𝐡̂_*^(l))^T]∈ℝ^M× M are performed in the simple DRGP-(V)SS case with

    𝐄_q_f^(l)_*[f^(l)_*]     = Ψ_1∗^(l)𝐦^(l),
    𝐕_q_f^(l)_*[f^(l)_*]     = (𝐦^(l))^T(Ψ_2∗^(l)-(Ψ_1∗^(l))^TΨ_1∗^(l))𝐦^(l)+ (𝐬^(l)Ψ_2∗^(l)),

where q_f^(l)_*=∫ p_f^(l)_*|𝐚^(l),𝐙^(l),𝐔^(l),𝐡̂_*^(l)q(𝐚^(l))q(𝐙^(l))q(𝐡̂^(l))d𝐚^(l)d𝐙^(l)d𝐡̂^(l) and for the optimal distribution case for a^(l) we have with A^(l) = Ψ_2^(l) + (_noise^(l))^2 I_M

    𝐦^(l)opt= (A^(l))^-1(Ψ_1^(l))^T^(l)_H_𝐱+1:,  𝐬^(l)opt=( _noise^(l))^2(A^(l))^-1,

for 1,…,L, and fully analog for l=L+1 by replacing ^(l)_H_𝐱+1: with 𝐲_H_𝐱+1:.

In the DRGP-(V)SS-IP-1/2-opt case we make predictions for each layer l and new 𝐡̂_*^(l) with Ψ_0* = tr(𝐄_q_𝐡̂_*^(l)[K_**])= N _power^2∈ℝ, Ψ_reg* = 𝐄_q_𝐡̂_*^(l)[K_M*K_*M]∈ℝ^M× M through

    𝐄_q_f^(l)_*[f^(l)_*]     = Ψ_1∗^(l)Λ^(l),
    𝐕_q_f^(l)_*[f^(l)_*]   = (Λ^(l))^T(Ψ_2∗^(l)-(Ψ_1∗^(l))^TΨ_1∗^(l))Λ^(l)
       +Ψ_0* - ((K_MM^(l))^-1Ψ_reg∗^(l)-(_noise^(l))^2(A^(l))^-1Ψ_2∗^(l)),

where A^(l) = Ψ_2^(l) + (_noise^(l))^2 I_M for 𝐦_1 and  A^(l) = Ψ_2^(l) + (_noise^(l))^2 K_MM^(l) for 𝐦_2 and

    Λ^(l)opt= (A^(l))^-1(Ψ_1^(l))^T^(l)_H_𝐱+1:,

for 1,…,L, and fully analog for l=L+1 by replacing ^(l)_H_𝐱+1: with 𝐲_H_𝐱+1:.


Prediction for a new input in all cases has time complexity 𝒪((L+1)M^3), which comes from the iterative prediction through all GP layers and the calculation of the statistics, see Appendix <ref>.


  §.§.§ Distributed variational inference for DRGP-(V)SS(-IP-1/2)



We refer to <cit.>, Equation (4.3), for a comparison.

Calculating the optimal REVARB-(V)SS and REVARB-(V)SS(-IP-1/2) requires 𝒪(NM^2Q_max(L+1)), where Q_max = max_l=1…,L+1Q^(l), Q^(l)def=(𝐡̂_̂î^(l)) and 𝐡̂_i^(l) is coming from the equations in (<ref>) 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-1/2) setting with distributed inference to 𝒪(M^3) if the number of cores scales suitably with the number of training-data. We show this for L^.8REVARB_.8(V)SS-IP-opt-2, because this is the most complex bound, and all other bounds reduce to special cases of this.

L^.8REVARB_.8(V)SS-IP-opt-2 in Appendix <ref>, separated for each hidden layer and the output layer (_H_𝐱+1: and _H_𝐱+1: replaced by 𝐲_H_𝐱+1:), can be written as L^.8REVARB_.8(V)SS-IP-opt-2= ∑_l=1^L+1B_l,
where we have the KL terms

    𝐊𝐋(q_Z^(l)||p_Z^(l)) =   1/2∑_m=1^Mtr(_m^(l)+_m^(l)(_m^(l))^T) - 1/2∑_m=1^Mln(|_m^(l)|) - MQ^(l)/2,

 for l=1,…, L+1, the terms

    -N/2+∑_i=1+H_𝐱-H_h^Nln(2π_i^(l))/2-∑_i=1+H_𝐱-H_h^H_hln(2π)/2 - (_i^(l) + (_i^(l))^2)/2,

for l=1,…, L and

    B_l     =-N̂-M/2ln((_noise^(l))^2)-(∑_i=H_𝐱+1^N_i^(l))+ ^T_H_𝐱+1:_H_𝐱+1:/2(_noise^(l))^2+ln(|K_MM^(L+1)|)/2
        + _H_𝐱+1:^TΨ_1^(l)(A^(l))^-1(Ψ_1^(l))^T_H_𝐱+1:/2(_noise^(l))^2+ln(|(A^(l))^-1|)/2-Ψ_0^(l)+tr((K_MM^(l))^-1Ψ_reg^(l))/2(_noise^(l))^2
       -N̂/2ln(2π),

for all l=1,…, L+1 and which can be separated further into B_l= ∑_i=H_𝐱+1^NB_li, a sum of N̂=N-H_𝐱 independent terms, extracting ln(|(A^(l))^-1|)/2, ln(|K_MM^(l)|)/2 and tr((K_MM^(l))^-1Ψ_reg^(l))/2(_noise^(l))^2, by

    B_li   = -1-M/N̂/2ln((_noise^(l))^2)-(_power^(l))^2+_i^(l) +_i^(l)_i^(l)/2(_noise^(l))^2
        + _i^(l)(Ψ_1^(l))_ i (A^(l))^-1(Ψ_1^(l))_ i^T_i^(l)/2(_noise^(l))^2 +ln(2π_i^(l))/2 -ln(2π)/2 ,

where (Ψ_1^(l))_ i means, taking the i-th column of Ψ_1^(l) for l=1,…,L+1, i=H_𝐱+1…,N. We further inspect 

    A^(l) = Ψ_2^(l) + (_noise^(l))^2 K_MM^(l) = ∑_i=H_𝐱+1^N(Ψ_2^i)^(l)+ (_noise^(l))^2 K_MM^(l),

and

    Ψ_reg^(l) = ∑_i=H_𝐱+1^N(Ψ_reg^i)^(l).

These terms and the sums of Ψ_2^(l) and Ψ_reg^(l) can be computed on different cores in parallel without communication. Only the 3(L+1) inversions and determinants of A^(l) and K_MM^(l) 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 𝒪(M^3) per single iteration with N̂(L+1) cores.



