Robust and Computation-Aware Gaussian Processes

Thank you for your thoughtful feedback. We are grateful for your recognition of the extended theoretical results and the improvements in uncertainty quantification.

Q1: 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?

A1: This question is reminiscent of Q4 from 8Bte, we paste the answer here:

In our work, we considered two formulations for $m$:

• A constant mean with a learnable constant.
• An expert-driven prior mean, detailed in Section 5 and Appendix G.

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.

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.

We apologize for any confusion caused by inconsistent labeling and will revise the manuscript for clarity.

Q2: How is $\beta$ (Equation 12) chosen?

A2: We answered the question regarding $\beta$ in Q2 to reviewer jPaV.

Q3: How is the projection dimensionality chosen.

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

Q4 (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].

A4: 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:

Proposition: 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}$:

$\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}]]{+}$

where $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}$:

where $\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$.

We 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.

Q5 (4 in Weaknesses): Proposition 4.2 and Equation (22) are hard to understand. (I guess there is some type in (22)?)

A5: There is indeed a typo in Equation (22): a semicolon and a comma have been wrongly inserted. Thank you for noticing that.

Intuition 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).

We have updated our SVGP with relevance pursuit formulation. Specifically, we update the trace term of the ELBO formulation since it also depends on $\mathbf{D}_{\boldsymbol{\rho} + \sigma^2}$.

In the proof, we show that the trace difference is scaled by

$\mathbf{K}_{ii} - [\mathbf{Q}_{\mathbf{f}}]_{ii}$

Furthermore, the second root of $\beta_i$ will be negative or zero, depending on

$\mathbf{K}_{ii} - [\mathbf{Q}_{\mathbf{f}}]_{ii}$

In both cases, we reject this beta, leaving the positive case only. Next, we report the results for SVGP with relevance pursuit. We choose 100 inducing points, exactly the same as RCSVGP and SVGP. For the $\mathcal{K}$ schedule, we choose $\{ 0.05n, 0.05n \}$, where $n$ denotes the number of observations and there is $0.1$ fraction of outliers. For fair comparison, we optimize $\boldsymbol{\rho}$ with the Adam optimizer with 50 iterations using $0.01$ learning rate. Other settings exactly follow the main experiment on UCI dataset.

$\frac{1}{\sqrt{\mathbf{K}_{ii} - [\mathbf{Q_f}]_{ii}}}$

These results, to be compared with the left part of Table 1 (asymmetric outliers), show that RCaGP outperforms SVGP with relevance pursuit on all datasets, except Energy.

Thank you so much for the clarifications and the additional comparison to [Ament, 2024]. Finally, could you please comment on the computational costs as I wrote in my previous review: "I guess the projection dimensionality directly affects the computational costs. Since 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. In particular, in Table 1, currently the proposed method is almost always more computationally expensive than RCSVGP."

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.

Q1: 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?

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

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

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

The 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.

Q3: 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?

A3: BO outcomes results for the Hartmann6D function with 20 initializations are currently running and will be displayed soon.

Q4 (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.

A4: In our work, we considered two formulations for $m$:

• A constant mean with a learnable constant.

• An expert-driven prior mean, detailed in Section 5 and Appendix G.

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, even with a constant mean.

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.

We apologize for any confusion caused by inconsistent labeling and will revise the manuscript for clarity.

That 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.

[1] Laplante et al. Robust and Conjugate Spatio-Temporal Gaussian Processes.

Q5 (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

A5: We implemented another robust baseline: SVGP with a Student-t likelihood, using Laplace approximation for inference. The 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.

Q6 (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.

A6: We answered the question regarding hyperparameters $c$ and $\beta$ in Q2 to reviewer jPaV. Regarding 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.

From 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.

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.

The ``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.

Q1: 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?

A1: 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. We 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.

Q2: 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?

A2: 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:

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.

These results suggest that $c$ can be selected based on the expected proportion of outliers.

We set $\beta=1$. An ablation study to explore the impact of varying $\beta$ on the yacht dataset with asymmetric outliers:

Both 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.

We will explicitly acknowledge hyperparameter tuning as a current limitation of the method.

Q3: 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?

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

In 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.

Q4: 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?

A4: 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. Concerning 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.

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

A5: Thank you for this suggestion. We will add the results of a paired sample test to the manuscript.

Q6: 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?

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

Q7: 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)

A7: 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. We 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.

Q8: Could the authors report results for BO without corruption across all tasks, not just Hartmann and Lunar?

A8: We provide the results for BO without corruption on the Rover 60D experiment. Lower is better.

Thank you for your answers to my questions and for all the clarifications.

"A4: [...] Concerning 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."

While I understand the value of controlled experiments with synthetic corruptions, the absence of evaluation on real-world data with naturally occurring misspecification and outliers remains a critical limitation of the paper in its current form. Real-world robustness is essential for demonstrating the practical utility of the proposed method. Simply acknowledging this limitation in the manuscript is insufficient to establish the method's effectiveness in realistic scenarios where such robustness properties are most needed.

"A5: Thank you for this suggestion. We will add the results of a paired sample test to the manuscript."

To facilitate an informed assessment during the discussion period, could the authors provide the paired sample test results? This would enable to evaluate the statistical significance of the reported improvements and potentially update my review accordingly.

"A7: 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. We 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."

Thank you for the clarification. However, this reveals an experimental design issue that significantly undermines the validity of your high-dimensional results (Rover and DNA tasks). Recent literature [1,2,3] has consistently demonstrated that careful lengthscale initialization is critical for Bayesian optimization performance in high-dimensional spaces. The default GPyTorch initialization is inappropriate for high-dimensional problems because it fails to account for the curse of dimensionality, specifically the increased Euclidean distances that render stationary covariance functions ineffective at capturing meaningful correlations.

In high-dimensional discrete spaces, this improper initialization typically reduces the Gaussian process surrogate to essentially random point selection, as the model cannot extrapolate meaningfully to nearby points. Given that the Rover and DNA tasks operate in 60D and 180D spaces respectively, the reported results on these benchmarks are likely unrepresentative of the true performance of the compared methods and may be dominated by this initialization artifact rather than reflecting the actual methodological differences.

This is not merely a minor technical detail but a core experimental validity issue that calls into question the conclusions drawn from the high-dimensional experiments. I strongly recommend either re-running these experiments with appropriate dimension-aware lengthscale initialization (for example initializing as per [2]) or removing the high-dimensional results entirely, as they currently provide little insight into the relative merits of the proposed approach.

[1] Vanilla Bayesian Optimization Performs Great in High Dimensions, Hvarfner et al. 2024.
[2] Standard Gaussian Process is All You Need for High-Dimensional Bayesian Optimization, Xu et al. 2025.
[3] Understanding High-Dimensional Bayesian Optimization, Papenmeier et al. 2025

Thank you for this extensive answer.

We agree on the absolute necessity to run Bayesian Optimization experiments again with a correct lengthscale initialization strategy, specifically for high-dimensional examples like Rover 60D or Lasso DNA 180D. These experiments are running at present, using the initialization proposed by Hvarfner et al. Results will be posted as soon as possible. This being said, it is worth noticing that the recent works you shared all deal with exact GPs and not custom models like SVGPs. While we expect these findings to largely transfer to approximate GP models such as SVGPs, some differences may arise due to the structural and inference-specific distinctions between exact and sparse variational GPs. Furthermore, the existing works on high-dimensional BO do not consider the presence of outliers. Finally, we consider a non-trivial EULBO acquisition function, where we need to simultaneously optimize the query, variational parameters, and the hyperparameters, which may lead to different behaviors.

Next, we present in the table below the statistical significance between methods. We used Wilcoxon signed-rank test, computed across 20 train/test split folds. This table is associated with the left part of Table 1, on UCI regression datasets with asymmetric outliers.

Finally, experiments on real-world data will also follow, as soon as we can.

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.

Q: 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_{noise} 2 \frac{d}{dg(x)}\log w(x, g(x))$?

The 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.

A: 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 $\nabla_y \log(w^2) \in \mathcal{H}_k$

This holds when $\nabla_y \log(w^2)$ has a finite RKHS norm; that is, $\| \nabla_y \log(w^2) \|_{\mathcal{H}_k}^2 < \infty$.

The 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$.

In 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.

A 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.

We will include a discussion of this condition in the revised manuscript.

For space concerns, we use the space available with the rebuttal answer to post experimental results requested in Q1 by reviewer 8Bte:

Q1 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?

A1 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.

Thank you for your response. It appears to me that the assumption in Proposition 4.2 is significantly stronger than the assumptions made by Wenger et al. (2024) and Kanagawa et al. (2018). In the case of the CaGP's results presented in Wenger et al. (2024), the predictive variance represents the worst-case error of the posterior mean, under the condition that the true function lies within the RKHS of the data-generating process's covariance. This seems like a reasonable assumption when the model is well-specified.

However, for the RCaGP, the predictive variance represents the worst-case error when the true function, minus the shrinkage term, belongs to the RKHS of the data-generating process's covariance. I find this somewhat puzzling because RCaGP, as a generalized Bayesian method, is predicated on the idea that standard Bayesian models are misspecified. This assumption makes it difficult to contextualize these results in practice, as it implies that the RCaGP itself is somehow well-specified. Given this context, what would you consider to be the key takeaway from Proposition 4.2?