Robust and Computation-Aware Gaussian Processes
Authors: Marshal Arijona Sinaga
Julien Martinelli
Samuel Kaski
Venue: NeurIPS 2025 poster ()
Keywords: Gaussian Processes
Bayesian optimization
Probabilist Numerics
Primary area: probabilistic_methods
TL;DR: RCaGP is a general-purpose Gaussian Process framework that combines robustness to outliers and approximation-aware uncertainty—two coupled limitations of sparse GPs that are critical to tackle for challenging tasks like high-throughput BayesOpt.

Abstract: Gaussian processes (GPs) are widely used for regression and optimization tasks such as Bayesian optimization (BO) due to their expressiveness and principled uncertainty estimates. However, in settings with large datasets corrupted by outliers, standard GPs and their sparse approximations struggle with computational tractability and robustness.
We introduce Robust Computation-aware Gaussian Process (RCaGP), a novel GP model that jointly addresses these challenges by combining a principled treatment of approximation-induced uncertainty with robust generalized Bayesian updating. The key insight is that robustness and approximation-awareness are not orthogonal but intertwined: approximations can exacerbate the impact of outliers, and mitigating one without the other is insufficient.
Unlike previous work that focuses narrowly on either robustness or approximation quality, RCaGP combines both in a principled and scalable framework, thus effectively managing both outliers and computational uncertainties introduced by approximations such as low-rank matrix multiplications.
Our model ensures more conservative and reliable uncertainty estimates, a property we rigorously demonstrate. Additionally, we establish a robustness property and show that the mean function is key to preserving it, motivating a tailored model selection scheme for robust mean functions. Empirical results confirm that solving these challenges jointly leads to superior performance across both clean and outlier-contaminated settings, both on regression and high-throughput Bayesian optimization benchmarks.

Robust and Computation-Aware Gaussian Processes
Marshal Sinaga
Julien Martinelli
ELLIS Institute Finland, Aalto University marshal marshal.sinagaaaaall sinagadaalto..fi
ELLIS Institute Finland, Aalto University
julien martinel 1 @aalta..o............
Samuel Kaski
ELLIS Institute Finland, Aalto University, University of Manchester samuel. kaski faalto.fi
Abstract
Gaussian processes (GPs) are widely used for regression and optimization tasks such Bayesian optimization (BO) due their expressiveness and principled uncertainty estimates. However, in settings with large datasets corrupted by outliers, standard GPs and their sparse approximations struggle with computational tractabil- ity and robustness. We introduce Robust Computation-awant aware Gaussian Process (RCaGP), a novel GP model that jointly addresses these challenges by combining a principled treatment of approximation-indueee uncertainty with robust generalized Bayesian updating. The key insight is that robustness and approximation caamaatioon-warensss.. are not orthogonal but intertwined: approximations can exacerbate the impact of outliers, and mitigating one without the other insufficient. Unlike previous work that focuses narrowly on either robustness or approximation quality, RCaGP combines both in a principled and scalable framework, thus effectively managing both outliers and computational uncertainties introduced by approximations such as low-rank matrix multiplications. Our model ensures more conservative and reliable uncertainty estimates, property we rigorously demonstrate. Addition- ally, we establish robustness property and show that the mean function is key preserving it, motivating tailored model selection scheme for robust mean functions. Empirical results confirm that solving these challenges jointly leads to superior performance across both clean and
settings, both on regression and high-throughput Bayesian optimization benchmarks.
Introduction
Gaussian Processes (GPs) are foundational tool in probabilistic machine learning, offering non- parametric modeling with principled uncertainty estimates [36]. Their use spans diverse domains such as time series forecasting, spatial modeling, regression on structured scientific data. One particularly impactful area where GPs have become integral is Bayesian Optimization (BO), cample-efficien framework for optimizing expensive black-box functions, where the surrogate model's 's predictive accuracy and uncertainty calibration directly influence decision-making [1]]. Many scientific and industrial tasks are high-throughput optimization problems in which numerous experiments simulations are performed, possibly in parallel. These can be effectively tackled by BO and include material design [15], drug discovery with extensive chemical space search [18], and fine-tuning gene expression in synthetic biology [6].
Yet, major challenges limit the reliability of GPs in both regression and BO settings: scalability observations. Sparse Variational to large datasets and robustness to corrupted or scalability via inducing point approximations, Gaussian Processes (SVGPs, [17, 41, 42]) address
yielding low-rank approximation that allows large-scale regression and high-throughput BO.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Acquisiion Function
XX
True functon
XCaGP
Octliees
Cati
SVCP
CatP
SVGP
Data
RCSVGP
EI-RCUOP
EI-RCSVGP
Global Optiman
Figure 1: Overview of proposed RCaGP against concurrent baselines (Left) SVGP fails to fit observed data contaminated by outliers, whereas RCaGP successfully fits high-density data regions while preserving higher variance due to the presence of outliers. (Middle- left) While enhancing robustness compared to SVGP, RCSVGP deviates more significantly from the true function than RCaGP. (Middle-right) Even if CaGP displays increased posterior variance near outliers, RCaGP provides superior posterior mean prediction. (Right) As result, the acquisition function landscape for RCaGP better prioritizes the true global optima.
1D example.
Despite their computational efficiency, significant challenges hinder the adoption SVGPs for high-throughput BO tasks, due the impact of surrogate posterior quality data acquisition. For starters, SVGP is not robust against outliers: upon contamination of the observed data, the posterior mean deviate significantly from the true latent function. Such deviation will misguide most acquisition strategies, resulting in suboptimal performance. Outliers are an almost unavoidable aspect of many real-world applications and datasets, specifically in high- tigh-throughput settings. They can arise from various factors, such as faulty measurements, malfunctioning sensors, or genetic anomalies. Next, SVGPs suffer from overconfident predictions, well-known pathology caused by the error incurred by the posterior approximation. Overconfident predictions occur in regions where the inducing points are far from observed data. Accurate uncertainty quantification essential for sequentially identifying the optimal data points in BO [45].
Contributions. To address these challenges, propose novel approach that combines robust- conjugate SVGP [2] with computation-aware GP [47]. This method, referred to as the Robust Computation-Awart GP (RCaGP), jointly tackles outlier robustness and approximation-induee uncertainty, two challenges that are deeply intertwined in sparse GP models. Computational approxi- mations such low-rank matrix factorizations one without addressing the other leads to suboptimal performance. RCaGP can be used as a surrogate in high-throughput BO, and its design is motivated by the need address these problems in tandem. This can be seen in Figure 1. In summary, this work
exacerbate the effects of outliers, and mitigating
1. (Methodological) introduces a novel GP model (RCaGP) that jointly addresses approximation- awareness and outlier robustness through a principled probabilistic framework, going beyond previous work that tackles these issues in isolation.
2. (Theoretical) establishes RCaGP's reliability through conservative uncertainty estimates and formal robustness property. Our analysis reveals the central role of the mean function in preserving robustness, motivating tailored model selection scheme.
3. (Empirical) demonstrates that jointly solving these two problems yields superior performance: RCaGP consistently outperforms existing state-of-the-art methods such as CaGP and RCSVGP on UCI regression tasks and BO settings, with and without outliers.
2 Background
Robust conjugate GP. The Gaussian noise assumption in GPs often leads to misspecification. e.g., due to outliers. To address this problem, recent work proposed robust conjugate GP (RCGP), which employs robust loss function through generalized Bayesian inference [2]. Consider noisy dataset D,, = 05)1, = (X We specify GP prior for the finite function U3)
evaluations f [/(xi). ... / /(x,)) p(f) = (f,,m, with m = [m{(x), ...xx))..... .., and
(.....,,(((x.))) vn(x)))] and
K = [k(x, (x)xx))]sssss Both m X functions, encoding prior beliefs about the function's structure and smoothness [36]. Conditioning
X - IR denote the mean and kernel

---

observations, the posterior predictive distribution is Gaussian with mean p. (x) and variance hs (x,x).
(x) = m(x) + k (K + T J w) (y mw). k. (x.x) = x(x,x) k2 (K + mosiso -
(2)
where ky = - [ E R The key difference with vanilla GPs is the introduction of weight function w, used by RCGP to down-weight the influence of potential outliers. Concretely, RCGP induces corrected mean and posterior variance by replacing (y ()) with the shrinkage term (y - mw). with mw m + o Vy (w²) and the identity matrix I with a diagonal matrix
Jw - diag Ciinte w w-2 for w = [mo(x), 31 ................................ w(xn; Mn)]".
Nevertheless, RCGP inherits the cubic inference complexity of vanilla GPs, prohibiting its usage large-scale datasets like those encountered in nigh-throughput BO. To remedy this issue, RCSVGP were introduced, enabling outlier-robust regression large-scale data [2]. This model introduces
inducing points u = {ur for r < n, with inducing locations z = {zr {zı such that u = f(4).. Given u, RCSVGP constructs an approximate posterior quu) (M En, with variational parameters Mu and En, obtained by maximizing the evidence lower bound (ELBO) criterion:
RCCSVr = (P) [log (y, f)] KL(g(u)lp(u)),
ELBO
where (y, f) denotes the robust loss function and ((f) = p(flu)q(u) is the approximate posterior over the function values conditioned on
While RCSVGP reduces the time complexity to (nr²), like other approximate GP methods, also exhibits overconfidence in regions where inducing points are far from the observations. This overconfidence can harm the exploration/explootttonn balance in BO, leading suboptimal results.
Computation-aware GP. [46, 47] proposed computation-aware GP (CaGP) to address the SVGP overconfidence. CaGP mitigates this issue by ensuring its posterior variance exceeds the vanilla GP variance. This corrected variance is obtained by adding computational uncertainty, associated with the approximation error in the representer weights v. (K + pusise-)) -yy -m). To obtain computational uncertainty, CaGP places a prior p(v.) Nvoo Co) and updates its posterior distribution (V,,: E,) using a probabilistic linear solver based on linear Gaussian identities:
V; C; (y m),
(4)
Σ K - C2
(5)
where C, 1SS K-1 is a low-rank approximation and S; E R the action matrix associated with the chosen approximation technique. In the case of inducing points, S; &(X,2 2.). Given a new data point x, the predictive distribution of f(x) is Gaussian, characterized by the mean (x(x) and variance (;(x, x), given by:
[x(x)= = m(() + kK Vs: 1x(x, x) = k((, x) kg - + E,k, = k(x, C,Kx
(6)
where K K R I and oz (x) = K E, denote the computational uncertainty. The posterior can be obtained O(?²t) time complexity. As n, CaGP recovers both the predictive mean and variance of the exact GP.
3 Robust and Computation-Aware Gaussian Processes
To tackle the intertwined challenges of outliers and approximation-inee indduced overconfidence, we introduce RCaGP a principled integration of robust-conjugate GPs with computation-awes ap- proximations. These issues must be addressed jointly: approximations can amplify outlier effects, while robustness alone does not correct the biases they introduce. RCaGP unifies both aspects to improve predictive reliability and decision- Section 3.1 details its composite inference scheme, Section 3.2 derives tailored optimization criterion, and Section 3.3 extends the framework to BO for joint model and query selection.
3

---

3.1 Robust Computation-aware Inference
Inspired by CaGP, we apply a probabilistic treatment to the representer weights & (K noise Jw > (y- mww) of RCGP in Equation (1). This enables us to quantify uncertainty in ap- proximating v-an aspect most approximate GPs ignore. RCaGP later incorporates this uncertainty into its predictive variance. We begin by placing a prior on v, i.e., & NV(vo 0, Σg K-1),
where K = K KKo noisee We then update this belief via the linear Gaussian identity [46], yielding Gaussian posterior with mean v; and variance Σ;:
Vi C, (y mw):
(8)
Σ = - C,,
(9)
where S. KS.) ( and S; exxi are low-rank approximations of K and the actions, respectively. The updates can be done sequentially [46, Algorithm I] or in batches [47, Algorithm S2].
Recent work generalized the action matrix S, E g by treating it as variable that can be optimized
during model selection [47]. Given new data point x X, the predictive distribution of f f(x) follows a normal distribution with mean ju: (x) and variance xx;
(x) = m(xx kx Vi.
(10)
((x,x))((xxx)) x) - k K, KK x(x,x) - kx Cik2.
(11)
In Equation (I1), RCaGP incorporates the computational uncertainty kg Σ, Kx, corresponding the variance of the representer weights Σ. resulting in a combined uncertainty (x, x). This computa- tional uncertainty guarantees that the predictive posterior variance of RCaGP is larger than RCGP, preventing RCaGP from overconfident predictions. Such behavior referred to as a conservative uncertainty estimate. The computational complexity of RCaGP can be found in Appendix B.
The robustness of RCaGP hinges on the weight function w, which we choose to incorporate into the representer weights V. To that end, we follow [2] and let:
w(x, y) = 3 (1 + (y - m(x))²/c)) 1/2
(12)
for 3 and O the learning rate and soft threshold.
The weight function assigns smaller values to observations y that deviate significantly from the mean prior m(x), treating them potential outliers. Lower weights reduce the influence of such points during inference, lcading to a more robust posterior. Additionally, the weight enters the noise term Jw diag([o" [ /2w 21), effectively inflating the noise for outliers and further limiting their impact. wxx, y) depends critically on the mean prior: for an outlier j. if m(x) is close to J. the weight increases. However, choosing an informative mean prior is difficult without domain knowledge.
3.2 Model hyperparameters optimization
We use the evidence lower bound (ELBO) as a loss function to optimize the kernel hyperpa- rameters A E TRP. Following [47], this enables RCaGP to scale to large datasets while avoid- ing overconfidence. The variational family is defined using the RCaGP posterior (f 0), i.e., (f:(F) = (ia (X), (x),,.(X (X,x))) E Rx parameterized by action S.. We replace the evidence term log pyy (0) with the robust loss log (y, from (2]. We then formulate ELBO
RRCaGP FLBO
Ey, log T (y (y,f)) (f])
(13)
This loss learns the hyperparameters if maximizing Ep(r/y on [log 1 (y, while minimizing the computational uncertainty associated with the approximation error. The derivation of the expected loss Ex(r) [log Toth (y, f) follows [2], while the derivation of the KI term [L[g(() [[p()) draws an analogy to CaGP. For the closed form and detailed derivation of ELBO, we refer the readers to Proposition C. We obtain the optimal hyperparameters 0" by minimizing the negative ELBO.
3.3 Joint model parameters and design selection strategy
To fully take advantage of RCaGP BO, we leverage the expected lower bound utility (EULBO) framework introduced by [29]. At acquisition time, instead of optimizing the surrogate's inducing

---

points, variational parameters, and hyperparameters, only then to find a design maximizing the acquisition function, global criterion maximized in an end-to-ent manner.
RRCGPP RRCGGP + (r) (log u(x, f: D:)]. EUULBO = ELBO
(14)
where u : X' x J + IR is the utility function. Here. we consider the expected improvement (EI) AF, reformulated as an expectation of the following utility:
(x, Pi) softplus(f(x) yf),
(15)
where i/ denotes the best evaluation observed so far and softplus x x log(1 exp(x)) This reformulation replaces the commonly used ReLU function with a softplus function, guaranteeing the utility function remains strictly positive whenever f(x) 2 y² Moreover, EULBO can be extended to support batch BO using Monte Carlo batch mode [5, 48]. Given set of candidates X ( the expected utility function corresponding the q g-improvement utility is:
Ellog uFI (X, x(JD,))
1 S
Σ
max softplus(r.), 1.....
(16)
with - (x.)) +c, Di and e, No, 1). Finally, when considering RCaGP, the action matrices S, are optimized, leading to the joint maximization problem
xi. S 02
argmax X1...X,,.S.00 S,0
RCCaGP EEULBO-
(17)
For efficiency, we impose the sparse block structure the action matrices S, E R following [47]. Specifically, we enforce each block to be column vector s, E R with k = n/i entries, so that the number of trainable parameters k X n (training data size). is worth noticing that unlike EULBO with SVGP, which maximizes variational parameters and inducing locations, our approach directly optimizes the action matrices S,
Maximizing the ELBO independently of the posterior-expected utility function can result in subopti- mal data acquisition decisions, the ELBO is primarily designed to model observed data [28]. In contrast, EULBO considers how the surrogate performs when selecting the next query, effectively guiding the solution of ELBO optimization towards high utility regions [29]. Our RCaGP model further balances such aggressive behavior by incorporating the combined uncertainty, preventing overconfidence during candidate query selection. Moreover, RCaGP enhances robustness in scenarios where function evaluations occasionally produce outliers.
4 Theoretical Analysis
This section presents theoretical analysis of RCaGP. Proposition 4.1 establishes its robustness outliers, while Proposition 4.2 links its uncertainty estimates to worst-case error, ensuring conservative predictions, thus addressing the challenges outlined in Section 1.
We demonstrate the robustness of RCaGP through the posterior influence function (PIF), similarly to [2, Proposition 3.2]. The PIF is a criterion proving robustness to misspecification in observation error [12, 27J. For this purpose, we define the dataset D = {(x), yj)} and the corresponding contaminated dataset D = (D \ (xm. ymm)) U (Xm= Vf. indexed by m E .......n.. The impact of vfr on inference is measured through PIF, expressed as Kullback- Leibler (KL) divergence between RCaGP's contaminated posterior p(f (..) and its uncontaminated counterpart p(f(D):
(I D) =
(18)
PIF views the KL-divergence as a function of lyf. ymml- In principle, we can consider any divergence that is not uniformly bounded. Here, we choose KL--divergence since it provides closed-form expression for multivariate Gaussian. A posterior robust sup [PIF D)I < oo [2]. We then show that RCaGP's PIF is bounded through the following proposition:
Proposition 4.1. Let f ~ GP(m, denote the RCaGP prior, and let E O,.... n represent the number of actions in RCaGP. Define constants C₁ and Cz. which are independent of yin For any i and assuming supx.y w(x, (1xx y) < oc, the PIF of RCaGP is given by:
PIF RCAGP (ym, D) = C4 (xm Vn.) y C.
(19)
Thus, if sup y wx y)² 0o, then RCaGP is robust since supys. |IIFCCGP(1//. D)1 < co.
5

---

The constraint s w(x,))
20 ensures no observation has infinite weight, and supxy V, w (x,x))) 3) < oo guarantees that w down-weights observations at least rate 1/y. The proof sketch shows that each term in Equation (18) is bounded. We leverage the pseudo-inverse property of positive semi-definite matrices and matrix norm bounds on the low-rank approximation to relate RCGP and CaGP. The full proof is in Appendix D. Crucially, the PIF of CaGP is unbounded, indicating a lack of robustness to outliers (see Appendix F.I).
Next, we show that RCaGP's uncertainty captures the worst-case error over all latent functions, paralleling [46, Theorem 2]. We assume the difference between the latent function g and the shrinkage term my lies an RKHS. This holds when , log nx(x v) < co. In the case of our robust weight function, the correction term is nonlinear rational function. Generally, such functions do not lie in the RKHS associated with common kernels (e.g., RBF, Matérn), unless specific and uncommon conditions are met-such as the RKHS being closed under composition with the given nonlinearity. A practical approach to ensure it belongs to the RKHS is to project the correction term onto the RKHS. Alternatively, explicitly accommodates this class of nonlinear transformations. Under this assumption, we establish the link between RCaGP's uncertainty estimates and the worst-case error through the following proposition:
could use design kernel whose RKHS
Proposition 4.2. Let :a((.-) = l.,)) + ōz((;-) be the combined uncertainty of RCaGP with zero-mean prior m. Then, for any new x E X we have that
sup (h(x) (x)) = Rx, x) + noise
(20)
[ 51
sup (1](x) (x)))) (x,,xx
(21)
I [1blegg $1
where (X,XKK(KX))))) is the RCGP's posterior and ()) K(, (gCC m (X) RCaGP's posterior mean for a function h E Hz = g(x) -mw (x) with latent function g X - IR
and the shrinkage function ma X x y + R.
The first equation in Proposition 4.2 shows that RCaGP's uncertainty ky(,,.) captures the worst-case error between the latent function g and the posterior mean (i.(.). Instead of pretending the model is fully correct, we allow a structured form of misspecification, namely, the part captured by mw- Proposition 4.2 then guarantees that any remaining discrepancy is still bounded by our variance. Even when classical Bayesian inference may fail to offer meaningful uncertainty due to model misspecification outliers, RCaGP maintains rigorous worst-case interpretation of its predictive variance, so long as the residual is sufficiently regular. This mirrors the CaGP assumption in spirit: just as CaGP analyzes worst-case error over functions in the RKHS, RCaGP does the same for the residual, which can be viewed as the portion of the signal that the robust update leaves uncorrected. The second shows that the computational uncertainty of captures the worst-case error between the RCGP and RCaGP posterior means, induced by approximation actions such as inducing points. The proof follows [46], leveraging [21, Proposition 3.9] and the RKHS associated with RCGP. Full details are in Appendix E. Lastly, the convergence of RCaGP in terms of mean function in RKHS norm is established in Appendix F.2.
5 Expert-guided robust mean prior definition
Our proposed RCaGP enforces outlier robustness through the weight function w, which as discussed in Section 3.1 and Proposition 4.2 criticltyyllyy depends on the prior mean m. Motivated by these insights, we propose defining a robust and informative mean prior using domain expert feedback.
We assume the expert identify subset of outliers in D, and provide corrections. Instead of discarding these outliers, we use them to inform the prior: removing them would increase evaluation cost, while RCaGP can still benefit from them. Since expert corrections are imperfect, model them probabilistically and define the mean via the inferred posterior.
Let {or J be binary outlier labels for dataset D,, {(x,, 2/ 1y= where o, 1 indicates an outlier. Each oj follows Bernoulli likelihood p(o; (6,) with latent probability 6; B(a,, 8y). We infer p(0 lo;) via Bayes' rule. Next, let ó {0.)a= and ý (2o J denote the identified outliers, and y 18 1 the expert-provided corrections. Each correction Do is drawn from N A o )or.)

---

P
Mo
Mo
10
1.0 0.0
1.5
griee nnfocmattiv
tidormativen
Figure 2: Expert-guided robust mean prior (Left) Graphical model of expert's 's feedback generative process for identified outliers. (Middle) User graphical model for expert outlier corrections. (Right) Toy example: with expert corrections, the expert gueedd prior better captures the true function than constant mean prior under contamination.
with latent to, (Po: To). Posterior inference yields p(ñoloo)- We define the mean prior by combining the expected labels and corrected values using a kernel-weighted average:
m(x)
16,(, 15y) By) sx)ūū)((xxxx)) Xy),
(22)
This corresponds averaging over the á nearest identified outliers, weighted by kernel and the expected confidence in each correction. In practice, we average over all identified outliers and use only posterior expectations, prioritizing computational efficiency over full Bayesian inference. The full generative process and toy illustration appear in Figure (left and middle). Additional inference details provided in Appendix G.
Related Work
Outlier-robust Gaussian processes fall into three main categories. The first replaces Gaussian noise with heavy-tailed alternatives e.g.. Student-t[20, 35], Huber[1], Laplace [23]. data-dependent [13], mixture models [8, 25]. These approaches break conjugacy and require approximate inference. The second class removes suspected outliers before applying standard inference [24, 32, 4], but this often unreliable and computationally expensive in high dimensions. The third approach, RCGP [2], builds on generalized Bayesian inference [7, 19, 22] to down-weight outliers without discarding data. We adopt RCGP the foundation for RCaGP due its strong performance in regression with outliers, its preserved conjugacy, and greater computational efficiency-enabllee enabling more accurate posterior inference and better predictive performance.
e..g..
SVGPs are the standard approach for scaling Gaussian processes to large-scale Bayesian optimiza- tion [14, 38, 44, 43]. [31] improved inducing point selection by optimizing their weighted Gram matrix via determinantal point process. [29] formulated acquisition as expected utility, enabling joint optimization of queries, variational parameters, hyperparameters, and inducing points in end-to-end fashion. [26] proposed an online conditioning strategy that updates SVGPs with new data without reoptimizing variational parameters, enabling look-ahead acquisition in BO.
Computation-aware Gaussian processes are novel class GPs that accounts for computational
uncertainty, due to approximations in matrix-vector multiplications [46]. These methods fall under the domain of probabilistic numerics, which quantify additional uncertainties arising from approximation errors [33, 40, 34]. In that respect, [47] introduced batch update algorithm for CaGP and proposed an ELBO objective to facilitate model selection. [16] extended CaGPs by proposing a calibrated probabilistic linear solver, resulting in calibrated computation-awee -wware GP model.
7 Experimental results
We evaluate the performance of our proposed RCaGP against concurrent baselines range of regression datasets, including a real-world dataset, followed by several tigh-throughhut Bayesian Optimization tasks. We conclude by presenting insights gained after carrying out dedicated ablation studies. Our implementation is available at https:: //t///gttttbbb//////// com/MarshaaAArijnnaa /CaaPP.
GP baselines. We compare RCaGP against several baseline methods: RCSVGP [2], SVGP [17]. and CaGP [47]. For the UCI regression benchmarks, we further include comparisons with SVGP models using Student-t likelihood [20] and SVGP with relevance pursuit (RRP), adapted from [3]
7

---

and formally derived in Section F.3. We set the number of actions (for RCaGP and CaGP) or inducing points (for SVGPs and RCSVGP) 25. All baselines employ a Matérn-5/2 kernel and are implemented in GPyTorch [10]. Model hyperparameters are optimized by maximizing the BLBO (Equation 13).
Outlier contamination protocol. In our experiments, we follow the settings described by [2]. Unless stated otherwise in dedicated ablation studies, for regression datasets, we sample uniformly at random 10% of the training dataset input-output pairs (x,, Mi), and replace the yi 's by asymmetric outliers, i.e., via subtraction of noise sampled from a uniform distribution ((3 96), with ō being the standard deviation of the original observations. For BO tasks, while running the optimization loop. each evaluation has a 25% chance of returning an asymmetric outlier, obtained by addition of noise sampled from a uniform distribution 2(155, 25), with ō being is the estimated standard deviation of the normalized evaluation function. Finally, for the weight function w (Equation 12), we set the soft threshold c = Qn - m|, c) with 0.2, with On the (1 -uquunie of ly -m/.
Bayesian Optimization settings. As described during Section 3.3, we integrate the baselines into EULBO, building on the implementation of [29]. We initialize the optimization process for all baselines with 250 data points sampled uniformly across the search space. All surrogates employ expected improvement (EI) as the acquisition function. Standard BO is applied across all tasks, while trust-region BO (TuRBO [9]) is additionally employed for high- nigh-dimensional tasks.
7.1 Regression on UCI datasets
We evaluate GP regression on four UCI datasets with asymmetric outliers (details in Section H) results are shown in the left part of Table 1. RCaGP delivers the overall best performance, combining strong predictive accuracy and well-calibrated uncertainty with acceptable runtimes, and being only surpassed in the Yacht dataset by RRP, much slower baseline. Notably, it significantly outperforms both CaGP and RCSVGP in this setting, indicating that the integration of robustness and is effective.
approximation-awaareee
To test robustness across contamination types. we further consider two additional outlier scenarios, uniform and focused outliers, inspired by [2] (definitions in Section H.3). The corresponding results are shown in Table SI and Table S2. RCaGP remains the top performer overall, particularly in uniform outliers for MAE, and focused outliers for NLL. These findings support its robustness across contamination types. More importantly, because outlier presence is not always known in advance, we also evaluate performance in the absence of outliers (right side of Table 1). RCaGP remains superior, except for the Energy dataset, where CaGP marginally outperforms RCaGP. This suggests that even in the absence of outliers, the weight function used by RCaGP (Fquation 12) proves useful. An alternative explanation would be that these datasets already contain outliers from the start.
Overall, these results demonstrate that RCaGP is more than the sum of its parts: the combination of RCGP and CaGP is not merely additive, but synergistic, consistently outperforming either component alone, both in outlier-rich and clean regimes.
Table UCT Regression datasets results with asymmetric outliers and without outliers. Average test set mean absolute error, negative ligelikelihood, and clock-time (in seconds), with std, for 20 train-test splits. Bolded results refer to the best baseline. Lower is better.
Asymmetric ontllers
No outllers
Boston
Energy
Yacht
Parkimsons
Boston
Energy
Yacht
Parkinsons
SVGP
MAE
1 10062
0.607 0.066
0.837 a 0.103
0.766
1 0.048
0-107
0.672
0.701 0.078
NLL
1,442 1 0.037
1,298 0.036
1,400 a 0,086
1,475
0.062
1,335 a 0.036
1,207
0.017
1.358
0.042
1 0.063
chokk-ime
1,400
+ 0.260
2.770 e 0,070
0.960 = 0.070
0.730 0.170
1,570
+ 370
2.930
0.310
0.950
0.030
0.520 = 0,030
CROP
MAE
0,738
1 10059
0,562 0.035
0.775 1 0.078
0.798
0.0S
0,158
1 1001
0.334
0.024
0.313 0.01
0.675 100TT
NLL
0.0357
0.027
1.395
1.278 0.040
1.106
1.263 0.042
1.362
clock time
2 =
0.720 a
= 0.040
3,710
0.730
0,840
RCSVGP
MAR
0.332
0.366 176600.00 0.30
0.477 = (.043
0.627
0.374
0.05
0.612 0.034
0.646 0.076
0.64K +000T
NLL
1.298 + 0.038
1 = 0.034
1.241 = 0.035
1.346
0.070
1,332
+ 0.048
1.324
+ 0.033
1.322 = 0.947
1,403 + 0.07%
time
1,500
0.000
2,860 # 0.060
1.030
0.750
0.020
1,560
0.040
2,970
0.120
1.010=0.030
0.510
RCSGP
MAE
= 0.846
0,380 0J013
0.435 =
0.586
0,084
0.42N
01303+
01,43 0.014
0,538
NLL
1,272 a: 0.036
1.162 0.012
1,253 a 0.043
1.317
0.067
1,261 a 0.049
1,160 0.008
1.251 0.040
1
chokkfliee
2,450
01000
3,790 + 0.150
0.870 = 0.060
1.250
+ 0.000
2.580
+ 150
3,780 e 0.190
0.910 = 0.070
0.920 +: 0.020
Student-
MAE
0.962
0.007
1.015 2 0.028
0,907
1 0.061
1,025
0.079
0.8%6
1 1.001
0,970
0.025
0.815 0.071
0.965 10071
NLL
1.564 0.067
0.085
1.572
0.026
1.53%
1.627
0.002
cloek-time
1,210 0.240
1,800 0.400
1,350 e 0.260
0.990
0.000
1,330 + 0.050
1.330
D.000
03890
0.150
0.730 0.020
RRF
MAE
0.550 1 0.050
0.597 2 0.112
0,356 2 0.051
0.894
0.100
0.815
1.15500.110
0.115
0,708
2 0.035
00.00
NLL
1,147 = 0.050
1.240 = 0.117
0.896 a 0.033
1.306
0.12
1.527
+ 194
0:214
1.179
= 0.045/
1,884 0.177
clock
time
5.120=00.100
3,840 0.090
7,200 0.160 3,010
0.080
5.260
+ 0.120 3.730
0.100
6.980
0.000
8

---

1.0
mWI M
0.8
0.6
DJA
0.4
Fre
0.2
10.00
11d0
1380
13.00
14.00
16.00
10.00
11.00
12.00
13,00
18.00
18ve
10:00
Figure 3: Left: Predictive mean of various approximate GP models on the DJIA index with constant mean prior. Outliers affect the baselines, whereas RCaGP maintains robustness. Right: learned weight functions wxx, y) for RCaGP and RCSVGP, illustrating their influence on the regression.
7.2 Twitter Flash Crash
We illustrate the woaldworld applicability of RCaGP by evaluating its performance and other baselines on the Dow Jones Industrial Average (DJIA) index for April 17, 2013. On this day. the Associated Press Twitter account was hacked, posting false report of an explosion the White House and an injury to the U.S. President. The resulting panic triggered sudden market sell-off, followed by rapid rebound, creating a brief interval during which the DJIA did not reflect the true economic state of the U.S. stock market. This scenario provides unique stress test for evaluating model robustness under sudden, anomalous market conditions. Observations during this period can reasonably be treated as outliers. Figure 3 (left) shows the DJIA index alongside the predictive means of RCaGP and other baselines. The predictive means of SVGP, RCSVGP, and CaGP exhibit noticeable deviations before or after PM, while RCaGP yields more reliable estimate, with only minor discrepancies around 2 PM. These results suggest that insufficient robustness or computational uncertainty in the baseline models contributes to their degraded performance under sudden market shocks. In the right panel of Figure 3, we plot the weight functions of RCaGP and RCSVGP. Although their weight functions exhibit similar profiles, they yield markedly different mean predictions. This supports our claim that robustness alone, without an appropriate treatment of computational uncertainty, can degrade the performance of approximate Gaussian process models.
7.3 High-Throughput Bayesian Optimization
We consider 4 tasks: Hartmann6DD, Lunar12D, Rover60D and Lasso-DNA180D (Section H.2). Figure 4 reports the best value found during BO, averaged over 20 trials, with =1 standard deviation. To ensure fair comparisons despite outlier contamination, we report the uncontaminated mean best value. In the top row (columns 1-4), using EULBO-EI and under asymmetric outliers, RCaGP clearly outperforms all baselines on Lasso DNA 180D, and ranks first or on par with CaGP on Hartmann 6D, Lunar 12D, and Rover60D, SVGP and RCSVGP display the same performance. The bottom row (columns 1-4) shows results with Determinantal Point Process (DPP)-based inducing point selection [31], still in the case of asymmetric outliers. RCaGP remains the top performer on DNA and performs similarly to CaGP on other tasks, with RCSVGP slightly ahead on Hartmann. Substituting El with Thompson Sampling in DPP-BO leads to equal performance across baselines (Figure $1).
Next, Figure S2 presents results with TurBO, a BO variant tailored to high dimensions [9]. RCaGP exhibits strong average performance, and notably, it is the only method that never ranks last in statistically significant manner. Lastly, as one might not know beforehand whether outliers will affect the trial, column 5 in Figure 4 shows results for EULBO El on Hartmann and Lunar without outliers. RCaGP maintains a clear advantage on Hartmann, while all baselines perform similarly on Lunar.
Together, these results highlight RCaGP's versatility and robustness across inducing point allocation strategies (EULBO, DPP), acquisition functions (EI, TS), BO behavior (global BO, trust-region BO), and noise regimes (with or without outliers). Its consistent superiority supports the central claim: combining robustness from RCGP with the uncertainty calibration of CaGP leads to a synergistic improvement, not just an additive one, across wide range of BO challenges.

---

Figure 4: High-throughput Bayesian Optimization task. Each panel shows the best value found each iteration found so far, averaged across 20 repetitions + I std. Columns 1-4 report results under asymmetric outliers contamination, with the 1st row using the BULBO-EI acquisition function, the 2nd row using DPPBO-EL. Column 5 features results without outliers using EULBO-EI.
Figure 5: Ablation studies. (Left) Varying in weight function (Equation 12) UCI regression datasets. MAE and NLL have been normalized for each dataset; the metries displayed represent the average. (Right) BO results for RCaGP and RCSVGP using expert-driven prior mean function in w. Results are averaged over across 20 repetitions, with standard deviation being shown.
7.4 Ablation studies
Effect of the soft threshold in the weight function. We evaluate RCaGP's sensitivity to the threshold e = Qn ml, <), defined the (1 - -q-anntlee of ly m| and controlling how many values are treated as outliers, by varying e. Figure 5 (left) shows MAE and NLL averaged over 4 UCI datasets (normalized; per-dataset results in Figure S3). RCaGP consistently achieves lower етгог across all with best results at 0.15, suggesting treating roughly 15% of the data as potential outliers is optimal. This aligns well with the naturally present and 10% injected asymmetric outliers.
Effect of the mean m the weight function. Given the central role of the mean function m in w (Equation 12), compare constant prior to the expert-driven one from Section 5, assuming perfect outlier correction. This yields four variants: RCaGP/RCSVGP with either constant (const.) or expert (exp.) mean. As shown in Figure (right), the informed prior greatly boosts RCSVGP across BO tasks, and modestly accelerates RCaGP on Lunar I.D. See Appendix G for further analysis.
8 Conclusion
This paper introduced Robust Computation-Aware Gaussian Processes (RCaGP), a principled frame- work unifying robust inference and approximation-awarr modeling coherently rather than tackling these challenges independently. Our theoretical and empirical results show that RCaGP yields conser- vative uncertainty estimates, provably robust outliers, and consistently outperforms methods addressing only one challenge. Additionally, we proposed an sxpert-guided prior mean to further en- hance robustness in practice. RCaGP represents a meaningful step forward probabilistic numerics, enabling robust and trustworthy inference for complex, large-scale regression and optimization tasks.
Limitations. While RCaGP demonstrates strong performance across diverse regression and BO tasks, it inherits certain practical limitations. Like CaGP, it requires full batching, leading to increased GPU memory demands compared to fully stochastic approaches like SVGP. Its performance also depends on manually selecting the projection dimensionality, which cannot be tuned via standard GP model selection techniques. Lastly, broader evaluations across more complex domains remain an avenue for future work.
10

---

Acknowledgments. This research was supported by the Research Council of Finland (flagship programme: Finnish Center for Artificial Intelligence, FCAI grants 358958. 345604, and 341763), and the UKRI Turing Al World-Leading Researcher Fellowship. EP/W002973/1. We also acknowledge the computational resources provided by the Aalto Science-IT Project.
References
[1] Pooja Algikar and Lamine Mili. Robust gaussian process regression with huber likelihood and projection pursuit. In Joint European Conference on Machine Learning and Knowledge Discovery Databases, 2025.
[2] Matias Altamirano, François- Xavier Briol, and Jeremias Knoblauch. Robust and conjugate gaussian process regression. In Proceedinggs of the 41st International Conference on Machine Learning 2024.
[3] Sebastian Ament, Elizabeth Santorella, David Eriksson, Ben Letham, Maximilian Balandat, and Eytan Bakshy. Robust gaussian processes via relevance pursuit. Advances in Neural Information Processing Systems, 37:61700-6111 734, 2024.
[4] Daniel Andrade and Akiko Takeda. Robust gaussian process regression with the trimmed marginal likelihood. In Uncertainty Artificial Intelligence, 2023.
[5] Maximilian Balandat, Brian Karrer, Daniel Jiang, Samuel Daulton, Ben Letham, Andrew G Wil- son, and Eytan Bakshy. Botorch: A framework for efficient monte-carlo bayesian optimization. In Advances in neural information processing systems, 2020.
[6] David Belanger, Suhani Vora, Zelda Mariet, Ramya Deshpande, David Dohan, Christof Anger- mueller, Kevin Murphy, Olivier Chapelle, and Lucy Colwell. Biological sequence design using batched bayesian optimization. In NeurIPS 2019 Workshop on Machine Learning and the Physical Sciences, 2019.
[7] Pier Giovanni Bissiri, Chris C Holmes, and Stephen G Walker. A general framework for updating belief distributions. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2016.
[8] Atefeh Daemi, Yousef Alipouri, and Biao Huang. Identification of robust gaussian process regression with noisy input using em algorithm. Chemometrics and Intelligent Laboratory Systems, 2019.
[9] David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local bayesian optimization. Advances in neural information processing systems, 32, 2019.
[10] Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. Advances in neural information processing systems, 2018.
[11] Roman Gareet.. Bayesian optimization. Cambridge University Press, 2023.
[12] Abhik Ghosh and Ayanendranath Basu. Robust bayes estimation using the density power divergence. Annals of the Institute of Statistical Mathematics, 2016.
[13] Paul Goldberg, Christopher Williams, and Christopher Bishop. Regression with input-dependent noise: gaussian process treatment. In Advances in neural information processing systems, 1997.
[14] Ryan-Rhy- Griffiths and José Miguel Hernández-Lobato. Constrained bayesian optimization for automatic chemical design using variational autoencoders. Chemical science. 2020.
[15] Trevor Hastings, Mrinalini Mulukutla, Danial Khatamsaz, Daniel Salas, Wenle Xu, Daniel Lewis, Nicole Person, Matthew Skokan, Braden Miller, James Paramore, et al. An interoperable multi objective batch bayesian optimization framework for high throughput materials discovery. arXiv preprint arXiv:2405. 240.:0000500 2024.

---

[16] Disha Hegde, Mohamed Adil, and Jon Cockayne. Calibrated computation-await gaussian processes. arXiv preprint arXiv:2410.08796 2024.
[17] James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data, In Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence, 2013,
[18] José Miguel Hemúndeez-oobtoo Edward Pyzer-Knapp. Alan Aspuru- Guzik, and Ryan P Adams. Distributed thompson sampling for large-scale accelerated exploration of chemical space. In NIPS Workshop on Bayesian Optimization, 2016.
[19] Jack Jewson, Jim Q Smith, and Chris Holmes. Principles of bayesian inference using general divergence criteria. Entropy, 2018.
[20] Pasi Jylänki, Jarno Vanhatalo, and Aki Vehtari. Robust gaussian process regression with student-t likelihood. Journal of Machine Learning Research, 2011.
[21] Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K Sriperumbudur. Gaus- sian processes and kernel methods: A review connections and equivalences. arXiv preprint arXiv:1807.02582, 2018.
[22] Jeremias Knoblauch, Jack Jewson, and Theodoros Damoulas. An optinization-centicc view on Bayes' rule: Reviewing and generalizing variational inference. Journal of Machine Learning Research, 2022.
[23] Malte Kuss. Gaussian process models for robust regression, classification, and reinforcement learning. PhD thesis, echnische Universität Darmstadt Darmstadt, Germany, 2006.
[24] Zhao-Zhou Li, Lu Li, and Zhengyi Shao. Robust gaussian process regression based on iterative trimming. Astronomy and Computing, 2021.
[25] Yifan Lu, Jiayi Ma, Leyuan Fang, Xin Tian, and Junjun Jiang. Robust and scalable gaussian process regression and its applications. In Proceedings of the IEEF/CVF Conference on Computer Vision and Pattern Recognition, 2023.
[26] Wesley Maddox, Samuel Stanton, and Andrew Wilson. Conditioning sparse variational gaussian processes for online gecision-making In Advances in Neural Information Processing Systems, 2021.
[27] Takuo Matsubara, Jeremias Knoblauch, François- Xavier Briol, and Chris J Oates. Robust generalised bayesian inference for intractable likelihoods. Journal of the Royal Statistical Society Series B: Statistical Methodology, 2022.
[28] Alexander G de G Matthews, James Hensman, Richard Turner, and Zoubin Ghahramani. On sparse variational methods and the kullback-leibler divergence between stochastic processes. In Artificial Intelligence and Statistics, 2016.
[29] Natalie Maus, Kyurae Kim, David Eriksson, Geoff Pleiss, John Patrick Cunningham, and Jacob R. Gardner. Approximation aware bayesian optimization. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024.
aware
[30] Natalie Maus, Kaiwen Wu, David Eriksson, and Jacob Gardner. Discovering many diverse solutions with bayesian optimization, In Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, 2023.
[31] Henry B Moss, Sebastian W Ober, and Victor Picheny. Inducing point allocation for sparse gaussian processes in high-throughput bayesian optimisation. International Conference on Artificial Intelligence and Statistics, 2023.
[32] Chiwoo Park. David Borth, Nicholas S Wilson, Chad N Hunter, and Fritz J Friedersdorf. Robust gaussian process regression with a bias model. Pattern Recognition, 2022.
[33] Marvin Pförtner, Ingo Steinwart, Philipp Hennig, and Jonathan Wenger. Physies-informed gaussian process regression generalizes linear pde solvers. arXiv preprint 12Xiv:2212. 12474, 2022.
12

---

[34] Marvin Pförtner, Jonathan Wenger, Jon Cockayne, and Philipp Hennig, Computation-awe aware kalman filtering and smoothing. In Proceedings of The 28th International Conference on Artificial Intelligence and Statistics, 2025.
[35] Rishik Ranjan, Biao Huang, and Alireza Fatehi. Robust gaussian process modeling using algorithm. Journal of Process Control, 2016.
[36] C. Rasmussen and C. Williams. Gaussian Processes for Machine Learning. MIT Press, 2006.
[37] Bernhard Schölkopf, Ralf Herbrich, and Alex Smola. generalized representer theorem. In International conference on computational learning theory, 2001.
[38] Samuel Stanton, Wesley Maddox, Nate Gruver, Phillip Maffettone, Emily Delaney, Peyton Greenside, and Andrew Gordon Wilson. Accelerating bayesian optimization for biological se- quence design with denoising autoencoders. In International Conference on Machine Learning, 2022.
[39] Sonja Surjanovic and Derek Bingham. Virtual library of simulation experiments: Test functions and datasets. https://wvi /ww.. sfu. saaa/sssrrssssseaaaaaaaseesss ano/ hwtnnn//ptiiiiii aptimiiaation html, 2020.
[40] Lukas Tatzel, Jonathan Wenger, Frank Schneider, and Philipp Hennig. Accelerating generalized linear models by trading off computation for uncertainty. CoRR, 2023.
[41] Michalis Titsias Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics, 2009,
[42] Michalis Titsias and Miguel AzzaroGGredille Doubly stochastic variational bayes for non- conjugate inference. In International conference on machine learning, 2014.
[43] Austin Tripp. Erik Daxberger, and José Miguel HernándezzLobato. Sample-efficient optimiza- tion in the latent space of deep generative models via weighted retraining. In Advances in Neural Information Processing Systems, 2020.
[44] Sattar Vakili, Henry Moss, Artem Artemev, Vincent Dutordoir, and Victor Picheny. Scalable thompson sampling using sparse gaussian process models. In Advances in neurul information processing systems, 2021.
[45] Zi Wang, Clement Gehring, Pushmeet Kohli, and Stefanie Jegelka. Batched scage bayesian optimization -ighdiimessionll spaces. In International Conference on Artificial Intelligence and Statistics, 2018.
[46] Jonathan Wenger, Geoff Pleiss, Marvin Pförtner, Philipp Hennig, and John P Cunningham. Pos- terior and computational uncertainty in gaussian processes. In Advances Neural Information Processing Systems, 2022.
[47] Jonathan Wenger, Kaiwen Wu, Philipp Hennig, Jacob R. Gardner, Geoff Pleiss, and John Patrick Cunningham. Computation-awarr gaussian processes: Model selection and linear-time inference. In The eiihth Annual Conference Neural Information Processing Systems, 2024.
[48] James Wilson, Frank Hutter, and Marc Deisenroth. Maximizing acquisition functions for bayesian optimization. In Advances in neural information processing systems, 2018.
[49] Kenan Šehic, Alexandre Gramfort, Joseph Salmon, and Luigi Nardi. Lassobench: A high- dimensional hyperparameter optimization benchmark suite for lasso. In Proceedings of the 1st International Conference on Automated Machine Learning (AutoML), 2022.
13

---

Supplementary Material Robust and Computation-Aware Gaussian Processes
The appendix is organized as follows:
Appendix A gives the results of additional experiments mentioned in the main text. These include experiments on UCI regression datasets with uniform and focused outliers (Table $1 and Table S2), BO results for DPP-BO-TE OO-TS (Figure S1), TuRBO (Figure S2), and dataset- specific results for the ablation study conducted Section 7.4, where c varied.
Appendix B gives the computational complexity of RCaGP.
Appendix C formally derives the ELBO for RCaGP introduced in Section 3.2.
Appendix D contains the proof for the robustness property presented in Proposition 4.1.
Appendix E proves Proposition 4.2, showing that RCaGP's uncertainty estimates capture the worse-case error over all latent functions.
Appendix F contains additional theoretical results: proposition showing the convergence mean of RCaGP RKHS norm (Section F.2), and proposition showing the lack of robustness of CaGP (Section F.1).
Appendix G provides further details regarding the expert-driven prior mean introduced in Section 5 and employed during the ablation study in Section 7.4.
Appendix H provides additional details about the UCI regression datasets (Section H.1), the test functions employed in BO (Section H.2), and the outlier-contamination protocols studied in this work (Section H.3).
Appendix provides tables that report the hyperparmeters used for UCI regression and high- throughput BO experiments (Section I.2) and the description of the computing resources used for our experiments (Section I.1).
A Additional experiments
The results of the UCI regression experiments in the presence of uniform and focused outliers are presented in Table S1 and Table S2, respectively. The MAE and NLL metrics indicate that RCaGP consistently outperforms the baseline methods across most datasets. These findings demonstrate robustness and versatility of RCaGP in handling different types of outliers. We further conducted experiments using the radial basis function (RBF) kernel on both clean and contaminated datasets. While RCaGP achieves competitive MAE and NLL performance compared to RCSVGP on clean data, substantially outperforms all baselines on contaminated data. However, this performance improvement comes the cost of increased computational time compared to the baselines.
The results of Bayesian optimization (BO) experiments using the Thompson Sampling (TS) ac- quisition function (AF) and Trust-Region BO (TuRBO) are presented in Figure S1 and Figure S2, respectively. In experiments utilizing the TF acquisition function, no model demonstrates a statisti- cally significant performance advantage. Furthermore, the best values identified by these models are generally lower than those obtained with the BO-EI method. We hypothesize that this performance degradation arises from the reliance of the TS approach on posterior sampling, which may reduce the effectiveness of the robustness properties inherent in RCaGP. Subsequent BO experiments using the TuRBO framework show improved performance across all models and except for the Lunar task. In this case, TuRBO yields lower results for all models. This suggests that the presence of outliers in the Lunar task may impede TuRBO's ability to effectively explore and identify the global optimum.
The MAE and NLL plots for varying values of in the weight function are shown in Figure S3. Across all datasets, RCaGP consistently yields lower MAE and NLL values than RCSVGP, indicating superior predictive performance. Furthermore, both metrics improve when is set to 0.1 or 0.,5, but begin deteriorate once e exceeds 0.2. This trend suggests that aligning with the true proportion of outliers in the data enhances model robustness and overall performance.
14

---

Table S1: UCI Regression datasets results in the presence of uniform outliers with 0.1. Average test set mean absolute error, negative log-likelihood, and clock-time (in seconds), with 1 std, for 20 train-test splits. Bolded results refer to the best baseline. Lower is better.
Boston
Energy
Yacht
Parkinsons
SVGP
MAE
0.799 = 0.054
0.918 = 0.070 0.899 0.099
0.837 = 0.093
NLL
1.511 = 0.047
1.486 = 0.042
1.558 + 0.095
1.523 + 0.057
clock-time
1,580 + 0.050
2.680 + 0.100
0.080
0.800 + 0.010
CaGP
MAE
0.782 0.042
0.908
=.861000099
0.826 + 0.087
NLL
1,483 = 0.047 1.477
0.009
1.542
1.479 + 0.073
clock-time
2.570 + 0.040 3.640 = 0.220
0.810
0.280
0.960 = 0.040
RCSVGP (const.)
MAE
0.763 # 0.065 0.921 # 0.082 0.799
0,8567 0.095 0.10
NLL
1.482 # 0.066 1.484 # 0.059
1.496 = 0.089
1.5494 = 0.0753
clock-time
1.620 + 0.040 2.700 = 0.060
1.050 = 0.030
0.830 = 0.020
RCaGP (const.)
MAE
+774340005
0.762 0.097
0.811 = 0.077
NLL
1.475 # 0.062
1.490 = 0.036
1.533 # 0.099
1.469 + 0.075
clock-time
2.630 = 0.040 3,770 0.110 0.960 = 0.050
1.080 + 0.020
Table S2: UCI Regression datasets results in the presence of focused outliers with e 0.1.
Average test set mean absolute error, negative log-likelihodd and clock-time (in seconds), with std, for 20 train-test splits. Bolded results refer to the best baseline. Lower is better.
Boston
Energy
Yacht
Parkinsons
SVGP
MAE
0.627 0.051
0.337 0.021
0.630 0.053
0.672 # 0.078
NLL
1,405 # 0.031
1.239 + 0.009
1.255 # 0.057
1.433 + 0.053
clock-time
1.500 + 0.050
2.800
0.080
1.000 # 0.060
0.770 # 0.020
CaGP
MAE
0.677 0.061
0,436
0.525 # 0.050 0.750 # 0.073
NLL
1.425 + 0.055
1.181 + 0.03
1.277 = 0.048 1.462 = 0.064
clock-time
2,460 = 0.060
3.670
0.830 = 0.290 0.930 = 0.010
RCSVGP (const.)
MAE
0,709 0.100 0.070
+ 0.551 = 0.079
NLL
1.089
1.496
1.321
1.392 0.087
clock-time
0.040
2.850
1.040
0.040
0.030
RCaGP (const.)
MAE
0.505 + 0.069
0.090
0.437
0.053
0
NLL
1.323 2 0.073
L.163
1,268
0.051
1 0.082
clock-time 2.570 0.050 3.780 0.110
0.930 0.060 1.010 + 0.040
Table S3: UCI Regression datasets results with RBF kernel. Average test set mean absolute error, negative log- ligeliiel,oodd, and clock-time (in seconds), with 1 std, for 20 train-test splits. Bolded results refer to the best baseline. Lower is better.
Asymmetric outliers
No outliers
Boston
Energy
Yacht
Purkinsons
Boston
Energy
Yacht
Parkinsons
SVGP
MAB
1.145 0.052
1.138 0.040
1,235 a 0.044 1,002 a 0,100 0,731 + 0.052
0,908 0.026 0.767 0.063 0.811 = 0.078
1.691
D.043
1.736= 0.037
1.732+
1,590
1.433
0,053
1,420
1,448
clock-time
0.930
0,280
1.170 = 0.070
0.650 + 0.050
0.720
+ 0.260
0.940
+ 0.030
1,230
+ 0.150
0,870
+ 0.240
0.670 = 0.130
CSGP
MAB
1,060 0.005
0971 0.090
1.133
0.955
0.110
0.48%
+ 0.442
0.334
0.024
0,512
+ D.02
0.674 = 0.077
NEL
1.611
0.086
1.535 0.071
1.679 0.088
1,556
0.088
1.278
1,106
0.011
1,263
0.042
1,362 0.067
clockttime
2.470 0.060
3,540 0,080
0.640 0.100
0.870
a 0.080
2410000000
3.540
160
0,770
0.060
0.830 0,100
RCSVOP
MAB
0.080 2 0.122
0.456
0.36600044 0,041
0.735
0.840
2: 0.157
0.31240.074
0.555
2 0.035
0.640
0.074
0.540 0.08.
NLL
1,364 0.080
1,200 10000000 0.026
1,381 0.084
1,518
0.116
1,302 1 0.063
1.30010.003
1.318
0.017
1.417 1 0.065
clock-time
0.090
2.560
0.970
1,530
0.040
0.720
RCGGP
MAR
0.753
0.353 0J053
0.651 0.095
0.913
01
0.379 a: 0.050
0.326
0.023
01366
0.005
0.534 = 01.7%
NLL
1.405 + 0.118
1.173 + 0.024
1,336 + 0,061
1,531 + 0.152
1.241 + 0.046
1.144 1:44++0.000
1.226 0.004
1,307 = 0.066
clock-time
0.060
0.090
0,900
0.030
0,940
=
Stndent-
MAE
0.562 0.037
1.015 0.028
0.907 0,061
1,225 1 1.0099
0.396 1 0.001
0.970
0.025
0.815
1 0.07]
0.365 2 0.077
NLL
1.619
0.039
1.606
0.024
1.564 1 0.067
16658
1 0.085
1.595
1 0.063
1.572
0.026
1,559
1.100
1.627 1 0.092
clock-time
1,030 a 0.040
1,350 =000=000 0.070
0.750 0.060
0,730++0.04)
1.000 +: 0.030
1,460
0.000
01920
0.050
0,730 : 0.010
B Computational Complexity
Assuming a constant mean prior, RCaGP exhibits the same time and memory complexity as CaGP, as both rely on the same batch update algorithm. Note that the weight function requires (n)) for time
15

---

and space complexity, thus it does not affect the complexity of the algorithm. The time complexity of RCaGP is O(ni i max(i, k)), where k n/é is the number of non-zero entries for each column S,. In addition, RCaGP and CaGP have linear memory complexity O(ni).
C Evidence Lower Bound for RCaGP
Proposition C.I. Define the variational family
Q _ (q:(f) = (X) E [=*1]
(S1)
and robust loss function
L (f. X)= = 2n
Tromis Jw' 2F V C(x,y, y, y. mnosse)
(S2)
where noise J (y
mm) and C(x,y, y. = ooois diag( ( - the evidence lower bound of RCaGP is given by
w?. Then,
RCaGP ELAO
(F) (ylf)] (f)
(S3)
1
tr noise Jw
1/2 X, X 122
2 J J (x)) 1/2 (x)) (x)
v
2
1
J( (noise)
[.
S KS. + log det(S KS.) log(o (nnisee)
2
det S,) - llog -r((s) KS.) S KS,))
(S4)
where v; (S - (y mw) and (y
X)) are the projected
representer weights and the
respectively.
Proof:
The ELBO given by
RCaGP ELBO
(r,, [logp" (L(((())p()))
(S5)
We first compute the expected loss function term:
/ log
g (F) p" (y/f)]
(F)
(S6)
/
1
(f TJ J f 2f v + C(x, Raisa Da (f)dd
(S7)
2
/
1 f' Tomm
mniso J I (f) df + (Q
2
C(x, y, o noise) (f)
(S8)
2
J
r noiss -
f. y: (f) df (n [PT vl C(x, y, o
(S9)
Next, we apply the identity
TT a noise
J= If ur(o (r(o 1/2 T
noisa J=//2 (see (2]):
noise
1 (n) trro J p/ σ Tuolee J 1/2 En if v) C(x, y, 2 moise)
(oogr)oopp (y(f)]
2
(S10)
Since the expectation of the trace is equal to the trace of the expectation, we find that
1
(r [log 1 " (ylf)]
- Srr(F (E (P) (o moooise J 11 "uoisee Jw 1/21) + (e) [f v] C(x. y, o2 moise)
2
(S11)
tr tr(EE (F) mnoisee J Jw12 1/2
Uoiise 1 + Eg [f] C(x, y, noise)
2
(S12)
16

---

By definition of mean and variance, we have that Ep(x) Xxx 1 Applying this identity, we obtain
V Vp(x) [x] + Bp(x) [x] Bx(x)
2 tr moiso Jw 1/2 ( [ff 1o noise Jw//2) (F) [f] 1 Cxx, y, moisse)
(log))loopp ((yl))
(S13)
1 1
2 noise Jw 1/2 (.)) + B8 (n )o
noise 1/3
+ -Bg(()) [T] v _ (x,, y, mnosse)
(S14)
2
tr(o noise Jw 1/2 + (X)) )o 11 JT + +X) Tv
2
1
C(x, noise))
(S15)
2
2 noise J 1/2 (X, (X.
tr(om
tr noise (X)); (
1/2)
+ (X) ((x, noise)
(S16)
2
- 1r(o moise 1//2 X, X)J
tr(o f T (i((X))
+ +x(XX) y moise)
(S17)
Next, we compute the KL term, where both (f) and p(f) are multivariate Gaussian:
1
det(K)
(f) ([ ||p((]] --(((()) - m(X)) (a;((X) m(X)) + log
(xe((k(Xxx)))
tr(K + - m)
(S18)
2 ((KK, mw (X))) (-)))Kk(K (KCily--y mm log deggde-K111 Ri
1
XX)
+ tr(Inxn - Č,K) -n)
(S19)
1 ((r mw (X)) T,JCC,K,Y--- C.KCC(y-m(XX) det -Č,K)
2
+ tr((Ixxn C,K) nn)
(S20)
1 [v.SS.
KS, (V..S.SSVV
log det (I Ixn Č,K) + tr (r[nn C,K) -n)
(S21)
2
1
v. S KS.V; - log det Č,K) - tr((S, tr((S,KK,) KS,)- (.S.))))
2
(S22)
Here, we apply the Weinstein-Aronszaji identity:
v S KS.V. - log det - (S. - KS,) ((s -ss.
(S23)
2
S...S. KS - let de, KS,) Ks, S. KS.)) - tr KS.) - KS.)
2
(S24)
1
(v.s.) S. KS. - log ((, 2 S.T J. (r Ks.) -1. KS.)
(S25)
2
V. S. KS. v, + (et KS.) log dot noise S J. tr((S) KS.) - KS.)
2
(S26)
1
2 (v. S. KS.V, v + log de. KS.. log(o log degdet(Jw)) - log de! (S.....s. S.)
(s kS,) - KS.))
(S27)
17

---

D Robustness Property
The following lemma contributes the proof of Proposition 4.1.
Lemma D.I. For an arbitrary matrice $ E 12mx and positive semidefinite matrice B R we have that
(SBST)= - S B- B-1
(S28)
where define G - 112 ( S-S and denotes the Moore-Penrose inverse.
Proof:
The whole proof is derived from answer to question posted on the Mathematics Stack Exchange Forums, which we write here for conciseness.
Denote Ô = -- 1-SS S S and H(a) = (S(al + +B ---
We also note that
($BS 1 lim H(a)
(S29)
By applying Woodbury matrix identity, rewrite H(o:) as follows:
-1
H(a)=
= SS .B-- 1/2 I+ -B1
1
(I+
B-1/22T
(S30)
a
α
a
a
Using the fact that SST is invertible and applying the Woodbury matrix identity for the second time, we obtain
((a) o(SS - 1/2
1/2
(- + B 1 --15 ($$ > BB "B
(S31)
= =a($ŠT)-1 + + ( 1/2 1/2 - 1/2)-
a
B
(S32)
We note that
S STSST)) S
(S33)
I STST(SS) - = Ô
(S34)
Then, rewrite H(a)) follows:
H(a) = (a + S 1/2
[+ T+- -B - 12200--11
B B-1/2
(S35)
Applying the Woodbury matrix identity for the third time gives
H(a) = (SST + S 1-1
B 1 ôB - 1-2--00-)-----11 1/2
(S36)
Since the Moore-Penreee inverse of matrice A is a limit:
lim A T A + al) AA lim a=0
(AAT + aI)
(S37)
a-
We can take the limit of (() as a - 0 and apply the limit relation above to obtain the following result:
( S+T 1/2 (I -1- B- - 1/26)++ -
112S++
(S38)
&
PIF RCaGP We now prove Proposition 4.1: Let f GP(m, k) denote the RCaGP prior, and let i E 0, .7 represent the number of actions in RCaGP. Assume the observation noise is
18

---

~ (0,, moise D). Define constants C2 E R; = 1,2.3, which are independent of you. For any i and assuming supx.y w(x, y) oc, the PIF of RCaGP given by:
PIF PIFRC (v D) C2(w(Xm v y C3-
(S39)
Therefore, if supx,y y w(x y)² suPas. RCaGP (y D)1 < xO.
00,
RCaGP regression
robust
since
Proof:
Without loss of generality, we aim to prove the bound for m - n.. We can extend proof for an arbitrary m E {l...., n.. Let (f = (f: a K,) and (f (f(;; p, KŞ) be the uncontaminated and contaminated computation-awar RCGP, respectively. Here,
A = m + KÖ
(S40)
K, k Joiso Jw
(S41)
i + KC,V
(S42)
K; KČO mnise we
(S43)
where w [w((1,31)) y1),................... ..... w(x,, v.ll The PIF has the following form
1
det(K;)
PIF DIFRCCGP(y D,i)
Tr( Tr((KK))- n (AC A.)) (KG)
(p
i.,) + In
det((K))
(1)
(2)
(3)
(S44)
We first derive the bound for (1):
(1) Tr((K;))
(S45)
Tr (KČčO Thisso we) (koooio
(S46)
Tr(o (i) Cio J (w) - 71
(S47)
Tr(o uoise
J we ( (rr( Cio Jw)
(S48)
Tr(o noise Ju (T,) C)Tr(C Pn
(S49)
The first and second inequality come from the fact that Tr(AF) Tr(A)Ir(F) for two positive semidefinite matrices A and F. Since Tr(C we can write C₁ TriCCoo Let (S K*S such that C S. BST. Observe that matrice B is positive semidefinite. Thus, we can apply Lemma D.1 to obtain the bound of Tr((C,)-)):
does contain the contamination term,
Tr((C)))))) -Tr(SS BS!)
(S50)
= Tr(S 1/2 1/2GB-1/23 1/2S;)
(S51)
< Tr(S S (B
(S52)
Furthermore, we derive the bound of Tr(G) follows:
Tr(G) = Tr(1 1/2 S. SB 1/2
(S53)
(S54)
1/2 < n 1/2
(S55)
The inequality S52 stems from the trace circular property and the properties of the product of two positive semidefinite matrices. Notably, observe that Tr(G) n since B 1/2 112 and (I (I) S S in S55 are positive semidefinite matrices. By definition, the trace of
19

---

positive semidefinite matrix is non-negative. We then apply the trace circular property for the second time obtain
) s s
(S56)
1 =nTr(S, ) JT(())
(S57)
= =C2Tr(K + o JooseeJ
(S58)
where we define C2 = TTr(S (TTT(S We then plug S58 into S49 to obtain (1) S Tr(o noise w- Tr(K on noise Jw (210)C0221 -n
(S59)
(40 noise w0 (x
Kik
nolsel 1 w 2 Wh)) (X5t)
C,O2 -n
(S60)
k=1
ssu moisee 1 w 2= (xx + Σ Kit )) ci0s-n==0 k=1
x
n sup 40 4 noise w²
(x,
n
x.y
x.
(S61)
Next, we derive the bound for the term (2). Following [2], we have that
(2) Amax ((K;)) fa,11::
(S62)
where Amax( (K;) 1) is the maximum eigemvalue of (KÇ)-1 We expand Amax (Kç)
1) to derive
the following bound:
Amax (K;9 Xmax (o JJ (C:)
(S63)
S Amax (5 "oiise J Jwd Amax ( -1 (Ca))"))mmx(KF Xmax
(S64)
= Amax (o noise J w Amin Ci) Amax (K 1K-
(S65)
Amax (o noise Jwh Aminn (K²) (Ko)- Danax (K-1)
(S66)
Amax (o nooiso w! ) (Amax (K) + Amax (o wc )) (ammaxxX)) (K-1)
(S67)
The first inequality follows from the property of the maximum eigenvalue of the product of two positive semidefinite matrices. Subsequently, the second equality is derived from the fact that a matrice's maximum eigenvalue equals its inverse's minimum eigenvalue, By the definition C (Ka)-1-3 11-1--17 S, and the fact that Öj (K-)-1 and E, are positive semidefinite matrices, the minimum eigenvalue of Ke is less or equal than of C; resulting in second inequality. The last inequality follows from the equivalence of the maximum eigenvalue and the addition property of the maximum eigenvalue of two positive semidefinite matrices.
Since w = diag( mntee (w)) and supx.y w(x, v) < oo, it holds that Amax (o (o Jw!) = Co < 2
to0 and Amax ( ( 2 J w = O3 < +00, such that
Amax (KK--)) S (K) C5) (oAAmaxX--0 (K-1) Co
(S68)
We substitute Co into (2) to obtain
(2) < CollAf - 11
(S69)
Co [Im KVØ) (m KV:)l
(S70)
Co [K(Č)(y" mw=) C,(y )l
(S71)
$ C. K|| || (y - mw:) - C,(y - mw
(S72)
S q Co |K| P ((K9)
mw=) -(K) (y mw 11?
(S73)
Co |K| [K + o Jw -1(y mwo) (K + 2 Jw) - mw)l
(S74)
for a constant y < +oc. The second equality comes from the definition of v. and the fourth row follows the Cauchy-Schwar inequality. Finally, the last inequality holds since K; (C, Σ.)
20

---

Applying results from [2], we obtain
(2) - q Co | [K + o Jwo)
mw:) - (K + ww) mw)lf
(S75)
S c K 2((C; + Cx) + (Co + C10) (oo))(w((())))) un) y()
(S76)
- C11 + C12( ( vii?
(S77)
where Cn = qC' K || Cs) and C12 [Coo|K]] The terms C7, Cs. Co. Ci0 correspond to Co C₁2 T7 Cg in [2].
The term (3) can be written as follows:
33 In
(S78)
dde((k))
dde( MeetCCjo maisc wo)
In
(S79)
det((C,o Ioosso
In(det(o noiso I C det Ö (eet)eeeeee (oosse
(S80)
Observe that we can write C13 = In(det(onois noise J C since it does not contain the contimation term. Furthermore, we obtain
(3) In(Cia (C)) del Jws))
(S81)
IeeC (k)) det (o moisoot we)
(S82)
det inoiooo we)
In C13 det(tK KK + mooie wc)
(S83)
det(o mooss ws)
Iin Ci3 det(K) + det(o noise 2 w
(S84)
The first inequality follows from the determinant property of positive semidefinite matrices. The last inequality leverages the fact that det(A F) > det(A) + det for A and F are positive semidefinite matrices. Since (e)( do (noiis w=) 2 0, we find that
(o nn
In det(K) + dot noise Jw.)
S
(S85)
leading to the following inequality:
(3) S (n(s))) C14-
(S86)
Finally, putting the three terms together, we obtain the following bound:
PIF (vfm D,1) I C3 + Cn + -C22(u(an. vn Ci4
(S87)
= C/ (w(En, vn v,o? + +C2
(S88)
where =²= and C2 = O3 + C11+
E RCaGP Capture the Worst-case Errors
Proposition 4.2. Let kel..-) (.) + o.(:.-) be the combined uncertainty of RCaGP with zero-mean prior m. Then, for any new x E X we have that
sup (h(x) (x))) = k,,xx, x) + to
(20)
[Allage S1
sup (p((x) - (x)( (x,
(21)
[|lllgg S1
where p=())=
is the RCGP's posterior and (2) (-) = kl, ()C,(((( m (X))
RCaGP's posterior mean for a function Hgr = g(x) - mx with latent function g TR. and the shrinkage function mao X x + R
21

---

Proof:
Let cj (X)]; for j = I,..., n, where we define k. (..-) k(.,.)
Sugius ónl(
w (x, v)
x' and D
8w 8 (x,,x)==
2
x' and D
(S89)
0
x +
where D denotes the observation set. Then, applying [21, Lemma 3.9] provides
2
sup
h(x) - --- 11{(x))
sup M Man
h(x)
Σ
c
(S90)
|| $1
j=1
= [1,"(,))) k(x, x) C) 113m
(S91)
(,,),, k" (1 x)){0 22% 2(k (,,x), k(x, XX,,X,X,,,,CCC 1)
+ (k(x, ()(xxX)CC C, (X. k(x, (XC,
(S92)
By reproducing property, we have
sup h(x) -1((x) = (x, 2x, X)C: + xx"xxxx">xx>>>"">>>>>""">>>">>>> (x,,x, X)C.A"(x,,", (X,x) [Mllige 51
(S93)
since * Xy, holds that kn (x, X) k(x, X). By definition, have k" (""(XXX K K and following in [46, Eq. (S42)], it holds that Č,KC, C,. Therefore, obtain
Oigine J w.
sup [Mll
M(x) 1P(x))
k(x,x) + onoic 2k(x, xXx + KJCK,CCKKK
(S94)
= =(x;x) + Guoise
k(x, X)C,((X
(S95)
xx(xx)) x)
-noise
(S96)
For the last result, we analogously choose Cj = [[K-1 -C,)," (X,,),
Then, we obtain
sup (x(x)) - 11(())
sup
Σgh(x,))
(897)
[Alllnm S1
IAllik $1 J-U
= -|A(x, X)(K C.) I1
(S98)
= (x, XjK KK "X, 2x" (x, XIKK
(S99)
1xxx,, X)(K C.)x(Xxx)
(S100)
(x,,
(S101)
The consequence of Proposition 4.2 is then formalized through the following corollary:
Corollary E.1. Asxume the conditions of Proposition 4.2 hold. Let μ. the corresponding mathematical RCGP posterior mean and fas the RCaGP posterior mean. Then, it holds that
[h(x)= p=((x)]
S (x,,xx) o2 moise
(S102)
|1 |
j(xx) (x()
S (xxxxx
(S103)
|
The proof follows immediately from Proposition 4.2 and notice that h h/]l/]]| Hew has unit norm. Corollary E.1 implies that the computational uncertainty provides the pointwise bound relative to RCGP's mathematical posterior, and the combined uncertainty is the pointwise bound relative the shifted latent function.
22

---

Additional Theoretical Results
F11
CaGP Lacks Robustness
PIF for the CaGP. We aim to prove the following statement: CaGP regression has the PIF for some constant C3 E R
PIF Ca D,i) = Ca((um Vom) and is not robust: PIFcuup (y D,i) + 00 Wial + 00.
(S104)
Proef:
Let p([[D) Mu, K.) and p([[D;.. computation-awes respectively. Here,
NL K.) be the uncontaminated and contaminated
11 Kv:
(S105)
K. iC
(S106)
1. m Kv;
(S107)
K; KC
(S108)
Note that both K, and Ki share the same matrice C. Then, the PIF has the following form:
PIF cacP (yms D,1) (TT(K; - (4) - -4,) (K)) (4f -1,) + In
det((K))
det(K.)
(S109)
Based [2]. the PIF leads the following form:
1 PIFCaGe (ym: D.1) = (14 m1) (K:)- (KK)) : m.))
(S110)
Notice that the term μa - 14; can be written
n. μ4 (m Kv;) (m Kv.)
(S111)
-K(v; -i)
(S112)
KyC.(y" - -m) C.(y m))
(S113)
K(C,(y" -y)
(S114)
Substituting the RHS of Eq. (S114) to i μi in Eq. (S110), we obtain
PIFcaGr(Uj D,1) = CC, - y)) " (KC K(C.ly- -))
1
(S115)
2
nosse (y y) K(y y) -
(S116)
Note that y and y" have only one exception for the m-th element. Thus, we have
PIFcaGr(yi. D, = C. Ko "I mmmmmUi ym
1
(S117)
F.2
Mean convergence of RCaGP
F.2.1
Empirical-risk minimization problem of RCGP.
We first show the corresponding empirical-risk minimization problem of RCGP whenever = 0. Following [2, proof of Proposition 3.1], rewrite L" and formulate the RCGP objective follows:
1 Ex-lj= 28 JJ (y mw)
f argmingeH, 2n
|
(S118)
L
for À > O. The regularization constant À controls the smoothness of the estimator to avoid overfitting, where the larger the is, the smoother the resulting estimator f' becomes. Next, we show the unique solution to $118 through the following lemma:
23

---

Lemma F.1. X > and the kernel & is invertible. the solution to S118 is unique, and is given by
(xx) = kx (K + AJw) ' (y mw) Σ
(S119)
==1
where
(di..., (K ) mw) RD
(S120)
Proof:
The optimization problem in S118 allows us to apply the representer theorem [37]. It implies that the solution of $118 can be written a weighted sum, i.c.,
f Σ x,)(,,))
(S121)
for Q1, , On IR. Let â (āi...... 6n) T E IR.. Substituting S121 into S118 provides
argmin
(x-'a à KJ kā 2 KJ (y mw)) + & Ká
(S122)
2n
Taking the differentiation of the objective w.r.t. setting it equal to zero, and arranging the result yields the following equation:
K(K + +nJJ)á K(y mw)
(S (S12)
Since the objective in $122 convex function a, we find that (K + AJ mw) provides minimum of the objective ($118 and $122). Furthermore, we can verify that T. is a convex function w.r.t. f. Therefore, conclude that â (K + AJw)) ) (y mw) provides the unique solution to $118. As remark, Proposition F.1 closely connects with [21, Theorem 3.4].
F.2.2 Proof of Proposition F.3
Relative bound errors. We first provide the equivalence [46, Proposition 2]:
Proposition F22. For any choice of actions relative bound error p(1) || volle ((1) Vll is given by
p(i) = (1 & Amax C, $
(S124)
where l/llVl
The proof is by directly substituting C., K, v. in [46] with Č,, K, v, respectively. Next, establish the convergence theorem of RCaGP through the following proposition:
Proposition F.3. Let Hss be the RKHS w.r.t kernel k, Joois > 0 and let js, Hr be the unique solution following empirical risk minimization problem
argmin JeHk 2n
P
J
28 J mw)
+
(S125)
which is equivalent to the mathematical RCGP mean posterior shifted by prior mean m. Then for the number of
actions i [0.....n] the RCaGP posterior mean P. satisfies:
liا - Pull, S p(1) c(Jw) l f. -m0m-m/le m||
(S126)
where is the relative bound corresponding the number actions and the constant c(Jw)
Auaa (Jw) + as Amaze (Jw)= + O.
Amia [K)
Proof:
Lemma F.1 implies there exists unique solution to the corresponding RCGP risk minimization problem.
Choosing p(i) described in Proposition F.2, we have that // Valle
p(i)]l vo where vo 0.

---

Then, for E (O.....2) we find that
v.l 1 (1) voll || - v.lli - ||
(S127)
[1 vollic + Arnaxe (Jw)
0000²
A Amin (K)]l8 110 vollē
(S128)
Amin (K)
[1 vollis + Amnu (onx((ww)
S (())
1/6 No liss
(S129)
Amnio (K)
Amnre (Jw)
S p(i)²
11
|| vollic
(S130)
Amin (K)
The third inequality stems from definition of Jw and matrix is the largest component of its diagonal. Assuming (-)() = m(.) + -;= (v.),A(,,)) = m(.) +
maximum eigenvalue a diagonal
le., Cy mmw and applying result from [46], we have that
||v villic = I p
(S131)
Combining both results and defining c(Jw) = -(1
2max (Jwc
Amin (K) ). we obtain
Aixe 24 II ||| S ((t)c(Jww)]l wlll vol.k p())c(Jw)ll. ml/
(S132)
F.3 SVGP with Relevance Pursuit
Proposition E4. Let Di, = (x3335) : $ 1}, Z 2m is the inducing points, p = Pie + pies where P.Pg R. [p,]e and e; is the ith canonical basis vector. Then keeping Pry Piy fixed,
By (Di)) - V. (el [ VK,, Qrla (u.
pi = argmax,, p. + =
(S133)
where yı [xx:) + es. The quantities can be expressed as a function of A = (Qr + D
Puiie +p)
E W: [A- y)-\A\./AA-'l and Valyi Drd 1/[A 'Ja
(S134)
where D 2 is diagonal matrix whose entries are oooie + p.
Proof: Following [41], the evidence lower bound (ELBO) of GP can be written as
-2ELBO(0) = (2m) log (A) + y A-'y tr(D
KK Qr)).
(S135)
p-o hoise
where Qr pta2 moise and Qr Kxa K K
Following [3]. partition matrix A to separate the effect of and use Schur's complement:
A + ú3,ú
(S136)
B.
where
A; [Qr) + D Pretto
(S137)
û = A. K.a K2 NoT (xi, Z),
(S138)
8. (Ik(xi, x:) + oniio + pi) -k(xi, Z)K, "IK,KKJ.A KrrK-th TT (x:, 2))
(S139)
Quadratic term:
y (A p+o
ATA ) --y==
B, Tu -
(S140)
noter
Determinant term:
det(A) (et((rr D
P. +o 2uuin
= det(A
(S141)
Following [3], formulate ELBO difference as 2(ELBO(p) ELBOO(P.)) )) - -3,Yy y1)² - 1og(3,
p. (.(Ka. (Qraa).
(8142)
The derivative of the difference ELBO w.r.t. Pt. is
0p (FFLBO()))))) ELBB((p))) (yT,ū v.)2 B.((K, -Qel:))
(S143)
25

---

Using the quadratic formula, obtain the solution $1,2 1 1 (- (Qrlie] Note that (y
2(x "-
v1)²² ; 0 and K - [Qr 2 0. Therefore, the smaller root P2 is negative when (K,, - (Qela) > 0 and zero when (K.: - [Qr]s]
Based on the / formulation (predictive variance cannot zero), we reject &2 Similarly, reject A2 0, following [3]. By ignoring the constant term in 82 and following the approach of [3], we solve 3- 1/ (K (v,) yi) for and projecting to non-negative half-line, get
(r v.) (k(x,, + onoise k(x;, 22 AxaA, KxnK, k" (x: Z))
(K.. (Qr]a)
(S144)
Lastly, we note that
N.u y (A KxrK, K,,K (x:, 2)) y = Ey(ly./Drd] yi:
(8145)
((x,, x;) + ouooso k(x,, ' A KxrKa'k (xx, 2) -.dl.
(S146)
Following [36], we have that
(D -A |AT")\/(AAT))
(S147)
V br [ 1/[A 'Ja.
(S148)
which follow the LOO predictive values.
G Details of Expert-guided Robust Mean Prior
In this section, we provide the posterior inference details of our probabilistic user model for defining the informative mean prior. We first discuss the posterior inference of latent variables 8,, governing human expert decisions in identifying outliers. For every outlier candidate (x., (o), we have that
p(3.)= B(aa an, B = =.)
prior,
(S149)
p(ōo do) = 6.0 + (1 -1- (1-b.)(1 - 5o)
likelihood,
(S150)
p(..)0,) x p(5.) p(ō.)
posterior
Bio
B (o) (5.) + (1-5.)(1 5.)) Bia ao oo, 3 = 30 + 0o),
(S151)
where we define the hyperparameters og and 3o as follows:
ak. 2m
and
8. 0.
The posterior (o.. exhibits conjugacy, allowing us to derive a closed form expression, i.c., Beta distribution. In this framework, we define a, using the z-score of the outlier candidates. The z-score information effectively recognizes outliers by comparing each data point the sample 9% close zero, reflecting that the expert confident that do is an outlier. Given the posterior p(6,) ( the expected value (3.1)_))) [8%] is expressed by:
and standard deviation. Next, set
a ao + 0o a. Bo +
Epic a. 16.]
(S152)
+
Next, we outline the posterior inference details for the outlier corrections provided by the human expert:
p(u)) = NNH no. 0², =;=)² p(oolTe) = N(u= = Too 0² o2or.)
prior,
(S153)
likelihood,
(S154)
p(RolTo) p(Fo) p(o)
posterior
(S155)
= -N(u Hoo
= = N = noo 0², = =
(S156)
=N To Ho + ooon Bo o² (To + 1.orenn.) )-1)
(S157)
+ 1.
where define the hyperparameters Ho and To follows:
P.
yj and
(yj Pa)
(NNN)(M))
26

---

Here. (x.) denotes set consisting of J closest neighbors to Specifically, we compute kernel function k(x,, D \ D) E IP (x--), where D represents inliers. Then, we retrieve (Xi....,x.) xx) corresponding with the largest k(x., D) E -1'(==0), o) set
our experiments, ensuring each outlier has four neighbors. This posterior also maintains conjugacy. yielding normally distributed posterior. For each identified oudier we define the prior for the latent variable oi using the sample mean and variance of its neighbors. Intuitively, this prior reflects the assumption that human expert's correction for a particularly identified outlier will be close its neighbors. This prior resembles prohahilistic model inspired by case-based reasoning. We then obtain the expected value Ep(Ao18)) as
To Mo + Feorr. Yo
Bp( (Bo)
(S158)
To L.FFoorr.
2
which equals the posterior mean. The algorithm implement the informative mean prior is presented in Algorithm
Algorithm 1 Expert-guided robust mean prior algorithm
1: Input: D. D. o, y, ozo..
2: j + D
Σ 1= yj
[P]
3: ō + V DD-
1
Σ 1ーー1 [1/5 xi)2
4: for 1,..., do
5:
a. + (9.- - p)/ol
6:
Bo + 0
7: No t Ho = = - (x),)) V/j
8:
To + J EyEyEN (V) mo02
to 9: Ep(5.1a.) 15o] f- a.t3.TI
Ep(RolB.)((o to oagee.
tattenir.
11: end for
12: m(x) E. Ep(5, 1o)l (p(
Output: mean prior mxx)
H Benchmark details
H.1 UCI regression datasets
Boston The dataset consists n 506 observations, each representing a suburban or town area in Boston. It encompasses d 13 features containing data like the average number of rooms in dwellings, pupil-teacher ratios, and per capita crime rates. We try to predict the median price of homes residents own (excluding rented properties). The dataset can be found at https: :///////wwww cs toronto... .eeaant..euuuaa// wwedeve/daaaa/beeedeeeeede bostonDetail.. html.
Energy The dataset describes the energy efficiency of buildings by correlating their heating and cooling load requirements with various building parameters. It consists of 768 data samples, each characterised by d = distinct features, with the ultimate goal of predicting a single continuous response variable found in the last column. The dataset can be found at attpss://arde eff efficiency.
ics uci www.uuu/ii.a............ dataset//222 /2e2////r////// enerryy+
The dataset's 's main focus is on predicting the residuary resistance of sailing yachts during their initial design phase, a critical aspect in evaluating vessel's performance and estimating the essential propulsive power required. This prediction relies d 6 primary input parameters, which include the fundamental hull dimensions and boat velocity. The dataset contains 308 observations. The dataset can be found ics.uci. edu/ http///aaastt///////caa/////arra
Yacht
https://aachie
The Parkinsons Telemonitoring dataset contains 197 data samples and 22 input features derived from voice measures. The dataset comprises collection biomedical voice measurements from 31 individuals, including 23 diagnosed with Parkinson's 's disease (PD). Each row represents one of 195 voice recordings, while each column corresponds to specific vocal feature.. The dataset can be found at https://archive. ics, uci edu/ htt//daaaset//7// wttpaaatttet////aaatooooo
Parkinsons
27

---

H.2
High-throughput Bayesian Optimization problems
Hartmann 6D.
The widely used Hartmann benchmark function [39].
Lunar Lander.
The goal of this task is to find an optimal 12-dimensional control policy that allows an autonomous lunar lander to consistently land without crashing. The final objective value we optimize is the reward obtained by the policy averaged over set of 50 random landing terrains. For this task, the controller setup used by [9].
The rover trajectory optimization task introduced by [45] consists of finding 60 60-dimensional policy that allows a to move along some trajectory while avoiding set of obstacles. We use the same obstacle set up in [30].
Rover.
Lasso DNA. We optimize the 180 dime-iinaii DNA task from the LassoBench library [49] of benchmarks based on weighted LASSO regression.
H.3 Additional outlier contamination protocols on UCI regression datasets
Our outlier contamination protocols follow the same settings as [2]. We recall the asymmetric outlier case described in the main text in Section
Asymmetric We sample uniformly at random 10% of the training dataset input-output pairs (x, 31), and replace the yi's asymmetric outliers, i.e., via subtraction of noise sampled from uniform distribution 44(3,997) 97), with ā being the standard deviation the original observations.
Next, Section 7.4 considers two additional kind outliers:
Uniform Instead of always subtracting noise as was the case for asymmetric outliers, the set of training data points selected for outlier contamination is divided in 2 subsets: half of the selected subset is contaminated by adding 2 U(3 9c)), while the other half contaminated by subtracting U33 9c).
Focused In this outlier generation process, we ramdomly select and remove subset of data points, which will be replaced by outliers. For these outliers, we deterministically choose their values in X. To do we calculate the median value for each input data dimension j. However, dircetly. Instead, we replace the removed input values by (m, + 81, m2 82 ...,, ma + 5a) where m, is the median in the j-th input data dimension, and 6, = u, where 0y is the median absolute deviation the j-th data dimension times 0. and u ~ (0, 1). Simultancously, the outlerr values on y are obtained by subtracting three times the standard deviation of the median of the observations M,. To not have the same value for every outlier position, also add small perturbation 6, = U, where ay the median absolute deviation of times 0, ~ (0, 1).
do place the outliers this median position
I Experimental Details
L1 Hardware
For the UCI regression experiments, all models didells-including our proposed method and the baselines were executed on a compute cluster consisting of two machines: one equipped with dual AMD EPYC 7713 processors (64 cores each, 2.0 GHz), and another with dual Intel Xeon Gold 6148 processors (20 cores each, 2.4 GHz). The longest-running UCI task, based the Energy dataset, required approximately 15 minutes complete, including both our model and all baseline methods. Bayesian Optimization experiments were conducted on a cluster with four NVIDIA V100 GPUs, each with 32 GB of memory. Among these tasks, the Hartmann 6D benchmark was the fastest, completing in roughly 20 minutes, while the most computationally intensive task DNA reeuure up 20 hours of runtime.
28

---

1.2 Hyperparameters
UCI regression
Optimizer
ADAM
Learning rate
0.01
Minibatch size
n-data
Number of iterations for optimizing ELBO
50
Proportion of test set
0.2
0.2
p-outliers
0.1
n-inducing points (SVGP)
100
Projection-dim (RCaGP)
5
Mean- Mean-prior m (RCaGP and RCSVGP)
Σj 9/
Mean-prior m (SVGP)
0
B
1.0
o 2
1.0
High- high-thrroghhut BO
Tig (data initialization)
250
ADAM step size for query x
0.001
ADAM step size for RCaGP parameters
0.01
The number of expert corrections
Ticrations /20
20
Minibatch size
250
-
1.0
p-outliers
0.25
Mean-prior m (all models)
0.0
Projection-dim (RCaGP)
25
n-inducing points (SVGP)
25
Table S4: Hyperparameter settings used for RCaGP. Most of them remain consistent across all the tasks.
(D 1S Hermmnn 60
80 (D00 (9) Limer 12b
BO (DDe (S) Bover 6/1
RCSVGP
CaGP
RCaGP (ours)
Figure $1: High-throughput Bayesian Optimization task under asymmetric outliers using DPP. BO and the Thompson Sampling (TS) acquisition strategy. Each panel shows the best value found each iteration found so far, averaged across 20 repetitions + 1 std.
29

---

CaGP
SVGP
RCaGP (ours)
TarB0 (DPP Luaer 120
TusBo (DF? EI) Bower 600
TuR8o (DPP DNA 1800
0311
RCSVGP
SVGP
RCaGP (ours)
Figure S2: High-throughput Bayesian Optimization task using TurBO under asymmetric outliers. Each panel shows the best value found each iteration found far, averaged across 20 repetitions + std. Top row: EULBO-TURBO-EL Bottom row: DPP-TURBO-EL
RGGGP TEa0 NLL
RCVVC?
RESVIG NLL
1bcht
Parkinson
10
NLL
MAE
MAE
4.5
Figure S3: Dataset-specific results for the ablation study conducted in Section 7.4, varying in weight function (Equation 12).
30

---

NeurIPS Paper Checklist
1. Claims
Question: the main claims made the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
The claims made in the abstract and introduction are theoretical, methodological, and empirical, and these are all backed in dedicated sections of the paper.
Justification:
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 results can he expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
A limitation paragraph has been added to the conclusion.
Guidelines:
The answer NA mcans that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors encouraged to create separate "Limitations" section in their paper.
The paper should point out strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well- speciflcciiii.nnnn asymptotic approximations only holding locally). The authors should reflect how these assumptions might he violated practice what the implications would
The authors should reflect on scope the claims made, e.g., 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, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to text system might not be used reliably to provide closed captions
online lectures because fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
applicable, the authors should discuss possible limitations 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 It acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play important role in developing norms that preserve the integrity the community. Reviewers will be specifically instructed not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: Yes, each proposition clearly states its theoretical assumptions. Proof's are deferred to the appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and refeeeeffeee
31

---

All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof' can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Companion code along with running instructions will be provided in the supplementary materials.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible important, regardless whether the code and data are provided or not.
the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can accomplished various ways. example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of a large language model), releasing of model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a way to access this model for reproducing the results a way to reproduce model (e.g.. with an open-source dataset or instructions for how to construct the dataset).
model (e.g.. large language model), then there should either be
(d) We recognize that reproducibility may be tricky in some cases, in which case authors welcome to deseribe the particular way they provide for reprodueibility. In the case of closed-soure source models, it may be that access to the model is limited some way (c.g., to registered users), but it should be possible for other rescarchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: All benchmarks datasets employed for regression tasks are known, publicly available online benchmarks. URLs are provided in the Appendix. The test functions employed for Bayesian Optimization have been introduced in earlier works introduced in the Appendix, and will be available in the code submit in Supplementary materials.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
/ http:/iissccccpubul/ cc/puiicc/ Please see the NeurIPS code and data submission guidelines (https: : http///iiss.aaa
www.bbmmsssoooooiioiy.o for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
32

---

The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (kttps://nipg) http:/iiiaa.cpppplii..c. guides/Code helpbbisssoooooiiioio for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for new proposed method and baselines. If only a subset experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release amonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g.. data splits, hyperparameters. how they were chosen, type of optimizer, etc.) necessary understand results?
Answer: [Yes]
Justification: We provide experimental details in Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error hars suitably and correctly defined or other appropriate informa- tion about the statistical significance of the experiments?
Answer: [Yes]
Justification: We report the mean and standard deviation for every numerical experiment computed across large number train/test split in regression tasks, or across repetitions with different seeds in Bayesian Optimization tasks.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
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, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to library function, bootstrap, etc.)
The assumptions made should given (c.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
is to report I-sigma error bars, one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CL. if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures
symmetric eтror bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments Compute Resources
Question: For each experiment. does the paper provide sufficient information the computer resources (type of compute workers, memory, time execution) needed reproduce the experiments?
Answer: [Yes]
Iustification:
We describe our computing resources in Appendix.
Guidelines:
33

---

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 well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (c.g., preliminary or failed experiments that didn' make into the paper).
9. Code Of Ethics
Question: Dnes the research conducted in the paper conform, in every respect, with the NeurIPS Code of Fthies https:/ httpe///werrsp...pppliiicc/////e wttppa/ce/aacaaauiieeeeleelleele
Answer: [Yes]
Justification:
We have read all the regulations and make to follow them.
Guidelines:
The answer means that authors have not reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain special circumstances that require deviation from Code of Ethics.
authors should make sure preserve anonymity (e-g., if there is a special consideration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This work foundational in nature and proposes general-purree ppreprsep Gaussian Process framework. It does not involve sensitive data, individual-leves prediction, or deployment decision- critical settings. As such, it does not pose direct societal risks or impacts.
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 unintended uses (e.g.. disinformation, generating fake profiles, surveillance), fairness considerations (e-g., deploy- ment 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 generate deepfakes for disinformation. On the other hand, it not needed point out that generic algorithm 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 intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative socictal impacts, the authors could also discuss possible mitigation strategies (c.g., gated release models, providing defenses in addition to attacks, mechanisms monitor- ing misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g.. pretrained language models, image generators, scraped datasets)?
Question:
Answer: [No]
Justification: Our work focuses on Gaussian Process methods and does not involve the release of models or datasets that pose a high risk of misuse. Therefore, safeguards for responsible release are not applicable.
34

---

Guidelines:
The answer NA means that the paper poses no such risks.
high risk 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 restrictions to access the model or implementing safety filters.
Released models that
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 encourage authors take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the ereators or original owners assets (e.g., code, data, models), used in the paper, properly credited and are license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
We cite the original creators of the codebase that we use for the experiments in the paper.
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, possible, include a URL.
The name of the license (c.g., ., CC-BY 4.0) should be included for each asset. scraped data from particular source (e.g.. website), the copyright terms of service of that source should be provided.
assets are released, the license, copyright information, and terms of use the package should he provided. For popular datasets, paperawi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors encouraged to reach out to the asset's creators.
13. New Assets
Question: Are new assets introduced the paper well documented and is the documentation provided alongside the assets?
Answer: [No]
Justification:
We do not introduce any new asset.
Guidelines:
The answer NA means that the paper does not release assets.
Researchers should communicate the details of the dataset/codemmodel as part of their sub- missions via structured templates. This includes details about training. license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset used.
At submission time, remember anonymize your (if applicable). You can either create an amonymized URL or include an anonymized zip file.
14. 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: [No]
Justification:
We not conduct rescarch with human subjects or crowdsourcing.
Guidelines:
The answer means that 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.
35

---

According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed 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?
Question:
Answer: [No]
Justification: Our research does not involve 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 expect authors adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
initial submissions, ble), such as institution conducting the review.
include any information that would break anonymity applica-
36


Reviews:
[{"review_id": "gVgppAdsep", "reviewer": "Reviewer_uEzz", "summary": "The authors extend the computationally aware Gaussian process (CaGP) model by integrating concepts from robust and conjugate Gaussian processes (RCGP) [Altamirano, Briol & Knoblauch, 2024]. They begin by establishing the model with an RCGP prior and then apply CaGP's variational inference approach to the representer weights of the conjugate posterior.\n\nIn addition to this core contribution, the authors:\n\n* Propose a prior mean function based on outlier detection from an expert (Section 5).\n* Demonstrate that the impact of adding contaminated data into the posterior computation, specifically the posterior influence function (PIF), is bounded in terms of the contamination (Proposition 4.1, Appendix D), an approach analogous to Altamirano, Briol & Knoblauch (2024).\n* Extend results from Wenger et al. (2024), based on Kanagawa et al. (2018), to directly link RCaGP's posterior variance with the supremum of the squared difference between the predictive centered mean and all $h(x)$ with an RKHS unit norm or less. This assumes the dataset is generated by a latent function $g:\\\\mathcal{X}\\\\to\\\\mathbb{R}$ and its centered version $h(x) = (g(x)-m\\_w(x))$ belongs to the RKHS of the data generating process (Proposition 4.2, Appendix E).\n* Empirically evaluate their model on four well-known small-scale ($n<1000$) UCI regression datasets and Bayesian optimization tasks, artificially introducing outliers.", "questions": "Could the authors clarify the definition of $m\\_w(x)$ in Proposition 4.2? While the meaning of $m\\_w(\\\\mathbf{X})$ is not explicitly clear, we can infer that it also depends on the vector $\\\\mathbf{y}=g(\\\\mathbf{X})$. However, for the assumption that $(g(x)-m\\_w(x))$ belongs to the RKHS of the data-generating kernel, are we meant to assume that $m\\_w(x)$ in this context signifies $m\\_w(x,g(x))=m(x)-\\\\sigma^2\\_{\\\\text{noise}}2\\\\frac{\\\\mathrm{d}}{\\\\mathrm{d}g(x)}\\\\log w(x,g(x))$?\n\nThe assumption that $(g(x)-m\\_w(x))$ belongs to an RKHS is similar to results from Wenger et al. (2024) based on Kanagawa et al. (2018). This assumption appears straightforward to accept in the standard Gaussian Process Regression (GPR) case, as the GP posterior mean is typically a member of the RKHS when the mean function is considered. In the current context, the shrinkage term alters this interpretation, as it effectively influences what would otherwise be considered the prior mean of the GP. I believe the authors should clarify this assumption and discuss the conditions under which it might or might not hold, thereby providing better context for Proposition 4.2.\n\n**Possible Typo:** In the inline equation on line 95 defining $\\\\mathbf{C}\\_i$, the inverse should be placed outside the parenthesis, rather than being attached to $\\\\hat{\\\\mathbf{K}}$.", "limitations": "Yes.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "I find this paper to be well-rounded in its contributions, translating expected results from its main source papers into this new setting. While some theoretical results might appear to be straightforward extensions of previous work, the proposal for the expert-guided prior mean and the comprehensive experimental evaluation provide interesting and valuable additions to the use of Gaussian Processes for robust regression, which should satisfy most readers.\n\nNonetheless, while the authors compare their method to the original CaGP and RCGP approaches, I believe readers might struggle to contextualize this proposal in terms of uncertainty quality and runtime when compared against non-conjugate GP regression methods, such as those employing MCMC or Sparse Variational Gaussian Processes (SVGP).", "was_revised": true, "final_justification": "I am keeping my score because I believe the authors have adequately responded to the concerns raised by myself and the other reviewers, especially with the addition of the new baselines. I recommend that the authors take our discussion into account when explaining their theoretical results. This will improve their clarity for a wider audience without compromising the precision of what is being conveyed", "rebuttal": "Thank you for your thoughtful review and for appreciating the overall contributions and experimental depth of our work on robust GP regression. Regarding the typo, thank you, we will update the manuscript accordingly.\n\n\nQ: Could the authors clarify the definition of $m_w(x)$ in Proposition 4.2? While the meaning of $m_w(\\mathbf{X})$ is not explicitly clear, we can infer that it also depends on the vector $\\mathbf{y} = g(\\mathbf{X})$. However, for the assumption that $(g(x) - m_w(x))$ belongs to the RKHS of the data-generating kernel, are we meant to assume that $m_w(x)$ in this context signifies  \n\n$m_w(x, g(x)) = m(x) - \\sigma^2_{noise} 2 \\frac{d}{dg(x)}\\log w(x, g(x))$?\n\nThe assumption that $(g(x) - m_w(x))$\n belongs to an RKHS is similar to results from Wenger et al. (2024) based on Kanagawa et al. (2018). This assumption appears straightforward to accept in the standard Gaussian Process Regression (GPR) case, as the GP posterior mean is typically a member of the RKHS when the mean function is considered. In the current context, the shrinkage term alters this interpretation, as it effectively influences what would otherwise be considered the prior mean of the GP. I believe the authors should clarify this assumption and discuss the conditions under which it might or might not hold, thereby providing better context for Proposition 4.2.\n\nA: We thank the reviewer for this insightful comment. As noted, the difference $g(\\mathbf{x}) - m_w(\\mathbf{x})$ belongs to the RKHS $\\mathcal{H}_k$ of the data-generating kernel if \n$$\\nabla_y \\log(w^2) \\in \\mathcal{H}_k$$\n\nThis holds when $\\nabla_y \\log(w^2)$ has a finite RKHS norm; that is, $\\| \\nabla_y \\log(w^2) \\|_{\\mathcal{H}_k}^2 < \\infty$.\n\nThe specific smoothness and regularity conditions required for membership in an RKHS depend on the choice of kernel. For instance, the squared exponential (RBF) kernel induces an RKHS of highly smooth functions, whereas Matérn kernels admit only a finite degree of differentiability. If $\\nabla_y \\log(w^2)$ is too irregular—e.g., non-smooth or not square-integrable under the kernel—it will not belong to $\\mathcal{H}_k$.\n\nIn the case of our robust weight function, the correction term is a nonlinear, rational function of $m(\\mathbf{x})$. Generally, such functions do not lie in the RKHSs associated with common kernels (e.g., RBF, Matérn), unless specific and uncommon conditions are met—such as the RKHS being closed under composition with the given nonlinearity.\n\nA practical approach to ensure $m_w(\\mathbf{x})$ belongs to the RKHS is to project the correction term onto the RKHS. Alternatively, one could use or design a kernel whose RKHS explicitly accommodates this class of nonlinear transformations, although this may come at the expense of flexibility or interpretability.\n\n__We will include a discussion of this condition in the revised manuscript.__\n\n\n--------\nFor space concerns, we use the space available with the rebuttal answer to post experimental results requested in Q1 by reviewer 8Bte:\n\nQ1 8Bte: The effectiveness of the algorithm does not seem to be limited to the Matérn kernel function. It would be beneficial to include some discussion and experiments regarding other common GP kernel functions?\n\nA1 8Bte: We generated a table similar to Table 1 but, experimenting with the RBF kernel in place of the Matern kernel. The results demonstrate the overall superiority of RCaGP both in terms of MAE and NLL, as our proposed method is only outperformed by RCSVGP on 2 out of 8 cases: Boston and Parkinsons with asymmetric outliers. \n\n| Method     | Metric       | Boston (A)         | Energy (A)         | Yacht (A)          | Parkinsons (A)     | Boston (N)         | Energy (N)         | Yacht (N)          | Parkinsons (N)     |\n|------------|--------------|--------------------|--------------------|--------------------|--------------------|--------------------|--------------------|--------------------|--------------------|\n| SVGP       | MAE          | 1.145 ± 0.052      | 1.138 ± 0.040      | 1.235 ± 0.044      | 1.002 ± 0.100      | 0.731 ± 0.052      | 0.908 ± 0.026      | 0.767 ± 0.063      | 0.811 ± 0.078      |\n|            | NLL          | 1.691 ± 0.043      | 1.736 ± 0.037      | 1.732 ± 0.041      | 1.590 ± 0.079      | 1.433 ± 0.053      | 1.430 ± 0.023      | 1.420 ± 0.078      | 1.448 ± 0.070      |\n|            | clock-time   | 0.930 ± 0.280      | 1.170 ± 0.070      | 0.650 ± 0.060      | 0.720 ± 0.260      | 0.940 ± 0.030      | 1.230 ± 0.150      | 0.870 ± 0.240      | 0.670 ± 0.130      |\n| CaGP       | MAE          | 1.060 ± 0.095      | 0.971 ± 0.090      | 1.133 ± 0.092      | 0.955 ± 0.110      | 0.488 ± 0.042      | 0.334 ± 0.024      | 0.512 ± 0.042      | 0.674 ± 0.077      |\n|            | NLL          | 1.611 ± 0.086      | 1.535 ± 0.071      | 1.679 ± 0.088      | 1.556 ± 0.080      | 1.278 ± 0.040      | 1.106 ± 0.011      | 1.263 ± 0.042      | 1.362 ± 0.067      |\n|            | clock-time   | 2.470 ± 0.060      | 3.540 ± 0.080      | 0.640 ± 0.100      | 0.870 ± 0.080      | 2.410 ± 0.060      | 3.540 ± 0.160      | 0.770 ± 0.060      | 0.830 ± 0.100      |\n| RCSVGP     | MAE          | 0.680 ± 0.122      | 0.456 ± 0.041      | 0.733 ± 0.114      | 0.840 ± 0.187      | 0.512 ± 0.074      | 0.558 ± 0.035      | 0.640 ± 0.074      | 0.540 ± 0.083      |\n|            | NLL          | 1.364 ± 0.080      | 1.200 ± 0.026      | 1.381 ± 0.084      | 1.518 ± 0.116      | 1.302 ± 0.063      | 1.300 ± 0.034      | 1.318 ± 0.047      | 1.417 ± 0.065      |\n|            | clock-time   | 1.500 ± 0.090      | 2.560 ± 0.100      | 0.970 ± 0.090      | 0.970 ± 0.090      | 1.530 ± 0.040      | 2.780 ± 0.120      | 1.130 ± 0.080      | 0.720 ± 0.060      |\n| RCaGP      | MAE          | 0.733 ± 0.156      | 0.353 ± 0.053      | 0.651 ± 0.095      | 0.913 ± 0.198      | 0.379 ± 0.050      | 0.326 ± 0.023      | 0.365 ± 0.035      | 0.534 ± 0.078      |\n|            | NLL          | 1.405 ± 0.118      | 1.173 ± 0.024      | 1.336 ± 0.061      | 1.531 ± 0.152      | 1.241 ± 0.046      | 1.144 ± 0.007      | 1.226 ± 0.034      | 1.307 ± 0.068      |\n|            | clock-time   | 2.550 ± 0.060      | 3.590 ± 0.090      | 0.790 ± 0.060      | 0.900 ± 0.030      | 2.470 ± 0.030      | 3.740 ± 0.170      | 0.940 ± 0.070      | 1.070 ± 0.240      |\n| Student-t  | MAE          | 0.962 ± 0.037      | 1.015 ± 0.028      | 0.907 ± 0.061      | 1.025 ± 0.079      | 0.896 ± 0.041      | 0.970 ± 0.025      | 0.815 ± 0.071      | 0.965 ± 0.074      |\n|            | NLL          | 1.619 ± 0.039      | 1.606 ± 0.024      | 1.564 ± 0.067      | 1.658 ± 0.085      | 1.593 ± 0.063      | 1.572 ± 0.026      | 1.559 ± 0.100      | 1.627 ± 0.092      |\n|            | clock-time   | 1.030 ± 0.040      | 1.350 ± 0.070      | 0.790 ± 0.060      | 0.730 ± 0.040      | 1.090 ± 0.030      | 1.460 ± 0.060      | 0.920 ± 0.050      | 0.730 ± 0.010      |"}, {"review_id": "odJTHAvmk2", "reviewer": "Reviewer_jPaV", "summary": "This paper proposes a generalized Gaussian process (GP) model that jointly addresses two practical limitations of GPs in real-world settings: robustness to outliers and computational tractability. The approach leverages generalized Bayesian inference with a robust loss to downweight corrupted observations, and incorporates uncertainty stemming from low-rank approximations to the posterior covariance (termed computational awareness). A prior over potentially corrupted data points is introduced. Hyperparameters are optimised using a generalised evidence lower bound. The method is supported by theoretically results on robustness and conservative uncertainty estimates. Experiments are conducted on four UCI regression tasks and four Bayesian optimization (BO) tasks under various data corruption scenarios.", "questions": "+ **Prior over corrupted observations:** The method relies on a prior over which points are corrupted. Can the authors clarify how this would be constructed in practical scenarios? What are real-world applications where such a prior is meaningful or available?\n+ **Choice of weight function parameters:** The method depends on the parameters $\\beta$ and $c$ in the weighting function $w(x,y)$. How sensitive is the performance to these values, and how should practitioners choose them?\n+ **Robustness definition:** The robustness guarantee in line 192 appears weak, as it allows for any corruption inducing a non-infinite KL divergence. Can the authors clarify the intuition and implications of this result, and compare it to bounded worst-case divergence?\n+ **Experimental realism:** The experimental setup, especially for BO, appears somewhat artificial. Can the authors include additional BO experiments without synthetic corruption and on higher-dimensional or larger-scale tasks where computational limitations are more pressing?\n+ **Statistical significance:** It would improve interpretability if Table 1 clearly indicated which methods are not significantly worse than the best. A paired sample test could help here.\n+ **Figure 1:** The predictive mean of RCaGP appears to poorly fit data at high input values. Can the authors explain this behavior? Is it due to initialization, optimization, or hyperparameters?\n+ **GP hyperparameters for BO:** How were the GP lengthscales initialized for BO? This is particularly important in high-dimensional tasks like Rover and DNA [1]\n+ **Additional BO results:** Could the authors report results for BO without corruption across all tasks, not just Hartmann and Lunar?\n\n[1] Vanilla Bayesian Optimization Performs Great in High Dimensions, Hvarfner et al. 2024", "limitations": "The authors acknowledge the need for a broader evaluation campaign in the conclusion. I see no obvious societal impact of this work that would have been missed by the authors.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 2, "originality": 3, "strengths_and_weaknesses": "### **Quality**\n\nThe paper is of high technical quality. It builds upon and synthesizes recent ideas from the GP literature on robust inference and computational uncertainty aware approximations. The methodology is sound, and the theoretical contributions are well presented. I expect this work to be of interest to the GP community. \n\n+ **Strengths:** Sound methodology and interesting integration of existing techniques.\n+ **Weaknesses:** Experimental evaluation is limited in both scope (only 4 UCI datatsets and 4 BO tasks) and realism (artifical corruption of real-world data). The proposed model is evaluated only under synthetic corruption regimes, raising concerns about the practical relevance of the results. In the BO experiments in particular, the small data regime and low input dimensionality seem insufficient to highlight the method’s strengths in realistic computational or model mismatch scenarios.\n\n### **Clarity**\n\nThe paper is clearly written, well organized, and enjoyable to read. Background and related work are well covered, and the key contributions are clearly articulated. Figures and tables are generally helpful, although there are a few missing items that would improve clarity (e.g., what \"std\" refers to in Table 1, explicit definition of the robust loss in the main text).\n\n### **Significance**\n\nThe paper tackles an important challenge by aiming to make Gaussian processes more robust to outliers and more reliable under computational constraints. These are meaningful goals, especially for real-world applications where data quality and compute resources are often limiting factors. However, I have concerns regarding the practical applicability of the proposed approach in its current form.\n\nMethodologically, the reliance on a prior over corrupted observations assumes some knowledge about which datapoints are likely to be outliers. In many practical scenarios (especially in Bayesian optimization) such prior knowledge may be unavailable or difficult to specify. A clearer discussion of real-world use cases where this modeling assumption is justified would help assess the method’s relevance. Similarly, the need to set additional hyperparameters in the weight function (e.g., $\\beta$, $c$) raises concerns: how sensitive is the method to these choices in practice and how would they be selected?\n\nExperimentally, the current evaluation setup appears somewhat artificial. The benefits of the method are demonstrated mainly through controlled corruptions of existing datasets. While the method indeed performs well under these scenarios, such setups are closely aligned with the assumptions it was designed around, making the improvements less surprising. In contrast, the paper does not convincingly demonstrate performance gains in settings that more closely reflect real-world noise or model mismatch. The BO tasks used are relatively high-dimensional and involve small datasets ($\\leq 10250$ samples), where misspecification and computational approximations of the posterior are unlikely to be a limiting factor. This undermines the motivation for computational awareness in this context.\n\nUltimately, the method may indeed offer substantial benefits in the presence of real-world corruptions or large-scale settings where exact inference is infeasible. However, the current experimental design does not convincingly establish such relevance. A stronger empirical case would require identifying or constructing realistic scenarios where the proposed robustness and approximation-aware inference offer clear advantages over standard GPs. \n\nAdditionally, reporting statistical significance (e.g., by bolding methods that are not significantly worse than the top performer in Table 1) would help provide a more nuanced view of the method's effectiveness.\n\n### **Originality**\n\nThis paper offers a combination of existing ideas that have not yet been jointly explored. The integration of outlier-aware inference with posterior approximation-aware uncertainty modeling is elegant and fills a gap in the GP literature. The contribution is more incremental than novel, but still worthwhile.", "was_revised": true, "final_justification": "I thank the authors for actively participating in the discussions.\n\n**Main criticism:**\n* *Limited experimental setup*: The paper only contains two experiments (regression on 4 UCI datasets and BO on 4 tasks). While the method performs well in these cases, the experiments are rather small-scale (768 samples for UCI and 10,000 for BO) and the corruptions are mostly synthetic, which does not make a strong case for the added computational uncertainty and robustness the method proposes. While the authors made a substantial effort to provide additional results to address this assessment, I believe that the empirical evaluation remains too weak.\n* *Sensitivity to hyperparameters:* The method seems very sensitive to hyperparameter choices, which can be problematic in practice. This needs to be more thoroughly addressed.\n\n**Score:** In light of the above remarks, I maintain my score of 4.\n\n**Recommendations:** I recommend the authors extend the empirical evaluation to considerably larger datasets (such as the power dataset in Wenger et al. [45]) where computational uncertainty is more of an issue, and include cases where misspecification on real-world data is more problematic to demonstrate the usefulness of the additional robustness.", "rebuttal": "Thank you for your detailed and thoughtful review. We appreciate your recognition of the paper’s technical quality, clear presentation, and the relevance of our contributions to robust and computationally aware Gaussian Process modeling.\n\nThe ``1 std'' reported in Table 1 refers to the standard deviation computed for the metrics across 20 train-test splits. The robust loss equation is presented in Equation S2. We will bring this equation to the main text from the Appendix.\n\nQ1: The method relies on a prior over which points are corrupted. Can the authors clarify how this would be constructed in practical scenarios? What are real-world applications where such a prior is meaningful or available?\n\nA1: Thank you for raising this important point. To clarify, the method does not assume a prior distribution over which data points are corrupted. Instead, it relies on a ``centering function'' $m(\\mathbf{x})$, and data points with large deviations from this function—i.e., where $\\lvert y - m(\\mathbf{x})\\rvert$ is high—are downweighted by the model. We considered two choices for $m(\\mathbf{x})$: the GP prior mean, which is a constant mean with a learnable scalar, and an expert-driven prior mean, which is introduced in Section 5 and Appendix G. The latter assumes that a domain expert can flag potential outliers and provide corrections. These corrected points are then used to learn a function $m(\\mathbf{x})$, which informs the weight function $w(\\mathbf{x}, y)$. In this context, the notion of a “prior over corrupted observations” is appropriate, though we emphasize that such an expert-driven mean is not required for the method to be applicable. We believe this setup is relevant in settings where expert annotation is available or can be integrated into the modeling process—such as in clinical data analysis or fault detection in manufacturing. However, as demonstrated in our experiments, the method performs robustly even when using simpler, data-agnostic choices for $m(\\mathbf{x})$, such as a constant or GP prior mean.\nWe will revise the manuscript to include this discussion clarifying the role and interpretation of $m(\\mathbf{x})$, and to better distinguish between optional expert involvement and baseline usage.\n\nQ2: The method depends on the parameters $\\beta$ and $c$ in the weighting function $w(\\mathbf{x},y)$. How sensitive is the performance to these values, and how should practitioners choose them?\n\nA2: In the first paragraph of Section 7.3, we conducted an ablation study on the UCI regression datasets to assess sensitivity to $c$. To recap the conclusions written in the manuscript:\n\n>We evaluate RCaGP's sensitivity to the threshold $c = Q_n(\\vert \\mathbf{y} - \\mathbf{m} \\vert, 1 - \\epsilon)$, defined as the $(1 - \\epsilon)$-quantile of $\\lvert \\mathbf{y} - \\mathbf{m} \\rvert$ and controlling how many values are treated as outliers, by varying $\\epsilon$. Figure 4 (left) shows MAE and NLL averaged over 4 UCI datasets (normalized; per-dataset results in Figure S3. RCaGP consistently achieves lower error across all $\\epsilon$, with best results at $0.15$, suggesting treating roughly 15\\% of the data as potential outliers is optimal. This aligns well with the naturally present and 10\\% injected asymmetric outliers.\n\nThese results suggest that $c$ can be selected based on the expected proportion of outliers.\n\nWe set $\\beta=1$. An ablation study to explore the impact of varying $\\beta$ on the yacht dataset with asymmetric outliers: \n\n| Method | Metric | **β = 0.1**       | **β = 0.5**       | **β = 1.0**       | **β = 1.5**       | β = 2.0         |\n|--------|--------|------------------|------------------|------------------|------------------|----------------|\n| RCaGP  | MAE    | 0.768 ± 0.075    | 0.421 ± 0.072    | 0.435 ± 0.052    | 0.472 ± 0.039    | 0.498 ± 0.035  |\n|        | NLL    | 1.439 ± 0.061    | 1.294 ± 0.055    | 1.253 ± 0.043    | 1.247 ± 0.031    | 1.250 ± 0.026  |\n\nBoth MAE and NLL decrease rapidly as $\\beta$ increases up to 0.5, indicating a significant performance improvement. After, the performance begins to degrade, but the increase in MAE and NLL is more gradual.\n\n__We will explicitly acknowledge hyperparameter tuning as a current limitation of the method.__\n\nQ3: The robustness guarantee in line 192 appears weak, as it allows for any corruption inducing a non-infinite KL divergence. Can the authors clarify the intuition and implications of this result, and compare it to bounded worst-case divergence?\n\nA3: Our theoretical result guarantees that the posterior under corrupted data remains close to the clean-data posterior whenever the corruption leads to a finite KL divergence. While this does not assume a worst-case, norm-bounded corruption, it encompasses a broad and realistic range of outliers that are not well captured by strict norm constraints.\n\nIn contrast to worst-case robustness guarantees—which often rely on conservative bounds under adversarial assumptions—our result shows that as the KL divergence of the corruption decreases (regardless of whether it is norm-bounded), the posterior transitions smoothly toward the clean posterior. This provides a more flexible and practically relevant information-theoretic notion of robustness. We will add a concise summary of these points as a discussion.\n\nQ4: The experimental setup, especially for BO, appears somewhat artificial. Can the authors include additional BO experiments without synthetic corruption and on higher-dimensional or larger-scale tasks where computational limitations are more pressing?\n\nA4: We would like to clarify that the BO tasks considered in our study already fall within the high-dimensional regime. Specifically, the 60D Rover and 180D Lasso DNA benchmarks are commonly used in the evaluation of high-dimensional Bayesian optimization methods. While high dimensionality was not the central focus of our work, we believe these tasks are representative of realistic and computationally demanding BO scenarios. Regarding the scale of the experiments, we interpret \"larger-scale\" as referring to a greater number of optimization steps or oracle evaluations. Extending the experiments in this direction is indeed valuable, but unfortunately not feasible within the constraints of the rebuttal period. That said, we agree this is an important direction for future work and will note it as such in the revised manuscript.\nConcerning the realism of the experimental setup: we acknowledge that the use of synthetic corruption introduces a controlled setting that may differ from real-world noise patterns. However, our goal was to isolate and evaluate the robustness properties of the proposed model under well-defined, repeatable conditions. We agree that incorporating naturally noisy or real-world BO tasks would strengthen the empirical evaluation, and we will make this limitation explicit in the manuscript.\n\nQ5: It would improve interpretability if Table 1 clearly indicated which methods are not significantly worse than the best. A paired sample test could help here.\n\nA5: Thank you for this suggestion. __We will add the results of a paired sample test to the manuscript.__\n\nQ6: The predictive mean of RCaGP appears to poorly fit data at high input values. Can the authors explain this behavior? Is it due to initialization, optimization, or hyperparameters?\n\nA6: It is due to the choice of the mean function. We set the mean function relatively far below the average outliers. In a low-input region, a sufficient amount of data overrides the penalty, resulting in a fitted curve. On the other hand, the high-input region is scarce in data, allowing the mean to act more dominantly. RCSVGP exhibits a similar trend, but more pronounced—likely due to its reliance on variational approximations, which are more sensitive to prior assumptions in data-scarce regions. In contrast, the computational uncertainty in RCaGP appears to buffer the impact of the mean prior, preventing it from dominating excessively in those regions.   \n\nQ7: How were the GP lengthscales initialized for BO? This is particularly important in high-dimensional tasks like Rover and DNA (Vanilla Bayesian Optimization Performs Great in High Dimensions, Hvarfner et al. 2024)\n\nA7: In our experiments, we used the default lengthscale initialization from GPyTorch, which sets a common value of approximately $\\ell \\approx 0.69$. The lengthscale is the same across all dimensions. While BoTorch includes strategies for dimension-aware initialization—as discussed in Hvarfner et al. (2024)—these are not part of GPyTorch’s default behavior, and we did not explicitly implement them in our setup.\nWe agree that evaluating the impact of dimension-aware lengthscale initialization in the context of high-dimensional tasks such as Rover and Lasso DNA would be a valuable addition. However, given the large number of experiments already being run and the time constraints of the rebuttal period, we may not be able to include this analysis at this stage. We will note this point in the revised manuscript as a direction for future work.\n\nQ8: Could the authors report results for BO without corruption across all tasks, not just Hartmann and Lunar?\n\nA8: We provide the results for BO without corruption on the Rover 60D experiment. Lower is better.\n\n| Method   | T = 50         | **T = 100**     | **T = 200**     | **d = 500**     | **T = 1000**     | T = 2000       |\n|----------|----------------|----------------|----------------|----------------|------------------|----------------|\n| CaGP     | -14.54         | -14.24         | -11.20         | -8.95          | -1.22            | -1.04          |\n| RCaGP    | -4.95          | -4.65          | -4.65          | -4.57          | -3.45            | -1.55          |\n| SVGP     | -6.88          | -5.53          | -4.27          | -2.90          | -2.03            | 0.05           |\n| RCSVGP   | -15.21         | -8.69          | -5.14          | -2.86          | -1.37            | 1.1146         |"}, {"review_id": "tPhRiaceXK", "reviewer": "Reviewer_8Bte", "summary": "This paper addresses the intertwined challenges of computational scalability and model robustness in Gaussian Processes when applied to large datasets corrupted by outliers. The authors introduce the RCaGP, whose core contribution is the principled integration of two key ideas: (1) leverages the mechanism of RCGPs to down-weight the influence of outliers, (2) simultaneously incorporating the approach of CaGPs to quantify and compensate for uncertainty arising from computational approximations, thereby preventing overconfident predictions. The paper also provides theoretical guarantees for RCaGP's robustness and conservative uncertainty estimates.", "questions": "(1) The effectiveness of the algorithm does not seem to be limited to the Matérn kernel function. It would be beneficial to include some discussion and experiments regarding other common GP kernel functions?\n(2) Outliers are caused by noise. Is there any discussion regarding the outliers generated by different magnitudes of noise? It seems that GPs perform differently under varying noise levels.\n(3) In the BO experiments, while using 250 points as the initial sample size ensures a fair comparison, this number is relatively large. For instance, 20 points are often sufficient for a 6-dimensional problem. This large initial set might mask the differences in sample efficiency among the models during the early exploration phase?", "limitations": "(1) The paper's core mechanism for achieving robustness through the weight function w(x,y) is critically dependent on a high-quality prior mean m(x). Algorithm identifies outliers by measuring the deviation of an observation y from m(x). This means that if m(x) is poorly specified, the model may erroneously down-weight normal data points or fail to identify true outliers, leading to a degradation in the performance of the entire robustness framework.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n(1) The paper is generally well written and well presented. \n(2) The proposed RCaGP model is not merely a simple superposition of existing methods (RCGP and CaGP), but rather a principled integration within a unified probabilistic framework. This approach of jointly solving the problem surpasses previous work that has tackled these challenges in isolation. \n(3) The paper provides theoretical support for the reliability of RCaGP. The model's robustness is proven through a bounded Posterior Influence Function (PIF), which stands in sharp contrast to the unbounded PIF of CaGP.\n\nWeaknesses:\n(1) Parameters such as the soft threshold c, the low-rank approximation dimension, and the size of the approximate action matrix require manual setting or cross-validation by the user, which could otherwise affect the balance between robustness and approximation quality. \n(2) The paper highlights the model's critical dependence on the prior mean m(x). The success of this mechanism hinges entirely on the quality of m(x). The use of an expert-guided prior is intended to directly address this weakness in model's robustness mechanism by providing a high-quality \"initial guess\" to enhance its core ability to identify and handle outliers. While the proposed expert-guided method is a good solution, it also intro a potential weakness of the model: in the complete absence of domain knowledge, how to set a reasonable prior mean remains an open question. If the domain expert's annotations are biased or insufficient, the constructed mean prior could potentially introduce misleading information. \n(3) Although the paper mentions that only posterior expectations are used for computational efficiency, requiring experts to review and correct a certain proportion of potential outliers in large-scale datasets is a process that is inherently difficult to scale, which seems to conflict with the paper's highlight on High-Throughput applications. \n(4) In the \"Related Work\" section, the paper mentions the first category of robust GP methods, which replace the Gaussian noise assumption with heavy-tailed distributions. Although the authors note that these methods often break conjugacy and require approximate inference, a direct comparison of RCaGP with a representative model from this category (e.g., an SVGP based on a Student-t likelihood) would more forcefully demonstrate whether RCaGP, while maintaining its advantage of conjugacy, still leads in terms of performance.", "was_revised": true, "final_justification": "I am raising my score to 5 as I am satisfied with the authors’ response and the additional experiments conducted. After reviewing the discussion with other reviewers, I believe the paper now presents sufficient content. I encourage the authors to integrate the discussed revisions into the revised manuscript or appendix to further strengthen the work.", "rebuttal": "Thank you for your kind and constructive review. We sincerely appreciate your recognition of the unified design of RCaGP, its theoretical grounding, and the clarity of the presentation.\n\nQ1: The effectiveness of the algorithm does not seem to be limited to the Matérn kernel function. It would be beneficial to include some discussion and experiments regarding other common GP kernel functions?\n\nA1: We generated a table similar to Table 1 but, experimenting with the RBF kernel in place of the Matern kernel. The results demonstrate the overall superiority of RCaGP both in terms of MAE and NLL, as our proposed method is only outperformed by RCSVGP on 2 out of 8 cases: Boston and Parkinsons with asymmetric outliers. The table spans 3000 characters and cannot be posted within this comment, it has been posted in the rebuttal answer associated with reviewer uEzz.\n\nQ2: Outliers are caused by noise. Is there any discussion regarding the outliers generated by different magnitudes of noise? It seems that GPs perform differently under varying noise levels.\n\nA2: We conducted an ablation study by running experiments on the Yacht dataset. We then vary the magnitude of noise by changing the parameter distribution of outliers. Specifically, we vary the upper bound of the uniform distribution and report the MAE and NLL. The settings follow the experiments reported in the main manuscript.\n\n| Method | Metric | **ub = 2**       | **ub = 5**       | **ub = 9**       | **ub = 15**      | ub = 20         |\n|--------|--------|------------------|------------------|------------------|------------------|-----------------|\n| RCaGP  | MAE    | 0.429 ± 0.051    | 0.430 ± 0.051    | 0.435 ± 0.052    | 0.440 ± 0.052    | 0.441 ± 0.056   |\n|        | NLL    | 1.257 ± 0.039    | 1.256 ± 0.042    | 1.253 ± 0.043    | 1.248 ± 0.040    | 1.244 ± 0.040   |\n\n\nThe results indicate that as the magnitude of outliers increases, RCaGP experiences a gradual degradation in performance. Importantly, this degradation is relatively modest, demonstrating the method’s resilience to increasingly severe noise.\n\nQ3: In the BO experiments, while using 250 points as the initial sample size ensures a fair comparison, this number is relatively large. For instance, 20 points are often sufficient for a 6-dimensional problem. This large initial set might mask the differences in sample efficiency among the models during the early exploration phase?\n\nA3: BO outcomes results for the Hartmann6D function with 20 initializations are currently running and will be displayed soon.\n\nQ4 (Limitation): The paper's core mechanism for achieving robustness through the weight function $w(\\mathbf{x},y)$ is critically dependent on a high-quality prior mean $m(\\mathbf{x})$. Algorithm identifies outliers by measuring the deviation of an observation $y$ from $m(\\mathbf{x})$. This means that if $m(\\mathbf{x})$ is poorly specified, the model may erroneously down-weight normal data points or fail to identify true outliers, leading to a degradation in the performance of the entire robustness framework.\n\nA4: In our work, we considered two formulations for $m$: \n\n- A constant mean with a learnable constant.\n\n- An expert-driven prior mean, detailed in Section 5 and Appendix G.\n\nIn Table 1, Figure 3, Figure 4 (left), Table S1, Table S2, Figure S1, Figure S2, and Figure S3, we use only the constant mean for the RCSVGP and RCaGP baselines. These are labeled either simply as “RCSVGP” and “RCaGP” (Table 1, Figure 3, Figure S1–S3), or as “RCSVGP (const.)” and “RCaGP (const.)” (Table S1 and S2). As shown in Sections 7.1 and 7.2, our proposed method RCaGP (const.) yields the best overall performance, even with a constant mean.\n\nThe expert-driven prior mean is evaluated only in a dedicated ablation study (Section 7.3, paragraph two), with results shown in the right panel of Figure 4. There, the baselines are labeled “RCSVGP (exp.)” and “RCaGP (exp.)”. These results show that the expert-driven prior can improve performance in BO. However, this experiment assumes perfect expert correction of all outliers, which is unrealistic in practice. This ablation study should therefore be seen as a proof of concept. In future work, we aim to extend this strategy to scenarios with partial and noisy expert supervision.\n\n__We apologize for any confusion caused by inconsistent labeling and will revise the manuscript for clarity.__\n\nThat said, we agree that using a constant prior mean can lead to degraded results, especially when the target function lies far from the prior. This limitation was discussed in [1], where the authors propose replacing the constant mean with the posterior predictive mean obtained via Generalized Bayes and filtering. Figure 4 in [1] illustrates the significant gains from this modification. We believe that integrating this idea into RCaGP is a promising and natural direction for future improvement.\n\n[1] Laplante et al. Robust and Conjugate Spatio-Temporal Gaussian Processes.\n\nQ5 (4 in Weaknesses): In the \"Related Work\" section, the paper mentions the first category of robust GP methods, which replace the Gaussian noise assumption with heavy-tailed distributions. Although the authors note that these methods often break conjugacy and require approximate inference, a direct comparison of RCaGP with a representative model from this category (e.g., an SVGP based on a Student-t likelihood) would more forcefully demonstrate whether RCaGP, while maintaining its advantage of conjugacy, still leads in terms of performance\n\nA5: We implemented another robust baseline: SVGP with a Student-t likelihood, using Laplace approximation for inference.\nThe results for the UCI regression datasets with asymmetric outliers and without outliers are displayed in a table below. These can be compared to the results provided in the main text (Table 1). In neither the corrupted nor non-corrupted setting is the Student-t baseline competitive with the other methods. \n\n| Method     | Metric       | Boston (A)         | Energy (A)         | Yacht (A)          | Parkinsons (A)     | Boston (N)         | Energy (N)         | Yacht (N)          | Parkinsons (N)     |\n|------------|--------------|--------------------|--------------------|--------------------|--------------------|--------------------|--------------------|--------------------|--------------------|\n| Student-t  | MAE          | 0.962 ± 0.037      | 1.015 ± 0.028      | 0.907 ± 0.061      | 1.025 ± 0.079      | 0.896 ± 0.041      | 0.970 ± 0.025      | 0.815 ± 0.071      | 0.965 ± 0.074      |\n|            | NLL          | 1.619 ± 0.039      | 1.606 ± 0.024      | 1.564 ± 0.067      | 1.658 ± 0.085      | 1.593 ± 0.063      | 1.572 ± 0.026      | 1.559 ± 0.100      | 1.627 ± 0.092      |\n|            | clock-time   | 1.210 ± 0.240      | 1.800 ± 0.400      | 1.350 ± 0.260      | 0.990 ± 0.090      | 1.330 ± 0.050      | 1.330 ± 0.080      | 0.890 ± 0.150      | 0.730 ± 0.020      |\n\n\n\nQ6 (1 in Weaknesses):  Parameters such as the soft threshold $c$, the low-rank approximation dimension, and the size of the approximate action matrix require manual setting or cross-validation by the user, which could otherwise affect the balance between robustness and approximation quality.\n\nA6: We answered the question regarding hyperparameters $c$ and $\\beta$ in Q2 to reviewer jPaV.\nRegarding the alternative hyperparameters of the method, such as the projection dimension, we conducted an additional ablation study to assess the impact of these hyperparameters on our results. We run the experiment on the Yacht dataset using the same settings as the UCI experiment, except now we vary the projection dimensionality for the CaGP and RCaGP models.\n\n| Method | Metric     | **d = 2**       | **d = 5**       | **d = 10**      | **d = 15**      | d = 20         | d = 25         |\n|--------|------------|------------------|------------------|------------------|------------------|------------------|------------------|\n| CaGP   | MAE        | 0.761 ± 0.089    | 0.775 ± 0.078    | 0.785 ± 0.082    | 0.795 ± 0.079    | 0.803 ± 0.080    | 0.815 ± 0.081    |\n|        | NLL        | 1.398 ± 0.073    | 1.395 ± 0.064    | 1.395 ± 0.062    | 1.398 ± 0.060    | 1.401 ± 0.059    | 1.409 ± 0.061    |\n|        | clock-time | 1.210 ± 0.240    | 0.720 ± 0.040    | 1.250 ± 0.110    | 1.430 ± 0.110    | 0.950 ± 0.180    | 1.070 ± 0.030    |\n| RCaGP  | MAE        | 0.426 ± 0.051    | 0.435 ± 0.052    | 0.438 ± 0.059    | 0.443 ± 0.055    | 0.436 ± 0.053    | 0.433 ± 0.047    |\n|        | NLL        | 1.246 ± 0.040    | 1.253 ± 0.043    | 1.250 ± 0.046    | 1.249 ± 0.046    | 1.246 ± 0.044    | 1.247 ± 0.040    |\n|        | clock-time | 1.700 ± 1.610    | 0.870 ± 0.060    | 1.190 ± 0.180    | 1.170 ± 0.060    | 1.230 ± 0.180    | 1.130 ± 0.180    |\n\nFrom the table, we observe that increasing the projection dimensionality does not necessarily reduce the MAE and NLL of RCaGP. However, CaGP gets higher MAE and NLL as we increase the projection dimensionality. This suggests that RCaGP is more robust to changes in projection dimension compared to CaGP."}, {"review_id": "zSepNoL150", "reviewer": "Reviewer_Q5xx", "summary": "The work proposes to combine the idea of computation-aware gaussian processes (CaGP) with robust conjugate GP (RCGP). This way the proposed method allows for conservative uncertainty estimates \nand robustness to outliers. \nExperiments on standard regression tasks and Bayesian optimization tasks suggest that the proposed method leads to better uncertainty estimate with and without outliers when compared to CaGP and RCGP.", "questions": "- It seems the weight function requires that an expert labels the outliers.\nHow are these labels exactly chosen in the experiments of Table 1 and Figure 3?\n\n- How is $\\beta$ of Equation (12) chosen?\n\n- How is the projection dimensionality chosen.\n\n- Clarification of the above questions and additional results with the method from [Ament, 2024] (see Weakness section) would affect my scoring.", "limitations": "no, the choice of all the hyper-parameters is not discussed as a limitation.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strength:\n- Extends theoretical results from CaGP and RCGP to their proposed method.\n- Empirical evaluation showing that proposed method improves uncertainty quantification compared to CaGP and RCGP\n\nWeakness:\n- The proposed method inherits from CaGP and RCGP various hyper-parameters like the projection dimension, choice of mean function, choice of weighting function (= choice of $\\beta$ and $c$).\nFor choosing the mean function, the work assumes that (noisy) expert labels are available that flag potential outliers. This seems to be a rather strong assumption.\nIt is also unclear how the other hyper-parameters are chosen (except $c$ which is explained in the experiments section)\n- RCGP is the only other robust method that is compared.\nHowever, like their proposed method RCGP is known for its sensitivity to mean prior.\nI would find it more convincing, if the work compared also to the method of [Ament, 2024].\n- I guess the projection dimensionality directly affects the computational costs.\nSince one of the main claim is computational awareness, I would like to see how the results of the proposed method change, if less computational resources are used.\nIn particular, in Table 1, currently the proposed method is almost always more computationally expensive than RCSVGP.\n- Proposition 4.2 and Equation (22) are hard to understand. (I guess there is some type in (22)?)\n\nReferences:\n[Ament, 2024]\nAment, Sebastian, et al. \"Robust Gaussian processes via relevance pursuit.\" Advances in Neural Information Processing Systems 37 (2024): 61700-61734.", "was_revised": true, "final_justification": "The main reason for my decision are:\nThe paper is technically sound and several additional experiments including the comparison to [Ament, 2024] were presented in response to my questions, leading me to increase my score.", "rebuttal": "Thank you for your thoughtful feedback. We are grateful for your recognition of the extended theoretical results and the improvements in uncertainty quantification.\n\nQ1: It seems the weight function requires that an expert labels the outliers. How are these labels exactly chosen in the experiments of Table 1 and Figure 3?\n\nA1: This question is reminiscent of Q4 from 8Bte, we paste the answer here:\n>In our work, we considered two formulations for $m$: \n>- A constant mean with a learnable constant.\n>- An expert-driven prior mean, detailed in Section 5 and Appendix G.\n\n>In Table 1, Figure 3, Figure 4 (left), Table S1, Table S2, Figure S1, Figure S2, and Figure S3, we use only the constant mean for the RCSVGP and RCaGP baselines. These are labeled either simply as “RCSVGP” and “RCaGP” (Table 1, Figure 3, Figure S1–S3), or as “RCSVGP (const.)” and “RCaGP (const.)” (Table S1 and S2). As shown in Sections 7.1 and 7.2, our proposed method RCaGP (const.) yields the best overall performance.\n\n>The expert-driven prior mean is evaluated only in a dedicated ablation study (Section 7.3, paragraph two), with results shown in the right panel of Figure 4. There, the baselines are labeled “RCSVGP (exp.)” and “RCaGP (exp.)”. These results show that the expert-driven prior can improve performance in BO. However, this experiment assumes perfect expert correction of all outliers, which is unrealistic in practice. This ablation study should therefore be seen as a proof of concept. In future work, we aim to extend this strategy to scenarios with partial and noisy expert supervision.\n\n>__We apologize for any confusion caused by inconsistent labeling and will revise the manuscript for clarity.__\n\nQ2: How is $\\beta$ (Equation 12) chosen?\n\nA2: We answered the question regarding $\\beta$ in Q2 to reviewer jPaV.\n\nQ3: How is the projection dimensionality chosen.\n\nA3: For both UCI and BO experiments, we first tried several candidates on one task and chose the one that provided the best results. Then, we fixed this dimensionality across all tasks. We acknowledge that selecting the optimum dimensionality remains a problem for the computation-aware GP. in A6 to 8Bte, we also varied the projection dimensionality in a dedicated ablation study.\n\nQ4 (2 in Weaknesses): RCGP is the only other robust method that is compared. However, like their proposed method RCGP is known for its sensitivity to mean prior. I would find it more convincing, if the work compared also to the method of [Ament, 2024].\n\nA4: We did our best to integrate the method proposed by Ament et al. as an additional baseline. Given that all baselines employ a form of approximate inference, the most straightforward way to achieve this is to derive an SVGP model with relevance pursuit, as proposed by the following lemma:\n\nProposition: Let $\\mathcal{D}{\\setminus i} = {(\\mathbf{x}j, y_j): j \\ne i}$ be the dataset excluding point $i$. Let $\\mathbf{Z} = {z_m}$ be the inducing points, and define $\\boldsymbol{\\rho} = \\boldsymbol{\\rho}{\\setminus i} + \\rho_i \\mathbf{e}i$, where $\\boldsymbol{\\rho}, \\boldsymbol{\\rho}{\\setminus i} \\in \\mathbb{R}+^n$, $[\\boldsymbol{\\rho}_{\\setminus i}]_i = 0$, and $\\mathbf{e}i$ is the $i$-th canonical basis vector. Fixing $\\boldsymbol{\\rho}{\\setminus i}$:\n\n$$\\rho_i^\\ast = \\arg\\max_{\\rho_i}  \\mathrm{ELBO}(\\boldsymbol{\\rho}_{\\setminus i} + \\rho_i \\mathbf{e}_i)  = [(y_i - \\mathbb{E}_q[y_i | \\mathcal{D}{\\setminus i} ])^2 - \\mathbb{V}_q[y_i \\vert \\mathcal{D}{\\setminus i}]]{+}$$\n\nwhere $y_i = f(x_i) + \\epsilon_i$. The quantities can be expressed as a function of $\\mathbf{A}^{-1} = (\\mathbf{Q}_\\mathbf{f} + \\mathbf{D}(\\sigma^2 + \\boldsymbol{\\rho}))^{-1}$:\n\n$$\\mathbb{E}_q[y_i \\vert \\mathcal{D}{\\setminus i}] = y_i - [\\mathbf{A}^{-1} \\mathbf{y}]_i / [\\mathbf{A}^{-1}]{ii} \\quad \\text{and} \\quad \\mathbf{V}_q[y_i \\vert \\mathcal{D}{\\setminus i}] = 1 / [\\mathbf{A}^{-1}]{ii}\n$$\n\nwhere $\\mathbf{D}_{\\sigma^2 + \\boldsymbol{\\rho}}$ is a diagonal matrix whose entries are $\\sigma^2 + \\boldsymbol{\\rho}$, and $\\mathbf{Q}\\mathbf{f} = \\mathbf{K}{\\mathbf{x} \\mathbf{z}} \\mathbf{K}\\mathbf{z}^{-1} \\mathbf{K}{\\mathbf{x} \\mathbf{z}}^\\top$.\n\nWe do not include the proof due to space limitations. Indeed, we have updated our manuscript with the full proof of this new class of SVGP. The idea is to consider ELBO by Titsias (2009), which has a closed form. The additional term trace between K and Q does not affect rho during model selection; thus, the derivation follows naturally. Experiments involving this novel baseline are currently running. We will post them as soon as possible.\n\nQ5 (4 in Weaknesses): Proposition 4.2 and Equation (22) are hard to understand. (I guess there is some type in (22)?)\n\nA5: There is indeed a typo in Equation (22): a semicolon and a comma have been wrongly inserted. Thank you for noticing that.\n\nIntuition behind Proposition 4.2: Just as a standard Gaussian Process (GP) accounts for the discrepancy between the true function and its posterior mean via posterior variance—reflecting uncertainty due to limited data—RCaGP introduces an additional layer. Specifically, RCaGP considers the discrepancy between its own posterior mean and that of the exact GP by quantifying computational uncertainty, which arises from approximations made in resource-constrained (limited computation) settings. Consequently, RCaGP captures the total discrepancy between the true function and its posterior mean through a combined uncertainty: the sum of mathematical uncertainty (from limited data) and computational uncertainty (from approximate inference)."}]

Author rebuttal:

Meta-review: This paper presents a technically sound and valuable contribution that merits publication. The reviewers are in unanimous agreement regarding acceptance.

The paper's core strength is its novel and principled integration of robustness to outliers with awareness of computational approximation uncertainty in Gaussian processes, which enables scalability. Jointly addressing these two problems, intertwined in the practical application of GPs, is a significant contribution to the field.

The rebuttal was critical for acceptance, in particular the addition of comparisons of the proposed method against other robust GP approaches (such as Student-t SVGP and a relevance pursuit-based model). In the final version of the paper, the authors are encouraged incorporate these as well as other feedback from reviewers (larger datasets where computational uncertainty is more relevant, cases where misspecification on real-world data is more problematic, etc).
Decision: Accept (poster)