Error Estimates For Multilevel Approximation Using

Polyharmonic Splines^{∗}

S.J. Hales ^{†} J. Levesley Department of Mathematics and Computer Science University of Leicester Leicester LE1 7RH, UK

January 5, 2001

Short title: Multilevel polyharmonic approximation

Abstract

Polyharmonic splines are used to interpolate data in a stationary multilevel iterative refinement scheme. By using such functions the necessary tools are provided to obtain simple pointwise error bounds on the approximation. Linear convergence between levels is shown for regular data on a scaled multiinteger grid, and a multilevel domain decomposition method.

Keywords: multilevel approximation, polyharmonic splines, convergence. AMS (MOS) Classification: 41A15, 41A65

1 Introduction

Consider data f_{i} ∈ IR located at distinct points x_{i} ∈ IR^{d} for i = 1, …, N. Let the interpolant be a linear combination of translated basis functions, appended by some polynomial thus

∗Technical Report No. 2001/01, Department of Mathematics & Computer Science, University of Leicester, University Road, Leicester LE1 7RH, England. \begin{matrix}c \\ ◯ \end{matrix} University of Leicester, 2000.

†This research has been supported by EPSRC.

s(x) = ∑_{i=1}^{N} λ_{i}Φ(x − x_{i}) + p_{m−1}(x),  (1)

where p_{m−1} ∈ P_{m−1}^{d}, and s(x_{i}) = f_{i} for i = 1, …, N, with conditions

∑_{i=1}^{N} λ_{i}q(x_{i}) = 0, for all q ∈ P_{m−1}^{d}.  (2)

Here P_{m−1}^{d} is the space of all d−variate polynomials of degree (m − 1) or less, and is of dimension M, where

M = ( \begin{matrix}m − 1 + d \\ d \end{matrix} ) .

Then

p_{m−1} = ∑_{j=1}^{M} α_{j}t_{j}  (3)

where t_{1}, …, t_{M} is a basis for P_{m−1}^{d}. The function Φ(x) = ϕ(∥x∥) is conditionally positive definite of order m, written CPD_{m}. Common choices for ϕ(r) in this setting are given in [10] and [13],

r^{β},  β ∉ 2IN, m ≥ ⌈β/2⌉, polyharmonic spline,
r^{β} log r, β ∈ 2IN, m > β/2, polyharmonic spline,
(c^{2} + r^{2})^{β/2}, β ∉ 2IN, m ≥ ⌈β/2⌉, generalised multiquadric.

The weights in (1), (2) and (3) form the solution vector to the following system

M ( \begin{matrix}λ \\ α \end{matrix} ) = [ \begin{matrix}A & P \\ P^{T} & 0 \end{matrix} ] ( \begin{matrix}λ \\ α \end{matrix} ) = ( \begin{matrix}f \\ 0 \end{matrix} ) .

Interpolating data in this way has been used in a variety of applications; see [1] and [3].

For large values of N the interpolation matrix M becomes ill–conditioned, causing computational difficulties. A method of circumventing this problem is to use a multilevel process of iterative refinement, as will be described in Section 2. A theoretical foundation for such an approach is provided by [9], and numerical viability has been demonstrated in [2] and [4]. As will be seen in Section 3, polyharmonic splines have properties which allow simple pointwise error bounds to be obtained for multilevel schemes. Such estimates are valid for regular data on a scaled multiinteger grid, and a domain decomposition method as described in [4].

2 Iterative Refinement

Although attention will later be restricted to these specialised cases, it is beneficial to explain iterative refinement in the general setting. Consider a sequence of subsets X_{k} = {x_{1}^{(k)}, x_{2}^{(k)}, …, x_{N_{k}}^{(k)}} of X, where X = {x_{1}, x_{2}, …, x_{N}} and ·^{(k)} denotes a quantity pertaining to level k. At the k^{th} level, a refinement to the approximation is made by interpolating the error of the previous level f − (s_{1} + s_{2} + ⋯ + s_{k−1}) at the points in X_{k} with

s_{k}(x) = ∑_{i=1}^{N_{k}} λ_{i}^{(k)}Φ (Δ^{(k)}(x − x_{i}^{(k)})) + p_{m−1}^{(k)}(x),  (4)

for some scale factor Δ^{(k)}. Formally

s_{1}|_{X_{1}} = f|_{X_{1}}
s_{2}|_{X_{2}} = (f − s_{1})|_{X_{2}}
⋮ ⋮
s_{k}|_{X_{k}} = (f − ∑_{i=1}^{k−1} s_{i})|_{X_{k}}.

For clarity, let the error at each level be

e_{k} = f − ∑_{j=1}^{k} s_{j}.

Then s_{k} interpolates e_{k−1} on X_{k}. Therefore the errors at different levels are related via

e_{1}(x) = f(x) − s_{1}(x),
e_{k}(x) = e_{k−1}(x) − s_{k}(x) for k ≥ 2.  (5)

The philosophy of this approach is to use tighter basis functions as the levels increase, so Δ^{(k+1)} > Δ^{(k)}. The main features of the data f are captured at the outset, and finer details are resolved as the levels progress.

The use of different basis functions at each level causes complications in the error analysis. By scaling the points in X_{k} it is possible to obtain an equivalent interpolation problem using a fixed basis. By decomposing the domain we can also keep the number of interpolation points inside each subdomain fixed. The benefit of scaling the points and not the basis function is that this allows the analysis to remain in a fixed function space.

For ease consider the case when Δ^{(k)} = δ^{k−1} for δ > 1, and define

y_{j}^{(k)} = u + δ^{k−1}x_{j}^{(k)}, for j = 1, 2, …, N_{k}, u ∈ IR^{d}.

Further, let y = u + xδ^{k−1} and Y_{k} = u + δ^{k−1}X_{k} = {y_{1}^{(k)}, y_{2}^{(k)}, …, y_{N_{k}}^{(k)}}. We include a translation here because of the application we have in mind in Section 3.

At the k^{th} level, the previous error e_{k−1} is scaled by a factor of δ^{k−1} and interpolated using Φ(y) at the points in Y_{k}. The result is a scaled version of the original interpolant. By introducing the affine scaling y = u + xδ^{k−1}, the interpolant in (4) becomes

s_{k}(x) = ∑_{i=1}^{N_{k}} λ_{i}^{(k)}Φ (δ^{k−1}(x − x_{i}^{(k)})) + p_{m−1}^{(k)}(x)
= ∑_{i=1}^{N_{k}} λ_{i}^{(k)}Φ (y − y_{i}^{(k)}) + p_{m−1}^{(k)} ((y − u)/δ^{k−1})
= ∑_{i=1}^{N_{k}} λ_{i}^{(k)}Φ (y − y_{i}^{(k)}) + q_{m−1}^{(k)}(y) =: t_{k}(y),

where q_{m−1}^{(k)} ∈ P_{m−1}^{d}. For clarity let the scaled error at each level be given by

g_{k}(y) = e_{k}(x) = e_{k} ((y − u)/δ^{k−1}) .  (6)

From (5) the errors at different levels are related via

e_{k}(x) = e_{k−1}(x) − s_{k}(x).

With scaling this becomes

e_{k} ((y − u)/δ^{k−1}) = e_{k−1} ((y − u)/δ^{k−1}) − s_{k}(x),

and from the definition of g_{k} and t_{k} this yields

g_{1} = f − t_{1},
g_{k}(y) = g_{k−1}(y/δ) − t_{k}(y),  k ≥ 2.  (7)

So as s_{k}(x) interpolates e_{k−1}(x) on X_{k}, then t_{k}(y) interpolates the scaled error g_{k−1}(y/δ) on Y_{k}, and is of the same form as (1).

The key aspect is that the interpolant obtained on the scaled set of points Y_{k} using Φ(y) is simply a scaled version of the solution to the problem on X_{k} which uses tighter basis functions Φ(xδ^{k−1}). From (6) the maximum pointwise error of the original and scaled problem is the same, therefore proof of convergence for ∥g_{k}∥_{∞} is also valid for ∥e_{k}∥_{∞}. It will be shown in Section 3 that the error at each level decreases linearly.

Simplification of the error analysis is possible when the points used in Y_{k} are the same as those in Y_{k−1}. Examples of data configurations for which this is satisfied are given.

2.1 Data Configuration

Let X be a scaled multiinteger grid where X = ℤ^{d}. Consider subsets of X given by X_{k} = δ^{1−k}ℤ^{d}, then the scaling y = δ^{k−1}x gives Y_{k} = ℤ^{d}. Therefore the interpolation problem on Y_{k} uses the same points for all k.

A more practical application o this phenomenon occurs in the domain decomposition method presented in [4]. For some fixed K ∈ IN let Ω = {x ∈ IR^{d} : ∥x∥_{∞} ≤ K}, and the associated grid X = {z ∈ ZZ^{d} : z ∈ Ω}. Let S = {−1, 1}^{d}. Then we can construct a refinement of the grid X using the sets

Ω_{σ_{1}} = \frac{\sqrt{d}}{2}σ_{1} + \frac{1}{2}Ω, σ_{1} ∈ S :

X_{1} = ⋃_{σ_{1}∈S} X_{σ_{1}},

where X_{σ_{1}} = {z/2 ∈ Ω_{σ_{1}} : z ∈ ZZ^{d}}. The cardinality of X_{σ_{1}} is exactly that of X.

At the kth level we have the decomposition

Ω_{σ_{1},...,σ_{k−1}} = ⋃_{σ_{k}∈S} Ω_{σ_{1},...,σ_{k}},

where

Ω_{σ_{1},...,σ_{k}} = \frac{\sqrt{d}}{2^{m}}σ_{k} + \frac{1}{2}Ω_{σ_{1},...,σ_{k−1}}, σ_{k} ∈ S.

The associated local grid refinement is X_{σ_{1},...,σ_{k}} = {z/2^{k} ∈ Ω_{σ_{1},...,σ_{k}} : z ∈ ZZ^{d}}, with global grid

X_{k} = ⋃_{σ_{1},...,σ_{k}∈S^{k}} X_{σ_{1},...,σ_{k}}.

Then, each sub-grid X_{σ_{1},...,σ_{k}} is just an affine scaling of the original grid:

X_{σ_{1},...,σ_{k}} = u_{σ_{1},...,σ_{k}} + 2^{−k}X,  (8)

where

u_{σ_{1},...,σ_{k}} = \sqrt{d} ∑_{i=1}^{k} 2^{−i}σ_{i}.

The residual interpolation problem on the grid X_{σ_{1},...,σ_{k}} is solved using the basis function Φ(2^{k}·) with support limited to the given sub-domain. This means that the algorithm is highly parallelizable because there is no communication between sub-domains. After the affine scaling prescribed by (8) described above we have Y_{k} = X for all k and the basis function used in the scaled problem is Φ.

3 Error Analysis

In this section we shall concentrate on the error analysis for the multilevel domain decomposition algorithm, using the refinement described in the previous section. It is straightforward to adapt this analysis to the infinite grid case.

Before addressing the multilevel scheme, let s_{f} be an interpolant of the form (1) to f at points in X. It is necessary to bound the size of the error (f − s_{f}) measured in some norm. Let Φ̂ be the Fourier transform of Φ given by

Φ̂(w) = ∫_{IR_{d}} Φ(x)e^{−ix^{T} w}dx.

The native space of Φ, denoted by G_{Φ}, contains all functions f such that

f̂/\sqrt{Φ̂} ∈ L_{2}(IR^{d}).

If Φ ∈CPD_{m}, then

∥f∥_{Φ} = ∥f̂/\sqrt{Φ̂}∥_{2}  (9)

is a semi-norm on G_{Φ} with kernel P_{m−1}. Then there exists a pointwise error bound of the form

|f(x) − s_{f}(x)| ≤ P (x)∥f∥_{Φ},  (10)

where P is called the power function; see [13]. The following proposition is a direct consequence of Theorem 4.3 of Light and Wayne [7], and gives the form of the power function which we shall use in the sequel.

Proposition 1 Let x_{1}, …, x_{N} be ordered so that x_{1}, …, x_{M} are unisolvent with respect to P_{m−1}. Let p_{1}, …, p_{M} satisfy p_{i}(x_{j}) = δ_{ij}, i, j = 1, …, M. Then,

P (x) = {Φ(0) − 2 ∑_{j=1}^{M} p_{j}(x)Φ(x − x_{j})
+ ∑_{i,j=1}^{M} p_{i}(x)p_{j}(x)Φ(x_{i} − x_{j})}^{1/2} .  (11)

Let P_{δ} be the power function obtained using Φ(·/δ) as a basis function. What we see immediately is that P_{δ} is invariant under affine scalings, where δ is the scaling parameter. In other words, if y_{j} = u + δx_{j}, j = 1, …, N, u ∈ IR^{d}, δ > 0, and y = u + δx, then

P_{δ}(y) = {Φ(0/δ) − 2 ∑_{j=1}^{M} p_{j}^{δ}(y)Φ((y − y_{j})/δ)
+ ∑_{i,j=1}^{M} p_{i}^{δ}(y)p_{j}^{δ}(y)Φ((y_{i} − y_{j})/δ)}^{1/2}
= P (x)  (12)

where p_{i}^{δ}(y_{j}) = δ_{ij}, i, j = 1, …, M. Above, we have used the (straightforward to establish) fact that p_{i}^{δ}(y) = p_{i}(x). Therefore, unlike the case for the usual interpolation problem (see [7]), we will not analyse the power function to infer convergence in the multilevel algorithm. In order to prove pointwise convergence it will be then necessary to show that the error in the ∥ · ∥_{Φ} seminorm decreases with successive levels.

For simplicity of notation we will consider the decrease in error at the first level of refinement, when we decompose Ω with the associated grid X and consider one of the 2^{d} ‘quadrants’ Ω_{σ} with grid X_{σ} for some σ ∈ S. The basis function for interpolation on X is Φ and on X_{σ} is Φ(2·). The residual on Ω is e and on Ω_{σ} is e_{σ}. Thus, if s_{σ} is the interpolant to e on X_{σ}, then e_{σ} = e−s_{σ}. For x ∈ X_{σ}, if we put y = 2x−\sqrt{d}σ, g(y) = e(x), g_{σ}(y) = e_{σ}(x), and t_{σ}(y) = s_{σ}(x), then we can rewrite (7) in the form

g_{σ}(y) = g(y/2) − t_{σ}(y).  (13)

Since t_{σ} is the best approximation to g(·/2) using Φ in G_{Φ} the following holds

∥g(·/2)∥_{Φ}^{2} = ∥g(·/2) − t_{σ}∥_{Φ}^{2} + ∥t_{σ}∥_{Φ}^{2}
= ∥g_{σ}∥_{Φ}^{2} + ∥t_{σ}∥_{Φ}^{2};  (14)

see [8]. Thus,

∥g_{σ}∥_{Φ} ≤ ∥g(·/2)∥_{Φ}.

What we wish to show is that

∥g_{σ}∥_{Φ} ≤ C · ∥g∥_{Φ},  (15)

for C < 1. This inequality can then be applied recursively, so that if g_{σ_{1},...,σ_{k}} is the scaled residual from Ω_{σ_{1},...,σ_{k}}, then

∥g_{σ_{1},...,σ_{k}}∥_{Φ} ≤ C^{k}∥g∥_{Φ}.

Since g is the residual when interpolating f using ϕ we can use (14) again to see that ∥g∥_{Φ} ≤ ∥f∥_{Φ}, so that

∥g_{σ_{1},...,σ_{k}}∥_{Φ} ≤ C^{k}∥f∥_{Φ},  (16)

as long as (15) is satisfied. We can now use (10) to see that, for x ∈ Ω_{σ_{1},...,σ_{k+1}},

|g_{σ_{1},...,σ_{k+1}}(x)| ≤ P_{2^{−k−1}}(x) C^{k}∥f∥_{Φ}
= P (y)C^{k} ∥f∥_{Φ},

where

y = 2^{k+1}(u_{σ_{1},...,σ_{k+1}} + x);

see (8). Since y ∈ Ω which is compact, and P is a continuous function, if (15) holds then

|g_{σ_{1},...,σ_{k+1}}(x)| ≤ AC^{k}∥f∥_{Φ},  (17)

for some constant A which is independent of k and x.

It now remains to show that C < 1 in (15).

3.1 Polyharmonic Splines and the Constant C

Semi-norms generated by some basis functions are able to be written explicitly in terms of the function itself rather than its Fourier transform; see [6]. Let Φ be a polyharmonic spline which is conditionally positive definite of order m. In order that G_{Φ} ⊂ C(IR^{d}), so that interpolation makes sense, it is necessary that m > d/2; see [11]. Then G_{Φ} coincides with the Beppo–Levi space of order m on IR^{d}. In these cases (see for instance [13])

∥f∥_{Φ} = ( ∑_{|α|=m} c_{α}∥D^{α}f∥^{2})^{\frac{1}{2}}  (18)

and the coefficients c_{α} are chosen such that

∑_{|α|=m} c_{α}x^{α} = |x|^{2m}.

Now, if g = f(·/δ) then

∥D^{α}g∥_{2}^{2} = ∫_{IR_{d}} |D^{α}g|^{2}
= δ^{−2m} ∫_{IR_{d}} |D^{α}f(·/δ)|^{2}
= δ^{d−2m} ∫_{IR_{d}} |D^{α}f|^{2}
= δ^{d−2m}∥D^{α}f∥_{2}^{2}.

Hence

∥f(·/δ)∥_{Φ} = δ^{\frac{d}{2} −m}∥f∥_{Φ}.  (19)

Since m > d/2, then by setting C = δ^{\frac{d}{2} −m} this proves (15). It is interesting to note that condition on m to obtain (18) is exactly required to guarantee that C < 1.

Thus, with (17), we have

Theorem 1 Let e_{k} be the error at level k approximating f by the multilevel domain decomposition method, using polyharmonic splines. Then

∥e_{k}∥_{∞,Ω} ≤ AC^{k−1}∥f∥_{Φ},

where A is independent of k and f, and C = 2^{\frac{d}{2} −m}.

4 Numerical Results

As a consequence of the error estimate in Theorem 1, one would expect to see the pointwise error decay linearly at a rate of C = δ^{\frac{d}{2} −m}. A number of examples are given where a test point has been chosen to examine the behaviour of ∥g∥_{∞}.

The multilevel domain decomposition scheme as described in Section 2.1 is implemented to interpolate a variety of datasets, where d = 1, given by the functions

F_{1}(x) = sin(10x) + 3x^{2},
F_{2}(x) = \frac{sin(16x)}{16(x + 1)},
F_{3}(x) = x exp(−x^{2}).

The initial grid X consists of 20 regularly spaced points on [0, 1], and basis functions of the form ∥x∥^{β} log(∥x∥) for β = 2, 4 are used. The scale factor is chosen such that δ = 2. Figures 1 & 2 show the error at an arbitrary test

Figure 1: Convergence of Multilevel Method for r^{2} log(r) in IR

point y = \frac{1}{2}, for initial data F_{3} in both cases. Table 1 contains approximate convergence rates, and the appropriate values of \sqrt{C} = 2^{− \frac{1}{2} (β+1)} for comparison.

In IR^{2} the initial grid X consists of a 10 × 10 regular grid on [0, 1]^{2}. Again the scale factor is chosen such that δ = 2, and the test point is located at y = (\frac{1}{2}, \frac{1}{2}). The initial data F_{4} is given by Franke’s function

F_{4}(x, y) = 0.75e^{−0.25(9x−2)^{2}−0.25(9y−2)^{2}} + 0.75e^{−(9x−2)^{2}/49−(9y−2)^{2}/10}
+0.5e^{−0.25(9x−7)^{2}−0.25(9y−3)^{2}} − 0.2e^{−(9x−4)^{2}−(9y−7)^{2}}.

In IR^{3} then initial grid X consists of a 6 × 6 × 6 regular grid on [0, 1]^{3}. The test point is located at y = (\frac{1}{2}, \frac{1}{2}, \frac{1}{2}), and the data F_{5} is given by

F_{5}(x, y, z) = x^{2}y^{3} cos(10z).

Tables 2 and 3 contain the approximate convergence rate, and the predicted values of C for the two and three dimensional cases.

Theorem 1 predicts the convergence rate well in the one dimensional case. In the two and three dimensional cases however, it becomes clear that the actual pointwise error converges faster than the error bound predicted by the theory. Interestingly, the error appears to decrease at the same rate for the thin plate spline in both IR, IR^{2}, and IR^{3}.

Figure 2: Convergence of Multilevel Method for r^{4} log(r) in IR 12

Dataset Basis Function
β = 2 β = 4
F_{1}  0.353 0.177
F_{2}  0.353 0.175
F_{3}  0.353 0.174
C  0.354 0.177

Table 1: Convergence Rates of Multilevel Method in IR

Dataset Basis Function
β = 2 β = 4
F_{4}  0.344 0.15
C 0.5 0.25

Table 2: Convergence Rates of Multilevel Method in IR^{2}

Dataset Basis Function
β = 2 β = 4
F_{5}  0.3281 0.086
C  0.707 0.354

Table 3: Convergence Rates of Multilevel Method in IR^{3}

References

[1] M.R.Dubal, Construction of three-dimensional black–hole initial data via multiquadrics. Physical Review D 45 (1992) 1178-1186.

[2] M.S. Floater and A. Iske, Multistep scattered data interpolation using compactly supported radial basis functions, J. Comp. Appl. Math. 73 (1996) 65-78.

[3] R.E. Hagan and E.J. Kansa, Studies of the R parameter in the multiquadric function applied to groundwater pumping. Journal of Applied Science and Computations 1 (1994) 266-282.

[4] S.J. Hales and J. Levesley, Multilevel approximation to scattered data using inverse multiquadrics, in Curve and Surface Fitting: Saint-Malo 1999, Albert Cohen, Christophe Rabut, and Larry L. Schumaker (eds.), Vanderbilt University Press, Nashville, 2000.

[5] E.J. Kansa and R.E. Carlson, Radial basis functions: A class of grid–free, scattered data approximations. Computational Fluid Dynamics Journal 3 (1995) 479-496.

[6] J. Levesley and W. Light, Direct form seminorms arising in the theory of interpolation by translates of a basis function. Adv. Comput. Math. 11 (1999) 161–182.

[7] W. Light and H. Wayne, Spaces of distributions, interpolation by translates of a basis function and error estimates, Numer. Math. 81 (1999), 415–450.

[8] W.R. Madych and S.A. Nelson, Multivariate interpolation and conditionally positive definite functions. Approx. Theory Appl. 4 (1988), 77–89.

[9] F.J. Narcowich, R. Schaback and J.D. Ward, Multilevel interpolation and approximation. Appl. Comput. Harmonic Anal. 7 (1999) 243-261.

[10] M.J.D. Powell, The theory of radial basis function approximation in 1990, in Advances in Numerical Analysis, Volume 2: Wavelets, Subdivision Algorithms and Radial Basis Functions, 105–210, W.A. Light (ed.), Oxford University Press, 1992.

[11] M. Renardy and R.C. Rogers, An Introduction to Partial Differential Equations. Texts in Applied Mathematics 13, Springer-Verlag, 1992.

[12] R. Schaback, Multivariate interpolation and approximation by translates of a basis function, Approximation Theory VIII 1 (College Station, TX, 1995) 491–514.

[13] R. Schaback, Reconstruction of Multivariate Functions from Scattered Data, University of Göttingen, 1997.

[14] H. Wendland, Error Estimates For Interpolation By Compactly Supported Radial Basis Functions Of Minimal Degree, J. Approx. Theory 93 (1998) 258–272.