Robust Gaussian Processes via Relevance Pursuit
Authors: Sebastian Ament
Elizabeth Santorella
David Eriksson
Benjamin Letham
Maximilian Balandat
Eytan Bakshy
Venue: NeurIPS 2024 poster (2024)
Keywords: Gaussian process
robust regression
Bayesian optimization
submodular
Primary area: probabilistic_methods
TL;DR: Robust Gaussian Processes with strong empirical performance, and surprising convexity and approximation guarantees.

Abstract: Gaussian processes (GPs) are non-parametric probabilistic regression models that are popular due to their flexibility, data efficiency, and well-calibrated uncertainty estimates. However, standard GP models assume homoskedastic Gaussian noise, while many real-world applications are subject to non-Gaussian corruptions. Variants of GPs that are more robust to alternative noise models have been proposed, and entail significant trade-offs between accuracy and robustness, and between computational requirements and theoretical guarantees. In this work, we propose and study a GP model that achieves robustness against sparse outliers by inferring data-point-specific noise levels with a sequential selection procedure maximizing the log marginal likelihood that we refer to as relevance pursuit. We show, surprisingly, that the model can be parameterized such that the associated log marginal likelihood is strongly concave in the data-point-specific noise variances, a property rarely found in either robust regression objectives or GP marginal likelihoods. This in turn implies the weak submodularity of the corresponding subset selection problem, and thereby proves approximation guarantees for the proposed algorithm. We compare the model’s performance relative to other approaches on diverse regression and Bayesian optimization tasks, including the challenging but common setting of sparse corruptions of the labels within or close to the function range.

Robust Gaussian Processes via Relevance Pursuit
Sebastian Ament
Elizabeth Santorella Meta santorella@meta, com
David Eriksson Meta
Meta ament@meta. com
deriksson@meta.com
Ben Letham
Maximilian Balandat Meta balladattmeta..om com
Eytan Bakshy Meta
Meta bletham@meta.com com
ebakshy@meta.com
Abstract
Gaussian processes (GPs) are non-parametric probabilistic regression models that are popular due to their flexibility, data efficiency, and well-calibrated uncertainty estimates. However, standard GP models assume homoskedastic Gaussian noise, while many real-world applications are subject to non- Gaussaa corruptions. Vari- ants of GPs that are more robust to alternative noise models have been proposed, and entail significant trade-offs between accuracy and robustness, and between computational requirements and theoretical guarantees. In this work, we propose and study a GP model that achieves robustness against sparse outliers by inferring spieseii noise levels with a sequential selection procedure maximizing the log marginal likelihood that we refer to as relevance pursuit. We show, surpris- ingly, that the model can be parameterized such that the associated log marginal likelihood is strongly concave in the data-point-specifie noise variances, property rarely found in either robust regression objectives or GP marginal likelihoods. This in turn implies the weak submodularity of the corresponding subset selection problem, and thereby proves approximation guarantees for the proposed algorithm. We compare the model's performance relative to other approaches on diverse gression and Bayesian optimization tasks, including the challenging but common setting of sparse corruptions of the labels within or close to the function range.
Introduction
Probabilistic models have long been central part of machine learning, and Gaussian process (GP) models are a key workhorse for many important tasks [54), especially in the small-data regime. GPs flexible, non-parametric predictive models known for their high data efficiency and well-calibrated uncertainty estimates, making them a popular choice for regression, uncertainty quantification, and
downstream applications such as Bayesian optimization (BO) g 24 26
and active learning 6 55
GPs flexibly model distribution over functions, but assume a particular observation model. The standard formulation assumes i.i.d Gaussian observation noise, i.e., y(x) f(x) + where f(x) is the true (latent) function value at point x and c ~ (0, o²), implying a homoskedastic Gaussian likelihood. While mathematically convenient, this assumption can be limitation in practice, since noise distributions are often heavy-tailed or observations may be corrupted due to issues such as sensor failures, data processing eπrоrs, software bugs. Using standard GP model in such settings can result in poor predictive performance.
A number of robust GP modeling approaches have been proposed to remedy this shortcoming, most of which fall into the following broad categories: data pre-processing (c.g., Winsorizing), modified likelihood functions (e.g.. Student-t), and model-based data selection and down-weighting procedures.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

These approaches offer different trade-off's between model accuracy, degree of robustness, broad applicability, computational requirements. and theoretical guarantees.
In this paper, we propose a simple yet effective implicit weightiigg approach that endows GPs with high degree of robustness to challenging label corruptions. Our approach is flexible and can be used with arbitrary kernels, is efficient to compute, and yields provable approximation guarantees. Our main contributions are as follows:
1. We propose modification to the standard GP model that introduces learnable data- lata-point-specifie noise variances.
2. We introduce novel greedy sequential selection procedure for maximizing the model's marginal log-likelihood (MLL) that we refer to as relevance pursuit.
3. We prove that, under particular parameterization, the MLL strongly concave in the data- point-specific noise variances, and derive approximation guarantees for our algorithm.
4. We demonstrate that our approach, Robust Gaussian Processes via Relevance Pursuit (RRP), performs favorably compared to alternative methods across various benchmarks, including challenging settings of sparse label corruptions within the function's ''s range, see e.g. Figure
2 Preliminaries
Posterior Predictive Distributions
We aim to model a function f X + R over some domain X Rd. With standard Gaus- sian noise model, for x; E X we obtain observa- tions y/- = fx,) + Eg, where & ~ NN0, o²) are i..d. draws from a Gaussian random variable. |I-II denotes the Euclidean norm unless indicated otherwise.
Gaussian
Truth
Student-t
Compttonn
1.5
Relevance Pursuit
Observations
1.0
0.5
0.0
2.1 Gaussian Processes
-05
GP f GP(p(-), k)( is fully defined by its mean function X + I and covariance or kernel function KØ X x X + E. which is parameterized by O. Without loss of general- ity, we will assume that μ = 0. Suppose we have collected data D = ((x,, (1) where == {vi )= Let Me E S
-1.0
-15
0.0
0.8
10
Figure 1: Comparison of RRP standard GP and a variational GP with a Student-t likelihood on regression example. While the other models are led astray by the corrupted observations, RRP successfully identifies the corruptions (red) and thus achieves a much better fit to the ground truth.
X : { denote the covariance matrix of the data set, i.e. [Eali = (x,,, xj) + 81j0², where d;j is the Kronecker delta. The negative marginal log- likelihood (NMLL) is given by
-2C(0) = lologgyy y E + log det Eg + n log 2n.
In the following, we will suppress the explicit dependence of the kernel matrix on for brevity of notation. For comprehensive background on GPs, we refer to Rasmussen et al. 54
2.2 Noise Models
Additive, heavy-tailed noise Instead of assuming the noise term E; in the observation model be Gaussian, other noise models consider zero- mean perturbations drawn from distributions with heavier tails, such as the Student-t [32], Laplace (39] or atSble distributions. These types of errors are common in applications such as finance, geophysics, and epidemiology [18). Robust regression models utilizing Student-t errors are commonly used combat heavy-tailed noise and outliers.
Sparse corruptions In practice, often a small number of labels are corrupted. We will refer to these as "outliers," " though emphasize that the corrupted values may fall within the range of normal outputs. Sparse corruptions are captured by model of the form y Zxf(x, + (1 Z,)W,, where

---

Z; E {0,1) and W E R is random variable. Note that W; need not have (and rarely has) fx;) as its mean. For instance, consider a faulty sensor that with some probability p reports random value within the sensor range [yt, yb]. In this case Z ~ Ber(p) and W; ~ U[y,, y/h. Software bugs, such as those found in MI. training procedures, or errors in logging data can result in sparse corruptions.
3 Related Work
Data pre-processing Data pre-procsssing can be an effective technique for handling simple forms of data corruption, such values that fall outside valid range of outputs. With such pre-processing outliers are handled upstream of the regression model. Common techniques include the power transformations [15], trimming, and winsorization. These methods add substantial bias if not used carefully, and generally do not handle data corruptions that occur within the normal range of the process to be modeled. See [17] for a review on data cleaning.
Heavy-tailed likelihoods One class of robust methods uses additive heavy-tailed noise likelihoods for GPs. particularly Student-/ [32]. Laplace (53). and Huber (. and could be extended with a-Stable distributions, which follow a generalized central limit theorem [5]. These models are less sensitive to outliers, but they lose efficiency when the outliers
sparse subset of the observations, opposed to global heavy-tailed noise. Furthermore, model inference is no longer analytic, necessitating the use of approximate inference approaches such MCMC [49]. Laplace approximation [64], expectation propagation (EP) 32). Expectation Maximization [53] or variational inference [61]. Shah et al. [56] take related approach using Student-t process prior in the place of the GP prior. Unfortunately, the Student-t process is not closed under addition and lacks the tractability that makes GPs so versatile. Alternative noise specifications include hierarchical mixture of Gaussians [19 and a "twinned" GP
model 488 that uses a two-component noise model to allow outlier behavior to depend on the inputs. This method is suited for settings where outliers are not totally stochastic, but generally is not able to differentiate "inliers" from outliers when they can occur with similar inputs.
Outlier classification Awasthi et al. 11 introduces the Trimmed MLE approach, which identifies the subset of data points (of pre-specified size) under which the marginal likelihood is maximized. Andrade and Takeda [10] fit GPs using the trimmed MLE by applying projected gradient method to an approximation of the marginal likelihood. The associated theory only guarantees convergence to a stationary point, with guarantee
quality. When no outliers present, this method can be worse than a standard GP. Li et al. [40] propose heuristic iterative procedure of removing those data points with the largest residuals after fitting standard GP, with subsequent reweighting. The method shows favorable empirical performance but has no theoretical guarantees, and fails the largest residual is not associated with an outlier. Park et al. [52] consider a model of the form yı 6; + fxx)) + 6i where outliers are regarded as data with a large bias 6;- Their random bias model is related to our model in that it also introduces learnable, tata-point-specffe variances. However, inference is done in one step by optimizing the NMLL with an inverse-gamma prior on the 5, ks,, which - in contrast to the method proposed herein - generally does not lead exactly sparse 0:'s.
Sample re-weighting Altamirano al. (2 propose robust and conjugate GPs (RCGP) based on a modification to the Gaussian likelihood function that is equivalent to standard GP inference, where the covariance of the noise o²I is replaced by o² diag(w 2) and the prior mean m is replaced by mw = + o² V, (w"). The authors advocate for the use of the inverse multi-quadratic weight function wxx, y) = B(1 + (y m(x))² /c2))//2 which introduces two additional :ppee-parameter::: the soft threshold c, and the "learning rate" 3. Importantly, the weights w are defined a-priori as a function of the prior mean mxx) and the targets y, thereby necessitating the weights identify the correct outliers without access to a model. This is generally only realistic if the outlier data points are clearly separated in the input or output spaces rather than randomly interspersed.
4 Robust Gaussian Process Regression via Relevance Pursuit
Our method adaptively identifies a sparse set of outlying data points that are corrupted by a mechanism that is not captured by the other components of the model. This is in contrast to many other approaches to robust regression that non-adaptively apply a heavy-tailed likelihood to all observations, which can be suboptimal if many observations are of high quality.
3

---

Relevance Pursuit Model Trace
Model Evidence
0.08
Prior
100
Posteroo
Postettor Marginal L.kellhund
0.06
0.5
0.0
0.04
-0.5
0.02
18
[5]
[5|=8
[S]
-1.0
18
[S] [S]
0.0 0.2
0.8
20
Supgoct Size (Price m Expomessuall/ncannnoo
Figure 2: Left: Evolution of model posterior during Relevance Pursuit, as the number of data-point specific variances |S|] increases (from light colors dark). Red points indicate corruptions that were generated by uniformly sampling from the function's 's range. Right: Comparison of posterior marginal likelihoods as a function of a model's [S]. The maximizer boxed in black is the preferred model.
4.1 The Extended Likelihood Model
We extend the standard GP observation noise variance a² with data-point-specific noise variances {= so that the i-th data point is distributed as
yi X ~ N (x(x)),
This is similar to Sparse Bayesian Learning 660 in which weight-specific prior variances control feature's 's degree of influence on a model's predictions. The marginal likelihood optimization of p in (2) gives rise to an automatic mechanism for the detection and weighting of outliers. The effect of on the estimate of f vanishes P; + similar to the effect of the latent varibales h in Bodin et [1]'s extended GP model fxx h), though h requires MCMC for inference. While many heteroskedastic GP likelihoods model noise an input-dependent process [28 351, our formulation does not require such assumptions, and thus suitable for corruptions that are not spatially correlated.
An elegant consequence of our modeling assumption is that we can compute individual marginal- likelihood maximizing P. in closed form when keeping all Pj for o i fixed. In particular,
Wj) j * i), P = pvr + pee where P. Pi E 2q, (pyli O, and e; is the ith canonical basis vector. Then keeping P fixed.
Lemma 1. /Optimal Robust Variances] Lei D.a = {(xj,,;)
p = arg max C(p,, + p.e.) = [(y: - E[y(x))) -ly(x))
(3)
-- where y(x,) = =5 [xx,) +e.. These quantities can expressed as functions of 12-1- = (K+Do2+p)=
- | and V(y(x;) Dr = 1/ |
where +p is diagonal matrix whose entries are o² + p.
The first component the i-th data point, i.e., the leave-one-out (LOO) cross-validation error (54). The second component (x,(x,) + (uy is the LOO predictive variance. The optimal solution to ps is only non-zero for those observations whose squared LOO error is larger than the LOO predictive variance at that point.
+e: of 3 is the empirical error to y of the model trained without
4.2 Optimization with a Maximum Number of Outliers
Without additional structure, inference of the noise variances P: does not yield desirable models, as the marginal likelihood can be improved by increasing the prior variance P: of any data point where Eq. (3) is greater than zero, even if that is due to regular (non-outlier) measurement noise. To avoid this, we constrain the number of non-zero Ps, that is, ||p|| = [[0 < p.) < < n. While this sparsity constraint mitigates over-flexibility, it gives rise to formidably challenging optimization problem, as there are combinatorial number of sparse outlier sets to consider. Even if the number of outliers n. were known, exhaustive search would still require considering n--hoose-n,, possibilities.

---

For tractability, we iteratively add data points to a set of potential "outliers" by setting their associated P; to be nonzero, using the closed-fome form expression for the optimal individual p; variances in Lemmal [ As the algorithm seeks to identify the most "relevant" data points (as measured by L) upon completion, we refer to it as Relevance Pursuit. This is Algorithm|I [ with useBayesian
as false.
Specifically, this is the "forward" variant; Algorithm| in the Appendix presents an alternative "backward" variant that we found to work well the number of corrupted data points is large.
Crucial to the performance of the optimizer, it never removes data from consideration completely; a data point is only down-wen wown-weighted/ if is apparently an outlier. This allows the down-weighte to be reversed if data point appears "inlying" after having down-weijs lown-weighted other data points, improving the method's robustness and performance. This is in contrast to Andrade and Takeda [[O]'s 's greedy algorithm, in which the exclusion of a data point can both increase or decrease the associated marginal likelihood. This means that their objective not monotonic, a necessary condition to provide constant-factor submodular approximation guarantees for greedy algorithms, see Section 5
Algorithm 1 Relevance Pursuit (Forward Algorithm)
Require: X, y, schedule K = (k), h2 .2....... .......,,),. kix),
useBayesianMooillssectiioon (boolean)
Initialize So S (1...., n (typically So = ()
for i in (I...., [K]] do
Optimize MLL: Ps. < arg maxps, L (ps.). where Ps. {pp Py = py:
Vj E S.}.
Expand Support:
Compute Ai( +- maxp; L(ps, + p,e,) ((ps.) for all j $ S via Lemma
A; < : Ju, such that (,) 2 (j'")) for all E A, and j' & (A, U U.).
S,+1 t S,UA
if use SoweessinMddelSeeeeettee Compute the marginal likelihood p(D ~ pDDSS, Ps.)
then
S* + arg maxs, ((S).))((S.)
else
Sk. Return S", Pg.
4.3 Automatic Outlier Detection via Bayesian Model Selection
In practice, it is often impossible to set hard threshold on the number of outliers for particular data set. For example, a sensor might have a known failure rate, but how many outliers it produces will depend on the specific application of the sensor. Thus, is often more natural to specify prior distribution p(5) over the number of outliers, rather than fix the number a priori. We leverage the Bayesian model selection framework [66] [5) determine the most probable number of outliers in data- and model-dependent way, aiming to maximize ((S) This gives rise to Algorithm I with Bawaaeeeeeeeeeeeeeeeeeeaaeeeeeee
true.
Computationally, we start by iteratively adding outliers up to the maximal support of the prior, similar to the procedure described in Section 4.2] We store trace of models generated at each iteration, then approximate the model posterior (() x p(D) (S )p(S)) at each point in the trace. As the exact f p p(D)S,, Ps, ddps p(D)S,, Pss)
posterior is intractable, we approximate it with p(D)S,)
Finally, we select the model from the model trace {8.); that attains the highest model posterior likelihood. Imposing a prior on the number of outliers differs notably from most sparsity-inducing priors, which are instead defined on the parameter values, like norm regularization. In practice, p(S) can be informed by empirical distributions of outliers. For our experiments, we use an exponential prior on [S] to encourage the selection of models that fit as much of the data tightly as possible.
Regarding the schedule IC in Algorithm I the most natural choice simply to add one data point at a time, i.e. K = 1,..), but this can be slow for large n. In practice, we recommend schedules that test fixed set of outlier fractions, e.g. K. = (0.05n, 0.05n....).
5 Theoretical Analysis
We now provide theoretical analysis of our approach. We first propose re-parameterization of the P. that maps the optimization problem to compact domain. Surprisingly, the re-parameterized
5

---

problem exhibits strong convexity and smoothness when the base covariance matrix (excluding the p;) is well-conditioned. We connect the convexity and smoothness with existing results that yield approximation guarantees for sequential greedy algorithms, implying a constant-factor approximation guarantee to the optimal achievable NMLL value for generalized orthogonal matching pursuit (OMP), greedy algorithm that is closely related to Algorithm
5.1 Preliminaries for Sparse Optimization
The optimization of linear models with respect to least-squares objectives in the presence of sparsity constraints has been richly studied in statistics [59]. compressed sensing f 62], and machine learning 2 [67. Of central importance the theoretical study of this problem class the eigenvalues of sub-matrices of the feature matrix, corresponding to sparse feature selections and so often referred to as sparse eigenvalues. The restricted isometry property (RIP) formalizes this.
Definition 2 (Restricted Isometry Property). An (n x mm-matrixx A satisfies the r-restricted isometry property (RIP) with constant 8, E (0, 1) for every submatrix As with |S] S m columns,
(1 6,) | < | || S (1 + 6-]||
where xs E IR", This is equivalent to 6, I (AA 1 (1 -,)L.
The RIP has been proven to lead to exact recovery guarantees [16] as well as approximation guarantees (2)_ Elenberg al. [23 generalized the RIP to non-linear models and other data likelihoods, using the notion of restricted strong convexity (RSC) and restricted smoothness.
Definition 3 (Restricted Strong Convexity and Smoothness). A function f R.d + R is m,.restricted strong convex and My-restricted smooth for all (x, x') in the domain D, C (IR x )R
m - 12/2 S f(x')
v[f](x)
(x' S M, -x -x|2/2 12/
In the context of sparse optimization, we let D, be the set of tuples of r-sparse vectors whose difference is also at most r-spars.. In particular, D, = {(x,x' x') s.f. ||x||o, ||x/|lo, ||x - xllo r).
44 is a greedy algorithm that keeps track of a support set S of non-zero coefficients, and expands the support based on the largest gradient magni- tudes, applied to the marginal liklihood optimization problem, S S U arg maxjes |V [p((p)] Algorithm I Ieeneralizes OMP (62 by allowing more general support expansion schedules K. and specializes the support expansion criterion using the special problem structure exposed by Lemmal
Generalized orthogonal matching pursuit (OMP) 4
5.2 The Convex Parameterization
The NMLL L of a GP TT is the sum of a convex function (-) and concave function log (.) of K, and is therefore not generally convex as function of the hyper-parameters 0, including the robust variances p. Here, we propose a re- re-parameterization that allows us to prove strong convexity guarantees of the associated NMLL. In particular, we let p(s) diag (Ko) C ((1 (1-s)-- - 1), where
K, := k(X, X) + I'1 and the inverse is element-wise. Note that ((s) is a diffeomorphism that maps from the compact domain s [0, 1|" to the entire range of p [0, oo]"
Henceforth, we refer to the original p as the canonical or p-parameterization and the newly proposed p(s) as the convex or s-parameterization. Lemma shows the Hessian of the s- s-parameterization.
Lemma 4. /Reparameterized Hessian] Let K2 k(X,X) + I + Dp(s)- Kg
diag (.) 1,, diag(K, 112 and à = -
Then
H [ = D [2 aaà o (K (K-1-1)) r)) diag(( (K -- - の o
D
Based on this representation, we now derive conditions on the eigenvalues of K that imply the m-strong convexity and M-smoothness of the NMLL. Lemma 5. [Strong Convexity via Eigenvalue Condition] Let K2 as Lemma Then Hg > m if
(2)_1x x-2 m)
Amin m 2(1 Amin// Ammax)
|y
(4)
where Aminm.ax (resp. )min.max) are the smallest and largest eigenvalues of Kg. respectively Kg.

---

12
Marginal 1 Likelihood as Function of Rbo
100
Second Derivative of the Marginal Likellhood
3.0
Hevsie Condiion Namber
10
Value
Optinsal
Scooed Denivatvve (west Rho)
M-000
B 0.53
Margical Hithhbood
Gradien
Conves Region
0.25
b 0.875
104
Hessian
Concove Region
50
075
b 1.0
2.0
bl 0.5
Optinal
2.5
bl 0.625
Non-comee
1.5
2.5
0.5
-50
Marginal Likelihood as Function of Rholo)
Second Derivative of the Marginal Likelihood
Hoasio Condition Number
12
10
3.0
Marginal Litebbood
10
Gradient (w.x1.))
2.5
10
Hessian (W.L5.5
Optimal
Desivalive
Bhol))2
2,0
Convex Region
1.5
1.0
0.5
0.0
1.0 0.0
Rbo(s]
Figure 3: Top: The behavior of the - log C(p) with respect to the canonical parameterization of p. Bottom: The behavior of log L(p(s)), highlighting the convexity property. Left: The value, and first two derivatives of - log L for a 1d example. Center: The second derivatives of 1d - log L as a function of |v/. The s-parameterization is everwhere convex for all considered yyl. while the canonical o-parameterization only convex around the origin and only for [y| > 0.5. Right: The heatmaps highlight that the original parameterization is non-convex (red) for larger values of p, and quickly becomes ill-conditioned, whereas the parameterization p(s) is convex and much better conditioned.
The behavior Lemma predicts is surprising and validated in Fig. Notably, the denominator "blows up" as K becomes close to unitary, making the inequality more likely to be satisfied, an indication that the convexity property of the NMLL intimately linked the RIP (Def, [ Note that Lemma 5 is a condition for par---sppttorstriiiee the approximation guarantees that rely on restricted convexity (Def. 3 However, sparse eigenvalues are generally difficult to compute exactly. Fortunately, covariance matrices of GPs naturally tend to exhibit a property that facilitates a different sufficient condition for convexity for all s [0, 1)".
convexity, which is stronger than is necessary for
Definition 6 (Diagonal Dominance). A matrix A is said to be 8-diagonally dominant if the elements asş satisfy Eity [ats] 8/au for all i.
Intuitively, the (s(s) that selected to be non-zero by the greedy algorithm take on large values, further encouraging the diagonal dominance of the sub-matrix K associated with the support of p. For this reason, the following condition on K, is sufficient guarantee convexity for all E [0, 1)a Lemma /Strong Convexity via Diagonal Dominance] Let m > 0 and Ko 8-diagonally dominant
with ((5 - m.) 255 9m + 17) /4 S (5 /17 /
Rs 0.44 and
Amin (Ko) 5)2 5)2
2) + -0)-1-(1 5) m 2 1lr12
2( 2(1
8)/(1 8))
Then the NMLL is m-strongly convex for all 10, -."",... i.e. p(s) E [0]"
We attain similar results for M M-smoothness, see Lemma 13 and Lemma 14 in the Appendix. Having proven m-noveexity and M M-smoothness conditions, we appeal to the results of Elenberg et [2]
Theorem 8. [Approximation Guarantee] Let Kn k(X, X) + I²I be 5-diagonally dominani,
Smax > O be an upper bound ||s|| sor and suppose lyll, satisfy the boundy of Lemmas 7 and 14 guaranteeing m-convexity and M -ssmoothness of the NMLLfor some m > O, M > 1/(1 Smax) Let SOMP(T (r) be the r-sparse vector attained by OMP on the NMLL objective for r steps, and let SOPT (r) = arg max||>]l =r.l=l= Stmas L(p(s)) be the optimal r-sparse vector. Then any 2r S n,
È (p(somp(r)) 2 (1 - -- mm/aa) E (p(SOPT(r)))
where 2(.)= L(-)) L()) is normalized that maxes [(ss) 2 for any support S.
7

---

A limitation of the theory is that it assumes the other hyper parameters of the GP model to be constant, as doing otherwise would introduce the non-convexity that common most marginal likelihood optimization problems. In practice, we typically optimize jointly with the other hyper-parameters of the model each iteration of RRP, as this yields improved performance, see App. D.5 for details.
6 Empirical Results
We evaluate the empirical performance of RRP against various baselines number of regression and Bayesian Optimization problems. Specifically, we compare against a standard GP with a Matern- 5/2 kernel ("Standard GP"), data pre-processing through Ax's adaptive winsorization procedure ("Adapt. Wins.") [12]. and power transformation ("Power Transf.") [15]. Further, we also consider Student-t likelihood model from Jylänki al. [32] ("Student-t"), the trimmed marginal likelihood model from Andrade and Takeda [10] ("Trimmed MLL"), and the RCGP model from Altamirano al. 2] Unless stated otherwise, all models are implemented in GPyTorch (25 and experiments in this section use 32 replications. See Appendix D for additional details.
6.1 Regression Problems
Synthetic We first consider the popular Friedman and Hartmann6 [21] test functions from the literature. We use two data generating processes: uniform noise, extreme outliers at some fixed value, and heavy-tailed (Student-/) noise function values. In these experiments, compare the performance predictive log-iikeldoodd The results are shown in Fig.
Hartmann6, 5% Constant Corruptions
Hartmann6, 20% Constant Comuptions
Hartmann6, Student-t Corruptions
0
Predictive Log Likelibood
44
Friedmam10, 10, 10% Uniform Corruptions
Friedman 10, 25% Uniform Coruptions
Friedman10. 10, Student-c Corruptions
Sudeni-1
0
RRF
Stndent
Fower Tranff
Adopt Wim.
Trimned MLLL
RRF
Studant-
Power Traisf
Standed GP
Adopt Wim.
Triounad MLL
0.2
0.5
1.0
Standad
Ourlier Fraction
Distribution predictive test-set log likelihood for various methods. Methods ommitted those that performed substantially worse. Right: Predictive log likelihood function of the corruption probability for GP model with the Student-t likelihood only starts outperforming RRP the corruption probability increases beyond 40%, and exhibits large variance in outcomes, which shrinks as the proportion of corruptions increases. All methods not shown were inferior to either RRP or Student-t.
Figure 4: Left:
corrpptions with two degrees of freedom. The
Twitter Flash Crash In Fig. we report comparison to Altamirano et al. [2]s RCGP on data from the Dow Jones Industrial Average (DJIA) index on April 22-23 2013, which includes sharp drop at 13:10 on the 23rd. The top panels shows that RCGP exhibits higher robustness than the standard GP. but still affected by the outliers, when trained on data from the 23rd. RRP is virtually unaffected. Notably, RCGP relies on an a-priori weighting of data points based on the target values' proximity to their median, which can be counter-pron produci-r when the outliers are not a-priori separated in range. To highlight this, we included the previous trading day into the training data for the bottom panels, leading RCGP to assign the highest weight to the outlying data points due to their proximity to the target values' median, thereby leading RCGP "trust" the outliers more than any inlier,

---

resulting in it being less robust than standard GP in this scenario. See Appendix D.6 for additional
comparisons RCGP, on data sets from the UCI machine learning repository 34
(a) 14725
(b) 14725
14700
14700
DJIA Index
DJIA index
14675
14675
14650
14650
RCGP [GPR Noise)
14625
Relevance
14625
Priur
14600
x
Sel off
14600
1.0
1.0
Wena
0.5
0.5
10:00
12:00
14:00
16:00
12:30
12:50
13:10
13:30
13:50
Time
Time
(c) 14800
GP
(d)
(GPR Noise!
14700
14700
Relevance Pursult
index
index
14600
14650
14500
14600
14400
1.00
1.00
Wima
0.75
0.75
12
14
16
18 08
10
12
Time
Time
Figure 5: Results on intra-day data from the Dow Jones Industrial Average (DJIA) index on April 22-23 2013, which includes sharp drop at 13:10 the 23rd, see (b) for detailed view. The accompanying panels labeled Wing show the function that Altamirano et al. 2 RCGP uses to down-weight weight data points. Top: RCGP, exhibits higher robustness than the standard GP, but still affected by the outliers. The RRP model is virtually unaffected. Bottom: Including the previous trading day into the training data in (c), leads RCGP assign the highest weight Wing to the outlying data points due to their proximity to the target values' median, thereby leading RCGP to be even more affected than standard GP, see (d) for detailed view of the results the data of April 23.
6.2 Robust Bayesian Optimization
GPs are commonly used for Bayesian optimization (BO), which is popular approach to sample- efficient black-box optimization (24). However, many of GP models used for BO are sensitive to outliers and may not perform well in settings where such outliers occur. While Martinez-Cantin al. [6 consider the use Student-/ likelihood for BO with outliers, the use of other robust GP models has not been thoroughly studied in the literature.
Experimental setup We use Ament et al. variant of the LogEl family of acquisition functions, 32 replications, and initialize all methods with the same quasi-random Sobol batch for each replication. We follow Hvarfner et al. (31] and plot the true value of the best in-sample point according the GP model posterior at each iteration. We also include Sobol and an "Oracle", which is a Standard GP that always observes the uncorrupted value, and consider the backward canonical version of relevance pursuit, denoted by RRP, for these experiments. The plots show the mean performance with a bootstrapped 90% confidence interval.
(qLogNEI), a
Synthetic problems We consider the popular 6-dimensional Hartmann test function with three different corruption settings: (1) a constant value of 100, (2) U [1[3,,3 3] distributed value, (3) the ob- jective value for randomly chosen point in the domain. The results for 10% corruption probability shown in Fig. 6 We also include results for 20% corruption probability in Appendix D.3

---

Constani outleers value 100 (10%)
Uniformi 3) outliers (10%)
Uniform input outliers (10%)
of Bee Inker
-1.
Trimmeed
Sobol
-20
Standard GP
-25
Student-4
-30
RRP
Oracle
15
20
25
30
35
40
15 20 25
35 40
40
Number of Evaluations
Number of Evaluation
Namber of Evaluations
Figure 6: BO results for Hartmann6: Left: Relevance pursuit performs well in the of constant outliers of value 100, almost as well as the oracle. Middle: Relevance pursuit performs the best followed by the Student-t likelihood in the case of U1-3, 3]. Right: Similar the middle plot, this setting hides the corruptions within the range of the function, making it a challenging task.
Real-world problems We include three real-world problems: 3D SVM problem, 5D CNN problem, and a 20D rover trajectory planning problem, see the App. A for details. For SVM and CNN, we simulate random corruptions corresponding to I/O error, which causes the corresponding ML model to be trained using only a small subset of the training data. For the rover planning problem follow the setup in 47 with the main difference that we consider 20D trajectory, and the corruptions are generated randomly, causing the rover to break down at an arbitrary point along its trajectory. In most cases, this results in a smaller reward than the reward of the full trajectory.
emoo with prbb.
CNN, emor with prob. Rover, crashes with prob
Value Beg Inferee
2a18
0.9
o.17
Trimmedd
16
Sohol
1015
Standard GP
0.14
04
Student-t
RRP
9.15
10
15
20
25
30
0.2
10
15
25
40
75 100 125 150
Orncle
Number of Evalustions
Number of Evalnations
Number of Evaluations
Figure 7: BO results for three real-world problems: Left: RRP is competitive with the oracle on 3D SVM problem. Middle; The power transform performs best on 5D CNN problem, outperforming RRP well as the Oracle. Right: RRP performs well the 20D Rover problem.
Conclusion and Future Work
Contributions Robust Gaussian Processes via Relevance Pursuit (RRP) provides novel and prin- cipled way to perform robust regression. It permits efficient and robust inference, performs well across a variety of label corruption settings, retains good performance in the absence of corruptions, and is fiexible, e.g., can used with any mean or kernel function. Our method can be readily applied to both robust regression problems well applications such Bayesian optimization and is available Etrough[BoToros
13]. Importantly, it also provides thcoretical approximation guarantees.
Limitations As our approach does not explicitly consider the locations of the data points in the out- lier identification, may be outperformed by other methods if the underlying noise is heteroskedastic and
On the other hand, those methods generally do not perform well in the presence of sparse, ieaapennnnddeeedennt data corruptions.
Promising extensions of this work include performing Bayesian model averaging, i.e., average the predictions of the different possible sparsity models according to their likelihoods instead of using a MAP estimate, applying RRP to specialized models such as Lin et al. [41's scalable learning-curve model for AutoML applications, and Ament [ model for sustainable conerete. On a higher level, the approach of combining greedy optimization algorithms with Bayesian model selection and leveraging convex parameterization achieve approximation guarantees might apply to other parameters that are optimized using the MLL objective: length-scales of stationary kernels, coefficients additive kernels, inducing inputs, and even related model classes like Tipping (60] Sparse Bayesian Learning (SBL), which seeks to identify sparse linear models and is intimately linked greedy matching pursuits ( Overall, the approach has the potential to lead to theoretical guarantees, insights, and performance improvements widely-adopted Bayesian models.
Extensions
10

---

References
[1] Algikar and L. Mili Robust gaussian process regression with huber likelihood, 2023. URL
https://arxiv.. org/ awtp/3/11.111 07858
[2] M. Altamirano, F.-X. Briol, and J. Knoblauch. Robust and conjugate gaussian process regression. arXiv preprint arXiv:2311.0046. 2023.
[3] S. Ament and C. Gomes. On the optimality of backward regression: Sparse recovery and subset selection. In ICASSP 2021 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 5599-5603 2021. doi: 10.1 I0.10091CAASP33922.000 9415082.
[4] Ament and C. Gomes. Generalized matching pursuits for the sparse optimization of separable objectives. In ICASSP 2022 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 9747510.
2022. doi: 10. (0.1009/CAASP44922..00000
[5] S. Ament and M. O'Neil. Accurate and efficient numerical calculation of stable densities via optimized quadrature and asymptoties. Statistics and Computing, 28:171-183 2018.
[6] S. Ament, M. Amsler, D. R. Sutherland, M.-C. Chang, D. Guevarra, A. B. Connolly, J. M. Gregoire, M. O. Thompson, C. P. Gomes, and R. B. van Dover. Autonomous materials synthesis via hierarchical active learning of nonequilibrium phase diagrams. Science Advances, 7(51); eabg4930, 2021. doi: 10.1 16/./c/aaaa.. URL https: wttp///wwwww.ce..e/.//////// org/ hogg/dii/aa /bs/ 10. www1://////.v...ooe
[7] S. Ament, S. Daulton, D. Eriksson, M. Balandat, and E. Bakshy. Unexpected im- provements to expected improvement for bayesian optimization. In A. Oh, T. Nau- mann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neu- ral Information Processing Systems, volume 36, pages 20577-20612. Curran Associates, Inc., 2023. URt.jhttpp::/prrcceediag http:://wwoocceeeeeoeeeereeeeccs www..is..ccppppppeppppeeeeeeeeee haaplaapprellsseeeeeee files/ee//ee Lttle///ppprr222///110 0187722865aa62181111336622222222 Porereeee....eee..e.cccereeeeee. pdff
[8] S. Ament, A. Witte, N. Garg, and Kusuma. Sustainable concrete via bayesian optimization, 2023. URLLhttps::/arriv..oo
18288
[9] E. Ament and C. Gomes. Sparse bayesian learning via stepwise regression. In M. Meila and T. Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 264-274, PMLR, 18-24 Jul 2021. URLhtttps http:/ts//ppocceeingg..... mlr. aitess/v139/amentn
html
[10] D. Andrade and A. Takeda. Robust Gaussian process regression with the trimmed marginal likelihood. In Evans and Shpitser, editors, Proceedings of the Thirty-Ninth Conference on Uncertainty in Artificial Intelligence, volume 216 of Proceedings of Machine Learning Research, pages 67-76. PMLR, 31 Jul-04 Aug 2023. URL ptttps://proceeeee press/v216/ andrade222 1ttmmmmmm html
mll
[11] Awasthi, A. Das, W. Kong, and R. Sen. Trimmed maximum likelihood estimation for robust generalized linear model. In A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, editors, Advances in Neural Information Processing Systems, 2022.
http://wwwwwwweeeteeeeeeeee net/forum?
id= #d=VHmdFPy4U.u
[12] E. Bakshy, L. Dworkin, B. Karrer, K. Kashin, B. Letham, A. Murthy, and S. Singh. Ae: A domain-agnostic platform for adaptive experimentation. In NeurIPS Workshop on Systems for Machine Learning, 2018.
[13] M. Balandat, B. Karrer, D. R. Jiang, S. Daulton, B. Letham, A. G. Wilson, and E. Bakshy. BoTorch: A Framework for Efficient Monte-Care Carlo Bayesian Optimization. In Advances in Neural Information Processing Systems 33, 2020.
[14] E. Bodin, M. Kaiser, Kazlauskaite, Dai, N. Campbell, and C. H. Ek. Modulating surrogates for Bayesian optimization. In H. D. III and A. Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 970-979. PMLR, 13-18 Jul 2020. URI. https // /proceedings mlr. wwwpppreeiiiiii saress//119/bbdin20e html
[15] G. F. Box and D. R. Cox. An analysis of transformations. Journal of the Royal Statistical Society Series B: Statistical Methodology, 26(2):211-24, 1964.

---

[16] E. J. Candes. The restricted isometry property and its implications for compressed sensing. Comptes rendus. Mathematique, 346(9-10):589 2008.
[17] X. Chu, F. Ilyas, S. Krishnan, and J. Wang. Data cleaning: Overview and emerging challenges. In Proceedings of the 2016 international conference management of data, pages 2201-2206, 2016.
[18] P. Congdon. Representing spatial dependence and spatial discontinuity in ecological epidemiol- ogy: scale mixture approach. Stochastic Environmental Research and Risk Assessment, 31(2): 291-304, 2017.
[19] A. Daemi, H. Kodamana, and B. Huang. Gaussian process modelling with gaussian mixture likelihood. Journal of Process Control, 81:209-222 220, 2019.
[20] A. Das and D. Kempe. Approximate submodularity and its applications: Subset selection, sparse approximation and dictionary selection. Journal of Machine Learning Research, 19(3): 1-34, 2018.
[21] L. C. Dixon. The global optimization problem: an introduction. Towards Global Optimiation
2. pages 1-15, 1978.
[22] K. Dong, D. Eriksson, H. Nickisch, D. Bindel, and A. G. Wilson. Scalable log determinants for gaussian process kernel learning. Advances in Neural Information Processing Systems, 30,
2017.
[23] E. R. Elenberg, R. Khanna, A. G. Dimakis, and S. Negahban. Restricted strong convexity implies weak submodularity. The Annals of Statistics, 46(66):3339 3568, 2018.
[24] P. I. Frazier. A tutorial on bayesian optimization. arXiv preprint
2018.
[25] J. R. Gardner, G. Pleiss, D. Bindel, K. Q. Weinberger, and A. G. Wilson. GPyTorch: Blackbox Matrix-matrix Gaussian Process Inference with GPU Acceleration. In Proceedings of the 32Nd International Conference on Neural Information Processing Systems, NIPS' 18, pages 7587-7597, USA, 2018. Curran Associates Inc.
[26] R. Gameet. Bayesian Optimization. Cambridge University Press, 2023. to appear.
[27] J. Gerritsma, R. Onnink, and A. Versluis. Yacht hydrodynamics. UCI Machine Learning Repository, 1981. DOI: https://doi.org/ http://di.orr/1..243442332 2.344444/2427..
[28] W. Goldberg. C. K. Williams, and C. Bishop. Regression with input-dependent noise: A gaussian process treatment. In Advances in Neural Information Processing Systems, volume 10, January 1998.
[29] J. Harrison, David and D. L. Rubinfeld. Hedonic prices and the demand for clean air. Journal of Environmental Economics and Management, 5:81-102, 1978.
[30] R. A. Horn and C. R. Johnson. Matrix analysis. Cambridge university press, 2012.
[31] C. Hvarfner, E. Hellsten, F. Hutter, and L. Nardi. Self-correcting bayesian optimization through bayesian active learning. Advances in Neural Information Processing Systems, 36, 2024.
[32] Jylänki, J. Vanhatalo, and A. Vehtari. Robust gaussian process regression with a student-t likelihood. Journal of Machine Learning Research, 12(11). 2011.
[33] O. Karaca, D. Tihanyi, and M. Kamgarpour. Performance guarantees of forward and reverse greedy algorithms for minimizing nonsupermodular nonsubmodular functions on a matroid. Operations Research Letters, (96)):85--8611 2021.
[34] M. Kelly, R. Longjohn, and K. Nottingham. The uci machine learning repository. URL https://archive. ics. uci. edu, 2023.
[35] K. Kersting, C. Plagemann, P. Pfaff, and W. Burgard. Most likely heteroscedastic gaussian process regression. In Proceedings of the 24th International Conference on Machine Learning, ICML '07, pages 393-400, 2007.
[36] A. Krause and D. Golovin. Submodular function maximization. In L. Bordeaux, Y. Hamadi, and P. Kohli, editors, Tractability, pages 71-104. Cambridge University Press, 2014. ISBN 9781 13911991778 URI, http:://bb/ uni trier dwtpdd/dddbbbbeebeeeee
htm #0001G14
12

---

[37] A. Krause and C. Guestrin. Near-optimal nonmyopic value of information in graphical models. In Proceedings of the Twenty-Firss First Conference on Uncertainty in Artificial Intelligence, UAT05,, page 3224-331, Arlington, Virginia, USA, 2005. AUAI Press. ISBN 0974903914.
[38] A. Krause, A. Singh, and C. Guestrin, Near-optimal sensor placements in gaussian processes: Theory, efficient algorithms and empirical studies. Journal of Machine Learning Research, 9(2), 2008.
[39] M. Kuss. Gaussian process models for robust regressilon, classification, and reinforcement learning. PhD thesis, Technische Universität Darmstadt, 2006.
[40] Z.-Z. Li, L. Li, and Z. Shao. Robust gaussian process regression based on iterative trimming. Astronomy and Computing, 36:100483, 2021. ISSN 2213-1337. doi: https///doi.org/10..10 1016/ j.ascom.2022 2021..0011...... UtI.htps::/ w/wp:////wwwww..//o wwwp://wwwwwccccccccccooo
$2213133721003778
[41] J. A. Lin, S. Ament, M. Balandat, and F. Bakshy. Scaling gaussian processes for learning curve prediction via latent kronecker structure, 2024. URL httpt::/tp//trtiii org/abs/ a////4/////100 09239
[42] D. C. Liu and Nocedal. On the limited memory bfgs method for large scale optimization. Math. Program., 45(1-3):5003
Aug. 1989. ISSN 0025-5610.
[43] F. Locatello, R. Khanna, M. Tschannen, and M. Jaggi. A Unified Optimization View on Generalized Matching Pursuit and Trank-Woffs In A. Singh and J. Zhu, editors, Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, volume 54 of Proceedings of Machine Learning Research, pages 860-865 PMLR, 20-22 Apr 2017. URL https://proceedinng
mlr press/
[44] F. Locatello, A. Raj. S. P. Karimireddy, G. Raetsch, B. Schölkopf, S. Stich, and M. Jaggi. On matching pursuit and coordinate descent. In Dy and A. Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 3198-3311 3207. PMLR, 10-15 Jul 2018. URLLhttps:/pprocedding html
mlr. ppresss/80/oocate///ooo
[45] S. Lotfi, P. Izmailov, G. Benton, M. Goldblum, and A. G. Wilson. Bayesian model selection, the marginal likelihood, and generalization, 2023. URL httpp::::///arti http:/aaxivv...g/bbe./// org/aba/// 2202. 11678
[46] R. Martinez-Cantin, M. McCourt, and K. Tee. Robust bayesian optimization with student-t likelihood, 2017.
[47] N. Maus, K. Wu, D. Eriksson, and Gardner. Discovering many diverse solutions with Bayesian optimization. In Proceedings of the International Conference on Artificial Intelligence and Statistics, volume 206, pages 1779-1798. PMLR, Apr. 2023.
[48] A. Naish-Guzman and S. Holden. Robust regression with twinned gaussian processes. In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Syy- lems, volume 20. Curran Associates, Inc., 2007. URL whtp:hhtt/////momcc. httpps:/pproceeednnggg neurips. cc/ paper_filest /tiles//aperr2007//11e/////
pdf
[49] R. M. Neal. Monte Carlo Implementation Gaussian Process Models for Bayesian Regression and Classification. Technical Report 9702, Department of Statistics. University of Toronto, 1997.
[50] G. Nemhauser, A. Wolsey, and M. L. Fisher. An analysis of approximations for maximizing submodular set functions-2 i. Mathematical programming. 14:265 294, 1978.
[51] R. K. Pace and R. Barry. Sparse spatial autoregressions. Statistics & Probability Letters, 33(3): 291-297, 1997.
[52] C. Park, D. J. Borth, N. S. Wilson, C. N. Hunter, and F. J. Friedersdorf. Robust gaussian process regression with a bias model. Pattern Recognition, 124:108444, 2022.
[53] R. Ranjan, B. Huang, and A. Fatchi. Robust gaussian process modeling using em algorithm. Journal of Process Control, 42: 125-136, 2016.
[54] C. E. Rasmussen, C. K. Williams, et al. Gaussian processes for machine learning, volume 1. Springer, 2006.
[55] C. Riis, F. Antunes, F. Hüttel, C. Lima Azevedo, and F. Pereira. Bayesian active learning with fully bayesian gaussian processes. Advances Neural Information Processing Systems, 35: 12141-12153, 2022.
13

---

[56] A. Shah, A. Wilson, and Z. Ghahramani. Student-t processes alternatives to gaussian processes. In Artificial intelligence and statistics, pages 877-885. PMLR, 2014.
[57] N. Srinivas, A. Krause, S. Kakade, and M. Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. In Proceedings of the 27th International Conference on International Conference Machine Learning, ICML'10, page 1015-1022, Madison, WI, USA, 2010. Omnipress. ISBN 978 16005890077.
[58] S. Sundararajan and S. Keerthi. Predictive approaches for choosing hyperparameters in gaussian processes. In S. Solla, T. Leen, and K. Müller, editors, Advances in Neural Information Pro- cessing Systems, volume 12. MIT Press, 1999. URL http:hhcps//p/roeeeeeiia wwteer_files/pppppe
neurips cc/ Paper pdf [59] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):222 (8)6)2278) 1996.
[60] M. E. Tipping. Sparse bayesian learning and the relevance vector machine. Journal of machine learning research, 1(Jun):211-22) 2001.
[61] M. E. Tipping and N. D. Lawrence. Variational inference for student-t models: Robust bayesian interpolation and generalised component analysis. Neurocomputing, 69(1):123-14) 2005.
[62] J. A. Tropp. Greed is good: Algorithmic results for sparse approximation. IEEE Transactions on Information theory, 50(10))2231-2244 2004.
[63] A. Tsanas and Xifara. Energy Efficiency. UCI Machine Learning Repository, 2012. DOI:
https://doi.org/ https:/ddiiorr//0.2445///// 10/1114///42/513000
[64] Vanhatalo, P. Jylänki, and A. Vehtari. Gaussian process regression with student-t likelihood.
In Advances in Neural Information Processing Systems, volume 22, 2009.
[65] 7. Wang, C. Gehring, P. Kohli, and S. Jegelka. Batched large-scale bayesian optimization in high-dimensional spaces. Proceedings of the International Conference on Artificial Intelligence and Statistics, volume 84 of PMLR, pages 745-754. JMLR, Mar. 2018.
[66] L. Wasserman. Bayesian model selection and model averaging. Journal of mathematical psychology, 44(1):92-107, 2000.
[67] D. Wipf and S. Nagarajan. new view of automatic relevance determination. In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Sys- tems, volume 20. Curran Associates, Inc., 2007. URL https:://procceeiinnngge//////e wtwpeeeeddddnnnnnnnnneeeeee newrips...... cc/
[68] I.-C. Yeh, Concrete Compressive Strength. UCI Machine Learning Repository, 1998. DOI: https//0oo.or//002433////CK
14

---

A Additional Details on the Model
Algorithm bellow is the "backward" variant of Algorithm T from See 4 As its name suggests, the main difference compared to the "forward" variant is that rather than building up set of "outliers", it starts from a (typically large) set of "outliers" and iteratively removes those data points from the set that have the smallest inferred data-point-dependen noise variance Pt-
While we have not derived theoretical guarantees for this "backward" version, we have found it to generally behave similarly to the "forward" version in terms of performance and robustness. One empirical observation from studies is that while the "forward" version tends perform slightly better than the "backward" version if there are only few outliers, the opposite is true if the outlier frequency very high. This behavior is rather intuitive and illustrates that relevance pursuit is particularly well-suited to identify sparse, low-cardinality subsets (note that in the "backward" variant under large corruptions, the uncorrupted data points can be viewed as the sparse subset that needs to be identified).
Algorithm 2 Relevance Pursuit (Backward Algorithm)
Require: X. y, schedule K = (k1, ......) Initialize Si S {1. {I..... (typically Sj (1.....mI)
for k; in IC do
Optimize MI.: Pīş. +- aag max psy ((Ps;), where Ps: = (p: Pj =0 v & S Compute the set R2 containing the k; elements of St with smallest inferred variance:
Rí !.t where j! E St such that Pis (11) P'sc (il') for I < I
S + S5 R₂
S; + (1,... n) S for each k; Compute the marginal likelihood y) = p(S.. Ps.. ()
S* + arg maxs, p(S, y))
Return S* P's..
B Additional Background on the Theory
B.1 Submodular Functions
Krause and Golovin 3d provides survey on the maximization of general submodular functions. Here, we focus on applications of submodularity to sparse regression and Gaussian process models.
Sparse Regression Das and Kempe [20] showed that the subset selection problem of regression features with an R² objective satisfies a weak submodularity property. which can be invoked to prove approximation guarantees for the greedy maximization of the objective. Elenberg et al. 2 generalized this work by proving that any log likelihood function exhibiting restricted strong concavity gives rise to the weak submodularity of the associated subset selection problem, which can be invoked to prove approximation guarantees for the greedy algorithm. Karaca et al. [33 contains a guarantee for the backward algorithm applied to the maximization of submodular set functions.
Gaussian Processes Submodularity has also found application to Gaussian process models. For sensor placement problem, Krause et al. (38 proved that the mutual information (MI) criterion, capturing the reduction in uncertainty in the entire search space, can be a submodular function. In this case, MI is not monotonic everywhere, but monotonic for small sets (2k) of sensors, which is sufficient to apply Nemhauser's guarantee for sparse sets of sensors up to size & [50]. Relatedly, the "myopic" joint entropy of set of observables is unconditionally submodular as consequence of the "information never hurts" principle (37). but generally leads to lower-quality sensor placements than the MI criterion. Srinivas et al. [57] used the submodularity of the joint entropy in order to prove regret bounds for the convergence of a Gasbased BO algorithm using the upper-confidence bound acquisition function.
Elenberg et al. 23 proved that any log likelihood function exhibiting restricted strong concavity and smoothness implies the weak submodularity of the associated subset selection problem.
15

---

Definition 9 (Submodularity Ratios (23). Let A. C [n] be two disjoint sets, and f 2( - R. The submodularity ratio of B with respect to is defined by
7B.A E (it) f(B)) ///////////AA/ f(())
VEA
The submodularity ratio of set C with respect to integer k is defined by AnB 0,
7C.k min 78,4 B.A
such that
B BCC,
and
| S k.
Then given Y > a function 7--eakly submodular set C with respect if,k 2 %-
Theorem 10 (Weak Submodularity via RSC (23)). The submodularity ratio 't8.k can be bound below using the restricted convexity and smoothness parameters missi and Mi
78, 2 mis++ + / MiSI+K-
Theorem 11 (OMP Approximation Guarantee (23)) Let (rMP be the r-sparse vector selected by OMP, and xopr(r) arg mi |o fxx) be the optimal r-sparse vector. Then
(xOMP(r) 11-e-ma.Ma.. [xOOrT(r))
where mar, M2, are the restricted strong convexity and smoothness parameters of f, respectively.
C Theoretical Results and Proofs
Lemma 1. /Optimal Robust Variances] Let D {(x, Mj) : + i}. Pye + pee. where
P. P/. E RT. (p O, and e, is the ith canonical basis vector. Then keeping P., fixed,
P arg max C(p,s + pre.) [(y: (Dy( Vlu(x:)
(3)
P
where y(x,) = =f fxx)) te,. These quantities can be expressed as functions of -- = (K (K+D)>+p)
and
Viy(x;) D = 1/ "m
where Dø²+p is a diagonal matrix whose entries are o² +p.
Proof. First, we partition the covariance matrix + o Dp to separate the effect of p and use the block matrix inverse
Σ + +uuu
-13,
E = =(kK xK+D)t=>)
(5)
-u
B,
where
Σ, 1(X,,, X(:) Pitto
ET k(X,,, (x)x,
and
(6)
8, (1((),,,) x4) 0² + pi] Xx,,, X (an_)((Xn k(X,,; xi)
Quadratic Term of the log likelihood
With the expression for the inverse of Σ above, we can write the quadratic term
yT (K Dpto:) y - u3 Jy 2y + v
-y y + B,y y:)
Determinant Term
The determinant of block matrix is given by
[
B =
D
CA -B.
(8)
D
Applying this identity Σ (K
D Dpto2), we get
K
D
p+oo = 13
(9)
16

---

Log Marginal Likelihood
2(((p) - = - y1)² - log(3,"),
(10)
Noting that 0, = B/ the derivative of the difference in log marginal likelihood w.r.t. Pi, is
pp. (((P) - ((),)) (y y Bc.
(11)
While 3, = 0 is root of the derivative, we ignore this solution since 3; is never zero when o² > 0. Therefore, the remaining stationary point is 3,- = (y,,"" y.)². Since we constrain / 2 0, this point might not always be attainable. However, because there is only single stationary point with respect 8, when o² 0, and 3, is a strictly decreasing function of pt. follows that the marginal likelihood is monotonic as function of to both the left and the right of the stationary point. Therefore, the optimal constraint pj simply the optimal unconstrained value, projected into the feasible space. In particular, solving 3 (y u y1)² for pt and projecting to the non-negative half-line, we get
P: (y vi)² -Ax(x,,x)++ + - 1xx, X M ( (x.
(12)
Lastly, note
y u y/ vT(((( xx) +D p ) 11(1 1/ x;) - y/ E(y(x; 1D y/:
(13)
x,x,,,;,) k(x;; X ) 1: ( xxxx + + V(b(x;) |Dyll=
As stated by Rasmussen et al.
54 (P. 117, Eq. 5.12), originally shown by Sundararajan and Keerthi 58 these quantities can be expressed as simple functions of Σ Σ-
E (ly((;) | y ΣΣ y]./ Σ
1m
(14)
(P 1/ | :s
Therefore, all LOO predictive values can be computed in (nn)) or faster, if an inducing point method used for K.
The following is a preliminary result for our analysis of the log marginal likelihood w.r.t. p.
Lemma 12. The gradient and Hessian of the log marginal likelihood C with respect to p are given by and -HH = (2xxx² T K K-1
-2Vp(4) = diag(K -axax
where K Ko for some base covariance matrix Ko and a K 'y.
Proof. Let = K- Regarding the gradient, note that 3, K e where e, is the canonical basis vector with a one as the ith element, and based on Equation 5.9 of Rasmussen et al. [54).
-20p, [4] = tr (KK 1 axx 0, K)
Lr(( axc eeee.)
=e_ (K-1 aa T) ei
-KK-1- -xxx it
Regarding the second derivatives, according to Dong al.
[22]
0o Oo, [log K|| = tr(kk-rkk))- 100, do K K [0,, K) (e0KK K])
0n, 0の, jy K K-ty| 2a" [0, K|K 100,, K|a - a 10p, Oo, K|a.
Therefore.
-20, 0, L = tr (KK ax )1,, 0, K| - (K-1 - -2a T)(0,, KIK Kl))
tr ((2xx T - K-1 ) ( 11,0;0,]
tr (e] (2xx - -K
KK -K-lee; )e:) KK (1()
[2axa T - 11- 1is
[(200 - -K-1) 3 K-
is
since (3p. O,, K] O. The third equality due to the invariance of the trace to circular shifts of its argument. The forth equality is due to the symmetry of the matrix in brackets.
17

---

C.1 Strong Convexity and Smoothness of the Reparameterized Robust Marginal Likelihood
Here, we re- re-parameterize p(s) diag Ku) O ((1-)) -1), and attain strong convexity for all inputs s, if conditions on the cigenvalues of the covariance matrix and the norm of the data vector | are met. The convexity result is surprising two ways: the negative log marginal likelihood of GPs is generally a non-convex function, and in addition, the negative log likelihoods of many alternative robust regression methods like the Student-/ likelihood or a-stable likelihoods are non-convex, and even Huber's proposal is non-strongly-conner
Lemma [Reparameterized Hessian] Let Ks x(X, + o²I D Dp(s) Ks
diag ( )11a 1, diag(K, (Kg) and â = K₂¹ diag(Ks) Then
H [L[[2[(p(s) = D (2 (x o (K-1 - r)) + diag(K< - - O O
Proof. Using the chain-rule, the Hessian Ha [4] can be expressed as a function of the Jacobian J. [[(ss] Dap(s). which is diagonal since p(s) is element-wise function, and the second derivatives (s;)) Then
-2H, = Jslp] H2l2 [],p[p] + Dv 2cj
D ( K ) K 'D' p'(s) diag(K (BoKKKKKDpp
DD" 1D"--ppp -Dp (KKKKK''- (K- ..-o--- 1DD
where substituted the relevant expressions from Lemma 12 Further substituting p(s);
[Ko]:](( (1 - 8s)_ 1), (((s); [o] - 8;)-², and p (8) = 2[Ko]a(11 8;)= noting
that K = Ko- (ing(Kn)(( -88 manipulation finish the proof.
(Ko (llag(Ko) ) diag(Kn)( - -s)-1,
1,)-1,a0 and algebraic
Lemma 5. [Strong Convexity via Eigenvalue Condition] Let Kg as in Lemma
Then Hg m if
(2xmLx 2-2
m)
Amin min
min
ly 12.
2(1 - Aminn/ >maxx)
where maix (resp. Amiin,max) are the smallest and largest eigenvalues of Kg, respectively Ky,
Proof. We seek to lower-bound the smallest eigenvalue of the Hessian matrix, which-for twice- aantimoously-difffeentaaa problems-is equivalent lower and upper bounds of the problem's Hessian matrix. Starting with the result of Lemma
H. H.(-2((p(s) = D 2 (ara" ② (K-11-1) r)) + 2diia((k))) - (iK'kokk) D
z 2 aa T o (K-1-1)) + diag((K)) (K- + k-
Now, we bound each of the three additive terms independently from below.
Term 1:
diag( -1-
1)) > 2Amin (K- 22 -K)
The first inequality comes from K being positive definite. and the absolute value of the diagonal of matrix, which is already positive for positive definite matrices, being lower bounded by the minimum eigenvalue of the matrix. The last steps a basic consequence of the eigenvalues of inverses matrices. Note that the eigenvalues of K can be further bound by the eigenvalues of the original matrix K:
Amin (K) = Amin (diag K) 1/2 K diag (K) (ii)
Amin (diag (K) (K) 11/2/3 Amin (K)
Amiu (diag (K) ) (K)
Amnin
Amin 1-11)1111)) Amii (K)
Amia (K) Amaxx (K).
In a similar way, we can show that Amax (aaxx()a. S Amaxx (K) Amin (K), which implies Amax Amax((K)1 (K)-12 Amin(K) (K)/ /ammmxx (K).
18

---

Term 2: Next, we have
O -K
mmax (K-1)21= -
Amin (K) II
which is due to the Hadamard product being sub-matrix of the Kronecker product of the same matrices, the largest eigenvalue of the former are bounded by the largest eigenvalue of the latter, which is the product of the largest cigenvalues of the constituent matrices.
Term 3: Lastly, note that
2 (ara" @k- 1)) = 2D -IID
= 2D (Amin (K I)Da
2(Amin (K 1) 1)D2
2(mmax 1 1)D
> 2(mmiu (K) / (K) 1)D2
= 2(Amin (K)/Amax(() (K)
1)
where the last inequality comes from the second last lower bound of Term 3 being non-positive, and therefore, being able to lower bound it with the largest magnitude entry of &
Term + 2 + 3: Putting together the inequalities for all terms, get
Amin (H,)= 2 2Amax (KK-1- Amin (K) -2 + 2(Amin (K) /max/ (K) 1)]
where in slight abuse of notation, we let H3 = Hal-2C] be the Hessian of the negative log likelihood.
Using the bound || < ||x]|2 S Amin 1|)-1 where y = diag(K)-1/2y Therefore, || S Amin (K)-1 Ami (K)-1/2|[y]]) I Amax (K) (K)-3/21ly]]>>
Finally, lower bounding the current lower bound by m > 0 yields sufficient condition for the convexity at S.
2) 2)maxx (K)-1- Amin (K)-2
Amin(K)) + 2(Amin (K) Aumax (K) - 1) (K)-22 Amii (K)-"|Y> > m.
Re-arranging, we attain
Amin 82 (2)max 2-2 -m) m)
11 12.
(15)
2(1
Aminn /mmax)
which is non-trivial guarantee when 2i-1
min >> O.
Lemma 13. /Smoothness via Eigenvalue Condition] Let K in Lemma 8max and
Suppose that ||s||>o S
(M(1 +max) )2 + min 28-11
Amin min
1
(16)
2(Amax. /miin
where Amin, (resp. Amin.max
m ) are the smallest and largest eigenvalues of K (resp. K) respectively.
Then H2 x M. This is non-trivial guarantee when M smax) )2 + 2-2
2 > 0.
Proof. We now derive an equivalent upper bound for the largest eigenvalue of the Hessian. Starting with the result of Lemma
H₂| (al-2[(p(s)] D 2 (àà T o (K- -K) 2daag(K (K--1 の k- D
(1 s) '13 [2 (àā o (K (K-- -P)-1) 1)) + diag (K- (K (K- K ')]
Therefore, it becomes immediately apparent that we will need to introduce upper bound on [[(1 - -) 12. which restriction on the domain that p can take. Proceeding similar way as above, bound the three terms in square brackets, now from above.
19

---

Term 1:
diag(K- 1 2Amax M (K- 2Amin 2mmin(K)-11 (I)
Term 2:
o ) 1 -mmin (K 1)2 -Amax (K) - - (K ② 1
Term 3:
2 (a @ (K-1-11)
1)) = & -I)Da
= 2D (Amax (K 1 11 1 I)Da
2mmmx (K - 1)DD3
2( KK) -1 11011111110 1)D2
3 2(Amax (K) (1a)minnn-nn 1)D3
3 2(2
(K) - I
where now the last inequality follows because the second last expression is always non-negative. Term + 2+ + 3: Putting together the inequalities for all terms, we get
Amax (H₂) S 22min (K)-1 - Amax (K) (K)- 2(Amax (K) Amin (K) d)
Finally, upper bounding the current upper bound by M/Il(1 s) 1| > 0 yields sufficient condition for the M-smoothness result,
Using the same bound for | derived for the convexity
2Amin (K) (K)-1 Amax (K) + 2((2((())))) (K) Amin (K) - 1) (K)- (K) -1yl22 < M.
Re-arranging, we attain
Amiu 2 (M(1 Smaxx )2 + 2-2 -imixx
1ly
(17)
2( 2(mmax Admia - 1)
which is non-trivial guarantee when M(1 - Smax) )² + i 2-2
22m
max > 0.
Lemma 7. [Strong Convexity via Diagonal Dominance] Let m > 0 and Ko be 8-diagonally dominant with < ((5 - m) - 25 - 9m + 17) /4 S (5 (17) F2 0.44 and
Amin (Ko) - -1)) 2(1 01-11/1111111 d)) m N 1
Then the NMLL m-strongly convex for all [0, 1)" iee. p(s) E o0)
Proof. Fist, Gershgorin's Disk Theorem implies that the eigenvalues of diag (A)-1/2 A diag(A (A) -1/2
lie in (1 - 8, + d) for 8-diagonally dominant matrix A. Further, the condition number of A is
bounded above by r((A) = Amax (amax(AA) (A) Am)mmin (A) I (1+ d) - 5). See Horn and Johnson [30 for more background matrix analysis. Plugging these bounds into the results of Lemma 5 vields
(2)_1x 2-2 m)
Amin (1) 1011-87 8)2 (2(1 8) (1 8) m)
Amin (2) min 2(1 Aunin/ A
2(1 (1 8)/(1 + 8))
Lower bounding the last expression by 11yl1 implies m-strong convexity. This gives rise non-trivial guarantee whenever the numerator is larger than zero. In particular,
2(1 + 5)-1 (1 5) m >
Expanding, noting that 0 < 8 < implies 83 < 82 in order to reduce power in the resulting expression, and collecting terms with like powers, we attain the following sufficient condition
28 + (m 5)6 + (1 m) > 0.
Note that for this condition to hold if = O, we need to have m < Fortunately, this is a quadratic in 5 whose smallest positive root is
1 (55 m)
8_ 4
(5 m)² 8(1 m)
20

---

In particular, for m = O, this reduces to
8_ = 1 V1 s 0.458447718799116977
Lastly, note that if Ko is 5-diagonally dominant, then so is K = Ku + Dp(s). since the robust variances add the diagonal, making it more dominant. Therefore, the convexity guarantee holds for all p(s), if it holds for the base covariance matrix Kg. Note that Amin (K) 2 o².
Similarly, we can prove similar statement relating diagonal dominance with M -ssmoothness. Lemma 14. [Smoothness via Diagonal Dominance] Suppose Kg is 6-diagonally dominant covari- ance matrix and suppose we constrain ||s||n0 S 8max I VI/M. Then
( (1)) 5) M(1 8max) - 2 |lyl
Amin
2((1 + 6)/(1 8) 1)
implies that the NLML is m-strongly convex for all E [0, 1]", i.e. p(s) oo)"
Proof. We proceed in similar way as for Lemma] 7 but with Lemma 13 as the starting point.
Amin (K) (Amin) (M(1 8 Smax) )2 + x-2 27_
2 Amin (K) 8)² M(1 Smax) ) 1 2((1 8)/(1 - 8) 1)²
min
2(Mmax/ /miin 1)
where we used (i-2
min - 2i 2 -1, which is tight when K is unitary. Lower bounding the last expression by llyl1 implies M -smoothsssss. This gives rise to a non-trivial guarantee whenever the numerator is larger than zero. In particular,
- >> O, which implies Smax I 11-v /I/M or equivalently. M > 1/(1 - $m Lastly, note that if Ko is ö-diagonally dominant, then so is Kg = Ko + Dp(s), since the robust variances only add to the diagonal. Therefore, if the incquality holds for the base covariance matrix Ko, the smoothness guarantee holds for all p(s) such that
8 S Smax < 1 VT/M. Note also that Amin (K) 2 0².
Theorem 8. [Approximation Guarantee] Let Ko k(X,X) g²I be 5-diagonally dominant, Smax > O be an upper bound on [|s|| and suppose ly|l. 5 satisfy the bounds of Lemmas [and Land[4 14 guaranteeing m-convexity and M ussmoothness of the NMLL some m > 0, M > 1/(1 Smax] ) Let (T) be the r-sparse vector attained by OMP on the NMLL objective for r steps, and let be the optimal r-sparse vector. Then for any 2r S n,
sOPT (r) = arg [aloo Semax C(p(s))
È (p(SOMP(r)) 2 (1-e-m/bl E (p(sop((r)).
where Z(-)= = C(-) ((0) is normalized that maxag [(ss) 2 0 for any support S.
Proof: The result is direct consequence of meeting the m-convexity and M-smoothness conditions of Lemmas 7 abooee and the OMP approximation guarantee of Theorem [T] due Elenberg et al. [23] Note that the condition on 2r < comes from the RSC condition in Theorem II being required for subsets size 2r. As we proved bounds for the m-convexity of the full Hessian of size n, has to be smaller than n for the assumptions of the theorem to hold. Regarding the upper bound 5max on s, we note that the constraint is convex and therefore doesn't 't change the convexity property of the optimization problem.
Further, note that maxps C(ps) S maxpsus C(psu:). since the additional non-zero zero element could stay at O, if the marginal likelihood does not improve with p; increasing. That is, the subset selection problem is monotonic. As a consequence, we can normalize the MLL by Z(-) = L(-) L()), which then only takes positive values for any s = arg maxx L(p(ss)), i.e. maxs; L(ss) 2 0. This normalization is required for the constant factor approximation guarantee to apply, similar to the original work of Nemhauser.
This theoretical approach could lead to approximation guarantees for Tipping 660 's Sparse Bayesian Learning (SBI.) model, for which Ament and Gomes [9] show that greedily optimizing the associated NMLL is equivalent to stepwise regression in the limit of - 0, proving exact recovery guarantees.
21

---

D Additional Detail on the Experiments
Our benchmark uses a modified version of the code from Andrade and Takeda [00, available at btttp://gitthte http://www...aaaaaaaaaeaaaaaaaaa GPLv2 license.
andrade-
under the GNU
D.1 Synthetic Regression Problems
Model fitting runtimes Fig. summarizes the fitting times for the different models on the different scenarios from Section 6. We observe that the outlier type and the fraction of outliers both have relatively limited effect on the fitting times for all of the models. Fig. therefore provides a more compact view of the same data (aggregated across outlier types and outhier fractions), giving a better sense of the distribution of the fitting times. Unsurprisingly, the baselines that simply single GP model ("vanilla", "winsorizze" approaches. While all of the robust models show similar fitting times on the Hartmann problem, fitting our RRP is significantly faster (note the logarithmic scale of the y-axis) than fitting the Student-f and trimmed MLE models on the 5-dim and 10-dim Friedman problems. The trimmed MLE model in particular ends up being quite slow, especially on the 5d Friedman function.
(taverrrtnnformm) are substantially faster than any of the robust
D.2 BO experiments, additional details
For Hartmann6, we consider the standard domain of [0, 116.
SVM The SVM problem, the goal is to optimize the test RMSE of an SVM regression model trained on 100 features from the CT slice UCI dataset. We tune the following three parameters: C [le-2, lc2], e E [le-2, 1], and f E [le-3, le-1]. All parameters are tuned in log-scale. Corruptions simulate 1/O failures in which case only train U [100, 1000] training points out of the available 50, 000 training observations.
CNN For the 5D CNN problem the goal is to optimize the test accuracy of CNN classifier trained on the MNIST dataset. We tune the following 5 parameters: learning rate in the interval [lc-4, le-1], momentum in the interval [O. I]. weight decay in the interval [O. 1], step size in the interval [1, 100), and [0. 1I. Similarly to the SVM problem, we only train on U1100, 1000] of the available training batches when an I/O failure occurs.
Rover The rover trajectory planning problem originally proposed in Wang al. 65 The
goal is to tune the trajectory of a rover in order to maximize the reward of its final trajectory. We the same obstacle locations and trajectory parameterization main difference being that we parameterize the trajectory using 10 points, resulting in 20 tunable parameters. When corruption occurs, the rover will stop uniformly random point along its trajectory, generally resulting in lower reward than the original trajectory.
in Maus al. (47). with the
D.3 Additional Synthetic BO Results
In addition to the results Fig. 6 we also include results when the corruption probability is 20% in Fig. I0 We observe that the results are similar with 10% corruption probability, but that the performance of several baselines regresses significantly.
D.4 Computational Setup and Requirements
Robust GP regression is very data-efficient, focuses on the small-data regime, and runs fast (faster than competing baselines studied in this paper). Therefore, each individual run required very limited compute resources (this includes the baseline methods). To produce statistically meaningful results, however, we ran large number of replications for both our regression and Bayesian optimization benchmarks on a proprietary cluster. We estimate the amount of compute spent on these experiments to be around 2 CPU years in total, using standard (Intel Xeon) CPU hardware. The amount of compute spent on exploratory investigations as part of this work was negligible (this was ad-hoc exploratory and development work on a single CPU machine).
22

---

friedman10.
friedman 10, student-t
friedman (0, uniform
10
10a
fit tittime (8)
10
friedman5,
friedman5, student-t
friedman5, uniform
10
107
fit time (8)
101
hartmann6, constant
hartmann6, student-t
hartmann6, uniform
10
10a
fit time (5)
101
0.05
0.10
0.15
0.25
0.05
0.10
0,20
0.25
0.05
0.20
025
Outlier Praction
Outlier Fraction
Outlier Fraction
Studeat-t
Standard
Trimmed
RRP (hkod)
RRP (hkwd, can)
Power Transf
Adapt.
Figure 8: Fitting times of the different robust GP modeling approaches. The plots show means and one standard deviation. Here "bkwd" indicates the backward variant of RRP (Algorithm 2 and "can" indicates the canonical (non-convex) parameterization.
D.5 Impact of Convex Parameterization on Joint Optimization of Hyper-Paraameeers
limitation of Lemma is that it only guarantees convexity for the sub-problem of optimizing the p. s. practice, we jointly optimize all GP byper-parameters, including length scales. In this case, the theory guarantees the positive-definitenes of the submatrix of the Hessian corresponding to p(s). and we expect this improve the quality of the results of the numerical optimization routines.
Indeed, positive-deffniteness is beneficial for quasi-Newton optimization algorithms like L-BFGS 42 which restarts its approximation the Hessian whenever it encounters non-convex regions, because the associated updates to the Hessian approximation are not positive-define definite. This leads the algorithm momentarily revert back gradient descent, with associated slower convergence rate.
To quantify the impact of this, we ran convergence analyses using the data from Fig. D allowing all p to be optimized jointly with other hyper- parameers (length scale, kemnll variance and noise
23

---

friedman5
10
fit time (8)
10a
10
10P
friedman 10
10
10
fit time
10
10c
hartmann6
100
fit time (5)
10
10
10c
RRP
RRP (bkwd, can)
Standard GP
RRP (bkwd)
Student-t Power Transf.
Adapt. Wins.
RRP (can)
Trimmed MILL
Figure 9: Fitting times of the different robust GP modeling approaches the regression tasks from Section Results are aggregated across outlier types and outlier fractions as those do not affect fitting times much (see Fig. 8).
Vilue of Best bstered Reir
Constaur cutliers value 100 (20%)
UniformE-3, 3) outlicus (20%)
Uniform imput outiees (20%)
-0.5
-1.0
Trimmied
-1.5
Sobol
Standard GP
-2.5
Student-t
-3.0
RRP
Ornele
Number
Figure 10: BO results for Hartmann6: Left: We see that Relevance pursuit performs well in the case of constant outliers of value 100 and almost performs just well as the oracle. Middle: Relevance pursuit performs the best followed by the Student-t likelihood from [46] in the case of U[-3, 3). No method performs as well as the oracle when the outlier probability is 20%. Right: Similarly the middle column, this setting hides the outliers within the range of the outliers making it difficult to match the performance of the oracle.
variance), recording the achieved negative log marginal likelihood (NLML) as a function of the tolerance parameter ftol of the L-BFGS optimizer. The results reported in Table T.5 and indicate that the optimizer terminates with much better NLML using the convex parameterization with the same convergence tolerance.
There are also settings which we do not actually jointly optimize the hyper-parameters, particularly when we have access to data from the same data generating process that has been manually labeled

---

ftol
Canonical p
Convex p(s)
Table 1: Comparison of the negative log marginal likelihood achieved after numeri- cal optimization of the canonical and con- vex parameterization of p with L-BFGS. No- tably, the convex parameterization leads to improved marginal likelihood values for for a given convergence tolerance ftol.
1e-3
-4.37
-14.18
le-4
-4.37
-93.00
le-5
-4.37
-93.01
le-6
-4.37
-135.05
1c-7 -98.62
-518.68
1e-8
-97.52
-1139.29
by domain experts as outlier-free, Then we can estimate the model hyper- parameters on that data, and fix them for the RRP on new data sets that do not know to be outlier-free.
D.6 Comparison Robust and Conjugate Gaussian Processes (RCGP)
We report extended empirical comparisons with Altamirano et al. 12's 's RCGP method, using their experimental setup and method iipllementtton GPFLow Including GPFlow our own bench- marking setup and compute resources proved difficult. To circumvent this, we wrote wrappers for both Torch's standard GP and RRP, which also accounts for any orthogonal implementation differences between the two frameworks, and ran the benchmarks locally
M-series MacBook.
See Tables 2 and B for the mean absolute error and negative log predictive density, respectively. The tables include the empirical mean and standard deviation over 20 replications corrupted version of the following base data sets: 1) Synthetic, which is generated as a draw of GP with a exponentiated quadratic kemel,, and four data sets available on the UCI machine learning repository [4]) in particular, 2) Boston (29, 3) Concrete [68], 4) Energy [63] and 5) Yacht (27]. The benchmark considers no corruptions (o Outliers"), "Asymmetric Outliers", which are uniform in are shifted negatively in y, "Uniform Outliers", which shift y in both directions (positively and negatively), and "Focused Outliers", which form concentrated clusters in both x and Any bold entry in the table signifies results that is within one standard- -error of the best result's one standard-errol confidence bound.
Table 2: Mean absolute error (MAE) using Altamirano et al. 2 experimental setup in GPFlow. RRP is always competitive with the other methods, and outperforms them significantly for uniform and asymmetric outliers.
Standard GP (CPPH)o)
Studens-t (GPFLow))
RCGP (GPH)ow)
Stundard (Belbech)
(BoToroh)
No Outliens
Synthetic
881e-2 (2.09e-)
(2.08e-3)
8,90e-2 (2.19e-3)
8.90l=2 (2.19e-3)
Boston
2:24e=1 (7. (7-1le-3)
2.24e-1 50e-3)
(7.29-6)))
2.08e-1 (5.53e-3)
2.08e-1 (5 (5.75e-3)
Concrese
2.00e 1 (3:72e 3)
1999-1 (3.73e
1.93e-1 (4.05e-3)
1.92e-1 (3.80e-3)
1.93e-1 (3.66e-3)
Energy
3.01e:2 (3.65e-4)
3.M4-22 (L.16e-3)
2.42e-2 (1.01e-3)
3.09e-2 (9.4le-4)
3,03e 2 (8 (8.88c
(2.09e3)
2.02e-2 (1.29e-3)
(4.6603)))
1.49e22 (1.95e-3)
1.42e22 (1.5He-3)
Unifoem Outliers
Synthetic
3:45e=1 (1. (1.48e-2)
2.93e (1.23e-2)
2.15e-1 (9.30e-)
3:48e-1 (1.62e:2)
8.99e-2 (2.37e-3)
Boston
4.97e1 (2 80e-2)
3.86e-1 (1 (1.52e-2)
4,94p-1 (1.81e-2)
6.81e=1 (2.25e-2)
2.83e-1 (3.16)-2)
Cinnrrtt
3.80e-1 (8 36c-3)
3.54G-1 98c-3)
4.40o1 (7.976-3)
1 1.80e-1 (7.77e-3)
2.03e-1 (3.98e-3)
Energy
2.73e=1 (2.40e-2)
(3.55e2)
5.600e11 (1.47e-2
3.06e=2
Yachi
3.23e-1 (1 (8e-2)
2.29-11 (1.54e-2)
4.18c-1 (3.4Be-2)
6.81e-1 (3.21e-2)
1.20e22 (1.04e-3)
Asymmetric Outlies
Synthetic
1.11es0 (2.24e2)
1.0eed (2.34e-2)
7.726-1 (3.91e-2)
1.7ea0 (2.55e-2)
2.71e-1 (9, (9.96e22
Boston
699e-1 (1 40e-))
5 300e-1 (1 (244-2
5.31e-1 (2.33e-2)
8.55p- (3 (3.04e-2
3,35e- (4.33e-2)
Conceate
6.60-1 (1-21o-2
5.550-1 (8o-2)
491c-1 (1.1le))
6.60e-1 (1.25e-2)
(4.45e-3)
Energy
603e-1 (1:1262)
4.74-1 (9.993-5)
(3.94e-2)
7.05e-1 (1.78e-2
3.07e-2 (8.33e-4)
Yachi
590e-1 (1 60e22)
4 58c-1 (9.31o-3)
4.326-1 (3.30e-2)
8.36e-1 (3.76e-2)
2,70e-2 (1 (1.35e22)
Fecssed Cuttiers
Syuthetic
1.89c-1 (1.50e-2)
192c-1 (1.56e-2)
Ló4e-1 (1.29e-2)
1880c-1 (1.46e-2)
2.00c-1 (1.62e22)
Boston
1.44e-1 (8.48e-3)
1.0e-1 (1.30e:2)
2.49e-1 (1.04e-2)
2.49e-1 (8.37e-3)
2.48e-1 (8.38e-3)
Concrcte
2.40e-1 (5 546-3)
2.35--1 (5.63623)
2.386-1 (1.10-2)
2.24c-1 (6.20c-3)
2.25e-1 (5.M6e-3)
Energy
8.33e22 15 46e-2)
3.24e22 (1.36e()
2.92e-2 (8.95c4)
3.03e-2 (8.95e4)
3.03e-2 (8.95e44
Yachi
2.43o-1 (6 (9)-))
1.106-1 (2.0762)
1.9362 (1.86-3)
1.76-2 (2.96e3)
1..06-1(2.20v...
D.7 Comparison Robust Gaussian Process with Huber Likelihood
In the following, we compare our method with additional variational GP baselines with Laplace and Huber likelihoods, and translated the Matlab code of the "projection statisties" Algikar and Mili
I to Porch.. We then combined the projection-stattstics
weighting of the Huber loss with
25

---

Table 3: Negative log predictive density (NLPD) using Altamirano et al. 2 experimental setup in GPFlow. RRP is generally competitive, and outperforms other methods significantly for uniform and asymmetric outliers.
Standard (IP (CiFFlow)
Student-d (GPFL.ow)
RCGP (CIPFlow)
Standard (ReToech)
RRP (RuTorch)
Outleers
Synthetic
.821e.1 (2.16e-2)
-8.00-11 (2.16e-2)
-8.33e11 (2.21e-2)
-805c-1 (2.25e-2)
(2.24e-2)
Buston
2,08c 1 (4.26c 2)
1.99c-1 (3.76e 2)
1.95e-1 (8.69e 2)
9.24c-2 (2.82e-2)
9.43e-2 (2.82e-2)
Concrete
1,48c-1 (2.33e-2)
1.29c-1 (2.31e-2)
1,09c-1 (3.19e-3)
1.11c-1 (2.73e-2)
L.6ee1 (2.68le2)
Eieeryy
-1.22eel (3.83e22
(3.46e-2)
1.96ee0 (4.37e-2)
-1,77ee0 (3.08c-2
1.09e00 (4.0022))
-1.79pe40 (3.23c-1)
-2.05e-00 (6.61e-2)
00le00 (3.27e-1)
-2.47ea0 (1.06e-1)
22.23es0 (1.846-1)
Unifomm Outliess
Synthetic
(1.41e0
1.4ree0 (1.55e
L (1.39e22)
1.58e0 (1.40e22
7.99t.1 (2.26e-2)
Baston
(7.5-0 3)
1.5cce0 (8.90 3..56e33 3)
(1.2-001.22 2)
181c=0 (5.95e 2)
4.52e-1 (1.34e-1)
Concrete
1.64p00 (6.83--3 83c-3)
1,51e40 36e-3)
1.6e400 (7.57e-)
1.63e00 (6.72e-3)
1.47e-1 (2.51e-2)
Energy
166ee00 (2.75e-3)
1.30ce00 (1.76e-2)
1.1e=0 (9.74e-2)
1.70e=0 (8.58e-3)
(3,72e-2)
Yarha
1.61e=0 (9.50e 3)
1.46c+0 (1.19e 2)
1,78c=0 (4.57c
-2.32e00 (1.97e-1)
Asy mmmerri Outliers
Synthetic
1,94e-0 (9. 39c-5)
190c40 (9 96c-3)
1.89--0 (1:13-22)
1.93e=0 (1 (066-2)
-.70e-1 (2.30e-1)
Baston
166e=0 (1.01e
1,56c00 (0le 2)
1.62=00 (9.51e:3)
2.52e=0 (2.26e-)
6.1ee-1 (1.19e-1)
(7,43e:3)
(3e5e33
162e-0
1.65=00 (7:55e3)
(2.86622)
Energy
1.62c=0 (8, 33
1.49ce0 (9.70e-3
1.62+00 (1.82e:2)
1.71e=0 (9.79-33
1,0ee00 (3,55e-2)
Yachu
1.59e=0 (9.35e-3)
1.45e+0 (1.01e-2)
1.5ee+0 (8.15e-3)
1.75e=0 (1.41e
2.17e-1 (1.82++0)
Focused Ouiliers
Synthetic
6.78e-1 (6.59e-2)
6,11e-1 (8-68e22)
5.69e-1 (4.07ee2)
6,66e11 (6.60e-2)
676740 (2.24cx0)
Baston
3.21e-1 (6.40e-2)
(3.04e-2)
(3.06e-2)
Cancrete
2.77e-1 (2. (2.5e-2)
2.72e=1 (2.81e-2)
2.43e-1 (4.36e-2)
2.22e-1 (2.87e-2)
2.28e-1 (2.73e-2)
Energy
2.81044 (2. (1e444)
-1.67-40 (4.13e-2)
-1.7ps40 (4. (4.996-2)
-1.70e40 (4.59e-2)
-1 -1,70c40 (4.59--2)
Yacht
9.38c+1 (5,09ce41
-4,60e-1 (2.91e-1)
2.47740 (8.09e22
-263e+0 (5.55e-2)
-2.1ee00 (S.66e-2)
variational (referred to Huber-Projections to get as close as possible to a direct comparison to Algikar and Mili 1 without access to a Matlab license.
Tables and 5 shows the root mean square error and negative log predictive density on the Neal, Friedman 5 and Friedman 10 test functions, as well as the Yacht Hydrodynamics [27] and California Housing (51) datasets from the UCI database (34), where 15% of the training data sets of the models were corrupted. Tables 6 and below were generated in similar way, but 100% of the data were subject to heavier-tailed noise, either Student-t or Laplace.
In summary, the Relevance Pursuit model generally outperforms the variational GPs with heavy-tailed likelihoods when the corruptions are a sparse subset of all observations. Unsurprisingly, the GPs with heavy-tailed likelihoods perform best when all observations are subject heavy-tailed noise. While such uniformly heavy-tailed noise does exist in practice, we stress that this is a distinct setting to the common setting where datasets contain subset of a-priori unknown outliers, while dominant fraction of the data can be considered inliers that, once they are identified, can be used to train model without additional treatment.
26

---

Table 4: Comparison with Huber GP: Root mean square error with 15% Corruptions
Data
Standard
Relevance Pursut
Student-t 1
Laplace
Huber
Huber + Projection
Neal
Uniform
3.87e=1 (2.40e-2)
3.79e-2 (1.10e2)
4.18e-1 (8,55-33
4.18e1 (d.72e-3)
8. 1 (1.76e: 3)
4. B.1B (d.76e-3)
Constant
1.1 18e-1)
4.37--2 (1.87e-2)
7:260-1 (1.38g-1)
4.75e=1 (1 (1.28e-2)
4.89c-1 (662-))
4.89e-1 (1 (2)-))
Studeas+
2.74e-1 (2.70e-1) 1)
6.71e-2 (2.29e-22
3.90e-1 (3.59e-3)
4.07e-1 (0.79e-3)
4.0% 1 (5.00 (5.02e 3)
4,09c-1 (5.02z-3)
Laphace
7.19e-1 (6.17e-2)
6,27e-2 (1.50e-2)
4,01e=1 (4,42e-3)
4,20e- (4,34e-3) 3)
4,20g-1 (1.00e-3)
4,20e- (4,00e-3)
Friedman 5
Uniform
20 (8. Ste-2)
9,40c-2 (1.23e-2)
7.01e-1 (1.57e-3)
7.30c-1 (7.35e-2)
6. 6.78e-1 (6 (6.67e-2)
2.12 (7.50e-3)
Cunstant
1.8c=1
7.29 (5.40-3)
8.98c-1 (1.40e-1)
2.37e=0 19e-2)
2.50==0 (B-2))
2.50c-0 (7.1ee2)
Suudent-
7.9 (8, 90c-1)
1.17e-1 (2.34e-2)
5.26e=1 (1 (1.26e-1)
1.27es0 (1499-7)
1.36-00 (9 34e-2)
1.31==0 (9.73e-2)
Taplace
7,59c-40 (4 (8--1)
9.63e-1 (1.15e-2)
6.62e-1 (1.34c-1)
1.51es0 (7 (7.78c-2)
1. 1.63e00 (6 (90--))
1.68e-0 (6.80e-2)
Priedman 10
Uniform
Laile+0 (8.94t 2)
4.84e (6.21e:3)
5,92e-1 (1.36e- 11
6.82c- 1 (1.03e-11
6,74g 1 (1.0% 1)
5.93e- 1 (9.43c 2)
Constaut
1.34c+1 (2.97e-1)
5.13a-1 (2,69c-3)
8.65e-1 (1.39e-1)
2.10cs0 (0e-2)
2.18es0 (3.01e-3)
2.18e-0 (5.01e22
Studear-t
7,93e40 (1/75ee0)
8,45e-2 (2.43e-2)
5.30e-1 (1.30e-11
1.40ca0 (9.10e-2)
1.40l-00 (9. 122-2)
1.40c=0 (9.44e-2)
I Taplsce
6.42e40 (3 (3-42e-1)
5.70x-2 (8.16e-3)
7.25e-1 (1.30e-1)
1.64c40 (4.98 98e-2)
1. 1,72e=0 (1 (66e-))
1,72e-0 (1 (.6e-2)
Yachi
Unaform
6,72c40 (5 (5.08e-1)
8.56e-1 (4.99)-22
1.62e1 (1.81e-1)
1.2ce1 (2.05e-1)
L-lee1 (2.04c-1)
LGle=1 (1.9e-13
Condint
3,50e=1 (1.51e+0)
8.27e-1 (5.78e-2)
1.0kll (2.11e:1)
1.59c]] (2.02e: 1)
1.58e=1 (1.73c
1.58=11 (1 (1.95e:1)
Studear-t
2.1lea1 (3.53eat)
1.13c+00 (2.03e-1)
1.6ke-1 (2.176-1)
1.67ss] (2.15e-1)
1.6ee41 (2.08e-1)
1.6e==1 (2.02e-1)
Laplace
1.lcel (1.97ee0)
1.15e+0 (1.54e-1)
1.63e11 (2.10e-1)
1.68ce] (4.00e-1)
1.63e1 (2.05e-1)
L.60e11 (1 (1.99e-1)
CA Housing
Uniform
7.10e=1 (7.58 (7.58e-3)
7.59e-1 (1.75e:2)
LI6c+0 (1.79- (1799-3)
7.7e++ (3
L.77-0 (2.93e-3) 3)
LIBc-O (6. (7:17)3)
Constant
2,28e=0 (5: 122-2)
6.35e-1 (4.38--3
1,170=0 (2.30e-3)
1.16e=0 (1 (7e-3)
1.1e-00 (1 (8e-3)
1.17e-0 (5.5le-3)
Student
1,34++0 (1 (1.88-
6.56e-1 (5.46e:3
L.8c+0 (3.56e
1.19e00 (9.13e-3)
1:18e00 3.833:33
1.19c-0 (6.24e
Laplace
1,00++0 (5.00005..0le)) (5.0le-2)
6.5le-1 (4.71e-3)
11
1.1ee=0 (3.91e:3)
1. 1.18e=0 (3:57e-3)
1.18--0 (5.30e-3)
Table 5: Comparison with Huber GP: Negative log predictive density with 15% Corruptions
Data
Standand
Relevance Puruit
Student-t
Laphace
Huber
Huber Projection
Neal
Unifoem
6.87es0 (5ead))
3,80c+0 (4.01cst))
185e+0 (9.93e-2)
1644s00 (1.20e-1)
1.70e00 (1.23e-1)
1.70cs0 (1.23e-1)
Constant
1.32e22 (1 29)
-2,00e+0 (2.90 (2.99c-1)
2,03e+0 (1 1..44e-1 1)
8.67e-1 (4.32e 2)
8.84e:1 (3.57e:2)
8.84e-1 I (3.87e-2)
Student-
1.98c40 (9
1.99--1 (2.29cxl)
1,7e=0 (1.17c-1)
1,25e40 (8 (8.32e-2)
1.25c40 (K0836-2
1.28440 (8013e22
Laplace
3.07e92 (5ee2)
2.02e10 (2.83est))
1.73cx0
1.20ev0 (8.59c-2)
1,21c10 (8.99e-2)
1.21es0 (8.99e-2)
Friedman
Unifsem
2.0e+00 (1.88e 2)
-1.35e00 (7.222-2)
(1 (1.0le-1) I)
8:29e (9. 19.ile22
8.88e (1.08e-]
Constant
4:28eell (2.25e-2)
-1.13c+0 (1.82e:2)
8.24e-1 (3.37e-1)
2.28e10 (2 (2.63e-2)
2.54ea0 (3. 33113-2)
2:34c00 (3.13e-2)
Student é
3.42e40 1)
-1.11e00 (1.M4-1)
1.78e 2 (2,86e 1)
1.55e+0 (1 1.44 I)
1.00e00 (1.1le 1)
1. 1.59++0 (1.10e 1)
Laplace
3,55e+0 (7.84e-2)
-1,22e+0 (5,89e-2)
3:47e=1 (3.13-11
1.75e+0 (8:46e 2)
1,90+00 (7.36e-2)
1 1.90++0 (7.36e
Friedmam 10
Unifoem
1.36e00 (3.04e 2)
-182e=0 (3.97e:2)
6.78e-2 (3.85e
9:20e=1 (1.99c-1
9.39e 12.01e:1)
7.65e: (1.80e-
Constant
4.21e+0(1.3 (1.39e-2)
-1,27e=0 (1. (1.10e:2)
8.85e-1 (3.79e-1) 1)
2.20e+0 () (5e-)))
2:27e+0 (1.51e-2)
2:23e+0 (1. (1.5de-2)
Student-
3.34400 (1.24g-1) 1)
-1.40c-0 (1.17e-1)
-6 80e-2 (3 13.65e-1 1)
1,96++0 (1 15---1)
1.95e40 11236-11 I
1 1664+0 (1 (1.22e-1)
Laplace
3.46e+0 (6.68e-2)
-1.44c-0 (5.91e-2)
4.89c- 1 (3:94c i)
2.13er0 (6:35e:2) 2)
2:21cH0 (2.01e:2)
221e+0 (2.01e 2)
Yacht
Uniform
345e+0 142--1)
2.37e+4 (2.79e-1)
L (1 ()
7.07e+1 (5 17040)
7.6ca1 38e+0)
500+2
Constant
70ce0t (1.70u 70e-1)
1.7eal0 (2.436-1)
7.16a=1 (7.68e40)
2.53e+1 (1 (.0eel))
1.95el 12640)
80ell (3 (3.46e-0)
Student-t
4 4.46440 (3.0% 3.08k-1) 1)
2.42+00 (3.38e1
1.32++2 (1 (1.05e+1)
8.20c=1 (5 (5.07ee+t)
8.13e+l (5 17e40)
1.55e+2 (1. (1.38e+1)
Laplace
4 35e+0 (2.88-1)
2.27e40 (377ee1
1. 1444-2 (9. (9.86e+0)
6.85c+1 (4 (4.720+0)
6. 6.12c+ (3.87040)
1. 1.17042 (8 (8.89)40)
CA Housing
Unitoem
3.14es0 (1.44c-1)
2.92es0 (2.06-11)
5.85e11 (9.5%--1)
6.77es1 (2.4lesl)
6.00es1 94eell)
854es1 (2.45e41)
Constant
3.57e40 (2.26e-1) I)
4.50e+0 (2. (e-1)
4.2o=1 (1 (1.50++0)
1.83+l 58c-1)
1.750+] (7003-1)
2.18 (A. 19c40)
Studen--
1.77e40 (6.656-2)
3.15call (1.676-1
5 5.95=4 (1 (1.44)0)
5.20c=1 (2. 19c4tl)
4.95+4 (1 (1.75e4t)
6:00e41 (2. (2.02e+1)
Laplace
161cs0 (4.21e-2)
3.5test (1.90)(111. (1.93e-1)
5.600e1 (1.58ex0)
4,23ce1 (1.72e40)
4,06es (1.47esl)
5.62es1 (1. (1.66e+1)
Table 6: Comparison with Huber GP: Root mean squared error for 100% Laplace noise
Data
Standard
Relevance Pursui
Studens it
Laplace
Huber
Huber + Projection
Neal
Student-t
1.96<=0 (3.64e-1)
I. 1.56+0 (2:48e- I)
7. 7.39e-1 (4.06e-2)
7,54e-1 (3.17e-2)
7.48e-1 (3.00e-2)
7.48e-1 (3.01e-2)
Laplace
1.51ce0 (1.06e 1)
2,40c+0 (2.72e
1.1er00 (9.18e 2)
1.13e+0 (9,90e22)
1.12c+0 (9.50e-2)
1.12++0 (9.50e-2)
Friedman
Studem-t
1.67e=1 (1 83e+0)
9, 0,30c+0 (3. ()055-1) 1.37e=1 (6.36e D)
5,26e+0 (1 (1.30e-1)
4.57e+0 (1.60e-1)
4.61e=0 (1 (1.55e-I)
4.61e+0 (1.55e-1)
Laplace
1.39ce1 (5.95e D)
8.33e+0 (3,69 I)
7,34e+0 (3.50e-1)
7.40c=0 (3.45e-1)
7.40e+0 (3-45e-1)
Friedman
Studeas-t
2.10e=1 (3)42ee0)
7,80e+0 (2X0e-1) 1.27er1 (3.99 L)
4.72c+0 (1.34e-D
4.02e+0 (9.57e:2)
4.05e=0 (9.34e-2)
4.05c+0 (9.34:22)
Laplace
1.30e=1 (3.77e-1)
7.24e+0 (1.99e L)
6.07e+0 (1.91e-1)
6,26e=0 (1.71e-1)
6.26e=0 (1.71e-1)
Yachi
Studeas-t
5.89c11 (1.57ex1)
2,44ee1 (1.77e+0)
1.57e11 (1 (1-36e11)
1.57ee1 (1.5te-1)
1.57ee1 (1.51e:1)
1.60e= (3.16e-1)
Laplace
2.57e=1 (1.17c+0)
4.75c+1 (3.51e+0)
1.64+11 (3.18e11
1.1ee1 (2.53e-1)
1.62e1 (2.60e-1)
1.67e1 (3.55e-1)
CA Housing
Student-t
2.92ce0 (5 30e-1)
1.34e+0 (6.17e-2)
1.17++0 (Ad55e-3)
1.18e+0 (4.42e-3)
1.17ev0 (4.14e-3
1
Laplace
1,57c+0 (8:55e2)
2.05e+00 (1:44e I)
1.21e+0 (1.08e-2)
1.20++0 (9.39--3)
1.20=00 (9,69-33
1.23e+0 (2.05e:2)
27

---

Table 7: Comparison with Huber GP: Negative log predictive density with 100% Laplace noise
Standand
Suuden-F
Hiher
Haber
Neal
Studem-&
1955a00 (1.22c-1)
(2.00e11)
(6.08e-2)
(5.71e-2)
(5.15e-2)
Laplace
1.89e00 (1.14p-1) 1)
1.38e22 13 13991e11)
1.57es0 (9.22e-2)
1.55er0 (9.20e:2)
1.55er0 (9.12e-2)
1,55er0 (9.12e:2)
Friedman
Studen-t
4.43es0
196-2)
3.93ev0 (2.19
3.11es0 (2.6992))
2.96c (3.49e-2)
2.97es0 (3.27e22
2.77e40 (327e-2)
Laplace
4.55e40 (1811-))
5.56=+0 886-2)
3.646+0 (4.03-2)
3,46e40 (4.45e22)
3,45c+0 (4.31-2)
3.48c40 (4.3e-2)
Friedman 10
Studem-t
4.00ss0 (8.54p-2)
5.89 (1.36e2)
5.04e10 (2.45e-2)
2.83esti (2.30e-2)
2.84c (2 (2.22e-2)
2.84e10 (2.22e-2)
(1.25o-2)
(1.276-2)
3.0000 (2:53e22)
(2.61o22)
Student-é
5.14ee0 (1 85e-1) 1)
4.60le00 (8.15e-2)
(.30e+0(1.40a
6.79ed (1.32e-1)
6.70es0 (1 (.322-1)
1.28ee1 17.75e-1)
Eaplace
4.74ex0 (4.45o-2)
5:30s+0 (8.31e-2)
4.91e40 (8.87-2)
5.10040 (1.176-1)
5.12040 (1 (46-1)
9,50600 44e-1) (444411)
Housing
Student-E
2.0000 Ole-1)
1.70ee0 (6.41le-2)
7.4le+0 1)
50e40 (2 54e-1)
6:41ee0 07e-1)
6.8e+00 (3 28e-1)
Laglace
1.88e00 (5 (5.25e-2)
2.24e=0 (9 (2.41e2)
2.97e40 (6.01e2)
3.53ea0 (7 (7777-2)
3.27e40 846-2)
6.64e40 (3.28e40)
28

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Everything that is claimed in the abstract and introduction discussed in detail in the rest of the paper.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Limitations of the work are mentioned throughout the paper as well as summa- rized the conclusion.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiitiinn asymptotic approximations only holding locally). ). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
29

---

Answer: [Yes]
Justification: All theoretical results clearly state the required assumptions; proofs are provided in Appendix C
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Code for our method as well for the benchmarking setup is provided in the supplementary material, and core models and algorithms will be open-sourced upon publication.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access a hosted model (e.g.. in the of large language model), releasing of a model checkpoint, or other means that are appropriate the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
30

---

reproducible, they 5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Question:
Answer: [Yes]
Justification: Code for reproducing the results the paper is included in the supplementary material. All data sets are publicly available and referenced in the work and the code.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.co
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/nips. wwtpcc/pccpccpcccceoeeoooop httpa//gggudeee/deuuuuubbiiiiiii
wwwsssssssiieeeeeeeoo for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Details on the experiments are provided in the experimental section as well as the Appendix. Exact implementation details are contained in the included code submission.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Variance in the results is provided throughout in the form of confidence intervals and/or violin plots.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emoor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
31

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Total compute resources spent are described in Section D.4 and estimated at around 2 CPU years in total.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics http:////eeerrpsssccplieileeeeee inneeiieeeeeiccsccdieeiieeee
Answer: [Yes]
Justification:
The research fully conforms to the NeurIPS Code of ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there a special consid- eration due to laws regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
The paper makes foundational methodological contributions to robust prob- abilistic regression, and we do not see any direct socictal impacts of the work that would require specific discussion.
Justification:
Guidelines:
32

---

The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negatíve societal impacts include potential malicious or unintended uses (e.g.. disinformation, generating fake profiles, surveillance), fairness considerations (c.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
All creators original owners of the code and data sets that were used in this paper are properly cited and credited. Licenses for these assets are mentioned in Sec. D and our use complies with those licenses.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
33

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paperssi thcode. cowc datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is available online, the authors encouraged to reach out the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: The only asset introduced in this paper is the code for method and benchmarks; it does not introduce any other assets (data sets or models). The code is provided anonymized zip file for the reviewers. Upon publication, the code released as open source and documented according to community standards.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
34

---

The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
35


Reviews:
[{"review_id": "QHl5pJxnNw", "reviewer": "Reviewer_CPiv", "summary": "The authors propose a robust Gaussian process model in the sparse heteroscedastic noise setting. Orthogonal Matching Pursuit- like algorithm is proposed to infer data point specific noise levels. The negative log marginal likelihood function to be optimized is claimed to be strongly convex by reparametrizing the additive noise levels. The method is validated on synthetic as well as real-world experiments.", "questions": "1. Why Matern kernel was chosen in one of the baselines?\n2. What is the difference between RGP-RP*-BW and RGP-RP-BW. I see that the former is canonical one. What is the difference between canonical and non-canonical reparameterization of $\\mathbf{s}$?\n3. What is $\\mathbf{K}_{0}$ on line 218?", "limitations": "Limitations are discussed as I hoped they would be.", "rating": 7, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. The authors formulated a convex $\\mathcal{L}$ in the Gaussian process regression model by reparametrizing the noise levels, which is very desirable for optimization perspective. The representation of the paper is good, however, I found the section 4.3 difficult to follow.\n2. The relevance pursuit is a novel approach to the heteroscedastic noise case of robust GP regression.", "weaknesses": "1. Much more work on heavy tailed outliers is done. So, more references and discussion need to be added in related work. Kersting et. al. 2007 proposed GPs for heteroscedastic noise setting. They also learned the homoscedastic part of the variance aiming to learn the remaining part using a second GP. The work is related and at least needs a mention. It would have been nice to see comparison with this model. \n\n2. While this work employs greedy algorithms for GPs, it is not the first to do so, so I would not consider it exceptionally novel.\n\n3. The experimentation is adequate but less discussed.\n \nRef: Most Likely Heteroscedastic Gaussian Process Regression\n\n\nMinor corrections: \n\n1. Consider adding references for these three categories on line 34. No reference is given for down weighting procedures even in related work.\n2. Figure 1 and 2 are not discussed at all. What regression example? \n3. broken line 150 \n4. unclosed braces figure 3\n5. no need of repeating the abbreviation RIP on line 200\n6. Broken line 264\n7. Is the x axis of right figure in Figure 4 percentage of corruptions or probability?\n8. Different notations need to be used for $\\delta$ on lines 67 and 112.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> Kersting et. al. 2007:\n\nThis is indeed related work, and we will discuss it in the updated manuscript. The main difference to our method is that their approach is unlikely to work well with outliers. In particular, while they do account for heteroskedastic observation noise, they model the logarithm of the standard deviation of the noise noise process as another GP. This not just restricts the observation to be (heteroskedastic) Gaussian, but, importantly, imposes a smoothness assumption on the variance of the observations as a function of the inputs. In contrast, our approach infers individual additional noise terms, and thus allows us to go beyond the setting in Kersting et. al. 2007 to handle even gross outliers (which require a noise process that is discontinuous in the inputs). Another downside of the approach in Kersting et. al. 2007 is that they require an EM algorithm to estimate the parameters of the two GP models, which does not have any theoretical guarantees and in practice often has convergence issues (see e.g. discussion here).\n\n> “I found the section 4.3 difficult to follow”\n\nThank you for your feedback. We will improve the clarity of our writing in this section in particular by introducing the general framework of Bayesian model selection before stating its instantiation in our work and providing further background references. Please let us know if there are additional edits you would like us to do.\n\n> “While this work employs greedy algorithms for GPs, it is not the first to do so”\n\nWhile greedy algorithms have been applied for many problems in general (outside of GPs), their application to the sparse optimization of noise variances for robustness, which we show permits strong theoretical guarantees, has to our knowledge not been studied and proposed before. Please let us know if you would like us to add any additional reference for related work.\n\n> “The experimentation is adequate but less discussed”\n\nWe produced a significant number of additional results with new baselines and test problems for this rebuttal, and will improve discussion of the entirety of the empirical results (see general response).\n\n> “Figure 1 and 2 are not discussed at all. What regression example?”\n\nThank you for catching this. The illustrative regression example is based on a synthetic one-dimensional modified sine function with a few outliers. We will detail the exact setup for this in the appendix.\n\n> “Is the x axis of right figure in Figure 4 percentage of corruptions or probability?”\n\nYes, this is stated in the figure’s caption. We will add an axis label as well in order to improve clarity. \n\n> \"Why Matern kernel was chosen in one of the baselines?”\n\nNo particular reason other than that the Matern-5/2 kernel is widely used in practice and the default in many GP / BO libraries (e.g. botorch, trieste). We expect the results to be quite similar if another common kernel, say an RBF kernel, were used.\n \n> “What is the difference between RGP-RP*-BW and RGP-RP-BW. I see that the former is canonical one. What is the difference between canonical and non-canonical reparameterization of s?”\n\nThis difference is described at the beginning of Section 5.2. The “canonical” parameterization is that which directly uses the robust variances $\\boldsymbol \\rho$ in the parameterization. The non-canonical parameterization is the convex parameterization we introduce in Section 5.2, i.e., that where ${\\boldsymbol \\rho}(\\mathbf s) = \\text{diag}({\\bf K_0}) \\odot ((1 - {\\mathbf s})^{-1} - 1)$ and the new parameter is $\\mathbf{s}$. The method label coding is as follows: A “*” indicates using the canonical (not necessarily convex) parameterization, and “BW” indicates using the backward algorithm for relevance pursuit (Algorithm 2 in Appendix A). We will improve the clarity of the exposition here and clearly explain the method labels. The convex parameterization is a key ingredient that allows us to prove strong theoretical guarantees for the result of the greedy algorithm, and also has beneficial effects on the convergence of the hyper-parameter optimization algorithms, as we demonstrate in the other responses.\n\n> “What is $\\bf K_0$ on line 218?”\n\nThanks for catching this forward reference, $\\mathbf{K}_0$ is defined in Theorem 8: $\\mathbf{K}_0 = k(\\mathbf{X}, \\mathbf{X}) + \\sigma^2 \\bf I$. We’ll move the definition up in the text.", "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": "", "extra_scores": {}}, {"review_id": "0Dns4c3xdN", "reviewer": "Reviewer_aaaZ", "summary": "The paper proposes a robust Gaussian Process regression by inferring data-point-specific noise levels with a sequential selection procedure maximising the log marginal likelihood. The authors show the good performance of their method in a mix of synthetic and real-world regression tasks, including Bayesian optimisation.", "questions": "- What's the difference between RGP-RP and RGP-RP* in the experiment section?\n- What's the schedule $\\mathcal{K}$ selected for the experiments? How does this affect the performance and fitting time of the method?\n- While the Friedman and Hartmann6 test functions are well-known, the authors should provide a reference for these functions.\n- Why use Hartmann6 as a test function for regression? Hartmann6 is typically used to test optimisation due to its properties. Why not use standard UCI benchmarks?\n- Do the authors know why the student-t GP performs so poorly in the Hartmann6 function, even when the noise comes from a student-t distribution?", "limitations": "No obvious negative societal impact", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-written and easy to follow. The method is clearly motivated and presented, with a good amount of related work discussed. The authors address a problem—robust regression in the context of Gaussian processes—of sufficiently broad interest to merit publication. The authors also provide a theoretical analysis to support their method.", "weaknesses": "In my opinion, the weakest part of the paper is the regression problems in the experiment section. It only benchmarks against two test functions, one of which is not a standard benchmark for regression; I would like to see more standard UCI benchmarks. Additionally, the proposed method does not clearly outperform other methods, particularly in the Hartmann6 test, where adaptive Winsorization outperforms RGP-RP. The authors do not discuss the results of the regression problems, so a discussion providing guidelines on the cases where the proposed method will perform better is needed.\n\nFinally, the authors claim their method \"permits fast and robust inference\"; however, in Appendix D.1, the method is slower than the Student-t GP for the Hartmann6. It would be useful if the authors discussed this in the main paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> “The weakest part of the paper is the regression problems in the experiment section”\n\nAs part of our rebuttal, we have produced a significant number of additional results on UCI benchmarks, the benchmarks and real-world Twitter crash example of Altamirano et al. 2024, and included an additional baseline (see general response and attached pdf).\n\n> “Additionally, the proposed method does not clearly outperform other methods, particularly in the Hartmann6 test, where adaptive Winsorization outperforms RGP-RP.”\n\nWe don’t necessarily expect our method to outperform all other baselines in all settings. But with our additional empirical results in our general response, we do show that it (i) works much better in situations in which outliers cannot be easily separated from the rest of the data, (ii) is competitive in most other settings. We also note that while the adaptive winsorization heuristic indeed performs slightly better than RGP-RP on Hartmann6 with constant, well-separated outliers, it completely fails on the Friedmann10 function and is barely better than a non-robust GP. In general, winsorization does not perform well if the outliers are not clearly and a-priori separated from the rest of the data.\n\n> “The authors do not discuss the results of the regression problems, so a discussion providing guidelines on the cases where the proposed method will perform better is needed”\n\nWe agree that the discussion can be improved and we will include this into the paper (see general response for discussion).\n\n> “the authors claim their method \"permits fast and robust inference\"; however, in Appendix D.1, the method is slower than the Student-t GP for the Hartmann6. It would be useful if the authors discussed this in the main paper”\n\nThe Student-t likelihood performs very poorly on Hartmann6 and has very high variance. We suspect that there are numerical issues in the optimization of the variational objective that results in the fitting terminating early and at suboptimal points, thus producing faster fit times. Note that the Student-t log likelihood is also non-convex, which can give rise to highly suboptimal local minima. We plan to investigate this in more detail.\n\n> “What's the difference between RGP-RP and RGP-RP* in the experiment section?”\n\nSee [our response to reviewer CPiv](https://openreview.net/forum?id=5FATPIlWUJ&noteId=JGfjtKQW5T).\n\n> “What's the schedule K selected for the experiments? How does this affect the performance and fitting time of the method?”\n\nThe default schedule used throughout the experiments is `[0%, 5%, 10%, 15%, 20%, 30%, 40%, 50%, 75%, 100%]` of the number of observations (traversed backwards by the backward algorithm). The fitting time is linear in the number of steps in the schedule, as a new model fit is carried out for each step of the schedule. \n\nThe advantage of a schedule with a finer granularity is that the model could be more data efficient, by introducing additional noise variances for outlying data points only, so that all other data is \"trusted\" up to the homoskedastic noise variance $\\sigma^2$. For most practical applications, it is probably sufficient to test for single-digit percentages of outliers, the default schedule here is designed to exhibit high generality, as is evidenced by RP's performance on the variety of benchmarks, including its outperformance on the RCGP benchmarks and Twitter crash example of Altamirano et al. 2024 included in the general response, using the same schedule.\n\nIf timing was a particular concern, the optimization could likely be accelerated substantially by keeping track of the approximation of the Hessian matrix generated by each L-BFGS call, which could be used to warm-start the optimization of each step and lead to significantly accelerated convergence, as the consecutive optimization problems tend to be similar. \n\n> “While the Friedman and Hartmann6 test functions are well-known, the authors should provide a reference for these functions.”\n\nThanks for the callout, we will provide the respective references.\n\n> “Why use Hartmann6 as a test function for regression? Hartmann6 is typically used to test optimisation due to its properties. Why not use standard UCI benchmarks?”\n\nWe have produced additional results on a number of UCI benchmarks as part of this rebuttal (see general response). \n\n> “Do the authors know why the student-t GP performs so poorly in the Hartmann6 function, even when the noise comes from a student-t distribution?”\n\nWe believe you are referring to the right column of Figure 4, whose x-axis corresponds to the fraction of data points that are corrupted by the outlier process. First, the axis label is missing in the submission, and we will correct this. \n\nWhat the figure shows is that the Student-t GP only starts to dominate the other methods once more than 40%-50% of data points are corrupted by Student-t noise. We indeed expect the Student-t GP to be the best model as the fraction of Student-t-distributed noise increases, as it becomes the correct noise assumption when 100% of data points are “corrupted” by Student-t noise. However, the results also show: if the main goal is to make a GP robust to a _sparse_ set of outliers – rather than a uniformly heavy tailed observation noise – then RP significantly outperforms the Student-t GP by leveraging this sparse structure.", "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": "", "extra_scores": {}}, {"review_id": "jx76G9U8VS", "reviewer": "Reviewer_geQY", "summary": "This work proposes a new way to perform heteroskedastic regression using Gaussian processes via data-point-specific noise levels. These noise levels are inferred using a sequential selection procedure maximizing the log-marginal likelihood. The authors show that under a specific parametrization the log-marginal likelihood is strongly concave in the noise variances and give approximation guarantees for their proposed algorithm. Experiments on regression and Bayesian optimization tasks demonstrate the benefits of this approach.", "questions": "- Is your assumption about the data corruption different from Huber's $\\varepsilon$-contamination model (see e.g. https://arxiv.org/pdf/1511.04144) or is it identical? Equation 2 suggests they are different, but the text and experiments often refer to only a fraction of the data being corrupted (e.g. in Section 4.3 and Section 6).\n- How does your approach compare to a vanilla GP on a dataset with no corruption? Does it \"fail gracefully\"?\n- Do you see better performance when *jointly* optimizing the hyperparameters for the convex reparametrization of the data-point-specific noise variances than if not reparametrizing?\n- How did you set the schedule to test data points for being outliers in the experiments (c.f. lines 182-184)?", "limitations": "I would have liked to see a bit more discussion on the limitations of the approach in Section 7. For example, the fact that the theoretical result does not apply when jointly optimizing the hyperparameters as is the case in the experiments, and a discussion of the wall-clock runtime.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Making Gaussian processes robust to some degree of model misspecification is a worthwhile goal and this paper provides a nice addition to the existing toolbox of methods. The paper is well-written and presents largely well-executed experiments for regression and Bayesian optimization. In particular, in tasks where Bayesian optimization is applied, it seems plausible that outliers appear in the way they are synthetically generated in the evaluation. I also appreciated that the authors submitted their code.", "weaknesses": "### Theory\nWhile correct as far as I can tell, I question the value of the theoretical result, if in practice you optimize the hyperparameters of the covariance function and $\\rho$ jointly. Is there reason to believe, that the convexity in $\\rho$ for fixed hyperparameters is beneficial given this choice? It would be informative to see an experiment that compares the non-convex parametrization to the convex one where all hyperparameters are optimized jointly. If the convex parametrization improves convergence, this would suggest that convexity is beneficial even when optimizing all hyperparameters jointly.\n\n### Related Work\nEven though this is quite recent work, I would have liked to see a discussion of Altamirano et al. (2024) in the related work section. They introduce data-point-specific observation noise based on a generalized Bayesian inference perspective. I think a comparison to their approach in the experiments would improve the paper, but I would not expect this given the recency of the work.\n\n- Altamirano, Matias, Francois-Xavier Briol, and Jeremias Knoblauch. \"Robust and Conjugate Gaussian Process Regression.\" International Conference on Machine Learning (ICML), 2024. URL: https://arxiv.org/abs/2311.00463\n\n### Experiments\nAs the authors laudably acknowledge, the fit times were 100x slower (see Figure 8) than the simple baseline approaches. How does that impact the Bayesian optimization results if plotted as a function of wall-clock time? It seems that eventually some of the other approaches catch up to your method. If the fit time is 100x that of the baseline approaches, the reward curves might look different if plotted as a function of wall-clock time. Or is there an implicit assumption that each evaluation is significantly more expensive than fitting the GP?\nThe experiments on regression problems could be expanded in my opinion, e.g. with benchmark datasets from the UCI repository. The test problems that were considered are arguably designed to test optimization algorithms (hence their usage in BO).\nFrom looking at the code the \"real-world problems\" in Section 6 have synthetically generated corruptions. I think the paper could be made stronger if there were at least one experiment where the corruption is not synthetically generated, but simply part of the observation process. \n\n### Other Weaknesses and Improvements\n- Almost all plots are way too small and barely readable. If you reformat these, then I am willing to raise my \"Presentation\" score.\n- Link to the proofs in the appendix after each theoretical result, otherwise they are hard to find.\n- When you reference details in the appendix, link to the appropriate section (e.g. line 286).\n- Citations in lines 32 to 34 for the methods that are referenced are missing\n- Typos (l.211, l. 286, l. 362)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> “While correct as far as I can tell, I question the value of the theoretical result, if in practice you optimize the hyperparameters of the covariance function and ρ jointly. Is there reason to believe, that the convexity in ρ for fixed hyperparameters is beneficial given this choice?” “It would be informative to see an experiment that compares the non-convex parametrization to the convex one where all hyperparameters are optimized jointly”\n\nWe would expect the overall problem to be “more convex” than with the canonical parameterization, and thus improve optimization performance. More formally, the convexity guarantee will guarantee the positive-definiteness of the submatrix of the Hessian corresponding to the rhos, which is better than not guaranteeing positive definiteness at all.\n\nPositive-definiteness is also beneficial for quasi-Newton optimization algorithms like L-BFGS, which restarts the Hessian-approximation whenever it encounters non-convex regions, because the associated updates to the Hessian approximation are not positive-definite. This leads the algorithm to momentarily revert back to gradient descent, with an associated slower convergence rate.\n\nPractically, we ran convergence analyses using the data from Figure 1, allowing all rhos to be optimized jointly with other hyper-parameters (length scale, kernel variance and noise variance), recording the achieved negative log marginal likelihood (NLML) as a function of the tolerance parameter `ftol` of the L-BFGS optimizer. The results indicate that the optimizer terminates with a much better _NLML_ using the convex parameterization with the same convergence tolerance:\n\nftol | Canonical | Convex |\n----|----|--------|\n1e-03 | -4.37 | -14.18 |\n1e-04 | -4.37 | -93.00 |\n1e-05 | -4.37 | -93.01 |\n1e-06 | -4.37 | -135.05 |\n1e-07 | -98.62 | -518.68 |\n1e-08 | -97.52 | -1139.29 |\n\nThere are also settings in which we do not actually jointly optimize the hyperparameters. For instance, consider a situation in which we have access to data from the same data generating process that has been manually labeled by domain experts as outlier-free. Then we can estimate the model hyperparameters (of the non-robust GP) on that data, and fix those for the RGP-RP on the new data set that we do not know to be outlier-free (and that we cannot label due to cost or other practical reasons).\n\n> “Even though this is quite recent work, I would have liked to see a discussion of Altamirano et al. (2024) in the related work section”\n\nThanks for pointing out this relevant work. We have added a discussion of and comparison to this work in our general response, and will add this to the paper.\n\n> “Or is there an implicit assumption that each evaluation is significantly more expensive than fitting the GP?”\n\nThis is indeed the case; GPs are commonly applied to Bayesian optimization or active learning tasks where evaluation times are on the order of hours or days, so the fitting time of the GP is generally not much of a concern. We will make this assumption explicit. \n\n> “The experiments on regression problems could be expanded in my opinion, e.g. with benchmark datasets from the UCI repository”\n\nWe have done so, see our general response. Please let us know if there are any additional datasets you would like us to consider for the CR if accepted.\n\n> “I think the paper could be made stronger if there were at least one experiment where the corruption is not synthetically generated [...]”\n\nThanks for this feedback. We have added the “Twitter Flash Crash” example from Altamirano et al. (2024) and an additional variant of it to our results (see general response). If you have any other examples of real-world data we would be happy to consider this for inclusion in the CR.\n\n> “Almost all plots are way too small and barely readable”\n\nWe have adjusted the size of our plots (in particular the axis labels and legends), see attached pdf. We will utilize part of the additional space for the CR submission to increase plot size throughout.  \n\n> “Is your assumption about the data corruption different from Huber's ε-contamination model or is it identical? Equation 2 suggests they are different, but the text and experiments often refer to only a fraction of the data being corrupted [...]”\n\nA significant high-level difference is that Huber’s model is homoskedastic while ours is not. Citing Huber’s original article, he states: “Let $x_1,⋯, x_n$ be independent random variables with common distribution function $F$”, i.e. Huber’s model assumes that every data point is equally likely to be an outlier. \n\nIn contrast, our model introduces data-point specific variances that are adapted using marginal likelihood optimization, leading the resulting likelihood function to be different for each datapoint with a distinct rho. \n\nThe article often refers to the “fraction of data being corrupted” because the inference of this fraction is an important problem that determines the performance of RP, which requires choosing a number of outlier variances ($\\rho_i$). We solve this problem by leveraging Bayesian model selection.\n\n> “How does your approach compare to a vanilla GP on a dataset with no corruption? Does it \"fail gracefully”?”\n\nYes, and the results in our general response demonstrate this. We expect this to hold as long as the prior on the occurrence probability is sufficiently uninformative and puts some weight on no occurrences. Throughout our experiments, we used a geometric prior on the support size with a mean equal to 20% of the number of observations. We will add the results and discussion to the manuscript.\n\n> “How did you set the schedule to test data points for being outliers in the experiments (c.f. lines 182-184)?”\n\nSee [our response to reviewer aaaZ](https://openreview.net/forum?id=5FATPIlWUJ&noteId=OjPuKthMyu).\n\n> “I would have liked to see a bit more discussion on the limitations of the approach in Section 7 [...]”\n\nGood point, we will include these in the discussion in Section 7.", "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": "", "extra_scores": {}}]

Author rebuttal:
We thank the reviewers for their detailed comments and valuable suggestions. We are glad to see that the reviewers found our "novel approach” to be “clearly motivated and presented” and of “broad interest”, and that our work “provides a nice addition to the existing toolbox of methods”. 

While reviewers in general found our work “well-written” and “easy to follow”, there were some requests for improving clarity (e.g. Section 4.3), more readable figures, and a more comprehensive discussion of the regression results. We will gladly make these adjustments, and will utilize the additional page to improve legibility and include the additional experiments suggested by the reviewers and added here.

## RCGP Benchmarks
In response to the desire for additional regression results, we ran a comprehensive suite of benchmarks that include various UCI datasets, and included the RCGP method from the contemporaneous work Altamirano et al. (ICML 2024) as an additional baseline, using their publicly available benchmarking suite. The results are presented in Table 1 in the attached pdf. The main takeaways are:
- Overall, RGP-RP outperforms the baselines in terms of Mean Absolute Error (MAE) - the metric chosen by Altamirano et al - in 50% of the test cases, is tied with the standard GP in 40% of cases, and the Student-T GP (t-GP) is best in 10% cases. There is not a single case where RCGP or standard GP outperforms RGP-RP in a statistically significant way (standard errors in parenthesis).

- For uncorrupted data, all methods perform comparably. Notably, RGP-RP shows predictive performance indistinguishable from the standard GP, an indication that the RGP-RP’s Bayesian model selection correctly chooses the outlier-free model. 

- For “Uniform” and “Asymmetric” Outliers, as defined in Altamirano et al. 2024, RGP-RP outperforms the baselines for all test cases, showing that our results generalize to more problems.

- RGP-RP takes longer to fit than other baselines but achieves superior predictive performance. Fitting times in these benchmarks are on the order of seconds, and so this should not be a concern for most practical applications. Relative timings that each method took to complete the entire benchmark are:
    - GPR (GPFlow): 1.9x,
    - t-GP (GPFlow): 20.6x
    - RCGPR (GPFlow): 5.7x
    - Standard (BoTorch): 1.0x
    - RGP-RP (BoTorch): 36.2x

_Note: The original RCGP benchmarks from Altamirano et al. (2024) [hard-code the seed that controls the random train-test-splits](https://github.com/maltamiranomontero/RCGP/blob/aff281a39a6be6eefc15d96db8acba6c49224d28/experiments/uci/dataset_api.py#L321), therefore leading the results to only use a single train-test split, contrary to what is indicated in the paper. We edited the code to use different seeds for each replication._

## Twitter Flash Crash
As an additional real-world example, we consider the “twitter flash crash” from Altamirano et al. (2024), see Figure 1 in the attached pdf for results.

_Note: The original notebook sets the RCGP noise variance to the one inferred by GPR, and turns off the learning of the variance. In contrast, RCGP jointly learns the noise variance in the original benchmarks._

In the top row of Figure 1, we show results for 
- “RCGP (GPR Noise)”, the original method that forces the variance to the GPR value.
- “RCGP (Trained Noise)” sets the variance to be trainable.
- “RCGP (Fixed Noise)” forces the noise variance to 0.25.

The top row shows that RCGP with the original setting (GPR Noise), and (Fixed Noise) are less affected by the outliers than GPR, but still substantially so. Surprisingly, RCGP (Trained Noise) is even more affected. Further, the RGP-RP (“Relevance Pursuit”) is virtually unaffected by the outlying data points while modeling other high-frequency components of the time series closely.

## Twitter Flash Crash with Additional Training Day

A key shortcoming of RCGP is that its weighting function only uses the distribution of outcomes Y relative to the prior mean $m({\bf x})$. For this to be effective, the outliers need to be separable from the marginal data distribution, like, e.g., for winsorization. The reason RCGP shows improvements on the single-day example is that the outliers at the sell-off time are a-priori separated from most of the data. 

If this is not the case, RCGP can up-weight outliers, while down-weighting real data. We illustrate this failure in the bottom row of Figure 1, generated by included data from the preceding day. Here, RCGP is not robust to the outliers - in fact, the weighting function assigns a very large weight to the outlier, while RP-GP exhibits consistently strong performance.

## Additional Results: UCI Data using RP Setup

We also added UCI datasets (yacht, energy, CA housing) into our own benchmarking suite. Figure 2 in the pdf contains the results for the yacht data with “uniform” outliers. Note that RCGP’s “uniform” outlier process for Table 1 is different from ours:

- RCGP’s “uniform” outliers: generated by adding or subtracting a random value uniformly sampled between 3 and 9 standard deviations of the original data, thereby separating outliers from the data.

- RP’s “uniform” outliers: generated by uniformly sampling inside the range of the uncorrupted data.
This leads to the differences between the results reported in Figure 2 and Table 1 on the yacht data.

The results show that RP generally outperforms the baselines in terms of MAE and LL over a range of corruption percentages, and is faster than the robust Student-T and Trimmed-MLL approaches. The results on the CA housing and energy data are qualitatively similar. We will include them in the manuscript.

## An Ask
With these additional results, we believe that the main weaknesses of the paper cited by the reviewers – the regression experiments and comparison to additional baselines (RCGP) – were addressed comprehensively. Based on this, we would like to ask the reviewers to consider raising their score.

Meta-review: All reviewers recommend acceptance. Authors: when preparing the final version of your manuscript, please follow the recommendations and guidance provided by the reviewers. In particular, Section 3 should be improved by including further relevant references from the heteroscedastic GPs literature (e.g., "Heteroscedastic Gaussian Process Regression", "Variational Heteroscedastic Gaussian Process Regression", "Most likely heteroscedastic Gaussian process regression") and discussing how they connect with your method, in particular the lack of an assumption of noise variance smoothly varying with the input.
Decision: Accept (poster)