Optimal kernel regression bounds under energy-bounded noise

We thank the Reviewer for the careful reading of the paper.
On modeling noise sequences as members of deterministic function spaces: We agree that the choice of giving a stochastic or deterministic representation of the noise depends on the problem at hand. There are indeed situations in which stochastic modeling is preferable -- e.g., when samples are i.i.d. with known distribution. However, the choice of modeling noise sequences as energy-bounded members of a RKHS (which is also performed in $Kanagawa+2018$ ) allows us to overcome two main difficulties of stochastic representations, namely dealing with biases and correlation among samples. Additionally, such a set-up allows us to deal with adversarial noises (in fact, we are computing the bound accounting for the worst-case noise realization), and rich-enough hypotheses classes for the noise allow us to compensate for model misspecification of the latent function. Therefore, the presented modeling choice has practical advantages that go beyond tractability. Finally, we would like to point out a further connection between our set-up and the stochastic one (cf. footnote 2, p.8): if one considers sub-Gaussian and i.i.d. noise, it is possible to find a certain probability level for which the norm-bounds in Assumption 1 hold -- but the converse is not true.
On energy-boundedness: While the reviewer is correct that an RKHS-norm assumption is significantly more restrictive than a bound on the $\mathscr{L}^2$ norm, we would like to clarify that we are never invoking energy bounds or $\mathscr{L}^2$ norms on continuous functions. When introducing the RKHS condition in Assumption 1, we directly reference Section 4.1 regarding its interpretation as energy-bounded noise. In particular, for Dirac noise kernels, Equation (16) in Section 4.1 shows that the condition is equivalent to bounded energy of the finite noise sequence affecting the measurements, i.e., $\sum_{i=1}^N w_i^2 \leq \Gamma_w^2$. Such energy-bounded noise assumptions are common in classical robust regression works such as

$Fogel 1979$ , $Bertsekas+2003$ and many others. We see that the word "energy" might also be interpreted as a $\mathscr{L}^2$-energy bound and we plan to use a more precise wording to highlight that this is a bound on the finite and discrete noise sequence that affects our finite and discrete measured data $(x_i,y_i)$, $i=1,\dots N$. Our general theoretical results are valid for any positive definite and uniformly bounded kernel $k^w$ characterizing the noise. This general modeling framework ensures that we can simultaneously recover results from classical energy bounded regression $Fogel 1979$ , connect our results to to recent kernel error bounds $Kanagawa+2018$ , and provide an interpretation of the bound in terms of mean and covariance of a GP (cf. Section 4).
On the claim in the contributions: As the Reviewer correctly points out, our sentence "we derive an analytical solution that exactly characterizes the worst-case realization within the function hypothesis class" is not precise. Indeed, given the non-parametric nature of the problem, the worst-case function $f^\star$ depends on the query point $x_{N+1}$ -- thus, the whole function $\bar{f}(x_{N+1})$, which takes the evaluation of different worst-case functions at different query points, is no longer in the hypothesis space. What we mean is that, for any given query point $x_{N+1}$, our procedure returns a worst-case latent function $f^\star$ in the RKHS space, cf. its parametrization in terms of the weight vector $\theta^\star$ in Eqs. (C.3), (C.16) for Cases 1 and 2. The evaluation of this function at the point $x_{N+1}$ then returns the error bound $\bar{f}(x_{N+1})$. We plan to better convey this idea by moving some of the expressions of the worst-case function $f^\star$ from the Appendix to the main text.
On the benefits of the proposed approach and application on a downstream task: First, we would like to emphasize that if the distribution of the noise is known, probabilistic bounds are the most natural choice to quantify the uncertainty -- and this would be the case for Gaussian i.i.d. noise. However, as soon as the distribution is unknown (e.g., biased and correlated), applying GP bounds would return wrong results, possibly endangering safety for downstream tasks. This displays one of the strengths of the proposed approach: it is distribution-free, and is also capable of dealing with adversarial noises. We defer to footnote 2 for a discussion on the connection between the proposed bound and the probabilistic ones, but we plan on expanding on this point in the paper.
Regarding the optimization of $\sigma$, we highlight that this is only needed to obtain the least conservative error bound, while any value $\sigma$ provides a valid estimate that can be safely used for downstream tasks, as discussed at the end of Section 3.2. In other words, optimizing $\sigma$ leads to a refinement of the uncertainty envelope, and the termination of iterative optimization procedures can be dictated by the downstream task at hand.
To highlight these aspects and also display the performance of the proposed approach, we conducted an experiment where we apply the proposed bounds to a downstream task. Specifically, we consider safety-critical control of an uncertain dynamical systems, where we demonstrate significant benefits of the reduced conservatism and manageable computation times when jointly optimizing $\sigma$. Additionally, we highlight another key property of the proposed bound -- they still hold when considering a subset of data -- which has a significantly bigger impact on the computational demand compare to the scalar $\sigma$.
Experimental details: We compare the proposed robust bound and a probabilistic bound akin to the setup in Section 4.3. The downstream application of the bound is thereby given by an optimization-based controller (specifically, a control barrier function) that ensures safety (constraint satisfaction) of an unknown dynamical system modeled using kernel regression with $10^2$ data points. The only difference in the application of the probabilistic bound, denoted by (a), compared to the proposed robust bound, (b), is the simultaneous optimization over the noise parameter $\sigma$ (one additional optimization variable), as well as the different scaling factor $\beta_\sigma$ of the covariance matrix, see also the discussion in Section 4.3. Additionally, we applied the proposed bound to a subset of data containing the 10 nearest training input locations and we denote this approach by (c). Notably, the derived robust bounds remain valid under arbitrary data selection rules, while the probabilistic approach relies on martingale bounds, which would fail under such data-selection rules. The following table compares the three approaches in terms of their computation time to solve the optimization problem to convergence, as well as the conservativeness of the bound, indicated by the share of the input domain for which a safety-preserving (feasible) input is computed:

The larger feasible domain of the proposed bound (b) illustrates the reduced conservatism of the bound. While, for the same number of data points, the computation time for optimizing the robust bound (b) is visibly larger than for the probabilistic one (a), the subset selection mechanism (c) leads to a similarly large feasible domain as the probabilistic bound at significantly reduced computation times.
On the sensitivity of the bounds to the choice of the RKHS: Choosing the RKHS is a widely investigated topic in the statistical learning literature, and connects with results obtained in the realm of Gaussian process regression, cf.  $Rasmussen+2006, Sec. 5.4$ . Our set-up falls also into this framework: indeed, a key result of our work is that the worst-case error bound has the same formula as the error bound for standard Gaussian process regression for a suitably chosen noise $\sigma$ when choosing a Dirac kernel for the noise, see Equation (8). Although less prevalent, we note that also general (non-Dirac) noise kernels $k^w$ are leveraged for Gaussian process regression

$Rasmussen+2006, Sec. 5.4.3$ and the equivalence (Eq. (8)) equally applies in this case. Hence, the sensitivity of our bound to the choice of the RKHS is comparable to the sensitivity of general Gaussian process regression to the choice of the kernel. For this reason, we see no need to separately investigate the sensitivity to related hyperparameters.

References:
  $Kanagawa+2018$ M. Kanagawa, P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur, "Gaussian Processes and Kernel Methods: A Review on Connections and Equivalences," arXiv:1807.02582, 2018;
  $Fogel 1979$ E. Fogel, "System identification via membership set constraints with energy constrained noise," IEEE Transactions on Automatic Control, vol. 24, no. 5, pp. 752--758, 1979;
  $Bertsekas+2003$ : D. Bertsekas and I. Rhodes, "Recursive state estimation for a set-membership description of uncertainty," IEEE Transactions on Automatic Control, vol. 16, no. 2, pp. 117--128, 1971;
  $Rasmussen+2006$ : C. E. Rasmussen and C. K. I. Williams, Gaussian processes for machine learning. in Adaptive computation and machine learning. Cambridge, Massachusetts: MIT Press, 2006.

We thank the Reviewer for the careful reading of the paper.
On knowing the RKHS norm bound: This assumption is common with other state-of-the-art bounds (e.g.

$Srinivas+2012,Chowdury+2017,Fiedler+2021$ ) that have been successful in downstream tasks such as safe reinforcement learning and Bayesian optimization. Given the practical impact of such results, there exist works that estimate RKHS norms: see e.g. $Tokmak+2025$ and references therein for extrapolation-based approaches, and $Csaji+2025$ for the case of Paley-Wiener kernels of band-limited functions using the uniformly-randomized Hoeffding's inequality or the empirical Bernstein inequality. Another strategy for estimating the RKHS bounds can be developed by expanding on ideas from the robust regression literature, such as $Calliess+2020$ : one starts from the smallest bound that is not falsified by the data (i.e., such that the problem remains feasible) and increases it with a small proportional buffer according to the so-called "lazy adaptation" technique; then, such a bound can be combined with recent RKHS-norm extrapolation heuristics  $Tokmak+2025$ . However, a formal analysis of this technique to the considered setting and intensive numerical validations would be required. We plan to expand on this point in the limitation section of the paper.
Sensitivity of result to norm-bound misspecification: A thorough quantitative analysis depends on the adopted estimation strategy and will be deferred to future works. Nevertheless, we can mention the following cases:

• if the norm bounds are too small, this could be detected by checking feasibility of the constraints in the problem formulation (C.1);

• if the problem (2) remains feasible despite $\Gamma_f$ being too small, then the noise function given by a positive definite noise kernel $k^w$ can generally compensate for the inadequacy of the constraint on the latent function norm.

Finally, we would like to point out that taking conservative (too large) norm bounds for both the noise and the latent function does not have a dramatic impact on the resulting bound. In fact, as displayed in Equation (7), conservative values for $\Gamma_f^2$ and $\Gamma_w^2$ lead to a just sublinear inflation of the uncertainty envelope.
Validity of proposed bound over multiple prediction points: Our deterministic guarantees are derived for any query point $x_{N+1}$ given the training data -- this is also what is done in our experiments (cf. Figure 1 and 2), where the bound is computed for any input in the given domain.
Deployment in downstream tasks: We agree that quantifying the performance of the proposed approach in a downstream task would better highlight some of the benefits. Therefore, we conducted the following experiment. We are given a dynamical system with unknown (scalar) dynamics, which is estimated from $10^2$ data points using kernel-based estimation. The estimated dynamic model and its uncertainty bound is utilized in a control barrier function (CBF) to ensure that the optimized control input satisfies safety-critical constraints. Through the size of the uncertainty envelope, the conservatism of the employed bound determines the feasible domain for which safe inputs can be computed, with smaller confidence sets leading to a larger feasible region. Similar to the numerical comparison in Section 4.3, we compare three methods for computing the uncertainty bound: (a) the probabilistic bound by $Abbasi-Yadkori 2013$ ; (b) the derived robust regression; (c) the derived robust regression using only the closest $10^1$ data points for each test point. Apart from the number of data points used, the corresponding optimization problems only differ in the number of optimization variables (for the robust bounds, we additionally optimize over the free noise parameter $\sigma$), as well as the definition of the scaling factor $\beta$, cf. Eq. (7) for the robust case and

$Fiedler+2024, Eq. (7)$ for the probabilistic case. Crucially, benchmark (c) also highlights the fact that the derived bound remains valid for any subset of data, independent of how it is chosen (in contrast to stochastic error bounds). The following table compares the robust and probabilistic bounds in terms of the computational complexity, as well as their conservatism in terms of the portion of the domain for which the size of the bound retains feasibility of the optimization problem:

While optimization over the noise parameter in (b) leads to larger computation times for the same number of data points, it leads to a larger feasible domain due to its reduced conservatism even in the zero-mean i.i.d. setting. The subset-of-data approach (c) highlights that the method can be computationally more efficient, while also providing rigorous guarantees for safety-critical applicants in downstream tasks.
On solving (9): The main difficulty when optimizing (9) lies in the fact that the optimizer $\sigma$ can tend towards the boundaries of $(0, \infty)$, which can be addressed numerically by using adequately permissive box constraints, for example by enforcing $\sigma^2 \in [10^{-6}, 10^6]$. To avoid scaling issues, in our experiments we have found that optimizing over $\log(\sigma^2)$ improves numerical stability; alternatively, squashing functions can be employed to optimize in the large parameter space. Additionally, we have added a small constant "jitter" $\sigma_{\mathrm{jitter}}^2 = 10^{-6}$ to the kernel matrix, i.e., used $K^f + (\sigma^2 + \sigma_{\mathrm{jitter}}^2) I$ in the implementation, as is commonly done to avoid numerically indefinite kernel matrices. Using these techniques, we have not experienced numerical issues in the limiting case $\sigma \rightarrow 0$, which we observed in our experiments for some query input locations $x_{N+1}$ that are close to a training input location. This is expected since, as shown in Proposition 2, the optimized bound attains a finite value for $\sigma \rightarrow 0$.

References:
  $Srinivas+2012$ N. Srinivas, A. Krause, S. M. Kakade, and M. W. Seeger, "Information-Theoretic Regret Bounds for Gaussian Process Optimization in the Bandit Setting," IEEE Transactions on Information Theory, vol. 58, no. 5, pp. 3250--3265, 2012;
  $Chowdury+2017$ S. R. Chowdhury and A. Gopalan, "On Kernelized Multi-armed Bandits," in Proceedings of the 34th International Conference on Machine Learning, PMLR, 2017, pp. 844--853;
  $Fiedler+2021$ C. Fiedler, C. W. Scherer, and S. Trimpe, "Practical and Rigorous Uncertainty Bounds for Gaussian Process Regression," Proceedings of the AAAI Conference on Artificial Intelligence, vol. 35, no. 8, Art. no. 8, 2021;
  $Tokmak+2025$ A. Tokmak, K. G. Krishnan, T. B. Schön, and D. Baumann, "Safe exploration in reproducing kernel Hilbert spaces", 2025, arXiv: arXiv:2503.10352;
  $Csaji+2025$ B. C. Csáji and B. Horváth, "Derandomizing Simultaneous Confidence Regions for Band-Limited Functions by Improved Norm Bounds and Majority-Voting Schemes", 2025, arXiv: arXiv:2506.17764;
  $Calliess+2020$ : J.-P. Calliess, S. J. Roberts, C. E. Rasmussen, and J. Maciejowski, "Lazily Adapted Constant Kinky Inference for nonparametric regression and model-reference adaptive control," Automatica, vol. 122, p. 109216, 2020;
  $Fiedler+2024$ : C. Fiedler, J. Menn, L. Kreisköther, and S. Trimpe, "On Safety in Safe Bayesian Optimization", 2024, arXiv: arXiv:2403.12948;
  $Abbasi-Yadkori 2013$ : Y. Abbasi-Yadkori, "Online learning for linearly parametrized control problems," University of Alberta, Edmonton, Alberta, 2013.

I thank the authors for their answers. The sublinear dependence on the norms in (7) is a key point which the authors point out, which does mitigate the issue of choosing these bounds, to some extent.

Regarding the multiple prediction points - yes, of course, the paper is about recovery of the latent function through regression. I was thinking more in the context of uncertainty (e.g. GP) where the predictive posterior across multiple points is correlated. I see this is not relevant here.

The dynamic system numerical experiment would be a very nice addition to the paper - thanks for adding this.

Thanks for addressing my questions.

We thank the Reviewer for the insightful comments.
On the energy bounds and adversarial noises: Assuming to have known RKHS-norm bound for the unknown latent function is quite standard in the literature $Srinivas+2012,Chowdhury+2017,Fiedler+2021$ (see the paragraph "Practical strategies to estimate RKHS-norm bounds" below for further discussion and existing techniques to estimate the RKHS norm). Therefore, given adequate RKHS-norm bounds, our approach is capable of addressing the case of adversarial noise, possibly with non-zero mean: by construction, the proposed bound accounts for the worst-case adversarial noise, regardless of its distribution. Note that this is typically not the case for standard regression results, which leverage the assumption that noises are zero-mean i.i.d. and cannot handle adversarial noise.
Discussion on the performance as the data-set size increases: We would like to clarify that the derived bound is tight, i.e., given the posed assumptions, it provides the exact description of uncertainty envelope for the latent function given the data. Under the assumption that the noise is zero-mean i.i.d., one can usually guarantee to learn the unknown function up to any tolerance. With the proposed method, we can also generate monotonically improving estimates of the unknown function; however, due to the adversarial noise assumption, there is no guarantee that the bounds will asymptotically shrink to zero (cf. Figure 2 in the paper). Note that standard stochastic regression bounds working with i.i.d. and zero-mean noises cannot handle adversarial (or just slightly biased) noises, and their well-known behavior of returning vanishingly small uncertainty bounds would lead to an incorrect bound in such a set-up.
Testing of parameter sensitivity and validation under noise correlation: Modeling noise sequences as energy-bounded members of some deterministic function space allows us by construction to handle biased and correlated noise. An informed choice of the noise kernel can thereby lead to improved estimates; in particular, by using the Dirac kernel $k^w$, the derived bounds are valid for any noise sequence satisfying $\sum_{i=1}^N w_i^2< \Gamma_w^2$ (see Equation (16)), even if the noise is correlated.
Practical strategies to estimate RKHS-norm bounds: Many state-of-the-art uncertainty bounds for kernel-based methods rely on RKHS-norm bounds $Srinivas+2012,Chowdury+2017,Fiedler2021$ , which are of high relevance in downstream tasks, such as safe reinforcement learning  $Berkenkamp+2023$ and Bayesian optimization  $Sui+2018$ . This has motivated works on estimating the RKHS-norm bound. For instance, $Tokmak+2025$ and its preceding works elaborate an extrapolation-based technique to estimate such a quantity, and

$Csaji+2025$ derives bounds for Paley-Wiener kernels of band-limited functions using a uniformly-randomized Hoeffding's inequality or an empirical version of Bernstein's inequality. Additionally, one would also be capable of detecting the inadequacy of the norm bounds (i.e., them being too small) by checking feasibility of the constraints in the formulation given in Equation (C.1): this would allow to rule out certain configurations, as they cannot explain the data.
What conditions would yield convexity in the $\sigma$-optimization problem? This is a very interesting point. We would like to highlight that we never claim that the problem is convex and application of the derived bound to downstream tasks also does not require this. Instead, the key point regarding Problem (9) is that any variable $\sigma\in(0,\infty)$ provides a correct error bound. Hence, for decision making, we can simply treat $(x_{N+1},\sigma)$ as a decision variable, and Theorem 1 ensures that we always obtain correct error estimates and that the optimal choice of $\sigma$ also ensures that this bound is sharp. As both the ground truth function $f^{\mathrm{tr}}$ and the upper bound $\bar{f}^\sigma(x)$ are typically non-convex, convexity in $\sigma$ plays no crucial role. Moreover, convexity is not required to achieve global optimality using simple gradient descent algorithms. In this regard, we expect that Problem (9) satisfies the Polyak-Łojasiewicz (PL)-conditions, although we currently have no formal proof.
Sensitivity of the bounds to parameter misspecification, especially when the RKHS does not contain the true function: The proposed method can deal with the presented misspecification, because we do not treat the noise sequence as a stochastic object, but as a member of a deterministic function space. As long as the latter is rich enough, the noise is therefore capable of compensating for the misspecification of the latent-function RKHS. In particular, if data are generated by $y_i=h(x_i)$ with some function $h$ that is not in the RKHS, then we can decompose $h(x)=f(x)+w(x)$, where $f$ lies in the RKHS and $\Gamma_f$ bounds its RKHS norm, while $w$ is the residual. By picking $k^w$ as the Dirac kernel, Assumption 1 holds with $\sum_{i=1}^N |h(x_i)-f(x_i)|^2<\Gamma_w^2$ (see Equation (16)): for the noise to "compensate" the RKHS misspecification, its bound on the RKHS norm of the noise has to be larger than the point-wise squared distance between the true function and the function that lies in the RKHS.

References:
  $Srinivas+2012$ N. Srinivas, A. Krause, S. M. Kakade, and M. W. Seeger, "Information-Theoretic Regret Bounds for Gaussian Process Optimization in the Bandit Setting," IEEE Transactions on Information Theory, vol. 58, no. 5, pp. 3250--3265, 2012;
  $Chowdury+2017$ S. R. Chowdhury and A. Gopalan, "On Kernelized Multi-armed Bandits," in Proceedings of the 34th International Conference on Machine Learning, PMLR, 2017, pp. 844--853;
  $Fiedler+2021$ C. Fiedler, C. W. Scherer, and S. Trimpe, "Practical and Rigorous Uncertainty Bounds for Gaussian Process Regression," Proceedings of the AAAI Conference on Artificial Intelligence, vol. 35, no. 8, Art. no. 8, 2021;
  $Tokmak+2025$ A. Tokmak, K. G. Krishnan, T. B. Schön, and D. Baumann, "Safe exploration in reproducing kernel Hilbert spaces", 2025, arXiv: arXiv:2503.10352;
  $Csaji+2025$ B. C. Csáji and B. Horváth, "Derandomizing Simultaneous Confidence Regions for Band-Limited Functions by Improved Norm Bounds and Majority-Voting Schemes", 2025, arXiv: arXiv:2506.17764;   $Berkenkamp+2023$ : F. Berkenkamp, A. Krause, and A. P. Schoellig, "Bayesian optimization with safety constraints: safe and automatic parameter tuning in robotics," Mach Learn, vol. 112, no. 10, pp. 3713--3747,2023;
  $Sui+2018$ : Y. Sui, V. Zhuang, J. Burdick, and Y. Yue, "Stagewise Safe Bayesian Optimization with Gaussian Processes," in Proceedings of the 35th International Conference on Machine Learning, PMLR, 2018, pp. 4781--4789.

We thank the Reviewer for the insightful comments.
On the claims in line 91: We definitely agree that the problem stated in (2) becomes finite-dimensional after applying the Representer Theorem presented in Lemma A.1: nevertheless, the original problem is stated as searching for functions in infinite-dimensional spaces, and without exploiting the RKHS structure it might appear to be intractable. We plan to use a more precise wording to clarify this point. Also, we will make sure not to use "analytical solution" for the general outcome of Theorem 1: indeed, analytical (closed-form) solutions are available only for the two limit-cases $\sigma \to 0$ and $\sigma \to \infty$. The general result of Theorem 1, encompassing both limit cases $\sigma \to 0$ and $\sigma \to \infty$ as well as the scenario $0 < \sigma < \infty$, reduces (2) to a scalar, unconstrained optimization problem. This formulation is amenable for efficient iterative optimization and each iterate returns a valid (iteratively improving) bound for the unknown latent function.
On Theorem 1 being standard in constrained convex optimization: One of the core contributions of our work is to reduce the challenging problem of non-parametric function regression under non-i.i.d. noise to the analysis of a convex optimization problem. This solution is only obvious in hindsight, as we directly posed the problem in terms of an optimization problem in Equation (2). Existing results, especially

$Kanagawa+2018, Prop. 3.8$ , also tried to address the same problem; however, their approach does not utilize a connection with convex optimization and, as we discuss in Section 4.2, fails to determine tight bounds for the robust regression problem. We note that the exact characterization of the optimal solution to this convex optimization problem is more involved than it may seem, as it cannot be simply assumed that both constraints are active -- instead, a careful handling of four different cases, corresponding to the possible sets of active constraints, is required. These are detailed in the supplementary material (Appendix C, pages 17-26).
On the relevance of the case $\sigma \to 0$: As we discuss in Section 3.3, degenerate kernel matrices leading to the case $\sigma\to 0$ occur in two practical situations:

1. when the hypothesis space is determined by a finite set of $r<N$ basis functions, making it finite-dimensional and thus returning a rank-deficient kernel matrix. This scenario occurs in many approximations of kernel-based regression that aim at reducing its $\mathcal{O}(N^3)$-computational complexity: these methods rely on providing the eigen-decomposition of the kernel operator thanks to Mercer's Theorem, and on considering a finite number of eigenfunction as the basis (approximately) representing the hypothesis space $Zhu+1997$ ; in the stochastic setting, a paradigmatic approach leveraging this rationale in frequency domain is that of $Lázaro-Gredilla+2010$ .

2. when the test input $x_{N+1}$ belongs to the training data-set, which can happen then evaluating the bound over the whole domain of interest (cf. Figure 1 in the paper).

Thus, the case $\sigma \to 0$ is also of practical relevance, and our analysis sheds some new light on regression with degenerate kernels. However, it is a byproduct of the general analysis in Theorem 1, which is capable of capturing and generalizing known results; see also the case $\sigma \to \infty$ and the discussion in Section 4.1 of the paper.
On Corollary 1: Thank you for this suggestion. In the current exposition, Corollary 1 follows directly from the more general claim of Proposition 2, which provides the analytic solution for the case $\sigma\to 0$ occurring when the kernel matrix is rank-deficient. From this perspective, the scenario in which the test-input belongs to the training data-set is one of the particular situations leading to a drop in the rank of the kernel matrix. However, we agree that Corollary 1 could potentially be established analogously to Proposition 1 -- an insight that we are happy to convey in the paper.

References:
  $Zhu+1997$ H. Zhu, C. K. Williams, R. Rohwer, and M. Morciniec, "Gaussian regression and optimal finite dimensional linear models," Aston University, NCRG/97/011, 1997;
  $Lázaro-Gredilla+2010$ M. Lázaro-Gredilla, J. Quinonero-Candela, C. E. Rasmussen, and A. R. Figueiras-Vidal, "Sparse spectrum Gaussian process regression," The Journal of Machine Learning Research, vol. 11, pp. 1865--1881, 2010;
  $Kanagawa+2018$ M. Kanagawa, P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur, "Gaussian Processes and Kernel Methods: A Review on Connections and Equivalences," arXiv:1807.02582, 2018;
  $Fogel 1979$ E. Fogel, "System identification via membership set constraints with energy constrained noise," IEEE Transactions on Automatic Control, vol. 24, no. 5, pp. 752--758, 1979.