On thin plate spline interpolation

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


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) Here, the seminorm |v| H m (R d ) is induced by the bilinear form v, w m : For m > d/2 and under very mild conditions on the point distribution, a unique minimizer I f exists. The name "thin plate splines" originates from the fact in the simplest case m = d = 2, I f can be represented in terms of translates of the fundamental solution of the biharmonic equation. For general m the interpolant I f can be expressed in terms fundamental solutions of ∆ m : There are constants c i ∈ R, i = 1, . . . , N, and a polynomial π ∈ P m−1 of degree m − 1 such that (with the Euclidean norm · 2 on R d ) where φ m is given explicitly by φ m (r) = r 2m−d log r d even r 2m−d d odd. (4) The representation (3) allows one to reformulate (1) as the problem of finding the coefficients c i and the polynomial π 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 Ω ⊂ R d and points X N = {x i | i = 1, . . . , N} ⊂ Ω, the fill distance h(X N ) is given by Starting with the seminal papers by J. Duchon [12,11] the error f − I f on Ω is controlled in terms of h and the regularity properties of f (on Ω): here, E Ω 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 It is worth noting that the interpolation operator I is a projection so that I( f − I f ) = 0. Proposition 1.1 applied to the function f − I f therefore yields 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 ∈ B m+1/2 2,1 (Ω) and ∂ Ω 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 ∂ Ω in [17,18] to Lipschitz continuity as discussed in more detail in Remark 2.8. Our main result therefore is a simpler proof of: Let Ω ⊂ R d be a bounded Lipschitz domain with sufficiently smooth boundary. Then there are constants h 0 , C 1 , C δ > 0 that depend solely on Ω, m, d, and δ such that for any collection X = {x 1 , . . . , x N } ⊂ Ω with fill distance h := h(X N ) ≤ h 0 there holds In particular, therefore, the estimates of [11,Prop. 3 , 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

Tools
The precise formulation of the minimization problem (1) is based on the classical Beppo-Levi space BL m (R d ), which is defined as We refer to [10] and [29,Sec. 10.5] for more properties of the space BL m (R d ); in particular, [29,Thm. 10.40] for the precise notion). We also need the minimum norm extension E Ω : H m (Ω) → BL m (R d ) given by The minimization property in (8) implies the orthogonality The connection with elliptic PDE theory arises from the fact that E Ω U satisfies an elliptic PDE in Ω c := R d \ Ω: It will be convenient to decompose B(u, v) The trace mapping is continuous H 1/2+ε (Ω) → H ε (∂ Ω) for ε ∈ (0, 1/2]; however, the limiting case ε = 0 is not true; it is true if the Sobolev space H 1/2 (Ω) is replaced with the slightly smaller Besov space B 1/2 2,1 (Ω): Then there exists C > 0 such that the multiplicative estimate u 2 L 2 (∂ Ω) ≤ C u L 2 (Ω) u H 1 (Ω) holds as well as Proof. The case k ≥ 1 in (11)

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 (Ω).
Then there exists a constant C > 0 that is independent of ε such that Proof. We start with the special case n = 1 and we abbreviate θ = θ 1 . Let f ∈ X θ . By definition of the K-functional we may choose f ∈ X 1 with Using the linearity of l, we can bound We now consider the general case n > 1. We choose f as in (12) and proceed as above to get In order to treat the terms involving f X θ i for i ≥ 2, we use the reiteration theorem to infer Inserting this result in (13), we get together with (12) Reasoning as in the case n = 1 now allows us to conclude the argument.

Elliptic regularity Lemma 2.3
Let Ω ⊂ R d be a bounded Lipschitz domain with a smooth boundary. Let m ∈ N and k ∈ N 0 . Then there is C Ω,m,k depending only on Ω, m, k such that the following is true: If g ∈ H −m+k (Ω) and u is the (variational) solution of the Dirichlet problem 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] However, for smooth ∂ Ω, it has additional mapping properties: Proof. We write Ω := B R (0) \ Ω.

PDE-based proof of Proposition 1.3 Lemma 2.5 Let Ω be a Lipschitz domain. Then
To that end, let again E be the universal extension of operator of [25,Chap. VI,3]. We write , where we used integration by parts and that δ | Ω ≡ 0; the integration by parts does not produce any terms "at infinity" since C ∞ 0 (R d ) is dense in BL m (R d ) (in the sense described in [29,Thm. 10.40]) and thus δ can be approximated by such compactly supported functions. The continuity of E implies and the reduction to a seminorm follows from the Deny-Lions Lemma and fact that I reproduces polynomials of degree m − 1.

We conclude that the linear functional
|B (Ω) = (H m (Ω), H m+1 (Ω)) 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 Ω f satisfies the homogeneous differential equation ∆ m E Ω f = 0 in Ω c .

Lemma 2.7
Let Ω be a bounded Lipschitz domain with a sufficiently smooth boundary. Then: Proof. Let f ∈ H 2m (Ω). By Corollary 2.4, we have The integration by parts does not produce any terms "at infinity" since C ∞ 0 (R d ) is dense in BL m (R d ) (in the sense described in [29,Thm. 10.40]) and thus E Ω f − I f ∈ BL m (R d ) can be approximated by such compactly supported functions.
Since ∇ j E Ω f = ∇ j f on ∂ Ω for j = 0, . . . , m − 1, we use again the multiplicative trace inequality and Corollary 1.2 to get We reduce the regularity requirement on f by applying Lemma 2.
We observe that the reiteration theorem of interpolation allows us to identify  and [28,23,9] show a shift theorem by 1/2 in the sense that for f ∈ B m+1/2 (∂ Ω), one can control ∇ j u| ∂ Ω for j = 0, . . . , m. This together with careful integration by parts arguments for the treatment of B Ω c allow for an extension of the proof of Proposition 1.3 to Lipschitz domain and will be given in [20].

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 H -matrices, [14,15]. 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 The matrix P N×3 is obtained by selecting a basis {b 1 , b 2 , b 3 } of P 1 (e.g., {1, x, y}) and setting P i, j = b j (x i ). The vector f ∈ R N collects the values f (x i ), the vector c ∈ R N the sought coefficients c i , and the vector λ ∈ R 3 is the Lagrange multiplier for the constrained problem (3). The function φ 2 (z) = z 2 log z is smooth for z > 0. Lemma 3.1 below shows that the function (x, y) → φ 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 G can in fact be approximated as a blockwise low-rank matrix, in particular in the form of an 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  where the vectors c 1 , f 1 ∈ R 3 and c 2 , f 2 ∈ R N−3 result from the permutations. The Schur complement is SPD. We computed an (approximate) Cholesky factorization of S using the library HLib [5]. This factorization can be employed as a preconditioner for a CG iteration. The H -matrix structure of 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 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 S is a rank-3 update of the matrix G 22 , which can also be realized in HLib.
Lemma 3.1 Let η > 0 be given. For any (closed) axiparallel boxes σ , τ ⊂ R 2 and a polynomial degree p ∈ N 0 denote by I Cheb p : C(σ × τ) → Q p the tensor product Chebyshev interpolation operator associated with σ × τ. Then there are constants C, b > 0 depending only on η such that under the condition max{diam(σ ), diam(τ)} ≤ η dist(σ , τ) there holds sup (x,y)∈σ ×τ Proof. The proof follows with the tool developed in [6]. Consider Q := ∏ n i=1 [a i , b i ] ⊂ R n and a function f ∈ C(Q; C). Denote by Λ p the Lebesgue constant for univariate Chebyshev interpolation (note that Λ p = O(log p)). Introduce, for each x ∈ Q and each i ∈ {1, . . . , n}, the univariate function f x,i : Then, standard tensor product arguments [6,Lemma 3.3] show that the tensor product Chebyshev interpolation error is bounded by The best approximation problems inf π∈P p f x,i − π L ∞ (−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 , t → d − tp 2 is holomorphic on U r := ∪ t∈[−1,1] B r (t) with r = dist(σ , τ)/ p 2 ≥ 2/η and maps into the left half plane C + = {z ∈ C | Re z > 0}. We note that sup z∈U r |n(z)| ≤ d 2 + r p 2 ≤ (2 + η) dist(σ , τ). In view of φ 2 (z) = z 2 log z, we conclude sup z∈U r | f x,i (z)| ≤ C(dist(σ , τ)) 2 (1 + | log dist(σ , τ)|) for a constant C > 0 that depends solely on η. We finish the proof by observing that there is ρ > 1 (depending only on r and thus on η) such that U r contains the Bernstein ellipse E ρ (see [6,Lemma 3.12]). A classical polynomial approximation result (see, e.g., [6,Lemma 3.11]) concludes the proof.

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 ∂ Ω. Without proof, we announce the following result: 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.