Heteroscedastic Gaussian Processes and Random Features: Scalable Motion Primitives with Guarantees
Authors: Edoardo Caldarelli
Antoine Chatalic
Adrià Colomé
Lorenzo Rosasco
Carme Torras
Venue: CoRL 2023 Poster (2023)
Keywords: Gaussian process regression
random features
motion primitives

Abstract: Heteroscedastic Gaussian processes (HGPs) are kernel-based, non-parametric models that can be used to infer nonlinear functions with time-varying noise. In robotics, they can be employed for learning from demonstration as motion primitives, i.e. as a model of the trajectories to be executed by the robot. HGPs provide variance estimates around the reference signal modeling the trajectory, capturing both the predictive uncertainty and the motion variability. However, similarly to standard Gaussian processes they suffer from a cubic complexity in the number of training points, due to the inversion of the kernel matrix. The uncertainty can  be leveraged for more complex learning tasks, such as inferring the variable impedance profile required from a robotic manipulator. However, suitable approximations are needed to make HGPs scalable, at the price of potentially worsening the posterior mean and variance profiles. Motivated by these observations, we study the combination of HGPs and random features, which are a popular, data-independent approximation strategy of kernel functions. In a theoretical analysis, we provide novel guarantees on the approximation error of the HGP posterior due to random features. Moreover, we validate this scalable motion primitive on real robot data, related to the problem of variable impedance learning. In this way, we show that random features offer a viable and theoretically sound alternative for speeding up the trajectory processing, without sacrificing accuracy.

Heteroscedastic Gaussian Processes and Random Features: Scalable Motion Primitives with Guarantees
Edoardo Caldarelli' Antoine Chatalic² Adrià Colomé' Lorenzo Rosasco 3,4 Carme Torras
Institut de Robôtica Informàtica Industrial, CSIC UPC, Barcelona, Spain MMaLGa Center DIBRIS Università di Genova, Genoa, Italy
CBMM Massachusets Institute of Technology. Cambridge, MA, USA *Istituto Italiano di Tecnologia, Genoa, Italy Correspondence to: ecaldarelli@irt www.iiii.ip...... www.uuc..... edu
Abstract: Heteroscedastic Gaussian processes (HGPs) are kernel-based, non- parametric models that can be used to infer nonlinear functions with varying noise. In robotics, they can employed for learning from demonstration as mo- tion primitives, i.e. model of the trajectories HGPs provide variance estimates around the reference signal modeling the jectory, capturing both the predictive uncertainty and the motion variability. How- ever, similarly to standard Gaussian processes they suffer from cubic complexity in the number of training points, due to the inversion of the kernel matrix. The un- certainty can be leveraged for more complex learning tasks, such as inferring the variable impedance profile required from a robotic manipulator. However, suit- able approximations are needed make HGPs scalable, at the price of potentially worsening the posterior mean and variance profiles. Motivated by these obser- vations, study the combination of HGPs and random features, which are popular, data-independent approximation strategy of kernel functions. In a the- oretical analysis, we provide novel guarantees on the approximation error of the HGP posterior due to random features. Moreover, we validate this scalable motion primitive on real robot data, related to the problem of variable impedance learning. In this way, we show that random features offer a viable and theoretically sound alternative for speeding up the trajectory processing, without sacrificing accuracy.
executed by the robot.
Keywords: Gaussian process regression, random features, motion primitives.
1 Introduction
Learning from demonstration (L.fD) is broadly used technique transfer skills from humans to robots in flexible and intuitive way [1]. Within the context of robotics manipulation, LíD consists of recording the movement of an arm performing a specific task multiple times. This data is then used to fit a model of the trajectory, called a motion primitive in this context, which allows the robot to reproduce the skill of interest. One popular way to achieve this objective is by means of the so-called Gaussian process (GP) regression [2]. Being a Bayesian model, a GP naturally provides a time-varying reference signal to be followed by the robot (the GP posterior mean), and an uncertainty quantification (the GP posterior variance) around the reference signal.
When used for LfD, GPs are usually heteroscedastic (HGP), i.e., the variance of the noise corrupting the recorded trajectories is not constant [3, 4]. As discussed by Arduengo et al. [5], such a time- varying noise variance is key asset for modeling motion variability when HGPs are used as motion primitives. The inconsistency of the human demonstrations, captured by the time-dependent noise variance, is form of epistemic uncertainty that is intrinsic the task, i.c. cannot be resolved by increasing the number of human demonstrations. Moreover, the HGP posterior variance quantifies the input-dependent uncertainty on the predictions, which is relevant when the testing points are far from the training data due to gaps in the human demonstrations.
7th Conference Robot Learning (CoRL. 2023), Atlanta, USA.

---

Figure I: The tasks involved our experiments. In the first two pictures, the human is guiding the robot perform: free motion above metal piece; the insertion of the end-effector in the trail inside the piece. the last picture, the robot is pulling a bed sheet to remove the wrinkles.
While being attractive for the aforementioned reasons, GP regression suffers from cubic scaling in the number of training samples as it requires an inversion of the kernel matrix. For this reason, several approximation techniques have been deployed, that allow for linear scaling in the number of observations. In particular, these approximations can be grouped into variational ones and the spectral ones. The former methods are data-dependent and rely on variational approximation of the posterior distribution of the GP [6, 8]. Spectral approximations, on the other hand, are most often inteeiinde,nn,nt,, and can rely, e.g-, on polynomial approximations 10], truncation of the kernel's spectrum [11], or randomized approximation of an integral representation of the kernel function. The latter technique includes the well known random fandomfeatuues (RFs) [12], that have successfully been used for GPs approximation [13, 14, 15, 16, 17]. In order to achieve a satisfactory model of the trajectory of interest, these approximations should not deteriorate the posterior mean and variance as these quantities fully specify the robot's motion; when such guarantees can be rigorously obtained, approximate GPs become an appealing motion primitive [11, 18, 19].
Contributions In this work, we study the combination of vector-valued HGPs with RFs for fast trajectory processing in LfD. In order assess the reliability of the approximation scheme, we perform theoretical analysis and propose novel bounds on the approximation error of the posterior mean and variance of the HGP. This is achieved by leveraging the deep connection between approx- imate GP regression and kernel ridge regression with RFs [20]. Moreover, given the well-studied relationship between the HGP posterior variance and the time-dependent robot stiffiness employed in variable impedance control (VIC) [21, 22, 23, 24), we further assess the quality of the approx- imation for the three VIC tasks described by Caldarelli et al. [24], and shown in Figure Overall, we demonstrate that the combination of RFs and HGPs is theoretically sound, and significantly improves the running time of the motion primitive fitting, without sacrificing accuracy.
2 Related Work
Motion primitives One of the most popular trajectory representations appearing in literature is given by the so-called probabilistic movement primitives (ProMPs) [25, 26, 27]. ProMPs model the trajectory at each time-step as a parametric weighted combination of basis functions, and are therefore different from HGPs, which are non-parametric by definition. ProMPs can be connected with another class of trajectory models, namely dynamic movement primitves, as shown by Li al. [28]. On the other hand, GP-based models exhibit strong analogies with the kernelized movement primitive model introduced by Huang al. [29]. In general, these approaches differ from GP regression as they also capture correlation between the different degrees of freedom (DOFs) of the trajectory. However, GPs can also be extended to account for output correlations by leveraging suitable coregionalization methods [30], based, e.g., on Gaussian mixture models [31].
Spectral approximations Some theoretical properties of spectral approximations for ho- moscedastic, scalar-valued GPs have been explored in previous works. Särkkä and Piché [9] propose a feature approximation based on different polynomial approximations of the RBF kernel function. They are able to prove uniform convergence of the kernel, mean and variance, but without conver- gence rates. Moreover, Solin and Särkkä [10] propose to compute the features based on the Fourier transform of the Laplace operator. They prove convergence of the kernel function values and con- vergence of the posterior mean and variance, but with dependency the size of the domain of the

---

expansion of the Laplace operator. RFs are also widely used in the Bayesian optimization setting, as they allow to provide scalable algorithms with theoretical guarantees [11, 19). In this scenario, the uniform convergence of the approximate kernel function is of interest as the optimization algorithms may require sampling points in the whole domain of the function being optimized. Compared to our work, as we will discuss in Section 4, these bounds have worse dependency on the number n of training samples (which is not due to the uniformity of the bounds). Although approximate GPs are known to suffer from the so-called variance starvation phenomenon, leading poor estimates of the posterior variance when moving far from the training data [32]. we stress that this is not prob- lematic in LfD where one is mainly interested in the interpolation of densely sampled trajectories, whose boundaries are delimited by the task duration.
Variational methods Multiple methods based on variational approximations of the posterior distribution of the GP have been proposed. They typically rely on the choice of inducing points, i.. of a subset of the data summarizing the whole training set [6, 7, 8]. One fundamental convergence result for the Kullback-Leibler (KL) divergence between the variational and exact GP posterior was proven by Burt et al. [18]. Moreover, Burt et al. [33, Proposition 1] showed that the error on the KI. divergence implies convergence of the variational posterior means and variances. However, their bound depends on the value of the exact posterior variance, which makes the comparison with our re- sult difficult, as it will be clear from Section 4. Hence, recent approaches such as variational Fourier features (VFFs) by Hensman et al. [8], which have shown good performance in practice, especially in the case of billions of datapoints, do not enjoy strong theoretical guarantees as RFs. Moreover, VFFs are limited to the Matérn class of kernels (and the work by Dutordoir et al. [34] provides an extension to stationary kernels on the sphere), while RFs be applied to any stationary kernel.
3 Background: Approximating HPGs with Random Features
HGP regression for trajectory encoding Let X' be an input space, and x: X - R be a scalar-valued function. LfD, the function x depends on time, i.e., X Rzo. A GP [2] specifies a prior distribution over the function z, which depends on mean function μ X - R and a kernel k: X x X R. We say that follows GP with mean fe and kernel function k if for any vector of time-steps t E R"" the vector x(t) of evaluations of x at t follows the multivariate normal dis- tribution x(t) ~ N (Jtt),, K), where K E Raxx is defined as Kú,- k(4,, The mean function is usually assumed to be O, as the function values are standardized. The kernel function depends on a set of hyperparameters, 6, which can be fixed or inferred from data. In LfD, the time-dependent trajectory to be encoded consists of DOFs. If all DOFs are fully observed, each of them is usually modeled by an independent GP [35]. This type of estimation can be linked Hilbert space (RKHS) of vector-valued functions [30]. We assume have access y R", vector of noisy measurements of x(t). As shown by Arduengo et al. [35]. assuming constant noise variance in LiD might severely limit the quality of the posterior prediction. Therefore, the noise is dime-dependent, turning the GP into an HGP. The posterior distribution of an HGP, conditioned on the noise variance values at the training points and the function's 's observations, is E R". Lastly,
reproducing kernel
analytically available. Let t" be a testing point, and k := [k((t,,t),.... "it,,tt))t,...t,,,,,, k(tn, t') let Emsise E R=x= be the noise variance at the training points, and 0² poise,t- be the noise variance at the testing point. The posterior mean and variance of the associated HGP are, respectively, [35]
Mpooa (1*) = (K + Emoise) 'y, op (1*) = k((", + 22 noise.t" k2- (K + Emcise) -k..
(2)
Random features Let (2, A, be a probability space over the sample space Ω and let s Ω x X' + R. Random features (RFs) are class of randomized methods that can be used to approximate a kernel function admitting an integral representation of the form
k t') I (Mw t) (oo(w,t))cc (w)
(3)
by dicretizing it. This is the case for standard kernels such as the ones belonging to the Matémn family or the radial-bas fasis-function (RBF) kernel, discussed by Rahimi and Recht
3

---

[12]. The expectation above is approximated by sampling m vectors (w))ısssm TW). For 6(t) _ (w),,)),,.... (wm,, t)]², the kernel is approximated as k(t, == 2(t) D(t). For instance for the case of a positive definite stationary kernel k(t, = k(t - t¹), Bochner's theorem ensures that &o has a non-negative Fourier transform, which can thus be used in place of т while defining 1 to be trigonometric function [12]. Rahimi and Recht [12] show that, for the RBF kernel k(t,t") RBF exp (-Ilt= the values of w are s.t. Wi (2 exp (-12w1?/2). Then, for b uniformly sampled from [0, 2r), they define ()w,() = V2RRAF cos(wl + b). To retrieve closed-form form expressions for the posterior mean and variance of an HGP approximated with RFs (in the following referred to as RF-HGP), we can define the matrix K E TRax with entries Ki,- (kt, {z), For any testing point t" we define k,- let Emoise and d noise, 2* be as in the previous paragraph. The posterior mean and variance of the associated HGP are, respectively [35]
ym
kit, (²}]" E R". Lastly,
Hpos(( = K. (K + Epoise) ''y,
(4)
(5) We can observe that these equations are structurally the same as Equations (1) and (2), with the approximated kernel & replacing the exact kernel k. Fixing the dimensionality of the RF vector, RF-HGPs allow for linear complexity in the number of samples when computing the posterior dis- tribution of the HGP, as the posterior expressions in Equations (4) and (5) can be simplified by the well-known Woodbury identity [36]. Such an identity allows to perform matrix inversion in ((m)), m < n being the dimension of the RF vector, as shown in Appendix A. Moreover, the kernel ap- proximation is inteeindde,n,,t,, as the functions (Ww),)) :) can be computed prior to seeing any data.
opoat (t"a k(t", (t) + o moise.t" K.f- (K' + -mniise
Quantification of the Approximation Errors
We now analyse the proposed approximate HGPs in two different settings. In the first one, the hy- perparameters 0 of the kernel function (which typically contains at least a lengthscale parameter). as well as the noise variance values at the training points, Enoise, and at a testing point, o mose,, are fixed priori (oracle setup). In the second setting, 0, Lnoise. and r are directly inferred from data (heuristic setup). In the first case, the error introduced by the RF approximation can be rigorously quantified, as we will show in the remainder of this section. In the second case, convergence to the exact HGP can be attained empirically, and suitable inference algorithm needs to be used, as we will discuss in Section Both in the oracle and in the heuristic setup, the RF-HGP's posterior dis- tribution is computed at set of testing time-steps T = { ........ t. which are fixed a priori based on the desired task duration. In this section, we derive error bounds on posterior mean and variance for the RF-HGP, at the test points in T, within an oracle setup. These results guarantee that RFs are a suitable approximation strategy, that can be used in robotics scenarios. The proofs, reported in Appendix B.2 and Appendix B.3, rely on techniques developed by Rudi and Rosasco [20].
4.1 Assumptions
In the following, we set (-) = vilw, -). where the latter function is defined in Section 3. More- over, we make the following assumptions.
Assumption 4.1. The kernel is bounded, i.e., k(t,t') S Moreover, the kernel admits the integral representation of Equation (3), in terms of suitable function (ww,-).
Assumption 4.2. The entries of the RF vector are bounded, i.c., (w) S a, Vj E {l...., m). Vt. Assumption 4.3. The noise variance values in the diagonal matrix Enoise are bounded from below, i.e. By E (0, 1] s.t. (Frmise)iz 2 7M, Vi E n.).
Assumption 4.1 is standard, and holds, e.g., for the RBF kernel (k = R Assumption 4.2 holds, e.g., for the RFFs described in Section 3 (with a (2rRuF). Assumption 4.3 is obtained by replac- ing, c.g, K + Enoise in Equations (1) and (2) with K (Eioise + ynl).. common practice aimed at avoiding numerical instabilities in matrix inversion. Lastly, for the RBF kernel and RFFs, a²,
4.2 Concentration of Posterior Mean and Variance
Under these assumptions, derive the two main results of this paper, bounding the deviation of the posterior mean and variance of our RF-HGP w.r.t. its exact counterpart.

---

RF-EGP
Boxd
RF-ROP
Roxt HGP
Pood Ene [+]
RP-HCP
Poue ICP
Prod tax
0.11
0.30
0.1
Poo.
0.1
0.05
50
100
1.50
200
Numbrr
Nomber
RPF
Nimber of RPFS
Figure 2: Time taken to compute the posterior distribution with an exact HGP and an RF-HGP, for free motion (left), assembly (center) and bed-making (right) tasks, in an oracle setup. Median, the 15th and 85th percentiles across all DOFs, 50 seeds.
Theorem 4.4. Let E (0, and m. be the dimension of the RF vector. Let m 2 8 ( a² a Mpost and Poon denote its exact and ft--aproxximted posterior mean of Equations (1) and (4). Let v maxd v. where y, denotes the observations DOF i. Lastly, let T be set of testing points. Under Assumptions 4.1 to 4.3, with probability at least - dô, Vtt E T, it holds that
log( (388)))) and consider a vector-valued HGP with d independent components. Let
2d 2d 2
V2dKv VT
Suc 2 log
a1+02/17
2 log
8x2 48"
(toor(t)) (t )l12 Sva?
76
+ VaV
myy
3m
ym
Theorem 4.5. Under the assumptions of Theorem 4.4, denoting pou and o² post respectively the exact and RF RF-appooximated posterior variances of Equations (2) and (5), holds with probability at least 1 - dö, vt* E TT,
2171 2d
+ ))
2d log 2 8
looporl()) - ppoan (1))2 /
m
a. 22
8x2 2 log
(1 a² //7)
log Sk2
78
+ =av
70
v7
3m
7mm
Our bounds show that both the mean and the variance errors for vector-valued RF-HGP scale as (m-11/2). Our results cover the RF-based approximation of the homoscedastic GP special case taking Moise = 31 for some 3 O. Our results ensure that RFs are viable data-independent approximation method, and match the rates of other data data-dependent strategies, such as the well- studied Nyström method [37] for GPs, as shown by Lu et al. [38]. Furthermore, by inspecting the proofs, and in particular the results in Appendix C, one can observe that the assumptions of the LiD framework (fixed testing points, scalar domain) translate in a dependency of the bounds on [T]. Nonetheless, our proofs can easily be adapted to the case of dense, multidimensional domains by replacing Corollary C.2 with the uniform convergence result for RFs [39, 40] (rather than using point-wise bound and union bound), which would leave the overall error in O(m -1/2).
As mentioned in Section 2, our bounds strictly improve over the state-of-the art bounds by Mutny and Krause [11]. For instance, Mutny and Krause [11, Proof of Theorem 5 and Proposition 1] in the context of homoscedastic GP approximated with deterministic features provide bounds of the form sup;- [Mpos((")) po (omma((())
(1*) ] S e and supt- lo (1*) op S en -2 where denotes the accuracy of a uniform bound on the kernel approximation, i.e. supx.y [x(2, y) ie(z, y) S In our setting, when d = 1 and in the homoscedatic setting Emoise = 021 and 7 = /n (cf. Assumption 4.3) we obtain dependence in n of order ((13/2 log both for the posterior mean and variance which strictly improves on these bounds. This result can be derived, e.g., for the mean by observing that the rate with respect to y translates to an error upper bound with complexity O(nyTogn) Onn3/2 l) + O(nylogn)
5 Empirical Evaluation
In this section, we empirically validate the RF-HGP. The implementation built on the package gpflow by Matthews et al. [41]'. We choose RFFs to approximate the RBF kernel, as reported
'The open-soure code for the experiments is available at https: http//ggttttttttete//e httpcchww..ooo.///.ffffetteeeeee
5

---

Free motion
Assembly
Ded-making
Free motion
Assembly
Ded-making
RMSEE var:
102
HMSES
104
10
100
50
100
150
200
250
50
100
150
200
250
Number of RFFs
Number of RFFs
Figure 3: RMSE between the posterior means (left) and posterior variances (right) of an HGP and an RF-HGP, in oracle setup, for three different tasks. Median, 15th and 85th percentiles across all DOFs, 50 seeds. The purple curves are the theoretical rates including the dependency on v, y and m, and show that the overall rate in 1/Vm matches experimental result. The lowest purple curve is for the bed making task, while the rates for assembly and proof-of-concept coincide.
3T-HGP
Nyström- HGP
RMSSE: moan
RF-HGP
Nystcim HGP
RF DE-HGP
Nystetun HGP
10P
RMSSE var.
RMS
200
20-
20
10
60
90
20
10
90
80
20
10
60
80
Namaber of features
Number of Seatures
Number of featurees
Figure 4: Incremental learning with missing chunks of data, in an oracle setup, for free motion (left), assembly (center) and bed-making (right) tasks. RF-HGP vs Nyström approximation with fixed centers, to achieve the same complexity as RFs. Median, 15th and 85th percentiles across all DOFs, 50 seeds.
in Section 3. In this work, we do not train the RFF parameters w's, as it is prone to overfitting [13, 42], and annuls the Monte-Carlo interpretation of RFFs introduced in Section 3. To assess the correctness of the GP approximation, we process real demonstrations of different robotic tasks, obtained by means of kinesthetic teaching with a 7-DOF Barrett WAM manipulator. The trajectories are recorded while performing a movement in free Cartesian space, an assembly task, and a bed- making skill [24]. These trajectories are partieularly interesting from the VIC point of view, and summarized in Figure 1. Since the time-depender impedance, stiffness, of manipulator can be tuned based on the variance of the human demonstrations [21, 23, 22, 24], it is important that the RF approximation does not deteriorate the quality of the posterior distribution of the HGP. The first task exhibits varying boundaries in which the manipulator can move, while the latter two involve physical constraints on the robot's motion. The assembly task requires a motion with low stiffness when the pieces being mounted are contact. On the other hand, the bed- bed-making skill requires the robot to stiffen up remove wrinkles from the sheet's surface, in spite of it opposing the motion of the robot's end-effector. Each task uses different number of human demonstrations, namely 6, 7 and 5. While the number of demonstration is relatively small, the total number of training points is 1286 for the -ffecccet-ppttt experiment, 1222 for the assembly task, and 864 for the bed-making skill, justifying the need for scalable GPs. Moreover, the number of samples might increase for longer or denser trajectories. We consider an oracle setup in Section 5.1, where all hyperparameters and noise variances are assumed be known, and then discuss and evaluate Section 5.2 the heuristic case of unknown hyperparameters and noise variance.
5.1 Oracle Setup
As start, we consider the case in which the GP kemel hyperparameters 0 are known, and the noise variance values Enoise, o are provided by an oracle, as discussed in Section 4. This experiment is useful to assess how the prediction capability of RF-HGP deteriorates due to the kernel spectral approximation, in view of the theoretical analysis carried out Section 4. The oracle is given by an exact HGP trained with the EM algorithm by Kersting al. [4]. for 15 iterations per DOF.
Offline learning In this setup, we are interested in comparing the prediction times (Figure 2) and accuracy (Figure 3) of exact and approximated GP regression. By prediction time, we refer here the time taken compute the posterior means and posterior variances at the testing points. The human demonstrations are temporally aligned on the same interval, and processed all together, in an offline fashion. The prediction times, as shown in Figure 2, indicate that the RF-HGP is more

---

10030 toom
RFF HGG
SVOP EEP
N00 1900
BFF-CPP
SVCP HSP
1000
REF HGP
SVOP HGP
10"
100
50
150
Figure 5: RMSE between the posterior means of an exact and an approximated HGP, in heuristic setup, for free motion (left), assembly (center) and bed-making (right) tasks. Median, 15th and 85th percentiles, across all DOFs and 5 seeds.
REF-HGP
SUGINFGP
RFF-IGP
SUGP-RGP
RFF-RGP
SUGP-HCP
10000g Mr.
Figure 6: RMSE between the posterior variances of an exact and an approximated HGP, in heuristic setup, for free motion (left), assembly (center) and bed-making (right) tasks. Median, 15th and 85th percentiles, across all DOFs and 5 seeds.
efficient, as expected. Moreover, to asses the predictive accuracy, we compute the normalized root- mean-squared-error (RMSE) between the posterior means or variances of the HGP and RF-HGP:
EUeT(Wrost,ll - ppost, (())² Ever :p)2,
RMSE Venan
100.
(6)
The RMSE& on the variance is computed in the same way. This experiment allows us to observe that the error rate follows the expected rate in 0mm2)/2), as shown in Figure 3, indicating that RFs can be used speed up the inference without sacrificing accuracy in the posterior's calculation. We can also observe that although all errors decrease with increasing number of features, the assembly task exhibits the largest eтror among the three tasks. This is likely to be due to the RF approximation being challenged by the sharp changes in the task variability.
Incremental learning RFs are data-independent, and they can be used to perform incremental learning of the posterior distribution of the HGP [43]. In an oracle setup, the posterior distribution can be updated every time a new demonstration is gathered (see Appendix D). This approach also possible, e.g., with the Nyström approximation of the kernel function, provided that the Nyström centers are sampled from the first demonstration and not updated afterwards. Doing so would imply re-computing the correlation matrices between the centers and the whole set of training points, which has linear complexity in the size of the training Nyström centers is not issue in general for LID, the domain is scalar and all demonstrations are temporally aligned on the time interval [0, Ij. and thus the first demonstration already covers the whole domain. However, this fact does not hold if, for any reason, portions of the human demonstrations are missing or invalid. In this case, the initialization of the centers might be poor, and the posterior calculation spoiled. This issue shown the examples of Figure 4, where chunks of 60 observations, sampled uniformly at random, were removed from each demonstration. RFFs offer a more reliable option, such a scenario.
Keeping fixed set of
5.2 Heuristic Setup
After validating our theoretical results, we now consider the more realistic setting in which the kernel hyperparameters and noise variance are unknown and need to be heuristically estimated from data. Note that this step can serve as preliminary stage with a few demonstrations before applying, e.g., the incremental learning described in the previous section. Here, we train the RF-HGP as follows.
Workflow with hyperparameters tuning Here, we consider the poneectttionmmxiiizzntte training proposed by Kersting et al. [4], adapted to the RF-HGP. The algorithm comprises the follow- ing steps, which are performed independently for each DOF (all the GPs involved take time as input);
(EM)
1. train a first homoscedastic GP (GP1), approximated with RFs (RF-GP), on the demonstrations' data with a maximum likelihood estimate (MIF),, to retrieve 0:
7

---

RP HGP
SVOP HOP
2xre HSP
Opi Fie H
RE HGP
SVOP HOP
Exsot HGP
RI HGP
SVGP HGP
Beact HGP
Opt. time
Ogt one
100 TMPPs
L50
Number
Nemror
RtF
Number RTDs
Figure 7: Time taken to perform an MLE with an exact HGP and an RF-HGP, with demonstrations of free motion (left), assembly (center) and bed-making (right) tasks, Median, the 15th and 85th percentiles, across all DOFs and 5 seeds.
compute the mean-squared-error (MSE) between the posterior mean of GP and the demonstra- tions' data, at each time-step;
3. train second RF-GP (GP2) by means of MLE, using the MSE training data; GP2 is surrogate model of the time-dependent noise variance function;
4. compute the posterior mean and variance of a new RF-HGP (GP;) as in Equations (4) and (5), with the kernel hyperparameters from step and the noise variance values in Enose and o2 moise,," given by the posterior mean of gP2, evaluated at the training and testing time-steps respectively; compute MSE between the posterior mean GP3 and the demonstrations, at each time-step; repeat from step 3 until convergence or until the maximum number of iterations is reached.
Sampling first whole noise variance profiles from the posterior of GP2. and then whole trajectories from the posterior specified by Equations (1) and (2) conditioned on the noise variance, yields trajectories that are strongly correlated in the regions of low variance, and vice-versa.
Results Figure 7 reports the time taken to complete an MLE for retrieving the GP hyperparame- ters, both in the exact and in the RFF setup. Considering the errors reported in Figures 5 and 6, we observe that convergence to the exact HGP posterior can be heuristically attained in the worst possi- ble scenario of having no knowledge about the necessary GP priors. Again, the largest errors are at- tained by the assembly task, as discussed for the oracle setup. Concerning variational methods, such as the sparse variational GP (SVGP) by Hensman et al. [44], Figure 7 shows that SVGP, employed in the EM training algorithm, requires many more iterations per step to converge, in the hyperparam- eter training, due to the greater complexity of the optimization problem being solved. To overcome this issue, we set the maximum number of iterations per step to 100. However, as we can observe from Figures 5 and 6, this choice hinders the convergence through the EM training process, and our complexity trade-off.
proposed feature-based approach has a strictly better
Limitations
The assumptions behind our theoretical analysis were stated in Section 4.1. Moreover, as discussed in Section 4, theoretical analysis focuses on an oracle setup. This type of scenario is standard in kernel theory [45, 20, 38], as it allows to decouple the error due the kernel approximation from the uncertainty surrounding the GP hyperparameters. If this oracle setup does not hold, the HGP training method is heuristic, and uses two approximate GPs iteratively in an EM fashion. The errors displayed in the heuristic experiments of Section 5.2 may change with different training algorithm.
Conclusion
In this work, we have studied the combination of heteroscedastic Gaussian processes and random features, used as scalable motion primitives in the context of learning from demonstration. In theoretical analysis, we derived novel upper bounds on the approximation error, induced by random features, on the posterior mean and variance of a Gaussian process. Moreover, we have validated this approximate motion primitive w.r.t. relevant tasks for variable impedance control of robotic manipulators, namely, a motion in free Cartesian space, an assembly task, and a bed-making skill. Our theoretical and empirical results demonstrate that random features are theoretically sound approximation method, that can be used to speed up the motion primitive fitting without sacrificing accuracy. Moreover, have shown that random features are well suited incremental learning from demonstration, thanks to their fata-independent nature.
8

---

Acknowledgments
E. Caldarelli, A. Colomé and C. Torras acknowledge support from the project CLOTHILDE ("CLOTH manlpulation Learning from DEmonstrations") funded by the European Research Coun- cil (ERC) under the European Union's ' Horizon 2020 research and innovation programme (Advanced Grant agreement No 741930). E. Caldarelli acknowledges travel support from ELISE (GA No 951847). L. Rosasco acknowledges the financial support of the European Research Council (grant SLING 819789), the AFOSR projects FA9550-18-1-7000 FA9550- 17-1-0390 and BAA-AFRL- AFOSR-2016-0007 (European Office of Aerospace Research and Development), the EU H2020- MSCA-RISE project NoMADS DLV-777826, and the Center for Brains, Minds and Machines (CBMM), funded by NSF STC award CCF-1231216.
References
[1] A. G. Billard, S. Calinon, and Dillmann. Learning from humans. robotics, pages 1995-2014, 2016.
Springer handbook of
[2] C. K. Williams and C. E. Rasmussen. Gaussian processes for machine learning, volume 2.
MIT press Cambridge, MA, 2006.
[3] P. Goldberg, C. Williams, and C. Bishop. Regression with input-dependent noise: Gaussian
process treatment. Advances in Neural Information Processing Systems, 10, 1997.
[4] K. Kersting, C. Plagemann, P. Pfaff, and W. Burgard. Most likely heteroscedastic Gaussian process regression. In International Conference Machine Learning (ICML), pages 393- 400. PMLR, 2007.
[5] M. Arduengo, A. Colomé, J. Lobo-Prat, L. Sentis, and C. Torras. Gaussian-process based robot learning from demonstration. Journal of Ambient Intelligence and Humanized Computing. 2023.
[6] H. Liu, Y.S.. Ong, X. Shen, and J. Cai. When Gaussian process meets big data: A review of scalable GPs. IEEE Transactions on Neural Networks and Learning Systems, 31(11):4405- 4423, 2020.
[7] M. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Interna- tional Conference on Artificial Intelligence and Statistics (AISTATS), pages 567-574. PMLR, 2009.
[8] J. Hensman, N. Durrande, A. Solin, et al. Variational Fourier features for Gaussian processes. The Journal of Machine Learning Research, 18(1):5537-5588, 2017.
[9] S. Särkkä and R. Piché. On convergence and accuracy of state-space approximations squared exponential covariance functions. In IEEE International Workshop Machine Learning for Signal Processing (MLSP), 2014.
[10] A. Solin and S. Särkkä. Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Computing, 30(2):419-446 2020.
[11] M. Mutny and A. Krause. Efficient high dimensional Bayesian optimization with additivity and quadrature Fourier features. Advances in Neural Information Processing Systems, 31, 2018.
[12] A. Rahimi and B. Recht. Random features for large-scale kemel machines. Advances in Neural Information Processing Systems, 20, 2007.
[13] M. Lázaro-Gredilla, Quinonero-Candela, C. E. Rasmussen, and A. R. Figueiras-Viche Sparse spectrum Gaussian process regression. The Journal of Machine Learning Research, 18:1865-1881, 2010.

---

[14] G. V. Karanikolas, Q. Lu, and G. B. Giannakis. Online unsupervised learning using ensemble Gaussian processes with random features. In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3190-3194, 2021.
[15] Y. Pan, X. Yan, E. A. Theodorou, and B. Boots. Prediction under uncertainty in sparse spectrum Gaussian processes with applications to filtering and control. In International Conference on Machine Learning (ICML), pages 2760-2768. PMLR, 2017.
[16] K. D. Polyzos, Q. Lu, and G. B. Giannakis. Graph-adaptive incremental learning using an ensemble of Gaussian process experts. In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 5220-5224, 2021.
[17] K. Cutajar, E. V. Bonilla, P. Michiardi, and M. Filippone. Random feature expansions for deep Gaussian processes. In International Conference Machine Learning (ICML), pages 984-893, PMLR, 2017.
[18] D. Burt, C. E. Rasmussen, and M. Van Der Wilk. Rates of convergence for sparse variational Gaussian process regression. In International Conference on Machine Learning (ICML), pages 862-871. PMLR, 2019.
[19] A. Dubey, No-regret algorithms for private Gaussian process bandit optimization. In In- ternational Conference on Artificial Intelligence and Statisties (AISTATS), pages 2062-2070. PMLR, 2021.
[20] A. Rudi and L. Rosasco. Generalization properties of learning with random features. Advances in Neural Information Processing Systems, 30, 2017.
[21] S. Calinon, Sardellitti, and D. G. Caldwell. Learning- based control strategy for safe human- robot interaction exploiting task and robot redundancies. In IEEE International Conference on Intelligent Robots and Systems (IROS), pages 249-254, 2010.
[22] L. Peternel, T. Petrič, and J. Babič. Robotic assembly solution by human-in-the-loop teaching method based on real-time stiffness modulation. Autonomous Robots, 42:1-17, 2018.
[23] D. Parent, A. Colomé, and C. Torras Variable impedance control in cartesian latent space while avoiding obstacles in null space. In IEEE International Conference on Robotics and Automation (ICRA), pages 9888-9894, 2020.
[24] E. Caldarelli, A. Colomé, and C. Torras. Perturbation- stiffieess inference in variable impedance control. IEEE Robotics and Automation Letters, 7(4):8823-88) -8300 2022.
[25] A. Paraschos, C. Daniel, J. R. Peters, and G. Neumann. Probabilistic movement primitives. Advances in Neural Information Processing Systems, 26, 2013.
[26] G. Maeda, G. Neumann, M. Ewerton, R. Lioutikov, O. Kroemer, and J. Peters. Proba- bilistic movement primitives for coordination of multiple human-robor collaborative tasks. Autonomous Robots, 41:593-612, 2017.
[27] S. Gomez-Gonzalez, G. Neumann, B. Schölkopf, and Peters. Adaptation and robust learning of probabilistic movement primitives. IEEE Transactions on Robotics, 36(2):366 2020.
[28] G. Li, Z. Jin, M. Volpp, F. Otto, R. Lioutikov, and G. Neumann. Prodmp: A unified perspective on dynamic and probabilistic movement primitives. IEEE Robotics and Automation Letters, (4):2325-233) 2023.
[29] Y. Huang, L. Rozo, J. Silvério, and D. G. Caldwell. Kernelized movement primitives. The International Journal of Robotics Research, 3837):833-852, 2019.
[30] M. A. Alvarez, L. Rosasco, N. D. Lawrence, al. Kernels for veltor-valued functions: A review. Foundations and Trends in Machine Learning, 4(3):195-266 2012.
10

---

[31] N. Jaquier, D. Ginsbourger, and S. Calinon. Learning from demonstration with model-based Gaussian process. In Conference on Robot Learning, pages 247-257. PMLR, 2020.
[32] Wilson, V. Borovitskiy, A. Terenin, P. Mostowsky, and M. Deisenroth. Efficiently sampling functions from Gaussian process posteriors. In International Conference Machine Learning, pages 10292-10302. PMLR, 2020.
[33] D. R. Burt, C. E. Rasmussen, and M. Van Der Wilk. Convergence of sparse variational in- ference in Gaussian processes regression. The Journal of Machine Learning Research, 21(1): 5120-5182, 2020.
[34] V. Dutordoir, N. Durrande, and Hensman. Sparse Gaussian processes with spherical har- monic features. In International Conference on Machine Learning, pages 2793-2802, PMLR, 2020.
[35] M. Arduengo, A. Colomé, J. Borràs, L. Sentis, and C. Torras. Task-adaptive robot learning from demonstration with Gaussian process models under replication. IEEE Robotics and Au- tomation Letters, 6(2):966 2021.
[36] M. A. Woodbury. Inverting modified matrices. Statistical Research Group, 1950.
[37] C. Williams and M. Seeger. Using the Nyström method speed up kernel machines. Advances in Neural Information Processing Systems, 13, 2000.
[38] X. Lu, A. Rudi, E. Borgonovo, and L. Rosasco. Faster kriging: Facing high-dimensional simulators. Operations Research, 68(1):233-249, 2020.
[39] B. Sriperumbudur and Z. Szabó. Optimal rates for random Fourier features. Advances in Neural Information Processing Systems, 28, 2015.
[40] F. Liu, X. Huang, Y. Chen, and J. A. Suykens. Random features for kemel approximation: A survey on algorithms, theory, and beyond. IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(10):7128-77 7148, 2021.
[41] A. G. d. G. Matthews, M. Van Der Wilk, T. Nickson, K. Fujii, Boukouvalas, León- Villagrá, Z. Ghahramani, and Hensman. Gpflow: A Gaussian process library using tensor- flow. The Journal Machine Learning Research, 18(40):1-6, 2017.
[42] Y. Gal and R. Turner. Improving the Gaussian process sparse spectrum approximation by rep- resenting uncertainty in frequency inputs. In International Conference on Machine Learning (ICML), pages 655-664. PMLR, 2015.
[43] R. Camoriano, Traversaro, Rosasco, G. Metta, and F. Nori. Incremental semiparametric inverse dynamies learning. In IEEE International Conference on Robotics and Automation (ICRA), pages 544-550, 2016.
[44] J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. In Uncertainty in Artificial Intelligence, page 282. Citeseer, 2013.
[45] A. Rudi, R. Camoriano, and L. Rosasco. Less is more: Nyström computational regularization. Advances in Neural Information Processing Systems, 28, 2015.

---

Table I: Summary of notations.
Variable
Meaning
Number of training samples
Number of testing samples
d
Degrees of freedom (one HGP for each)
ar
We (L) S a (ef. Assumption 4.2)
k(t, t') S (cf. Assumption 4.1)
max1 ما V
Emnice E RRax
Diagonal time-varying noise variance at training points
R = Enoise/" n
Normalized noise variance at training points
quise,t" E R
Noise variance testing point /*
7. noise, n/n
Normalized noise variance at testing point t*
0 0<<<1
R₁ > Sis
S: R2n
Sampling operator with normalization n 112 Normalized Gram matrix with exact kernel
== SS = K/nn E g=
LR E RRax
LR = L
/ (t) (/y.tt)
Element of approximate feat. vector Approximate feat. vector
@(t = m 11 (t),.... ..... e (t))²
Sm : 1.m + 1Rn
Sampling operator with normalization n 1/2 Normalized Gram matrix with RF kernel LmnRR Lun R
Lm = SmS = K k/nc E R
LmM.R Raxn
A HGP Posterior Equations Revisited
In this section, we will rewrite the exact and approximated HGP posterior equations from Sec- tion 3 in terms of standard linear operators used in RKHS theory. For a linear operator A, denote its adjoint by A". Let H be the RKHS associated to the kernel of interest. In denote S H + Rn the sampling operator de- order to retrieve suitable expression, Moreover, the adjoint of the sampling operator is de- fined Sʃ = (. ((mn))
fined as 5* : Rd - H S*a to E a;k((4,)), :), as being the i-th entry of a. Now, let R + R SS Note that K nL. Let = E let rz- a
Lo2 -iise, Lastly,
let (m denote the inner product of n-dimensional vectors. With this notation, let us consider single DOF of the trajectory be processed. The posterior mean of the associated exact HGP from Equation (1) is
Mpuaa (t*) -1 - (L++R) , vā
R"
Moreover, the posterior variance from Equation (2) is given by
po (1*) = k((",, (") nrr- (Sk(1" ) (L + Ri) Ski(1, ,))=..
(8)
Considering RFs, we can define the operator Sm : m - Sm = v (t((t) (T, and
Lm R" + R", Lm SmS With this notation, let us consider a single DOF of the trajectory to be processed. The RF-based posterior mean of the associated HGP from Equation (4) can be rewritten as
Ppou(t") RR
y
(9)
On the other hand, the RF-based posterior variance from Equation (5) is given by
opose (t*) = kett. t*) + nrt- - (Smm(tt) (Lm + R) -1 Sm ((t)))
(10)
A summary of the main operators and constants that will appear in the proof's can be found in Table 1.
12

---

Fast matrix inversion By defintion, the operators Lm and Sm are matrices. The inversion of the matrix I'm + R appearing in Equations (9) and (10) can be performed by means of Woodbury identity [36], as follows:
L L st S5 St + R))
(11)
= - Sm (1 + S
R-
(12)
The latter expression involves inverting posterior calculation if m < n.
matrix, which boosts the speed of the HGP
B Proofs of the Main Results
In this appendix, we report the proofs of the two main theoretical results of our paper, along with some technical propositions that will be extensively used. In the following, we denote by Ar the operator A + R. with R diagonal positive definite matrix, and by A, the operator A+yI. Moreover, in the remainder, ||-|| denotes the operator norm, while I|-I12 denotes the Euclidean norm a vector.
B.1 Useful Propositions
In this part, we report three propositions that will be useful in the proofs. Proposition B.1 (Proposition 8 of [20]). Let H. be a separable Hilbert space, A, B be two bounded self-adjoint positive linear operators on H. and X 0. Then
|| 1/2 /2 - |
1/2 B
1 (1 -)//22
(13)
where
B 1/ (B- (B-)) 1/2 17
8= Amaa
(14)
Proposition B.2. Let Sm 1 of the RF vectors are bounded, that is, [t)( a, Vj E fl..... m) Then,
R" Sm 1 (1((1t)) .... (t.n)) and assume that the entries
v
|| a.
(15)
Proof. The result follows from the definition of operator norm:
||Sml| =
sup aEle__lallsS1
||m
(16)
sup Dallesll Vn V (5(11), +.... (6(L)))
(17)
Vna= = a,
(18)
reported in the statement.
Proposition B.3. Let A be bounded positive semi-definite operator, and let AR := A+ R, with R diagonal positive definite and A, A + T.. Lastly, assume all entries in R are greater or equal to
7,. Then.
||
1/2 1//2 11/1111111 51
(19)
Proof. Noting that AR Ay (R- 71) /) x by hypothesis, it holds A, & Ag and thus
/2 A!/2122 = || 1/2 A 1/2 < |4| =1
(20)
13

---

B.2 Proof of Theorem 4.4 (Deviation of Approximate Posterior Mean)
We report here the proof of Theorem 4.4. We start by considering a single DOF, and generalize to a d-valued GP at the end of this section. We begin by proving a lemma that will be used to retrieve the main result.
Lemma B.4. Let m ( } + ( log(
80 70
and (0, 11. Then, the following bound holds, with
probability at least 8,
V2x V7
log
8x² 76
(1 + (1+a2/1) 02/22)
Scc2 2 log 4p a²
[1 Ln))s(((",
(21)
3m
ym
Proof. In order to bound the term of interest, we can use the fact that, for any invertible matrices A
and B, A B- = AI - AB-1) = ( - A)B-1, and Proposition B.3, as follows:
[1m,,1 [. )Sk(tt" : 1)12
(22)
= -1 12 ( RR L.m R LR Sk((t" 1ll
(23)
1 1/2 1-12 1
mm, R
mm, mm.
1/2
mm.7
1/ 1//2/1110 b.n 1/2 Imm) 1 Lp 5k((", 1)
1/2 1,
(24)
1/1/
1/2
1/2 /
1/2 - .mm) (k((",)) 1))
(25)
mm.y
1/2 || /
(26)
We can now proceed to bound each of the three factors. To start off, us consider |r.p. |/.-1//51, term can be bounded by using the polar decomposition of the bounded linear operator S, as follows. Let S (SS*)1/2 /?U, where U is a partial isometry. By Proposition B.3, the definition of polar decomposition, and by considering that T. & T.,, by definition,
S||. This
I R
1/2 /|| || 1/2 (SS*))
(27)
- || 1/2 /|
(28)
< ||
||
(29)
$
(30)
Now, we can move on to bound definition,
||- ||2 - Im L 1/21 11
To do so, we can observe that, by
Lm Sr
(31)
1 m (t()t))
(t1lt
x
(32)
72 m
(t,)
(tx))
Moreover, due to linearity of expectation,
E [Tmm] = L.
(33)
We can therefore apply Proposition C.4, with m. Q L, and Qp = L.m Note that Tr L is the trace of the normalized Gram matrix 1 K and hence is smaller or equal to x² under Assumption 4.1.
72 Lastly, the value of the constant Fox (7) in Proposition C.4 can be computed as follows:
(tt)1)
1
a
S
(34)
Vn
(t,))
Vn
(t..)
R
Thus, obtain, with probability least 1- 8.
log
SN² (1 + a² a2/77)
2 log 70
Sic
a
|L,1/2(L ,,
1/211
70
(35)
3m
7mm
14

---

To conclude the proof, we can bound |mmm. 1//21 By Proposition B.I, we have that
/4 1
where 3 = Amax
( L
(36)
(1 3)1/2'
According to Equation (33), we can apply Proposition C.4 and see that with probability least 1-
2 log Sx.2
2 log Sez2 a²
B <
S 0.5
(37)
3m
ym
provided that m 2 a²) ))
log(
8ca 78
Proof of Theorem 4.4 In order to retrieve the main concentration result, we can consider the fol- lowing decomposition of the error on the posterior mean. By Cauchy-Schwarl inequality and Equa- tions (7) and (9),
( 1,-1 m. Smm( - L. Sk(t", ))
1/ (1*) 1'posa (t.) =
(38)
vn
R2
1|1
Lm :)) L L,, L Sk(l"
)
(39)
S v|| R 0(0) - Sk((t, 1))l 1 + v ( mn L )Sk(tt
(40)
/// Shi(1, 111 ~ || L. S)"("
(41)
Now, upper bound the two norms appearing in the expression above. The first addend can be directly bounded by applying Corollary C.3. The second addend in Equation (41) can be bounded by Lemma B.4. Hence, we obtain the following bound with probability at least
|Apost (t") Mposs (f') S // 1/>1lSm Sk((t, ,)122 + vil(L L )S)((" 1))22
(42)
2v2a4 log my?
2Tn V2Ku V1
2 log Sx² (1+ a²/7)
Sec 2 log
a²
8
+
76
3m
mm
(43)
The final result for the vector-valued GP can be obtained by applying a union bound.
B.3 Proof of Theorem 4.5 (Deviation of Approximate Posterior Variance)
In this section, we prove our result related to the concentration the approximate posterior variance. Again, we begin by stating some lemmas that will be used the proof.
Lemma B.5. Let (0, 1). Then, the following bound holds, with probability at least
2n log
27n 8
(S)(1",)) - Sm1 Sk((t,(t S
(44)
ym
Proof. By Cauchy-Schwarz.
((Sk(t",)) Smm(tt), Sk(t", ))2= [|S((t")) Sm Smt sk(t", (45)
15

---

By using the polar decomposition of S, for suitable partial isometry operator U, and according to Propositions B.I and B.3
| Sk(t", 12 - (
(46)
(47)
'n 1/2 T 11
(48)
| 1/2 11 2 1 1/2 1/2
(49)
14 1/22 1/
(50)
(51)
VT
To conclude the proof, we can observe that, according to Corollary C.3,
2x log nTn 6
-.m S
(52)
Lemma B.6. Let (0, 1]. Then, the following bound holds, with probability at least
a²
2a log 2Tn
(Smp((1)) I (Sk(tt,) (tm)(("))))))
(53)
m
Proof. By Cauchy-Schwarz inequality and Proposition B.2,
((Smp(t)) I. (Sk((,)) -mn s
- Sm)((")) (54)
1s ( S (55)
a² [SK(((,))) S. (11))11
(56)
Now, we can again observe that, according Corollary C.3,
2a log 2Tn
|[Sk(t*.)) Sm 5mt)
(57)
which concludes the proof.
L
Lemma B.7. Let m 2 8 + 2 log( So² probability least 8. 78 ) and ā = =(0, 1). Then, the following bound holds, with
22/2
log
Sk² 78
(1 + a² a2/7)
Exc2 2 2 76
a²
V7
3m
7m
(58)
16

---

Firstly, we can observe that, by Cauchy-Schward inequality, Propositions B.2 and B.3, the polar decomposition of Sm, and the fact that T'm & I'm. by definition, we have that
Proof.
KSmm (L - R) ((t)))
((Sm 1_ I. Imm In Sm ((1))())))))
(59)
m Sm 2 (t) L, m
_ (I Lm L R Sm (t)))))) ))z=
(60)
< 1 Sm (1 || 1//2 I. Ti Sm (1 1|2
(61)
m,1
3 | m, (Sm Smm //2 1/22 1/2 Lm (L - Lm LR Sm (12)2
1/2
(62)
$ 3 L mm,
1/2
1/2
L
1/1111| -,
1/2 Lm) 1/2 - 1/ LR
1/2
|
1/2
m.
m7.7
(63)
2
1 | Lm) )
(64)
m.7
Now, we can bound the factors. According Propositions B. and C.4, with probability least
1 - 8, for d E (0, 1] and 2 8
2/) log( 8a²2 ) we have that
log
8x² 78
(1 + a22
2 log See2 a²
|| L 1////11///1111 (I Lm S v2
+
(65)
3m
ym
concluding the proof.
Proof of Theorem 4.5 We are now ready to prove Theorem 4.5. According to Equations (8)
and (10), and similarly to what we did for the posterior mean, we can decompose the error on the variance a single DOF as follows:
lo (t") oposa (t*) = k(t", ,t) - (Sk(t", LR Sk(t", "S' ))2=
-k(t,tt)) (t)) + (Smp((t)) I SmmpSm@@tt"
(66)
<ik(t",t) (t) h(t",t))
+ +((SK(1" ") A.LSSk(( Sk((', (5mp(1") S (())))))
(67)
sia(t"t) k((t"tt))
+ +((Sk(t", ,) -mm L Sk(t", ,))=
( Lm S. T( (()))))) < (((,tt)) - [t(tt,t")]
(68)
Sk(t", ,))=" + (((")"")),_="(((((( -)- Sm " ))) ((8mm(1") (L - Lm,,)) ((" [z)))
+I(Sk(t",)
(69)
Now, we can upper bound the four addends appearing in the decomposition above. The first addend can by directly bounded by Corollary C.2. The second addend of the decomposition in Equation (69) can be bounded by Lemma B.5. The third addend in Equation (69) can be bounded by Lemma B.6. The last addend in Equation (69) can be bounded by Lemma B.7. In this way, retrieve the result of Theorem 4.5, obtaining the following bound holding with probability at least 1 --5. Having defined
17

---

C
2 log
22
2 log
24m
20 log 23
V2
2ng
2og
m
7m
m
v9
3m
7/77
post (1))
post (1*) )1 S (), ()) ( ( ((m)))
+ (Sk((t, - -mmm(()))) L
Sk(t", ;k((t;,))R='
S (( + (( (L - L) (ttn((tt)))) )2=)
+ (() Ln Ln
(70)
SC.
(71)
The final result for the vector-valued GP can be obtained by applying a union bound.
C Concentration Results
We first provide few lemmas for the concentration of the approximate kernel functions that derive from Hoeffding inequality, and then lemma for the concentration of random operators that derives from Bernstein inequality. Again, we denote by A, the operator A + 7/. ||-Il denotes the operator norm, while |-1|2 denotes the Euclidean norm of a vector.
C.1 Approximation of the Kernel Function
Note that if uniform convergence of the RF-HGP posterior is seeked w.r.t. the domain of the func- tion modelled with the HGP, our proofs could be adapted by replacing the following Corollary C.2 with a uniform convergence result. For instance, in the case RFFs, such a result can be found in [39, Theorem and Remark 1].
Lemma C.1. Let 8 = (0,1) 1]. Then, for any (t1, t2), with probability at least - 8, it holds
2a log 2
((t1)) k(tz, b2) s
(72)
Proof. To upper bound the quantity of interest, we can use Hoeffding's inequality for bounded random variables. Let (t1,t2) == bhy (t) (f1)))) (t2) E., (w, , f (w, t2). Since -a² <
(11))))))(2)) S a² according to Assumption 4.2, by Hoeffding inequality, we have that
Pr
Σ
(2,),,
2
t
2e
Tmm
(73)
Therefore, by setting the above upper bound smaller than d, for (0, 1]. we get that with proba- bility at least 1 --6 5
1 m
2a log 2
((t1) k(t,,
Σ
t21,
$
(74)
m
Corollary C.2. Let 8 (0, 1]. Then with probability at least 8, it holds
4 log
21T1
((()))(()) k(tt, (*) S
Vt" TT.
(75)
m
Proof. We apply Lemma C.I each element of T with 8 == 8/T. The claimed result then follows using a union bound.
18

---

Corollary C.3. Let = (0, 1]. Then with probability at least 1 - 5,
2a log 2
[|mm)((1)) (000U)))S(((()))) S
11 ET.
(76)
m
Proof. It holds
S -
Σ
((t)) - k(t,, t2)
(77)
The result thus follows from applying nT times Lemma C.I on the pairs ((t,, t" t' with probability : 6/ and using a union bound.
C.2 Concentration of the Kernel matrix
The following result derives from the Bernstein inequality for sums of random operators separa- ble Hilbert spaces in operator norm.
Proposition C.4 (Proposition 6 and Remark 10 of [20]). Let VLy..., Vp with p N 1, be independent and identically distributed random vectors on separable Hilbert spaces H such that Q = Ev @ v
is trace-class, and for any À >> O there exists constant (A) < oo such that (v, (Q (Q+ S Foo (X) almost everywhere. Let Qp = 7 Ei=1 V, @ V; and take O < X |QQ]- Then for any 8 2 0, the following holds with probability at least 6:
2w + Fol)
(2-
<
(78)
3p
p
where w = log
8T Q 26
Moreover, with the same probability,
Amax
Qx
Q2) -2)))00)
1/2
2w
Fwo
(79)
3p
P
Moreover, for any & E (0, [vll| S a, we have that, with probability at least 1 - 6,
Amax Q (Q - Qp) 1/2
S.
(80)
provided that p 2V
+ Foe (7) log 8a and Q[2].
D Efficient Matrix Inversion and Online Updates
In this section, show how the expression of the posterior mean and variance easily be updated when adding new samples to the dataset.
We recall that the operators S. and Lm,R are matrices and are defined in Appendix A. As discussed in Appendix A, the inversion of Lm,R. involved the posteriors Equations (9) and (10), can be simplified by applying Woodbury identity [36], as follows:
Lm R
=S St +R)
RRR-
(81)
= =R - R-'RRS Sm (1 + S. R.
Rm) - S
(82)
The posterior mean of the HGP in Equation (9) becomes:
Alpost (t") {[R R-R--s - Sm + St R. ( (5m)-'sp R R1 (mt)(t) 1
y)z
(83)
Vn"
= =0(tt)7' [1 St. R. Sm (I + St R Sm))
St. R. -y
(84)
n
- [-- B(I B) - A
(85)
19

---

where A == Sm E R and B == R Sm E Em.m. Moreover, the only term in the expression of the posterior variance of Equation (10)
opost (t') (,((",, + nrr (())(()), (Im R)
(86)
which varies with n is
(Sm.(tt),, (Lm + (m
(87)
= =Ø" (t") S R Sm S R Sm (I S Sm) R Sm x(t"),
(88)
=t=(t") B B(I + B) B] x(t").
(89)
When new human demonstration is gathered, the training is enlarged by adding fnew training points. This means that the matrix Sm is updated by adding Rnew rows (and renormalized), contain- ing the RF embeddings of the training points. The same happens vector y and to the diagonal matrix R, which is enlarged by adding Rnew rows and columns. This means that matrices A and B support online updates. In particular, after initializing A and B to the null matrix, having collected the new embeddings in Sm,new IR "Ree Xm (with normalization n and the new noise variance values in Rnew E "nee XPaee (with normalization neo the updates
follows:
+ A+ Vinaww Se. Rnew new
(90)
B t B S m,neww S Sm,new-
(91)
Having computed the updates, the matrices appearing in the posterior mean and variance can be computed in constant time w.r.t. the current size of the training set during the data streaming.
20


Reviews:
[{"review_id": "t0ML5k2Ga0b", "reviewer": "Reviewer_kgu9", "summary": "", "questions": "", "limitations": "", "rating": 3, "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": "This paper proposes a sparse spectrum Gaussian process model for settings with heteroscedastic noise in learning from demonstration and provides theoretical guarantees for the approximation error with respect to an exact Gaussian process (GP) model. The analysis is based on previous results in the literature on random Fourier features, especially Rudi and Rosasco [25]. An algorithm is also provided to learn the hyper-parameters of the model via expectation maximisation, adapting the method proposed by Kersting et al. [4]. Experimental results with data collected from tasks executed with a robotic arm show the evolution of the error in the approximation and the computation time with respect to the number of Fourier features.", "summary_of_recommendation": "My concerns are mainly due some claims in the theoretical analysis and issues with the experimental comparisons. I'll be looking forward to clarifications during the author discussion phase.", "technical_quality": "Very Good:  Convincing.", "clarity_of_presentation": "Very Good:  Well organized and clear.", "potential_impact": "", "robotics_focus": "Sufficient demonstration on hardware", "quality_of_limitations": "Limitations are addressed clearly", "questions_for_rebuttal": "1. Sec. 4.2, **Feature approximations.** The error bounds in Mutny and Krause [16] deal with possibly continuous function domains and data that might be not i.i.d. The results in Thr. 4.4 and 4.5, however, hold for a finite domain $\\mathcal{T}$, and the size of the domain $|\\mathcal{T}|$ is a component of the error bound. So I am not sure if the claimed improvement in convergence rate is applicable.\n\n2. Sec. 4.2, **convergence rate.** I could not see how a convergence rate of $\\mathcal{O}(n^{3/2}\\log n)$ comes out of the results in Thr. 4.4 and 4.5, when the only terms explicitly dependent on $n$ are $\\mathcal{O}(\\sqrt{\\log n})$.\n\n3. **Experimental comparisons.** There are no experimental comparisons against other sparse GP approximations that can deal with heteroscedastic noise, such as variational approximations based on inducing points. The only baseline is an exact GP. I believe comparisons against another sparse GP baseline should strengthen the paper and provide better insight into the significance of it contribution, especially if there are gains with respect to runtime complexity or approximation errors.\n\n4. **Application to learning from demonstration.** The paper also seems not to apply the proposed RF-HGP model to the generation of motion primitives within robotic control or motion planning tasks. So it is hard to assess the impact of such modelling framework within the robotic tasks it has been proposed for.\n\nMinor issues:\n* When dealing with sparse spectrum GPs, we usually do not resort to eq. 4 and 5 to compute the posterior mean and variance, but to alternative versions obtained by the Woodbury identities, as mentioned in the paper. I believe it would be useful to have such versions explicitly available somewhere in the main paper.\n\n* What are the total number of data points for each task in the experiments? I've only found mention that they are in the order of $10^3$, but not an explicit number.", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The main contribution of the paper seems to be the theoretical results in Sec. 4, which bound the approximation error of the posterior mean and variance of the proposed model with respect to the posterior of an exact GP model. I haven't checked the proofs, but the results seem to be novel and are should be beneficial to applications beyond the robotics domain. However, I have a few concerns regarding a few claims in the paper and the experimental results, which I detail in my questions.", "extra_scores": {"impact": "4: The work contains interesting new ideas that will potentially have major impact in robotics or machine learning."}, "extra_text": {}}, {"review_id": "fG6EKVFJQwz", "reviewer": "Reviewer_qoam", "summary": "", "questions": "", "limitations": "", "rating": 3, "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": "The authors motivate heteroscedastic Gaussian process regression for learning from demonstrations, with random Fourier features to approximate the kernel.", "summary_of_recommendation": "I advocate for acceptance based on the strong theoretical contribution, robot experiments and relevance to CoRL. ", "technical_quality": "Very Good:  Convincing.", "clarity_of_presentation": "Very Good:  Well organized and clear.", "potential_impact": "", "robotics_focus": "Sufficient demonstration on hardware", "quality_of_limitations": "Limitations are addressed clearly", "questions_for_rebuttal": "1. Since you have a Bayesian linear regression model with the random features, could it not be easier to just fit the predictive distribution over trajectories directly rather than the EM approach? E.g.\n\nParametric Gaussian Process Regressors\nMartin Jankowiak, Geoff Pleiss, Jacob Gardner ICML 2020\n\nIt would be rather high-dimensional, but the training would be a lot simpler and more inline with behaviour cloning. \n\n2. See earlier question about function space sampling recovering the temporal correlations of the demonstrations. \n\n3. Why aren't inducing points used as a baseline? Line 285 says they are more complex, but I think the inducing points in the time domain have a nice intuitive interpretation as via-points . \n\n4. Could it be beneficial to use both RFFs and RBF features for the approximation? Since RBFs are good at modelling the datapoints. C.f.\n\nEfficiently sampling functions from Gaussian process posteriors\nJames Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, Marc Deisenroth, ICML 2020", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The use of GPs for regularized learning of demonstrations enjoys a long history, and some rigorous theoretical analysis is welcome. \nWhen reading the paper, it sometimes felt it could be better suited as a ML paper on heteroscedastic GPs, as it's never made clear what makes learning from demonstrations special as a heteroscedastic regression problem. Usually in movement primitive research there is a focus on incorporating different inductive biases, such as periodicity, and blending lower-level primitives through composability.  \n\n**Why are demonstrations heteroscedastic?** The central heteroscedasticity assumption here is motivated by the demonstration data / expert being inconsistent and 'corrupted'. However, while there is some motor noise at play, the inconsistency is often more due to ambiguity. For example, for a 2D go-to task, the expert trajectory could be a straight line, or an arc above or below it. If you fit this data into a heteroscedastic model and sample predictions IID, you will not recover these trajectories samples but rather rough heteroscedastic noise. It could be argued that the demonstration ambiguity is a form of permanent task-induced epistemic uncertainty, rather than aleatoric uncertainty. The benefit of fitting the predictive distribution of a single GP over the whole trajectory sequence (which I discuss in the question section) is that sampling the model in function space should yield trajectories like the demonstrations. Since this model has two GPs for the mean and variance, sampling both in function space and combining should yield the correlated trajectories one would expect? This feature is important (and perhaps unique to LfD) because if these learned primitives are used to initialize an RL finetuning step, you want to sample trajectories that explore 'correctly' by reproducing the correlations of the demonstrations. I think a discussion regarding whether we are modelling aleatoric or epistemic uncertainty, as well as the function-space view of the trajectories data, is a useful and necessary addition.\n\n**Scholarship.** The authors don't cite many robotics papers and omit discussing many similar prior works on the intersection of GPs and movement primitives, for example\n\nProbabilistic movement primitives / Using Probabilistic Movement Primitives in Robotics\nA Paraschos, C Daniel, J Peters NeurIPS 2013 / AuRo\n\nKernelized movement primitives\nY Huang et al IJRR 2019\n\nProbabilistic movement primitives for coordination of multiple human–robot collaborative tasks\nMaeda et al. Autonomous robots 2017\n\nLatent spaces for dynamic movement primitives\nS Bitzer, S Vijayakuma ICRA 2009\n\nLearning from demonstration with model-based Gaussian process\nN Jaquier et al Conference on Robot Learning 2020 \n\nDifferentiable Gaussian Process Motion Planning\nBhardwaj et al ICRA 2020\n\nAdaptation and Robust Learning of Probabilistic Movement Primitives\nGomez-Gonzalez et al. TRO 2020\n\nBayesian Disturbance Injection: Robust Imitation Learning of Flexible Policies\nOh et al, ICRA 2021\n\nProDMP: A Unified Perspective on Dynamic and Probabilistic Movement Primitives\nLi et al RAL IEEE 2023\n\nThere was also a paper at CoRL last year using GPs as action priors for MPC and looked at random Fourier features\n\nInferring Smooth Control: Monte Carlo Posterior Policy Iteration with Gaussian Processes\nJ Watson, J Peters - Conference on Robot Learning, 2023\n\nWhile probabilistic movement primitives are not sold as GPs, if you look at the math they are essentially Bayesian linear regression models with RBF features, so almost like a sparse inducing point GP approximation.  The difference is that the weight distribution is fit in weight-space, by fitting a primitives per demonstration first, rather than prediction space over all demonstrations. This is in order to not make the IID assumption.\n\n**Quadrature Random features.** Lines 52 - 55 seems to separate QRFFs from RFFs, but QRFFs are using Gauss Hermite quadrature to approximate the Gaussian integral that RFFs use Monte Carlo samples for. Whats nice about QRFFs is that when the input is 1 dimensional (like time in this paper) the curse of dimensionality associated with Gauss Hermite quadrature is not so bad and you can achieve a very good approximation of the integral. The authors are correct that this approximation truncates the spectrum, but this is true for RFFs too in practice based on whatever the largest frequency sample is. \n ", "extra_scores": {"impact": "3: The work is incremental and unlikely to have much impact, even though it may be technically correct and well executed."}, "extra_text": {}}, {"review_id": "5wMhgC-XHap", "reviewer": "Reviewer_zk8k", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 1, "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": "This paper learns a vector-values heteroscedastic GP using Random Fourier features (RFF) for modeling motion primitives. The authors additionally provide bounds on the approximation error of the posterior mean and variance of the HGP approximated through  RFF .", "summary_of_recommendation": "I did not find the paper clear and cohesive. The motivation of the paper was especially not clear to me. \n\n-----------POST REBUTTAL ------------\n\nBased on the rebuttal provided by the authors, I am happy to raise my score to accept", "technical_quality": "Good: Generally convincing, but some additional experiments or analyses needed.", "clarity_of_presentation": "Fair: Major revision needed to improve organization or clarity.", "potential_impact": "", "robotics_focus": "Sufficient demonstration on hardware", "quality_of_limitations": "Limitations are addressed clearly", "questions_for_rebuttal": "1. What is 'd'? If 'd' is the output dimensionality, it seems that d>1 in the experiments, and the strict improvement in the error bounds does not hold.", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "I didn’t find the paper to very well organized even though the individual sections were clearer to read. The was not a clear section on Related Works (some subsections of Section 4 seemed to refer to related approaches, but works related to ). Furthermore, the paper covers several different concepts – learning a motion primitive, learning a heteroscedastic GP with RFF approximation and bounds on this approximation -  that are not presented cohesively. The presented approximation bounds for RFF seemed significant. But it seems to be a strict improvement over previously presented bounds under some assumptions that seem to be too conservative.\n\n\n**Strengths:**\n- theoretical results\n- results on hardware\n\n**Weaknesses:**\n- Paper didn’t seem to be organized very well\n- Some clarity issues: for example it is not clear to me what ’d’ is. Is it the dimensionality of the outputs?\n- There isn't much reporting on related works. This work is not compared to other approaches such as probabilistic movement primitives and other dynamic movement primitives that can also be used for LfD. \n", "extra_scores": {"impact": "2: Unsure, difficult to assess."}, "extra_text": {}}, {"review_id": "4cUyp9NEWd", "reviewer": "Reviewer_JwFW", "summary": "", "questions": "", "limitations": "", "rating": 2, "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": "This paper provides bounds for the approximation error for the mean and variance of heteroscedastic Gaussian processes (HGs) whose kernel is approximated via random features (RFs). The interest of such approximation is to reduce the computational complexity of the posterior computation by reducing the cost of the training kernel matrix inversion. The resulting HGP are evaluated in the context of variable impedance learning from demonstrations.", "summary_of_recommendation": "Although the provided bounds seems theoretically sound, I believe that the paper generally lacks significance and novelty. Moreover, it does not consider the latest state of the art in terms of efficient GPs and its experimental evaluation is limited. Therefore, I recommend rejecting the paper.\n\n----\nPost-rebuttal:\n\nI updated my score from strong to weak reject in appreciation for the theoretical bound presented in Section 4. However, I still advocate for rejection as the rest of the paper generally lacks novelty, and the impact of Section 4 in (robotics) applications beyond VIL remains unclear.", "technical_quality": "Good: Generally convincing, but some additional experiments or analyses needed.", "clarity_of_presentation": "Very Good:  Well organized and clear.", "potential_impact": "", "robotics_focus": "Sufficient demonstration on hardware", "quality_of_limitations": "Additional details required", "questions_for_rebuttal": "- How does the considered RF approximation related to recent efficient variational and/or feature-based approximations in terms of computational cost and approximation error? \n- Are Theorems 4.4 and 4.5 limited to input domain in R or are they valid for any input domain? ", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The main strength of the paper is the theoretical analysis leading to the bounds for the approximation error of HGP with RF approximation. The main interest of these bounds is that they apply to the heteroscedastic case. In addition, the paper is clearly written.\n\nDespite the novel bound, I am generally concerned with the originality and significance of the paper, as well as with the application scope in which the contributions are phrased. \n\nFirst, HGPs have been widely studied and already applied to robot learning, see e.g., [4, 5]. Second, a large body of work is concerned with reducing the computational cost of (H)GPs with respect to the number of training samples, including, as discussed in the paper, variational methods based on inducing points and feature-based kernel approximations. Therefore, the RF approximation considered in the paper is not new and has been widely studied in the literature. In addition, recent works such as [Hensman, 2018], [Dutordoir, 2020], [Cunningham, 2023] (to cite a few) proposed efficient variational and/or feature-based approximations beyond RF approximations allowing the consideration of millions of datapoints with some order of magnitude faster than traditional (sparse) GPs. In other words, there exist other, more performant approximations that considerably reduce the computational cost of GPs and that could be considered in robot learning as well. \n\nConcerning the scope, the paper seems to limit the input space of the considered GPs to time (i.e., R), while GP are generally interesting because they are able to deal with high-dimensional input spaces. I am not sure if the main theorems 4.4 and 4.5 depend on this assumption. If not, it would be good to precise it. Moreover, the paper phrases its main contributions in the context of variable impedance learning (VIL). It seems to me that HGPs with RF approximations are applicable to many other robotic problems, including, e.g., learning from demonstrations of general trajectories, Bayesian optimization, etc. Therefore, reducing the scope of the paper to VIL generally undermines the contributions of the paper. \n\nFinally, the experimental evaluation validates the fact that RF approximations are computationally cheaper than computing the HGP without approximation and that they provide a good approximation of the mean and variance. However, I would have liked to see a clear link between the empirical evaluation and the bounds provided in Section 4. For instance, the bound may be represented in Figures 3 and 6. Moreover, the performance of the RF approximation is not evaluated against other efficient variational and/or feature-based approximations. This would provide the reader with evident to choose one of these approximation when using GPs or HGPs in the context of robot learning from demonstrations. In particular, a comparison against classical variational GPs would be straightforward. Moreover, as recent approaches have been shown to be very efficient to approximate GPs with many training points, they should be explicitly considered in the paper. \n\n\n[Hensman, 2018] Hensman, Durrande, & Solin, \"Variational Fourier Features for Gaussian Processes\", JMLR, 2018.\n\n[Dutordoir, 2020] Dutordoir, Durrande, & Hensman, \"Sparse Gaussian Processes with Spherical Harmonic Features\", ICML, 2020.\n\n[Cunningham, 2023] Cunningham, Augusto de Souza, Takao, van der Wilk, & Deisenroth, \"Actually Sparse Variational Gaussian Processes\", AISTATS, 2023.\n\n\n--------\nPost-rebuttal comments:\n\nAlthough I appreciate the efforts of the authors to answer the reviewers during the rebuttal, my general concerns about the paper mostly remain. \n\nThe main contribution of the paper is the presented theoretical analysis leading to the bounds for the approximation error of HGP with RF approximation. However, the contributions of the paper remain very limited from a robotics perspective as (i) HGPs have been widely studied and already applied to robot learning, and (ii) the paper shows only a single case of application (using the tasks described in [10]) where the derived bound may be useful. Despite that I appreciate the presented theoretical derivation, I find that its impact on (robotics) applications is not well reflected by the paper. In this sense, I agree with Reviewer qoam that the theoretical contributions of the paper may be better suited for publication within the ML / GP community due to its potential applications beyond the robotics domain.\nIf there is interest to publish the paper for robotics, I would expect to see a clear discussion on the impact of its results for robotics applications. \n\nLastly, I appreciate the discussion provided during the rebuttal about RF approximations and variational approaches. However, I believe that the newly-added comparison lacks significance, and I think that the discussion (and eventual comparisons) should consider more recent works and implementations.", "extra_scores": {"impact": "2: Unsure, difficult to assess."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper combines Gaussian processes with random Fourier features to derive a variable observational noise model (aka heteroskedastic) used for learning from demonstration. Theoretical guarantees for the approximation are provided as well as an algorithm to learn the hyperparameters of the model. Experiments with a robotic arm demonstrate the benefits of the approach on assembly and bed-making tasks.

Reviewers were positive about the theoretical contributions of the work but were concerned about the presentation and the overall connection and relevance of these results to robotics. In the rebuttal and in the new version, authors have addressed most of these issues but the significance to robotics remains unclear. All in all, I believe the theoretical contributions are important and can be applied to other problems authors have yet to consider.
Decision: Accept (Poster)