Learning Nonparametric Volterra Kernels with Gaussian Processes
Authors: Magnus Ross
Michael Thomas Smith
Mauricio A Álvarez
Venue: NeurIPS 2021 Poster (2021)
Keywords: Gaussian Processes

Abstract: This paper introduces a method for the nonparametric Bayesian learning of nonlinear operators, through the use of the Volterra series with kernels represented using Gaussian processes (GPs), which we term the nonparametric Volterra kernels model (NVKM). When the input function to the operator is unobserved and has a GP prior, the NVKM constitutes a powerful method for both single and multiple output regression, and can be viewed as a nonlinear and nonparametric latent force model. When the input function is observed, the NVKM can be used to perform Bayesian system identification. We use recent advances in efficient sampling of explicit functions from GPs to map process realisations through the Volterra series without resorting to numerical integration, allowing scalability through doubly stochastic variational inference, and avoiding the need for Gaussian approximations of the output processes. We demonstrate the performance of the model for both multiple output regression and system identification using standard benchmarks.

Learning Nonparametric Volterra Kernels with Gaussian Processes
Magnus Ross
Michael T. Smith
Department of Computer Science University of Sheffield, UK mross1@sheff ield. ac. uk
Department of Computer Science
University Sheffield, UK
m.t. smi tw@sheffiellleeeeeeeeeeeeeeeeefe ac. uk
Mauricio A. Álvarez
Department of Computer Science University of Sheffield, UK
mauricio .wraree0eeheffiildded ac.
Abstract
This paper introduces method for the nonparametric Bayesian learning of non- linear operators, through the use of the Volterra series with kernels represented using Gaussian processes (GPs), which we term the nonparametric Volterra kernels model (NVKM). When the input function to the operator is unobserved and has a GP prior, the NVKM constitutes a powerful method for both single and multiple output regression, and can be viewed as a nonlinear and nonparametric latent force model. When the input function is observed, the NVKM can be used perform Bayesian system identification. We use recent advances in efficient sampling of explicit functions from GPs to map process realisations through the Volterra series without resorting to numerical integration, allowing scalability through doubly stochastic variational inference, and avoiding the need for Gaussian approximations of the output processes. We demonstrate the performance of the model for both multiple output regression and system identification using standard benchmarks.
Introduction
Gaussian processes (GPs) constitute general method for placing prior distributions over functions, with the properties of samples from the distribution being controlled primarily by the form of the covariance function [22]. Process convolutions (PCs) are one powerful method for building such
covariance functions 4 S In the PC framework, the function we wish to model is assumed be generated by the application of some convolution operator to base GP with some simple covariance, and since linear operators applied to GPs result in GPs, the result another GP with a covariance we deem desirable. PCs allow models for multiple correlated output functions to be built with ease, by assuming each output is generated by different operator applied to the same base function, or set of functions 30 [3].
The PC framework unifies number different ideas in the GP literature. Latent force models (LFMs) [1] use PCs to include physics based inductive biases in multiple output GP (MOGP) models by using the Green's 's function of linear differential operator as the kernel of the convolution. This leads to the interpretation of each output as having been generated by inputting a random latent force into a linear system, with physical properties described by the differential operator. The Gaussian process convolutin model (GPCM) of Tobar et al. (29) treats the convolution kernel itself as an unknown function be inferred from data, and places a GP prior over it. Linear systems are entirely
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

described by their Green's 's function, so one can interpret the GPCM as a nonparametric, linear LFM, in which the form of the system itself inferred from data.
In the physical world, nonlinear systems are the norm, and linearity is an approximation. As a consequence it is desirable when dealing with physical data to have models that can incorporate nonlinearity naturally, Often however, given a certain set of data, it is not clear exactly what form this nonlinearity takes, and so introducing specific parametric nonlinear operators can be overly restrictive. Álvarez et al. [2] present model known as the nonlinear convolved MOGP (NCMOGP) which introduces nonlinearity to MOGPs via the Volterra series [8] a nonlinear series expansion used widely for systems identification, whose properties are controlled by a set of square integrable functions of increasing dimensionality known as the Volterra kernels (VKs). The NCMOGP assumes these functions are both separable and homogeneous, and of parametric Gaussian form, it also approximates the outputs as GPs in order to make inference tractable.
The present work introduces new model which drops the separabilityy and homogeneity assumptions on the VKs, allows their form to be learned directly from data, and makes no approximation on the distribution of the outputs. We refer to it the nonparametric Volterra kernels model (NVKM). We develop a fast sampling method for the NVKM which leverages the recent results of Wilson al. 335 the sampling of explicit functions from GPs to analytically map function realisations through the Volterra series, avoiding the need for computationally expensive and inaccurate high dimensional numerical integration. Fast sampling allows for the application of doubly stochastic variational inference (DSVI) (28) for scalable learning.
The NVKM is well suited to both single and multiple output regression problems, and can be thought of as an extension of the GPCM to both nonlinear systems and multiple outputs. The NVKM can also be interpreted as nonlinear LFM in which the operator is learned directly from data. We additionally present a variation of the NVKM that can be used for Bayesian systems identification, where the task is to learn operator mappings between observed input and output data, and show that it allows for considerably better quantification of uncertainty than competing methods which use recurrence [8]
2 BBackrroudd
In this section we give a brief introduction to the mathematical background of PCs and the Volterra series.
Process convolutions In the PC framework, the set of output functions (fa(t)) with D being the number of outputs, is generated by the application of some set of linear operators, specifically convolution operators, to a latent function u represented by a Gaussian process, fa(() S Ga(t T)u(T)dT, where T is the domain of integration, and the function G3 is known variously as the convolutional kernel, smoothing kernel, impulse response or Green's 's function, depending on context. Assuming that the input is bounded, the function Ga must be square integrable to ensure the output finite. linear operator acting on GP produces another GP 22 and so we obtain D distinct GPs. Since the latent function u is shared across the outputs, these D GPs are correlated, allowing joint variations to be captured, whilst the convolution with Ga adapts to each output. Álvarez et al. [3] show that many MOGP models can be recast in terms of the PC framework by particular choices of Gd and u. In LFMs, Ga is taken to be the Green's 's function of some differential operator. We can then interpret each output as resulting from a shared random force being fed into a distinct linear system, represented by the differential operator. The smoothing kernels are usually taken to have parametric, often Gaussian form. Tobar et (29 use the PC framework for a single output, and make the smoothing kernel itself GP. giving rise to the GPCM. Bruinsma 7 extends the GPCM to the multiple output case, although the model was not applied to data.
Volterra Series The PC framework can be extended to represent broader class of output functions by instead considering the outputs {fa as being the result of some nonlinear system, acting on the latent function u. The Volterra series is a series approximation for nonlinear systems that is widely used in the field of system identification [8. It is given by,
/ f(1) - Σ T
Golt T...., To)
ul u(T dr,
(I)
2

---

where f is the system output. Ge is the ct order VK, is the system input, and C' is the order of the approximation. We can think of the Volterra series as an extension of the well known Taylor expansion, that allows f to have memory of the past values of u, that is to say I' depends on at all values & T. The Volterra series can approximate broad class of nonlinear operators, however unable to represent certain properties of general nonlimear systems, for example chaos. Álvarez et al. 12 use Equation (I to construct the NCMOGP, applying the Volterra series with d distinct sets shared latent GP input u, to produce outputs {ff(tt) Since the Volterra series is nonlinear operator, the output becomes an intractable, non-Gaussian process. The authors perform inference by approximating the outputs as GPs and using the first and second moments of the output process to form its mean and covariance function. To enable to computation of these moments, the authors restrict the set of VKs to those which are both separable and homogeneous, i.c. Ga.o(ti...., te) I Gt(t.). Additionally, since the moment computation requires analytically solving number of non-trivial convolution integrals, the authors only consider a Gaussian form for the VKs.
of VKs {Ga,e)om 1e
3 The nonparametric Volterra kernels model
The NVKM relaxes the restrictions of separability and homogeneity which are placed on the VKs in the NCMOGP, and represents these kernels as independent GPs, allowing their form and uncertainty to be inferred directly from data. The generative process for the NVKM can be stated as
v(t) ~ GP(O, k (1,1)),
(t. ~ GP(O, (Ga.)) (G,aett, (t,t)), t tER, Ve E I,...,C, Ad 1, I,...,,,
(2)
C fa(t) Σ Ga,olt - -I.... Te) II u(T)))dT,
where (:(()(t t') is the covariance function for the imput process, and (GGa, t') is the covariance function for the th VK of the thh output. We follow Tobar et al. [29] in using the decaying square ponential (DSE) covariance for the VKs, which is a modification to the ubiquitous square exponential (SE) covariance that ensures the samples have finite energy [ The DSE covariance has the form
( t') = 0² exp(-a(|lt|)--- ||t'[]2) - - +²))²)
where || || is the or e norm, a is the amplitude, controls the rate at which the samples decay away from the origin, and y is related to the length scale / of the samples by 7 #. A diagram of the generative process for the model is shown in Figure Obtaining an exact distribution over the outputs fu intractable, since involves integration over nonlinear combinations of infinite dimensional stochastic processes. In order to sample from the model, and perform inference, approximations must be introduced. In particular, we employ the results of Wilson et al. [35] to sample linear
time, which enables efficient learning through the use of variational inducing points (27 with doubly stochastic variational inference (DSVI) 28
3.1 Sampling
One could sample from the model by drawing from the input and filter GPs at some finite set of locations, and then using the samples with some method for numerical integration find the output. As the dimensionality of the filters increases, however, many points would be needed to obtain an accurate answer, this quickly becomes computationally intractable, since sampling exactly from GP has cubic time complexity with respect to the number of points requested. We can sidestep this problem, and avoid the need for any numerical integration. by representing samples from the GPs explicitly as functions. Using the results from Wilson et al. (35), and following their notation, we can write sample from a GP f R - IR, with covariance function ke(t, t'), given M inducing variables v E IR M with corresponding inducing inputs {z,
IM 13=1, with z, E R", as
(flv)(t)= Σ widi(() + Σ gtk(t,z))
'This shown in Appendix A of the supplemental material.
3

---

Gaite
Gg)lfi.
1u(9)
ymn d-1
*
-ymm
/ d=2
Now
/ (outt-rju(rddf
Gealt - a.t (..t--))a(nn(eeee ()))())((s)))ddr
fa(t)
Figure 1:
A diagram of the generative process for the NVKM with C 2 and D showing the stages of computation Volterra series for three outputs, shown in the rows, with the star representing convolution. The 1D contribution from the second order term is obtained by taking the diagonal of the result of the 2D convolution.
the first order (on the left side) and second order (in the center) terms the
where f is random Fourier basis, with No being the number of basis functions in the approximation, w E R&N, with entries wy ~ (0, 1), and are the entries of the vector = (v kw), where K E RMX with elements Km,n = k(zn, Zm). is the covariance matrix of the inducing points, and & E Mr is a feature matrix, with each basis function being evaluated at each inducing location. The random Fourier basis is obtained by first sampling |3, (O, 2 where U is the uniform distribution, and then sampling 0; ~ (K))., where FT is the Fourier transform of the covariance function, which is the spectral density of the process. The basis functions are then given by (t) 22/N, cos(O, (0, t B:). We can see Equation (4) as consisting of an approximate GP prior, using a random Fourier features approximation [21], with a correction term which uses Matheron's update rule account for the inducing points. By using Equation 4 samples can be obtained in linear time with respect to the number of requested points. It should be noted that Equation (4) only applies GPs with stationary kernels. however, the DSE covariance required for Ged, is non-stationary. It can be shown that the process
(l'' (tt)) has the DSE covariance if G²,dd has the SE covariance. We can then write sample from the output of the NVKM
as,
(fallv
Vilveelleivv (t)
1- E/e - (Gael -........ Te) Σ
(5)
II u (T;))
which can computed analytically, assuming that u also has the SE, covariance, by representing the Fourier basis in complex form, and factorising the integrals, leading to combinations of sums and products of single dimensional integrals of the form f exp(-ax² + bz)dx VT/aexp(b") (app)b)/)))) See Appendix B in the supplemental material for details of the computation.
3.2 Inference
Learning with the NYKM implies making inference of the input process u, along with all VKs C,D
(G.,e) from observed output data {ya) with Ya E Ma, which are the functions ( evaluated at points (ta)2_ with t, E Na corrupted by some i.i.d Gaussian noise. That is to say ya. fa(ta,) + ca.: with e ~ N(O,, o. Let v Gz denote the inducing points for the VKs, and " u(z") denote the inducing points for the input. The joint distribution over these

---

inducing points and the latent functions then has the following form
plyvalan (Gd.o: v "e U, u,,,,, v²)=
D.Na
D.C II (GGaovddp(yy) d,c=1
II p(yas)Sa((a.))
(6)
where fa(ta,)) depends on the VKs and input through Equation ( the likelihood is ((yaslla((a..)) N Nya, fa(tas), ) p(Ga,alvā,d) and p(ulv") are GP posterior distributions, and p(v") and p(vV.c) are the prior distributions over the inducing points. The dependency structure of the model is described in Figure 2 We form an approximate variational distribution, in similar way to Tobar et al. (29), using structured mean field approximation. That is to say, we mirror the form of the true joint distribution, and replace the prior distributions over the inducing points with variational distributions, q and q(v"), leading to the variational distribution
dee)
D.C
MiGd, V.G. G.D d,c u, v")= = II d,c=1
p(Ga,el (
Given the assumed factorisation of the variational distibution, the optimal form of the variational posteriors are multivariate Gaussians, g(v")= =Nv";; ,4", Σ") and g(v =N(y Si fd,e Eg,) where the mean vectors and covariance matrices of these distributions are variational parameters. This form of variational approximation leads to a variational lower bound,
D.Na
D.C
F E.loop(uuslfatta.)))
- Σ KL(g(v (llp(v)))) KL[g(v") Ip(v(") (8)
where KL[-II-] represents the Kullback -iebler (KL) divergence, for details Appendix C in the supplemental material. The expression above is optimised using gradient descent. The KL divergences have closed form. The derivation of the bound and KL divergences material. The expectation of the log likelihood of the outputs, given in the first term, is intractable, due to the nature of nonlinearity introduced by the Volterra series. We instead compute stochastic estimate of the log likelihood by sampling from the model and using,
given in the supplementary
B.y pogg
log oogalaallaaall (lai)))),,
S
where vG and " are first sampled from their re- spective variational distributions, and then used in Equation (5) to generate a sample from fa. To make the inference scheme scalable, we compute the bound on randomly sub-sampled mini batches of the data set, which alone is known as stochastic variational inference [14. When this source of stochasticity combined with the stochastic estimate of the expected log likelihood, we have DSVI 28
Gá.e
ya,e
ta
L....
fa
V
In the standard NVKM model, the input process u a latent function with no observed data associated with it. There are many situations in which, instead of learning a distribution over some output functions alone, we wish to learn an operator mapping between an input function and an output function, functions. That is to say, in addition to observing data {ya)2 we also observed the input process u at locations t. E R= corrupted with i.i.d noise, which we denote
Figure 2: A graphical model for the NVKM, where the dashed elements added to form the IO-NVKM. Note that nodes u, fa, and Gd.e not random variables, they ran- dom processes, but the distinction is not made in the diagram for the sake of clarity.
x, so I;; u((T) + ee with e (0, az) We apply a simple modification to the inference scheme
of the NVKM, in order to form new model which we term the input/output NVKM (IO-NVKM).
5

---

For the IO-NVKM, we pick up an additional likelihood term in Equation ( with the variational distribution remaining unchanged, and we obtain a new bound,
N₂ Fi0 = = Σ Eg Eo
(10)
where p(ry (((()()))) (ary; 0, 22), F is given by Equation (8) and approximate the expecta-
tion as in Equation (9) The relationship between the NVKM and IO-NVKM is illustrated in Figure
where the dashed sections are added to form the IO model.
4 Related Work
In the this section, we give a brief overview of existing ideas in the literature that have some connection to the NVKM.
Nonparametric covariances Covariance function design and selection is key to achieving good performance with GPs. Much work has been done on automatically determining covariance functions from data, for example by building complex covariances by composing simpler ones together [10], or by using deep neural networks to warp the input space in complex ways before the application of simpler kernel [3 Flexible parametric covariances, can also be designed in frequency space [3 Alternatively, some efforts have been made learn covariance functions using GPs themselves. In addition the the GPCM, another model that learns covariances nonparametrically is due to Benton et al. (5) who use GPs to represent the log power spectral density, and then apply Bochner's theorem to convert this to a representation of the covariance function. GPs are fully specified by their first two moments, so by learning the covariance function and mean function, one knows there is to know about the process. The present work uses the formalism of the Volterra series to learn the properties of more complex, non- non-Gaussian processes, nonparametrically. We can think of this as implicitly learning not just the first and second order moments of the process. but also the higher moments, depending on the value of C. In the case of = 1 and = 1, the NVKM and the GPCM are the same, except for the fact the GPCM uses the white noise the input process, whereas we use SE GP.
LFMs and MOGPs As discussed in Section 2 can interpret the first order filter function as the Green's ' function of some lincar operator or system, and so by placing a GP prior over it, we implicitly place a prior over some set of linear systems [29]. Since standard LFMs use an operator of fixed form, we can interpret the NVKM in the case C = 1, 2 1 as being an LFM in which the generating differential equation itself is learned from data. LFMs can be extended cases in which the differential operator is nonlinear. Hartikainen et al. [12] recast a specific nonlinear LFM in terms of a state space formalism allowing for inference in linear time. Lawrence al. [17 use Markov chain Monte Carlo to infer the parameters of a specific nonlinear ODE describing gene regulation, using a GP as latent input function. Ward et al. 31 use black box VI with inverse auto-regressive flows to infer parameters of the same ODE. When C > 1, the NVKM can be interpreted as a nonlinear nonparametric LFM. Álvarez et al. 12 use the fixed, parametric VKs to build an MOGP model. In contrast to the NVKM, they approximate the outputs as GPs and use analytical expressions for the moments to perform exact GP inference.
Nonlinear system identification The IO-NVKM falls into the class of models which aim to perform system identification. A key concern of systems identification is determining how the output of certain systems, often represented by differential operators, respond to a given input. GPs have long been used for the identification of both linear and nonlinear systems. Many models exist which use GPs to recurrently map between previous and future states, including GP-NARX models [6). various state space models [26] 25) and recurrent GPs (RGPs) [9] The thesis of Mattos [88 gives summary of these methods. Worden al. (36] detail a method for the combination of the GP-NARX model with mechanistic model based on the physical properties of the system under study, leading to improved performance over purely data driven approaches. The IO-NVKM differs from these models in that instead of learning a mapping from the state of the system at given point to the next state, we use GPs to learn an operator that maps the whole input function to the whole output function.

---

M m
916
GPEN
NYRMIC (1
-6-
NVKM (C=2)
-8-
NVRMIC-
-10
-10
-20
-15
-10
10
15
20
-000
-9.5
-9.0
-8.5
-8.0
-55
-70
Figure 3: Model predictions the synthetic data, with crosses indicating training points and dots indicating test points. The right plot shows an enhanced view of the peak around t = -9, and the shaded regions show 2o confidence.
Experiments
For all the following experiments we place the inducing locations for both the input process and VKs on a fixed grid. The GP representing the e VK has input dimension c, which means that the number of inducing points required to fully characterise it scales exponentially with if they are placed on grid. For all experiments we use 15, 102 100, 6 216 and 4 256 inducing points, for each of the 1s to 4th order filters respectively, centered zero. We treat the range of the points of each VK as hyperparameter, and fix such that the decaying part of the DSE covariance causes samples to be near zero at the edge of the range. For we use approximately 1/10 of the number of inducing points as data points (average number data per output for multi-output problems). The VK GP length scales, VK GP amplitudes, and imput GP amplitude parameters by maximising the variational bound using gradient descent. For computational reasons, the input process length scale is fixed based the spacing of the input inducing points, and the noise hyperparameters are fixed to a small value whilst the other hyperparameters and variational parameters optimised, and then estimated afterwards by minimising the bound with all other variables fixed. The model is implemented using the Jax framework
optimised along with the variational
[6] For all experiments we Adam [15]. All models were trained a single Nvidia K80 GPU. The code is available github. com/ magnusross /osk/ nvkm In the following sections, details of the three main experiments given. An additional experiment, in which data is sampled from the model, and the aim is to recover the VKs that generated it, is given in Appendix of the supplemental material.
5.1 Synthetic data
To illustrate the advantage of including nonlinearity in the model we generate a synthetic single output regression prob- lem which includes both hard and soft nonlinearities by sampling y from an
Table 1: Performance on the synthetic data set, showing mean and standard deviation for 10 repeats.
Model
NMSE
NLPD
GPCM
0.199 +0.023
1,080 13.130
SE GP with length scale 2, computing fa(() = fe ha(T)g(t 7)dr for h1 = (6t), he(t) sin"(5t) and
NVKM (C =)
0.196 +0.047
2.084 #0.398
NVKM (C 2)
0.108 #0.065
0.638 #0.580
NVKM (C = 3)
0.055 #0.016
0.124 10.107
hy(t) cos(4t) by numerical integration.
NVKM (C = 4)
0.084 +0.087
0.149 20.331
then computing the output as,
(t) (0)2))))())())) 51)
(11)
with (0, 0.05?). We generate 1200 points in the range t. 20, 20) and random subset of a third for training and the rest for testing. Table [ shows the normalised mean square errors (NMSEs) and negative log probability densities (NLPDs) the test set for the NVKM with various values of C' as well as the GPCM, with repeats using different random train/test split, and different random seeds As we would expect, the NMSE values are very similar for the NVKM
Results generated using the implementation available at github. www.t....oo.aeeeeee. wttmm//seeeb/gpo

---

20-
10
10
0-
-500
250
250
500
-500
-250
250
500
-500
-250
250
500
(uipat(V)
7.5
20
50-
2.5
0.0
3000
4000
Timo
Figure 4: (Top) Diagonal of the inferred Volterra kernels for the the IO-NVKM with C' 3, showing 2o confidence region. (Bottom) The predicted output for the test set, with the dashed line showing the true values.
with C and the GPCM, since the models are nearly equivalent except for the prior on the input GP. Interestingly the NL.D values are better for the GPCM than the NVKM with C 1, likely due to the fact we do not optimise the noise jointly with the bound. As C' increases the performance of the NVKM improves until = 4. The fact performance does not improve after = 3 illustrates the difficulty of identifying higher order nonlinearities in relatively small training set, effect supported by the results of the Cascaded Tanks experiment the following section. Although the
C'= 4 model does have more capacity to represent nonlinearities, the optimisation procedure is challenging, illustrated by the high variance of the results. Plots of the predictions for the model can be seen in Figure B] We can that increasing the nonlinearity for the NVKMs allows the sharp spike and the finer grained features, as well as the hard nonlincarities, to be captured simultaneously.
5.2 Cascaded tanks
To demonstrate the IO-NVKM, we standard benchmark for nonlinear systems identification know Cascaded Tanks (24] The system comprises two vertically stacked tanks filled with water, with water being pumped from reservoir to the top tank, which then drains into the lower tank and finally back to the reservoir. The training data is two time series of 1024 points, one being the input to the system, which is the voltage fed into the pump, and the second being the output, which is the measured water level in the lower tank. For testing. an additional input signal, again of 1024 points, is provided, and the task is to predict the corresponding output water level. The system consid- ered challenging because it contains hard nonlinearities when the tanks reach maximum capacity and overflow (see the regions around 600s and 2900s in Figure 4 it
Table 2: Comparison of performance on the Cascaded Tanks dataset, with the last best models reported in (18). H indicates the
number of hidden layers in the RGP.
Model
RMSE
NLPD
IO-NVKM (C 1)
0.835
1.724
1O-NVKM (C 2)
0.716
1.311
IO-NVKM (C 3)
0.532
0.879
IO-NVKM (C' 4)
0.600
0.998
RGP (H = 1)
0.797
2.33
RGP (H =
0.308
7.79
GP-NARX
1.50
1080
Var. GP-NARX
0.504
119.3
has unobserved internal state, and has a relatively small training set. Table shows the predictive root mean square eтrors (RMSEs) and NLPDs for the IO-NVKM with various C, well as four other GP based models for system identification from [18]. For each C, five random settings of VK ranges were tested, and each training was repeated three times with different initialisations. The setting and initialisation with the lowest combined NLPD on the training input and output data is shown. Although the RGP with H = 2 provides the best RMSE of the model, this comes the cost of poor NLPD values. All 1O-NVKMs achieve considcrably better NLPD values than the alternatives indicating much better quantification of uncertainty. Of the IO-NVKMs, C 3 performs best in both metrics. Figure 4 show the predictions of the C' = model on the test set, as well as the inferred VKs. The uncertainty in the VKs increases with their order, which is natural given the difficulty of estimating higher order nonlinear effects from a small training set. We can see the first order
Available sites www.ggooo..oooooooo http://..c.ee/eeee///ooooeeeeeeo benchhaark

---

24
24
22
Test
22-
20
Chimet Temperatre
20-
18.
18.
16-
16-
14-
14-
10 11
12
13
14
15
10
11
12
13
14
15
(days)
Figure 5: Predictive means and 2o confidence regions from the NVKM with C' 3, the Cambermet (left) and Chimet (right) outputs from the Weather data set. Orange circles are the artificially removed test data, black pluses training data.
Table 3:
Comparison of performance on the Weather data set, for the NVKM mean and standard
deviation of three initialisation is shown, along with the best model from 2
Cambermet
Chimet
Model
NMSE
NLPD
NMSE
NLPD
NVKM (C' 1)
0.212*0.088
2.182+0.743
1.669#0052
7.148+0.111 1.14
NVKM (C = 2)
0.440+0.286
3.88442.380
0.939=0.216
NVKM (C 33
0.253+0.002
2.390=#.123
0.871#0.394
3.994+1.924
NCMOGP (C (C=3) 3)
0.44
2.33
0.43 2.18
kernel approximating the response of linearized system, and the third order allowing the sharp nonlinearities around 600s and 2900s to be captured. It should be noted that Worden et al. [36] achieve a much lower RMSE of 0.191 by using specific physics model of the system tandem with GP-NARX model, but since we are considering purely data driven approaches here, it is not directly comparable.
5.3 Weather data
To illustrate the utility of the NVKM for multiple output regression problems, we consider popular benchmark in the MOGP literature, consisting of multiple correlated time series air temperature measurements taken at four nearby locations the south coast of England, originally described by Nguyen et al. [20), which we refer to as Weather The four series are named Bramblemet, Sotonmet, Cambermet and Chimet, with 1425. 1097, 144I, and 1436 data points, respectively. Bramblemet and Sotonmet both contain regions truly missing data, 173 and 201 points in a continuous region artificially removed form Cambermet and Chimet with the task being to predict them based on the all the other data. Table 3 shows the performance of the multiple output NVKM on the Weather dataset, along with the best performing NCMOGP model of Álvarez et al. (2]. For each C, five random settings of VK ranges different initialisations, the setting with the best average NLPD value on the training data is shown, along with its standard deviation. All NVKM models show better equivalent performance than the NCMOGP on the Cambermet output, but all show worse performance on the Chimet output, although the Chimet output the variance between repeats is high. It should be noted that the LFM reported by Guarnizo and Álvarez [1]] achieves much lower scores, having NMSEs of 0.11 and 0.19 on Cambermet and Chimet respectively, but that model uses six latent functions as opposed to a single latent function for the NVKM and NCMOGP, Including multiple latent functions may lead to large performance improvements for the NVKM and is a promising direction for future work.
tested, with each training being repeated three times with
Available for download in convenient from using the wbml package, thub. com/veeseeb/vvbml

---

6 Discussion
Societal Impacts is possible that the NVKM model could incur some negative societal impacts. GP and MOGP models have long been applied to problems in robotics [9] 32]. Better inclusion of nonlinearities in these models may enhance the ability of robots, potentially leading loss of jobs and livelihoods to automation. On the other hand, the NVKM may have positive impact on society, its ability to model weather and climate data have been demonstrated in Section [ and accurate climate forecasting models will likely play a key role in any solution to the ongoing climate crisis
23
Future Work There are a number of extensions to both the NVKM and IO-NVKM that could lead to substantial improvements in performance. As briefly mentioned in Section S the number of inducing points required for the VKs scales exponentially with the order of the series, meaning it is difficult to represent complex features in the higher order terms, without using computationally intractable number of points. Whilst initially we saw the increased flexibility of non-separable VKs as a virtue, may be that introducing separability leads to more powerful models, since the number of points needed to specify separable VKs scales linearly. Currently the models do not support multidimensional inputs, but this could be easily added, requiring the computation of few extra integrals, with the complexity scaling linearly with the number of inputs. For the multiple output model, allowing a shared set of latent functions, with the input to each output's Volterra series being trainable linear combination, in a similar way to LFMs, is highly likely to improve performance especially for problems with a large number of outputs. Additionally it is likely that the inference method can be improved significantly, perhaps by going beyond the structured mean field approximation and using a more flexible variational distribution, or even by forgoing VI entirely and sampling from the posterior using Markov chain Monte Carlo methods.
Conclusions We have presented a novel model which uses Gaussian processes to learn the kernels of the Volterra series nonparametrically. allowing for the effective modeling of data with nonlinear properties. We have developed fast and scalable sampling and inference methods for the the model and shown its performance modification to the model was presented that achieves significantly better uncertainty quantification than competitors challenging benchmark for nonlinear systems identification.
single and multiple output regression problems. Additionally,
Acknowledgments and Disclosure of Funding
We thank Lizzy Cross, Felipe Tobar, and Carl Henrik Ek for helpful conversations. We would also like to thank Wessel Bruinsma for useful and insightful discussions and advice, well support in using the code for the GPCM. Magnus Ross and Michael T. Smith thank the Department of Computer Science at the University of Sheffield financial support. Mauricio A. Álvarez has been financed by the EPSRC Research Projects EP/R033303/1, EP/T000433/2 and EP/V029045/1.
References
[1] M. Alvarez, D. Luengo, and N. D. Lawrence. Latent force models. In Artificial Intelligence and Statistics, pages 9- 16. PMLR, 2009.
[2] M. Álvarez, W. Ward, and C. Guarnizo. Non-linear process convolutions for multi-output Gaussian processes. In The 22nd International Conference Artificial Intelligence and Statistics, pages 1 1969-1977, 2019.
[3] M. A. Álvarez, L. Rosasco, and N. D. Lawrence. Kernels for vector-valued functions: A review, Found. Trends Mach. Learn,, 4(3): 195-266, Mar. 2012. ISSN 1935-8237. doi: 10.1561/2200000336 URL nutps://doi... orggooo//// 1 2200002000
[4] R. P. Barry and J. M. Ver Hoef. Blackbox kriging: spatial prediction without specifying variogram models. Journal of Agricultural, Biological and Environmental Statistics, 133):297322, 1996.
[5] G. Benton, W. J. Maddox, J. Salkey, J. Albinati, and A. G. Wilson. Function-space distributions over kernels. In Advances in Neural Information Processing Systems, pages 14965-14976. 2019.
[6] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke,
J. VanderPlas. S. Wanderman Milne, and Q Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL httphtt::/// htatt:://gihhut com/ google/jax
10

---

[7] W. Bruinsma. The Generalised Gaussian Process Convolution Model. MPhil. thesis, University of Cambridge, 2016.
[8] C. Cheng, Z. Peng. W. Zhang, and G. Meng. Volterra-series-baaed nonlinear system modeling and its engincering applications: A of art review. Mechanical Systems and Signal Processing, 87: 340 364, 2017. ISSN 0888-3270. doi: http:///...opp..0100001ooommmppp
URL httpp
com/ sci http////aattcccc//////////222000
[9] M. P. Deisenroth, D. Fox, and C. E. Rasmussen. Gaussian processes for data-efficient learning in rohotics and control. IFFE transactions on pattern ondlysis and machine intelligence, 37(2):488 2013.
[10] D. K. Duvenaud, H. Nickisch, and C. Rasmussen. Additive Gaussian processes. In J. Shawe-Taylorr R. Zemel, P. Bartlett, Pereira, and K. Q. Weinberger, editors. Advances in Neural Information Processing Systems, volume 24, Curran Associates, Inc., 2011, http::://tppeeeeeeeeeeeeee
neurips. cctppper/
2011/fiee/<c5bde745 www.eeeddde4881111615611 1055617490200 37800 001-667774900007777888----P- Paper. pdf
[11] C. Guarnizo and M. A. Álvarez. Fast kernel approximations for latent force models and convolved multiple- output Gaussian processes. In Uncertainty in Artificial Intelligence: Proceedings of the Thirty-Fourth Conference (2018). AUAI Press, 2018.
[12] J. Hartikainen, M. Seppänen, and S. Särkkä. State-space inference for non lineir latent force models with application to satellite orbit prediction. In Proceedings of the 29th International Coference on International Conference on Machine Learning, pages 723-730, 2012.
[13] D. Higdon. Space and space-time modeling using pracess convolution. In Quantitative methods for current environmental issues, pages 37-56. Springer, 2002.
[14] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14(5), 2013.
[15] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In Y. Bengio and Y. LeCun, editors, 3rd International Conference on Learning Representationss ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015. URL http:://axii wrgg/bs//41/2/
[16]] Kocijan, A. Girard, B. Banko, and R. Murray- Smith. Dynamic systems identification with Gaussian processes. Mathematical and Computer Modelling of Dynamical Systems, 1 11(4):11 424, 2005. doi: 10. 10.1080/1387300 38733995000006666. URL. http:htttt//aaat org/ 10. 1et.100001110 13771/889399500085577
[17] N. Lawrence, M. Rattray, and M. Titsias. Efficient sampling for Gaussian process inference using control variables. In D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, editors, Advances in Neural Information Processing Systems, volume 21. Curran Associates, Inc., 2009. htRpssstttt// / htt.://220...0000000000000666553
/ttp//eeecceeee..cc/oe neurips pdf
[18] C. C. Mattos. Recurrent Gaussian processes and robust dynamical modeling. PhD thesis, 2017.
[19] C. C. Mattos, Z. Dai, A. Damianou, J. Forth, G. A. Barreto, and N. D. Lawrence. Recurrent Gaussian processes. In H. Larochelle, B. Kingsbury. and S. Bengio, editors, Proceedings of the International Conference Learning Representations. volume 3, Caribe Hotel, San Juan, PR, 2016. URL httpp //anverseprobability hwtpbiiiiittt.yoooouuaaeeaeeaeee
rew.oosseeruueet.16. html
[20] T. V. Nguyen, F. V. Bonilla, et al. Collaborative multi-output Gaussian processes. In UAI, pages 643-652. Citeseer, 2014.
[21] A. Rahimi and B. Recht. Random features for large-scale kernel machines. Advances in Neural Information Processing Systems, 3, 2007.
[22] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press, 2005. ISBN 026218253X.
[23] D. Rolnick, P. L. Donti, L. H. Kaack, K. Kochanski, A. Lacoste, K. Sankaran, A. S. Ross, N. Milojevic- Dupont, N. Jaques, A. Kording, Gomes, A. Ng, D. Hassabis. J. C. Platt, F. Creutzig, Chayes, and Y. Bengio. Tackling climate change with machine learning, 2019.
A. Luccioni, T. Maharaj. E. D. Sherwin, S. K. Mukkavilli, K. P.
[24] M. Schoukens and P. Noël. Three benchmarks addressing open challenges nonlinear system identifica- tion. IFAC- PFFC-PaapesOOnLine, 50(1):446-44 2017.
[25] A. Svensson and T. B. Schön. A flexible state-space model for learning nonlinear dynamical systems. Automatica, 80:1 1991199 2017.
[26] A. Svensson, A. Solin, S. Särkkä, and T. Schön. Computationally efficient bayesian learning of Gaussian process state space models. A. Gretton and C. Robert, editors, Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, volume 51 of Proceedings of Machine Learning Research, pages fress/v51/svenseee html
Cadiz, Spain, 09- 11 May 2016. PMLR. URL http::/procceediiiiee mlr

---

[27] M. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics, pages 567-574, 2009.
[28] M. Titsias and M. Lázaro-Gredilla. Doubly stochastic variational Bayes for non-conjugate inference. In International conference machine learning, pages 1971-1979 PMLR, 2014.
[29] F. Tobar, D. Bui, R. Turner. Learning stationary time series using Gaussian processes with nonparametrie kernels. In Advances Neural Information Processing Systems, pages 3501-3509, 2015.
[30] J. M. Ver Hoef and R. P. Barry. Constructing and fitting models for cokriging and multivariable spatial prediction. Journal of Statistical Plannig and Inference, 69:275-294, 2944, 1998.
[31] W. Ward, T. Ryder, D. Prangle, and M. Alvarez. Black-bor box inference for non conear latent force models. In International Conference Artificial Intelligence and Statistics, pages 3088 3098. PMLR, 2020.
[32] C. Williams, S. Klanke, S. Vijayakumar, and K. Chai Multi-task Gaussian process learning of robot inverse dynamics. In D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, editors, Advances in Neural Information Processing Systems, volume 21, Curran Associates, Inc., 2009. BtRphttps:/ fttt:s://prrceediggg neurrip cc/paper/2008 /tt///feeee6666899189797777/1//0
.pff
[33] A. Wilson and R. Adams. Gaussian process kemels for pattern discovery and extrapolation. In International conference machine learning, pages 1067-1075, 2013.
[34] A. Wilson, statistics. pages 370-378, 2016.
R. Salakhutdinov,
Xing. Deep kernel learning. Artificial intelligence and
[35] J. Wilson, Borovitskiy, A. Terenin, Mostowsky, and M. Deisenroth. Efficiently sampling functions from Gaussian process posteriors. H. III and Singh. editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 10292- 10302. PMLR, 13-18 Jul 2020. /tLL/ttt://pooceediiii
mlr prrssprssssee.e wwwps////11/////eeemmam htmmll00..mmmm11mm
[36] Worden, Barthorpe, Cross, N. Dervilis, Holmes, Manson, and Rogers. evolutionary system identification with applications to nonlinear benchmarks. Mechanical Systems and Signal Pro- cessing, 112:194-232, 2018. ISSN 0888-3270 3270. doi:
10..0 URL.
https://www sci sciencedirect. com sci http:////o...oo........o///////.
12


Reviews:
[{"review_id": "vU_waviemrl", "reviewer": "Reviewer_Aww8", "summary": "The paper expands previous work to learn convolutional operators that map a latent (unobserved) GP to observed data. The convolutional operators/kernels are more general as in previous works (non-factorizing and non-homogeneous kernels), and inference in made scalable via a combination of a variation approach and recent results to efficiently sample GPs. Furthermore, the authors provide an extension, where the latent input process is also associated with data.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I found the paper interesting and very well written. It is providing some methodological expansion on previous work, such that less assumption on the Gaussian process convolutional models (GPCM) are needed. While the expansions of the model are straightforward, the variational inference of the model needed to utilize some recent findings on GP-sampling to get an efficient approximation of the ELBO. To use the sampling technique from [Wilson et al, 2020] the authors transform a GP with non-stationary kernels, to a stationary one. Then, the nontractable parts of the ELBO can be computed efficiently, and standard stochastic VI methods can be employed.\nThe experiments seem to be done thoroughly. However, I’d have some suggestion for improvements (see below). \nOverall, the paper is not the most innovative one, since novelty mainly in combinations of recent findings and an interesting, but straightforward model extension. But ideas are executed neatly, and I consider it as some solid work, which could be published at this stage.\n\n## Update\n\nI thank the authors for their response. Considering the additional information and experiments, that will be added to the final publication, I am happy to increase my score to 7.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "+ l. 153: I don’t think, one can say that the optimal variational distributions are Gaussians. Optimal would be, that they turn out to be Gaussians under the factorizing assumption of Eq. (7). Here, however, I believe it is an assumption that needs to be made.\n+ Sec.5.1: While I agree, that the synthetic example is a difficult one, it would be insightful, if an inference example can be shown, that is coming from the generative model. Can the model then correctly identify C? Or does the fact, that the authors choose a decreasing number of inducing points for larger C’s corrupt the model selection?\n+ Also for the real experiments, it would be nice to choose an example, where learnt filters can be interpreted in a meaningful way, and potentially demonstrate, why it is important to assume non-homogenous and/or non-factorizing filters, instead of just referring to the likelihood.\n+ As discussed, I also consider the fact, that the input process u needs to be univariate the strongest restriction. I would appreciate, in the discussion, how many dimensions are likely to be feasible.\n+ I am missing the report of runtimes in the experiments. How does it scale with number of data points, C, etc. How does it compare to the more limited models? Does one pay a high price in terms of computation time for the proposed model?\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5h"}, "extra_text": {"ethical_concerns": "I found the comment on drones (l. 324ff) a bit out of place. Surely, it is not impossible, that the proposed algorithm could be used in such a malicious way, but the comment seems a bit farfetched, and any technology can be used in such away. Since neither methodology, nor experiments point towards this application, I do not consider this comment as crucial here."}}, {"review_id": "hPHHee2-rjG", "reviewer": "Reviewer_aBJf", "summary": "The paper places GP priors over the kernels of a Volterra series which is driven by an input process drawn from a GP, which may optionally have associated observations (the IO variant of the model). Sampling and mean field variational inference procedures are described. The method is demonstrated on two real data sets.", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper seemed to be inspired by [2]. However, it moved away from the full GP approximation to directly model the the output function via equation (5). This make the paper more original and significant. Sampling and inference procedure are clearly described. \n\n# COMMENTS\n1. It will be worthwhile to see if (Markov chain) Monte Carlo posterior sampling of the model can yield better results than currently reported. This will allow the reader to assess if the limitation is because of the model or because of the approximate inference. Section 3.1 is a good start to build this upon. In the least, a comment on this aspect will be useful for the reader.\n\n2. In the experiments, a number of runs via different random seeds is needed, and the best results are reported (though see point 3).  Does the various runs give drastically different results? What is the main resultant source of randomness: the gradient descent; the Monte Carlo estimate of eq (9); or the sampling procedure in eq (5)?\n\n3. Section 5.3 says the best results from five runs are picked, but Table 3 shows standard deviation across three runs. Please clarify or correct.\n\n4. Figure 5 is not very useful without the actual data also plotted in the test region; and it is not referred to at all the the text. Suggest to remove this, and use the space to expand on how one may interpret the kernels learnt in section 5.2.\n\n5. Another future work could be to move beyond mean field approximation.\n\nChanges after rebuttal\n=================\n4. Figure 5 does indeed have the actual data plotted. It was not very clear, and I missed it in the first reading. Suggest the authors to make this clearer in the revised version.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "# Societal Impacts\nSuggest to include some possible impact on the modelling aspect. For example, climate/weather modelling which is actually one of the examples.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "g4whwewBwf", "reviewer": "Reviewer_5eAE", "summary": "This paper proposes to model volterra kernels with Gaussian processes and develops an inference algorithm based on variational inference. ", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\n\n1. square integrability of greens functions (line 66).\n\nWhat is required of G to ensure finite output depends on G. That is, if u is bounded then it is sufficient that G is integrable (BIBO stability). \n\n2. Decaying square exponential (line 100-101).\n\nSquare integrability of G is sufficient for Tobar [29] because it is convolved with a white noise signal. In the present situation it is not clear that this is works. For example, take the second term in the sum of eq (2):  \n\nint  ( int G_{d,2}(t-s_1,t-s_2) u(s_1) d s_1` )  u(s_2) d s_2 \n\nwhich can be made finite if  u is in L_2 and  int G_{d,2}(t-s_1,t-s_2) u(s_1) d s_1` )   is an L_2 function for every t, but is this ensured by square integrability of G? \n\n3. Typo (Line 157). \n\nLieber -> Leibler.\n\n\n4. Number of inducing points (Line 240).\n\nThe number of inducing points assigned to the filters appears counterintuitive. You need few points to cover a small space and many points to cover a large space but the number of points assigned to the filters is directly opposite to this?\n\n5. Length scale and noise parameters (line 249). \n\nIt is not made clear how the length scales are fixed or how the noise hyperparameters are post-hoc estimated. \n\n\n\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "yes.", "review_text": "", "extra_scores": {"ethical_concerns": "no. ", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2.5"}, "extra_text": {}}, {"review_id": "QlXfDD2kcT8", "reviewer": "Reviewer_o4oR", "summary": "The article proposes a methodology for designing covariances kernels in a nonparametric way using Volterra series. The ideas proposed build on existing works on convolution models, Volterra kernels, SysID and MOGP. The paper substantially differ from previous works and the contribution is validated empirically. ", "questions": "", "limitations": "", "rating": 6, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The authors deliver a thorough and clear literature review, which makes it easy for the reader to identify the motivation and the contribution of this article. As other works in the last decade, this paper focuses on kernel design for GPs, in particular, this kernel takes a nonparametric perspective, a practice that has done in the past but only by few researchers. The main distinction of this work is (in addition to the nonarametirc kernel construction) its application for multiple channels and the inclusion of nonlinearities\n\nIn my opinion, this paper is appropriate for NeurIPS but I would appreciate if the authors referred to the following points \n\n- The authors use a fixed grid of the input and kernel, please elaborate why this does not constitute a discrete/parametric model. \n\n- The authors refer several times to the similarities between the proposed method and GPCM. Is the approximate inference the same for both? (structured VI) or the implementation of GPCM in this paper is \"enhanced\" wrt to the original one? \n\n- The approaches for nonparametric kernels are known by their non-trivial implementation, perhaps the authors would make their code available so that the public can use their method. \n\n- Please elaborate on the computational complexity of the method, hopefully quantitatively. Due to the involved high-dimensional elements, this might be worth considering.\n\n- The proposed method is designed for non-Gaussian data. Then, in the experiments, it seems that the authors implicitly assume that real-world data is non-Gaussian (it's most likely the case). Why not explain which nonlinearities are modelled by the Volterra approach?  \n\n\n\nminor: \nl77: Voterra -> Volterra\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "yes", "review_text": "", "extra_scores": {"ethical_concerns": "NA", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers found the paper interesting and well written, with a thorough and clear literature review, and methodologically sound. All reviewers recommend acceptance, and the remaining concerns are possible to address in the final camera-ready version. I recommend that the authors carefully go through the manuscript in light of the reviewer comments.

Decision: Accept (Poster)