A Uniform Error Bound for Stochastic Kriging: Properties and Implications on Simulation Experimental Design

In SK, the system output obtained at a point \({\bf x}\) in the d-dimensional input space \(\mathcal {X}\subseteq \mathbb {R}^d\) on the jth simulation replication \(\mathcal {Y}_{j}({\bf x})\) is modeled aswhere \(f({\bf x})\) denotes the unknown true mean response that we intend to estimate at \({\bf x}\in \mathcal {X}\). For ease of exposition, we assume that \(f: \mathcal {X}\mapsto \mathbb {R}\) represents a sample of a stationary mean-zero GP [37], in line with a vast array of GP literature [13, 32, 40, 44]. The spatial covariance between any two points in the GP is given by \(\Sigma _{\sf M}({\bf x}, {\bf x}^\prime):= \xi ^2 K({\bf x}, {\bf x}^\prime)\), where \(\Sigma _{\sf M}: \mathbb {R}^d \times \mathbb {R}^d \rightarrow \mathbb {R}_+\) denotes the spatial covariance function, \(\xi ^2\) is the spatial variance of the GP, and \(K({\bf x}, {\bf x}^\prime)\) denotes the kernel function satisfying \(K({\bf x}, {\bf x}^\prime)\le 1\) for all \({\bf x}, {\bf x}^\prime \in \mathcal {X}\). The simulation noise terms incurred at \({\bf x}\) on different replications, \(\varepsilon _j({\bf x})\)’s, are assumed to be independent and identically distributed (i.i.d.) normal random variables with mean zero and input-dependent variance \({\sf V}({\bf x}) := {\rm Var}(\varepsilon _j({\bf x}))\). The normality of \(\varepsilon _j({\bf x})\) could be anticipated, since in a discrete-event simulation, the output \(\mathcal {Y}_{j}({\bf x})\) typically represents the average of a large number of more basic random variables obtained on the jth simulation replication.

Given a fixed simulation budget B to expend for approximating the mean response surface via SK, an experimental design for performing the simulation runs can be given as \(\lbrace ({\bf x}_i, n_i)_{i=1}^k: \sum _{i=1}^k n_i = B\rbrace\), where k denotes the number of distinct design points selected from the input space \(\mathcal {X}\), \({\bf x}_1, {\bf x}_2,\ldots , {\bf x}_k\) denote the k design points in the design-point set \(\mathcal {D}:=\lbrace {\bf x}_i, i=1,2,\ldots ,k\rbrace\), and \(n_i\) represents the number of replications to apply at \({\bf x}_i\), \(i=1,2,\ldots ,k\). Based on the simulation dataset \(\mathbb {D}_{k,B}=\lbrace {\bf x}_i,\lbrace \mathcal {Y}_{j}({\bf x}_i)\rbrace _{j=1}^{n_i}, i=1,2,\ldots ,k: \sum _{i=1}^k n_i = B\rbrace\) obtained, one can obtain the SK predictor of \(f({\bf x}_0)\) at any \({\bf x}_0 \in \mathcal {X}\) as follows:and its corresponding predictive variance is given bywhere \({\bf X}= ({\bf x}_1^\top , {\bf x}_2^\top , \ldots , {\bf x}_k^\top)^\top\) denotes the \(k\times d\) design matrix; \(\bar{\mathcal {Y}}=(\bar{\mathcal {Y}} ({\bf x}_1), \bar{\mathcal {Y}} ({\bf x}_2), \ldots , \bar{\mathcal {Y}} ({\bf x}_k))^\top\) denotes the \(k\times 1\) vector of the sample averages of simulation outputs, with \(\bar{\mathcal {Y}} ({\bf x}_i)=f({\bf x}_i)+\bar{\mbox{ $\varepsilon $}}({\bf x}_i)\) and \(\bar{\mbox{ $\varepsilon $}}({\bf x}_i)\) denoting the average random noise incurred at \({\bf x}_i\), \(i=1,2,\ldots ,k\). With a slight abuse of notation, we use \(\Sigma _{\sf M}({\bf X},{\bf X})\) to denote the \(k\times k\) covariance matrix that records the spatial covariances across the k design points and \(\Sigma _{\sf M}({\bf x}_0,{\bf X})\) to represent the \(k\times 1\) vector that contains the spatial covariances between the k design points and the prediction point \({\bf x}_0\); and the kernel matrix \(K({\bf X},{\bf X})\) and the kernel vector \(K({\bf x}_0,{\bf X})\) can be understood in the same manner. The \(k\times k\) diagonal matrix \(\Sigma _{\varepsilon }\) denotes the variance-covariance matrix of the \(k\times 1\) vector of average random noise terms \(\bar{\mbox{ $\varepsilon $}}=\left(\bar{\mbox{ $\varepsilon $}}({\bf x}_1),\bar{\mbox{ $\varepsilon $}}({\bf x}_2),\ldots , \bar{\mbox{ $\varepsilon $}}({\bf x}_k)\right)^\top\) given by \(\Sigma _{\varepsilon }=\mbox{diag} \left({\sf V}({\bf x}_1)/n_1, {\sf V}({\bf x}_2)/n_2, \ldots , {\sf V}({\bf x}_k)/n_k\right)\).

A pointwise confidence interval (CI) of \(f({\bf x}_0)\) at a given prediction point \({\bf x}_0\in \mathcal {X}\) with a prescribed confidence level \((1-\alpha)\) can be constructed as \(\mu _{k,B}({\bf x}_0)\pm z_{1-\alpha /2}\sigma _{k,B}({\bf x}_0)\), where \(\alpha \in (0,1)\) and \(z_{1-\alpha /2}\) denotes the \((1-\alpha /2)\)-quantile of the standard normal distribution [25]. One way to construct a confidence region for \(f(\cdot)\) at N prediction points (say, \({\bf x}_{0,i}\in \mathcal {X}\), \(i=1,2,\ldots ,N\)) with confidence level \((1-\alpha)\) is to apply the Bonferroni correction and obtain a simultaneous confidence bound comprising N pointwise confidence intervals [26]: \(\mu _{k,B}({\bf x}_{0,i})\pm z_{1-\alpha /(2N)}\sigma _{k,B}({\bf x}_{0,i})\), \(i=1,2,\ldots ,N\). Nevertheless, the literature on constructing a joint confidence bound for SK over an arbitrary set of prediction points is sparse.

In the remainder of this work, we employ the following notation consistently. Given a function \(g: \mathcal {X}\mapsto \mathbb {R}\), let \(\Vert g\Vert _\infty\) denote its supremum norm. For a vector \({\bf v}\), \(\left\Vert {\bf v}\right\Vert\) denotes the Euclidean norm of \({\bf v}\). Let \({\bf I}_k\) denote the \(k\times k\) identity matrix. For any symmetric matrix \({\bf Q}\), the \(\ell _2\)-operator norm or spectral norm of \({\bf Q}\) is denoted as \(\Vert {\bf Q}\Vert\), defined as \(\Vert {\bf Q}\Vert := \lambda _{\max }({\bf Q})\), where \(\lambda _{\max }({\bf Q})\) and \(\lambda _{\min }({\bf Q})\), respectively, denote the maximum and minimum eigenvalues of \({\bf Q}\). Define \({\sf V}_{k,\max } :=\max \limits _{1\le i \le k}{\sf V}({\bf x}_i)/n_i\) and \({\sf V}_{k,\min } :=\min \limits _{1\le i \le k}{\sf V}({\bf x}_i)/n_i\) as the maximum and minimum noise variance at the k design points. Let \(\lesssim\) and \(\gtrsim\) denote inequalities up to a constant multiple and write \(a \asymp b\) when both \(a \lesssim b\) and \(a \gtrsim b\) hold.

This subsection performs an analysis in the fixed design setting when the kernel K has a finite smoothness level. We start with the following result, which provides an upper bound on the error of the regularized kernel based approximation \(s_\lambda\) given in Equation (14).

We note that the interior cone condition essentially states that the input space \(\mathcal {X}\) is sufficiently regular with no sharp corners or cusps; and the d-dimensional unit cube \(\mathcal {X}= [0,1]^d\) provides an example [41]. By setting \(f:=K_{{\bf x}_0}\), \(s_\lambda := \widehat{K}_{{\bf x}_0}\), \(j=0\), and \(q=\infty\) in Proposition 1, we obtain the following result, which reveals the optimal convergence rate of the maximum discrepancy between \(K_{{\bf x}_0}\) and its approximation \(\widehat{K}_{{\bf x}_0}\).

In the finite smoothness case under the fixed design setting, Corollary 1 indicates that \(\sigma _{k,B}({\bf x}_0) \lesssim k^{-\frac{\nu -d/2}{2d}}\) for any \({\bf x}_0\in \mathcal {X}\). The condition on the smoothness parameter \(\lambda\) stipulated in Corollary 1 can be expressed as \({\sf V}_{k,\max }\lesssim k^{1-\frac{2\nu }{d}}\) following the definition \(\lambda ={\sf V}_{k,\max }/\xi ^2\). This condition indicates that: (i) the upper bound of \({\sf V}_{k,\max }\) must decrease with the number of design points k; and (ii) given a fixed number of design points k, the higher the smoothness level \(\nu\) (respectively, the higher the input-space dimensionality d) a given problem has, the faster (respectively, the slower) the decreasing rate required on \({\sf V}_{k,\max }\). The above sheds light on the following important aspects of a desirable simulation experimental design to achieve the optimal convergence rate of \(\Vert \sigma ^2_{k,B}\Vert _{\infty }\). First and foremost, the design-point set must be quasi-uniform. Moreover, the numbers of simulation replications \(n_i\)’s must grow with the number of design points k, and the smoother the kernel K, the faster \(n_i\)’s must grow.

Given that the kernel has a finite smoothness level, the next result extends Theorem 2 by elucidating the order of the uniform error bound as \(k, B \rightarrow \infty\) in the fixed design setting.

When the kernel K has infinite smoothness, the corresponding RKHS \(\mathcal {H}\) can be continuously embedded into every classical Sobolev space \(W^{\nu }_p(\mathcal {X})\) for all \(\nu \in \mathbb {N}\) and some fixed \(p\in [1,\infty)\) [36]. By restricting our attention to two popular infinitely smooth kernels, respectively, the Gaussian kernel (abbreviated as G): \(K({\bf x},{\bf x}^\prime) = \exp (-\Vert {\bf x}-{\bf x}^\prime \Vert ^2/\theta)\) with \(\theta \gt 0\) and the inverse multiquadrics kernel (abbreviated as M): \(K({\bf x},{\bf x}^\prime) = (c^2+\Vert {\bf x}-{\bf x}^\prime \Vert ^2)^{-\xi },\) where \(c\gt 0\) and \(\xi \gt d/2\) (page 76 of Reference [48]), we reveal the convergence rate of \(\sigma ^2_{k,B}(\cdot)\) under the supremum norm and the asymptotic order of the uniform error bound as \(k, B \rightarrow \infty\) in the fixed design setting. We first state a result that provides an upper bound on the error of the regularized kernel-based approximation \(s_{\lambda }\) given in Equation (14).

Setting \(f:=K_{{\bf x}_0}\), \(s_\lambda := \widehat{K}_{{\bf x}_0}\), and \(q=\infty\) in Proposition 3 yields the following result on the optimal convergence rate of the maximum discrepancy between \(K_{{\bf x}_0}\) and its approximation \(\widehat{K}_{{\bf x}_0}\) in the infinite smoothness case:

In the infinite smoothness case under the fixed design setting, Corollary 2 indicates that \(\sigma _{k,B}({\bf x}_0) \lesssim \exp (-Ck^{\frac{1}{d}}/2)\) at any \({\bf x}_0\in \mathcal {X}\). The condition on the regularization parameter \(\lambda ={\sf V}_{k,\max }/\xi ^2\) stipulated in Corollary 2 can be expressed as \({\sf V}_{k,\max }\lesssim \exp (-(2C{d}^{-1}(\log k)+2\tilde{C})k^{1/d})\) and \({\sf V}_{k,\max }\lesssim \exp (-2(C+\tilde{C})k^{1/d})\) for a Gaussian kernel and an inverse multiquadrics kernel, or equivalently, \({\sf V}_{k,\max }\lesssim \exp (-\bar{C}k^{1/d})\) with \(\bar{C}=2(C+\tilde{C})\) if the logarithmic factor term is ignored. It implies that: (i) the upper bound of \({\sf V}_{k,\max }\) must decrease exponentially with \(k^{1/d}\); and (ii) given a fixed number of design points k, the higher the input-space dimensionality d of a given problem, the slower the decreasing rate required on \({\sf V}_{k,\max }\). The above also provides some insight into a desirable simulation experimental design to achieve the optimal convergence rate of \(\Vert \sigma ^2_{k,B}\Vert _{\infty }\). First, the design-point set must be quasi-uniform. Moreover, the numbers of simulation replications \(n_i\)’s must grow with the number of design points k at a much faster rate in the infinitely smoothness case compared to in the finite smoothness case.

Given that the kernel has infinite smoothness, the next result extends Theorem 2 by elucidating the order of the uniform error bound as \(k, B \rightarrow \infty\) in the fixed design setting.

We close this section with some practical insights into suitable experimental designs for SK metamodeling. To achieve a desirable convergence rate of \(\sigma _{k,B}(\cdot)\) under the supremum norm given a sufficiently large budget B to expend, an experimental design suitable for the infinite smoothness case is expected to use fewer design points but with more replications allocated to each design point to fulfill its condition on \({\sf V}_{k,\max }\), hence, a “sparse and deep” design is recommended. However, an experimental design suitable for the finite smoothness case can adopt relatively more design points but with fewer replications allocated at each point to fulfill its condition on \({\sf V}_{k,\max }\); hence, a “dense and shallow” design is recommended. Such insights echo those reported in other contexts in the stochastic simulation metamodeling literature; see, e.g., References [46, 55].

This subsection provides some further insight into the impact of experimental designs and budget allocation schemes on SK metamodeling in light of the results obtained in Section 4 for the fixed design setting.

Given a fixed number of design points k, it is not difficult to see that \({\sf V}_{k,\max } = \max \limits _{1\le i \le k}{\sf V}({\bf x}_i)/n_i\) is minimized by setting the number of replications at design point \({\bf x}_i\) toUnder the budget allocation in Equation (16), we have \(\lambda = {\sf V}_{k,\max }/\xi ^2 = {(B\xi ^2)}^{-1} \sum _{i=1}^k {\sf V}({\bf x}_i)\). When k becomes large, it follows thatwhere the last step follows from \(|\mathcal {X}|^{-1}\int _{\mathcal {X}} {\sf V}({\bf u}) {\tt d} {\bf u}\lt \infty\).

However, if an equal budget allocation rule is applied, i.e., \(n_i = B/k\) for \(i=1,2,\ldots ,k\), then it is easy to see that in this case \(\lambda = (B\xi ^2)^{-1} k\max _{1\le i\le k} {\sf V}({\bf x}_i)\). When k becomes large, we havewhere in the last step we have used \(\sup _{\mathcal {X}}{\sf V}({\bf x})\lt \infty\). In brief, when k is sufficiently large, different budget allocation rules lead to the same order of \(\lambda\), which only depends on the total budget B and the number of design points k.

The discussion above sheds some light on the relationship between the number of design points k and the total budget B to achieve the optimal decreasing rate of \(\Vert \sigma _{k,B}\Vert _{\infty }\). In light of the results for the fixed design setting in Section 4, we have the following remarks regarding the finite and the infinite smoothness cases. When the underlying kernel function has a smoothness level \(\nu \lt \infty\), the discussion above and Corollary 1 indicate that \(\lambda \asymp k/B \lesssim k^{1-\frac{2\nu }{d}}\), or equivalently, \(k\lesssim B^{\tfrac{d}{2\nu }}\), when the budget B is sufficiently large. Therefore, in addition to adopting a quasi-uniform design-point set, to achieve the optimal convergence rate of \(\sigma _{k,B}(\cdot)\), one should use a large number of design points k if the input-space dimensionality d of a given problem is high (respectively, the smoothness level \(\nu\) is low) and vice versa. However, when the kernel function has an infinite smoothness level, the discussion above and Corollary 2 indicate that \(\lambda \asymp k/B \lesssim \exp (-(2Cd^{-1}(\log k)+2\tilde{C})k^{1/d})\) (for a Gaussian kernel), or equivalently, \(k\lesssim (\log B)^d\) if the logarithmic factor terms are ignored when the budget B is sufficiently large. Hence, besides adopting a quasi-uniform design-point set, to achieve an optimal convergence rate of \(\Vert \sigma _{k,B}\Vert _{\infty }\) in this case, one should use a relatively large number of design points k if the input space dimensionality d of a given problem is high. However, the number of design points required is expected to be much smaller than that in the finite smoothness case.

Last but not least, we note that the root mean squared error of SK is upper bounded by the predictive standard deviation under the supremum norm:Therefore, the aforementioned desirable simulation experimental designs for achieving the optimal convergence rate of \(\Vert \sigma _{k,B}\Vert _{\infty }\) also apply if the target is to achieve an adequate global response surface fitting with a guaranteed RMSE convergence rate.

In this subsection, we conduct a theoretical investigation of the impact of noise variance estimation on the performance of the uniform error bound. Recall that, in light of Theorem 1, we can obtain a high-probability uniform bound as given in Expression (7) for \(f(\cdot)\), assuming the knowledge of true noise variances \({\sf V}({\bf x}_i)\) for \(i=1,2,\ldots ,k\). Refer to this bound as the nominal uniform bound. Denote \(\widehat{{\sf V}}({\bf x}_i)\) as the sample variance calculated based on \(n_i\ge 2\) replications simulated at \({\bf x}_i\) for \(i=1,2,\ldots ,k\). Define \(\widehat{LB}({\bf x}):= \widehat{\mu }_{k,B}({\bf x})-\sqrt {\beta (\tau)}\widehat{\sigma }_{k,B}({\bf x})-\widehat{\gamma }_B(\tau)\) and \(\widehat{UB}({\bf x}):= \widehat{\mu }_{k,B}({\bf x})+\sqrt {\beta (\tau)}\widehat{\sigma }_{k,B}({\bf x})+\widehat{\gamma }_B(\tau)\), which are constructed using the sample variance \(\widehat{{\sf V}}({\bf x}_i)\) in place of \({\sf V}({\bf x}_i)\) for \(i=1,2,\ldots ,k\). We refer to \([\widehat{LB}({\bf x}), \widehat{UB}({\bf x})]\) as the empirical uniform bound for \(f({\bf x})\) at all \({\bf x}\in \mathcal {X}\) and study its performance as k and B tend to infinity.

Define \(n^*_k : = \min _{1 \le i\le k} n_i\) as the minimum number of replications allocated across the k design points, and recall \({\sf V}_{k,\max } = \max \limits _{1\le i \le k}{\sf V}({\bf x}_i)/n_i\). The following result provides a high-probability bound of the maximum discrepancy between the diagonal entries in the true noise variance-covariance matrix \(\Sigma _{\varepsilon }\) and those in its estimate \(\widehat{\Sigma }_\varepsilon :=\mbox{diag} (\widehat{{\sf V}}({\bf x}_1)/n_1, \widehat{{\sf V}}({\bf x}_2)/n_2, \ldots , \widehat{{\sf V}}({\bf x}_k)/n_k)\), which holds for all \(k\ge 1\).

The proof of Proposition 5 is provided in Appendix A.1. Based on Proposition 5, we can analyze the difference between the empirical uniform bound and the nominal uniform bound as k and B tend to infinity.

The proof of Theorem 3 is deferred to Appendix A.2. We see from Theorem 3 that the difference between the empirical uniform bound and the nominal uniform bound becomes negligible so long as the numbers of replications at all design points grow sufficiently fast with the expansion of the design-point set \(\mathcal {D}\). While the conditions stipulated on the minimum number of replications \(n_k^*\) may seem stringent, the numerical results in Section 6 indicate that the empirical and nominal uniform bounds perform similarly in practice, unless the numbers of replications at some design points become very small. We note that the conservativeness in the conditions arises from the techniques used to upper bound various terms involving the inverse of \({\mathbf {\Sigma }}:=\left(\Sigma _{\sf M}+\Sigma _{\varepsilon }\right)\) in \(| \widehat{LB}({\bf x}) - LB({\bf x})|\) and \(| \widehat{UB}({\bf x}) - UB({\bf x})|\) and to ensure that these terms vanish faster than the nominal uniform bound’s half-width for all \({\bf x}\in \mathcal {X}\). Other proof techniques may provide less stringent conditions. The proof of Theorem 3 also reveals the importance of stabilizing the SK metamodeling process as the design-point set expands and the budget allocated grows. To be more precise, the condition of matrix \({\mathbf {\Sigma }}\) can become increasingly worse as \(k,B\rightarrow \infty\), which is reflected by the magnitude of \(\Vert {\mathbf {\Sigma }}^{-1}\Vert = 1/\lambda _{\min }({\mathbf {\Sigma }})\) in the proof of Theorem 3. This also leads us to briefly discuss the condition number of \({\mathbf {\Sigma }}\), \(\mbox{cond}({\mathbf {\Sigma }}) = \lambda _{\max }({\mathbf {\Sigma }})/ \lambda _{\min }({\mathbf {\Sigma }})\). It is known that the maximum eigenvalue typically incurs no problem and behaves nicely compared to the minimum eigenvalue (Chapter 12 of Reference [48]). Since \(\lambda _{\min }({\mathbf {\Sigma }}) \ge \lambda _{\min }(\Sigma _{\sf M}) + \lambda _{\min }(\Sigma _{\varepsilon }) \ge \max \lbrace \lambda _{\min }(\Sigma _{\varepsilon }), \lambda _{\min }(\Sigma _{\sf M})\rbrace\), we know more about the stability and accuracy of SK’s prediction if we know more about the lower bounds for \(\lambda _{\min }(\Sigma _{\sf M})\) and \(\lambda _{\min }(\Sigma _{\varepsilon })\). Lower bounds for \(\lambda _{\min }(\Sigma _{\sf M})\) are known to depend on the separation distance \(q_{\mathcal {D}}\). Recall from Definition 2 that \(q_{\mathcal {D}} \asymp h_{\mathcal {D}, \mathcal {X}}= \mathcal {O}(k^{-1/d})\) if the design-point set \(\mathcal {D}\) is quasi-uniform. In this case, it is known that \(\lambda _{\min }(\Sigma _{\sf M}) \gtrsim q_{\mathcal {D}}^{(2\nu -d)} = \mathcal {O}(k^{1-2\nu /d})\) in the finite smoothness case [51], whereas in the infinite smoothness case, \(\lambda _{\min }(\Sigma _{\sf M}) \gtrsim \exp (-C_0 q_{\mathcal {D}}^{-2})q_{\mathcal {D}}^{-d} = \mathcal {O}(\exp (-C_0 \cdot k^{2/d})k)\) for some \(C_0\gt 0\) (Section 12.2 of Reference [48]). To maintain SK’s numerical stability as \(k, B \rightarrow \infty\), \(\lambda _{\min }(\Sigma _{\varepsilon })= {\sf V}_{k,\min }\) should not diminish too rapidly. Since \(\lambda _{\min }({\mathbf {\Sigma }}) \ge \max \lbrace \lambda _{\min }(\Sigma _{\varepsilon }), \lambda _{\min }(\Sigma _{\sf M})\rbrace = \lambda _{\min }(\Sigma _{\varepsilon }) = {\sf V}_{k,\min }\) and \(\Vert {\mathbf {\Sigma }}^{-1}\Vert \le 1/\lambda _{\min }(\Sigma _{\varepsilon }) = {\sf V}^{-1}_{k,\min }\), the condition \({\sf V}_{k,\max }/{\sf V}_{k,\min }\lt \infty\) and that on the order of \({\sf V}_{k,\max }\) stipulated in Theorem 3 ensure the numerical stability of SK while also fulfilling the regularization parameter role played by the noise variances (manifested by \({\sf V}_{k,\max }\)) to achieve the best convergence rate of SK’s predictive variance. We will examine the impact of noise variance estimation on the performance of the uniform bound through numerical experiments in the next section.

Consider simulating an M/M/1 queue with service rate 1 and arrival rate x with \(x \in \mathcal {X}= [0.3, 0.9]\). In this example, the simulation output recorded on each replication at a given design point is the average number of customers in the system from time 0 to T time units. It is well known from queueing theory that the steady-state mean number of customers in the queue is \(f(x)=x/(1-x)\), which is the mean function we intend to estimate. The noise variance function is \({\sf V}(x)\approx 2x(1+x)/(T(1-x)^4)\) for T large [50]. When simulating the M/M/1 queue at each design point, we initialize each replication in steady state and set the run length of each replication T to 1,000. Besides the general experimental setup given at the beginning of Section 6, the specific setup for this example is as follows: The value of the grid constant \(\tau\) is set to ensure that \(\sqrt {\beta _k(\tau)}\sigma _{k,B}(x)\) dominates \(\gamma _{k,B}(\tau)\) in Equation (7). We consider using a total budget \(B \in \lbrace 2,560,25,600\rbrace\) and varying the number of distinct design points \(k\in \lbrace 16,32,128,512\rbrace\) when \(B=2,560\) and \(k\in \lbrace 32,128,512,2,048\rbrace\) when \(B=25,600\). Given each choice of k, we choose k design points according to a uniform design or a grid design, where the former comprises a sample of k independent and uniformly distributed points and the latter k equispaced grid points in \(\mathcal {X}\). The three budget allocation schemes mentioned at the beginning of Section 6 are applied given each combination of B, k, and the design-point sampling scheme. A grid of \(N = 1,000\) equispaced points is selected from \(\mathcal {X}\) for prediction. Summary of results. Tables 1 and 2 show the \(\widehat{\mbox{SCP}}\)’s achieved by the empirical uniform bound (\(\mbox{CI}_{{\mbox{eu}}}\)), the nominal uniform bound (\(\mbox{CI}_{{\mbox{u}}}\)), and the Bonferroni correction-based simultaneous confidence intervals (\(\mbox{CI}_{{\mbox{b}}}\)) under different experimental settings when SK is implemented with the Gaussian and Matérn kernels. The following common observations are made across all experimental settings: First, given a fixed budget B, the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{u}}}\) are typically higher than the target level 0.95, while the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{b}}}\) are not close to 0.95, except when the number of design points k is relatively small under the grid design. Second, the simultaneous coverage performance of \(\mbox{CI}_{{\mbox{eu}}}\) remains satisfactory as the number of design points k increases until it reaches a point where only a few replications are allocated to some design points, regardless of the kernel function and the design-point sampling scheme adopted; see, e.g., the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{eu}}}\) obtained under the equal allocation and unequal allocation 1 when \(k=512\) and \(B = 2,560\). However, the simultaneous coverage performance of \(\mbox{CI}_{{\mbox{b}}}\) typically deteriorates as k increases given a fixed budget B, regardless of the kernel function adopted. This is evidenced by examining the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{b}}}\) obtained across all experimental settings shown in Tables 1 and 2. Third, the \(\widehat{\mbox{SCP}}\)’s of all bounds increase with the budget B given that the number of design points k is fixed. Last, applying unequal budget allocation schemes typically help improve the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{b}}}\) when the number of design points k is not too large relative to the budget B given. Such an impact on the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{u}}}\) is not as obvious, however.

We next examine the maximum half-widths \(H_{\max }\) of \(\mbox{CI}_{{\mbox{eu}}}\), \(\mbox{CI}_{{\mbox{u}}}\), and \(\mbox{CI}_{{\mbox{b}}}\) obtained under different experimental settings. Figure 1 summarizes the magnitudes of \(H_{\max }\) (in a logarithmic scale) of \(\mbox{CI}_{{\mbox{eu}}}\), \(\mbox{CI}_{{\mbox{u}}}\), and \(\mbox{CI}_{{\mbox{b}}}\) obtained under the grid design with a total budget of \(B=2,560\) when SK is implemented with the Gaussian and Matérn kernels. The corresponding results obtained with \(B=25,600\) and those obtained with \(B \in \lbrace 2,560, 25,600\rbrace\) under the uniform design are shown in Figures 6–8 in Appendix D.1. We have the following observations from Figures 1 and 6: First, the \(H_{\max }\)’s of \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{u}}}\) are larger than those of \(\mbox{CI}_{{\mbox{b}}}\) under each experimental setting considered, regardless of the kernel function adopted. Second, the magnitudes of \(H_{\max }\) of \(\mbox{CI}_{{\mbox{eu}}}\), \(\mbox{CI}_{{\mbox{u}}}\), and \(\mbox{CI}_{{\mbox{b}}}\) decrease with the number of design points k when a large total budget \(B=25,600\) is applied under all budget allocation schemes. However, this trend may not hold for a small total budget. Third, for \(\mbox{CI}_{{\mbox{eu}}}\), \(\mbox{CI}_{{\mbox{u}}}\), and \(\mbox{CI}_{{\mbox{b}}}\), when a given budget B is allocated to a fixed number of design points, the magnitudes of \(H_{\max }\) produced under the equal allocation scheme are the highest. It is worth noting that, regarding all bounds, the two unequal budget allocation schemes lead to shorter half-widths than the equal budget allocation scheme without compromising the simultaneous coverage performance, as observed in Table 1. The magnitudes of \(H_{\max }\) produced under the two unequal allocation schemes are similar, with unequal allocation 2 leading to higher \(H_{\max }\) values than unequal allocation 1. Last, regarding the impact of the design-point sampling scheme, the uniform design tends to lead to more outliers of the \(H_{\max }\) values compared to the grid design; see Figures 7 and 8 in Appendix D.1 for details. Moreover, under the uniform design, the unequal allocation schemes yield shorter half-widths of all bounds compared to the equal allocation scheme, but the difference between the two unequal allocation schemes is slight.

We summarize the impacts of noise variance estimation, the design-point sampling scheme, and the choice of kernel function observed in the M/M/1 queueing example. First, comparisons based on the \(\widehat{\mbox{SCP}}\)’s and \(H_{\max }\)’s of \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{u}}}\) indicate that the performance of the empirical uniform bound is very similar to that of the nominal one, except when the number of replications is too few at some design points, since using the variance estimates can incur additional variability in the SK model fitting. This corroborates the theoretical results provided in Section 5.2. Second, with other experimental settings held constant, the performance of all bounds under the grid design that is quasi-uniform typically outperforms that under the uniform design. This aligns with the theoretical results provided in Section 4. Last, with other experimental settings held constant, the performance of all bounds built with the Matérn kernel is comparable to that built with the Gaussian kernel for the M/M/1 queueing example. However, the magnitudes of \(H_{\max }\) of all bounds built with the Matérn kernel tend to be higher than those with the Gaussian kernel; this is because a faster convergence rate of the predictive variance is achieved with the Gaussian kernel than with the Matérn kernel, corroborating Corollaries 1 and 2, along with the comments following the results. Moreover, when the number of replications at some design points is too small, using the Matérn kernel can result in fitting issues that render the SK model unstable and adversely affect simultaneous coverage performance; see the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{b}}}\) obtained under unequal allocation 1 when \(B=2,560\) and \(k=512\) shown in Table 1 and the corresponding \(H_{\max }\)’s shown in Figure 1(d) in this subsection and Figure 7(d) in Appendix D.1.

Consider the following dimension-flexible example adapted from the one studied by Reference [23], where \(\mathcal {X}= [-1,1]^d\) and d denotes the input-space dimension. The simulation output at design point \({\bf x}=(x_1,x_2,\ldots ,x_d)^\top\) on the jth replication is generated as \(\mathcal {Y}_{j}({\bf x}) = f({\bf x}) + \varepsilon _j({\bf x}),\) where \(\varepsilon _j({\bf x})\)’s are i.i.d. normally distributed with mean zero and variance \({\sf V}({\bf x})\). The noise variance function is \({\sf V}({\bf x}) = (2+\cos (\pi +\sum _{i=1}^d x_i/d))^2\) for any \({\bf x}\in \mathcal {X}\). The mean function of interest is given by \(f(x)=\sin (9x^2)\) if \(d=1\) and \(f({\bf x})=\sin (9x_1^2) + \sin ((\sum _{i=2}^d 3 x_i/(d-1))^2)\) if \(d\gt 1\). We describe the specific setup for this example next and refer the reader to the beginning of Section 6 for the general experimental setup. We consider two problem cases corresponding to \(d=1\) and \(d=5\). For \(d=1\), we vary the number of design points \(k \in \lbrace 16, 64, 256, 512\rbrace\) when \(B=2,560\) and \(k\in \lbrace 16,64,256,512,2,048\rbrace\) when \(B=10,240\). For \(d=5\), we vary \(k \in \lbrace 16, 64, 256, 512\rbrace\) when \(B=2,560\) and \(k \in \lbrace 64,256,512,1,024,2,048\rbrace\) when \(B=10,240\). Given each choice of k, a set of k design points is chosen according to a uniform design or a maximin Latin hypercube design (LHD). The three budget allocation schemes mentioned at the beginning of Section 6 are applied given each combination of B, k, and the design-point sampling scheme. On each macro-replication, we independently generate a set of \(N = 2,500\) points from \(\mathcal {X}\) according to a maximin LHD as the prediction points for performance evaluation.

Summary of results. We first look into the 1D case. Tables 3 and 4 present the \(\widehat{\mbox{SCP}}\)’s of the empirical uniform bound (\(\mbox{CI}_{{\mbox{eu}}}\)), the nominal uniform bound (\(\mbox{CI}_{{\mbox{u}}}\)), and the Bonferroni correction-based simultaneous confidence intervals (\(\mbox{CI}_{{\mbox{b}}}\)) under different experimental settings when SK is implemented with the Gaussian and Matérn kernels. Observations akin to those noted in the M/M/1 queueing example are summarized as follows: First, given a fixed budget B, the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{u}}}\) are typically higher than the target level 0.95. Moreover, the performance of \(\mbox{CI}_{{\mbox{eu}}}\) is very similar to that of the nominal one in this 1D case, except when the number of replications is too few at some design points. Second, the \(\widehat{\mbox{SCP}}\)’s of all bounds increase with the budget B given that the number of design points k is fixed. Third, applying unequal budget allocation schemes typically helps improve the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{b}}}\) when the number of design points k is not too large relative to the budget B given. Last, the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{b}}}\) deteriorate as the number of design points k becomes large, regardless of the budget allocation schemes adopted. We highlight that this 1D case differs from the M/M/1 example in the following aspects: First, the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{b}}}\) are not significantly lower than those of \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{u}}}\), so long as the number of design points k is not too large relative to the budget B given. Second, all bounds built with the Matérn kernel achieve comparable or higher \(\widehat{\mbox{SCP}}\)’s compared to those built with the Gaussian kernel, especially when the number of design points k is small.

We further delve into the maximum half-widths \(H_{\max }\) of \(\mbox{CI}_{{\mbox{eu}}}\), \(\mbox{CI}_{{\mbox{u}}}\), and \(\mbox{CI}_{{\mbox{b}}}\) obtained under different experimental settings. To save space, we concentrate on analyzing the results obtained under LHD with the equal allocation scheme. Similar observations apply to the results obtained under LHD with the two unequal budget allocation schemes and under the uniform design with all three budget allocations; see Figures 9 and 10 in Appendix D.2.1 for details. Figures 2 and 3 summarize the magnitudes of \(H_{\max }\) (in a logarithmic scale) of \(\mbox{CI}_{{\mbox{eu}}}\), \(\mbox{CI}_{{\mbox{u}}}\), and \(\mbox{CI}_{{\mbox{b}}}\) obtained under LHD with a total budget of \(B=2,560\) and 10,240. First, the \(H_{\max }\)’s of \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{u}}}\) are larger than those of \(\mbox{CI}_{{\mbox{b}}}\) under each experimental setting considered, regardless of the kernel function adopted. Second, the magnitudes of \(H_{\max }\) of \(\mbox{CI}_{{\mbox{eu}}}\), \(\mbox{CI}_{{\mbox{u}}}\), and \(\mbox{CI}_{{\mbox{b}}}\) decrease with the number of design points k when a large total budget \(B=10,240\) is applied. However, this trend is less discernible with a small budget. Third, the magnitudes of \(H_{\max }\) of \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{u}}}\) are similar, except when k becomes large relative to the given budget B so using variance estimates based on a few replications can incur additional variability in SK model fitting, which impacts \(\mbox{CI}_{{\mbox{eu}}}\); this observation holds regardless of the kernel function adopted. Fourth, applying unequal budget allocation schemes only slightly reduces the half-widths of all three bounds, as observed in Figures 9(a)–(d) and 10(a)–(d) in Appendix D.2.1. Last, the magnitudes of \(H_{\max }\) of all bounds built with the Matérn kernel tend to be higher than those with the Gaussian kernel, which explains the corresponding higher \(\widehat{\mbox{SCP}}\)’s observed in Table 3. We also highlight some differences in the 1D case from the M/M/1 example. Given a fixed total budget B, especially a small one, using a small number of design points k can lead to unstable SK model fitting when the Gaussian kernel is adopted, resulting in highly variable \(H_{\max }\) values for all bounds; see, e.g., the \(H_{\max }\)’s obtained for \(k=16\) in Figure 2(a). This is because the true mean function \(f(\cdot)\) in this 1D case is nonlinear, and its trend changes more frequently than that of the M/M/1 example across the input space, due to its combination of the trigonometric function and higher-order polynomial terms. It is known that the Gaussian kernel is not flexible enough for modeling a function with such features, and using a small number of design points exacerbates this deficiency in capturing the true mean function \(f(\cdot)\) in the 1D case.

We next investigate the 5D case. Compared to the 1D case, notice that the true mean function \(f(\cdot)\) in this 5D case is more complex and changes more rapidly across its input space. Tables 5 and 6 present the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{eu}}}\), \(\mbox{CI}_{{\mbox{u}}}\), and \(\mbox{CI}_{{\mbox{b}}}\) obtained under various experimental settings. We observe that, similar to the 1D case, the performance of \(\mbox{CI}_{{\mbox{b}}}\) built with the Gaussian kernel is notably weak, regardless of the budget size B. In contrast, \(\mbox{CI}_{{\mbox{b}}}\) built with the Matérn kernel delivers relatively robust performance under the same experimental settings. Furthermore, in this 5D case, an increase in the number of design points k is associated with higher \(\widehat{\mbox{SCP}}\)’s for \(\mbox{CI}_{{\mbox{b}}}\) especially when the Matérn kernel is adopted, underscoring the need of using a dense design with more design points for capturing the function in high-dimensional problems and enhancing the reliability of the confidence bound built. Regarding \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{u}}}\), the Matérn kernel also facilitates consistently strong simultaneous coverage performance across different experimental settings, highlighting its suitability for modeling highly nonlinear and rapidly changing response surfaces in high-dimensional problems. In contrast, due to the limitations of the Gaussian kernel in modeling the mean surface \(f(\cdot)\) in this 5D case, the \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{eu}}}\) and \(\mbox{CI}_{{\mbox{u}}}\) suffer especially when the number of design points k is small.

Figures 4 and 5 summarize the magnitudes of \(H_{\max }\) (in a logarithmic scale) of \(\mbox{CI}_{{\mbox{eu}}}\), \(\mbox{CI}_{{\mbox{u}}}\), and \(\mbox{CI}_{{\mbox{b}}}\), employing LHD with a total budget \(B=2,560\) and 10,240. For simplicity, the discussion below focuses on the results obtained under LHD with the equal allocation scheme. Similar observations apply to the results obtained under LHD with the two unequal budget allocation schemes and under the uniform design with all three budget allocations, which are detailed in Figures 11 and 12 provided in Appendix D.2.2. First and foremost, the variation exhibited by the \(H_{\max }\) values of all bounds built with the two kernels demonstrates the notable difference between the two kernels in this 5D case as the number of design points k increases, especially when the given budget B is small. Specifically, the Gaussian kernel leads to increased instability in modeling as the number of design points k grows, contrasting with the Matérn kernel, which maintains consistent modeling stability across the range of k values considered; compare Figures 4(a) and (b). Furthermore, given a large budget B, the \(H_{\max }\)’s of all bounds built using the Gaussian kernel show a decreasing trend as the number of design points k increases. However, such a trend seems absent when the Matérn kernel is adopted. This observation explains the corresponding low \(\widehat{\mbox{SCP}}\)’s of \(\mbox{CI}_{{\mbox{b}}}\) observed in Tables 5 and 6. An explanation is that the predictive variance (hence the confidence bound width) associated with the Gaussian kernel diminishes faster compared with the Matérn kernel, as indicated by Corollaries 1 and 2 and the comments immediately following them.

Finally, we summarize key insights obtained from the numerical evaluations. When applying SK for metamodeling purposes, selecting an appropriate kernel suitable for modeling the function of interest is crucial. The Gaussian kernel performs well in queueing applications like the M/M/1 example considered here [53]. In contrast, only the Matérn kernel with smoothness parameter \(\nu = 5/2+d/2\) was considered by Reference [23] from where the dimension-flexible example originated, which also showed its competence in our evaluations. Therefore, choosing an appropriate kernel function is example-dependent and should take into account the features of the function of interest. The performance of Bonferroni correction-based simultaneous confidence intervals is sensitive to this choice, and they can be unreliable even when built by adopting an appropriate experimental design. In contrast, the uniform error bound is more robust, especially in high-dimensional problems. Second, fixed designs can outperform random designs in facilitating stability in SK model fitting and simultaneous coverage performance. Last, unequal allocation schemes can enhance the performance of confidence bounds, but their impact is limited.

We are now in a position to state the assumption on the eigenfunctions of the kernel function K and the main result that will be useful for establishing the convergence rate of the predictive variance under the supremum norm in the random design setting.

Define \(\tilde{\kappa }^2 := \sup _{{\bf x}\in \mathcal {X}}\tilde{K}({\bf x},{\bf x})\) andwhich is often referred to as the effective dimensionality of kernel K with respect to \(L_2(\mathcal {X})\) [56]. It follows directly from Assumption 1 that \(\tilde{\kappa }^2 \le C_\psi ^2 \sum _{i=1}^\infty \nu _i \lesssim \gamma (\lambda)\).

Theorem 4 follows from applying Lemma 4 and Theorems 1 in Reference [30] to the noiseless KRR estimator \(\widehat{K}_{{\bf x}_0}\) given in Equation (A.33). For the sake of brevity, we omit the proof here. Notice that setting \(\delta = k^{-10}\) simplifies the bound in Expression (A.36) and facilitates further analysis. Such a choice is quite common in the nonparametric regression literature; see, e.g., References [31, 33, 54]. Theorem 4 enables us to further investigate the convergence rate of \(\Vert K_{{\bf x}_0} - \widehat{K}_{{\bf x}_0}\Vert _{\infty }\) for commonly used kernel functions.

Polynomially decaying kernels. Kernels with polynomially decaying eigenvalues include those that underlie the Sobolev spaces with different orders of smoothness. One popular kernel function belonging to this category is the Matérn kernel when \(\mathbb {P}\) is the uniform distribution on \(\mathcal {X}\) [40].

With respect to Corollary 3, one viable choice of \(\lambda\) is \(\lambda \asymp k^{-p}{(\log k)}^{b}\), where \(p\gt 0\) and \(b\in \mathbb {R}\). We can show that, if \(p=\nu /d\) and \(b\gt \nu /d\), then \(\tilde{\kappa }^2 \lesssim \gamma (\lambda) \lesssim {(k^{-p} {(\log k)}^{b})}^{-\tfrac{d}{2\nu }} \lesssim k^{\tfrac{1}{2}}{(\log k)}^{-\tfrac{d}{2\nu }b} = o(\sqrt {k/\log k})\). In this case, it follows that \(\Vert K_{{\bf x}_0} - \widehat{K}_{{\bf x}_0}\Vert _{\infty } \lesssim k^{-\tfrac{2\nu -d}{2d}}{(\log k)}^{\tfrac{2\nu -d}{2\nu }b}\). Hence, with a high probability, for any \({\bf x}_0 \in \mathcal {X}\), \(\sigma _{k,B}({\bf x}_0) \lesssim k^{-\tfrac{2\nu -d}{4d}}{(\log k)}^{\tfrac{2\nu -d}{4\nu }b}\). Recall the definition \(\lambda = {\sf V}_{k,\max }/(\xi ^2 k)\); it is not difficult to allocate the simulation budget to fulfill this choice of \(\lambda\). In case the budget allocation adopted makes \(\lambda\) decay to zero at a faster rate than the aforementioned case, Corollary 3 makes it clear that the corresponding upper bound on \(\Vert K_{{\bf x}_0} - \widehat{K}_{{\bf x}_0}\Vert _{\infty }\) and hence the upper bound on the order of \(\sigma _{k,B}(\cdot)\) still hold; however, it can be loose and does not imply the optimal convergence rate of \(\sigma _{k,B}(\cdot)\).

We can further extend Theorem 2 and study the order of the uniform error bound in the random design setting as \(k, B \rightarrow \infty\) when polynomially decaying kernels are adopted. We see from Expression (8) that both \(\gamma _{k,B}(\tau)\) and \(\sqrt {\beta _k(\tau)}\sup _{{\bf x}\in \mathcal {X}}\sigma _{k,B}({\bf x})\) must converge to zero as \(k,B \rightarrow \infty\) to guarantee a vanishing approximation error at \(\forall {\bf x}\in \mathcal {X}\). Given the aforementioned choice of \(\lambda \asymp k^{-p} {(\log k)}^{b}\) with \(p=\nu /d\) and \(b\gt \nu /d\), we have that \(\lim \limits _{k, B\rightarrow \infty } \gamma _{k,B}(\tau)=0\) can be achieved if the grid constant as a function of k, \(\tau (k)\), decreases sufficiently fast. Recall from Expression (12) that \(\gamma _{k,B}(\tau) \lesssim ((\beta _k(\tau) \cdot \tau (k) \cdot k {\sf V}_{k,\min }^{-1})^{\tfrac{1}{2}} + \tau (k) + \tau (k) \cdot k {\sf V}_{k,\min }^{-\tfrac{1}{2}}) .\) Assuming that \({\sf V}_{k,\max }/{\sf V}_{k, \min }\lt \infty\), we have \({\sf V}_{k,\min } \asymp {\sf V}_{k,\max } \asymp k^{1-p} (\log k)^{b}\). Setting \(\tau (k) = \mathcal {O}(k^{-a})\) as \(k, B\rightarrow \infty\) with \(a = 2\nu /d\) ensures that \(\lim \limits _{k, B\rightarrow \infty } \gamma _{k,B}(\tau)=0\) and moreover \(\gamma _{k,B}(\tau) = o(\sqrt {\beta _k(\tau)}\sup _{{\bf x}\in \mathcal {X}}\sigma _{k,B}({\bf x}))\) as \(k, B\rightarrow \infty\). The latter follows from the fact that \(\beta _k(\tau) = \mathcal {O} (\log k)\) due to Expression (9) and \(\sigma _{k,B}({\bf x}) = \mathcal {O} (k^{-\tfrac{2\nu -d}{4d}}{(\log k)}^{\tfrac{2\nu -d}{4\nu }b})\) at \(\forall {\bf x}\in \mathcal {X}\), thanks to Corollary 3. Therefore, \(\gamma _{k,B}(\tau) + \sqrt {\beta _k(\tau)}\sup _{{\bf x}\in \mathcal {X}}\sigma _{k,B}({\bf x}) \lesssim \sqrt {\beta _{k}(\tau)} \sup _{{\bf x}\in \mathcal {X}}\sigma _{k,B}({\bf x}) = \mathcal {O}(k^{-\tfrac{2\nu -d}{4d}}{(\log k)}^{\tfrac{2\nu -d}{4\nu }b+\tfrac{1}{2}})\) as \(k, B \rightarrow \infty\).

Exponentially decaying kernels. Consider kernels having exponentially decaying eigenvalues. Roughly speaking, infinitely smooth stationary kernels belong to this category. One important kernel function that falls into this category is the squared exponential or Gaussian kernel when \(\mathbb {P}\) is the uniform distribution or normal on \(\mathcal {X}\) [39, 40]: \(K({\bf x},{\bf x}^\prime) = \exp \left(-\Vert {\bf x}-{\bf x}^\prime \Vert ^2/\theta \right),\) where \(\theta \gt 0\) denotes the lengthscale parameter.

Some comments follow from Corollary 4. One viable choice of \(\lambda\) is \(\lambda \asymp \exp (-k^{\frac{1}{p}})\) for some \(p\gt 0\). By choosing \(p\gt 2d\), it follows that \(\tilde{\kappa }^2 \lesssim \gamma (\lambda) \lesssim k^{\frac{d}{p}} = o(\sqrt {k/\log k})\). The high probability upper bound given by Expression (A.40) in this case follows as \(\Vert K_{{\bf x}_0} - \widehat{K}_{{\bf x}_0}\Vert _{\infty } \lesssim \exp (-k^{\frac{1}{p}}) k^\frac{d}{p}\). Therefore, it follows that, with a high probability, for any \({\bf x}_0 \in \mathcal {X}\), \(\sigma _{k,B}({\bf x}_0) \lesssim \exp (-\frac{1}{2}k^{\frac{1}{p}}) k^\frac{d}{2p}\). Recall that \(\lambda := {\sf V}_{k,\max }/(\xi ^2 k)\); we see that it is possible to allocate the computational budget to fulfill this choice of \(\lambda\). In case the budget allocation adopted makes \(\lambda\) diminish at a faster rate than the case discussed above, the corresponding upper bound on \(\Vert K_{{\bf x}_0} - \widehat{K}_{{\bf x}_0}\Vert _{\infty }\) and hence the upper bound on the order of \(\sigma _{k,B}(\cdot)\) still hold; however, in this case, the upper bound is not as tight.

Finally, we further extend Theorem 2 by discussing the order of the uniform error bound in the random design setting as \(k, B \rightarrow \infty\) when exponentially decaying kernels are adopted. We see from Expression (8) that both \(\gamma _{k,B}(\tau)\) and \(\sqrt {\beta _k(\tau)}\sup _{{\bf x}\in \mathcal {X}}\sigma _{k,B}({\bf x})\) must converge to zero as \(k,B \rightarrow \infty\) to guarantee a vanishing approximation error at \(\forall {\bf x}\in \mathcal {X}\). Given the aforementioned choice of \(\lambda \asymp \exp (-k^{\frac{1}{p}})\) with \(p\gt 2d\), we have that \(\lim \limits _{k, B\rightarrow \infty } \gamma _{k,B}(\tau) = 0\) can be achieved if the grid constant as a function of k, \(\tau (k)\), decreases sufficiently fast. Recall from Expression (12) that \(\gamma _{k,B}(\tau) \lesssim ((\beta _k(\tau) \cdot \tau (k) \cdot k {\sf V}_{k,\min }^{-1})^{\tfrac{1}{2}} + \tau (k) + \tau (k) \cdot k {\sf V}_{k,\min }^{-\tfrac{1}{2}}) .\) Assuming that \({\sf V}_{k,\max }/{\sf V}_{k, \min }\lt \infty\), we have \({\sf V}_{k,\min } \asymp {\sf V}_{k,\max } \asymp k \exp (-k^{\tfrac{1}{p}})\). Setting \(\tau (k) = \mathcal {O}(\exp (-a k^{\tfrac{1}{p}}))\) where \(a \ge 2\) ensures that \(\lim \limits _{k, B\rightarrow \infty } \gamma _{k,B}(\tau)=0\) and moreover \(\gamma _{k,B}(\tau) = o(\sqrt {\beta _k(\tau)}\sup _{{\bf x}\in \mathcal {X}}\sigma _{k,B}({\bf x}))\) as \(k, B\rightarrow \infty\). The latter follows from the fact that \(\beta _k(\tau) = \mathcal {O} (k^{\tfrac{1}{p}})\) due to Expression (9) and \(\sigma _{k,B}({\bf x}) = \mathcal {O} (\exp (-\frac{1}{2}k^{\frac{1}{p}}) k^\frac{d}{2p})\) at \(\forall {\bf x}\in \mathcal {X}\), thanks to Corollary 4. Therefore, \(\gamma _{k,B}(\tau) + \sqrt {\beta _k(\tau)}\sup _{{\bf x}\in \mathcal {X}}\sigma _{k,B}({\bf x}) \lesssim \sqrt {\beta _{k}(\tau)} \sup _{{\bf x}\in \mathcal {X}}\sigma _{k,B}({\bf x}) = \mathcal {O}(\exp (-\frac{1}{2}k^{\frac{1}{p}}) k^\frac{1+d}{2p})\) as \(k, B \rightarrow \infty\).