pESTO: Iterative p-greedy kernel regression with risk-driven Early STOpping

Introduction

Given a set of measurements sampled at multivariate scattered points, a kernel-based greedy scheme partitions the data into a training set and a validation set. The training set is used to construct an initial coarse interpolant, which is then progressively refined by iteratively adding the most informative samples selected from the validation set according to a prescribed rule. Such selection rules may rely on residuals or on deterministic error indicators, such as the power function, that exclusively depends on the kernel and on the node locations. In the interpolation setting, the latter, namely the $p$-greedy scheme, provides a reliable and well-established criterion for a priori node selection. However, when observations are affected by noise, relying solely on the power function may lead to suboptimal performances.

Kernel regression in the presence of noise has been extensively studied (see, e.g., [1 Chapter 15]). Likewise, greedy techniques for the corresponding problem have been considered [2]. However, to the best of our knowledge, a stopping rule for the respective greedy algorithms is still missing. Of course, the availability of a stopping criterion is of interest not only from a theoretical perspective, but also for computational efficiency. Typical stopping rules for kernel methods in the absence of noise can be based on the analysis of errors or residuals [3]. In the general case of iterative methods with noisy data, one can also consider the discrepancy principle introduced by Morozov [4], Generalized Cross Validation (GCV) [5], and risk-based criteria [6]. Among these methods, the Stein’s unbiased risk estimator (SURE) [7] has been employed as a data-driven criterion for parameter selection [8]. In the present work, we propose to employ such an estimator to identify the optimal iteration of greedy methods applied to kernel interpolation in the presence of noise. This is of interest, for instance, in combination with inverse problems for the purpose of image reconstruction [9].

The goal of the present work is to study error estimates for noisy data that can be used as the target quantity of iterative greedy methods. In particular, we show that a statistically equivalent derivation of the power-function error bound in Reproducing Kernel Hilbert Spaces (RKHSs) needs to account for the trade-off between approximation accuracy and noise amplification. These non-sharp bounds, involving the power function, a Lebesgue-type term, and a noise contribution, are not computable in practice. Indeed, they require the knowledge of the norm of exact function, while only noisy observations are available. This motivates the use of SURE in this context as a proxy for the Mean Squared Error (MSE). Building upon this characterization, the $p$ESTO algorithm, proposed in this paper, provides a principled stopping rule for the iterative greedy selection by balancing approximation error and statistical noise effects.

The problem we deal with reads as follows, and we refer the reader to the monographs [1] [10, 11] for further details on kernel-based methods. Let $\kappa: \Omega\times\Omega\to\mathbb{R}$ be a strictly positive definite kernel and $\mathcal{H}_{\kappa}$ the corresponding RKHS. Let $X_n = \{x_1,\ldots,x_n\} \subset\Omega\subseteq\mathbb{R}^{d}$ be a set of nodes and let $f \in\mathcal{H}_{\kappa}(\Omega)$ be the unknown function to approximate for which noisy observations are available. We assume that the associated observations are realizations of a Gaussian random variable $\mathbf{y}_n$ with independent components, namely $y_i := (\mathbf{y}_n)_i = f(x_i) + \varepsilon_i$, $i = 1,\ldots,n$, where $\boldsymbol{\varepsilon}_n = (\varepsilon_1,\ldots,\varepsilon_n) \sim\mathcal{N}(0,\operatorname{diag}(\sigma_1^2,\ldots,\sigma_n^2))$. Consequently, $\mathbf{y}_n = (y_1,\ldots,y_n) \sim\mathcal{N}(\mathbf{f}_n,\operatorname{diag}(\sigma_1^2,\ldots,\sigma_n^2))$, where $\mathbf{f}_n = (f_1,\ldots,f_n)^{\mathsf{T}} = (f(x_1),\ldots,f(x_n))^{\mathsf{T}}$. With abuse of notation, we denote both the random variables and their realizations by the same boldface symbol. We further assume that the noise variance is constant, i.e., $\sigma_i^2 = \sigma^2$ for $i = 1,\ldots,n$.

In RKHSs, the problem consists in computing the regularized interpolant

whose solution is

where $\mathbf{c}_n=(c_1,\ldots,c_n)^{\mathsf{T}}=(K_n+\lambda I_n)^{-1}\mathbf{y}_n$, $(K_n)_{ij}:=\kappa(x_i,x_j)$, $i,j=1,\ldots,n$, and $I_n$ is the $n\times n$ identity matrix. Under the assumptions listed above, $K_n$ is positive definite and the existence and uniqueness of the regularized interpolant $s_{X_n,y_n}$ is guaranteed by the generalized representer theorem [12] for any $\lambda\ge0$. In the case when $K_n$ is only positive semidefinite, the uniqueness of the interpolant coefficients is no longer guaranteed for $\lambda= 0$ [2 Theorem 1].

Greedy methods construct surrogate models for $f$ by selecting a subset of the data available at $X_n$. At the $(j+1)$-th step of the iteration, a new node is selected as $x_{j+1} = \operatorname*{argmax}_{x \in X_n \setminus X_j} \eta(x)$, where $X_j \subseteq X_n$ and $\eta$ is a given error indicator, such as the power function; we refer the reader to [3] [13, 14]. The objective of this work is to develop a new error bound that balances approximation accuracy and noise amplification. Furthermore, as this bound is not computable in practice, we utilize the SURE capturing both contributions. The resulting method, referred to as pESTO, provides a reliable estimate of the MSE and yields an a priori stopping criterion for the $p$-greedy scheme [9].

The paper is organized as follows. § 2 presents the expected approximation error. § 3 highlights the connection with SURE. § 4 reports numerical results supporting our findings, and conclusions are drawn in § 5.

Expected interpolation error

To begin with, we analyze the expected value of the error in the presence of noisy data. The following results hold for any subset $X_j$ of the set of interpolation nodes $X_n$. In the subsequent sections, we will in fact propose a greedy strategy for selecting such a subset. The results presented in this section provide theoretical error estimates which may be hard to compute in practice. For this reason, in the next section we will introduce a statistically-based and computable proxy.

Letting $\mathbf{k}_j(x)^\mathsf{T} = \left(\kappa(x,x_1),\ldots,\kappa(x,x_j)\right)$ and $\mathbf{c}_j = (K_j+\lambda I_j)^{-1}y_j$, we obtain

where $\boldsymbol{\ell}_{\lambda,j}(x)^\mathsf{T} = \mathbf{k}_j(x)^\mathsf{T}(K_j+\lambda I_j)^{-1}$ are known as the regularized Lagrange or cardinal functions. The following proof follows the classical arguments used in the noise-free setting (see, e.g., [1 §9]), while incorporating the additional contribution due to the observation noise.

is the sum of the squared regularized Lagrange functions defined in [2 Equation (22)] and referred to as “$\ell_2$-Lebesgue function”.

Proof. Since $s_{X_j,f_j}(x)=\boldsymbol{\ell}_{\lambda,j}(x)^\mathsf{T}\mathbf{f}_j$, by [1 Equation (9.6)] one has

Hence, by the reproducing property and the Cauchy-Schwarz inequality

In particular, recalling that the power function $P_{X_j,\kappa}(x)$ satisfies

we have $\left|f(x)-s_{X_j,f_j}(x)\right|\leq\|f\|_{\mathcal{H}_\kappa}P_{X_j,\kappa}^{2}(x)$. Moreover, by the Cauchy-Schwarz inequality and basic properties of normal distributions

Since $(a+b)^2\leq2(a^2+b^2)$ for any $a,b\in\mathbb{R}$, and by applying the triangular inequality, we get the bound in (2), i.e.,

The above obtained (4) gives the inequality in (2). The remaining equality in (2) immediately follows observing that $\left|f(x)-s_{X_j,y_j}(x)\right|^2$ does not depend on the elements in $\mathbf{y}_n\setminus\mathbf{y}_j$.

Remark 2.2. Let $\mathcal{L}_{X_j,\kappa,\lambda}(x):=\sum_{h=1}^{j}|(\ell_{\lambda,j}(x))_h|$ be the usual Lebesgue function. Then, since its terms are positive, one clearly has

That is, the square of the Lebesgue constant $\sup_x\mathcal{L}_{X_j,\kappa,\lambda}(x)$ is a (non sharp) upper bound for $\Lambda_{X_j,\kappa,\lambda}^{2}(x)$. The role of the Lebesgue constant in kernel approximation is of independent interest [15].

Remark 2.3. In [2], the term $\Lambda_{X_j,\kappa,\lambda}$ is related with the fill distance and the regularizing parameter $\lambda$. In particular, in such a paper it is proven that $\Lambda_{X_j,\kappa,\lambda}$ is uniformly bounded by $F_\kappa(h_{X_n})/\sqrt{\lambda}$, being $F_\kappa$ a positive function (depending on the kernel $\kappa$) and $h_{X_n}$ the fill distance associated with the set $X_n$. This makes the last term of (4) uniformly bounded for $\lambda$ large enough, and large for small $\lambda$.

The right-hand side of Equation (2) contains two terms. The first one accounts for the usual approximation error, while the second encodes information on the amplification of noise. The balance between the two terms is given by the signal-to-noise ratio $\lVert f\rVert_{\mathcal{H}_{\kappa}}/\sigma$. The main practical limitation is that the RKHS norm $\lVert f\rVert_{\mathcal{H}_{\kappa}}$ is unknown and cannot be directly estimated from the available noisy observations. Moreover, estimating this norm is already a challenging task even in the noise-free setting (e.g. [16, 17]). Therefore, greedy approaches based on the minimization of (2) cannot be applied in practice. To overcome this issue, we propose pESTO: we iteratively select the nodes by minimizing the power function (3) until a probabilistic stopping criterion that takes into account noise amplification is met.

In the next Section we propose an approach levering on statistical risk to estimate the expected interpolation error on the entire set $X_n$ of the evaluation nodes. To this end, we shall need an immediate consequence of Proposition 2.1.

Risk estimator

In the noisy setting, evaluations of the interpolant on the set of nodes, i.e., $(s_{X_j,\boldsymbol{y}_j}(x_1),\ldots,s_{X_j,\boldsymbol{y}_j}(x_n))^{\mathsf{T}}$, depend on the random variable $\boldsymbol{y}_j$ and, therefore, form a random vector. Realizations of this random vector provide estimates of the unknown exact evaluations $\boldsymbol{f}_n=(f_1,\ldots,f_n)^{\mathsf{T}}$. For this reason, such a vector can be thought of as an estimator $\hat{\boldsymbol{f}}(\boldsymbol{y}_j)$ of $\boldsymbol{f}_n$ [18, pag. 286].

The estimator $(s_{X_j,\boldsymbol{y}_j}(x_1),\ldots,s_{X_j,\boldsymbol{y}_j}(x_n))^{\mathsf{T}}$ depends on the number $j$ of interpolation nodes. In practice, it may be convenient to construct the interpolant over a set of nodes of cardinality $j<n$. This subset can be selected by a greedy procedure. To determine the optimal number $j$ of nodes we propose the use of the Stein’s unbiased risk estimator (SURE) to estimate the (expected) distance of $(s_{X_j,\boldsymbol{y}_j}(x_1),\ldots,s_{X_j,\boldsymbol{y}_j}(x_n))^{\mathsf{T}}$ from $\boldsymbol{f}_n$ and terminate the greedy procedure at the iteration $j$ minimizing such a distance. In contrast with the theoretical error bounds derived in § 2, SURE is explicitly computable since it depends only on the available noisy data.

Given an estimator $\hat{\boldsymbol{f}}(\boldsymbol{y}_j)$ of $\boldsymbol{f}_n$, SURE is defined as

where, throughout the paper, $\partial_i$ indicates the partial derivative with respect to the $i$-th component of $\boldsymbol{y}_j$. SURE provides an estimate of the expected MSE distance between the unknown exact evaluations $f_n$ and $(s_{X_j,y_j}(x_1),\ldots,s_{X_j,y_j}(x_n))^\mathsf{T}$. In the statistical setting, this expectation is called the risk of $(s_{X_j,y_j}(x_1),\ldots,s_{X_j,y_j}(x_n))^\mathsf{T}$ with respect to the MSE loss.

SURE is an unbiased estimator of the expected value of the MSE,

SURE can therefore be computed to estimate the squared distance between the values of the interpolant when evaluated over the set of nodes $X_n$ and the noise-free data $f_n$ (i.e., the left-hand side of Equation (5)), provided that the noise variance $\sigma^2$ is available. This is a standard assumption as in many practical applications the noise level is in fact either characterized by the acquisition process or estimated from the observations [19]. The following holds true.

Proof. From

it follows

Plugging this into (6) one finds (8). Finally, Equation (9) follows from (7). $\square$

being $\mu_i$ the $i$-th eigenvalue of $D_j$ (and hence of $K_j$). As, further, $K_j$ is positive definite and $\lambda> 0$ one has $\mu_i/(\mu_i + \lambda) \le1$ for each $i$, and this yields the chain of inequalities

Numerical verifications show that, if the magnitude of $\lambda$ is smaller than the majority of the eigenvalues of $K_j$, the first bound in (10) gives an acceptable and sensibly cheaper approximation for the central term of SURE, while the second appears to be generally too rough.

Remark 3.3. There is an analogy between $\mathcal{E}(f,\kappa,\sigma,\lambda,j)$ and SURE. Indeed, the term involving the squared power function in (5) and the average of the squared residuals $\frac{1}{n}\sum_{i=1}^{n}(s_{X_j,y_j}(x_i)-y_i)^2$ represent the approximation error and decrease for increasing number of interpolation nodes $j$. Differently, the term involving the sum of the squared regularized Lagrange functions in (5) and $2\sigma^2\operatorname{trace}(K_j(K_j+\lambda I_j)^{-1})$ describe the noise amplification and increase over the iterations of the greedy approach. This heuristic interpretation of the constituting terms of both the expected interpolation error and SURE suggests that the greedy approach may semi-converge, in the sense that the quality of the interpolant improves up to some iteration; however, due to noise propagation, the interpolant diverges from the ground truth function after more interpolation nodes are added. Consistently, both the expected interpolation error and SURE reach their respective minimum at a number of iterations smaller than the total number of iterations, and attain such a minimum value when approximation error balances noise amplification. This is a reflection of the bias–variance tradeoff defined in statistical settings, even if the minima of the two quantities need not coincide (as discussed in the numerical § 4). With reference to Remark 3.2, we note that replacing the trace term in SURE with the constant $2\sigma^2$ (see Equation (10)) removes the dependence on the size $j$ of the interpolation set and prevents the typical increase of the tail of the estimator.

An example of the tradeoff and the respective tail behavior is shown in Figure ??, which refers to the numerical experiments subsequently discussed in § 4.

By definition, SURE only depends on the random variable $\mathbf{y}_n$ and can hence be computed without knowledge of $f_n$, which is unknown. Therefore, in practice, we compute SURE by considering the given realization of the data in place of the random variable. As we pointed out at the end of § 2, greedy approaches based on the minimization of $\mathcal{E}(f,\kappa,\sigma,\lambda,j)$ are not viable. Therefore, the proposed a hybrid greedy approach operates as follows: (1) the nodes are selected independently from the data by minimizing the power function $P_{X_j,\kappa}^{2}$, and (2) a data-driven rule based on the minimization of the SURE estimator is employed to identify the greedy iteration that provides the best trade-off between approximation accuracy and noise amplification. The idea is to select the nodes by minimizing only one of the two terms of the theoretical error bound given in (5) and to control how the second term increases by monitoring SURE. The resulting $p$ESTO procedure is summarized in Algorithm 1.

Algorithm 1 Pseudo-code of the $p$ESTO algorithm

1. Take an initial interpolation set $X_1=\{x_1\}$.

2. For $j=1,\ldots,n-1$: (i) Define

(ii) Set

(iii) Compute the regularized interpolant $s_{X_{j+1},y_{j+1}}$.

(iv) Compute $\operatorname{SURE}(s_{X_{j+1},y_{j+1}},y_n)$.

3. Select

4. Return $X_{j^*}$.

Remark 3.4. In the numerical experiments of § 4, we applied the $p$ESTO procedure as presented in Algorithm 1. However, owing to the typical behavior of SURE depicted in Figure 1, in practice, one can stop the execution of the algorithm when the value of SURE plateaus or increases for a given number of iterations.

Numerical experiments

For the numerical experiments, we consider two different kernels that are widely adopted in Radial Basis Function (RBF) theory: the Gaussian kernel and the inverse multiquadric kernel, which are respectively defined as

The shape parameter $r$ is set equal to 10 for both kernels. The optimization of the shape parameter is beyond our scope; however, for details regarding this aspect, we address the reader to, e.g., [1] [21].

To assess the effectiveness of the proposed method, we consider two different ground truth functions: the Franke’s function [22] and the product function [23] defined as

where $x=(x_{1},x_{2})$. For each kernel and ground truth function, we consider both the Halton interpolation nodes [24] and an equispaced grid in the domain $[0,1]\times[0,1]$. We test the performances of the proposed $\mathit{p}$ESTO using sets of nodes with different cardinality, equal to $13^{2}$, $17^{2}$, and $21^{2}$. We perturb the data with white Gaussian noise with standard deviation $\sigma$ equal to 0.05 or 0.1, and we generate 30 different realizations for each level of noise. We perform the proposed $p$-greedy approach from each different data realization. In particular, at each step of the greedy method, we construct the regularized interpolant using (1), where we empirically set the regularization parameter $\lambda$ equal to $10^{-6}$. This value of $\lambda$ makes the bound of Remark 3.2 tight.

To quantitatively evaluate the effectiveness of the proposed $\mathit{p}$ESTO approach, we utilize two different metrics: $\operatorname{MSE}(\mathcal{X}_{N})$, which represents the MSE between the regularized interpolant and the ground truth function computed over a fine equispaced grid $\mathcal{X}_{N}$ of $N=200\times200$ points in $[0,1]\times[0,1]$, and the Structural Similarity Index (SSIM) [25]. Table 1 and Table 2 report the metric values corresponding to the greedy iteration selected by $\mathit{p}$ESTO. For comparison, we also report the metric values for the interpolant computer over all the nodes (corresponding to the last greedy iteration, i.e., for $j=n$), and for the optimal greedy iteration. The latter is defined as the greedy iteration minimizing $\operatorname{MSE}(\mathcal{X}_{N})$, which is, in essence, the target quantity. For each metric, we report the median and the Median Absolute Deviation (MAD) across the 30 noisy data realizations. Such indicators are robust to the presence of outliers.

Kernel	Nodes	Num.	$\sigma$	MSE($\mathcal{X}_N$) pESTO	MSE($\mathcal{X}_N$) All	MSE($\mathcal{X}_N$) Opt.	SSIM pESTO	SSIM All	SSIM Opt.
Gaussian	Halton	$13^2$	0.05	$3.0\pm0.8$	$5.2\pm0.7$	$2.6\pm0.2$	$0.88\pm0.04$	$0.76\pm0.03$	$0.89\pm0.02$
Gaussian	Halton	$13^2$	0.1	$7.9\pm1.8$	$19.8\pm2.6$	$7.5\pm1.6$	$0.76\pm0.06$	$0.49\pm0.03$	$0.77\pm0.04$
Gaussian	Halton	$17^2$	0.05	$2.5\pm0.6$	$10.0\pm2.4$	$2.3\pm0.4$	$0.90\pm0.03$	$0.61\pm0.04$	$0.90\pm0.02$
Gaussian	Halton	$17^2$	0.1	$7.8\pm1.3$	$40.2\pm10.0$	$7.5\pm1.2$	$0.76\pm0.03$	$0.32\pm0.04$	$0.77\pm0.03$
Gaussian	Halton	$21^2$	0.05	$2.0\pm0.3$	$27.7\pm6.0$	$2.0\pm0.2$	$0.91\pm0.02$	$0.43\pm0.03$	$0.91\pm0.02$
Gaussian	Halton	$21^2$	0.1	$6.8\pm0.9$	$110.7\pm24.1$	$6.8\pm0.8$	$0.78\pm0.04$	$0.19\pm0.02$	$0.78\pm0.04$
Gaussian	Equispaced	$13^2$	0.05	$2.2\pm0.2$	$2.1\pm0.3$	$2.1\pm0.3$	$0.90\pm0.02$	$0.88\pm0.01$	$0.89\pm0.02$
Gaussian	Equispaced	$13^2$	0.1	$8.4\pm1.2$	$8.3\pm1.1$	$7.7\pm1.1$	$0.74\pm0.04$	$0.67\pm0.02$	$0.75\pm0.06$
Gaussian	Equispaced	$17^2$	0.05	$2.2\pm0.4$	$2.3\pm0.2$	$2.1\pm0.2$	$0.91\pm0.02$	$0.82\pm0.02$	$0.90\pm0.03$
Gaussian	Equispaced	$17^2$	0.1	$8.0\pm1.0$	$9.1\pm1.0$	$7.7\pm0.8$	$0.76\pm0.03$	$0.56\pm0.03$	$0.76\pm0.03$
Gaussian	Equispaced	$21^2$	0.05	$2.3\pm0.3$	$2.4\pm0.2$	$2.2\pm0.2$	$0.88\pm0.05$	$0.76\pm0.02$	$0.90\pm0.02$
Gaussian	Equispaced	$21^2$	0.1	$7.3\pm1.2$	$9.8\pm1.0$	$7.2\pm1.1$	$0.77\pm0.05$	$0.47\pm0.03$	$0.77\pm0.05$
Inverse multiquadric	Halton	$13^2$	0.05	$1.7\pm0.4$	$2.8\pm0.3$	$1.7\pm0.3$	$0.94\pm0.02$	$0.84\pm0.02$	$0.94\pm0.02$
Inverse multiquadric	Halton	$13^2$	0.1	$5.8\pm1.5$	$11.0\pm1.2$	$5.3\pm0.9$	$0.86\pm0.05$	$0.60\pm0.03$	$0.86\pm0.04$
Inverse multiquadric	Halton	$17^2$	0.05	$1.7\pm0.3$	$3.3\pm0.4$	$1.5\pm0.3$	$0.95\pm0.02$	$0.76\pm0.01$	$0.94\pm0.02$
Inverse multiquadric	Halton	$17^2$	0.1	$5.8\pm1.7$	$13.3\pm1.5$	$5.6\pm1.4$	$0.85\pm0.04$	$0.49\pm0.02$	$0.85\pm0.05$
Inverse multiquadric	Halton	$21^2$	0.05	$1.6\pm0.3$	$4.3\pm0.3$	$1.6\pm0.3$	$0.94\pm0.02$	$0.66\pm0.01$	$0.94\pm0.02$
Inverse multiquadric	Halton	$21^2$	0.1	$5.6\pm0.9$	$17.4\pm1.3$	$5.4\pm0.9$	$0.85\pm0.04$	$0.36\pm0.02$	$0.84\pm0.05$
Inverse multiquadric	Equispaced	$13^2$	0.05	$1.7\pm0.4$	$1.8\pm0.2$	$1.6\pm0.3$	$0.93\pm0.02$	$0.90\pm0.01$	$0.94\pm0.03$
Inverse multiquadric	Equispaced	$13^2$	0.1	$6.0\pm1.4$	$7.4\pm1.0$	$5.3\pm1.0$	$0.84\pm0.05$	$0.70\pm0.02$	$0.84\pm0.05$
Inverse multiquadric	Equispaced	$17^2$	0.05	$1.6\pm0.2$	$2.0\pm0.2$	$1.5\pm0.3$	$0.94\pm0.02$	$0.84\pm0.02$	$0.94\pm0.02$
Inverse multiquadric	Equispaced	$17^2$	0.1	$4.7\pm1.0$	$8.0\pm0.7$	$4.6\pm1.1$	$0.89\pm0.04$	$0.60\pm0.02$	$0.89\pm0.05$
Inverse multiquadric	Equispaced	$21^2$	0.05	$1.7\pm0.3$	$2.1\pm0.1$	$1.6\pm0.3$	$0.93\pm0.02$	$0.79\pm0.02$	$0.94\pm0.02$
Inverse multiquadric	Equispaced	$21^2$	0.1	$5.7\pm1.5$	$8.4\pm0.6$	$5.5\pm1.3$	$0.85\pm0.06$	$0.52\pm0.03$	$0.84\pm0.06$

Kernel	Nodes	Num.	$\sigma$	MSE($\mathcal{X}_N$) pESTO	MSE($\mathcal{X}_N$) All	MSE($\mathcal{X}_N$) Opt.	SSIM pESTO	SSIM All	SSIM Opt.
Gaussian	Halton	$13^2$	0.05	$3.2\pm0.9$	$5.1\pm0.8$	$2.8\pm0.4$	$0.85\pm0.05$	$0.75\pm0.03$	$0.88\pm0.01$
Gaussian	Halton	$13^2$	0.1	$9.5\pm2.6$	$20.0\pm3.1$	$8.2\pm1.4$	$0.70\pm0.07$	$0.48\pm0.03$	$0.73\pm0.04$
Gaussian	Halton	$17^2$	0.05	$2.6\pm0.6$	$10.2\pm2.6$	$2.4\pm0.4$	$0.87\pm0.04$	$0.59\pm0.04$	$0.88\pm0.02$
Gaussian	Halton	$17^2$	0.1	$8.3\pm1.8$	$40.7\pm10.4$	$7.9\pm1.7$	$0.74\pm0.04$	$0.31\pm0.04$	$0.75\pm0.03$
Gaussian	Halton	$21^2$	0.05	$2.3\pm0.4$	$27.7\pm6.1$	$2.2\pm0.3$	$0.88\pm0.02$	$0.42\pm0.03$	$0.88\pm0.02$
Gaussian	Halton	$21^2$	0.1	$7.5\pm1.1$	$110.7\pm24.4$	$7.4\pm1.2$	$0.75\pm0.03$	$0.18\pm0.02$	$0.76\pm0.03$
Gaussian	Equispaced	$13^2$	0.05	$2.3\pm0.2$	$2.1\pm0.3$	$2.1\pm0.3$	$0.87\pm0.02$	$0.87\pm0.01$	$0.87\pm0.02$
Gaussian	Equispaced	$13^2$	0.1	$8.5\pm1.8$	$8.3\pm1.0$	$8.1\pm1.1$	$0.71\pm0.05$	$0.66\pm0.02$	$0.69\pm0.05$
Gaussian	Equispaced	$17^2$	0.05	$2.3\pm0.2$	$2.3\pm0.2$	$2.2\pm0.2$	$0.89\pm0.03$	$0.81\pm0.02$	$0.88\pm0.05$
Gaussian	Equispaced	$17^2$	0.1	$8.0\pm0.9$	$9.2\pm1.0$	$7.8\pm0.7$	$0.75\pm0.03$	$0.55\pm0.03$	$0.75\pm0.03$
Gaussian	Equispaced	$21^2$	0.05	$2.3\pm0.4$	$2.4\pm0.2$	$2.2\pm0.3$	$0.87\pm0.03$	$0.76\pm0.02$	$0.87\pm0.03$
Gaussian	Equispaced	$21^2$	0.1	$7.8\pm1.4$	$9.8\pm0.9$	$7.7\pm1.1$	$0.74\pm0.04$	$0.46\pm0.03$	$0.75\pm0.04$
Inverse multiquadric	Halton	$13^2$	0.05	$2.2\pm0.4$	$3.0\pm0.3$	$2.0\pm0.2$	$0.91\pm0.01$	$0.83\pm0.02$	$0.91\pm0.01$
Inverse multiquadric	Halton	$13^2$	0.1	$6.7\pm1.4$	$11.2\pm1.2$	$6.2\pm1.0$	$0.83\pm0.03$	$0.60\pm0.03$	$0.84\pm0.03$
Inverse multiquadric	Halton	$17^2$	0.05	$1.8\pm0.3$	$3.3\pm0.3$	$1.7\pm0.2$	$0.92\pm0.02$	$0.75\pm0.02$	$0.92\pm0.02$
Inverse multiquadric	Halton	$17^2$	0.1	$6.1\pm1.3$	$13.2\pm1.3$	$5.9\pm1.2$	$0.84\pm0.04$	$0.46\pm0.02$	$0.84\pm0.04$
Inverse multiquadric	Halton	$21^2$	0.05	$1.7\pm0.4$	$4.4\pm0.3$	$1.6\pm0.3$	$0.92\pm0.02$	$0.65\pm0.01$	$0.92\pm0.01$
Inverse multiquadric	Halton	$21^2$	0.1	$5.6\pm1.1$	$17.4\pm1.3$	$5.5\pm1.2$	$0.84\pm0.03$	$0.34\pm0.01$	$0.83\pm0.04$
Inverse multiquadric	Equispaced	$13^2$	0.05	$1.9\pm0.4$	$1.9\pm0.2$	$1.7\pm0.2$	$0.91\pm0.01$	$0.88\pm0.01$	$0.92\pm0.02$
Inverse multiquadric	Equispaced	$13^2$	0.1	$6.3\pm1.1$	$7.4\pm1.0$	$5.8\pm0.8$	$0.82\pm0.05$	$0.69\pm0.02$	$0.84\pm0.03$
Inverse multiquadric	Equispaced	$17^2$	0.05	$1.7\pm0.3$	$2.0\pm0.2$	$1.6\pm0.3$	$0.92\pm0.02$	$0.84\pm0.01$	$0.92\pm0.02$
Inverse multiquadric	Equispaced	$17^2$	0.1	$5.9\pm1.4$	$8.1\pm0.8$	$5.4\pm1.6$	$0.84\pm0.05$	$0.59\pm0.03$	$0.85\pm0.04$
Inverse multiquadric	Equispaced	$21^2$	0.05	$1.7\pm0.4$	$2.1\pm0.1$	$1.7\pm0.4$	$0.92\pm0.02$	$0.79\pm0.01$	$0.92\pm0.01$
Inverse multiquadric	Equispaced	$21^2$	0.1	$6.0\pm1.4$	$8.4\pm0.6$	$5.8\pm1.2$	$0.85\pm0.04$	$0.51\pm0.02$	$0.83\pm0.04$

Results shown in Table 1 and Table 2 demonstrate that the subsets of nodes selected by $\mathit{p}$ESTO lead to performances comparable to those of the optimal subset of nodes. In fact, the median metric values corresponding to $\mathit{p}$ESTO are systematically consistent within MAD with those corresponding to the optimal subset of nodes. Further, the $\mathit{p}$ESTO stopping criterion almost systematically yields improved performances compared to the use of all nodes. The only exceptions occur in the case of interpolation nodes located on an equispaced grid, where sometimes the median $\operatorname{MSE}(\mathcal{X}_{N})$ values for the interpolant obtained with all nodes are slightly smaller (although comparable within MAD). However, for most of these exceptions, the SSIM metric of $\mathit{p}$ESTO is larger than that corresponding to all nodes, showing the $\mathit{p}$ESTO leads to comparable performances using a smaller number of nodes.

Table 3 compares the greedy iterations selected by the proposed $\mathit{p}$ESTO with the optimal ones (i.e., those minimizing the $\operatorname{MSE}(\mathcal{X}_{N})$ metric). Specifically, we report the median and MAD of the selected iterations computed across the 30 noise realizations. Results show that the median of the greedy iterations selected by $\mathit{p}$ESTO is consistent within MAD with the median of the optimal ones. Therefore, $\mathit{p}$ESTO represents a good proxy for $\operatorname{MSE}(\mathcal{X}_{N})$, although it does not make use of evaluations of the ground truth function. Only in three cases corresponding to the Gaussian kernel and $13^{2}$ equispaced interpolation nodes, the optimal iteration appears to be close to the last one. In this cases, $\mathit{p}$ESTO selects a sensibly smaller number of nodes, however, the corresponding metrics shown in Table 1 and Table 2 demonstrate that the performances of $\mathit{p}$ESTO are comparable to the optimal ones with fewer interpolation nodes.

Kernel	Nodes	Num.	$\sigma$	Franke $p$ESTO	Franke Opt.	Product function $p$ESTO	Product function Opt.
Gaussian	Halton	$13^2$	0.05	$100.0\pm25.9$	$102.0\pm10.4$	$108.0\pm35.6$	$96.5\pm20.0$
			0.1	$64.5\pm9.6$	$69.0\pm7.4$	$84.5\pm21.5$	$80.0\pm9.6$
	Halton	$17^2$	0.05	$98.5\pm27.4$	$88.5\pm12.6$	$123.0\pm34.1$	$111.5\pm21.5$
			0.1	$74.5\pm9.6$	$76.0\pm10.4$	$84.5\pm13.3$	$79.0\pm8.2$
	Halton	$21^2$	0.05	$91.0\pm24.5$	$91.0\pm21.5$	$123.5\pm41.5$	$125.5\pm26.7$
			0.1	$73.5\pm3.0$	$74.5\pm5.2$	$81.5\pm15.6$	$81.5\pm15.6$
	Equispaced	$13^2$	0.05	$119.0\pm29.7$	$167.0\pm3.0$	$114.0\pm23.7$	$168.5\pm0.7$
			0.1	$99.0\pm11.9$	$97.5\pm8.9$	$99.5\pm10.4$	$167.0\pm3.0$
	Equispaced	$17^2$	0.05	$93.0\pm8.9$	$103.0\pm21.5$	$98.0\pm18.5$	$115.5\pm35.6$
			0.1	$86.5\pm8.2$	$87.0\pm9.6$	$80.0\pm2.2$	$81.5\pm5.2$
	Equispaced	$21^2$	0.05	$111.0\pm22.2$	$114.0\pm15.6$	$129.5\pm28.2$	$125.5\pm31.9$
			0.1	$74.0\pm5.2$	$71.0\pm5.9$	$72.5\pm8.9$	$83.0\pm20.0$
Inverse multiquadric	Halton	$13^2$	0.05	$49.5\pm19.3$	$56.0\pm18.5$	$63.0\pm19.3$	$75.0\pm23.0$
			0.1	$35.0\pm8.9$	$35.0\pm9.6$	$41.0\pm8.2$	$38.0\pm7.4$
	Halton	$17^2$	0.05	$53.5\pm14.8$	$72.0\pm29.7$	$64.0\pm18.5$	$78.0\pm28.9$
			0.1	$44.5\pm4.4$	$49.0\pm17.0$	$45.5\pm11.9$	$52.0\pm23.0$
	Halton	$21^2$	0.05	$77.0\pm24.5$	$77.0\pm26.7$	$69.5\pm24.5$	$69.0\pm24.5$
			0.1	$48.0\pm14.1$	$48.5\pm15.6$	$48.0\pm14.8$	$48.0\pm14.8$
	Equispaced	$13^2$	0.05	$54.0\pm8.2$	$80.5\pm40.8$	$74.5\pm17.8$	$97.5\pm60.8$
			0.1	$48.0\pm18.5$	$50.0\pm18.5$	$44.0\pm17.0$	$46.0\pm14.1$
	Equispaced	$17^2$	0.05	$72.0\pm22.2$	$79.0\pm13.3$	$79.0\pm22.2$	$79.0\pm25.2$
			0.1	$38.0\pm17.8$	$39.0\pm19.3$	$51.0\pm17.8$	$59.5\pm24.5$
	Equispaced	$21^2$	0.05	$80.0\pm31.1$	$83.0\pm27.4$	$92.0\pm35.6$	$76.0\pm28.2$
			0.1	$42.0\pm21.5$	$42.5\pm23.0$	$48.5\pm15.6$	$53.5\pm23.0$

To visually appreciate the behaviour of SURE and compare it with the theoretical error bound derived in § 2, in Figure 2 and Figure 3 we display the sample mean of SURE, $\operatorname{MSE}(s_{X_{j},y_{j}},f_{n})$, and $\operatorname{MSE}(\mathcal{X}_{N})$ computed across the 30 noisy realizations, for the Franke’s and the product function, respectively. The shaded area corresponds to a region within one standard deviation around the mean of SURE. In all cases, the sample mean of $\operatorname{MSE}(s_{X_{j},y_{j}},f_{n})$ falls within the shaded area, thus showing that SURE is an unbiased estimator of $\operatorname{MSE}(s_{X_{j},y_{j}},f_{n})$. Further, consistent with the results presented in Table 3, the minimum of the sample mean of MSE($\mathcal{X}_N$) is close to the SURE one. It is interesting to observe that the theoretical error bound attains a minimum value at a generally larger greedy iteration. This suggests that considering the bound on the expected interpolation error as a stopping rule leads to a less computationally efficient strategy even assuming that $\lVert f\rVert_{\mathcal{H}_\kappa}$ is known.

Conclusion

In this work, we analyzed error bounds for kernel approximation in the presence of noisy data. In such a setting, usual error bounds based on the power function must be corrected by a statistical term. The balance of the two contributions is unknown, since the norm of the exact function is not available; as a consequence, bounds that separate approximation accuracy and noise amplification turn out to be ineffective. We proposed pESTO, a risk-driven stopping rule for iterative $p$-greedy kernel regression in the presence of noisy data, which leverages SURE to estimate the MSE without access to the true solution. Numerical results confirm that $p$ESTO effectively identifies the optimal greedy iteration. Future work will focus on extending the proposed framework to non-Gaussian noise models, in particular to Poisson noise [26].

Acknowledgements

MT and EP acknowledge the support of the Fondazione Compagnia di San Paolo within the framework of the Artificial Intelligence Call for Proposals, AIxtreme project (ID RoI: 71708). LBB and EP are further supported by the GOSSIP project (“Greedy Optimal Sampling for Solar Inverse Problems”) – funded by the Ministero dell’Università e della Ricerca (MUR) - within the PRIN 2022 program, CUP: E53C24002330001. LBB is also funded by the INδAM-GNCS project CUP: E53C25002010001. PM is supported by Swiss PRODEX grant for STIX.