A study of the uniform accuracy of univariate thin plate spline interpolation
aa r X i v : . [ m a t h . NA ] M a y A study of the uniform accuracy ofunivariate thin plate spline interpolation
Aurelian Bejancu and Simon Hubbert Abstract
The usual power function error estimates do not capture the trueorder of uniform accuracy for thin plate spline interpolation to smoothdata functions in one variable. In this paper we propose a new type ofpower function and we show, through numerical experiments, that theerror estimate based upon it does match the expected order. We alsostudy the relationship between the new power function and the Peanokernel for univariate thin plate spline interpolation.
For each γ >
0, define the basis function φ γ : IR → IR by φ γ ( x ) = (cid:26) | x | γ , if γ , | x | γ log | x | , if γ ∈ . Let m γ = ⌊ γ/ ⌋ be the integer part of γ/
2. For any integer n ≥ m γ andany set of values { f , . . . , f n } of a target function f prescribed at the set ofequi-spaced knots { , h, h, . . . , } , where h = 1 /n , Micchelli’s theory [7] ofconditionally positive definite radial basis functions guarantees the existenceof a unique function s h,γ of the form s h,γ ( x ) = n X k =0 a k φ γ ( x − hk ) + m γ X l =0 b l x l , x ∈ IR , (1.1)that satisfies the interpolation conditions s h,γ ( hi ) = f i , i = 0 , , . . . , n, (1.2) Corresponding author. E-mail: [email protected]
1s well as the ‘side’ conditions n X k =0 a k k l = 0 , l = 0 , , . . . , m γ . (1.3)In this paper we focus on the special case γ = 2 corresponding to the thin platespline (TPS) basis function φ ( x ) = | x | log | x | and we investigate the rate atwhich the interpolant (1.1) converges to f uniformly over [0 , h →
0. If f has a Lipschitz continuous third derivative on [0 , ,
1) inherits the maximalconvergence rate O ( h ) obtained by Powell [9] and Buhmann [5] for ‘cardinal’TPS interpolation on the infinite grid h ZZ . Due to boundary effects, however,the uniform norm of the error over the full interval [0 ,
1] decays at the muchslower rate O ( h / ), as illustrated numerically in [2, 11, 6].The usual method for error estimation in radial basis function interpolation,reviewed in the next section, delivers bounds of the form | f ( x ) − s h,γ ( x ) | ≤ c f,γ P h,γ ( x ) , x ∈ [0 , , where P h,γ is the so-called ‘power function’ associated with φ γ , while f belongsto the ‘native space’ generated by φ γ [15, 17]. It is well known that theoreticalconvergence rates based upon bounding P h,γ ( x ) uniformly for x ∈ [0 ,
1] do notmatch the actual rates of decay of the error achieved in numerical experimentsif f has sufficiently many continuous derivatives. This discrepancy was firstobserved by Powell [10] for the bivariate TPS interpolant.For γ = 2, in section 3, we obtain a new error bound which employs a ‘mixedpower function’ M h,µ defined by means of the basis functions φ and φ µ , for µ ∈ (0 , x ∈ [0 , M h,µ ( x ) as h → µ ∈ [3 , h / . This matches exactly the previously known numerical orderof uniform convergence of the error f − s h, on [0 , f . In section 4 we prove that, for µ = 3 and x ∈ [0 , M h, ( x ) is, up to a constant factor, the L -norm ofthe Peano kernel of the error functional at x . Moreover, we provide numericalevidence that the smaller L -norm of this Peano kernel does not in fact decayfaster than the mixed power function when measured uniformly over [0 , O ( h / ) for univariate TPS interpolation to sufficiently smooth targetfunctions. 2 Error estimates via the standard power function
In this section we review the power function technique to obtain error estimatesfor univariate interpolation with the radial basis function φ γ . A key role in thistechnique is played by the generalized or distributional Fourier transform of φ γ . Lemma 2.1. [15, section 8.3]
For each γ > , the generalized Fouriertransform of φ γ satisfies c φ γ ( t ) = A γ | t | γ , t ∈ IR \ { } , for some constant A γ such that ( − m γ +1 A γ > . As above, let m γ = ⌊ γ/ ⌋ , n ≥ m γ , and h = 1 /n. For each x ∈ IR which is notin the knot-set { , h, . . . , } , Micchelli’s theory implies that the quadratic form Q γ,n ( v ) := ( − m γ +1 n X j =0 n X k =0 v j v k φ γ ( hj − hk ) − n X j =0 v j φ γ ( x − hj ) , (2.1)is strictly positive whenever the non-zero vector v = ( v , . . . , v n ) T ∈ IR n +1 satisfies x l = n X j =0 v j ( hj ) l , l = 0 , , . . . , m γ . (2.2)Further, for each j ∈ { , , . . . , n } , let ℓ ( γ ) j,h be the unique function of the type(1.1)–(1.3) which satisfies the Lagrange interpolation conditions ℓ ( γ ) j,h ( ih ) = δ ij , i = 0 , , . . . , n. Then we have the Lagrange representation formula for the interpolant (1.1): s h,γ ( x ) = n X j =0 f j ℓ ( γ ) j,h ( x ) , x ∈ IR , (2.3)as well as the reproduction formula x l = n X j =0 ( hj ) l ℓ ( γ ) j,h ( x ) , x ∈ IR , l = 0 , , . . . , m γ . (2.4)3 roposition 2.2. [17] With the above notations, for each x ∈ IR , the vector v x = (cid:16) ℓ ( γ )0 ,h ( x ) , . . . , ℓ ( γ ) n,h ( x ) (cid:17) T ∈ IR n +1 has the property that it minimizes the quadratic form (2.1) among all non-zerovectors v ∈ IR n +1 that satisfy (2.2) . The minimum value of the quadratic form Q γ,n defines the square of the so-called ‘power function’ P h,γ : IR → [0 , ∞ ), namely P h,γ ( x ) := Q γ,n ( v x ) . (2.5)Note that P h,γ ( hj ) = 0, ∀ j ∈ { , , . . . , n } . Proposition 2.3. [17]
For each x ∈ IR , let Θ x,γ ( t ) := e ixt − n X j =0 ℓ ( γ ) j,h ( x ) e ihjt , t ∈ IR . (2.6) Then we have the absolutely convergent integral representation P h,γ ( x ) = | A γ | π Z IR | Θ x,γ ( t ) | | t | γ dt. (2.7) For each γ >
0, let κ γ = ⌈ γ/ ⌉ be the least integer that is greater thanor equal to γ/ f ∈ C κ γ [0 , f ∗ : IR → IR of f as follows (cf. [3]).By the Whitney extension theorem [16], there exists e f ∈ C κ γ (IR) such that e f ( x ) = f ( x ), for x ∈ [0 , ν be an infinitely differentiable cut-off functionwhich satisfies ν ( x ) = 1 for x ∈ [0 ,
1] and ν ( x ) = 0 for sufficiently large | x | ,and set f ∗ ( x ) := ν ( x ) e f ( x ) , ∀ x ∈ IR . Clearly f ∗ ∈ C κ γ (IR) is compactly supported and coincides with f on [0 , c f ∗ , defined as the continuous function c f ∗ ( t ) := Z IR e − ixt f ∗ ( x ) dx, (cid:12)(cid:12)(cid:12) t κ γ c f ∗ ( t ) (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) d ( f ∗ ) ( κ γ ) ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:13)(cid:13)(cid:13) ( f ∗ ) ( κ γ ) (cid:13)(cid:13)(cid:13) L (IR) , (2.8)for any t = 0. In particular, c f ∗ is integrable over IR, so f ∗ can be recoveredvia the Fourier inversion formula f ∗ ( x ) = 12 π Z IR e ixt c f ∗ ( t ) dt, x ∈ IR . (2.9)Next, let f j := f ( hj ) in (1.2) for j = 0 , . . . , n . Then (2.9) and (2.3) imply theerror representation f ( x ) − s h,γ ( x ) = 12 π Z IR c f ∗ ( t ) Θ x,γ ( t ) dt, x ∈ [0 , , (2.10)where Θ x,γ is given by (2.6). Moreover, as a consequence of (2.8) and thedefinition of κ γ , we have c f,γ := n π | A γ | Z IR | c f ∗ ( t ) | | t | γ dt o / < ∞ , i.e., f ∗ belongs to the so-called ‘native space’ generated by φ γ . Using (2.7) andthe Cauchy-Schwarz inequality in (2.10), we obtain the error bound | f ( x ) − s h,γ ( x ) | ≤ c f,γ P h,γ ( x ) , x ∈ [0 , . (2.11)Further, Wu and Schaback [17] showed that the variational characterization ofthe power function given in Proposition 2.2 impliesmax x ∈ [0 , P h,γ ( x ) ≤ B γ h γ/ , as h → , (2.12)for a constant B γ independent of h . On the other hand, Schaback and Wend-land [14] proved that the exponent γ/ O ( h γ/ ) forthe uniform norm of the error over [0 , In this section we focus on the TPS basis function φ ( x ) = | x | log | x | , i.e. γ = 2. According to (2.12), in this case the power function satisfiesmax x ∈ [0 , P h, ( x ) = O ( h ) , as h → . x ∈ [0 , | f ( x ) − s h, ( x ) | (3.1)is of the magnitude of h / for a sufficiently smooth target function f .To address this discrepancy, we start from the integral representation (2.7): P h, ( x ) = | A | π Z IR | Θ x, ( t ) | | t | dt. (3.2)Note that expression (2.6) satisfies | Θ x, ( t ) | = ( O (cid:16) | t | (cid:17) , as t → ,O (1) , as | t | → ∞ , (3.3)for each fixed x and h. Indeed, since m = 1, by (2.4) we have1 = n X j =0 ℓ (2) j,h ( x ) and x = n X j =0 hj · ℓ (2) j,h ( x ) , for x ∈ IR . (3.4)This fact, together with the series expansion of the exponential, provides thebound (3.3) for t →
0. The bound for | t | → ∞ follows from the triangleinequality.As a consequence of (3.3), the integral (3.2) is still well defined if | t | is replacedby | t | µ , for any µ ∈ (0 , µ = 2. We may thus define the mixed powerfunction M h,µ : IR → [0 , ∞ ) whose square is given by M h,µ ( x ) := | A µ | π Z IR | Θ x, ( t ) | | t | µ dt, x ∈ IR , µ ∈ (0 , . (3.5)Under this integral, the Lagrange functions entering in expression (2.6) ofΘ x, (and generated by the TPS basis function φ ) are combined with thegeneralized Fourier transform of the basis function φ µ (cf. Lemma 2.1).We now let κ = ⌈ µ/ ⌉ , f ∈ C κ [0 ,
1] and, as in subsection 2.2, consider thecompactly supported extension f ∗ ∈ C κ (IR) of f to the whole real axis. Thenthe error analysis of subsection 2.2 recast in terms of the mixed power functionimplies | f ( x ) − s h, ( x ) | ≤ c f,µ M h,µ ( x ) , x ∈ [0 , , µ ∈ (0 , , (3.6)6here c f,µ = n π | A µ | Z IR (cid:12)(cid:12)(cid:12)c f ∗ ( t ) (cid:12)(cid:12)(cid:12) | t | µ dt o / < ∞ . (3.7)This shows that, for a fixed µ ∈ (0 , M h,µ as h → Problem 3.1.
Given µ ∈ (0 , , µ = 2 , does there exist an algebraic decay rateof the mixed power function uniformly over [0 , , i.e., a largest value α µ > such that max x ∈ [0 , M h,µ ( x ) = O ( h α µ ) , as h →
0? (3.8)Before embarking on a numerical answer to this problem, a few remarks are inorder. Firstly, note that, due to (3.7), the above target function f ∈ C κ [0 , f ∗ in the native space generated by thebasis function φ µ , rather than the native space generated by the TPS basisfunction φ as in the standard estimate (2.11) for γ = 2. Thus, for a given µ ∈ (0 ,
4) ( µ = 2) , the resulting mixed power function error bound (3.6) isprecisely what we would expect if we measured the TPS interpolation error fortarget functions in the native space of φ µ ; this approach is investigated in [12]in the context of approximation rather than interpolation. In particular, themixed power function bound (3.6) applies to the smooth ( C ∞ ) target functionsemployed in the numerical experiments of [2, 6, 11].Secondly, recall the two equivalent expressions for the standard power function:the direct form (2.5) and the integral representation (2.7). Letting m = ⌊ µ/ ⌋ ,an application of Theorem 3 from [17] shows that the mixed power functioncan also be expressed as M h,µ ( x ) =( − m +1 n X j =0 n X k =0 ℓ (2) j,h ( x ) ℓ (2) k,h ( x ) φ µ ( hj − hk ) − n X j =0 ℓ (2) j,h ( x ) φ µ ( x − hj ) ! , x ∈ IR . (3.9)Thirdly, note that (3.4) implies that the TPS Lagrange functions ℓ (2) j,h satisfyconstraint (2.2) of the variational problem from Proposition 2.2 with µ in placeof γ . However, the solution to that problem is provided by the values of the7 = 1 / µ = 2 / µ = 1 h − M ( max ) h,µ α h,µ M ( max ) h,µ α h,µ M ( max ) h,µ α h,µ
128 4.774E-01 0.167 1.768E-01 0.333 6.342E-02 0.500256 4.253E-01 0.167 1.404E-01 0.333 4.485E-02 0.500512 3.789E-01 0.167 1.114E-01 0.333 3.171E-02 0.5001024 3.376E-01 0.167 8.842E-02 0.333 2.242E-02 0.5002048 3.007E-01 0.167 7.018E-02 0.333 1.586E-02 0.500 c µ µ ∈ (0 , φ µ . Asa result, the bounding technique [17] that leads to the estimate (2.12) for thestandard power function cannot be applied to obtain estimates on M h,µ .We now turn to a numerical investigation of the behaviour of the mixed powerfunction. For a fixed parameter µ ∈ (0 , µ = 2, we compute an approxima-tion M ( max ) h,µ of the left-hand side of (3.8) for h = 1 /n , starting from n = 128and proceeding as follows:1. For the current mesh-size h and each j ∈ { , , . . . , n } , express the TPSLagrange function ℓ (2) j,h in the form (1.1) and compute its coefficients bysolving the system (1.2)–(1.3), where f i = δ ij , i ∈ { , , . . . , n } .2. Use (3.9) to evaluate the mixed power function at the set of mid-points X eval,h = { h/ , h/ , . . . , − h/ } and determine its maximum value M ( max ) h,µ = max {M h,µ ( x ) : x ∈ X eval,h } .
3. Double n and repeat steps 1–2 as long as n ≤ µ , the values of M ( max ) h,µ satisfy M ( max ) h,µ = c µ h α h,µ , where c µ and α h,µ are also included in the tables. On the basis of these nu-merical results, we are led to the following conjecture.8 = 4 / µ = 5 / h − M ( max ) h,µ α h,µ M ( max ) h,µ α h,µ
128 2.143E-02 0.667 6.258E-03 0.833256 1.350E-02 0.667 3.512E-03 0.833512 8.503E-03 0.667 1.971E-03 0.8331024 5.356E-03 0.667 1.106E-03 0.8332048 3.374E-03 0.667 6.208E-04 0.833 c µ µ ∈ (1 , µ = 7 / µ = 8 / µ = 3 h − M ( max ) h,µ α h,µ M ( max ) h,µ α h,µ M ( max ) h,µ α h,µ
128 1.061E-03 1.167 6.221E-04 1.334 3.327E-04 1.507256 4.727E-04 1.167 2.473E-04 1.334 1.196E-04 1.503512 2.106E-04 1.167 9.828E-05 1.333 4.296E-05 1.5001024 9.381E-05 1.167 3.905E-05 1.333 1.543E-05 1.4982048 4.179E-05 1.167 1.551E-05 1.333 5.534E-06 1.496 c µ µ ∈ (2 , µ = 10 / µ = 11 / h − M ( max ) h,µ α h,µ M ( max ) h,µ α h,µ
128 2.032E-04 1.491 1.661E-04 1.497256 6.995E-05 1.497 5.808E-05 1.499512 2.419E-05 1.501 2.039E-05 1.5001024 8.402E-06 1.503 7.182E-06 1.5012048 2.919E-06 1.505 2.523E-06 1.502 c µ µ ∈ (3 , onjecture 3.2. The mixed power function satisfies the estimate (3.8) withthe algebraic decay rate α µ = ( µ , for µ ∈ (0 , \ { } , , for µ ∈ [3 , . µ = 3 Note that if Conjecture 3.2 can be established for a particular value µ ∈ [3 , f ∈ C [0 , x ∈ [0 , | f ( x ) − s h, ( x ) | = O (cid:16) h / (cid:17) , as h → . (4.1)This would provide a theoretical explanation of the numerical results reportedin [2, 6, 11].In this section, we investigate Conjecture 3.2 for the special case µ = 3. By(3.5) and (3.9), the square of the mixed power function M h, , combining TPSLagrange functions with the cubic basis function φ , is given by M h, ( x ) = A π Z IR | Θ x, ( t ) | | t | dt = n X j =0 n X k =0 ℓ (2) j,h ( x ) ℓ (2) k,h ( x ) | hj − hk | − n X j =0 ℓ (2) j,h ( x ) | x − hj | . (4.2)Figure 1 illustrates the decay of M h, as h →
0. It can be confirmed nu-merically that the decay rate O ( h / ) of M h, suggested by Table 3 appliesuniformly on [0 , M h, to a classical error analysismethod, namely the Peano kernel representation. Let E h,x ( f ) := f ( x ) − s h, ( x ) = f ( x ) − n X j =0 f ( hj ) ℓ (2) j,h ( x ) . For each x ∈ [0 , E h,x is a continuous linear functional on C [0 ,
1] with theusual max norm, and (3.4) implies that the linear polynomials are in the null10 −3 Figure 1: Plot of M h, ( x ) for h − = 16 ,
32 and 64 . space of E h,x . Then, for any f with an absolutely continuous first derivativeon [0 , f ( x ) − s h, ( x ) = Z K h,x ( u ) f ′′ ( u ) du, ∀ x ∈ [0 , , (4.3)where K h,x is the ‘Peano kernel’ given by K h,x ( u ) := ( x − u ) + − n X j =0 ℓ (2) j,h ( x )( hj − u ) + , u ∈ IR . Proposition 4.1.
For each x ∈ [0 , , the mixed power function value M h, ( x ) is a constant multiple of the L [0 , -norm of the Peano kernel K h,x .Proof. The reproduction property (3.4) implies that K h,x is compactly sup-ported on [0 , K h,x ( u ) = 12 | x − u | − n X j =0 ℓ (2) j,h ( x ) | hj − u | , u ∈ IR . Then, by Lemma 2.1, the Fourier transform of the kernel K h,x is the analyticand square integrable function d K h,x ( t ) = A x, ( − t ) t , t ∈ IR , x, is defined by (2.6). Therefore, using the first line of (4.2), theParseval-Plancherel formula, and the compact support of K h,x , we deduce2 πA M h, ( x ) = Z IR | Θ x, ( − t ) | t dt = 4 A Z IR | d K h,x ( t ) | dt = 8 πA Z | K h,x ( u ) | du, which is the required conclusion.As a consequence, we obtain an alternative way of bounding the error f − s h, interms of the mixed power function M h, , by using Cauchy-Schwarz directly inthe right-hand side of the Peano formula (4.3). The resulting bound applies toany f with an absolutely continuous first derivative on [0 ,
1] and f ′′ ∈ L [0 , f ∈ C [0 , µ = 3.Finally, a related question of interest is whether a sharper uniform error boundcan be obtained from (4.3) via H¨older’s inequality | f ( x ) − s h, ( x ) | ≤ B h ( x ) (cid:13)(cid:13) f ′′ (cid:13)(cid:13) L ∞ [0 , , x ∈ [0 , , where f ′′ ∈ L ∞ [0 ,
1] and B h ( x ) := Z | K x,h ( u ) | du. We note that this technique was used by Atkinson [1] in the late 1960s to inves-tigate the error behavior of natural cubic spline interpolant near the endpointsof the unit interval; see also Schaback [13] for a treatment that is closer to ourpresentation. In the case of the TPS interpolant, a numerical answer to thequestion is provided in Table 5, whose entries satisfy B h (cid:18) h (cid:19) = 0 . h β h , and B h (cid:18) − h (cid:19) = 0 . h σ h , i.e., B h decays approximately with the rate O ( h / ) near the endpoints of [0 , O ( h ) near the midpoint. Also, Figure 2 shows thatthe extreme peak value is well approximated by B h (cid:0) h (cid:1) , while all of the lower12eaks decay at the faster rate. Therefore estimating the L -norm of the Peanokernel leads to the same rate of decay O ( h / ) of the uniform error (3.1) asthat predicted in (4.1) by the mixed power function M h, .We conclude the paper by remarking that any theoretical proof of the uniformdecay rate O ( h / ) of M h, ( x ) or B h ( x ) for x ∈ [0 ,
1] will have to rely on specificproperties of the TPS Lagrange functions ℓ (2) j,h , j ∈ { , , . . . , n } . A potentiallyuseful such property is the special case of [4, Theorem 3.1] stating that theLebesgue-type constant max x ∈ [0 , n X j =0 h ℓ (2) j,h ( x ) i admits an upper bound independently of the mesh-size h . It remains an openquestion whether this or other properties of the TPS Lagrange functions canlead to further progress on the above conjectures. −3 Figure 2: Plot of B h ( x ) for h − = 16 ,
32 and 64 . Acknowledgements.
A.B. is grateful to Michael Johnson (Kuwait Univer-sity) for a suggestion which improved the presentation of the numerical resultsand for spotting a fallacious argument in a tentative proof of Conjecture 3.2for µ = 3. The final version has also benefited from the referees’ suggestions.13 − B h (cid:0) h (cid:1) β h B h (cid:0) − h (cid:1) σ h
64 1.024E-04 1.491 3.633E-05 2.001128 3.533E-05 1.498 9.098E-06 2.001256 1.228E-05 1.501 2.293E-06 1.999512 4.289E-06 1.503 5.694E-07 2.0001024 1.502E-06 1.504 1.434E-07 1.999Table 5: Decay of the L -norm of the Peano kernel. References [1] K. E. Atkinson,
On the Order of Convergence of Natural Cubic SplineInterpolation,
SIAM Journal on Numerical Analysis, No.1 (1969),89–101.[2] A. Bejancu,
Convergence Properties of Surface Spline Interpolation,
Ph.D. Dissertation, University of Cambridge, (1999).[3] A. Bejancu,
Local accuracy for radial basis function interpolation onfinite uniform grids,
J. Approx. Theory (1999), 242–257.[4] A. Bejancu, On the accuracy of surface spline approximation and in-terpolation to bump functions,
Proc. Edinburgh Math. Soc. (2001),225–239.[5] M.D. Buhmann, Multivariate cardinal interpolation with radial basisfunctions,
Constr. Approx. (1990), 225–255.[6] S. Hubbert and S. M¨uller, Thin plate spline interpolation on the unitinterval,
Numer. Alg. (2007), 167–177.[7] C.A. Micchelli, Interpolation of scattered data: Distance matrices andconditionally positive definite functions,
Constr. Approx. (1986), 11–22.[8] M.J.D. Powell, Approximation Theory and Methods . Cambridge Uni-versity Press (1981).[9] M.J.D. Powell,
The theory of radial basis function approximation in1990, in Advances in Numerical Analysis, vol. II: Wavelets, subdivision14lgorithms, and radial basis functions, W.A. Light (ed.), Clarendon,Oxford (1992), 105–210.[10] M.J.D. Powell,
The uniform convergence of thin plate spline interpo-lation in two dimensions,
Numer. Math. (1994), 107–128.[11] M.J.D. Powell, Recent research at Cambridge on radial basis functions, in New Developments in Approximation Theory, M. W. Mueller, M.D. Buhmann, D. H. Mache and M. Felten (eds), Birkh¨auser (1999),215–232.[12] R. Schaback,
Approximation by Radial Basis Functions with FinitelyMany Centers,
Constructive Approximation (1996), 331–340.[13] R. Schaback, Radial Basis Functions Viewed from Cubic Splines, in: Multivariate Approximation and Splines, G. N¨urnberger, J.W.Schmidt, and G. Walz (eds.), Birkh¨auser, (1997), 245–258.[14] R. Schaback and H. Wendland,
Inverse and saturation theorems forradial basis function interpolation,
Math. of Comp. (2002), 669–681.[15] H. Wendland, Scattered Data Approximation.
Cambridge UniversityPress (2005).[16] H. Whitney,
Analytic extensions of differentiable functions defined inclosed sets,
Trans. AMS (1934), 63–89.[17] Z. Wu and R. Schaback, Local error estimates for radial basis functioninterpolation of scattered data,