Radial basis and kernel interpolation

Fix an integer $d\ge 1$. An admitted domain $\mathrm{\Omega }\subset {\mathbb{R}}^d$ is a nonempty bounded Lipschitz domain (its boundary is locally the graph of a Lipschitz function, with $\mathrm{\Omega }$ on one side) with fixed $r_0>0$ and $0<\theta <\pi \mathrm{/}2$ such that for every $x\in \bar{\mathrm{\Omega }}$ there is an orthogonal map $Q_x$ for which $x+Q_x\{z\in {\mathbb{R}}^d:0<\mathrm{\mid }z\mathrm{\mid }<r_0,z_1>\mathrm{\mid }z\mathrm{\mid }\cos \theta \}\subset \mathrm{\Omega }\mathrm{.}$ An orthogonal map is a real linear map with $Q_x^{\mathsf{T}}{Q}_x=I$, equivalently one preserving the Euclidean inner product. The map $Q_x$ may depend on $x$; no simultaneous selection of the maps is part of this definition.

A real radial kernel on ${\mathbb{R}}^d$ has the form $K(x,y)=\varphi (\mathrm{\mid }x-y\mathrm{\mid })$ for a function $\varphi :[0,\mathrm{\infty })\to \mathbb{R}$. For a finite ordered list of distinct sites $X=(x_1,\dots ,x_N)$, its Gram matrix is $K_X=(K(x_i,x_j){)}_{i,j=1}^N$. The kernel is strictly positive definite if $c^{\mathsf{T}}{K}_{X}c>0$ for every such $X$ and every nonzero $c\in {\mathbb{R}}^N$; the nonstrict positive-semidefinite condition replaces $>$ by $\ge$. Only strict positivity is used for the ordinary interpolation solve below.

For an integer $m\ge 1$, let ${\mathrm{\Pi }}_{m-1}({\mathbb{R}}^d)$ be the real polynomials of total degree at most $m-1$, with an ordered basis $p_1,\dots ,p_J$, where $J=\dim {\mathrm{\Pi }}_{m-1}({\mathbb{R}}^d)$. The evaluation matrix on $X$ is $P_X=(p_j(x_i){)}_{i,j}$. The sites are $m$-unisolvent for this polynomial space when $P_X$ has linearly independent columns. A kernel is strictly conditionally positive definite of order $m$ when, for every finite distinct $X$ and this evaluation matrix, $c^{\mathsf{T}}{K}_{X}c>0$ for every nonzero $c\in {\mathbb{R}}^N$ satisfying $P_X^{\mathsf{T}}c=0$. Thus strict conditional positivity is a strict inequality on the moment-null subspace; it is stronger than the corresponding nonstrict inequality.

For data $y\in {\mathbb{R}}^N$, the strictly positive definite interpolation system is $K_{X}a=y$. For an order-$m$ strictly conditionally positive definite kernel and $m$-unisolvent sites, the polynomially augmented interpolation system is $\left[\begin{array}{cc}{K}_X& P_X\\ P_X^{\mathsf{T}}& 0\end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right]=\left[\begin{array}{c}y\\ 0\end{array}\right]\mathrm{.}$ Whenever the relevant system has a unique solution, its kernel-plus-polynomial function $x\mapsto {\sum }_i{a}_{i}K(x,x_i)+{\sum }_j{b}_j{p}_j(x)$ is called the interpolant. Invertibility for the stated strict hypotheses is supplied by Strict positive and conditionally positive augmented interpolation systems are invertible ↗.

For a finite nonempty site set $X\subset \mathrm{\Omega }$, the fill distance is $h_X={\sup }_{x\in \mathrm{\Omega }}{\min }_{x_i\in X}\mathrm{\mid }x-x_i\mathrm{\mid }$. If $N\ge 2$, the separation radius and mesh ratio are $q_X=\frac{1}{2}{\min }_{i\mathrm{\ne }j}\mathrm{\mid }{x}_i-x_j\mathrm{\mid }$ and ${\rho }_X=h_X\mathrm{/}{q}_X$, respectively. These are exact geometric quantities, independent of any floating-point solver.

Let $X=(x_1,\dots ,x_N)$ be a finite ordered list of distinct sites in ${\mathbb{R}}^d$. If $K$ is strictly positive definite, then $K_X$ is invertible and $K_{X}a=y$ has a unique solution for every $y\in {\mathbb{R}}^N$. If $K$ is strictly conditionally positive definite of order $m$ and $X$ is $m$-unisolvent for ${\mathrm{\Pi }}_{m-1}({\mathbb{R}}^d)$, then the augmented matrix in Radial kernels, polynomial augmentation, and site geometry is invertible and its system has a unique solution for every right-hand side.

Fix the dimension $d\ge 1$ of an admitted domain, a smoothness parameter $\nu >0$, and an inverse length $a>0$. Assume Countable Choice for the Lebesgue-integration suppliers used to justify this formula; this definition does not assume full Choice. Put $s=\nu +d\mathrm{/}2$ and $A=\frac{\mathrm{\Gamma }(s)(4\pi {)}^{d\mathrm{/}2}{a}^{2\nu }}{\mathrm{\Gamma }(\nu )},w(\xi )=A(a^2+4{\pi }^2\mathrm{\mid }\xi {\mathrm{\mid }}^2{)}^{-s}\mathrm{.}$ All integrals are Lebesgue integrals with respect to ${\lambda }_d$ (Lebesgue measurable sets, the family $\mathcal{L}({\mathbb{R}}^n)$, and the restricted set function ${\lambda }_n$), and $\mathrm{\Gamma }$ is the Euler integral of The real Gamma function by Euler's integral. With the negative-sign Fourier transform convention $\widehat{U}(\xi )=\int U(x)e^{-2\pi ix\cdot \xi }dx$, define the normalized stationary kernel by $k(z)={\int }_{{\mathbb{R}}^d}w(\xi )e^{2\pi iz\cdot \xi }d\xi ,K(x,y)=k(x-y)\mathrm{.}$ The corresponding length parameter in the usual Matérn notation is $\mathrm{\ell }=\sqrt{2\nu }\mathrm{/}a$. For $b>0$, define the modified Bessel function in its integral normalization by $K_{\nu }(b)=\frac{1}{2}{\left(\frac{b}{2}\right)}^{\nu }{\int }_0^{\mathrm{\infty }}{u}^{-\nu -1}{e}^{-u-b^2\mathrm{/}(4u)}du\mathrm{.}$ The classical Bessel expression for $k(z)$ when $z\mathrm{\ne }0$, its value $k(0)=1$, absolute convergence, continuity, and radial symmetry are established by The Matérn spectral integral, normalization, and radial formula ↗. This definition makes no finite-precision stability claim.

Assume the Axiom of Choice when using the well-definedness results below; it is inherited conservatively from the published $L^2$ completeness supplier. Let $d\ge 1$ and $s>d\mathrm{/}2$. Write $L_s^2({\mathbb{R}}^d)$ for the equivalence classes modulo Lebesgue-almost-everywhere equality of measurable complex functions $g$ satisfying $\mathrm{\parallel }g{\mathrm{\parallel }}_{L_s^2}^2={\int }_{{\mathbb{R}}^d}(1+\mathrm{\mid }\xi {\mathrm{\mid }}^2{)}^s\mathrm{\mid }g(\xi ){\mathrm{\mid }}^{2}d\xi <\mathrm{\infty }\mathrm{.}$ For such $g$, the inverse Fourier representative is the function $U_g(x)={\int }_{{\mathbb{R}}^d}g(\xi )e^{2\pi ix\cdot \xi }d\xi \mathrm{.}$ The whole-space Bessel-potential space $H^s({\mathbb{R}}^d)$ consists of these representatives, with norm $\mathrm{\parallel }{U}_g{\mathrm{\parallel }}_{H^s({\mathbb{R}}^d)}=\mathrm{\parallel }g{\mathrm{\parallel }}_{L_s^2}$. For the admitted domain $\mathrm{\Omega }$ of Radial kernels, polynomial augmentation, and site geometry, define $H^s(\mathrm{\Omega })$ to be the restrictions $U{\mathrm{\mid }}_{\mathrm{\Omega }}$ of whole-space members, with the quotient norm $\mathrm{\parallel }u{\mathrm{\parallel }}_{H^s(\mathrm{\Omega })}=\inf \{\mathrm{\parallel }U{\mathrm{\parallel }}_{H^s({\mathbb{R}}^d)}:U{\mathrm{\mid }}_{\mathrm{\Omega }}=u\}\mathrm{.}$ The whole-space inner product is $⟨U_g,U_h{⟩}_{H^s}=\int (1+\mathrm{\mid }\xi {\mathrm{\mid }}^2{)}^{s}g(\xi )\overline{h(\xi )}d\xi$, linear in the first variable as in Real and complex inner-product spaces and their induced length. The inverse Fourier map is injective, whole-space representatives are continuous, the whole-space model is Hilbert, and Schwartz functions are dense in that norm; these facts are supplied by The weighted Fourier model is an injective Hilbert representation ↗ and Schwartz functions are dense in the Fourier Bessel-potential norm ↗. The restriction quotient is a norm by bounded point evaluation; its Hilbert realization as the restricted native space is proved in The restricted Matérn native space is the Sobolev quotient space ↗. No extension operator or identification with an intrinsic fractional Sobolev norm on $\mathrm{\Omega }$ is part of this definition.

Let $\mathrm{\Omega }\subset {\mathbb{R}}^d$ be an admitted domain of Radial kernels, polynomial augmentation, and site geometry, let $X\subset \mathrm{\Omega }$ be finite and nonempty, and put $h=h_X$. Fix an integer $m\ge 1$ and let $J=\dim {\mathrm{\Pi }}_{m-1}({\mathbb{R}}^d)$. There are constants $R,B,h_0>0$, depending only on $d,m,r_0,\theta$, such that if $0<h\le h_0$, then for every $x\in \bar{\mathrm{\Omega }}$ there are distinct nodes $y_1,\dots ,y_J\in X$ and real weights ${\lambda }_1,\dots ,{\lambda }_J$ with $\mathrm{\mid }{y}_j-x\mathrm{\mid }\le Rh,{\sum }_{j=1}^J\mathrm{\mid }{\lambda }_j\mathrm{\mid }\le B,p(x)={\sum }_{j=1}^J{\lambda }_{j}p(y_j)(p\in {\mathrm{\Pi }}_{m-1}({\mathbb{R}}^d))\mathrm{.}$ No lower bound on site separation is assumed.

Assume Countable Choice, and fix $d\ge 1$, $\nu >0$, $a>0$, $s=\nu +d\mathrm{/}2$, and $A,w,k$ as in The normalized Matérn kernel in the 2π Fourier convention. Then $w>0$, $w\in L^1({\mathbb{R}}^d)$, and ${\int }_{{\mathbb{R}}^d}w(\xi )d\xi =1.$ The integral defining $k(z)$ is absolutely convergent for every $z$; $k$ is real, even, radial, and continuous, and $k(0)=1$. For every $z\in {\mathbb{R}}^d$, $k(z)=\frac{a^{2\nu }}{\mathrm{\Gamma }(\nu )}{\int }_0^{\mathrm{\infty }}{t}^{\nu -1}\exp \text{\hspace{0.17em}⁣}(-a^{2}t-\frac{\mathrm{\mid }z{\mathrm{\mid }}^2}{4t})dt\mathrm{.}$ For $z\mathrm{\ne }0$ this is equivalently $k(z)=\frac{2^{1-\nu }}{\mathrm{\Gamma }(\nu )}(a\mathrm{\mid }z\mathrm{\mid }{)}^{\nu }{K}_{\nu }(a\mathrm{\mid }z\mathrm{\mid }),$ where $K_{\nu }$ has the integral normalization in The normalized Matérn kernel in the 2π Fourier convention. Under the stated Fourier convention, the $d$-dimensional Gaussian identity is, for $t>0$, ${\int }_{{\mathbb{R}}^d}{e}^{-4{\pi }^{2}t\mathrm{\mid }\xi {\mathrm{\mid }}^2}{e}^{2\pi iz\cdot \xi }d\xi =(4\pi t{)}^{-d\mathrm{/}2}{e}^{-\mathrm{\mid }z{\mathrm{\mid }}^2\mathrm{/}(4t)}\mathrm{.}$

The exact augmented RBF finite-difference weight procedure has the following inputs: a dimension $d\ge 1$; an integer polynomial degree $r\ge 0$; a finite ordered stencil $Y=(y_1,\dots ,y_N)\subset {\mathbb{R}}^d$ of distinct sites that is $(r+1)$-unisolvent for ${\mathrm{\Pi }}_r({\mathbb{R}}^d)$; a real kernel $K$ strictly conditionally positive definite of order $r+1$; an ordered basis $p_1,\dots ,p_J$ of ${\mathrm{\Pi }}_r({\mathbb{R}}^d)$; a center $x\in {\mathbb{R}}^d$; a homogeneous constant-coefficient differential functional $L={\sum }_{\mathrm{\mid }\alpha \mathrm{\mid }=t}{c}_{\alpha }{D}^{\alpha }$ with $0\le t\le r$; and a neighborhood of $x$ on which every section $z\mapsto K(z,y_i)$ has continuous derivatives through order $t$. Here $\alpha \in {\mathbb{N}}_0^d$, $\mathrm{\mid }\alpha \mathrm{\mid }={\sum }_j{\alpha }_j$, and $D^{\alpha }$ is the corresponding partial derivative. All arithmetic in this procedure is exact real arithmetic.

The input is supplied with total, correct oracles for admissibility and exact evaluation. The admissibility oracle returns whether all the stated input conditions hold, including the polynomial-basis condition, global strict conditional positive definiteness, and continuous section differentiability on the supplied neighborhood. On an affirmative answer, the evaluation oracles return the exact real values $K(y_i,y_j)$, $p_j(y_i)$, $(L_{z}K(z,y_i)){\mathrm{\mid }}_{z=x}$, and $(L{p}_j)(x)$ for the finitely many requested indices. Oracle calls are primitive operations in this model; neither global kernel admissibility nor differentiation is claimed to be decidable by real arithmetic alone.

Query the admissibility oracle. If it returns false, return inadmissible-input. Otherwise use the evaluation oracles to form $P_{ij}=p_j(y_i)$, $b_i=(L_{z}K(z,y_i)){\mathrm{\mid }}_{z=x}$, ${\mathrm{\ell }}_j=(L{p}_j)(x)$, and $M=\left[\begin{array}{cc}{K}_Y& P\\ P^{\mathsf{T}}& 0\end{array}\right],g=\left[\begin{array}{c}b\\ ell\end{array}\right]\mathrm{.}$ Solve $Mz=g$ by exact Gaussian elimination with deterministic nonzero-pivot row selection: set $A=M$, $c=g$, and $n=N+J$. For $k=1,\dots ,n$, choose the least $i\ge k$ with $A_{ik}\mathrm{\ne }0$; if none exists, return singular-system. Swap rows $i$ and $k$ in both $A$ and $c$. For each $j=k+1,\dots ,n$ in increasing order, set $\mu =A_{jk}\mathrm{/}{A}_{kk}$, replace row $j$ by row $j-\mu {\mathrm{r}\mathrm{o}\mathrm{w}}_k$ in $A$ and $c$, and set $A_{jk}=0$. After all pivots, set $z_n=c_n\mathrm{/}{A}_{nn}$; then, for $k=n-1,n-2,\dots ,1$ in decreasing order, set $z_k=(c_k-{\sum }_{j=k+1}^n{A}_{kj}{z}_j)\mathrm{/}{A}_{kk}$. Return success with the first $N$ coordinates $w$ of $z$ and the remaining $J$ polynomial multipliers. The stated strict hypotheses imply that singular-system is unreachable, as proved by Strict positive and conditionally positive augmented interpolation systems are invertible.

The returned $w$ is an exact mathematical weight vector. The nonsingularity lemma Strict positive and conditionally positive augmented interpolation systems are invertible makes singular-system unreachable for every input that passes the stated strict hypotheses. This definition makes no floating-point stability or solver-error claim.

Assume the Axiom of Choice as a conservative cover of the published $L^2$ completeness dependency, which declares it in its transitive dependencies although its statement requires only Countable Choice. For every $s>d\mathrm{/}2$, each $g\in L_s^2({\mathbb{R}}^d)$ lies in $L^1({\mathbb{R}}^d)$, its inverse Fourier representative $U_g$ is bounded and continuous, the map $g\mapsto U_g$ is injective, and the whole-space Fourier Bessel-potential space of Fourier Bessel-potential and restriction-quotient Sobolev spaces is a complex Hilbert space. In particular, for a constant $C_{d,s}<\mathrm{\infty }$, $\mathrm{\mid }{U}_g(x)\mathrm{\mid }\le C_{d,s}\mathrm{\parallel }g{\mathrm{\parallel }}_{L_s^2}(x\in {\mathbb{R}}^d)\mathrm{.}$

Assume Countable Choice. Fix $d\ge 1$, $\nu >0$, $a>0$, $s=\nu +d\mathrm{/}2$, and the normalized Matérn kernel with spectral density $w(\xi )=A(a^2+4{\pi }^2\mathrm{\mid }\xi {\mathrm{\mid }}^2{)}^{-s}$ from The normalized Matérn kernel in the 2π Fourier convention. For a list of $N\ge 2$ distinct sites $X=(x_1,\dots ,x_N)\subset {\mathbb{R}}^d$, put $q=\frac{1}{2}{\min }_{i\mathrm{\ne }j}\mathrm{\mid }{x}_i-x_j\mathrm{\mid }$. There is $c=c(d,\nu ,a)>0$ such that if $0<q\le 1$, then ${\lambda }_{\min }(K_X)\ge c{q}^{2\nu },{\lambda }_{\max }(K_X)\le N,{\kappa }_2(K_X)\le c^{-1}N{q}^{-2\nu }\mathrm{.}$ The eigenvalues and condition number are those of the exact real symmetric positive-definite matrix. For $N=1$, $K_X=[1]$ and ${\kappa }_2(K_X)=1$.