Kernel-based interpolation at approximate Fekete points

We construct approximate Fekete point sets for kernel-based interpolation by maximising the determinant of a kernel Gram matrix obtained via truncation of an orthonormal expansion of the kernel. Uniform error estimates are proved for kernel interpolants at the resulting points. If the kernel is Gaussian we show that the approximate Fekete points in one dimension are the solution to a convex optimisation problem and that the interpolants convergence with a super-exponential rate. A numerical experiment is provided for the Gaussian kernel.


Introduction
Kernel-based methods are widely used in interpolation and approximation of functions (Wendland, 2005;Fasshauer, 2007;Fasshauer and McCourt, 2015). Let d ∈ N and Ω ⊂ R d be a compact set with a non-empty interior. Given evaluations of a function f : Ω → R at a scattered set of distinct points X n = {x x x 1 , . . . , x x x n } ⊂ Ω and a continuous positive-definite kernel K : Ω × Ω → R, the kernel interpolant s f is where the coefficients c k are uniquely determined by the interpolation conditions s f (x x x k ) = f (x x x k ) for every k = 1, . . . , n. The choice of the evaluation points X n can have a significant effect on the accuracy of the approximation s f (x x x) ≈ f (x x x) at x x x / ∈ X n . Popular methods for constructing "good" point sets include different types of greedy algorithms (Schaback and Wendland, 2000;De Marchi et al., 2005;Müller, 2009;Wirtz and Haasdonk, 2013;Santin and Haasdonk, 2017) that construct the next point x x x n+1 by maximising the power function. An alternative approach is to select n points concurrently by maximising det K Xn = det(K(x x x k , x x x m )) n k,m=1 , the determinant of the kernel Gram matrix, over all sets of n points X n ⊂ Ω. The resulting points are called Fekete points in an analogue to the classical Fekete points that maximise the Vandermonde determinant (Bos et al., 2010;Briani et al., 2012). The asymptotic distribution of these points for kernel-based interpolation in one dimension has been studied by Bos and Maier (2002) and Bos and De Marchi (2011).
Because maximisation of det K Xn is typically intractable, in this article we study approximate Fekete points that are obtained by maximising the determinant of the kernel matrix of a truncated version of the kernel. Let {ϕ } ∞ =1 be an orthonormal basis of H K (Ω), the reproducing kernel Hilbert space (RKHS) of K. Then the kernel can be written as K(x x x, y y y) = ∞ =1 ϕ (x x x)ϕ (y y y).
The approximate Fekete points X * n are then defined as any set of n points maximising . (1.1) This and related constructions have been recently suggested by Tanaka (2019) and in the context of numerical integration and sampling from determinantal point processes by Belhadji et al. (2019) and Gautier et al. (2019). Our construction differs slightly from the prior work in that we do not require the basis functions {ϕ } ∞ =1 to arise from Mercer's theorem, which significantly simplifies analysis and construction of the points, at least when the kernel is Gaussian. This article contains two main contributions: • Let f ∈ H K (Ω). In Section 3 we use a bound on the Lebesgue constant for interpolation with {ϕ } n =1 to prove that for interpolation at any approximate Fekete points.
• In Section 4 we show that for a certain simple orthonormal expansion (Minh, 2010) of the univariate Gaussian kernel with a scale parameter ε > 0 the objective function (1.1) is convex and has a unique maximiser. We then specialise the uniform error estimate (1.2) and some other results from Section 3 for the Gaussian kernel.
A numerical experiment for the Gaussian kernel is given in Section 5. We also discuss improved error estimates in subspaces of H K (Ω) and tensor product extensions of the univariate approximate Fekete points for anisotropic multivariate Gaussian kernels.

Background
This section reviews basic properties of kernel interpolants and defines the approximate Fekete points studied in the remainder of the article.

Kernel-based interpolation
Every positive-definite kernel K : Ω × Ω → R on a general domain Ω ⊂ R d induces a unique reproducing kernel Hilbert space H K (Ω), which is a Hilbert space consisting of real-valued functions defined on Ω. The RKHS is characterised by the properties that K(·, x x x) ∈ H K (Ω) for every x x x ∈ Ω and f, K(·, x x x) H K (Ω) = f (x x x) for every f ∈ H K (Ω) and x x x ∈ Ω, the latter of which is known as the reproducing property.
Given a set of n distinct points, X n = {x x x 1 , . . . , x x x n } ⊂ Ω, the kernel interpolant s f is the minimum-norm interpolant to a function f : Ω → R at these points: This definition implies that s f H K (Ω) ≤ f H K (Ω) . The main advantage in working in an RKHS as opposed to some different function space is that the minimum-norm interpolant has a simple algebraic form: where we denote c c c = (c 1 , . . . , c n ) ∈ R n and k k k is the positive-definite kernel Gram matrix. From this it follows that s f is the unique interpolant to f at X n in the span of {K(·, x x x k )} n k=1 . The interpolant can be written using the cardinal functions From the reproducing property and the Cauchy-Schwarz inequality it now follows that for any f ∈ H K (Ω) the interpolation error admits the bound where the non-negative power function P Xn is The latter form is the point-wise worst-case approximation error. The power function can be also written in a determinantal form (e.g. Schaback, 2005, Lemma 3) which suggests, via (2.3), that points X n that maximise det K Xn ought to provide to small approximation error. Numerous explicit bounds on the error f − s f in different norms and for different classes of kernels and functions within and without the RKHS can be found in Wendland (2005, Chapter 11); Wendland and Rieger (2005);Narcowich et al. (2006); Arcangéli et al. (2007) and Rieger and Zwicknagl (2010).

Approximate Fekete points
For the remainder of this article we assume that Ω is a compact subset of R d with a non-empty interior and that the positive-definite kernel K : Ω × Ω → R is continuous.
Let {ϕ } ∞ =1 be an orthonormal basis of H K (Ω). Then the kernel can be written as for all x x x, y y y ∈ Ω. It is easy to verify that K in (2.5) is the reproducing kernel: Any The Fekete points for interpolation with the kernel (2.5) are the points that maximise the determinant of the kernel matrix. As exact computation of the Fekete points is typically challenging, we fix an orthonormal basis {ϕ } ∞ =1 of H K (Ω), truncate the expansion (2.5) after n terms and consider maximisation of the resulting approximation of the objective function (2.6). Define the truncated kernel and its kernel matrix K Xn = ( K(x x x k , x x x m )) n k,m=1 ∈ R n×n . From (2.7) it is easy to see that The approximate Fekete points X * n = {x x x * 1 , . . . , x x x * n } ⊂ Ω are then any points such that Note that because {ϕ } n =1 are linearly independent, there exists X n ⊂ Ω such that det Φ Xn > 0. As Ω is compact and the continuity of K implies the continuity of the basis functions, there exist points X * n at which det Φ Xn attains a maximal value. Given a set X n of n previously selected points, the popular P -greedy algorithm (De Marchi et al., 2005;Santin and Haasdonk, 2017) which, using the block determinant identity and (2.4), can be written equivalently as That is, the P -greedy points can be interpreted as greedily computed Fekete points. Because it is known (Santin and Haasdonk, 2017) that the interpolation error of the P -greedy algorithm decays fast (in some cases with an optimal rate), it is reasonable to expect that these rates are inherited by interpolation at the Fekete points, and by extension perhaps by interpolation at the approximate Fekete points.

Error estimates
This section provides upper bounds on the error of approximating f ∈ H K (Ω) with the kernel interpolant s f when the interpolation points are the approximate Fekete points from Section 2.2.

Interpolation with basis functions and Lebesgue constants
For any f : Ω → R and any points From these requirements it follows that Alternatively, the interpolant can be written in the Lagrange form where u ϕ k are the Lagrange basis functions solved from for every x x x ∈ Ω. Define the Lebesgue constant A standard argument yields a conservative upper bound on the Lebesgue constant at approximate Fekete points (Bos et al., 2010).
Proposition 3.1. If X * n are any approximate Fekete points (2.8), then the Lebesgue Proof. Cramer's rule applied to (3.3) gives is obtained by replacing the kth row of the matrix Φ Xn with the row vector (ϕ 1 (x x x), . . . , ϕ n (x x x)) ∈ R n . Because any approximate Fekete points maximise det Φ Xn among all sets of n points within Ω and Φ k (2010) for bounds on the Lebesgue constant sup x x x∈Ω n k=1 |u k (x x x)| for kernel interpolation when the RKHS is a Sobolev space.

Uniform error estimates
In this section we consider the uniform interpolation error sup That is, H K (Ω) consists of functions having the form (3.7) such that their norm in (3.8) is finite.
Proof. The claim follows from the Cauchy-Schwarz inequality and (3.8): Then for any f ∈ H K (Ω), The first term on the right-hand side can be bounded with Lemma 3.2. The second term vanishes because g ∈ span{ϕ } n =1 and s ϕ g being the unique interpolant to g in span{ϕ } n =1 imply that s ϕ g = g. Finally, the Lagrange form (3.2) and Lemma 3.2 yield a bound on the third term: To obtain a bound on |f ( . . , n and s f H K (Ω) ≤ f H K (Ω) by the norm-minimality property (2.1), both terms on the right-hand side obey the bound (3.10). The claim follows.
Proposition 3.1 immediately yields an error estimate for any approximate Fekete points.
Corollary 3.4. Suppose that f ∈ H K (Ω) is interpolated at any approximate Fekete points (2.8). Then (3.11) Due to the presence of a supremum on the right-hand side of (3.9) and (3.11) it is difficult to make the bounds explicitly dependent on, for example, smoothness of the kernel as is usual in the error analysis of radial basis function interpolants (Wendland, 2005, Chapter 11).

Improved error estimates in subspaces
It is known that the rate of convergence of kernel interpolation can be improved if the function being interpolated lives in a subset of the RKHS. The existing results in Schaback (1999Schaback ( , 2000Schaback ( , 2018 and Wendland (2005, Section 11.5) are particularly interesting when the kernel is finitely smooth 1 . Roughly speaking, in this case a typical algebraic rate of convergence is "doubled" for sufficiently smooth elements of the RKHS. Specifically, let µ be a Borel measure on Ω that assigns positive measure to every open set and let {ψ } ∞ =1 and (λ ) ∞ =1 be the eigenfunctions and the positive decreasing eigenvalues of the integral operator T f (x x x) = Ω K(x x x, y y y)f (y y y) dµ(y y y). By Mercer's theorem (e.g. Sun, 2005), The standard improved error estimate states that for Because the range of T is the collection of functions for which (3.12) holds is a subset of the RKHS. Theorem 3.5 below is significantly more flexible than this result and does not require that the Mercer expansion be used. Let (α ) ∞ =1 be a positive, increasing, and divergent sequence and define the subspace It is easy to verify that H α K (Ω) is an RKHS and that its reproducing kernel is Theorem 3.5. If f ∈ H α K (Ω) is interpolated at any approximate Fekete points (2.8), then Proof. When f ∈ H α K (Ω), we replace the estimate of Lemma 3.2 with the following estimate: The proof of Theorem 3.3 then yields claimed uniform bound.

Gaussian kernel
For the most of this section we set d = 1 and consider the Gaussian kernel with a scale parameter ε > 0. This kernel has the orthonormal expansion 2 (4.2) and its RKHS is thus These results and other properties of the Gaussian kernel and its RKHS are studied in more detail in Steinwart et al. (2006) and Minh (2010). In Section 4.1 we show that, owing to the special structure of the above basis functions, the approximate Fekete points for the one-dimensional Gaussian kernel are solved from a convex optimisation problem. The rest of this section is then devoted to specialising the general error estimates of the previous section. Note that most prior work, such as Tanaka (2019) and Belhadji et al. (2019), uses a well-known Mercer expansion of the Gaussian kernel instead of (4.2). This expansion is where the eigenfunctions are orthonormal with respect to the Gaussian measure with variance σ 2 : The eigenfunctions and values are (Fasshauer and McCourt, 2012) where H is the th probabilists' Hermite polynomial and the constants are α = 1 √ 2σ , β = (1 + 8ε 2 σ 2 ) 1/4 and δ 2 = 1 4σ 2 (β 2 − 1).
The expansion (4.2) used in this article is evidently much simpler to work with.

Approximate Fekete points via convex optimisation
Let Observe that in this section we begin indexing of the expansion from zero to simplify notation.
be the truncation of the Gaussian kernel (4.1) and K Xn = ( K(x k , x m )) n k,m=1 ∈ R n×n the corresponding kernel matrix. Define the matrices the latter of which is the classical Vandermonde matrix. Because K Xn = Φ Xn Φ T Xn and the kth row of Φ Xn is that of V Xn multiplied by (2 k−1 ε 2(k−1) /(k − 1)!) 1/2 exp(−ε 2 x 2 k ), we have where the last equation uses the standard explicit expression for the Vandermonde determinant. This expression verifies that K Xn and Φ Xn are invertible whenever the points are distinct. Define The approximate Fekete points (2.8) for the Gaussian kernel are thus seen to be x where Q ε (x) = ε 2 x and N (x) = 1/ log |x|. To ensure that I is well-defined and nonuniqueness arising from ordering of the points is eliminated, define the simplex R n = (x 1 , . . . , x n ) ∈ Ω n : x 1 < x 2 < · · · < x n−1 < x n ⊂ Ω n and consider I as a function on R n . Adaptation of the proof of Theorem 3.3 in Tanaka and Sugihara (2019) shows that the objective function I is convex and that there exists a unique minimiser X * n ∈ R n .
If Ω ⊂ R is a closed interval, then the energy function I : R n → R is convex and has a unique minimiser.
Proof. The Hessian matrix ∇ 2 I of I is

Because both
N (x) = log 1 |x| and Q ε (x) = ε 2 x 2 are strictly convex on R\{0} and R, respectively, we have N > 0 and Q > 0. Therefore all diagonal elements of ∇ 2 I are always positive. Moreover, which verifies that the Hessian is diagonally dominant and hence positive-definite. That is, the energy function I is convex on R n . To verify that there is a unique minimiser in the non-closed set R n , consider the function J(X n ) = exp(−I(X n )) which is continuous on the closure of R n if we set J(X n ) = 0 for every X n = {x 1 , . . . , x n } ∈ Ω n such that x i = x i+1 for some i. Because it is positive on R n , any maximiser of J is in R n . As a maximiser of J is a minimser of I and I is convex it follows that I must have a unique minimiser in R n .
Remark 4.2. If we set ε = 0, the above optimisation problem becomes that of finding the Fekete points for polynomial interpolation. However, in this case the objective function I is no longer convex because Q ε (x i ) = 0 in (4.6). Our optimisation problem can be thus viewed as a regularised version of the standard Fekete problem. Based on this and the well-known convergence of kernel interpolants to polynomial interpolants at the so-called flat limit (Schaback, 2005;Lee et al., 2007;Karvonen and Särkkä, 2019) it may be expected that X * n converge to the polynomial Fekete points as ε → 0. We do not attempt to prove this.
Using the estimate of Lemma 4.3 in Corollary 3.4 yields an explicit error estimate for interpolation with the Gaussian kernel.
Proof. The claim follows from Corollary 3.4, Lemma 4.3 and the lower bound n! ≥ √ 2πn n+1/2 e −n of Stirling's approximation (Robbins, 1955): √ 2 εc Ω n e n/2 n n n/2+1/4 = 128 π Also Theorem 3.5 can be specialised, and in some cases the kernel of the subspace H α K (Ω) has an explicit form. For instance, set α = √ ! 2 ε 2 . Then which can be written in terms of I 0 , the modified Bessel function of the first kind: Theorem 4.5. Consider the Gaussian kernel with the orthonormal expansion (4.2) and suppose that Ω ⊂ R is a closed interval. If f ∈ H α K (Ω) is interpolated at the unique approximate Fekete points X * n defined in (4.5) and n ≥ 2ε 2 c 2 Ω , then where C 1 = (128/π) 1/4 ≈ 2.53 and C 2 = √ 2 e εc Ω .
If Ω = [a, b] ⊂ R is a closed interval, the standard fill-distance based bound (Rieger and Zwicknagl, 2010, Theorem 6.1) for interpolation error is when the fill-distance h Xn,Ω = sup x∈Ω min x k ∈Xn |x − x k | is sufficiently small. The constant satisfies C ≤ 1 8 min{(b − a)/6, 1}. 3 For the equispaced points X n = a, a + b − a n , . . . , b − b − a n , b , which have the minimal fill-distance h Xn,Ω = (b − a)/n, the bound (4.8) becomes where C/(b − a) ≤ 1 48 . Our bound (4.7) for points X * n , being essentially of order exp(− 1 2 n log n), is thus tighter when n is sufficiently large. However, a significant advantage of bounds of the type (4.8) is that they apply to nested point sets (i.e., X n ⊂ X n+1 for every n ≥ 1). It cannot be expected that the approximate Fekete point sets are nested. Further error estimates for Chebyshev-type nodes that cluster near the boundary are provided in Rieger and Zwicknagl (2014).
Remark 4.6. It is easy to see that in the Gaussian case the Lagrange basis functions in (3.3) can be expressed in terms of the classical polynomial Lagrange functions: Let Λ pol (X n ) = sup x∈Ω n k=1 |l k (x)| be the Lebesgue constant for polynomial interpolation. It follows easily from (4.9) and the boundedness of Ω that there exist C 1 , C 2 > 0 such that for any X n ⊂ Ω. This implies that in Theorem 3.3 the coefficient 1 + Λ ϕ (X n ) can be replaced with 1 + C 2 Λ pol (X n ), which means that convergence results are available if polynomial Lebesgue constants can be controlled (e.g., if X n are the Chebyshev points).

Tensor product algorithms
In this section we provide error estimates for interpolation with anisotropic Gaussian kernels in higher dimensions when the evaluation points are constructed as tensor products of the approximate Fekete points (4.5). Besides Beatson (2010) there does not appear to be much work on error estimates for general anisotropic kernels.  and Sloan and Woźniakowski (2018) analyse the L 2 -error of general linear algorithms for functions in the RKHS of an anisotropic Gaussian. Let be a hyper-rectangle and with ε i > 0 an anisotropic Gaussian kernel on Ω. Because the kernel (4.11) is a product of one-dimensional Gaussian kernels, its RKHS is the tensor product of their RKHSs: and Let n 1 , . . . , n d ∈ N and denote N = n 1 × · · · × n d . We take the point set to be a tensor product of approximate Fekete point sets (4.5) for Gaussian kernels K i on Ω i : where X * i,ni ⊂ Ω i stands for the set of n i approximate Fekete points for kernel K i on Ω i . Due to the tensor product structure of the point set and the RKHS, the kernel interpolant s f to f ∈ H K (Ω) can be written as where s i,fi is the kernel interpolant, based on K i , of f i ∈ H Ki (Ω i ) at the points X * i,ni . Theorem 4.7. Consider the multi-dimensional Gaussian kernel (4.11) and suppose that Ω ⊂ R d is a hyper-rectangle of the form (4.10). If f ∈ H K (Ω) is interpolated at the tensor product points X * N defined in (4.12) and n i ≥ 2ε 2 i c 2 Ωi for every i = 1, . . . , d, then where C 1 = (128/π) 1/4 ≈ 2.53 and C i,2 = √ 2 e ε i c Ωi .
Proof. The argument is standard. (4.14) where the notational convention f d+1:d (x x x) = 1 is used. By the reproducing property and the minimum-norm property (2.1), for any i ≤ d and x i ∈ Ω i . Inserting these estimates and the bound (4.7) into (4.14) yields In particular, if n 1 = · · · = n d = n (so that N = n d ) and all Ω i and ε i are equal, the bound of Theorem 4.7 becomes It would be straightforward to generalise Theorem 4.5 to the tensor product setting. 2) based on the Gaussian kernel (4.1). Right column: Auxiliary interpolants (3.1) constructed out of the basis functions (4.2). Also displayed is the theoretical rate (5.2) from Theorem 4.4.