Analysis of a new class of rational RBF expansions

We propose a new method, namely an eigen-rational kernel-based scheme, for multivariate interpolation via mesh-free methods. It consists of a fractional radial basis function (RBF) expansion, with the denominator depending on the eigenvector associated to the largest eigenvalue of the kernel matrix. Classical bounds in terms of Lebesgue constants and convergence rates with respect to the mesh size of the eigen-rational interpolant are indeed comparable with those of classical kernel-based methods. However, the proposed approach takes advantage of rescaling the classical RBF expansion providing more robust approximations. Theoretical analysis, numerical experiments and applications support our findings.


Introduction
Multivariate approximation is one of the most investigated topics in applied mathematics and finds applications in a wide variety of fields.Many successful methods, such as multivariate splines, meshfree or meshless approaches and finite elements (see e.g.De Boor, 1978;Brenner & Scott, 1994), have already been proven to be effective numerical tools.All these approaches have useful properties and we here concentrate on mesh-free methods (refer e.g. to Buhmann, 2003;Wendland, 2005;Fasshauer, 2007).Concerning approximation theory, their historical foundation lies in the concept of positive definite functions or, more generally, positive definite kernels.Their development can be traced back to the work of both Mercer (1909) and Mathias (1923).Nowadays many positive definite functions are classified as radial basis functions (RBFs), a new term that appeared for the first time (Dyn & Levin, 1983).
Such mesh-free methods, taking advantage of being easy to implement in any dimension, adapt to different applications.Their main branches deal with interpolation, collocation and quasi-interpolation (refer e.g. to Franke, 1982;Kansa, 1990;Buhmann & Dai, 2015).The kernel of this investigation lies in interpolation.The linear spaces spanned by RBFs provide useful properties for multivariate approximation.Nevertheless, it might be advantageous to study approximants in nonlinear spaces generated by RBFs.A few examples of these kinds of approaches already exist in the literature and are known as rational RBFs, introduced in Jakobsson et al. (2009) and further developed in De Marchi et al. (2019), Perracchione (2018) and Sarra & Bay (2018).The scheme essentially consists in considering quotient expansions of two kernel-based interpolants.The authors show that they provide more accurate results for functions with steep gradients, in analogy to rational polynomial approximation.However, the rational basis is constructed by means of function values and consequently it is not a data-dependent method.
To avoid this drawback, we propose a new approach, namely what we now call the eigen-rational method, which is extended to work with conditionally positive definite kernels and enables us to define a rational RBF expansion depending exclusively on kernel and data.Roughly speaking, it consists in rescaling the classical RBF interpolant, taking into account the eigenvector associated to the largest eigenvalue of the kernel matrix.Being rescaled, it shows strong similarities with the approach presented in De Marchi et al. (2017).However, the latter paper limits its attention to compactly supported RBFs (CSRBFs).
Moreover, despite the fact that from the analysis of Lebesgue functions and error bounds it turns out that the eigen-rational method behaves similarly to classical kernel-based interpolants, weighting the interpolant numerically provides more accurate approximations.Such a phenomenon is particularly evident for RBFs characterized by a fast decay, such as the Gaussian.This is shown by means of numerical evidence and applications to real data.
The paper is organized as follows.In Section 1 we review the main theoretical aspects of kernelbased interpolation.Section 2 is devoted to the presentation of the eigen-rational interpolant and to its theoretical analysis.Extensive numerical tests and applications to image registration are presented in Sections 3 and 4, respectively.Concluding remarks and a brief outline of future work are discussed in Section 5.

Review of kernel-based interpolation
In this section we review the main theoretical aspects of kernel-based interpolation methods.For further details we refer the reader to the books Buhmann (2003), Fasshauer (2007), Fasshauer & McCourt (2015) and Wendland (2005).The approximation problem we consider can be formulated as follows.Given X N = {x i , i = 1, . . ., N} ⊂ Ω, which should be a set of distinct data points (or data sites or nodes), arbitrarily distributed on a domain Ω ⊆ R d , with an associated set F N = {f i = f (x i ), i = 1, . . ., N} of data values (or measurements or function values), which are obtained by sampling some (unknown) function f : Ω −→ R at the nodes x i , the scattered data interpolation problem consists in finding a function P f : Ω −→ R from our vector space (e.g., spanned by shifts of our RBF; see below) such that (1.1)

Well-posed problems
We restrict our attention to conditionally positive definite radial kernels called K : Ω −→ R of order on R d ; see e.g.Wendland (2005, Def. 10.14, p. 141).Denoted by Π d −1 , the set of polynomials of total degree less than or equal to − 1 in d variables, conditionally positive definite radial kernels are so that for all c ∈ R d \ {0}, e −εr (15 112/45r 9/2 + 16/3r 7/2 − 7r 4 − 14/15r 2 + 1/9 0 for all p ∈ Π d −1 and To each conditionally positive definite kernel we can associate a univariate function φ : [0, ∞) → R (possibly depending on a shape parameter ε > 0), such that for all x, y ∈ Ω, Here, r = ||x − y|| 2 .For several examples of RBFs and their regularities we refer the reader to Table 1.Such functions will be used for the numerical experiments.Usually, since C ∞ functions might lead to instability (see e.g.Fasshauer, 2007), RBFs with finite regularities are strongly recommended in applications.Moreover, observe that the Buhmann functions (Buhmann, 2000;Zastavnyi, 2006;Porcu et al., 2017;Zastavnyi & Porcu, 2017) are independent of the shape parameter and, as the Wendland ones, belong to the class of CSRBFs.Indeed, the Buhmann kernels listed below are defined for 0 ≤ r ≤ 1.
An RBF interpolant P f : Ω −→ R assumes the form (Fasshauer, 2007) where p 1 , . . ., p L are a basis for the L-dimensional linear space Π d −1 , with Since the conditions (1) must be satisfied, the interpolation problem (3) leads to solving a linear system of the form where T , and the 0 denotes a zero vector of length L and, finally, O is an L × L zero matrix.
In order to establish sufficient conditions under which the interpolation problem is well posed (uniquely solvable), we give the following definition (cf.Fasshauer, 2007, Definition 6.1, p. 53).
data points is called an ( − 1)-unisolvent set if the only polynomial of total degree at most − 1 interpolating zero data on X N is the zero polynomial.
In this situation the following theorem provides suitable conditions under which the stated approximation problem admits a unique solution; we refer e.g. to Fasshauer (2007, Theorem 7.2, p. 64). 3) is conditionally positive definite of order on R d and the set X N = {x i , i = 1, . . ., N} ⊂ Ω of data points forms an ( − 1)-unisolvent set, then the system of linear equations (4) admits a unique solution.
Remark 1.3 If = 0, then we have conditionally positive definite functions of order zero, i.e. strictly positive definite functions.As a consequence, we can interpolate without adding the polynomial term in (3).

Error bounds
In this subsection, we report error estimates for the interpolation process based on conditionally positive definite kernels.To this aim, we need to introduce the so-called native spaces.For each conditionally positive definite kernel K with respect to Π d −1 , we define an associated real Hilbert space, the native space N K (Ω).Following Wendland (2005), we start by introducing the space Since H K (Ω) is only a pre-Hilbert space, we define a Hilbert space H K (Ω) associated to the kernel K as the completion of H K (Ω) with respect to the norm To get the native space N K (Ω), we consider the ( − 1)-unisolvent set E L = {ξ 1 , . . ., ξ L } ⊂ Ω and the Lagrange basis p j , j = 1, . . ., L of Π d −1 constructed on these centres.Then we consider The following definition and the successive theorem (cf. Wendland, 2005, Definition 10.16, p. 144 and Theorem 11.4, p. 176) provide a first error bound in terms of the power function P K,X N .
Definition 1.4 The native space corresponding to a symmetric kernel that is conditionally positive definite on Ω with respect to Π d −1 is defined by equipped with a semiinner product via A first important property for the interpolant in the native spaces is that it has the minimum Hilbert space norm among all functions interpolating the given data, i.e.
which in particular shows the boundedness of the interpolant.Other bounds in terms of L ∞ norms can be found in De Marchi & Schaback (2010) and Rieger & Wendland (2017).
Theorem 1.5 Let Ω ⊆ R d be open.Suppose that K ∈ C 2k (Ω × Ω) is a conditionally positive definite kernel of order .Assume that X N = {x i , i = 1, . . ., N} ⊂ Ω is an ( − 1)-unisolvent set.Then, for every α ∈ N d 0 with |α| ≤ k, the error between f ∈ N K (Ω) and its interpolant can be bounded by Theorem 1.5 bounds the error with respect to the power function, i.e. a term independent of f and dependent on K and X N .Moreover, to test the convergence of the method, we need to investigate error bounds in terms of mesh size.Since we deal with scattered data we consider the following quantity, known as the fill-distance; see Fasshauer (2007): It is an indicator of how Ω is filled out by points.We now report the following error bound in terms of the fill-distance (refer e.g. to Fasshauer, 2007, Theorem 14.6, p. 123).
Theorem 1.6 Suppose Ω ⊆ R d is open and bounded and satisfies an interior cone condition.Let K ∈ C 2k (Ω × Ω) be symmetric and conditionally positive definite of order .Fix α ∈ N d 0 with |α| ≤ k.
Then there exist positive constants h 0 and C independent of x, f and K, such that where Fasshauer (2007, Theorem 14.4, p. 120), and where Remark 1.7 Referring to Theorem 1.6 and to its proof (see also Wendland, 2005, Theorem 11.3, p. 182), we can see that the error bound is obtained by means of a Taylor expansion of order α of the power function about the origin.In particular, if |α| = 0, we have (refer also to Fasshauer, 2007, Theorem 14.5, p. 121) Theorem 1.6 provides convergence orders for a kernel-based interpolant.Essentially, as long as the fill-distance decreases and as long as the effect of ill-conditioning is negligible, the error decreases.Moreover, note that equations ( 6) and ( 7) provide an upper bound for the power function, indeed for all x ∈ Ω (see e.g.Fasshauer, 2007, p. 121 ), (1.9) Remark 1.8 We finally point out that if all the derivatives of order 2k are continuous on ( Ω × Ω) then C K (x) is uniformly bounded Ω; see Wendland (2005, p. 183).

Eigen-rational kernel-based interpolation
In this section we provide a new definition for a rational RBF expansion.This study is motivated by the analogy with polynomial rational approximation, which is well known to be particularly well performing.However, it is a mesh-dependent approach and, as a consequence, extending polynomial approximation in higher dimensions is quite hard (refer e.g. to Hu et al., 2002;Lehmensiek & Meyer, 2011).This is the main reason for which recent research focuses on rational RBFs.
A rational expansion for RBF-based methods was first introduced in Jakobsson et al. (2009).Later this scheme was extended to work with the partition of unity method (De Marchi et al., 2019) and for collocation problems (Sarra & Bay, 2018).This algorithm has already been proven to be an effective numerical tool for interpolation.However, the main drawback consists of the fact that it does not provide a basis independent of the function values.On the contrary, the method described in what follows enables us to provide stability and convergence analysis.

The eigen-rational interpolant
The new class of rational kernel-based interpolants we propose, namely eigen-rational interpolants, assumes the form defined for some function values g i , h i , i = 1, . . ., N. Roughly speaking, once we provide the function values P h (x i ) = h i , i = 1, . . ., N, we can construct P g in the standard way, i.e. such that it interpolates g = (f 1 h 1 , . . ., f N h N ) .Then, obviously, Pf interpolates the given function values F N at the nodes X N .Note that we allow the more general case of conditionally positive definite functions for constructing P g , but we need to restrict to the case of strictly positive kernels for P h .As evident from what follows, this enables us to make the eigen-rational interpolant well defined, i.e. such that P h (x) = 0 for all x ∈ Ω.Unfortunately, such a choice implies that we might consider two different RBFs in (10), i.e. two different kernel matrices.Indeed, let us suppose that K is a conditionally positive definite kernel; then we define K as its associated kernel, which turns out to be strictly positive definite.Since many conditionally positive definite functions φ(r) of order inherit their properties from the th derivative of φ( √ r), if the numerator is a conditionally positive definite function of order , then the associated kernel that we choose for the denominator is any nonzero multiple of For instance, if K is the GM of Table 1, then the associated K is the IM of Table 1.Finally, if K is strictly positive definite, we fix K = K so that we deal with the same kernel matrix for both numerator and denominator.
To show under which conditions the rational interpolant is well defined, i.e. such that P h (x) = 0 for all x ∈ Ω, we need to introduce Perron's theorem (Perron, 1907).
Definition 2.1 A square matrix is called positive when all entries are positive.
By means of this definition, the following statements find application in our setting (cf.Shapiro, 2016, Theorem 1.8, p. 13 and Theorem 1.10, p. 13).
Theorem 2.2 Every positive square matrix has an eigenpair with a positive eigenvalue and an eigenvector with all positive components.The eigenpair, eigenvalue and eigenvector, of the above theorem is known as Perron's eigenpair.
Theorem 2.3 Each positive square matrix possesses exactly one Perron eigenpair.Among all the (real or possibly complex) eigenvalues, the Perron eigenvalue has the largest modulus.
We are now able to discuss when the rational interpolant is well defined.The coefficients β = (β 1 , . . ., β N ) are selected so that β is the eigenvector associated to the largest eigenvalue of Φ K , i.e. such that it satisfies max where Φ K is the kernel matrix generated by K. Of course (11) reduces to finding the eigenvector β associated to the largest eigenvalue of Φ K .
Thus, in our case, if the kernel K is such that K(x, •) > 0, for all x ∈ Ω, by virtue of Theorem 2.3, we can conclude that P h (x) = 0 for all x ∈ Ω.This is, for instance, the case for all the globally defined functions listed in Table 1.We can also accommodate the CSRBFs under certain restrictions.In particular, for the Wendland functions, we know that there always exists a shape parameter ε m such that the kernel matrix generated by the Wendland functions is positive.However, working in this way, we lose the sparsity of the kernel matrix for CSRBFs, which might be essential for huge data sets.Therefore, in the case of CSRBFs one can simply define β as in (11) subject to the constraint that β i be positive for all i = 1, . . ., N.
Remark 2.4 The choice made in (11) follows from the fact that the aim consists in finding P h such that its native space norm relative to the function values is minimized.This would lead to finding min β β Φ −1 K β.Furthermore, dividing the interpolant by the eigenvector associated to the largest eigenvalue of Φ K should make the computation more accurate.Moreover, it enables us to give an eigenrational RBF expansion independent of the function values of the approximant and depending only on the kernel and on the data set.
Once we compute the vector β, we calculate the function values P h (x i ) = h i , i = 1, . . ., N.Then, assuming that K is a conditionally positive definite radial kernel of order on R d , we construct the eigen-rational interpolant Pf by solving (4), where the vector of function values is replaced by (g, 0) , with g = fh and 0 is a zero vector of length L.
In what follows, we point out that this eigen-rational interpolant is comparable in terms of accuracy and stability to the standard RBF approximation.However, since we consider in the denominator an interpolant constructed with the eigenvector associated to the largest eigenvalue, we expect a substantial improvement in terms of numerical stability.
Unlike Jakobsson et al. (2009), we here define the vector h considering only the interpolation matrix Φ K , i.e. depending exclusively on the nodes and on the RBF, and not on the function values.This enables us to construct a data-dependent rational basis.Indeed, let us consider for simplicity a strictly positive definite kernel K; then the eigen-rational interpolant Pf reads as follows: , where we used Since P h is not vanishing, the function is strictly positive definite (see Aronszajn, 1950;De Marchi et al., 2017).In the next subsections, for the eigen-rational interpolant, we investigate the theoretical behaviour of Lebesgue functions and we provide generic error bounds from above.The technique for computing error estimates is then extended to work with the eigen-rational approximant.

Stability analysis
Lebesgue constants are stability indicators that should be studied for a better comprehension of the error behaviour.To this aim we need to dispose of the cardinal form for our eigen-rational interpolant.
Theorem 2.5 Suppose K is a conditionally positive definite kernel of order in R d and K is the associated strictly positive definite kernel.Suppose Furthermore, if K = K is strictly positive definite, the ûk , k = 1, . . ., N form a partition of unity.
Proof.If K is a conditionally positive definite kernel with respect to Wendland (2005, Theorem 11.1, p. 173), we know that there exist functions If the kernel K is strictly positive definite, the same argument holds for functions ūk ∈ span{ K(•, x j ), j = 1, . . ., N}.Therefore, Thus, the resulting eigen-rational interpolant in cardinal form is given by and, furthermore, ûi = 1, as required.Once we have a cardinal form for the eigen-rational interpolant, we can define the Lebesgue functions and constants.Indeed, We note that, as in De Marchi et al. (2017), the cardinal functions of the eigen-rational interpolant are essentially a rescaled form of those of the standard RBF approximant.The same holds for the Lebesgue functions and so we expect that their behaviour is similar to that of the standard interpolants.
Remark 2.6 Given x ∈ Ω \ X N , we have and from Wendland (2005, Theorem 12.1, p. 208), we know that where ω is the smallest eigenvalue of the kernel matrix constructed on the node set X N ∪ {x}.Since the quadratic form ( 2) is positive, we necessarily have ω > 0. This, together with ( 5) and ( 9), shows an upper bound for the Lebesgue functions ûj .

Error analysis
To formulate error bounds, we need to think of g and h as function values obtained by sampling two functions g = fh and h belonging to N K (Ω) and N K (Ω), respectively.We will show that the error for the eigen-rational interpolant can be bounded in terms of both power function and fill-distance.
The following proposition bounds the error from above for the eigen-rational interpolant in terms of the power function.
Proposition 2.7 Let Ω ⊆ R d be open.Suppose K ∈ C(Ω × Ω) is a conditionally positive definite kernel of order and K the associated strictly positive definite kernel.Assume that X N = {x i , i = 1, . . ., N} ⊂ Ω is an ( − 1)-unisolvent set.Then, for x ∈ Ω, the error between f and its eigen-rational interpolant can be bounded by Proof.Let us consider x ∈ Ω; then we have Then by applying Theorem 1.5 we have as claimed.
The statement of Proposition 2.7 enables us to bound the error in terms of the power functions defined on X N associated to K and K.However, for testing the convergence of the method, we need to investigate how the errors depend on the fill-distance.
Proposition 2.8 Suppose Ω ⊆ R d is open, bounded and satisfies an interior cone condition.Further, let K ∈ C 2k (Ω × Ω) be symmetric and conditionally positive definite of order and K be the associated strictly positive definite kernel.Assume that X N = {x i , i = 1, . . ., N} ⊂ Ω is an ( − 1)-unisolvent set.Then there exist positive constants h 0 and C such that provided h X N ≤ h 0 and g, h ∈ N K (Ω) and N K (Ω), respectively, where and C K (x) and C K (x) are computed as in (8).
Proof.The proof is straightforward.Indeed, by using the same argument used for Proposition 2.7 and Theorem 1.5, for certain constants C1 and C2 , we have Therefore, by defining C = max( C1 , C2 ) the assertion follows.
Proposition 2.7 simplifies as soon as we consider strictly positive definite kernels for constructing P g .Indeed, we then interpolate with the same kernel, which in particular implies that P K,X N (x) = P K,X N (x).Furthermore, we will empirically verify the bound given by Proposition 2.8.Theoretically, the eigen-rational interpolant and the standard one show similar convergence rates.
Moreover, we point out that the error depends on the shape parameter of the basis function.Indeed, even if its dependence is not explicitly indicated, the kernel matrices are constructed by means of this scaling parameter.Usually, its value is selected via a priori error estimates.Thus, in what follows we propose a method to determine an optimal value for the shape parameter of an eigen-rational interpolant.We use the term optimal with abuse of notation.Indeed, error estimates only provide an approximation of the (true) optimal value, i.e. the one that can be found via trial and error for which the solution is supposed to be known.

Error estimates
Techniques allowing us to select a predicted optimal shape parameter via error estimates have already been designed.Here we focus on the well-known cross-validation algorithm.It was introduced in Allen (1964) and further developed in Golub & Von Matt (1997).An efficient variant of such a method for strictly positive definite functions, known in literature as leave one out cross validation (LOOCV), is detailed in Rippa (1999) and adapted to the eigen-rational interpolant in what follows.In fact, the polynomial term does not play a crucial role in the accuracy of the interpolant and therefore, only for computing error estimates, we omit the polynomial part.For completeness, we have to mention that there exist methods that are able to handle the polynomial term for computing error estimates (see e.g.Scheuerer, 2011;Biazar & Hosami, 2017).Further extensions of LOOCV for huge data sets are also available in the literature (Cavoretto et al., 2018).
To adapt LOOCV to our setting, let us first define for a classical RBF interpolant the following cost for a node x k (Rippa, 1999): where P k f is computed leaving out the kth node.From Rippa (1999) we know that, if the kernel is strictly positive definite, (12) consists of calculating (2.4) Note that this is the key step for having an efficient method.More specifically, this simplification implies that we compute the inverse of the interpolation matrix only once, while using ( 12) it would have to be computed N times, i.e. for each node left out.Therefore, even if ( 12) applies to the more general case of conditionally positive definite kernels, we prefer to simplify the problem and efficiently compute a measure of the error via (13).
Then for a classical interpolant we define the error estimate as in Rippa (1999) to be where at this point we fix ν = ∞.
For the eigen-rational interpolant we have .
Then the error estimate for the eigen-rational interpolant is defined as Since Ê depends on the shape parameter, i.e.Ê = Ê(ε), the criterion enables us to also select an a priori optimal value for it.To evaluate the proposed algorithm we need to compare the estimated optimal shape parameter, namely ε * , with the true one ε * T .

Numerical experiments
In this section we carry out numerical experiments for testing the behaviour of the eigen-rational interpolant.It is compared with the standard one.We analyse the behaviour of the Lebesgue functions and of the convergence rates.Moreover, we test the method for error estimates.We use the RBFs reported in Table 1.It should be clear that if we use the GM kernel, then the eigenrational interpolant is computed by means of the IM at the denominator.We take different kinds of data sets: grid (G), Halton (H), random (R) and Chebyshev (C) points.
The eigen-rational interpolants are evaluated on v = 40 d equally spaced points X v = {x i , i = 1, . . ., v}.To point out the accuracy we refer to the root mean square error ( R) and the maximum absolute error ( Â),

Experiments for Lebesgue functions and constants
In these tests we compare the Lebesgue constants Λ N and functions λ N of the standard RBF interpolant with the corresponding ones of the eigen-rational approximant, i.e.ΛN and λN .The direct computation of the cardinal functions is an unstable operation and thus we take simple examples to avoid the possible effect of the ill-conditioning.
The first experiment consists in computing the Lebesgue constants for different kernels on N = 10 equally spaced points on Ω = [−1, 1].The shape parameter, unless we consider the case of the Buhmann functions, affects the calculation: that is why we repeat the experiment for two values of ε.The results are reported in Table 2.We can note that, as expected from the stability analysis carried out in Section 2.2, the Lebesgue constants for the eigen-rational and standard interpolants are comparable.
To have a better understanding of the phenomenon, we calculate the Lebesgue functions also for different kinds of nodes (Halton, random and Chebyshev); refer to Fig. 1.Since the cardinal form of the eigen-rational interpolant is essentially a rescaled form of the standard one, as expected, we obtain Lebesgue functions for the eigen-rational interpolant that behave similarly to those of the standard method.Similar behaviour can be observed in the two-dimensional framework.In Fig. 2 we show a comparison of the Lebesgue functions for the W6 and B2 kernels.Therefore, from this study it comes out that the behaviour of the Lebesgue functions and constants are comparable.

Experiments for convergence and error estimates
In this section we compare the errors obtained via the classical RBF interpolation and the eigen-rational one.It will be evident that the eigen-rational interpolation gives more accurate results, especially for RBF characterized by fast decay.In particular, the aim on the one hand consists in comparing the error estimates described in Section 2.4 and on the other hand in checking the convergence rates reported in Section 2.2.Therefore, in the first part of this comparison, since we are interested in analysing how the error varies with respect to the shape parameter, we keep N fixed and let ε vary-whereas to test the convergence rates we essentially reverse this setting.
In our first examples we consider N = 81 Chebyshev and random points on Ω = [−1, 1] and the following test functions: In Figs 3 and 4 we plot the maximum absolute error A and Â for standard and eigen-rational interpolation, respectively, by varying ε in the range Ξ = [10 −2 , 10 2 ].In general, the eigen-rational one is more stable in the region close to the optimal shape parameter.Furthermore, especially for the GA function, it turns out to be really accurate.In other words, we numerically show that it enables us to partially overcome the trade-off between accuracy and stability, i.e. a conflict between theoretically achievable accuracy and numerical stability.
Furthermore, in Figs 3 and 4 we also always report the estimated error via LOOCV E and Ê for standard and eigen-rational approximation, as computed in Section 2.4.If we think of E and Ê as error  Table 3 Error estimates and maximum absolute errors with the approximated optimal and true optimal shape parameter (ε * and ε * T , respectively) for the standard interpolant.They are compared with those obtained via eigen-rational approximation (ε * and ε * T ).Results are computed in the same framework as Figs 3 and 4  functions depending on the value of ε, i.e.E = E(ε) and Ê = Ê(ε), we are able to define a priori optimal values for the shape parameter as follows: The estimated values of the optimal shape parameters ε * and ε * for standard and eigen-rational RBF interpolants, are compared with the true ones ε * T and ε * T , respectively.The latter are found via trial and error by minimizing the maximum absolute errors.In other words, we think of A and Â similarly, as error functions depending on ε, A = A(ε) and Â = Â(ε); refer to Table 3.For both eigen-rational and classical interpolants, the LOOCV provides suitable approximation of the optimal values.However, we might note that in correspondence with the optimal shape parameter we usually register higher absolute errors for the standard interpolant.This is particularly evident for the GA function.For comparing the convergence rates of eigen-rational and standard interpolants, we need to study the error by varying the number of nodes N and evaluate the empirical convergence rates where Rk is the root mean square error for the eigen-rational interpolant on the kth numerical experiment, and h X N k is the fill-distance of the kth computational mesh.Similarly, we define the convergence rates Γ k for the standard interpolant.Tests are carried out by considering several nested sets of Halton and grid points on Ω = [0, 1] 2 and the following test functions: The results of these experiments are reported in Tables 4 and 5. Of course, to evaluate the convergence, we need to keep ε fixed for all the data sets.We note that, as expected, the convergence rates of eigen-rational and standard interpolants are comparable.In particular, we recover spectral convergence rates for C ∞ RBFs, while for functions with finite regularities the convergence rates are limited by the order of smoothness of the kernel.Furthermore, note that despite such rates being comparable, the eigen-rational interpolant is revealed to be more accurate, especially for the GA function.For such a particular RBF it is also possible to see faster convergence rates.

Test with real data: application to image registration
Besides the experiments shown in the previous section, we test the method on real data.We consider the application to image registration.The purpose of image registration is to find a transformation so that the transformed form of the source image is similar to the target one.In what follows we focus on landmark-based image registration; see e.g.Cavoretto et al. (2015), Cavoretto & De Rossi (2018) and (MIPAV) for further details.
Let S = {s i ∈ R 2 , i = 1, . . ., N} be a given set of landmarks belonging to the source image and T = {t i ∈ R 2 , i = 1, . . ., N} the corresponding landmarks in the target image; an example taken by Cavoretto & De Rossi (2018) is shown in Fig. 6, where we take 21 landmarks.The problem consists in finding a transformation F : R 2 −→ R 2 such that F(s i ) = t i , i = 1, . . ., N. In particular, if the problem is solved separately for each component of F, the computational issue reduces to the one presented in Section 1.The results of the image registration for the eigen-rational interpolant via W2 and M2 are plotted in Fig. 5 and Fig. 6 .Since it is difficult to carry out a comparison only with the displayed images, we omit the results obtained by means of standard interpolation.However, we make a comparison by evaluating the mean error (M) of the standard interpolant, As usually done, we denote by M the mean error calculated with eigen-rational interpolants.A comparison of the mean errors vs. different values of the shape parameters belonging to Ξ = [10 −3 , 2] is plotted in Fig. 7. Once more, we note that the eigen-rational interpolants perform better than the classical ones.However, we point out that we did not investigated any topology preservation properties for which the classical RBF approach might take advantage; see e.g.(Cavoretto & De Rossi, 2018).It holds for the classical interpolant for several RBFs and in those cases the standard one might take advantage of this.

Conclusions and work in progress
We proposed a new numerical method, the eigen-rational interpolant.For such an approximant we provided a cardinal form and we defined the Lebesgue functions and constants.They are the standard ones rescaled by a factor depending on the largest eigenvalue of the kernel matrix.Moreover, its convergence was studied and numerically we observed that the eigen-rational interpolant provides more accurate approximations.
Nevertheless, further investigations from the scientific community in reconstructing functions with steep gradients and discontinuities are needed.Generally, in these cases the scheme proposed in Sarra & Bay (2018) turns out to be more robust than the eigen-rational approach.Such results call for further studies and comparisons with the method based on discontinuous kernels outlined in De Marchi et al. (2018).

Fig. 3 .
Fig. 3. Error estimates via LOOCV Ê and E and maximum absolute errors Â and A for eigen-rational and classical interpolants, respectively.Results are computed with f 1 and 81 Chebyshev points on Ω = [−1, 1] via GM (left) and GA (right) kernels.

Fig. 5 .
Fig. 5.The image considered in the example.Landmarks are plotted with squares on the source image (left) and with dots on the target image (left).

Fig. 7 .
Fig. 7. Mean error M by varying the shape parameter for the standard interpolant compared with those obtained with the eigenrational interpolant M. Results are computed in the same framework as of Figs 6-7.

Table 2
Lebesgue constants Λ N and ΛN for classical and eigen-rational interpolants, respectively.They are computed on N = 10 equally spaced points on Ω = [−1, 1].We remark that the Buhmann function is independent of ε

Table 4
Root mean square errors R, R, convergence rates Γ and Γ for standard and eigen-rational interpolants, respectively, computed with f 3 on five sets of grid data on