Stable interpolation with isotropic and anisotropic Gaussians using Hermite generating function

TL;DR

This paper introduces a new stabilization algorithm for multivariate Gaussian interpolation, leveraging Hermite polynomial generating functions, which improves accuracy and stability especially in high-precision regimes.

Contribution

It proposes a novel stabilization method based on Hermite generating functions and an analytic cut-off criterion for Gaussian interpolation, enhancing existing techniques.

Findings

New stabilization algorithm for Gaussian interpolation

Automatic adjustment of basis functions via analytic cut-off

Improved stability in high-accuracy regimes

Abstract

Gaussian kernels can be an efficient and accurate tool for multivariate interpolation. In practice, high accuracies are often achieved in the flat limit where the interpolation matrix becomes increasingly ill-conditioned. Stable evaluation algorithms for isotropic Gaussians (Gaussian radial basis functions) have been proposed based on a Chebyshev expansion of the Gaussians by Fornberg, Larsson & Flyer and based on a Mercer expansion with Hermite polynomials by Fasshauer & McCourt. In this paper, we propose a new stabilization algorithm for the multivariate interpolation with isotropic or anisotropic Gaussians derived from the generating function of the Hermite polynomials. We also derive and analyse a new analytic cut-off criterion for the generating function expansion that allows to automatically adjust the number of the stabilizing basis functions.

Equations

$$s(\mathbf{x})=\sum _{k=1}^N{\alpha }_k{\varphi }_k(\mathbf{x}),$$

$$s({\mathbf{x}}_i^{\mathrm{col}})=f_i\text{for}i=1,\dots ,N.$$

$${\Phi }^{\mathrm{col}}\mathbf{\alpha }=\mathbf{f},\text{with}{\Phi }_{ij}^{\mathrm{col}}={\varphi }_j({\mathbf{x}}_i^{\mathrm{col}}).$$

$${\varphi }_k(\mathbf{x})=\exp (-{\epsilon }^2\Vert \mathbf{x}-{\mathbf{x}}_k{\Vert }^2)$$

$${\varphi }_k(\mathbf{x})=\exp (-(\mathbf{x}-{\mathbf{x}}_k{)}^T{E}^{T}E(\mathbf{x}-{\mathbf{x}}_k))$$

$$h_{\ell +1}(x)=2x{h}_{\ell }(x)-2\ell h_{\ell -1}(x),h_0(x)=1,h_{-1}(x)=0.$$

$$|h_{\ell }(x)|\le e^{\frac{x^2}{2}}c\sqrt{2^{\ell }\ell !}.$$

$$\epsilon >0,\gamma >0,t\in (0,1),$$

$$H_{\ell }^{\gamma ,\epsilon ,t}(x)=\frac{t^{|\ell |/2}}{\sqrt{2^{\ell }\ell !}}{h}_{\ell }(\gamma x){\mathrm{e}}^{-{\epsilon }^2{x}^2},$$

$${\varphi }_k(x)={\mathrm{e}}^{-{\epsilon }^2{\left(x-x_k\right)}^2}=\exp \left({\epsilon }^2{\Delta }_k^2\left(\frac{{\epsilon }^2}{{\gamma }^2}-1\right)\right)\sum _{\ell \ge 0}\frac{{\epsilon }^{2\ell }\sqrt{2^{\ell }}}{{\gamma }^{\ell }\sqrt{t^{\ell }\ell !}}{\Delta }_k^{\ell }{H}_{\ell }^{\gamma ,\epsilon ,t}(x-x_0),$$

$$s(x)=\sum _{k=1}^N{\alpha }_k\exp \left({\epsilon }^2{\Delta }_k^2\left(\frac{{\epsilon }^2}{{\gamma }^2}-1\right)\right)\sum _{\ell \ge 0}\frac{{\epsilon }^{2\ell }\sqrt{2^{\ell }}}{{\gamma }^{\ell }\sqrt{t^{\ell }\ell !}}{\Delta }_k^{\ell }{H}_{\ell }^{\gamma ,\epsilon ,t}(x-x_0).$$

$${\mathrm{e}}^{2ba-a^2}=\sum _{\ell \ge 0}\frac{a^{\ell }}{\ell !}{h}_{\ell }(b),\forall a,b\in \mathbb{R}.$$

$$\sum _{\ell \ge 0}\frac{{\epsilon }^{2\ell }}{{\gamma }^{\ell }\ell !}{\Delta }_k^{\ell }{h}_{\ell }(\gamma (x-x_0))=\exp \left(2{\epsilon }^2{\Delta }_k(x-x_0)-\frac{{\epsilon }^4{\Delta }_k^2}{{\gamma }^2}\right)$$

$$\exp \left({\epsilon }^2{\Delta }_k^2\left(\frac{{\epsilon }^2}{{\gamma }^2}-1\right)\right)\sum _{\ell \ge 0}\frac{{\epsilon }^{2\ell }\sqrt{2^{\ell }}}{{\gamma }^{\ell }\sqrt{t^{\ell }\ell !}}{\Delta }_k^{\ell }{H}_{\ell }^{\gamma ,\epsilon ,t}(x-x_0)$$

$$=\exp \left({\epsilon }^2{\Delta }_k^2\left(\frac{{\epsilon }^2}{{\gamma }^2}-1\right)\right)\sum _{\ell \ge 0}\frac{{\epsilon }^{2\ell }}{{\gamma }^{\ell }\ell !}{\Delta }_k^{\ell }{h}_{\ell }(\gamma (x-x_0)){\mathrm{e}}^{-{\epsilon }^2(x-x_0{)}^2}\stackrel{\text{eqref}eq:hg{f}_{e}xpansio{n}_{p}roof}{=}{\mathrm{e}}^{-{\epsilon }^2(x-x_k{)}^2},$$

$$h_{\mathbf{\ell }}(\mathbf{x})=h_{{\ell }_1}(x_1)\cdots \cdot h_{{\ell }_d}(x_d),\mathbf{\ell }\in {\mathbb{N}}^d,\mathbf{x}\in {\mathbb{R}}^d.$$

$${\varphi }_{\mathbf{q}}(\mathbf{x})=\exp (-(\mathbf{x}-\mathbf{q}{)}^T{E}^{T}E(\mathbf{x}-\mathbf{q})),\mathbf{x},\mathbf{q}\in {\mathbb{R}}^d,$$

$$H_{\mathbf{\ell }}^{G,E,t}(\mathbf{x})=\frac{t^{|\mathbf{\ell }|/2}}{\sqrt{2^{|\mathbf{\ell }|}\mathbf{\ell }!}}{h}_{\mathbf{\ell }}(G^T\mathbf{x})\exp (-{\mathbf{x}}^T{E}^{T}E\mathbf{x}),$$

$${\varphi }_{\mathbf{q}}(\mathbf{x})=\exp ({\mathbf{\Delta }}_{\mathbf{q}}^T(\tilde{G}-E^{T}E){\mathbf{\Delta }}_{\mathbf{q}})\cdot \sum _{\mathbf{\ell }\in {\mathbb{N}}^d}\frac{(G^{-1}{E}^{T}E{\mathbf{\Delta }}_{\mathbf{q}}{)}^{\mathbf{\ell }}\sqrt{2^{|\mathbf{\ell }|}}}{\sqrt{t^{|\mathbf{\ell }|}\mathbf{\ell }!}}{H}_{\mathbf{\ell }}^{G,E,t}(\mathbf{x}-{\mathbf{x}}_0),$$

$$\sum _{\mathbf{\ell }\in {\mathbb{N}}^d}\frac{{\mathbf{a}}^{\mathbf{\ell }}}{\mathbf{\ell }!}{h}_{\mathbf{\ell }}(\mathbf{b})=\exp (2{\mathbf{b}}^T\mathbf{a}-{\mathbf{a}}^T\mathbf{a}),\forall \mathbf{a},\mathbf{b}\in {\mathbb{R}}^d.$$

$$\sum _{\mathbf{\ell }\in {\mathbb{N}}^d}\frac{(G^{-1}{E}^{T}E{\mathbf{\Delta }}_{\mathbf{q}}{)}^{\mathbf{\ell }}}{\mathbf{\ell }!}{h}_{\mathbf{\ell }}(G^T(\mathbf{x}-{\mathbf{x}}_0))=\exp (2(\mathbf{x}-{\mathbf{x}}_0{)}^T{E}^{T}E{\mathbf{\Delta }}_{\mathbf{q}}-{\mathbf{\Delta }}_{\mathbf{q}}^T\tilde{G}{\mathbf{\Delta }}_{\mathbf{q}}).$$

$$2{\mathbf{x}}^T{E}^{T}E\mathbf{q}$$

$$-{\mathbf{x}}^T{E}^{T}E\mathbf{x}$$

$$E_{\mathrm{isotropic}}=\epsilon {\mathrm{Id}}_d.$$

$$\sum _{|\mathbf{\ell }|=0}^{\infty }\frac{t^{|\mathbf{\ell }|}{h}_{\mathbf{\ell }}(\mathbf{x})h_{\mathbf{\ell }}(\mathbf{y})}{2^{|\mathbf{\ell }|}\mathbf{\ell }!}=\frac{\exp \left(\frac{t}{1-t^2}({\mathbf{x}}^T\mathbf{y}+{\mathbf{y}}^T\mathbf{x})-\frac{t^2}{1-t^2}(\Vert \mathbf{x}{\Vert }_2^2+\Vert \mathbf{y}{\Vert }_2^2)\right)}{(1-t^2{)}^{d/2}}.$$

$${\mathrm{e}}^{{\mathbf{x}}^T\mathbf{x}}={\pi }^{d/2}{\int }_{{\mathbb{R}}^d}{\mathrm{e}}^{2\pi \mathrm{i}{\mathbf{x}}^T\mathbf{\xi }}{\mathrm{e}}^{-{\pi }^2{\mathbf{\xi }}^T\mathbf{\xi }}\mathrm{d}\mathbf{\xi }={\pi }^{-d/2}{\int }_{{\mathbb{R}}^d}{\mathrm{e}}^{2\mathrm{i}{\mathbf{x}}^T\mathbf{\xi }}{\mathrm{e}}^{-{\mathbf{\xi }}^T\mathbf{\xi }}\mathrm{d}\mathbf{\xi }.$$

$$h_{\mathbf{\ell }}(\mathbf{x})={\mathrm{e}}^{{\mathbf{x}}^T\mathbf{x}}(-\nabla {)}^{\mathbf{\ell }}{\mathrm{e}}^{-{\mathbf{x}}^T\mathbf{x}}=\frac{(-2\mathrm{i}{)}^{|\mathbf{\ell }|}}{{\pi }^{d/2}}{\mathrm{e}}^{{\mathbf{x}}^T\mathbf{x}}{\int }_{{\mathbb{R}}^d}{\mathbf{\xi }}^{\mathbf{\ell }}{\mathrm{e}}^{2\mathrm{i}{\mathbf{x}}^T\mathbf{\xi }-{\mathbf{\xi }}^T\mathbf{\xi }}\mathrm{d}\mathbf{\xi },$$

$$\sum _{|\mathbf{\ell }|=0}^{\infty }\frac{t^{|\mathbf{\ell }|}{h}_{\mathbf{\ell }}(\mathbf{x})h_{\mathbf{\ell }}(\mathbf{y})}{2^{|\mathbf{\ell }|}\mathbf{\ell }!}=$$

$$=\frac{{\mathrm{e}}^{{\mathbf{x}}^T\mathbf{x}+{\mathbf{y}}^T\mathbf{y}}}{{\pi }^d}\sum _{|\mathbf{\ell }|=0}^{\infty }\frac{(-2t{)}^{|\mathbf{\ell }|}}{\mathbf{\ell }!}{\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}{\mathbf{\xi }}_{\mathbf{x}}^{\mathbf{\ell }}{\mathbf{\xi }}_{\mathbf{y}}^{\mathbf{\ell }}{\mathrm{e}}^{2\mathrm{i}({\mathbf{x}}^T{\mathbf{\xi }}_{\mathbf{x}}+{\mathbf{y}}^T{\mathbf{\xi }}_{\mathbf{y}})}{\mathrm{e}}^{-{\mathbf{\xi }}_{\mathbf{x}}^T{\mathbf{\xi }}_{\mathbf{x}}-{\mathbf{\xi }}_{\mathbf{y}}^T{\mathbf{\xi }}_{\mathbf{y}}}\mathrm{d}{\mathbf{\xi }}_{\mathbf{x}}\mathrm{d}{\mathbf{\xi }}_{\mathbf{y}}$$

$$=\frac{{\mathrm{e}}^{{\mathbf{x}}^T\mathbf{x}+{\mathbf{y}}^T\mathbf{y}}}{{\pi }^d}{\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}{\mathrm{e}}^{-2t{\mathbf{\xi }}_{\mathbf{x}}^T{\mathbf{\xi }}_{\mathbf{y}}}{\mathrm{e}}^{2\mathrm{i}({\mathbf{x}}^T{\mathbf{\xi }}_{\mathbf{x}}+{\mathbf{y}}^T{\mathbf{\xi }}_{\mathbf{y}})}{\mathrm{e}}^{-{\mathbf{\xi }}_{\mathbf{x}}^T{\mathbf{\xi }}_{\mathbf{x}}-{\mathbf{\xi }}_{\mathbf{y}}^T{\mathbf{\xi }}_{\mathbf{y}}}\mathrm{d}{\mathbf{\xi }}_{\mathbf{x}}\mathrm{d}{\mathbf{\xi }}_{\mathbf{y}},$$

$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{g}_{(x,y)}(u,v)\mathrm{d}v\mathrm{d}u=\frac{\pi {\mathrm{e}}^{-(x^2+y^2)/2}}{\sqrt{1-t^2}}\exp \left(\frac{x^2-y^2}{2}-\frac{(x-yt{)}^2}{1-t^2}\right).$$

$$g_{(x,y)}(u,v)=\exp (-u^2-2tuv-v^2+2\mathrm{i}xu+2\mathrm{i}yv).$$

Full text

1. Introduction

Multivariate interpolation is a topic that is relevant for a vast number of applications. Gaussian radial basis functions (Gaussian RBFs) are a class of functions for which interpolation generalizes to higher dimensions in a simple way while yielding spectral accuracy [6]. However, it is known that rather small values of the shape parameter $\varepsilon>0$ (the width of the Gaussians) are often required. In this case, the Gaussians become increasingly flat, and the interpolation matrix becomes ill-conditioned. This problem has been extensively studied in the literature (see the review [8] by Fornberg & Flyer and [19] for Tarwater’s description of this phenomenon in 1985). It has been quantified by Fornberg & Zuev [13], that the eigenvalues of the interpolation matrix are proportional to powers of the shape parameter, causing the notorious ill-conditioning in the flat limit regime $\varepsilon\to 0$.

A direct collocation solution of the interpolation problem, referred to as RBF-Direct in the literature, computes the expansion coefficients of the Gaussian interpolant by inverting the collocation matrix and then evaluating the expansion. Several algorithms have been proposed to stabilize this procedure in the flat limit regime, see [12, 11, 9, 7, 16, 5, 10]. A common idea of many of the stabilization algorithms—including the one proposed in this paper—is to evaluate the interpolant in a sequence of well-conditioned steps by a transformation to a different basis so that the ill-conditioning is isolated in a diagonal matrix that can be inverted analytically.

1.1. The new HermiteGF stabilization approach

In this paper, we propose a stabilizing expansion of isotropic Gaussian functions, later referred to as HermiteGF expansion, built on the exponential generating function of the classic Hermite polynomials. For certain classes of functions, anisotropic Gaussians yield improved accuracy as shown in [2]. To include these cases in our description, we use an anisotropic generating function recently obtained by Dietert, Keller, & Troppmann [4] as well as by Hagedorn [14] and generalize our HermiteGF expansion to the anisotropic case. We also propose and analyze a novel cut-off criterion, that contrary to the existing ones does not only account for the diagonal contributions of the stabilization but measures the impact of the full stabilization basis.

1.2. Previous stabilization approaches

The first stabilization method was the Contour–Padé approximation proposed by Fornberg & Wright for multiquadrics [12]. Later Fornberg & Piret [11] proposed the so-called RBF-QR method for stable interpolation with Gaussians on the sphere by expanding the Gaussians in spherical harmonics. The method has been extended to more general domains in one to three dimensions by Fornberg, Larsson, & Flyer [9]. This expansion is based on a combination of Chebyshev polynomials and spherical harmonics. This method will be referred to as Chebyshev-QR in this paper. The technique has also been used for the stable computation of RBF-generated finite differences by Larsson, Lehto, Heryudono, & Fornberg [16]. Fornberg, Lehto, & Powell [10] developed an alternative stabilization technique for the same problem. To treat complex domains, the Chebyshev-QR method has been combined with a partition-of-unity approach by Larsson, Shcherbakov, & Heryudono [17]. Fasshauer & McCourt [7] have developed another RBF-QR method, called Gauss-QR, that relies on a Mercer expansion of the Gaussian kernel. The basis transformation involves scaled Hermite polynomials. Compared to the Chebyshev-QR method by Fornberg et al. [9], the Gauss-QR algorithm extends to higher dimensions in a simpler way. On the other hand, the method introduces an additional parameter that needs to be hand-tuned. Our new basis is similar to the one in [7] with the difference that it can be extended to the interpolation with anisotropic Gaussians. Moreover, the generating function framework enables us to derive a new cut-off criterion that accounts for the full Hermite basis effect.

1.3. Organization of the paper

The paper is organized as follows: In section 2, we introduce our HermiteGF expansion of isotropic and anisotropic Gaussians and derive important properties of the HermiteGF basis. In section 3, we propose a new RBF-QR method based on the HermiteGF basis. A cut-off criterion for the new HermiteGF expansion is derived in section 4. Numerical results show the accuracy of our method in section 5 and finally, section 6 concludes the paper.

2. HermiteGF expansion

In this section, similarly to the earlier stabilization approaches [11, 9, 7], we propose an expansion of the anisotropic Gaussians in a “better” basis, that spans the same space, but avoids instabilities related to the flat limit. For the sake of simplicity, we first derive the expansion for Gaussians in 1D. We then extend our expansion to the case of multivariate anisotropic Gaussians.

2.1. Interpolation problem

Before introducing our expansion of the Gaussian basis, let us briefly define the interpolation problem.

Given a set $\{\phi_{k}(\mathbf{x})\}_{k=1}^{N}$ of expansion functions and the values $\boldsymbol{f}=\{f_{i}\}_{i=1}^{N}$ of the function $f:\mathbb{R}^{d}\to\mathbb{R}$ at collocation points $\{\mathbf{x}_{i}^{\mathrm{col}}\}_{i=1}^{N}$ we seek to find an interpolant of the form,

[TABLE]

such that it satisfies the $N$ collocation conditions,

[TABLE]

The direct approach is to find the coefficients $\boldsymbol{\alpha}=\{\alpha_{k}\}_{k=1}^{N}$ as a solution of the linear system,

[TABLE]

The matrix $\Phi^{\mathrm{col}}\in\mathbb{R}^{N\times N}$ is called collocation matrix. Then, after solving the linear system (2.2), the interpolant (2.1) can be evaluated at any point of the domain. In this paper, we consider Gaussian radial basis functions (isotropic Gaussians)

[TABLE]

with center points $X^{\mathrm{cen}}=\{\mathbf{x}_{k}\}_{k=1}^{N}$, shape parameter $\varepsilon>0$ and anisotropic Gaussians

[TABLE]

with invertible shape matrix $E\in\mathbb{R}^{d\times d}$.

2.2. HermiteGF expansion in 1D

Let $\{h_{\ell}\}_{\ell\geq 0}$ be the Hermite polynomials in the physicists’ version, that satisfy the recurrence relation,

[TABLE]

The following upper bound holds for the magnitude of the Hermite polynomials [1, Expr. 22.14.17]: There exists a positive constant $c\approx 1.086435$ such that

[TABLE]

for all $\ell\geq 0$ and all $x\in\mathbb{R}$. The combinatorial factors $\sqrt{2^{\ell}\ell!}$ grow very fast with $\ell$. Therefore, in order to avoid overflow for large $\ell$, we work with a scaled version of the Hermite polynomials, $\frac{1}{\sqrt{2^{\ell}\ell!}}h_{\ell}$. We introduce three parameters,

[TABLE]

that is, the usual shape parameter $\varepsilon>0$, a scaling parameter $\gamma>0$ for varying the evaluation domain of the Hermite polynomials to improve conditioning, and a truncation parameter $t\in(0,1)$ for controlling the truncation error of the stabilization expansion. We then define the following basis functions,

[TABLE]

that we refer to as HermiteGF functions.

Based on the generating function theory, we derive an infinite expansion of the one dimensional Gaussian RBFs in the new HermiteGF basis $\{H^{\gamma,\varepsilon,t}_{\ell}\}_{\ell\geq 0}$.

Proposition 2.1** (HermiteGF expansion (1D)).**

For all parameters $\varepsilon>0$, $\gamma>0$, $t\in(0,1)$, and for all shifts $x_{0}\in\mathbb{R}$, we have a pointwise expansion

[TABLE]

for all $k=1,\ldots,N$, where $\Delta_{k}=x_{k}-x_{0}$. The RBF interpolant $s(x)$ can then be pointwise computed as

[TABLE]

Proof.

The Hermite polynomial’s generating function is given by (see e.g. [1, Expr. 22.9.17]),

[TABLE]

Choosing $a=\frac{\varepsilon^{2}\Delta_{k}}{\gamma}$ and $b=\gamma(x-x_{0})$, we obtain

[TABLE]

Hence, we get

[TABLE]

which proves expansion (2.3). Using expansion (2.3) in the interpolant (2.1), we get the representation (2.4). ∎

Remark 2.1** (Basis centering).**

Hermite polynomials are symmetric with respect to the axis $x=0$. Due to its growth behavior, it is advantageous to have the basis centered around this point of symmetry, that is, to use the translation $x_{0}=\frac{B-A}{2}$, where $[A,B]$ is the interval of interest for evaluating the function $f$.

Remark 2.2** (The parameter $\gamma$).**

The parameter $\gamma>0$ in the basis $\{H^{\gamma,\varepsilon,t}_{\ell}\}_{\ell\geq 0}$ allows control over the evaluation domain of the Hermite polynomials. When choosing it, one has to consider two counteracting effects: For small values of $\gamma$, ill-conditioning appears since the values of the basis functions at the collocation points are too similar. On the other hand, Hermite polynomials take very large values on large domains which can lead to an overflow. An optimal balance depends on the particular function and the number of basis functions.

2.3. Multivariate HermiteGF expansion of anisotropic Gaussians

The HermiteGF expansion can be easily extended to higher dimensions for the case of isotropic Gaussians, using tensor products of 1D physicists’ Hermite polynomials,

[TABLE]

However, finding a stable interpolant for anisotropic Gaussian functions of the type

[TABLE]

is a more challenging task. A similar question was raised in [7, $\mathsection$ 8.5], however, without further investigation. McCourt & Fasshauer [18] considered anisotropic Gaussians with diagonal shape matrix $E$ using Mercer expansion theory, but this result has not been extended to the case of arbitrary $E$. Analogously to the 1D case, we define the multivariate version of our HermiteGF functions by

[TABLE]

where $G,E\in\mathbb{R}^{d\times d}$ are arbitrary invertible matrices.

Proposition 2.2** (HermiteGF expansion of anisotropic Gaussians).**

Let $\mathbf{q}\in\mathbb{R}^{d}$ and $E,G\in\mathbb{R}^{d\times d}$ invertible matrices. Consider the anisotropic Gaussian $\phi_{\mathbf{q}}(\mathbf{x})$. Then, for any shift $\mathbf{x}_{0}\in\mathbb{R}^{d}$ the following relation holds pointwise in $\mathbf{x}\in\mathbb{R}^{d}$:

[TABLE]

*where $\boldsymbol{\Delta}_{\mathbf{q}}=\mathbf{q}-\mathbf{x}_{0}$, $\tilde{G}=E^{T}EG^{-T}G^{-1}E^{T}E$. *

Proof.

We use the following anisotropic Hermite generating function (see [4, Lemma 5] or [14, Theorem 3.1] with $A=\mathrm{Id}_{d}$):

[TABLE]

We denote $\mathbf{b}=G^{T}(\mathbf{x}-\mathbf{x}_{0})$ and $\mathbf{a}=G^{-1}E^{T}E\boldsymbol{\Delta}_{\mathbf{q}}$. Then, using (2.7), we get

[TABLE]

We observe that

[TABLE]

Putting everything together we arrive to (2.6). ∎

This expansion provides a new powerful tool for dealing with anisotropic approximation. Note that the standard multidimensional isotropic Gaussian interpolation corresponds to the following matrix $E$:

[TABLE]

2.4. Mehler’s formula (nD)

In this section, we investigate the behavior of the HermiteGF basis functions with large index $\boldsymbol{\ell}$ approaching infinity. To be able to decide where to cut the HermiteGF expansion for numerical purposes, it is useful to quantify the size of the tail of the truncated expansion. For that, we take a look at the magnitude of the values of the HermiteGF basis. We consider the infinite-dimensional vector that contains the values of all basis functions $H_{\mathbb{\boldsymbol{\ell}}}^{G,E,t}(\mathbf{x})$ at a certain point $\mathbf{x}\in\mathbb{R}^{d}$. Its Euclidean norm can be computed analytically using a multivariate extension of Mehler’s formula for Hermite polynomials.

Theorem 2.1** (Bilinear generating function).**

The following relation holds for all $t\in(0,1)$ and all $\mathbf{x},\mathbf{y}\in\mathbb{R}^{d}$

[TABLE]

Proof.

We extend the 1D proof proposed by Watson in [20, p. 4] to multiple dimensions. Applying the inverse Fourier transform to the Fourier transform of a normal distribution with $\sigma=1/\sqrt{2}\mathrm{Id}_{d}$, we get

[TABLE]

Hence,

[TABLE]

where we used the Rodrigues formula for multivariate tensor-product Hermite polynomials (see [4, Expr. 11] with $M=\mathrm{Id}$). Then,

[TABLE]

where we used the Taylor series of the exponential function. Recall the following formula for a bivariate Gaussian Fourier integral (see [20, p. 4]):

[TABLE]

with

[TABLE]

Using (2.10) $d$ times for the integral (2.9) together with the algebraic identity

[TABLE]

yields (2.8).

∎

With the help of the multidimensional Mehler’s formula we can now compute the square of the Euclidean norm

[TABLE]

of the infinite vector containing the values of all basis functions $H_{\mathbb{\boldsymbol{\ell}}}^{G,E,t}(\mathbf{x})$ with $\boldsymbol{\ell}\in\mathbb{N}^{d}$ at some point $\mathbf{x}\in\mathbb{R}^{d}$. Indeed,

[TABLE]

We take a look if this value matches our numerical result. One can see in Figure 1 that for the square domain the limit value is matched for different values of $\varepsilon$, $t$.

3. Stabilization of the RBF interpolation

In this section, we derive a numerical stabilization algorithm for RBF interpolation based on the HermiteGF expansion. The main idea is to represent the RBF interpolant in the more stable basis $\{H^{G,E,t}_{\boldsymbol{\ell}}\}_{\boldsymbol{\ell}\geq 0}$. For appropriately chosen parameters $G$ and $t$, we expect it to be better conditioned. We now first derive the HermiteGF-QR algorithm in 1D and then generalize it to a multidimensional form.

3.1. HermiteGF-QR (1D)

In 1D we denote by

[TABLE]

the vector of Gaussians with center points $X^{\mathrm{cen}}$ evaluated at $x\in\mathbb{R}$ and write the stabilization expansion (2.3) as an infinite matrix-vector product

[TABLE]

where the vector

[TABLE]

contains all the elements of the polynomial basis $\{H^{\gamma,\varepsilon,t}_{\ell}\}_{\ell\geq 0}$ evaluated in the translated point $x-x_{0}$ and

[TABLE]

is an $\infty\times N$ matrix. The major part of the ill-conditioning is now confined in the matrix $B$. Since $B$ is independent of the point $x$ where the basis function is evaluated, both the evaluation and interpolation matrix can be expressed in the form (3.1) with the same matrix $B$.

We follow the RBF-QR approach and further split

[TABLE]

into a well-conditioned full $N\times\infty$ matrix $C$ and an infinite diagonal matrix $D$, where all harmful effects are confined in $D$. In the case of expansion (2.3), the following setup follows naturally from the Chebyshev-QR theory [9, $\mathsection$ 4.1.3],

[TABLE]

Here we also divide each coefficient by the radius of the domain $L$ containing the center points in order to avoid ill-conditioning in $C$ coming from taking high powers of $x_{k}$. That might be dangerous when the domain is too large, however, it still extends the range of domain diameters possible.

The goal is now to find a basis $\{\psi_{j}\}$ spanning the same space as $\{\phi_{k}\}$ but yielding a better conditioned collocation matrix. In particular, we need an invertible matrix $X$ such that $X^{-1}\Phi(x)^{T}$ is better conditioned. Let us perform a QR-decomposition on $C=QR$. Then we get,

[TABLE]

where $R_{1}$ and $D_{1}$ are $N\times N$ matrices containing the upper (left) block of the infinite matrices $R$ and $D$, respectively, while $R_{2}$ and $D_{2}$ assemble the remaining entries. Consider $X=QR_{1}D_{1}$. The new basis $\Psi(\mathbf{x})^{T}:=X^{-1}\Phi(x)^{T}$ can be formed as,

[TABLE]

To avoid under/overflow in the computation of $D_{1}^{-1}R_{1}^{-1}R_{2}D_{2}$, we form the two matrices $\tilde{R}=R_{1}^{-1}R_{2}$ and $\tilde{D}$ with elements

[TABLE]

and compute their Hadamard product, $\tilde{R}.*\tilde{D}$. Despite the harmful effects contained in $D$, the resulting term $D_{1}^{-1}R_{1}^{-1}R_{2}D_{2}=\tilde{R}.*\tilde{D}$ is then harmless.

3.2. Multivariate anisotropic HermiteGF-QR

In this section, we derive an analog of the HermiteGF-QR algorithm for the multivariate case. Since we are now dealing with matrices, in the general case, it is impossible to separate $E$ from $\mathbf{x}_{k}$ in $(E\mathbf{x}_{k})^{\boldsymbol{\ell}}$ as we did before. Therefore, the flow of the HermiteGF-QR does not apply directly. However, it is still possible to tackle some part of the ill-conditioning analytically. Consider the following splitting of the matrix $G^{-1}E^{T}E\in\mathbb{R}^{d\times d}$:

[TABLE]

where $\mathrm{Diag}$ is the $d\times d$ diagonal matrix containing the diagonal elements of $G^{-1}E^{T}E$ and $\mathrm{Rem}$ contains the remaining off-diagonal terms. Denote

[TABLE]

Then, it holds that

[TABLE]

where $\mathrm{Diag}^{-1}\mathrm{Rem}$ can be computed analytically. Denote $\mathbf{d_{\mathrm{vec}}}=\mathrm{diag}(\mathrm{Diag})$. Then, the generating function expansion (2.6) can be written as:

[TABLE]

As before, we can write the expansion above as the infinite matrix-vector product

[TABLE]

with

[TABLE]

As before, we write the transpose of the infinite matrix $B$ as a product $CD$ with

[TABLE]

The $d\times d$ matrix product $\mathrm{Diag}^{-1}\mathrm{Rem}$ contained in the vectors $\mathbf{v}_{k}$ can be computed analytically. The diagonal part of the matrix $G^{-1}E^{T}E$, that is now in the matrix $D$ can be handled in the exact same fashion as it has been done in 1D. In particular, we perform the QR-decomposition of the matrix $C=QR$, block decompose

[TABLE]

such that the entries related to the first $N$ stabilizing basis functions are contained in the $N\times N$ matrices $R_{1}$ and $D_{1}$, and consider the preconditioner $X=QR_{1}D_{1}$. Hence, analogously to the HermiteGF-QR case, the new basis can be formed as

[TABLE]

The action of $D_{1}^{-1}$ and $D_{2}$ can again be computed as the Hadamard product with

[TABLE]

Remark 3.1**.**

When the magnitude of $\mathbf{v}_{k}$ gets too large, the matrix $C$ can also become ill-conditioned. To avoid that, one can increase the magnitude of the elements of $G$. Alternatively, one can add a scaling in $C$ which should then be compensated in the matrix $D$, similarly to the scaling with the domain size $L$ in the 1D version.

3.2.1. Isotropic HermiteGF-QR

In the isotropic case, when $E=\varepsilon\mathrm{Id}_{d}$ and $G=\gamma\mathrm{Id}_{d}$, the expressions for the matrices $C$ and $D$ can be written in a simpler way. In particular, in this case, we have

[TABLE]

Hence, the diagonal vector simplifies to $\mathbf{d_{\mathrm{vec}}}=\gamma^{-1}\varepsilon^{2}(1,\ldots,1)$ and the elements of the matrices $C$ and $D$ take the form

[TABLE]

3.3. HermiteGF interpolant

In the new basis, we can write the equivalent formulation of the interpolant (2.1) as

[TABLE]

To compute the interpolant in the new formulation numerically, we have to cut the infinite expansion (2.6). We discuss the cut-off strategy in section 4.

3.4. Alternative splitting based on the Vandemonde matrix

Instead of using a QR-decomposition of the matrix $C\in\mathbb{R}^{N\times\infty}$ one could also split it as $C=\bar{E}W$, where $\bar{E}\in\mathbb{R}^{N\times N}$ is a diagonal matrix for the exponential part and $W\in\mathbb{R}^{N\times\infty}$ accounts for the polynomial contributions,

[TABLE]

We now decompose the original basis using this splitting,

[TABLE]

where $W_{1}\in\mathbb{R}^{N\times N}$ and $W_{2}\in\mathbb{R}^{N\times\infty}$. With the preconditioner $X_{V}=\bar{E}W_{1}D_{1}$, the new basis reads as

[TABLE]

One can see that

[TABLE]

where $C_{1}$ is the first $N\times N$ block of $C$. Therefore, in exact arithmetic $\Psi$ and $\Psi_{V}$ are the same, however, in floating-point arithmetic the values of the bases might differ. However, we will use this alternative splitting for the derivation and analysis of a suitable cut-off criterion for the HermiteGF basis in the next section.

4. Cut-off of the expansion

To make the RBF-QR methods usable for numerical computations, one has to cut the expansion (2.6) at a certain polynomial degree ${j_{\mathrm{max}}}\in\mathbb{N}$ which in 1D also corresponds to the number of stabilizing basis functions $M$. In the multivariate setting, the number of basis functions $M$ equals

[TABLE]

However, choosing an efficient cut-off degree ${j_{\mathrm{max}}}$ is not a trivial task. We first derive a criterion that is analogous to the state-of-the-art criteria for other RBF-QR methods ([9, Expr. 5.2], [7, Expr. 4.10]). However, especially in the HermiteGF case, it turns out to be inefficient, i.e. it overestimates the number ${j_{\mathrm{max}}}$. We then derive a new cut-off criterion based on the theoretical framework presented in the previous sections. This new criterion allows to directly control the approximation error of the stable basis which is more efficient while still being effective.

4.1. State-of-the-art criterion

Similarly to [9, Expr. 5.2], [7, Expr. 4.10] we take a look at the matrix $\tilde{D}$ that contains the effects of $D_{1}^{-1}$ and $D_{2}$. Recall that for the QR methods in section 3.1 the matrix $\tilde{D}$ is then multiplied element-wise with the matrix $\tilde{R}$. We stop once all elements of the new block of $\tilde{D}$ are below machine precision, i.e.

[TABLE]

The criterion (4.1) guarantees that all additional columns that could be added to $D_{2}$ would yield elements in $\tilde{D}$ that are below machine precision. We now take a look at the behavior of the elements of the matrix $D$. Here, for the Gauss-QR method we used $\alpha_{\mathrm{Gauss-QR}}=\gamma$. One can see in Figure 2 that for small values of $\varepsilon$ the behavior is very similar for all three methods. However, for large $\varepsilon$ the decay in $D$ is particularly bad in our formulation. This criterion also neglects the matrix $\tilde{R}$ and the effect of the polynomial vector $H^{\gamma,\varepsilon,t}(x-x_{0})$. In particular, we know from section 2.4 that the tail of the polynomial vector has some decay.

4.2. New HermiteGF cut-off criterion

In this section, we derive a more holistic criterion for the cut-off in the HermiteGF expansion. We use the Vandermonde formulation of the method since it provides an explicit expression for the elements of all matrices which simplifies the analytical study of the method. We cut the polynomial vector as

[TABLE]

with $\hat{H}^{G,E,t}\in\mathbb{R}^{1\times M}$, where the number of basis functions used is larger than the number of collocation points, that is, $M\geq N$. Analogously we cut the $N\times\infty$ Vandermonde matrix $W_{2}$ and the infinite diagonal matrix $D_{2}$,

[TABLE]

with $\hat{W}_{2}\in\mathbb{R}^{N\times(M-N)}$ and $\hat{D}_{2}\in\mathbb{R}^{(M-N)\times(M-N)}$. We note that the infinite matrix $W_{\infty}$ contains the columns of the full Vandermonde matrix $W$ from column $M+1$ onward, while the infinite matrix $D_{\infty}$ contains the diagonal entries of the full diagonal matrix $D$ starting from the entry $M+1$. We then rewrite the formulation of the method (see section 3.4) after the cut-off,

[TABLE]

We want to make sure that

[TABLE]

is small for all collocation points by choosing a sufficiently large but not too large truncation parameter ${j_{\mathrm{max}}}$. For estimating $\delta\Psi(\mathbf{x})$, we need the following lemma.

Lemma 4.1** (Exponential tail).**

Consider $\mathbf{y}\in\mathbb{R}^{d}$ with $y_{i}\geq 0$ for all $i=1,\ldots,d$ and ${j_{\mathrm{max}}}\in\mathbb{N}$. Then,

[TABLE]

where $\|\mathbf{y}\|_{1}=|y_{1}|+\ldots+|y_{d}|$.

Proof.

Consider the function

[TABLE]

We note that the estimated sum coincides with the remainder of the Taylor series for the function $f$ at the point $\mathbf{a}=0$. Then, according to the multivariate Taylor’s theorem, there exists $\boldsymbol{\xi}\in[0,\mathbf{y}]$ such that

[TABLE]

Noting that $\partial^{\boldsymbol{\ell}}f(\boldsymbol{\xi})=\exp(\mathbf{\boldsymbol{\xi}})\leq\exp(\|\mathbf{y}\|_{1})$ for $\mathbf{y}\geq 0$, we arrive at (4.2). ∎

Before proceeding to the estimation of the truncation error $\|\delta\Psi(\mathbf{x})\|_{2}$, we recall the definition of the vectors

[TABLE]

for $k=1,\ldots,N$. They will contribute to the upper bound of the following estimate. In the isotropic case, they have the particularly simple form

[TABLE]

Theorem 4.1** (Truncation estimate).**

For $k=1,\ldots,N$ we set

[TABLE]

where $W_{1}$ is the upper left $N\times N$ block of the infinite Vandermonde matrix $W=(\mathbf{v}_{k}^{\boldsymbol{\ell}})$. For ${j_{\mathrm{max}}}\in\mathbb{N}$ we denote

[TABLE]

Then, the truncation error $\delta\Psi$ satisfies for all $\mathbf{x}\in\mathbb{R}^{d}$,

[TABLE]

where $H^{G,E,t}_{\mathrm{lim}}(\mathbf{x}-\mathbf{x}_{0})$ can be evaluated via (2.11).

Proof.

We start by observing that

[TABLE]

Hence, it holds that

[TABLE]

and due to compatibility of Frobenius norm and the 2-norm

[TABLE]

We further consider the two norms on the right-hand side separately. We first take a look at the Frobenius norm. Recall that in the RBF-QR method we evaluate the effect of the impact of $D_{1}^{-1}\ldots D_{2}$ analytically. We can do the same here:

[TABLE]

where $.*$ denotes the Hadamard product and $\tilde{D}_{\infty}$ is constructed analogously to (3.3). We write the Frobenius norm as

[TABLE]

and estimate with the help of the Cauchy-Schwarz inequality,

[TABLE]

where $\boldsymbol{\ell}\in\mathbb{N}^{d}$ is the $\ell$-th multi-index corresponding to our basis enumeration. We used the explicit expression of $D_{\ell\ell}$ as defined in (3.2) and write

[TABLE]

We denote $\tilde{\mathbf{y}}_{i}=\begin{pmatrix}\frac{2}{t}(\mathbf{y}_{i})_{1}^{2}&\ldots&\frac{2}{t}(\mathbf{y}_{i})_{d}^{2}\end{pmatrix}$. Then, by lemma 4.1 we get

[TABLE]

Using the multinomial theorem and the fact that

[TABLE]

we arrive to estimate (4.3). ∎

The denominator $\mathbf{d}_{\mathrm{vec}}^{2\mathbf{k}}$ in expression (4.3) can take extremely small values that can lead to underflow. To avoid this, it can be combined with $\|\mathbf{y}_{i}\|_{2}^{2({j_{\mathrm{max}}}+1)}$. For this, we define the $d$-dimensional index $\mathbf{j}_{d}=\left(\frac{1}{d},\,...,\,\frac{1}{d}\right)$ and use the following transformation

[TABLE]

where $./$ denotes component-wise division. Pulling out $\mathrm{Diag}^{\mathbf{j}_{d}}$ the constant of the estimate (4.3) can be rewritten as

[TABLE]

Note, that in the isotropic case, one can simplify $\mathbf{d}_{\mathrm{vec}}^{-2\mathbf{k}+2({j_{\mathrm{max}}}+1)\mathbf{j}_{d}}=(\varepsilon^{4}/\gamma^{2})^{{j_{\mathrm{max}}}+1-|\mathbf{k}|}$. We are now ready to formulate our cut-off criterion.

Criterion 1**.**

We choose ${j_{\mathrm{max}}}$ for the HermiteGF-QR method such that

[TABLE]

where $\{\mathbf{x}_{k}\}_{k=1}^{N}$ are the collocation points.

Since we are looking at the relative error, the tolerance TOL need not be machine precision. The crucial difference to the state-of-the-art criterion (4.1) is that now the TOL directly controls the accuracy of the stable basis $\Psi$. Depending on the desired accuracy, the tolerance can be adjusted for the specific problem.

4.3. Automatic detection of $t$

One can use the criterion above also for determining the value of the parameter $t$. We scan the whole spectrum of the values of $t$ and detect the one that yields the minimum amount of basis functions

[TABLE]

Note that very small values of $t$ can cause cancellations and should be excluded (see LABEL:{sec:numerics_iso_N}). Even though this introduces additional computational cost in the determination of the suitable expansion, it could be profitable for the cases where the basis is used multiple times after having fixed the number ${j_{\mathrm{max}}}$ as e.g. in a time loop.

5. Numerical results

We have implemented the HermiteGF interpolation in MATLAB. The code is available for download at https://gitlab.mpcdf.mpg.de/clapp/hermiteGF. We compare the isotropic HermiteGF-based algorithm with the existing stabilization methods, the Chebyshev-QR method111Code downloaded from http://www.it.uu.se/research/scientific_computing/software/rbf_qr on September 10, 2018. and the Gauss-QR method222Code downloaded from http://math.iit.edu/~mccomic/gaussqr/ on September 5, 2018.. We evaluate the influence of different parameters, such as $\varepsilon$, $\gamma$, number of collocation points $N$ on the quality of the interpolation. For the Gauss-QR method, we take the free parameter $\alpha$ to be equal to our value of $\gamma$, i.e., $\alpha_{\mathrm{Gauss-QR}}=\gamma$. Since there are no stabilization methods available for fully anisotropic interpolation, we test the anisotropic HermiteGF-QR only against the direct algorithm to verify the correctness. To determine the cut-off degree in the HermiteGF method, we use the HermiteGF cut-off criterion with $\text{{TOL}}=10^{-6}$, unless stated otherwise. For this tolerance, the HermiteGF-QR method provides results that match the results from Chebyshev-QR and Gauss-QR. The parameter $t$ is detected automatically. In all tests we evaluate the interpolant at a set of evaluation points $\{\mathbf{z}_{k}\}_{k=1}^{N_{\rm ev}}$ and look at the average error of the form [7, $\mathsection$ 5.1, Expr. (5.2)]:

[TABLE]

5.1. 2D isotropic interpolation

In this section, we take a look at the two-dimensional interpolation with HermiteGF-QR. We take multiples of the identity for both $E$ and $G$. We look at a hyperbolic domain (see Figure 3) defined by the inequality

[TABLE]

with a boundary condition $x^{2}+y^{2}\leq 1$. The hyperbola of type (5.1) can then be parameterized as

[TABLE]

We run the tests for the following function ($f_{4}$ from [9, $\mathsection$ 6]):

[TABLE]

We investigate the behavior of the performance of the interpolation with respect to the parameters $\gamma$, $\varepsilon$, and number of functions $N$. We use $\gamma=3.5$, $\varepsilon=0.05$ and optimize $t$ from the set $\mathrm{tvec}=\texttt{linspace(0.1, 0.99, 10)}$, unless stated otherwise.

We sample the collocation points from Halton points that are clustered to the boundary and mapped to the hyperbolic domain. For all tests, we use $N_{\rm ev}=53^{2}$ evaluation points that are sampled similarly to the collocation points, but based on a uniform grid and without clustering. The nodes distribution is depicted in Figure 3. This domain and sampling strategy choice was inspired by [9, $\mathsection$ 6.1.2].

5.1.1. The number of nodes $N$

Let us first look at the behavior of the method for different numbers of nodes, $N$. We take the values of $N\in[100,410]$ of the form $\begin{pmatrix}P+2\\ 2\end{pmatrix}$ for some integer $P$, such that there are no same powers of $\varepsilon$ present in both $D_{1}$ and $D_{2}$. In Figure 4, we see that the error consistently decays for all the tested methods. Choosing the truncation parameter $t$ in the interval $t\in(0.1,1)$, the conditioning of the HermiteGF-QR method is slightly worse than for the other methods, since big powers of smaller values of $t$ yield cancellations. Indeed, limiting the range of $t$ to $t\in(0.3,1)$ brings the conditioning to the level of the other methods. Using all integers in the interval $[100,410]$ also provided consistent results for all three methods, however, the picture gets noisy. A snippet of that behavior can be seen in the zoomed regions in Figure 4. This can be related to the fact that for the values of $N$ of the form above the limit of the RBF interpolant in the flat limit $\varepsilon\rightarrow 0$ is a unique polynomial of degree $P$ [15, $\mathsection 4$ Theorem 4.1] whereas for other values the uniqueness is not guaranteed.

5.1.2. Sensitivity to $\gamma$

Let us take a look at the influence of the parameter $\gamma$ on the interpolation quality. We see in Figure 5 that for small $N$ the interpolation quality is not sensitive to the value of $\gamma$. However, for larger $N$ the parameter $\gamma$ has to be chosen with care. One can see in the Figure 5(b) that the conditioning is worse for small $\gamma$. However, one should be careful while increasing $\gamma$ since it also increases the evaluation domain of the Hermite polynomials, which take very large values on large domains which can lead to overflow. This effect becomes more pronounced as the degree of the Hermite polynomials increases. The optimal balance depends on the particular function and the number of basis functions.

5.1.3. Cut-off degree ${j_{\mathrm{max}}}$

Next, we look at the influence of the value of TOL on the quality of the interpolation. We compare the error only to the Gauss-QR method since the difference between the Chebyshev-QR and Gauss-QR results is down to machine precision. One can see in Figure 6 that for TOL $=10^{-6}$ the difference HermiteGF-QR and Gauss-QR is also down to machine precision. If we relax the tolerance to $10^{-2}$, the error is still small compared to the magnitude of the interpolation error, while having smaller ${j_{\mathrm{max}}}$, which yields an improved computational efficiency. Also, the figure shows a general trend that the expansion decays rather fast for small $\varepsilon$ while an increasing number of basis functions is needed for $\varepsilon$ close to 1.

5.2. 2D anisotropic interpolation

To test the performance of the HermiteGF expansion for anisotropic basis functions, we consider the function:

[TABLE]

As collocation points, we use Halton points clustered toward the boundaries. For the evaluation grid we use $53\times 53$ uniformly distributed points. We choose a non-diagonal matrix $G$. As for the shape matrix $E$, we check whether the off-diagonal elements influence the quality of the results. We fix $E$ and $G$ to be of the following form:

[TABLE]

We restrict $t$ to the interval $\mathrm{tvec}=\texttt{linspace(0.3, 1, 10)}$ in order to improve the stability of the computations. Let us take a look on how much the off-diagonal elements of the matrix $E$ influence the error. For our scan, we take 30 logarithmically distributed values of $\varepsilon\in[10^{-3},10^{0.1}]$. One can see from Figure 7 that properly chosen off-diagonal elements improve the quality of the interpolation. However, for larger $p$ slight instabilities occur which might be explained since $E$ becomes singular for $p\to 1$.

5.3. Multivariate interpolation

In this section, we consider an example of the usage of HermiteGF-QR in higher dimension. For all tests, we use the function

[TABLE]

where $|\mathbf{x}|=\sum_{i=1}^{d}x_{i}$. We use Halton collocation points and 1000 Halton points, excluding the ones used for collocation, for the evaluation grid. We fix $G=5\mathrm{Id}_{d}$, $\varepsilon=$ logspace(-3, 0.1, 30), and we again optimize the parameter $t$ over the set $\mathrm{tvec}=\texttt{linspace(0.3, 1, 10)}$. As before, we choose the tolerance $\texttt{TOL}=10^{-6}$. We first look whether the anisotropic HermiteGF-QR method converges to the results of the direct interpolation as $\varepsilon$ increases. We choose an arbitrary pattern for $E$ in order to verify that the stabilization works for a truly anisotropic interpolation,

[TABLE]

In Figure 8(a), we can see that HermiteGF-QR interpolation works stably even for very small values of $\varepsilon$ for different $N$. On the other hand, for larger values of $\varepsilon$ the result matches the direct anisotropic interpolation.

In order to validate the HermiteGF-QR method in higher dimensions against the existing methods, we compare the isotropic HermiteGF-QR method with the Gauss-QR method in 3-5D. For that, we fix the shape matrix $E$ and the number of interpolation points $N$ as

[TABLE]

where $d$ is the dimensionality. We choose the tolerance $\texttt{TOL}=10^{-2}$ since it is enough to meet the overall accuracy of the method. One can see in Figure 8(b) that the HermiteGF-QR method matches the reference Gauss-QR method in 3-5D.

6. Conclusion and Outlook

In this paper, we derive a new stabilization algorithm for Gaussian RBF interpolation in the flat limit ($\varepsilon\rightarrow 0$). The main idea of “isolating” the ill-conditioning in a special matrix is the same as in the previous approaches [11, 7, 9]. However, our new algorithmic framework draws from generating functions that naturally extend to the interpolation with anisotropic Gaussians. We introduce several parameters ($\varepsilon$/$E$, $\gamma$/$G$, $t$) for the new HermiteGF basis which have a distinct connotation: $\varepsilon$/$E$ is the original shape parameter of the Gaussian basis, $\gamma$/$G$ stands for the size of the evaluation domain of the Hermite polynomials, and the technical basis truncation parameter $t$ can be chosen automatically thanks to a novel analytically derived truncation criterion. The interpolation quality is not sensitive to the precise value of $\gamma$/$G$. The generic formulation of the method essentially provides an algorithm in any dimension, and we have reported results for up to five dimensions.

For the cut-off of the HermiteGF expansion, we derive a novel truncation criterion, generalizing Mehler’s theory for bilinear generating functions. In particular, we analytically estimate the truncation error $\delta\Psi$ in the stable HermiteGF basis $\Psi$. This allows adjusting the number of basis functions based on the desired accuracy in the basis $\Psi$.

The HermiteGF-QR algorithm has been implemented in MATLAB. For all isotropic test cases, the accuracy of the HermiteGF-QR method is consistent with the ones of established stabilization methods (Chebyshev-QR, Gauss-QR). For the anisotropic case, the results matches the RBF-Direct method, where the latter is applicable.

The stability of our algorithm could be further improved by employing special algorithms for the stable inversion of the Vandermonde type matrices. Based on the successful experience of the automatic detection of the truncation parameter $t$, it could be possible to develop an algorithm of choosing an optimal parameter matrix $G$ defining the effective evaluation domain. We believe that our anisotropic method could be also of interest in statistical data fitting [5, $\mathsection$ 17] or for continental size ice sheet simulations, see e.g. [3].

Acknowledgments

Fruitful discussions with Elisabeth Larsson (Uppsala University) and Michael McCourt (SigOpt) are gratefully acknowledged.