“While correct as far as I can tell, I question the value of the theoretical result, if in practice you optimize the hyperparameters of the covariance function and ρ jointly. Is there reason to believe, that the convexity in ρ for fixed hyperparameters is beneficial given this choice?” “It would be informative to see an experiment that compares the non-convex parametrization to the convex one where all hyperparameters are optimized jointly”

We would expect the overall problem to be “more convex” than with the canonical parameterization, and thus improve optimization performance. More formally, the convexity guarantee will guarantee the positive-definiteness of the submatrix of the Hessian corresponding to the rhos, which is better than not guaranteeing positive definiteness at all.

Positive-definiteness is also beneficial for quasi-Newton optimization algorithms like L-BFGS, which restarts the Hessian-approximation whenever it encounters non-convex regions, because the associated updates to the Hessian approximation are not positive-definite. This leads the algorithm to momentarily revert back to gradient descent, with an associated slower convergence rate.

Practically, we ran convergence analyses using the data from Figure 1, allowing all rhos to be optimized jointly with other hyper-parameters (length scale, kernel variance and noise variance), recording the achieved negative log marginal likelihood (NLML) as a function of the tolerance parameter ftol of the L-BFGS optimizer. The results indicate that the optimizer terminates with a much better NLML using the convex parameterization with the same convergence tolerance:

There are also settings in which we do not actually jointly optimize the hyperparameters. For instance, consider a situation in which we have access to data from the same data generating process that has been manually labeled by domain experts as outlier-free. Then we can estimate the model hyperparameters (of the non-robust GP) on that data, and fix those for the RGP-RP on the new data set that we do not know to be outlier-free (and that we cannot label due to cost or other practical reasons).

“Even though this is quite recent work, I would have liked to see a discussion of Altamirano et al. (2024) in the related work section”

Thanks for pointing out this relevant work. We have added a discussion of and comparison to this work in our general response, and will add this to the paper.

“Or is there an implicit assumption that each evaluation is significantly more expensive than fitting the GP?”

This is indeed the case; GPs are commonly applied to Bayesian optimization or active learning tasks where evaluation times are on the order of hours or days, so the fitting time of the GP is generally not much of a concern. We will make this assumption explicit.

“The experiments on regression problems could be expanded in my opinion, e.g. with benchmark datasets from the UCI repository”

We have done so, see our general response. Please let us know if there are any additional datasets you would like us to consider for the CR if accepted.

“I think the paper could be made stronger if there were at least one experiment where the corruption is not synthetically generated [...]”

Thanks for this feedback. We have added the “Twitter Flash Crash” example from Altamirano et al. (2024) and an additional variant of it to our results (see general response). If you have any other examples of real-world data we would be happy to consider this for inclusion in the CR.

“Almost all plots are way too small and barely readable”

We have adjusted the size of our plots (in particular the axis labels and legends), see attached pdf. We will utilize part of the additional space for the CR submission to increase plot size throughout.

“Is your assumption about the data corruption different from Huber's ε-contamination model or is it identical? Equation 2 suggests they are different, but the text and experiments often refer to only a fraction of the data being corrupted [...]”

A significant high-level difference is that Huber’s model is homoskedastic while ours is not. Citing Huber’s original article, he states: “Let $x_1,⋯, x_n$ be independent random variables with common distribution function $F$”, i.e. Huber’s model assumes that every data point is equally likely to be an outlier.

In contrast, our model introduces data-point specific variances that are adapted using marginal likelihood optimization, leading the resulting likelihood function to be different for each datapoint with a distinct rho.

The article often refers to the “fraction of data being corrupted” because the inference of this fraction is an important problem that determines the performance of RP, which requires choosing a number of outlier variances ($\rho_i$). We solve this problem by leveraging Bayesian model selection.

“How does your approach compare to a vanilla GP on a dataset with no corruption? Does it "fail gracefully”?”

Yes, and the results in our general response demonstrate this. We expect this to hold as long as the prior on the occurrence probability is sufficiently uninformative and puts some weight on no occurrences. Throughout our experiments, we used a geometric prior on the support size with a mean equal to 20% of the number of observations. We will add the results and discussion to the manuscript.

“How did you set the schedule to test data points for being outliers in the experiments (c.f. lines 182-184)?”

See our response to reviewer aaaZ.

“I would have liked to see a bit more discussion on the limitations of the approach in Section 7 [...]”

Good point, we will include these in the discussion in Section 7.

“The weakest part of the paper is the regression problems in the experiment section”

As part of our rebuttal, we have produced a significant number of additional results on UCI benchmarks, the benchmarks and real-world Twitter crash example of Altamirano et al. 2024, and included an additional baseline (see general response and attached pdf).

“Additionally, the proposed method does not clearly outperform other methods, particularly in the Hartmann6 test, where adaptive Winsorization outperforms RGP-RP.”

We don’t necessarily expect our method to outperform all other baselines in all settings. But with our additional empirical results in our general response, we do show that it (i) works much better in situations in which outliers cannot be easily separated from the rest of the data, (ii) is competitive in most other settings. We also note that while the adaptive winsorization heuristic indeed performs slightly better than RGP-RP on Hartmann6 with constant, well-separated outliers, it completely fails on the Friedmann10 function and is barely better than a non-robust GP. In general, winsorization does not perform well if the outliers are not clearly and a-priori separated from the rest of the data.

“The authors do not discuss the results of the regression problems, so a discussion providing guidelines on the cases where the proposed method will perform better is needed”

We agree that the discussion can be improved and we will include this into the paper (see general response for discussion).

“the authors claim their method "permits fast and robust inference"; however, in Appendix D.1, the method is slower than the Student-t GP for the Hartmann6. It would be useful if the authors discussed this in the main paper”

The Student-t likelihood performs very poorly on Hartmann6 and has very high variance. We suspect that there are numerical issues in the optimization of the variational objective that results in the fitting terminating early and at suboptimal points, thus producing faster fit times. Note that the Student-t log likelihood is also non-convex, which can give rise to highly suboptimal local minima. We plan to investigate this in more detail.

“What's the difference between RGP-RP and RGP-RP* in the experiment section?”

See our response to reviewer CPiv.

“What's the schedule K selected for the experiments? How does this affect the performance and fitting time of the method?”

The default schedule used throughout the experiments is [0%, 5%, 10%, 15%, 20%, 30%, 40%, 50%, 75%, 100%] of the number of observations (traversed backwards by the backward algorithm). The fitting time is linear in the number of steps in the schedule, as a new model fit is carried out for each step of the schedule.

The advantage of a schedule with a finer granularity is that the model could be more data efficient, by introducing additional noise variances for outlying data points only, so that all other data is "trusted" up to the homoskedastic noise variance $\sigma^2$. For most practical applications, it is probably sufficient to test for single-digit percentages of outliers, the default schedule here is designed to exhibit high generality, as is evidenced by RP's performance on the variety of benchmarks, including its outperformance on the RCGP benchmarks and Twitter crash example of Altamirano et al. 2024 included in the general response, using the same schedule.

If timing was a particular concern, the optimization could likely be accelerated substantially by keeping track of the approximation of the Hessian matrix generated by each L-BFGS call, which could be used to warm-start the optimization of each step and lead to significantly accelerated convergence, as the consecutive optimization problems tend to be similar.

“While the Friedman and Hartmann6 test functions are well-known, the authors should provide a reference for these functions.”

Thanks for the callout, we will provide the respective references.

“Why use Hartmann6 as a test function for regression? Hartmann6 is typically used to test optimisation due to its properties. Why not use standard UCI benchmarks?”

We have produced additional results on a number of UCI benchmarks as part of this rebuttal (see general response).

“Do the authors know why the student-t GP performs so poorly in the Hartmann6 function, even when the noise comes from a student-t distribution?”

We believe you are referring to the right column of Figure 4, whose x-axis corresponds to the fraction of data points that are corrupted by the outlier process. First, the axis label is missing in the submission, and we will correct this.

What the figure shows is that the Student-t GP only starts to dominate the other methods once more than 40%-50% of data points are corrupted by Student-t noise. We indeed expect the Student-t GP to be the best model as the fraction of Student-t-distributed noise increases, as it becomes the correct noise assumption when 100% of data points are “corrupted” by Student-t noise. However, the results also show: if the main goal is to make a GP robust to a sparse set of outliers – rather than a uniformly heavy tailed observation noise – then RP significantly outperforms the Student-t GP by leveraging this sparse structure.

Kersting et. al. 2007:

This is indeed related work, and we will discuss it in the updated manuscript. The main difference to our method is that their approach is unlikely to work well with outliers. In particular, while they do account for heteroskedastic observation noise, they model the logarithm of the standard deviation of the noise noise process as another GP. This not just restricts the observation to be (heteroskedastic) Gaussian, but, importantly, imposes a smoothness assumption on the variance of the observations as a function of the inputs. In contrast, our approach infers individual additional noise terms, and thus allows us to go beyond the setting in Kersting et. al. 2007 to handle even gross outliers (which require a noise process that is discontinuous in the inputs). Another downside of the approach in Kersting et. al. 2007 is that they require an EM algorithm to estimate the parameters of the two GP models, which does not have any theoretical guarantees and in practice often has convergence issues (see e.g. discussion here).

“I found the section 4.3 difficult to follow”

Thank you for your feedback. We will improve the clarity of our writing in this section in particular by introducing the general framework of Bayesian model selection before stating its instantiation in our work and providing further background references. Please let us know if there are additional edits you would like us to do.

“While this work employs greedy algorithms for GPs, it is not the first to do so”

While greedy algorithms have been applied for many problems in general (outside of GPs), their application to the sparse optimization of noise variances for robustness, which we show permits strong theoretical guarantees, has to our knowledge not been studied and proposed before. Please let us know if you would like us to add any additional reference for related work.

“The experimentation is adequate but less discussed”

We produced a significant number of additional results with new baselines and test problems for this rebuttal, and will improve discussion of the entirety of the empirical results (see general response).

“Figure 1 and 2 are not discussed at all. What regression example?”

Thank you for catching this. The illustrative regression example is based on a synthetic one-dimensional modified sine function with a few outliers. We will detail the exact setup for this in the appendix.

“Is the x axis of right figure in Figure 4 percentage of corruptions or probability?”

Yes, this is stated in the figure’s caption. We will add an axis label as well in order to improve clarity.

"Why Matern kernel was chosen in one of the baselines?”

No particular reason other than that the Matern-5/2 kernel is widely used in practice and the default in many GP / BO libraries (e.g. botorch, trieste). We expect the results to be quite similar if another common kernel, say an RBF kernel, were used.

“What is the difference between RGP-RP*-BW and RGP-RP-BW. I see that the former is canonical one. What is the difference between canonical and non-canonical reparameterization of s?”

This difference is described at the beginning of Section 5.2. The “canonical” parameterization is that which directly uses the robust variances $\boldsymbol \rho$ in the parameterization. The non-canonical parameterization is the convex parameterization we introduce in Section 5.2, i.e., that where ${\boldsymbol \rho}(\mathbf s) = \text{diag}({\bf K_0}) \odot ((1 - {\mathbf s})^{-1} - 1)$ and the new parameter is $\mathbf{s}$. The method label coding is as follows: A “*” indicates using the canonical (not necessarily convex) parameterization, and “BW” indicates using the backward algorithm for relevance pursuit (Algorithm 2 in Appendix A). We will improve the clarity of the exposition here and clearly explain the method labels. The convex parameterization is a key ingredient that allows us to prove strong theoretical guarantees for the result of the greedy algorithm, and also has beneficial effects on the convergence of the hyper-parameter optimization algorithms, as we demonstrate in the other responses.

“What is $\bf K_0$ on line 218?”

Thanks for catching this forward reference, $\mathbf{K}_0$ is defined in Theorem 8: $\mathbf{K}_0 = k(\mathbf{X}, \mathbf{X}) + \sigma^2 \bf I$. We’ll move the definition up in the text.

Thank you again for your suggestions. We present the additional results you suggested below and hope that this will allow you to fully recommend our work for acceptance.

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

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

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

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

RMSE

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

We highlight the main limitation of Relevance Pursuit versus the additional baselines: if 100% of the data points are subject to heavy-tailed noise, the methods based on uniformly heavy-tailed likelihoods, including Huber + Projection will perform best, as the following results demonstrate, similar to the setup in Table 1 of Algikar et.al., 2023.

100% Laplace Noise

RMSE

NLP

Last, we provide the code of the main parts of our translation of the Matlab code of Algikar et.al., 2023 in the following, for you to be able to check its correctness.

def projection_statistics(H: Tensor) -> Tensor:
    """

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

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

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


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

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

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

For the Huber loss, we use $\epsilon = 0.45$ and $b = 0.5$, identical to the setting of the Github repository of Algikar et.al., 2023. All reported results were generated with 32 independent replications, reporting the mean and standard error.

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

A word of thanks Thank you again for your valuable suggestions, we believe they make the paper stronger. Let us know if you have additional questions or if all your remaining concerns have been addressed.