Practical and Rigorous Extremal Bounds for Gaussian Process Regression via Chaining

Thank you very much for your comments. We address them below.

Weaknesses:

1. Response:

Thank you for your comment. We acknowledge that our initial writing may not have clearly conveyed the relevance of Theorem 1 to the GPR bound problem. In the revised paper, we have introduced a general bound (lines 205–215) applicable to all kernels immediately after Theorem 1. For specific kernels, such as RBF and Matérn, tighter estimates of $\mathbb{E}\left[ \sup\_{t \in T} |X\_t - X\_{t\_0}| \right]$ are derived, leading to tighter bounds, which we present along with their corresponding proofs in the subsequent subsections.

Theorem 2's general bound is proven using Theorem 1, highlighting its importance. The detailed proof is provided in the appendix (lines 985-1054). Due to space constraints, we could only include this explanation in the appendix but not the main paper. However, we would try adding it in the camera-ready version if extra space is provided.

2. Response:

Thank you for your comment. In the original manuscript, $t\_0$ was introduced in Section 4.4. In the revised paper, we have clarified its role earlier, directly following Theorem 2 (lines 204–210). Specifically, $t\_0$ is introduced to establish the bound

Due to the zero-mean property of the GP and the symmetry of the covariance function, $t\_0$ is chosen such that $X\_{t\_0}$ is close to zero; otherwise, $\mathbb{E}\left[ \sup\_{t \in T} |X\_t - X\_{t\_0}| \right]$ may be overestimated. A more detailed explanation of this choice and its implications is provided in the revised manuscript, Section 4.4 (lines 386–392).

3. Response:

Thank you for your comment. Theoretically, our method relies solely on the kernel function to compute distances, enabling it to handle both extrapolation (forecasting) and interpolation (infill) tasks. This is because the kernel quantifies similarities between points based on their relative positions, making the method independent of whether the test points lie within or outside the observed range.

We have clarified this in the revised paper (lines 861-891) and included an interpolation experiment using NOAA data to demonstrate this capability. In this experiment, the middle 70% of the data was used as the test set (interpolation points), while the leftmost and rightmost 15% were used as the training set. The results show that the predicted bounds successfully cover all test points, achieving a PICP of 1.0 and a CWC of 1.65, further supporting the effectiveness of our method for interpolation scenarios. For more details, please refer to Appendix B.5 (lines 861-891) of the revised paper.

A minor issue. Response:

Thank you for pointing this out. The font issue has been corrected, and the updates can be found in the revised manuscript at lines 444-447.

Questions:

1. Response:

Thank you for the question. For the general Matérn kernel with arbitrary smoothness, our proposed method provides an applicable bound. As mentioned in our response to W1, we have added a general bound after Theorem 1 (lines 204–210 of Theorem 2) in the revised manuscript, which is applicable to any kernel, including Matérn kernels with $\nu = 1/2$ and $\nu = 5/2$ that you mentioned. The bounds derived for RBF and Matérn $\nu = 3/2$ are tighter than the general bound in Theorem 1 and are specifically tailored for these kernels. For other cases, such as the popular $\nu = 1/2$ and $\nu = 5/2$ of Matérn kernel, the general bound can be applied.

2. Response:

Thank you for the question. Theoretically, the domain dimensionality does not directly affect the bound in our method, as GPR relies on kernel functions to compute distances. This is clarified in the revised paper (lines 254–256), where the formula for $d(s, t)$ is given:

Regardless of the dimensionality of the domain, the kernel function reduces these computations to a single scalar value, ensuring that the bound remains unaffected by the input dimension.

Thank you very much for your comments. We address them below.

Weaknesses:

1. Response:

Thank you for your thoughtful comment. We sincerely apologize for the lack of clarity in our explanation, which may have caused some misunderstanding.

To clarify, this paper focuses on extremal bounds (upper and lower bounds) rather than confidence intervals. Mathematically, extremal bounds are defined by the infimum and supremum, representing absolute limits that fully encapsulate all possible outcomes in the prediction space without any probabilistic interpretation.

While the baselines do not directly compute extremal bounds, they calculate error bounds, which are used to construct upper and lower bounds by adding and subtracting the error from the predictions. These are also deterministic bounds, providing absolute limits rather than probabilities.

Thus, the upper and lower bounds can be understood as corresponding to $\alpha = 0$ offering absolute guarantees instead of probabilistic coverage.

2. Response:

Thank you for your comment, and we apologize for the misunderstanding caused by unclear wording. Similar to the previous question, this paper focuses on extremal bounds (upper and lower bounds) rather than confidence intervals. Because these are deterministic bounds, the ideal outcome is for the prediction interval to cover all points, making a higher PICP preferable. As described on lines 444-447 of the revised paper, PICP directly reflects the proportion of points covered by the bounds. Additionally, as shown in the CWC formula on line 461, when $\mu = 1$, a PICP closer to 1 leads to a smaller CWC, demonstrating better uncertainty calibration.

3. Response:

Thank you for highlighting this point. We agree that uncertainty estimates and statistical significance are crucial for rigorous evaluation. To address this concern, we conducted paired t-tests to evaluate the statistical significance of differences in the CWC metric across three datasets: Boston House Price, NOAA Weather, and Sarcos. In each dataset, the training and testing sets were randomly sampled 100 times, and the models were trained and tested independently in each trial, resulting in 100 CWC values for each model.

As lower CWC values are preferable, negative t-statistics indicate that our model outperformed the compared baselines. The results confirm that our approach achieves statistically significant improvements over the baselines, with $p < 0.01$ in all comparisons. These analyses have been included in the revised manuscript and are detailed in the appendix (lines 813–859).

Table 1: Paired t-Test Comparisons of Our Method against Baselines on the Boston House Price Data. (** indicates $p < 0.01$; * indicates $p < 0.05$; negative t-statistics indicate that our model performs better than the compared model, as lower CWC values are preferable.)

Table 2: Paired t-Test Comparisons of Our Method against Baselines on the NOAA Weather Data. (** indicates $p < 0.01$; * indicates $p < 0.05$; negative t-statistics indicate that our model performs better than the compared model, as lower CWC values are preferable.)

Table 3: Paired t-Test Comparisons of Our Method against Baselines on the Sarcos Data. (** indicates $p < 0.01$; * indicates $p < 0.05$; negative t-statistics indicate that our model performs better than the compared model, as lower CWC values are preferable.)

(Continued in the next comment.)

Thank you for your insightful comments and valuable feedback. We have carefully considered the questions and concerns you raised and have addressed them in our revised submission and responses. We would greatly appreciate it if you could review our responses and let us know if they sufficiently address your points. Should you have any remaining questions or additional suggestions, we would be happy to address them to further enhance the paper.

Thank you very much for your comments. We address them below.

Weaknesses:

[W1]. Response:

Thank you for your positive feedback on the background and related work in our manuscript. We recognize that our earlier writing may have lacked coherence, potentially hindering readability and making it challenging for readers to follow the flow of our explanations.

In the revised paper (uploaded), we have focused on improving logical continuity. Specifically, after introducing Talagrand’s Theorem 1, we now present a general bound earlier (lines 200-215 in the revised manuscript). This ensures that the bounds are clearly introduced early in the discussion, improving readability. The general bound applies to all GPR kernels. However, when restricted to specific kernels such as RBF and Matérn, tighter estimates of $\mathbb{E}\left[ \sup\_{t \in T} |X\_t - X\_{t\_0}| \right]$ can be derived (as discussed in Sections 4.2 and 4.3), leading to tighter bounds for these cases.

Regarding the RBF and Matérn kernels, we agree that they are indeed well-known. However, we included a brief reintroduction to ensure clarity for readers who may not recall their details. Following this, we derived specific bounds for these kernels, which are tighter compared to the general bound, along with the corresponding proofs.

[W2]. Response:

We should have paid more attention to the tone of our writing. We will replace "stellar" with a more neutral "strong".

The wave-style behavior you mentioned in Figure 1c stems from the inherent fluctuations within the Sarcos dataset. As illustrated in the figure, all four methods exhibit similar fluctuations patterns, aligning closely with the black points that represent the actual test data. This indicates that our method is consistent with other SOTA methods in capturing the fluctuation of the dataset, rather than being different.

[W3]. Response:

Thank you for your insightful question. We will address this in two parts.

First, regarding the limitations of GPR, the computational cost of $O(N^3)$ is a fundamental constraint inherent in GPR training, and is not specific to our bounding method. This complexity stems from the necessity of decomposing and inverting the covariance matrix during training. It is important to emphasize that this limitation is intrinsic to the training of GPR models and thus is borne by all bounding methods that work with GPR, including $both$ the baselines and our approach.

Second, we consider the computational complexity of solving the bounds in our method, and the corresponding steps are outlined in the pseudocode in the revised paper (lines 400–420). Constructing $\{T\_n\}$ requires $\mathcal{O}(|D\_{\text{train}}|^2 \cdot \log \log |D\_{\text{train}}|)$ time, which is dominated by kernel distance computations. Computing bounds for $|D\_{\text{test}}|$ test points has a time complexity of $\mathcal{O}(|D\_{\text{test}}| \cdot |D\_{\text{train}}| \cdot \log \log |D\_{\text{train}}|)$. The total computational complexity depends on the relative sizes of the training and test sets. Since the sizes of the training and test sets can vary, the overall complexity is determined by the more computationally intensive step. Thus, the total time complexity is: $\mathcal{O}(\max(|D\_{\text{train}}|^2 \cdot \log \log |D\_{\text{train}}|, |D\_{\text{test}}| \cdot |D\_{\text{train}}| \cdot \log \log |D\_{\text{train}}|))$.

Please note that the papers of the baseline methods do not report their computational complexity. Thus, for the purpose of comparison, we conducted runtime experiments across datasets of varying sizes, with detailed results presented in the appendix of revised paper (lines 627–671).

(Continued in the next comment.)

Thank you for your insightful comments and valuable feedback. We have carefully considered the questions and concerns you raised and have addressed them in our revised submission and responses. We would greatly appreciate it if you could review our responses and let us know if they sufficiently address your points. Should you have any remaining questions or additional suggestions, we would be happy to address them to further enhance the paper.