Robust Gaussian Processes via Relevance Pursuit

Authors: sebastian  ament
elizabeth  santorella
david  eriksson
ben  letham
maximilian  balandat
eytan  bakshy
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 19:46:42

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.

Main paper: 
 introduction: 
 Probabilistic models have long been a central part of machine learning, and Gaussian process (GP) models are a key workhorse for many important tasks [49] , especially in the small-data regime. GPs are 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) [22, 24] and active learning [6, 50] .GPs flexibly model a 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 a point x and ϵ ∼ N (0, σ 2 ), implying a homoskedastic Gaussian likelihood. While mathematically convenient, this assumption can be a limitation in practice, since noise distributions are often heavy-tailed or observations may be corrupted due to issues such as sensor failures, data processing errors, or software bugs. Using a 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 (e.g., Winsorizing), modified likelihood functions (e.g., Student-t), and model-based data selection and down-weighting procedures.These approaches offer different trade-offs between model accuracy, degree of robustness, broad applicability, computational requirements, and theoretical guarantees.In this paper, we propose a simple yet effective implicit data-weighting approach that endows GPs with a 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 a modification to the standard GP model that introduces learnable data-point-specific noise variances.2. We introduce a 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 a particular parameterization, the MLL is strongly concave in the datapoint-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 range, see e.g. Figure 1 . Comparison of RRP to a standard GP and a variational GP with a Student-t likelihood on a 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. 
preliminaries: 
 We aim to model a function f : X → R over some domain X ⊂ R d . With a standard Gaussian noise model, for x i ∈ X we obtain observations y i = f (x i ) + ϵ i , where ϵ i ∼ N (0, σ 2 ) are i.i.d. draws from a Gaussian random variable. ∥•∥ denotes the Euclidean norm unless indicated otherwise.A GP f ∼ GP (µ(•), k θ (•, •)) is fully defined by its mean function µ : X → R and covariance or kernel function k θ : X × X → R, which is parameterized by θ. Without loss of generality, we will assume that µ ≡ 0. Suppose we have collected data D = {(x i , y i )} n i=1 where X := {x i } n i=1 , y := {y i } n i=1 . Let Σ θ ∈ S n ++ denote the covariance matrix of the data set, i.e., [Σ θ ] ij = k θ (x i , x j ) + δ ij σ 2 , where δ ij is the Kronecker delta. The negative marginal loglikelihood (NMLL) L is given by -2L(θ) := -2 log p(y|X, θ) = y ⊤ Σ -1 θ y + log det Σ θ + n log 2π.In the following, we will suppress the explicit dependence of the kernel matrix on θ for brevity of notation. For a comprehensive background on GPs, we refer to Rasmussen et al. [49] .Additive, heavy-tailed noise Instead of assuming the noise term ϵ i in the observation model to be Gaussian, other noise models consider zero-mean perturbations drawn from distributions with heavier tails, such as the Student-t [30] , Laplace [35] , or α-Stable [5] distributions. These types of errors are common in applications such as finance, geophysics, and epidemiology [16] . Robust regression models utilizing Student-t errors are commonly used to combat heavy-tailed noise and outliers.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 a model of the form y i = Z i f (x i ) + (1 -Z i )W i , where Z i ∈ {0, 1} and W i ∈ R is a random variable. Note that W i need not have (and rarely has) f (x i ) as its mean. For instance, consider a faulty sensor that with some probability p reports a random value within the sensor range [y l , y h ]. In this case Z i ∼ Ber(p) and W i ∼ U[y l , y h ]. Software bugs, such as those found in ML training procedures, or errors in logging data can result in sparse corruptions.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 summarized in the conclusion.Guidelines:• The answer NA means that the paper has no 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 (e.g., independence assumptions, noiseless settings, model well-specification, 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 a few datasets or with a 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, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to 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 fairness. • 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 acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
related work: 
 Data pre-processing Data pre-processing can be an effective technique for handling simple forms of data corruption, such as values that fall outside a valid range of outputs. With such pre-processing, outliers are handled upstream of the regression model. Common techniques include the power transformations [13] , trimming, and winsorization. These methods can 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 [15] for a review on data cleaning.Heavy-tailed likelihoods One class of robust methods uses additive heavy-tailed noise likelihoods for GPs, particularly Student-t [30] , Laplace [48] , and Huber [1] , and could be extended with α-Stable distributions, which follow a generalized central limit theorem [5] . These models are less sensitive to outliers, but they lose efficiency when the outliers are a sparse subset of the observations, as opposed to global heavy-tailed noise. Furthermore, model inference is no longer analytic, necessitating the use of approximate inference approaches such as MCMC [44] , Laplace approximation [59] , expectation propagation (EP) [30] , Expectation Maximization [48] , or variational inference [56] . Shah et al. [51] take a related approach using a 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 a hierarchical mixture of Gaussians [17] and a "twinned" GP model [43] 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. [10] 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 [9] fit GPs using the trimmed MLE by applying a projected gradient method to an approximation of the marginal likelihood. The associated theory only guarantees convergence to a stationary point, with no guarantee on quality. When no outliers are present, this method can be worse than a standard GP. Li et al. [36] propose a heuristic iterative procedure of removing those data points with the largest residuals after fitting a standard GP, with subsequent reweighting. The method shows favorable empirical performance but has no theoretical guarantees, and fails if the largest residual is not associated with an outlier. Park et al. [47] consider a model of the formy i = δ i + f (x i ) + ϵ i ,where outliers are regarded as data with a large bias δ i . Their random bias model is related to our model in that it also introduces learnable, data-point-specific variances. However, inference is done in one step by optimizing the NMLL with an inverse-gamma prior on the δ i 's, which -in contrast to the method proposed herein -generally does not lead to exactly sparse δ i 's .Sample re-weighting Altamirano et 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 σ 2 I is replaced by σ 2 diag(w -2 ) and the prior mean m is replaced bym w = m + σ 2 ∇ y log(w 2). The authors advocate for the use of the inverse multi-quadratic weight function w(x, y) = β(1 + (y -m(x)) 2 /c 2 ) -1/2, which introduces two additional hyper-parameters: the soft threshold c, and the "learning rate" β. Importantly, the weights w are defined a-priori as a function of the prior mean m(x) and the targets y, thereby necessitating the weights to 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.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?• 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 or the standard error of the mean. Justification: Total compute resources spent are described in Section D.4 and are estimated at around 2 CPU years in total.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs 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). 
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. We extend the standard GP observation noise variance σ 2 with data-point-specific noise variances ρ = {ρ i } n i=1 , so that the i-th data point is distributed asEQUATIONThis is similar to Sparse Bayesian Learning [55] in which weight-specific prior variances control a feature's degree of influence on a model's predictions. The likelihood in ( 2) is similar in that the marginal likelihood optimization of the ρ i gives rise to an automatic mechanism for outlier detection and down-weighting. Intuitively, the effect of y i on the estimate of f vanishes as ρ i → ∞. This model can be viewed as a type of GP with heteroskedastic noise model. However, while many heteroskedastic models model noise as an input-dependent process [26, 33] , our formulation does not make such assumptions, and is thus suitable for corruptions that are not spatially correlated.An elegant consequence of our modeling assumption is that we can compute individual marginallikelihood maximizing ρ i 's in closed form when keeping all ρ j for j ̸ = i fixed. In particular, Lemma 1.[Optimal Robust Variances] Let D \i = {(x j , y j ) : j ̸ = i}, ρ = ρ \i + ρ i e i , where ρ, ρ \i ∈ R n + , [ρ \i ] i = 0, and e i is the ith canonical basis vector. Then keeping ρ \i fixed,EQUATIONwhere y(x i ) = f (x i )+ϵ i . These quantities can be expressed as functions of Σ -1 = (K+D σ 2 +ρ ) -1 :E[y(x i )|D \i ] 2 = y i -Σ -1 y i Σ -1 ii , and V[y(x i )|D \i ] = 1 Σ -1ii , where D σ 2 +ρ is a diagonal matrix whose entries are σ 2 + ρ. 3) is the empirical error to y i of the model trained without the i-th data point, i.e., the leave-one-out (LOO) cross-validation error [49] . The second componentThe first component E[f (x i ) + ϵ i |D \i ] 2 of (V[f (x i ) + ϵ i |D \i ]is the LOO predictive variance. The optimal solution to ρ i is only non-zero for those observations whose squared LOO error is larger than the LOO predictive variance at that point.Without additional structure, inference of the noise variances ρ i does not yield desirable models, as the marginal likelihood can be improved by increasing the prior variance ρ i 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 ρ i , that is, ∥ρ∥ ∞ = |{0 < ρ i }| ≤ k < n. While this sparsity constraint mitigates over-flexibility, it gives rise to a formidably challenging optimization problem, as there are a combinatorial number of sparse outlier sets to consider. Even if the number of outliers n o were known, exhaustive search would still require considering n-choose-n o possibilities.For tractability, we iteratively add data points to a set of potential "outliers" by setting their associated ρ i to be nonzero, using the closed-form expression for the optimal individual ρ i variances in Lemma 1. 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 1 with useBayesianModelSelection as false. Specifically, this is the "forward" variant; Algorithm 2 in the Appendix presents an alternative "backward" variant that we found to work well if 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-weighted if it is apparently an outlier. This allows the down-weighting to be reversed if a data point appears "inlying" after having down-weighted other data points, improving the method's robustness and performance. This is in contrast to Andrade and Takeda [9] '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 is 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 1 , k 2 , . . . , k K ), useBayesianModelSelection (boolean)Initialize S 0 ⊆ {1, . . . , n} (typicallyS 0 = ∅) for i in (1, . . . , |K|) do Optimize MLL: ρ Si ← arg max ρ S i L ρ Si , where ρ Si = {ρ : ρ j = 0, ∀ j ̸ ∈ S i }. Expand Support: Compute ∆ i (j) ← max ρj L(ρ Si + ρ j e j ) -L(ρ Si ) for all j ̸ ∈ S i via Lemma 1 . A i ← {j 1 , . . . , j ki } such that ∆ i (j) ≥ ∆ i (j ′ ) for all j ∈ A i and j ′ ̸ ∈ (A i ∪ S i ). S i+1 ← S i ∪ A i if useBayesianModelSelection then Compute the marginal likelihood p(D|S i ) ≈ p(D|S i , ρ Si ) S * ← arg max Si p(D|S i )p(S i ). else S * = S K . Return S * , ρ S * .In practice, it is often impossible to set a hard threshold on the number of outliers for a 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 a prior distribution p(S) over the number of outliers, rather than fix the number a priori. We leverage the Bayesian model selection framework [61, 40] to determine the most probable number of outliers in a data-and model-dependent way, aiming to maximize p(S|D). This gives rise to Algorithm 1, with useBayesianModelSelction as 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 a trace of models generated at each iteration, then approximate the model posterior p(S i |D) ∝ p(D|S i )p(S i ) at each point in the trace. As the exact posterior is intractable, we approximate it with p(D|S i ) = p(D|S i , ρ Si )dρ Si ≈ p(D|S i , ρ * Si ). Finally, we select the model from the model trace {S i } i 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 l 1 -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 as tightly as possible.Regarding the schedule K in Algorithm 1, the most natural choice is simply to add one data point at a time, i.e. K = (1, 1, ...), but this can be slow for large n. In practice, we recommend schedules that test a fixed set of outlier fractions, e.g. K = (0.05n, 0.05n, . . . ). 
theoretical analysis: 
 We now provide a theoretical analysis of our approach. We first propose a re-parameterization of the ρ i that maps the optimization problem to a compact domain. Surprisingly, the re-parameterized problem exhibits strong convexity and smoothness when the base covariance matrix (excluding the ρ i ) 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), a greedy algorithm that is closely related to Algorithm 1.The optimization of linear models with respect to least-squares objectives in the presence of sparsity constraints has been richly studied in statistics [54] , compressed sensing [3, 57] , and machine learning [8, 62] . Of central importance to the theoretical study of this problem class are 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 × m)-matrix A satisfies the r-restricted isometry property (RIP) with constant δ r ∈ (0, 1) if for every submatrix A S with |S| = r ≤ m columns,(1 -δ r )∥x∥ ≤ ∥A S x S ∥ ≤ (1 + δ r )∥x∥, where x S ∈ R r . This is equivalent to (1 -δ r )I ⪯ (A * S A S ) ⪯ (1 -δ r )I.The RIP has been proven to lead to exact recovery guarantees [14] , as well as approximation guarantees [18] . Elenberg et al. [21] 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 r -restricted strong convex and M r -restricted smooth if for all (x, x ′ ) in the domain D r ⊂ (R d × R d ), m r ∥x ′ -x∥ 2 /2 ≤ f (x ′ ) -f (x) -∇[f ](x) ⊤ (x ′ -x) ≤ M r ∥x ′ -x∥ 2 /2.In the context of sparse optimization, we let D r be the set of tuples of r-sparse vectors whose difference is also at most r-sparse. In particular,D r = {(x, x ′ ) s.t. ∥x∥ 0 , ∥x ′ ∥ 0 , ∥x ′ -x∥ 0 ≤ r}.Generalized orthogonal matching pursuit (OMP) [4, 39] 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 magnitudes, applied to the marginal liklihood optimization problem, S i+1 = S i ∪ arg max j̸ ∈S |∇ ρ L(ρ)| j . Algorithm 1 generalizes OMP [57] by allowing more general support expansion schedules K, and specializes the support expansion criterion using the special problem structure exposed by Lemma 1.The NMLL L of a GP (1) is the sum of a convex function (•) -1 and a concave function log det(•) of K, and is therefore not generally convex as a function of the hyper-parameters θ, including the robust variances ρ. Here, we propose a re-parameterization that allows us to prove strong convexity guarantees of the associated NMLL. In particular, we let ρ(s) = diag(K 0 ) ⊙ ((1 -s) -1 -1), where K 0 := k(X, X) + σ 2 I and the inverse is element-wise. Note that ρ(s) is a diffeomorphism that maps s from the compact domain s∈ [0, 1] n to the entire range of ρ ∈ [0, ∞] n .Henceforth, we refer to the original ρ as the canonical or ρ-parameterization and the newly proposed ρ(s) as the convex or s-parameterization. Lemma 4 shows the Hessian of the s-parameterization.Lemma 4. [Reparameterized Hessian] Let K s = k(X, X) + σ 2 I + D ρ(s) , Ks = diag(K s ) -1/2 K s diag(K s ) -1/2 , and α = K-1 s diag(K s ) -1/2 y. Then H s [-2L(ρ(s)] = D -1 1-s 2 α α⊤ ⊙ ( K-1 -I) + 2 diag( K-1 ) -( K-1 ⊙ K-1 ) D -1 1-s .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.EQUATIONwhere λ min,max (resp. λmin,max ) are the smallest and largest eigenvalues of K s , respectively Ks . The behavior Lemma 5 predicts is surprising and validated in Fig. 3 . 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 is intimately linked to the RIP (Def. 2). Note that Lemma 5 is a condition for non-support-restricted convexity, which is stronger than is necessary for 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] n . Definition 6 (Diagonal Dominance). A matrix A is said to be δ-diagonally dominant if the elementsa ij satisfy i̸ =j |a ij | < δ|a ii | for all i.Intuitively, the ρ i (s) that are selected to be non-zero by the greedy algorithm take on large values, further encouraging the diagonal dominance of the sub-matrix of K associated with the support of ρ.For this reason, the following condition on K 0 is sufficient to guarantee convexity for all s ∈[0, 1] n . Lemma 7. [Strong Convexity via Diagonal Dominance] Let m > 0 and K 0 be δ-diagonally dominant with δ < (5 -m) - √ 25 -9m + 17 /4 ≤ (5 - √ 17)/4 ≈ 0.44 and λ min (K 0 )(1 -δ) 2 2(1 + δ) -1 -(1 -δ) -2 -m 2(1 -(1 -δ)/(1 + δ)) ≥ ∥y∥ 2 2 .Then the NMLL is m-strongly convex for all s∈ [0, 1] n , i.e. ρ(s) ∈ [0, ∞] n .We attain similar results for M -smoothness, see Lemma 13 and Lemma 14 in the Appendix. Having proven m-convexity and M -smoothness conditions, we appeal to the results of Elenberg et al. [21] . Theorem 8. [Approximation Guarantee] Let K 0 = k(X, X) + σ 2 I be δ-diagonally dominant, s max > 0 be an upper bound on ∥s∥ ∞ , and suppose ∥y∥, δ satisfy the bounds of Lemmas 7 and 14, guaranteeing m-convexity and M -smoothness of the NMLL for some m > 0, M > 1/(1 -s max ) 2 . Let s OMP (r) be the r-sparse vector attained by OMP on the NMLL objective for r steps, and let s OPT (r) = arg max ∥s∥0=r,∥s∥∞≤smax L(ρ(s)) be the optimal r-sparse vector. Then for any 2r ≤ n,L (ρ(s OMP (r))) ≥ 1 -e -m/M L (ρ(s OPT (r))) ,whereL(•) = L(•) -L(0) is normalized so that max s S L(s S ) ≥ 0 for any support S.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 is common to most marginal likelihood optimization problems. In practice, we typically optimize ρ jointly with the other hyper-parameters of the model in each iteration of RRP, as this yields improved performance, see App. D.5 for details. 
empirical results: 
 We evaluate the empirical performance of RRP against various baselines on a 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.") [11] , and a power transformation ("Power Transf.") [13] . Further, we also consider a Student-t likelihood model from Jylänki et al. [30] ("Student-t"), the trimmed marginal likelihood model from Andrade and Takeda [9] ("Trimmed MLL"), and the RCGP model from Altamirano et al. [2] . Unless stated otherwise, all models are implemented in GPyTorch [23] and all experiments in this section use 32 replications. See Appendix D for additional details.Synthetic We first consider the popular Friedman10 and Hartmann6 [19] test functions from the literature. We use two data generating processes: uniform noise, extreme outliers at some fixed value, and heavy-tailed (Student-t) noise at true function values. In these experiments, we compare the performance predictive log-likelihood. The results are shown in Fig. 4 . Twitter Flash Crash In Fig. 5 , we report a 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 a sharp drop at 13:10 on the 23rd. The top panels shows that RCGP exhibits higher robustness than the standard GP, but is 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-productive when the outliers are not a-priori separated in the 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 to "trust" the outliers more than any inlier, resulting in it being less robust than a standard GP in this scenario. See Appendix D.6 for additional comparisons to RCGP, on data sets from the UCI machine learning repository [32] . The accompanying panels labeled w imq show the function that Altamirano et al. [2] 's RCGP uses to down-weight data points. Top: RCGP, exhibits higher robustness than the standard GP, but is 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 to assign the highest weight w imq to the outlying data points due to their proximity to the target values' median, thereby leading RCGP to be even more affected than a standard GP, see (d) for a detailed view of the results on the data of April 23.GPs are commonly used for Bayesian optimization (BO), which is a popular approach to sampleefficient black-box optimization [22] . However, many of the GP models used for BO are sensitive to outliers and may not perform well in settings where such outliers occur. While Martinez-Cantin et al. [41] consider the use of a Student-t 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. [7] 's qLogNoisyExpectedImprovement (qLogNEI), a variant of the LogEI 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. [29] and plot the true value of the best in-sample point according to 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.We consider the popular 6-dimensional Hartmann test function with three different corruption settings: (1) a constant value of 100, (2) a U [-3, 3] distributed value, (3) the objective value for a randomly chosen point in the domain. The results for a 10% corruption probability are shown in Fig. 6 . We also include results for a 20% corruption probability in Appendix D.3. Real-world problems We include three real-world problems: A 3D SVM problem, a 5D CNN problem, and a 20D rover trajectory planning problem, see the App. D.2 for details. For SVM and CNN, we simulate random corruptions corresponding to an 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 we follow the setup in [42] with the main difference that we consider a 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.  
conclusion and future work: 
 Contributions Robust Gaussian Processes via Relevance Pursuit (RRP) provides a novel and principled way to perform robust GP 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 flexible, e.g., can be used with any mean or kernel function. Our method can be readily applied to both robust regression problems as well as applications such as Bayesian optimization and is available through BoTorch [12] . Importantly, it also provides theoretical approximation guarantees.Limitations As our approach does not explicitly consider the locations of the data points in the outlier identification, it may be outperformed by other methods if the underlying noise is heteroskedastic and location-dependent. On the other hand, those methods generally do not perform well in the presence of sparse, location-independent data corruptions.Extensions 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. [37] 's scalable learning-curve model for AutoML applications, and leveraging forward-backward greedy algorithms. On a higher level, the approach of combining greedy optimization algorithms with Bayesian model selection and leveraging a convex parameterization to achieve approximation guarantees might apply to other parameters that are optimized using the MLL objective: length-scales of stationary kernels, coefficients of additive kernels, inducing inputs, and even related model classes like Tipping [55] 's Sparse Bayesian Learning (SBL), which seeks to identify sparse linear models and is intimately linked to greedy matching pursuits [8] . Overall, the approach has the potential to lead to theoretical guarantees, new insights, and performance improvements to widely-adopted Bayesian models. 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA]Justification: The paper makes foundational methodological contributions to robust probabilistic regression, and we do not see any direct societal impacts of the work that would require specific discussion.Guidelines:• 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 negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.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 a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an 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 to 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 or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 a 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 or 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 or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a 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 not available online, the authors are encouraged to reach out to the asset's creators. 
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 new 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 as an 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 is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing 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 as details about compensation (if any)?Answer: [NA] Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
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 crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing 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 we expect authors to adhere to 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. 
 
 Appendix: Algorithm 2 below is the "backward" variant of Algorithm 1 from Sec 4. As its name suggests, the main difference compared to the "forward" variant is that rather than building up a 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-dependent noise variance ρ i .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 our studies is that while the "forward" version tends to perform slightly better than the "backward" version if there are only few outliers, the opposite is true if the outlier frequency is 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 = (k 1 , k 2 , . . . ) Initialize S c 0 ⊆ {1, . . . , n} (typically S c 0 = {1, . . . , n}) for k i in K do Optimize ML: ρ * S c i ← arg max ρ S c i L(ρ S c i ), where ρ S c i = {ρ : ρ j = 0, ∀ j / ∈ S c i } Compute the set R i containing the k i elements of S ci with smallest inferred variance: R i ← {j 1 i , . . . , j ki i } wherej l i ∈ S c i such that ρ * S c i (j l i ) ≤ ρ * S c i (j l ′ i ) for l < l ′ S c i+1 ← S c i \ R i S i ← {1, . . . , n} \ S c i for each k i Compute the marginal likelihood p(S i |X, y) ≈ p(S i , ρ * Si , |X, y) S * ← arg max Si p(S i |X, y)p(S i ). Return S * , ρ * S * .Das and Kempe [18] showed that the subset selection problem of regression features with an R 2 objective satisfies a weak submodularity property, which can be invoked to prove approximation guarantees for the greedy maximization of the objective. Elenberg et al. [21] 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. [31] contains a guarantee for the backward algorithm applied to the maximization of submodular set functions.Submodularity has also found application to Gaussian process models. For a sensor placement problem, Krause et al. [34] 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 k [45] . Srinivas et al. [52] used the submodularity of the joint entropy in order to prove regret bounds for the convergence of a GP-based BO algorithm using the upper-confidence bound acquisition function.Elenberg et al. [21] proved that any log likelihood function exhibiting restricted strong concavity and smoothness implies the weak submodularity of the associated subset selection problem.Definition 9 (Submodularity Ratios [21] ). Let A, B ⊂ [n] be two disjoint sets, and f : 2[n] → R.The submodularity ratio of B with respect to A is defined byγ B,A = i∈A (f (B ∪ {i}) -f (B)) (f (B ∪ A) -f (B))The submodularity ratio of a set C with respect to an integer k is defined byγ C,k = min B,A γ B,A such that A ∩ B = ∅, B ⊆ C, and |A| ≤ k.Then given γ > 0, a function is γ-weakly submodular at a set C with respect to k if γ C,k ≥ γ.Theorem 10 (Weak Submodularity via RSC [21] ). The submodularity ratio γ S,k can be bound below using the restricted convexity and smoothness parameters m |S|+k and M |S|+k ,γ S,k ≥ m |S|+k M |S|+k .Theorem 11 (OMP Approximation Guarantee [21] ). Let x OMP (r) be the r-sparse vector selected by OMP, and x OPT (r) = arg min ∥x∥0=r f (x) be the optimal r-sparse vector. Thenf (x OMP (r)) ≥ 1 -e -m2r/M2r f (x OPT (r)) ,where m 2r , M 2r are the restricted strong convexity and smoothness parameters of f , respectively.Lemma 1.[Optimal Robust Variances] Let D \i = {(x j , y j ) : j ̸ = i}, ρ = ρ \i + ρ i e i , where ρ, ρ \i ∈ R n + , [ρ \i ] i = 0, and e i is the ith canonical basis vector. Then keeping ρ \i fixed,EQUATIONwhere y(x i ) = f (x i )+ϵ i . These quantities can be expressed as functions of Σ -1 = (K+D σ 2 +ρ ) -1 :E[y(x i )|D \i ] 2 = y i -Σ -1 y i Σ -1 ii , and V[y(x i )|D \i ] = 1 Σ -1 ii ,where D σ 2 +ρ is a diagonal matrix whose entries are σ 2 + ρ.Proof. First, we partition the covariance matrix K + σ 2 I + D ρ to separate the effect of ρ i and use the block matrix inverseEQUATIONwhereEQUATIONEQUATIONQuadratic Term With the expression for the inverse of Σ above, we can write the quadratic term of the log likelihood asEQUATION)Determinant Term The determinant of a block matrix is given byEQUATIONApplying this identity to Σ = (K + D ρ+σ 2 ), we getEQUATIONEQUATIONNoting that ∂ ρi β i = -β 2 i , the derivative of the difference in log marginal likelihood w.r.t. ρ i , isEQUATIONWhile β i = 0 is a root of the derivative, we ignore this solution since β i is never zero when σ 2 > 0.Therefore, the remaining stationary point is β -1 i = (y ⊤ \i u -y i ) 2 . Since we constrain ρ ≥ 0, this point might not always be attainable. However, because there is only a single stationary point with respect to β i when σ 2 > 0, and β i is a strictly decreasing function of ρ i , it follows that the marginal likelihood is monotonic as a function of ρ i to both the left and the right of the stationary point. Therefore, the optimal constraint ρ i is simply the optimal unconstrained value, projected into the feasible space. In particular, solving β -1 i = (y ⊤ \i u -y i ) 2 for ρ i and projecting to the non-negative half-line, we getEQUATIONLastly, noteEQUATIONAs stated by Rasmussen et al. [49] (P. 117, Eq. 5.12), originally shown by Sundararajan and Keerthi [53] , these quantities can be expressed as simple functions of Σ -1 :EQUATIONTherefore, all LOO predictive values can be computed in O(n 3 ) or faster, if an inducing point method is used for K.The following is a preliminary result for our analysis of the log marginal likelihood w.r.t. ρ. Lemma 12. The gradient and Hessian of the log marginal likelihood L with respect to ρ are given byEQUATION-2H ρ [L] = (2αα ⊤ -K -1 ) ⊙ K -1, where K = K 0 + D ρ for some base covariance matrix K 0 and α = K -1 y.Proof. Let α = K -1 y. Regarding the gradient, note that ∂ ρi K = e i e ⊤ i , where e i is the canonical basis vector with a one as the ith element, and based on Equation 5.9 of Rasmussen et al. [49] ,-2∂ ρi [L] = tr K -1 -αα ⊤ ∂ ρi K = tr( K -1 -αα ⊤ e i e ⊤ i ) = e ⊤ i K -1 -αα ⊤ e i = K -1 -αα ⊤ii . Regarding the second derivatives, according to Dong et al. [20] ,∂ θi ∂ θj [log |K|] = tr(K -1 [∂ θi ∂ θj K] -K -1 [∂ θi K]K -1 [∂ θj K]) ∂ θi ∂ θj [y ⊤ K -1 y] = 2α T [∂ θi K]K -1 [∂ θj K]α -α ⊤ [∂ θi ∂ θj K]α. Therefore, -2∂ ρi ∂ ρj L = tr (K -1 -αα T )[∂ ρi ∂ ρj K] -(K -1 -2αα ⊤ )([∂ ρi K]K -1 [∂ ρj K]) = tr (2αα ⊤ -K -1 )([K -1 ] ij e i e ⊤ j ) = tr e ⊤ j (2αα ⊤ -K -1 )e i [K -1 ] ij = [2αα ⊤ -K -1 ] ij [K -1 ] ij = (2αα ⊤ -K -1 ) ⊙ K -1 ij , since [∂ ρi ∂ ρj K] = 0.The third equality is 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.Here, we re-parameterize ρ(s) = diag(K 0 ) ⊙ ((1 -s) -1 -1), and attain strong convexity for all inputs s, if conditions on the eigenvalues of the covariance matrix and the norm of the data vector ∥y∥ are met. The convexity result is surprising in 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-t likelihood or α-stable likelihoods are non-convex, and even Huber's proposal is non-strongly-convex.Lemma 4. [Reparameterized Hessian] Let K s = k(X, X) + σ 2 I + D ρ(s) , Ks = diag(K s ) -1/2 K s diag(K s ) -1/2 , and α = K-1 s diag(K s ) -1/2 y. Then H s [-2L(ρ(s)] = D -1 1-s 2 α α⊤ ⊙ ( K-1 -I) + 2 diag( K-1 ) -( K-1 ⊙ K-1 ) D -1 1-s .Proof. Using the chain-rule, the Hessian H s [L] can be expressed as a function of the Jacobian J s [ρ(s)] = D ∂ρ(s) , which is diagonal since ρ(s) is an element-wise function, and the second derivatives∂ 2 s ρ(s i ). Then -2H s [L] = -J s [ρ] ⊤ H ρ [2 log L]J s [ρ] + D ∇ρ[-2L] D ∂ 2 s ρ = D ρ ′ (s) [(2αα ⊤ -K -1 ) ⊙ K -1 ]D ρ ′ (s) + diag(K -1 )D ρ ′′ -D ρ ′ (K -1 • K -1 )D ρ ′ ,where we substituted the relevant expressions from Lemma 12. Further substituting ρ(s) i = [K 0 ] ii [(1 -s i ) -1 -1], ρ ′ (s) i = [K 0 ] ii (1 -s i ) -2 , and ρ ′′ (s) i = 2[K 0 ] ii (1 -s i ) -3 , noting that K = K 0 + diag(K 0 )[(1 -s) -1 -1] = (K 0 -diag(K 0 )) + diag(K 0 )(1 -s) -1, and algebraic manipulation finish the proof.EQUATIONwhere λ min,max (resp. λmin,max ) are the smallest and largest eigenvalues of K s , respectively Ks .Proof. We seek to lower-bound the smallest eigenvalue of the Hessian matrix, which-for twicecontinuously-differentiable problems-is equivalent to lower and upper bounds of the problem's Hessian matrix. Starting with the result of Lemma 4,H s [-2L(ρ(s)] = D -1 1-s 2 α α⊤ ⊙ ( K-1 -I) + 2 diag( K-1 ) -( K-1 ⊙ K-1 ) D -1 1-s ⪰ 2 α α⊤ ⊙ ( K-1 -I) + 2 diag( K-1 ) -( K-1 ⊙ K-1 ).Now, we bound each of the three additive terms independently from below.Term 1: 2 diag( K-1 ) ⪰ 2λ min ( K-1 )I = 2λ max ( K) -1 I. The first inequality comes from K being positive definite, and the absolute value of the diagonal of a matrix, which is already positive for positive definite matrices, being lower bounded by the minimum eigenvalue of the matrix. The last steps is 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:λ min ( K) = λ min (diag(K) -1/2 K diag(K) -1/2 ) ≥ λ min (diag(K) -1/2 ) 2 λ min (K) = λ min (diag(K) -1 )λ min (K) ≥ λ min (K -1 )λ min (K) = λ min (K)/λ max (K).In a similar way, we can show that λ max ( K) ≤ λ max (K)/λ min (K), which implies λ max ( K) -1 ≥ λ min (K)/λ max (K).Term 2: Next, we have-( K-1 ⊙ K-1 ) ⪰ -λ max ( K-1 ) 2 I = -λ min ( K) -2 I,which is due to the Hadamard product being a 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 eigenvalues of the constituent matrices.Term 3: Lastly, note that2 α α⊤ ⊙ ( K-1 -I) = 2D α( K-1 -I)D α ⪰ 2D α(λ min ( K-1 )I -I)D α = 2(λ min ( K-1 ) -1)D 2 α = 2(λ max ( K) -1 -1)D 2 α ⪰ 2(λ min (K)/λ max (K) -1)D 2 α ⪰ 2(λ min (K)/λ max (K) -1)∥ α∥ 2∞ I, where the last inequality comes from the second to last lower bound of Term 3 being non-positive, and therefore, being able to lower bound it with the largest magnitude entry of α.Term 1 + 2 + 3: Putting together the inequalities for all terms, we getλ min (H s ) ≥ 2λ max ( K) -1 -λ min ( K) -2 + 2(λ min (K)/λ max (K) -1)∥ α∥ 2∞ , where in slight abuse of notation, we let H s = H s [-2L] be the Hessian of the negative log likelihood.Using the bound ∥ α∥ ∞ ≤ ∥ α∥ 2 ≤ λ min ( K) -1 ∥ŷ∥, where ŷ = diag(K) -1/2 y. Therefore, ∥ α∥ ∞ ≤ λ min ( K) -1 λ min (K) -1/2 ∥y∥ 2 ≤ λ max (K)λ min (K) -3/2 ∥y∥ 2 .Finally, lower bounding the current lower bound by m > 0 yields a sufficient condition for the convexity at s.2λ max ( K) -1 -λ min ( K) -2 + 2(λ min (K)/λ max (K) -1)λ min ( K) -2 λ min (K) -1 ∥y∥ 2 2 > m.EQUATIONwhich is a non-trivial guarantee when 2 λ-1 max -λ-2 min -m > 0.Lemma 13. [Smoothness via Eigenvalue Condition] Let K as in Lemma 4. Suppose that ∥s∥ ∞ ≤ s max andEQUATION)where λ min,max (resp. λmin,max ) are the smallest and largest eigenvalues of K (resp. K) respectively. Then H s ≺ M . This is a non-trivial guarantee when M (1-s max ) 2 + λ-2 min -2 λ-1 max > 0.Proof. We now derive an equivalent upper bound for the largest eigenvalue of the Hessian. Starting with the result of Lemma 4,H s [-2L(ρ(s)] = D -1 1-s 2 α α⊤ ⊙ ( K-1 -I) + 2 diag( K-1 ) -( K-1 ⊙ K-1 ) D -1 1-s ⪯ ∥(1 -s) -1 ∥ 2 ∞ 2 α α⊤ ⊙ ( K-1 -I) + 2 diag( K-1 ) -( K-1 ⊙ K-1 ) .Therefore, it becomes immediately apparent that we will need to introduce an upper bound on∥(1 -s) -1 ∥ 2 ∞, which is a restriction on the domain that ρ can take. Proceeding in a similar way as above, we bound the three terms in square brackets, now from above.Term 1:2 diag( K-1 ) ⪯ 2λ max ( K-1 )I = 2λ min ( K) -1 I. Term 2: -( K-1 ⊙ K-1 ) ⪯ -λ min ( K-1 ) 2 I = -λ max ( K) -2 I. Term 3: 2 α α⊤ ⊙ ( K-1 -I) = 2D α( K-1 -I)D α ⪯ 2D α(λ max ( K-1 )I -I)D α = 2(λ max ( K-1 ) -1)D 2 α = 2(λ min ( K) -1 -1)D 2 α ⪯ 2(λ max (K)/λ min (K) -1)D 2 α ⪯ 2(λ max (K)/λ min (K) -1)∥ α∥ 2 ∞ I, where now the last inequality follows because the second to last expression is always non-negative. Term 1 + 2 + 3: Putting together the inequalities for all terms, we getλ max (H s ) ≤ 2λ min ( K) -1 -λ max ( K) -2 + 2(λ max (K)/λ min (K) -1)∥ α∥ 2∞ , Finally, upper bounding the current upper bound by M/∥(1 -s) -1 ∥ 2 ∞ > 0 yields a sufficient condition for the M -smoothness at s. Using the same bound for ∥ α∥ ∞ derived for the convexity result,2λ min ( K) -1 -λ max ( K) -2 + 2((λ max (K)/λ min (K) -1)λ min ( K) -2 λ min (K) -1 ∥y∥ 2 2 < M.EQUATIONwhich is a non-trivial guarantee when M (1 -s max ) 2 + λ-2 min -2 λ-1 max > 0. λ min (K 0 )(1 -δ) 2 2(1 + δ) -1 -(1 -δ) -2 -m 2(1 -(1 -δ)/(1 + δ)) ≥ ∥y∥ 2 2 .Then the NMLL is m-strongly convex for all s ∈ [0, 1] n , i.e. ρ(s) ∈ [0, ∞] n .Proof. Fist, Gershgorin's Disk Theorem implies that the eigenvalues of diag(A) -1/2 A diag(A) -1/2 lie in (1 -δ, 1 + δ) for a δ-diagonally dominant matrix A. Further, the condition number of A is bounded above by κ(A) = λ max (A)/λ min (A) ≤ (1 + δ)/(1 -δ). See Horn and Johnson [28] for more background on matrix analysis. Plugging these bounds into the results of Lemma 5 yieldsλ min (K) λ2 min (2 λ-1 max -λ-2 min -m) 2(1 -λ min /λ max ) ≥ λ min (K)(1 -δ) 2 (2(1 + δ) -1 -(1 -δ) -2 -m) 2(1 -(1 -δ)/(1 + δ)) .Lower bounding the last expression by ∥y∥ 2 2 implies m-strong convexity. This gives rise to a non-trivial guarantee whenever the numerator is larger than zero. In particular,2(1 + δ) -1 -(1 -δ) -2 -m > 0.Expanding, noting that 0 < δ < 1 implies δ 3 < δ 2 in order to reduce a power in the resulting expression, and collecting terms with like powers, we attain the following sufficient condition2δ 2 + (m -5)δ + (1 -m) > 0.Note that for this condition to hold if δ = 0, we need to have m < 1. Fortunately, this is a quadratic in δ whose smallest positive root isδ -= 1 4 (5 -m) -(5 -m) 2 -8(1 -m) .In particular, for m = 0, this reduces toδ -= 1 4 5 - √ 17 ≈ 0.4384471871911697.Lastly, note that if K 0 is δ-diagonally dominant, then so is K s = K 0 + D ρ (s), since the robust variances add to the diagonal, making it more dominant. Therefore, the convexity guarantee holds for all ρ(s), if it holds for the base covariance matrix K 0 . Note that λ min (K) ≥ σ 2 .Similarly, we can prove a similar statement relating diagonal dominance with M -smoothness. Lemma 14. [Smoothness via Diagonal Dominance] Suppose K 0 is a δ-diagonally dominant covariance matrix and suppose we constrain ∥s∥ ∞ ≤ s max ≤ 1 -1/M . Thenλ min (K 0 )(1 -δ) 2 M (1 -s max ) 2 -1 2((1 + δ)/(1 -δ) -1) ≥ ∥y∥ 2 2 ,implies that the NLML is m-strongly convex for all s ∈ [0, 1] n , i.e. ρ(s) ∈ [0, ∞] n .Proof. We proceed in a similar way as for Lemma 7, but with Lemma 13 as the starting point.λ min (K)( λmin ) 2 (M (1 -s max ) 2 + λ-2 min -2 λ-1 max ) 2(λ max /λ min -1) ≥ λ min (K)(1 -δ) 2 M (1 -s max ) 2 -1 2((1 + δ)/(1 -δ) -1) ,where we used ( λ-2 min -2 λ-1 max ) ≥ -1, which is tight when K is unitary. Lower bounding the last expression by ∥y∥ 2 2 implies M -smoothness. This gives rise to a non-trivial guarantee whenever the numerator is larger than zero. In particular, M (1-s max ) 2 -1 > 0, which implies s max ≤ 1-1/M or equivalently, M > 1/(1 -s max ) 2 . Lastly, note that if K 0 is δ-diagonally dominant, then so is K s = K 0 + D ρ (s), since the robust variances only add to the diagonal. Therefore, if the inequality holds for the base covariance matrix K 0 , the smoothness guarantee holds for all ρ(s) such that s ≤ s max ≤ 1 -1/M . Note also that λ min (K) ≥ σ 2 . Theorem 8. [Approximation Guarantee] Let K 0 = k(X, X) + σ 2 I be δ-diagonally dominant, s max > 0 be an upper bound on ∥s∥ ∞ , and suppose ∥y∥, δ satisfy the bounds of Lemmas 7 and 14, guaranteeing m-convexity and M -smoothness of the NMLL for some m > 0, M > 1/(1 -s max ) 2 . Let s OMP (r) be the r-sparse vector attained by OMP on the NMLL objective for r steps, and let s OPT (r) = arg max ∥s∥0=r,∥s∥∞≤smax L(ρ(s)) be the optimal r-sparse vector. Then for any 2r ≤ n,L (ρ(s OMP (r))) ≥ 1 -e -m/M L (ρ(s OPT (r))) ,where L(•) = L(•) -L(0) is normalized so that max s S L(s S ) ≥ 0 for any support S.Proof. The result is a direct consequence of meeting the m-convexity and M -smoothness conditions of Lemmas 7, 14 above, and the OMP approximation guarantee of Theorem 11 due to Elenberg et al. [21] . Note that the condition on 2r ≤ n comes from the RSC condition in Theorem 11 being required for subsets of size 2r. As we proved bounds for the m-convexity of the full Hessian of size n, r has to be smaller than n/2 for the assumptions of the theorem to hold. Regarding the upper bound s max on s, we note that the constraint is convex and therefore doesn't change the convexity property of the optimization problem.Further, note that max ρ S L(ρ S ) ≤ max ρ S∪i L(ρ S∪i ), since the additional non-zero element could stay at 0, if the marginal likelihood does not improve with ρ i increasing. That is, the subset selection problem is monotonic. As a consequence, we can normalize the MLL by L(•) = L(•) -L(0), which then only takes positive values for any s * S = arg max s \S =0 L(ρ(s S )), i.e. max s S L(s S ) ≥ 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 [55] 's Sparse Bayesian Learning (SBL) model, for which Ament and Gomes [8] show that greedily optimizing the associated NMLL is equivalent to stepwise regression in the limit of σ → 0, proving exact recovery guarantees.Our benchmark uses a modified version of the code from Andrade and Takeda [9] , available at https://github.com/andrade-stats/TrimmedMarginalLikelihoodGP under the GNU GPLv2 license.Model fitting runtimes Fig. 8 summarizes the fitting times for the different models on the different scenarios from Section 6.1. 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. 9 therefore provides a more compact view of the same data (aggregated across outlier types and outlier fractions), giving a better sense of the distribution of the fitting times. Unsurprisingly, the baselines that simply fit a single GP model ("vanilla", "winsorize", "power_transform") are substantially faster than any of the robust 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-t 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.For Hartmann6, we consider the standard domain of [0, 1] 6 .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: Rover The rover trajectory planning problem was originally proposed in Wang et al. [60] . The goal is to tune the trajectory of a rover in order to maximize the reward of its final trajectory. We use the same obstacle locations and trajectory parameterization as in Maus et al. [42] , with the main difference being that we parameterize the trajectory using 10 points, resulting in 20 tunable parameters. When a corruption occurs, the rover will stop at a uniformly random point along its trajectory, generally resulting in lower reward than the original trajectory.C ∈ [1e-2, 1e2], ϵ ∈ [1e-2,In addition to the results in Fig. 6 , we also include results when the corruption probability is 20% in Fig. 10 . We observe that the results are similar as with 10% corruption probability, but that the performance of several baselines regresses significantly.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 a 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). A limitation of Lemma 4 is that it only guarantees convexity for the sub-problem of optimizing the ρ i 's. In practice, we jointly optimize all GP hyper-parameters, including length scales. In this case, the theory guarantees the positive-definiteness of the submatrix of the Hessian corresponding to ρ(s), and we expect this to improve the quality of the results of the numerical optimization routines.Indeed, positive-definiteness is beneficial for quasi-Newton optimization algorithms like L-BFGS [38] , which restarts its approximation to the Hessian 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.To quantify the impact of this, we ran convergence analyses using the data from Fig. 1 , allowing all ρ to be optimized jointly with other hyper-parameters (length scale, kernel variance and noise We see that Relevance pursuit performs well in the case of constant outliers of value 100 and almost performs just as well as the oracle. Middle: Relevance pursuit performs the best followed by the Student-t likelihood from [41] in the case of U [-3 , 3] . No method performs as well as the oracle when the outlier probability is 20%. Right: Similarly to 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 are reported in Table D . 5 , and indicate that the optimizer terminates with a much better NLML using the convex parameterization with the same convergence tolerance.There are also settings in 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 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 we do not know to be outlier-free.We report extended empirical comparisons with Altamirano et al. [2] 's RCGP method, using their experimental setup and method implementation in GPFLow. Including GPFlow in our own benchmarking setup and compute resources proved difficult. To circumvent this, we wrote wrappers for both BoTorch's standard GP and RRP, which also accounts for any orthogonal implementation differences between the two frameworks, and ran the benchmarks locally on an M-series MacBook.See Tables 2 and 3 for the mean absolute error and negative log predictive density, respectively. The tables include the empirical mean and standard deviation over 20 replications on corrupted version of the following base data sets: 1) Synthetic, which is generated as a draw of a GP with a exponentiated quadratic kernel, and four data sets available on the UCI machine learning repository [32] , in particular, 2) Boston [27] , 3) Concrete [63] , 4) Energy [58] , and 5) Yacht [25] . The benchmark considers no corruptions ("No Outliers"), "Asymmetric Outliers", which are uniform in x 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 y. Any bold entry in the table signifies a results that is within one standard-error of the best result's one standard-error confidence bound. a variational (referred to as Huber-Projection) to get as close as possible to a direct comparison to Algikar and Mili [1] without access to a Matlab license.Tables 4 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 [25] and California Housing [46] datasets from the UCI database [32] , where 15% of the training data sets of the models were corrupted. Tables 6 and 7 below were generated in a 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 to 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 a subset of a-priori unknown outliers, while a dominant fraction of the data can be considered inliers that, once they are identified, can be used to train a model without additional treatment. Justification: Everything that is claimed in the abstract and introduction is 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 in 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 as long as it is clear that these goals are not attained by the paper.Justification: The research fully conforms to the NeurIPS Code of ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 


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

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 

2. While this work employs greedy algorithms for GPs, it is not the first to do so, so I would not consider it exceptionally novel.

3. The experimentation is adequate but less discussed.
 
Ref: Most Likely Heteroscedastic Gaussian Process Regression


Minor corrections: 

1. Consider adding references for these three categories on line 34. No reference is given for down weighting procedures even in related work.
2. Figure 1 and 2 are not discussed at all. What regression example? 
3. broken line 150 
4. unclosed braces figure 3
5. no need of repeating the abbreviation RIP on line 200
6. Broken line 264
7. Is the x axis of right figure in Figure 4 percentage of corruptions or probability?
8. Different notations need to be used for $\delta$ on lines 67 and 112.
questions: questions: 1. Why Matern kernel was chosen in one of the baselines?
2. 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}$?
3. What is $\mathbf{K}_{0}$ on line 218?
limitations: limitations: Limitations are discussed as I hoped they would be.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

Finally, 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.
questions: questions: - What's the difference between RGP-RP and RGP-RP* in the experiment section?
- What's the schedule $\mathcal{K}$ selected for the experiments? How does this affect the performance and fitting time of the method?
- While the Friedman and Hartmann6 test functions are well-known, the authors should provide a reference for these functions.
- 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?
- 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: limitations: No obvious negative societal impact
paper_summary: paper_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.
strength_weakness: strength_weakness: ### Theory
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 $\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.

### Related Work
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. 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.

- 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

### Experiments
As 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?
The 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).
From 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. 

### Other Weaknesses and Improvements
- Almost all plots are way too small and barely readable. If you reformat these, then I am willing to raise my "Presentation" score.
- Link to the proofs in the appendix after each theoretical result, otherwise they are hard to find.
- When you reference details in the appendix, link to the appropriate section (e.g. line 286).
- Citations in lines 32 to 34 for the methods that are referenced are missing
- Typos (l.211, l. 286, l. 362)
questions: 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).
- How does your approach compare to a vanilla GP on a dataset with no corruption? Does it "fail gracefully"?
- 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?
- How did you set the schedule to test data points for being outliers in the experiments (c.f. lines 182-184)?
limitations: 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.

Comments: title: A PyTorch implementation of Algikar et.al.'s projection statistics
Results for uniformly heavier-tailed noise
Results for Suggested Benchmarks
[time sensitive question]
comment: Thank you for providing the performance comparisons with additional models, especially on such short notice. I find the comparison to be very thorough and have adjusted my score accordingly.
Last, we provide the code of the main parts of our translation of the Matlab code of Algikar et.al., 2023 in the following, for you to be able to check its correctness.
```
def projection_statistics(H: Tensor) -> Tensor:
    """

    Args:
        H: (n x d)-dim Tensor, i.e. number of data points x dimensionality. NOTE:
            in the original code, this is taken to be X with an appended ones column.

    Returns:
        A (n)-dim Tensor of projection statistics.
    """
    dtype = H.dtype
    device = H.device
    m, n = H.shape
    M = torch.median(H, dim=0, keepdim=True).values  # (1 x d) row vector
    u = torch.zeros(m, n, dtype=dtype, device=device)  # i.e. (n x d) matrix
    v = torch.zeros(m, n, dtype=dtype, device=device)
    z = torch.zeros(m, 1, dtype=dtype, device=device)  # i.e. (n x 1) matrix
    P = torch.zeros(m, m, dtype=dtype, device=device)  # i.e. (n x n) matrix
    eps = 1e-6  # avoiding divide-by-zero issues
    for kk in range(m):
        u[kk, :] = H[kk, :] - M  # looping over data points
        v[kk, :] = u[kk, :] / max(torch.linalg.norm(u[kk, :]), eps)
        for ii in range(m):
            z[ii, :] = torch.dot(H[ii, :], v[kk, :])
        zmed = torch.median(z, dim=0, keepdim=True).values
        MAD = 1.4826 * (1 + (15 / (m))) * torch.median(torch.abs(z - zmed))
        for ii in range(m):
            P[kk, ii] = torch.abs(z[ii] - zmed) / max(MAD, eps)

    PS = torch.amax(P, dim=0)
    return PS


def _compute_projection_statistics_weights(H: Tensor, PS: Tensor) -> Tensor:
    """This computes the weights for the projection statistics, according to the procedure
    in the original Matlab code:
    https://github.com/apooja1/GP-Huber/blob/fee038963b471eb198d59b22d57b89e69f451d8c/Experiments/Friedman
    """
    niu = torch.sum(H != 0, dim=-1)
    cutoff_PS = torch.zeros_like(niu)
    for i in range(len(niu)):
        # for chi2, see this post for correspondence of ppf with invcdf:
        # https://stackoverflow.com/questions/53019080/chi2inv-in-python
        # the 0.975 was copied from the original matlab code:
        # https://github.com/apooja1/GP-Huber/blob/fee038963b471eb198d59b22d57b89e69f451d8c/Experiments/Friedman.m#L186
        cutoff_PS[i] = chi2.ppf(0.975, df=niu[i].item())

    weights = (cutoff_PS / PS.square()).clamp(max=1.0)
    return weights

def compute_projection_statistics_weights(X: Tensor) -> Tensor:
    # X: n x d
    n, d = X.shape
    H = torch.cat((X, torch.ones(n, 1)), dim=-1)  # n x (d + 1)
    PS = projection_statistics(H)  # (n x d) -> n
    return _compute_projection_statistics_weights(H, PS)
```

For the Huber loss, we use $\epsilon = 0.45$ and $b = 0.5$, identical to the setting of the [Github repository](https://github.com/apooja1/GP-Huber/blob/fee038963b471eb198d59b22d57b89e69f451d8c/basics/lik_huber.m#L191) of Algikar et.al., 2023.  All reported results were generated with 32 independent replications, reporting the mean and standard error.

An additional implementation detail: the `median` and `max` calls take the extra `dim=0` argument in PyTorch, as the Matlab function are defined to compute the row-wise maximum. Also, we added a small numerical constant `eps` in the code of `projection_statistics` to avoid division-by-zero errors. The rest of the code should be an almost verbatim translation. 

**A word of thanks** Thank you again for your valuable suggestions, we believe they make the paper stronger. Let us know if you have additional questions or if all your remaining concerns have been addressed.
**We highlight the main limitation** of Relevance Pursuit versus the additional baselines: if 100% of the data points are subject to heavy-tailed noise, the methods based on uniformly heavy-tailed likelihoods, including `Huber + Projection` will perform best, as the following results demonstrate, similar to the setup in Table 1 of Algikar et.al., 2023.

# 100% Laplace Noise 

## RMSE

| Data | Standard | Relevance Pursuit | Student-t | Laplace | Huber | Huber + Projection |
| --- | --- | --- | --- | --- | --- | --- |
neal + laplace | 1.51e+00 (1.06e-01) | 2.40e+00 (2.72e-01) | **1.16e+00 (9.18e-02)** | **1.12e+00 (9.90e-02)** | **1.12e+00 (9.80e-02)** | **1.12e+00 (9.80e-02)**
friedman5 + laplace | 1.39e+01 (5.95e-01) | 1.37e+01 (6.36e-01) | 8.33e+00 (3.69e-01) | **7.34e+00 (3.50e-01)** | **7.40e+00 (3.45e-01)** | **7.40e+00 (3.45e-01)**
friedman10 + laplace | 1.30e+01 (3.77e-01) | 1.27e+01 (3.99e-01) | 7.24e+00 (1.99e-01) | **6.07e+00 (1.91e-01)** | **6.26e+00 (1.71e-01)** | **6.26e+00 (1.71e-01)**
yacht_hydrodynamics + laplace | 2.57e+01 (1.17e+00) | 4.75e+01 (3.51e+00) | **1.64e+01 (3.18e-01)** | **1.61e+01 (2.52e-01)** | **1.62e+01 (2.60e-01)** | **1.67e+01 (3.55e-01)**
california_housing + laplace | 1.57e+00 (8.25e-02) | 2.06e+00 (1.44e-01) | **1.21e+00 (1.08e-02)** | **1.20e+00 (9.39e-03)** | **1.20e+00 (9.69e-03)** | 1.23e+00 (2.05e-02)

## NLP 

| Data | Standard | Relevance Pursuit | Student-t | Laplace | Huber | Huber + Projection |
| --- | --- | --- | --- | --- | --- | --- |
neal + laplace | 1.89e+00 (1.14e-01) | 1.38e+02 (3.98e+01) | **1.57e+00 (9.22e-02)** | **1.55e+00 (9.20e-02)** | **1.55e+00 (9.12e-02)** | **1.55e+00 (9.12e-02)**
friedman5 + laplace | 4.55e+00 (2.81e-02) | 4.56e+00 (1.88e-02) | 3.64e+00 (4.03e-02) | **3.46e+00 (4.45e-02)** | **3.48e+00 (4.31e-02)** | **3.48e+00 (4.31e-02)**
friedman10 + laplace | 4.58e+00 (1.25e-02) | 4.57e+00 (1.27e-02) | 3.60e+00 (2.53e-02) | **3.32e+00 (3.01e-02)** | **3.37e+00 (2.61e-02)** | **3.37e+00 (2.61e-02)**
yacht_hydrodynamics + laplace | **4.74e+00 (4.45e-02)** | 5.30e+00 (8.31e-02) | 4.91e+00 (8.87e-02) | 5.19e+00 (1.17e-01) | 5.12e+00 (1.04e-01) | 9.90e+00 (6.44e-01)
california_housing + laplace | **1.88e+00 (5.25e-02)** | 2.24e+00 (9.41e-02) | 2.92e+00 (6.01e-02) | 3.33e+00 (7.47e-02) | 3.27e+00 (5.84e-02) | 6.64e+00 (3.28e+00)
Thank you again for your suggestions. We present the additional results you suggested below and hope that this will allow you to fully recommend our work for acceptance.

**Summary** We added additional variational GP models with Laplace and Huber likelihoods, and translated the Matlab code of the "projection statistics" of Algikar et.al., 2023 ([Github](https://github.com/apooja1/GP-Huber/tree/main)) to PyTorch. We then combined the projection-statistics-based weighting of the Huber loss with a variational GP (referred to as `Huber-Projection` ) to get as close as possible to a direct comparison to Algikar et.al., 2023 without access to a Matlab license. 

## Additional Benchmark Results with 15% Corruptions (Negative Log Likelihood)

The following results include the Friedman 10 function used in Algikar et.al., 2023, and the UCI CA Housing data set, which - according to UCI - is a replacement for the deprecated UCI Boston Housing data, which was used in Algikar et.al., 2023.

In summary, for a sparse 15% of corrupted observations, Relevance Pursuit is able to almost uniformly achieve better predictive negative log likelihoods (NLL) for the datasets and outlier distributions we consider, the only exception being the CA Housing data, where the standard GP performs surprisingly well in terms of NLL, but not in terms of RMSE.

| Data | Standard | Relevance Pursuit | Student-t | Laplace | Huber | Huber + Projection |
| --- | --- | --- | --- | --- | --- | --- |
friedman5 + uniform | 2.07e+00 (4.88e-02) | **-1.35e+00 (7.22e-02)** | 4.67e-01 (3.00e-01) | 8.88e-01 (1.01e-01) | 8.29e-01 (9.14e-02) | 8.88e-01 (1.08e-01)
friedman5 + constant | 4.28e+00 (2.25e-02) | **-1.13e+00 (1.82e-02)** | 8.24e-01 (3.17e-01) | 2.28e+00 (2.63e-02) | 2.34e+00 (3.13e-02) | 2.34e+00 (3.13e-02)
friedman5 + student-t | 3.42e+00 (1.34e-01) | **-1.11e+00 (1.04e-01)** | 1.78e-02 (2.86e-01) | 1.55e+00 (1.14e-01) | 1.66e+00 (1.11e-01) | 1.59e+00 (1.10e-01)
friedman5 + laplace | 3.55e+00 (7.84e-02) | **-1.22e+00 (5.89e-02)** | 3.47e-01 (3.13e-01) | 1.75e+00 (8.46e-02) | 1.90e+00 (7.36e-02) | 1.90e+00 (7.36e-02)
friedman10 + uniform | 1.86e+00 (3.04e-02) | **-1.82e+00 (3.97e-02)** | 6.78e-02 (3.85e-01) | 9.20e-01 (1.99e-01) | 9.39e-01 (2.01e-01) | 7.65e-01 (1.80e-01)
friedman10 + constant | 4.21e+00 (1.39e-02) | **-1.27e+00 (1.10e-02)** | 8.85e-01 (3.79e-01) | 2.20e+00 (1.45e-02) | 2.23e+00 (1.54e-02) | 2.23e+00 (1.54e-02)
friedman10 + student-t | 3.34e+00 (1.24e-01) | **-1.40e+00 (1.17e-01)** | -6.80e-02 (3.65e-01) | 1.96e+00 (1.25e-01) | 1.95e+00 (1.23e-01) | 1.96e+00 (1.22e-01)
friedman10 + laplace | 3.46e+00 (6.68e-02) | **-1.44e+00 (5.91e-02)** | 4.89e-01 (3.94e-01) | 2.13e+00 (6.35e-02) | 2.21e+00 (2.01e-02) | 2.21e+00 (2.01e-02)
yacht_hydrodynamics + uniform | 3.45e+00 (1.42e-01) | **2.37e+00 (2.79e-01)** | 1.18e+02 (1.06e+01) | 7.07e+01 (5.17e+00) | 7.46e+01 (5.38e+00) | 1.89e+02 (1.64e+01)
yacht_hydrodynamics + constant | 5.29e+00 (1.70e-01) | **1.79e+00 (2.43e-01)** | 7.16e+01 (7.68e+00) | 2.32e+01 (1.59e+00) | 1.95e+01 (1.12e+00) | 3.89e+01 (3.46e+00)
yacht_hydrodynamics + student-t | 4.46e+00 (3.08e-01) | **2.42e+00 (3.38e-01)** | 1.32e+02 (1.05e+01) | 8.20e+01 (5.07e+00) | 8.13e+01 (5.17e+00) | 1.55e+02 (1.38e+01)
yacht_hydrodynamics + laplace | 4.35e+00 (2.88e-01) | **2.27e+00 (3.71e-01)** | 1.14e+02 (9.86e+00) | 6.85e+01 (4.72e+00) | 6.12e+01 (3.87e+00) | 1.17e+02 (8.89e+00)
california_housing + uniform | **3.14e+00 (1.44e-01)** | **2.92e+00 (2.06e-01)** | 5.85e+01 (9.59e-01) | 6.77e+01 (2.41e+00) | 6.69e+01 (1.94e+00) | 8.94e+01 (2.45e+01)
california_housing + constant | **3.57e+00 (2.26e-01)** | 4.51e+00 (2.11e-01) | 4.02e+01 (1.59e+00) | 1.83e+01 (6.58e-01) | 1.75e+01 (7.03e-01) | 2.18e+01 (4.19e+00)
california_housing + student-t | **1.77e+00 (6.65e-02)** | 3.15e+00 (1.67e-01) | 5.95e+01 (1.44e+00) | 5.20e+01 (2.19e+00) | 4.95e+01 (1.75e+00) | 6.90e+01 (2.02e+01)
california_housing + laplace | **1.61e+00 (4.21e-02)** | 3.51e+00 (1.93e-01) | 5.60e+01 (1.58e+00) | 4.23e+01 (1.72e+00) | 4.06e+01 (1.47e+00) | 5.62e+01 (1.66e+01)

## RMSE 

RP performs uniformly better w.r.t. RMSE, we only include CA Housing here due to space limitations.

| Data | Standard | Relevance Pursuit | Student-t | Laplace | Huber | Huber + Projection |
| --- | --- | --- | --- | --- | --- | --- |
california_housing + uniform | **7.10e-01 (7.58e-03)** | 7.39e-01 (1.75e-02) | 1.16e+00 (1.79e-03) | 1.17e+00 (3.04e-03) | 1.17e+00 (2.93e-03) | 1.18e+00 (6.17e-03)
california_housing + constant | 2.28e+00 (5.12e-02) | **6.35e-01 (4.38e-03)** | 1.17e+00 (2.39e-03) | 1.16e+00 (1.87e-03) | 1.16e+00 (1.88e-03) | 1.17e+00 (5.51e-03)
california_housing + student-t | 1.34e+00 (1.68e-01) | **6.56e-01 (5.46e-03)** | 1.18e+00 (3.56e-03) | 1.19e+00 (4.13e-03) | 1.18e+00 (3.83e-03) | 1.19e+00 (6.74e-03)
california_housing + laplace | 1.00e+00 (5.04e-02) | **6.51e-01 (4.71e-03)** | 1.18e+00 (3.41e-03) | 1.18e+00 (3.91e-03) | 1.18e+00 (3.67e-03) | 1.18e+00 (6.30e-03)
I don't intend to create last minute pressure of additional experimentation on you on the last day. I am glad to see that you offer the comparison with the requested models. I am okay with you adding them in the final submission of the paper. You got the references for the preprints right. I will change my score to 6.
Thank you for your response. I'll increase my score to 7.
We are glad we addressed your concerns in detail via our additional experiments which demonstrate the effectiveness of RGP.

The reviewer now suggests additional methods to compare against. At a minimum, we can discuss our work in the context of those papers, and promise to include these comparisons in the CR.   We hope that time permits us to run some additional simulations in the next day to meet your standards.

Before we can do this, can you please clarify the stated references in your follow-up response above?

Is "Kuss et. al 2006" Malte Kuß's PhD thesis https://pure.mpg.de/rest/items/item_1791134/component/file_3167621/content ?

Is "Algikar et.al., 2023" the unpublished manuscript of Pooja Algikar & Lamine Mili, Robust Gaussian Process Regression with Huber Likelihood?
(preprint at https://www.researchgate.net/publication/367280832_Robust_Gaussian_Process_Regression_with_Huber_Likelihood)
Thank you for addressing my concerns in detail. While the additional experimentation certainly demonstrates the effectiveness of the RGP, a more comprehensive comparison would benefit from including benchmarks that explicitly model heavy-tailed distributions, such as those employing Laplace (Kuss et. al., 2006) and Huber (Algikar et.al., 2023) likelihoods. Both methods address similar challenges – spatially uncorrelated gross outliers and input-independent noise – making them valuable benchmarks. Algikar et.al's method claims to not lose efficiency at Gaussian distribution using Huber loss and leverage weighting. My current score remains unchanged.
Thank you again for your detailed, actionable feedback, and consideration.
Thank you for your response. I believe the additional results you presented strengthen the paper. I have raised my score to reflect this.