On thin plate spline interpolation

TL;DR

This paper provides a PDE-based proof for improved error estimates in thin plate spline interpolation and demonstrates the effectiveness of ${ m f H}$-matrix techniques for large-scale problems.

Contribution

It introduces a new PDE-based proof for error estimate improvements and applies ${ m f H}$-matrix methods to efficiently solve large interpolation problems.

Findings

Error estimates can be improved by $h^{1/2}$

${ m f H}$-matrix techniques are effective for large problems

The proof simplifies understanding of error bounds in spline interpolation

Abstract

We present a simple, PDE-based proof of the result [M. Johnson, 2001] that the error estimates of [J. Duchon, 1978] for thin plate spline interpolation can be improved by $h^{1/2}$. We illustrate that ${\mathcal H}$-matrix techniques can successfully be employed to solve very large thin plate spline interpolation problems

Equations

$$\displaystyle\mbox{ minimize }\qquad|v|_{H^{m}({\mathbb{R}}^{d})}\qquad\mbox{ under the constraint $v(x_{i})=f(x_{i})$, $i=1,\ldots,N$.}$$

$$⟨v,w{⟩}_m:=\sum _{|\alpha |=m}\frac{m!}{\alpha !}{\int }_{{\mathbb{R}}^d}{D}^{\alpha }v{D}^{\alpha }wdx.$$

$$If(x)=\sum _{i=1}^N{c}_i{\varphi }_m(\Vert x-x_i{\Vert }_2)+{\pi }_{m-1}(x),\sum _{i=1}^N{c}_{i}q(x_i)=0\forall q\in {\mathbb{P}}_{m-1},$$

$${\varphi }_m(r)=\{\begin{array}{ll}{r}^{2m-d}\log r& d\text{mbox}even\\ r^{2m-d}& d\text{mbox}odd.\end{array}$$

$$h(X_N):=\underset{x\in \Omega }{\sup }\underset{i=1,\dots ,N}{\inf }\Vert x-x_i{\Vert }_2.$$

$$\sum _{|\alpha |=k}\Vert D^{\alpha }(f-If){\Vert }_{L^p(\Omega )}\le C_1{h}^{m-k-d/2+d/p}|E^{\Omega }f-If{|}_{H^m({\mathbb{R}}^d)}\le C_2{h}^{m-k-d/2+d/p}|f{|}_{H^m(\Omega )};$$

$$\sum _{|\alpha |=k}\Vert D^{\alpha }(f-If){\Vert }_{L^p(\Omega )}\le C_2{h}^{m-k-d/2+d/p}|f-If{|}_{H^m(\Omega )}.$$

$$|E^{\Omega }f-If{|}_{H^m({\mathbb{R}}^d)}$$

$$\text{mbox}{\mathrm{BL}}^m({\mathbb{R}}^d):=\{u\in {\mathcal{D}}^{\prime }|{\nabla }^{m}u\in L^2({\mathbb{R}}^d)\}.$$

$$E^{\Omega }U=\arg \min \{|u{|}_{H^m({\mathbb{R}}^d)}|u\in \text{mbox}{\mathrm{BL}}^m({\mathbb{R}}^d),u{|}_{\Omega }=U\}.$$

$$⟨E^{\Omega }U,v{⟩}_m=0\forall v\in \{v\in \text{mbox}{\mathrm{BL}}^m({\mathbb{R}}^d)|v{|}_{\Omega }=0\}.$$

$$\Delta^{m}E^{\Omega}U=0\qquad\mbox{ in ${\Omega^{c}}$}.$$

$$B_{\Omega }(u,v):=\sum _{|\alpha |=m}\frac{m!}{\alpha !}{\int }_{\Omega }{D}^{\alpha }u{D}^{\alpha }v,B_{{\Omega }^c}(u,v):=\sum _{|\alpha |=m}\frac{m!}{\alpha !}{\int }_{{\Omega }^c}{D}^{\alpha }u{D}^{\alpha }v.$$

$$\Vert u{\Vert }_{L^2(\partial \Omega )}\le C\Vert u{\Vert }_{B_{2,1}^{1/2}(\Omega )},\Vert {\nabla }^{k}u{\Vert }_{L^2(\partial \Omega )}\le C\Vert u{\Vert }_{B_{2,1}^{k+1/2}(\Omega )}.$$

$$|l(f)|$$

$$|l(f)|\le C{\epsilon }^{{\theta }_1}\Vert f{\Vert }_{X_{{\theta }_1}}\forall f\in X_{{\theta }_1}.$$

$$\Vert f-\tilde{f}{\Vert }_{X_0}+\epsilon \Vert \tilde{f}{\Vert }_{X_1}\le 2K(\epsilon ,f).$$

$$|l(f)|$$

$$|l(f)|\le C{\epsilon }^{\theta }\Vert f{\Vert }_{X_{\theta }}.$$

$$|l(f)|$$

$${\theta }_i={\theta }_1(1-s_i)+s_i.$$

$${\epsilon }^{{\theta }_i}\Vert \tilde{f}{\Vert }_{X_{{\theta }_i}}$$

$$=C\left[{\epsilon }^{{\theta }_1}\Vert \tilde{f}{\Vert }_{X_{{\theta }_1}}+\epsilon \Vert \tilde{f}{\Vert }_{X_1}\right].$$

$$|l(f)|\le C\left[K(\epsilon ,f)+{\epsilon }^{{\theta }_1}\Vert \tilde{f}{\Vert }_{X_{{\theta }_1}}\right].$$

$$\Delta^{m}u=g\quad\mbox{ in $\Omega$},\qquad u=\partial_{n}u=\cdots\partial_{n}^{m-1}u=0\quad\mbox{ on $\partial\Omega$},$$

$$\Vert u{\Vert }_{H^{m+k}(\Omega )}\le C_{\Omega ,m,k}\Vert g{\Vert }_{H^{-m+k}(\Omega )}.$$

$$|E^{\Omega }f{|}_{H^m({\mathbb{R}}^d)}\le C_{\Omega }\Vert f{\Vert }_{H^m(\Omega )}.$$

$$\Vert {\gamma }_0^c({\nabla }^{m+j}{E}^{\Omega }f){\Vert }_{L^2(\partial \Omega )}\le C_{j,\Omega }\Vert f{\Vert }_{B_{2,1}^{m+j+1/2}(\Omega )},$$

$$\Delta^{m}u=-\Delta^{m}Ef\in L^{2}(\widetilde{\Omega})\quad\mbox{ in $\widetilde{\Omega}$},\qquad u=\partial_{n}u=\cdots=\partial_{n}^{m-1}u=0\qquad\mbox{ on $\partial\widetilde{\Omega}$}.$$

$$\Vert E^{\Omega }f{\Vert }_{H^m(\tilde{\Omega })}\le C\Vert f{\Vert }_{H^m(\Omega )},\Vert E^{\Omega }f{\Vert }_{H^{2m}(\tilde{\Omega })}\le C\Vert f{\Vert }_{H^{2m}(\Omega )}.$$

Full text

1 Introduction and Main Results

Interpolation with so-called thin plate splines (also known as surface splines, $D^{m}$-splines, or polyharmonic splines) is a classical topic in spline theory. It is concerned with the following interpolation problem (1): Given a (sufficiently smooth) function $f$ and points $x_{i}\in{\mathbb{R}}^{d}$, $i=1,\ldots,N$, find the minimizer $If$ of the problem

[TABLE]

Here, the seminorm $|v|_{H^{m}({\mathbb{R}}^{d})}$ is induced by the bilinear form

[TABLE]

For $m>d/2$ and under very mild conditions on the point distribution, a unique minimizer $If$ exists. The name “thin plate splines” originates from the fact in the simplest case $m=d=2$, $If$ can be represented in terms of translates of the fundamental solution of the biharmonic equation. For general $m$ the interpolant $If$ can be expressed in terms fundamental solutions of $\Delta^{m}$: There are constants $c_{i}\in{\mathbb{R}}$, $i=1,\ldots,N$, and a polynomial $\pi\in{\mathbb{P}}_{m-1}$ of degree $m-1$ such that (with the Euclidean norm $\|\cdot\|_{2}$ on ${\mathbb{R}}^{d}$)

[TABLE]

where $\phi_{m}$ is given explicitly by

[TABLE]

The representation (3) allows one to reformulate (1) as the problem of finding the coefficients $c_{i}$ and the polynomial $\pi_{m-1}$ so that the (constrained) interpolation problem (3) is solved. The classical error analysis for (1) is formulated in terms fill-distance: For a bounded domain $\Omega\subset{\mathbb{R}}^{d}$ and points $X_{N}=\{x_{i}\,|\,i=1,\dots,N\}\subset\Omega$, the fill distance $h(X_{N})$ is given by

[TABLE]

Starting with the seminal papers by J. Duchon [12, 11] the error $f-If$ on $\Omega$ is controlled in terms of $h$ and the regularity properties of $f$ (on $\Omega$):

Proposition 1.1** ([11, Prop. 3])**

Let $\Omega\subset{\mathbb{R}}^{d}$ be a bounded Lipschitz domain. Let $m>d/2$, $k\in{\mathbb{N}}$, $p\in[2,\infty]$ be such that $H^{m}(\Omega)\subset W^{k,p}(\Omega)$. Then, there are constants $h_{0}$, $C_{1}$, $C_{2}>0$ depending only on $\Omega$, $m$, $d$ such that for any collection $X_{N}=\{x_{1},\ldots,x_{N}\}\subset\Omega$ with fill distance $h:=h(X_{N})\leq h_{0}$

[TABLE]

here, $E^{\Omega}f$ denotes the minimum norm extension of $f$ defined in (8).

In Proposition 1.1 and throughout the present note, we will use the standard notation for Sobolev spaces $W^{s,p}$ and Besov spaces $B^{s}_{2,q}$; we refer to [26] for their definition. Interpolation space will always be understood by the so-called “real method” (also known as “$K$-method”) as described, e.g., in [26, 27]. We will use extensively that the scales of Sobolev and Besov spaces are interpolation spaces. We will also use the notation $|\nabla^{j}f|^{2}=\sum_{|\alpha|=j}\frac{j!}{\alpha!}|D^{\alpha}f|^{2}$.

It is worth noting that the interpolation operator $I$ is a projection so that $I(f-If)=0$. Proposition 1.1 applied to the function $f-If$ therefore yields

Corollary 1.2

Under the assumptions of Proposition 1.1 there holds

[TABLE]

A natural question in connection with Proposition 1.1 is whether the convergence rate can be improved by requiring additional regularity of $f$. It turns out that boundary effects limits this. We mention that a doubling of the convergence rate is possible by imposing certain homogeneous boundary conditions on high order derivatives as shown in [22] and, more abstractly, in [24]. If this highly fortuitous setting is not given, then only a small further gain is possible as shown by M. Johnson, [17, 18]. For example, he showed that a gain of $h^{1/2}$ is possible if $f\in B^{m+1/2}_{2,1}(\Omega)$ and $\partial\Omega$ is sufficiently smooth. The purpose the present note is to give a short and simple proof of this result using different tools, namely, those from elliptic PDE theory. The techniques also open the door to reducing the smoothness assumptions on $\partial\Omega$ in [17, 18] to Lipschitz continuity as discussed in more detail in Remark 2.8. Our main result therefore is a simpler proof of:

Proposition 1.3** ([17])**

Let $\Omega\subset{\mathbb{R}}^{d}$ be a bounded Lipschitz domain with sufficiently smooth boundary. Then there are constants $h_{0}$, $C_{1}$, $C_{\delta}>0$ that depend solely on $\Omega$, $m$, $d$, and $\delta$ such that for any collection $X=\{x_{1},\ldots,x_{N}\}\subset\Omega$ with fill distance $h:=h(X_{N})\leq h_{0}$ there holds

[TABLE]

In particular, therefore, the estimates of [11, Prop. 3] (i.e., Prop. 1.1) can be improved by $h^{1/2}$ for $f\in B^{m+1/2}_{2,1}(\Omega)$ and by $h^{\delta}$ for $f\in H^{m+\delta}(\Omega)$.

Remark 1.4

*A common route to error estimates for $f-If$ is via the so-called “power function” $P(x)$. Indeed, classical pointwise estimates take the form $|f(x)-If(x)|\leq P(x)|E^{\Omega}f-If|_{H^{m}({\mathbb{R}}^{d})}$ (cf., e.g., [8, Prop. 5.3], [29, Thm. 11.4]) and $P$ is subsequently estimated in terms of the fill distance $h$. Thus, Proposition 1.3 allows for improving estimates in this setting.  *

We close this section by referring the reader to the monographs [29, 8] as well as [16] for further details on the approximation properties of radial basis functions, in particular, thin plate splines.

2 Proof of Proposition 1.3

2.1 Tools

The precise formulation of the minimization problem (1) is based on the classical Beppo-Levi space ${\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$, which is defined as

[TABLE]

We refer to [10] and [29, Sec. 10.5] for more properties of the space ${\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$; in particular, $C^{\infty}_{0}({\mathbb{R}}^{d})$ is dense in ${\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$ (see [29, Thm. 10.40] for the precise notion). We also need the minimum norm extension $E^{\Omega}:H^{m}(\Omega)\rightarrow{\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$ given by

[TABLE]

The minimization property in (8) implies the orthogonality

[TABLE]

The connection with elliptic PDE theory arises from the fact that $E^{\Omega}U$ satisfies an elliptic PDE in ${\Omega^{c}}:={\mathbb{R}}^{d}\setminus\overline{\Omega}$:

[TABLE]

It will be convenient to decompose $B(u,v):=\langle u,v\rangle_{m}=\sum_{|\alpha|=m}\frac{m!}{\alpha!}\int_{{\mathbb{R}}^{d}}D^{\alpha}uD^{\alpha}v$ as $B(u,v)=B_{\Omega}(u,v)+B_{{\Omega^{c}}}(u,v)$, where

[TABLE]

The trace mapping is continuous $H^{1/2+\varepsilon}(\Omega)\rightarrow H^{\varepsilon}(\partial\Omega)$ for $\varepsilon\in(0,1/2]$; however, the limiting case $\varepsilon=0$ is not true; it is true if the Sobolev space $H^{1/2}(\Omega)$ is replaced with the slightly smaller Besov space $B^{1/2}_{2,1}(\Omega)$:

Lemma 2.1** (Trace theorem)**

Let $\Omega\subset{\mathbb{R}}^{d}$ be a Lipschitz domain, $k\in{\mathbb{N}}_{0}$. Then there exists $C>0$ such that the multiplicative estimate $\|u\|^{2}_{L^{2}(\partial\Omega)}\leq C\|u\|_{L^{2}(\Omega)}\|u\|_{H^{1}(\Omega)}$ holds as well as

[TABLE]

Proof. The case $k\geq 1$ in (11) follows immediately from the case $k=0$. The case $k=0$ is discussed in [27, Thm. 2.9.3] for the case of a half-space. The generalization to Lipschitz domains can be found, for example, in [1, Lemma 1.10].

2.2 An interpolation argument

The following technical result, which is of independent interest, will be used to reduce regularity assumptions to $B^{m+1/2}_{2,1}(\Omega)$.

Lemma 2.2

Let $X_{1}\subset X_{0}$ be two Banach spaces with continuous embedding. Let $q\in[1,\infty]$, $\theta\in(0,1)$. Define (by the real method of interpolation) $X_{\theta}:=(X_{0},X_{1})_{\theta,q}$ for $\theta\in(0,1)$. Let $0<\theta_{1}<\theta_{2}<\cdots<\theta_{n}<1$ be fixed and assume that $l\in X_{0}^{\prime}$ satisfies for some $C_{0}$, $C_{1}$, $\varepsilon>0$

[TABLE]

Then there exists a constant $C>0$ that is independent of $\varepsilon$ such that

[TABLE]

Proof. We start with the special case $n=1$ and we abbreviate $\theta=\theta_{1}$. Let $f\in X_{\theta}$. By definition of the $K$-functional we may choose $\widetilde{f}\in X_{1}$ with

[TABLE]

Using the linearity of $l$, we can bound

[TABLE]

We now use the bound $\|f-\widetilde{f}\|_{X_{\theta}}\leq 3K(\varepsilon,f)$ from [7, eqn. (2.8)] and then $K(\varepsilon,f)\leq C\varepsilon^{\theta}\|f\|_{X_{\theta}}$ (see, e.g., [27, Thm. 1.3.3]) to conclude

[TABLE]

We now consider the general case $n>1$. We choose $\widetilde{f}$ as in (12) and proceed as above to get

[TABLE]

In order to treat the terms involving $\|\widetilde{f}\|_{X_{\theta_{i}}}$ for $i\geq 2$, we use the reiteration theorem to infer $X_{\theta_{i}}=(X_{\theta_{1}},X_{1})_{s_{i},q}$, where $s_{i}\in(0,1)$ is given by

[TABLE]

Next, the interpolation inequality $\|\widetilde{f}\|_{X_{\theta_{i}}}\leq C\|\widetilde{f}\|_{X_{\theta_{1}}}^{1-s_{i}}\|\widetilde{f}\|_{X_{1}}^{s_{i}}$ together with the elementary bound $ab\leq a^{p}+b^{q}$ ($a$, $b>0$, $p$, $q>1$ with $1/p+1/q=1$) gives

[TABLE]

Inserting this result in (13), we get together with (12)

[TABLE]

Reasoning as in the case $n=1$ now allows us to conclude the argument.

2.3 Elliptic regularity

Lemma 2.3

Let $\Omega\subset{\mathbb{R}}^{d}$ be a bounded Lipschitz domain with a smooth boundary. Let $m\in{\mathbb{N}}$ and $k\in{\mathbb{N}}_{0}$. Then there is $C_{\Omega,m,k}$ depending only on $\Omega$, $m$, $k$ such that the following is true: If $g\in H^{-m+k}(\Omega)$ and $u$ is the (variational) solution of the Dirichlet problem

[TABLE]

then $u\in H^{m+k}(\Omega)$ with the a priori bound

[TABLE]

Proof. This regularity result is a special case of a more general result for the regularity of solutions of elliptic systems, [2, 3]. Self-contained proofs of this result can also be found, for example, in [30, Sec. 20] and in [19, Chap. 2, Thm. 8.2].    The minimum norm extension $E^{\Omega}:H^{m}(\Omega)\rightarrow{\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$ satisfies

[TABLE]

However, for smooth $\partial\Omega$, it has additional mapping properties:

Corollary 2.4

Let $\Omega$ be a bounded Lipschitz domain with a smooth boundary and let $\overline{\Omega}$ be contained in the (open) ball $B_{R}(0)$ of radius $R$ centered at [math]. For each $j\in\{0,\ldots,m\}$ there is a constant $C_{j,\Omega}>0$ depending only on $j$, $\Omega$, and $R$ such that the following is true for the minimum norm extension $E^{\Omega}:H^{m}(\Omega)\rightarrow{\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$: It is also a bounded linear map $H^{m+j}(\Omega)\rightarrow H^{m+j}(B_{R}(0)\setminus\overline{\Omega})$ and, with $\gamma_{0}^{c}$ denoting the trace operator for $B_{R}(0)\setminus\overline{\Omega}$,

[TABLE]

Proof. We write $\widetilde{\Omega}:=B_{R}(0)\setminus\overline{\Omega}$. The operator $E^{\Omega}$ is clearly a bounded linear map $E^{\Omega}:H^{m}(\Omega)\rightarrow H^{m}(\widetilde{\Omega})$. From Lemma 2.3, we also see that $E^{\Omega}$ maps $H^{2m}(\Omega)$ boundedly into $H^{2m}(\widetilde{\Omega})$: We denote by $E$ the universal extension operator of [25, Chap. VI, 3], which we may choose such that $\operatorname*{supp}Ef\subset B_{R}(0)$. Next, we write $E^{\Omega}f$ in the form $E^{\Omega}f=Ef+u$, where $Ef\in H^{2m}(\widetilde{\Omega})$ (since $f\in H^{2m}(\Omega)$) and $u$ solves the differential equation

[TABLE]

Lemma 2.3 then gives $u\in H^{2m}(\widetilde{\Omega})$ with the a priori estimate $\|u\|_{H^{2m}(\widetilde{\Omega})}\leq C\|\Delta^{m}Ef\|_{L^{2}(\widetilde{\Omega})}\leq C\|Ef\|_{H^{2m}(\widetilde{\Omega})}\leq C\|f\|_{H^{2m}(\Omega)}$. We have thus obtained

[TABLE]

An interpolation argument then gives us

[TABLE]

By the trace theorem (Lemma 2.1), we arrive at $\|\nabla^{j+m}E^{\Omega}f\|_{L^{2}(\partial\Omega)}\leq C\|f\|_{B^{m+j+1/2}_{2,1}(\Omega)}$ for $j=0,\ldots,m-1$.

2.4 PDE-based proof of Proposition 1.3

Lemma 2.5

Let $\Omega$ be a Lipschitz domain. Then

[TABLE]

Proof. We exploit that $\Delta^{m}(E^{\Omega}f-If)=0$ in ${{\Omega^{c}}}$. To that end, let again $E$ be the universal extension of operator of [25, Chap. VI, 3]. We write $E^{\Omega}f-If=E(f-If)+\delta$ for some $\delta\in{\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$ with $\delta|_{\Omega}=0$. We get

[TABLE]

where we used integration by parts and that $\delta|_{\Omega}\equiv 0$; the integration by parts does not produce any terms “at infinity” since $C^{\infty}_{0}({\mathbb{R}}^{d})$ is dense in ${\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$ (in the sense described in [29, Thm. 10.40]) and thus $\delta$ can be approximated by such compactly supported functions. The continuity of $E$ implies

[TABLE]

and the reduction to a seminorm follows from the Deny-Lions Lemma and fact that $I$ reproduces polynomials of degree $m-1$.    The solution $If$ of the minimization problem (1) satisfies the orthogonality condition

[TABLE]

since $E^{\Omega}f-If\in{\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$ and $(E^{\Omega}f-If)(x_{i})=f(x_{i})-If(x_{i})=0$, $i=1,\ldots,N$. Therefore,

[TABLE]

These last two terms are treated separately in Lemmas 2.6, 2.7. Inserting (19), (21) in (18) we get

[TABLE]

which readily implies (6) of Proposition 1.3. The bound (7) follows from (6) and an interpolation argument since the reiteration theorem asserts for $0<\delta<1/2$ that $H^{m+\delta}(\Omega)=(H^{m}(\Omega),B^{m+1/2}_{2,1}(\Omega))_{2\delta,2}$ and $|E^{\Omega}f-If|_{H^{m}({\mathbb{R}}^{d})}\leq C\|f\|_{H^{m}(\Omega)}$, which follows from combining (17) and (14).

Lemma 2.6

Let $\Omega$ be a Lipschitz domain. Then:

[TABLE]

Proof. Let $\widetilde{f}\in H^{m+1}(\Omega)$. Integration by parts once gives

[TABLE]

The multiplicative trace inequality $\|z\|^{2}_{L^{2}(\partial\Omega)}\lesssim\|z\|_{L^{2}(\Omega)}\|z\|_{H^{1}(\Omega)}$, Corollary 1.2 with $k=m-1$, and the trace estimate $\|\nabla^{m}z\|_{L^{2}(\partial\Omega)}\lesssim\|z\|_{B^{m+1/2}_{2,1}(\Omega)}$ yield

[TABLE]

We conclude that the linear functional $\widetilde{f}\mapsto B_{\Omega}(f-If,\widetilde{f})$ satisfies

[TABLE]

since $B^{m+1/2}_{2,1}(\Omega)=(H^{m}(\Omega),H^{m+1}(\Omega))_{1/2,1}$ Lemma 2.2 implies the estimate (19).    We now turn to the second part of (20). The key step is to observe that the minimum norm extension $E^{\Omega}f$ satisfies the homogeneous differential equation $\Delta^{m}E^{\Omega}f=0$ in ${\Omega^{c}}$.

Lemma 2.7

Let $\Omega$ be a bounded Lipschitz domain with a sufficiently smooth boundary. Then:

[TABLE]

Proof. Let $\widetilde{f}\in H^{2m}(\Omega)$. By Corollary 2.4, we have $E^{\Omega}\widetilde{f}\in H^{2m}(B_{R}(0)\cap{\Omega^{c}})$ for every $R>0$ sufficiently large. Furthermore, $\Delta^{m}E^{\Omega}\widetilde{f}=0$ in ${\Omega^{c}}$. Next, $m$-fold integration by parts yields

[TABLE]

The integration by parts does not produce any terms “at infinity” since $C^{\infty}_{0}({\mathbb{R}}^{d})$ is dense in ${\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$ (in the sense described in [29, Thm. 10.40]) and thus $E^{\Omega}f-If\in{\mbox{\rm BL}^{m}({\mathbb{R}}^{d})}$ can be approximated by such compactly supported functions.

Since $\nabla^{j}E^{\Omega}f=\nabla^{j}f$ on $\partial\Omega$ for $j=0,\ldots,m-1$, we use again the multiplicative trace inequality and Corollary 1.2 to get

[TABLE]

We reduce the regularity requirement on $\widetilde{f}$ by applying Lemma 2.2 to $\widetilde{f}\mapsto B_{{\Omega^{c}}}(E^{\Omega}f-If,E^{\Omega}\widetilde{f})$: We observe that the reiteration theorem of interpolation allows us to identify

[TABLE]

hence, we get (21) from an application of Lemma 2.2 with $X_{0}=H^{m}(\Omega)$, $X_{1}=B^{2m-1/2}_{2,1}(\Omega)$ and $\varepsilon=h^{m-1/2}$ since we have additionally the stability bound $|B_{{\Omega^{c}}}(E^{\Omega}f-If,E^{\Omega}\widetilde{f})|\leq C|f-If|_{H^{m}(\Omega)}\|\widetilde{f}\|_{H^{m}(\Omega)}$ by Lemma 2.5 and (16).

Remark 2.8** (Generalization to Lipschitz domains)**

The proof Proposition 1.3 relies on three ingredients: a) integration by parts arguments to treat $B_{\Omega}$, b) the approximation properties given in [11] of the thin plate spline interpolation operator $I$, and c) regularity properties of $u:=E^{\Omega}f$. Ingredients a) and b) are already formulated for Lipschitz domains. However, the regularity properties of $u=E^{\Omega}f$ are delicate in their generalization to the case of Lipschitz domains. We note that $u$ solves in ${\Omega^{c}}$ the Dirichlet problem

[TABLE]

*and [28, 23, 9] show a shift theorem by $1/2$ in the sense that for $f\in B^{m+1/2}(\partial\Omega)$, one can control $\nabla^{j}u|_{\partial\Omega}$ for $j=0,\ldots,m$. This together with careful integration by parts arguments for the treatment of $B_{\Omega^{c}}$ allow for an extension of the proof of Proposition 1.3 to Lipschitz domain and will be given in [20].  *

3 Numerical example

We illustrate Proposition 1.3 for the case $m=d=2$, i.e., the classical thin plate splines. We employ uniformly distributed nodes on two geometries, the unit square $\Omega_{1}=(0,1)^{2}$ and the L-shaped domain $\Omega_{2}=(-1/2,1/2)^{2}\setminus[0,1/2]^{2}$. As usual, we denote $r:x\mapsto\|x\|_{2}$. We interpolate 4 functions with different characters: the functions $r^{1.05}$ and $r^{2.76}$, which are, for any $\varepsilon>0$, in $H^{2.05-\varepsilon}$ and $H^{3.76-\varepsilon}$, respectively, and the smooth functions $\exp(xy)$ and $F(x,y)$, where the so-called Franke function $F$ is given by

[TABLE]

The results are presented in Fig. 1 and corroborate the assertions of Proposition 1.3, which read, for $m=2$, $\|f-If\|_{L^{\infty}(\Omega)}\leq Ch^{1+\delta}\|f\|_{H^{2+\delta}(\Omega)}$ with $\delta\in[0,1/2)$ and $\|f-If\|_{L^{\infty}(\Omega)}\leq Ch^{3/2}\|f\|_{B^{5/2}_{2,1}(\Omega)}$. These numerical results were first presented at the conference [21].

3.1 ${\mathcal{H}}$-matrix techniques for solving the TPS interpolation problem

The numerical solution of the thin plate interpolation problem is numerically challenging since the system matrix is fully populated. Nevertheless, several approaches for fast solution techniques exist. For example, the matrix-vector multiplication can be realized in log-linear complexity using techniques from fast multipole methods. This leads to efficient solution strategies based on Krylov subspace methods provided suitable preconditioners are available. We refer to [29, Sec. 15], [8, Sec. 7.3] as starting points for a literature discussion. For our calculations, we employed related techniques based on the concept of ${\mathcal{H}}$-matrices, [14, 15]. ${\mathcal{H}}$-matrices come with an (approximate) factorization that can either be used as a solver (if the approximation is sufficiently accurate) or as a preconditioner in an iterative environment. The latter use has been advocated, in a different context, for example, in [4, 13].

For the case $m=2=d$, the interpolation problem (3) results in a linear system of equations of the form

[TABLE]

The matrix ${\mathbf{P}}^{N\times 3}$ is obtained by selecting a basis $\{b_{1},b_{2},b_{3}\}$ of ${\mathbb{P}}_{1}$ (e.g., $\{1,x,y\}$) and setting ${\mathbf{P}}_{i,j}=b_{j}(x_{i})$. The vector ${\mathbf{f}}\in{\mathbb{R}}^{N}$ collects the values $f(x_{i})$, the vector ${\mathbf{c}}\in{\mathbb{R}}^{N}$ the sought coefficients $c_{i}$, and the vector $\lambda\in{\mathbb{R}}^{3}$ is the Lagrange multiplier for the constrained problem (3). The function $\phi_{2}(z)=z^{2}\log z$ is smooth for $z>0$. Lemma 3.1 below shows that the function $(x,y)\mapsto\phi_{2}(\|x-y\|_{2})$ can be approximated by a polynomial, which is in particular a separable function, i.e. a short sum of products of functions of $x$ and $y$, only. This in turn implies that the fully populated matrix ${\mathbf{G}}$ can in fact be approximated as a blockwise low-rank matrix, in particular in the form of an ${\mathcal{H}}$-matrix, [14, 15].

By forming a Schur complement, the linear system of (24) can be transformed to SPD form. To that end, we select three points and rearrange the problem (24) as

[TABLE]

where the vectors ${\mathbf{c}}_{1}$, ${\mathbf{f}}_{1}\in{\mathbb{R}}^{3}$ and ${\mathbf{c}}_{2}$, ${\mathbf{f}}_{2}\in{\mathbb{R}}^{N-3}$ result from the permutations. The Schur complement

[TABLE]

is SPD. We computed an (approximate) Cholesky factorization of ${\mathbf{S}}$ using the library HLib [5]. This factorization can be employed as a preconditioner for a CG iteration. The ${\mathcal{H}}$-matrix structure of ${\mathbf{S}}$ was ensured by so-called geometric clustering of the interpolation points. Specifically, we used this hierarchical structure to set up $G_{22}$ by approximating its entries with the Chebyshev interpolant as described in Lemma 3.1. In the interest of efficiency, the thus obtained ${\mathcal{H}}$-matrix approximation of $G_{22}$ was further modified by using SVD-based compression of blocks as well as coarsing of the block structure (these tools are provided by HLib). The matrix ${\mathbf{S}}$ is a rank-$3$ update of the matrix $G_{22}$, which can also be realized in HLib.

Lemma 3.1

Let $\eta>0$ be given. For any (closed) axiparallel boxes $\sigma$, $\tau\subset{\mathbb{R}}^{2}$ and a polynomial degree $p\in{\mathbb{N}}_{0}$ denote by $I^{Cheb}_{p}:C(\sigma\times\tau)\rightarrow{\mathbb{Q}}_{p}$ the tensor product Chebyshev interpolation operator associated with $\sigma\times\tau$. Then there are constants $C$, $b>0$ depending only on $\eta$ such that under the condition $\max\{\operatorname*{diam}(\sigma),\operatorname*{diam}(\tau)\}\leq\eta\operatorname*{dist}(\sigma,\tau)$ there holds

[TABLE]

Proof. The proof follows with the tool developed in [6]. Consider $Q:=\prod_{i=1}^{n}[a_{i},b_{i}]\subset{\mathbb{R}}^{n}$ and a function $f\in C(Q;{\mathbb{C}})$. Denote by $\Lambda_{p}$ the Lebesgue constant for univariate Chebyshev interpolation (note that $\Lambda_{p}=O(\log p)$). Introduce, for each $x\in Q$ and each $i\in\{1,\ldots,n\}$, the univariate function $f_{x,i}:[-1,1]\rightarrow{\mathbb{C}}$ by $f_{x,i}(t):=f(x_{1},\ldots,x_{i-1},(a_{i}+b_{i})/2+t(b_{i}-a_{i})/2,x_{i+1},\ldots,x_{n})$. Then, standard tensor product arguments [6, Lemma 3.3] show that the tensor product Chebyshev interpolation error is bounded by

[TABLE]

The best approximation problems $\inf_{\pi\in{\mathbb{P}}_{p}}\|f_{x,i}-\pi\|_{L^{\infty}(-1,1)}$ in turn lead to exponentially small (in $p$) errors, provided the holomorphic extensions of the functions $f_{x,i}$ can be controlled. We show this for the case $f(x_{1},x_{2},x_{3},x_{4})=\phi_{2}(\|(x_{1},x_{2})-(x_{3},x_{4})\|_{2})$ under consideration here. Note that $f_{x,1}(t)=\phi_{2}(\|{\mathfrak{d}}-t{\mathfrak{p}}\|_{2})$, where ${\mathfrak{d}}=((a_{1}+b_{1})/2-x_{3},x_{2}-x_{4})^{\top}$ and ${\mathfrak{p}}=((a_{1}-b_{1})/2,0)^{\top}$. Note $\|{\mathfrak{d}}\|_{2}\leq(1+\eta)\operatorname*{dist}(\sigma,\tau)$ and $\|{\mathfrak{p}}\|_{2}\leq 1/2\max\{\operatorname*{diam}(\sigma),\operatorname*{diam}(\tau)\}\leq\eta/2\operatorname*{dist}(\sigma,\tau)$. As is shown in [6, Lemma 3.6, proof of Thm. 3.13], the holomorphic extension of the function ${\mathfrak{n}}:t\mapsto\|{\mathfrak{d}}-t{\mathfrak{p}}\|_{2}$ is holomorphic on $U_{r}:=\cup_{t\in[-1,1]}B_{r}(t)$ with $r=\operatorname*{dist}(\sigma,\tau)/\|{\mathfrak{p}}\|_{2}\geq 2/\eta$ and maps into the left half plane ${\mathbb{C}}_{+}=\{z\in{\mathbb{C}}\,|\,\operatorname*{Re}z>0\}$. We note that $\sup_{z\in U_{r}}|{\mathfrak{n}}(z)|\leq\|{\mathfrak{d}}\|_{2}+r\|{\mathfrak{p}}\|_{2}\leq(2+\eta)\operatorname*{dist}(\sigma,\tau)$. In view of $\phi_{2}(z)=z^{2}\log z$, we conclude $\sup_{z\in U_{r}}|f_{x,i}(z)|\leq C(\operatorname*{dist}(\sigma,\tau))^{2}(1+|\log\operatorname*{dist}(\sigma,\tau)|)$ for a constant $C>0$ that depends solely on $\eta$. We finish the proof by observing that there is $\rho>1$ (depending only on $r$ and thus on $\eta$) such that $U_{r}$ contains the Bernstein ellipse ${\mathcal{E}}_{\rho}$ (see [6, Lemma 3.12]). A classical polynomial approximation result (see, e.g., [6, Lemma 3.11]) concludes the proof.

3.2 Edge effects and concentrating points at the boundary

The convergence behavior of thin plate splines is limited by edge effects. Above, we mentioned that imposing certain boundary conditions on $f$ mitigates this effect. An alternative is to suitably concentrate points near $\partial\Omega$. Without proof, we announce the following result:

Proposition 3.2

Assume that the points $x_{i}$, $i=1,\ldots,N$, satisfy for a $\delta>0$ sufficiently small

[TABLE]

Then, for $f\in H^{m+1}(\Omega)$ there holds $|f-If|_{H^{m}(\Omega)}\leq C\left(h_{min}^{1/2}+h\right)|f|_{H^{m+1}(\Omega)}$.

Inserting the result of Proposition 3.2 in the estimates of Proposition 1.1 shows that a factor $h_{min}^{1/2}+h$ can be gained in the convergence estimates. Fig. 2 presents numerical examples for the square $\Omega_{1}$ and the functions given in Sec. 3. We selected $h_{min}=h^{2}$ and distributed the points so as ensure the condition

[TABLE]

For the present case $d=2$, it can then be shown that the number of points $N$ is $O(h^{-2})$, which is also illustrated in Fig. 2.