Unfortunately, some of the constants depend exponentially on the space dimension. Thus,

in higher dimensions and also for non-quasi-uniform data sets, kd-trees and their extension,

bd-trees, seem to have better properties than the ¬xed-grid strategy.

What about the stability of the process? Our results from Chapter 12 show that the condi-

tion number does not depend on the number of centers but on the separation distance. The

separation distance of the local center sets X j is in general of the same size as the separation

distance of the whole set of centers X . So we run into the same problem. Fortunately, for

basis functions with a ¬nite number of continuous derivatives the problems seem to give

much less dif¬culty than predicted. Moreover, if we are working with thin-plate splines

or any other basis function that satis¬es the conditions of Theorem 12.13 we know that in

fact we could instead work with the kernel κ(·, ·). Then, by Corollary 12.14 the condition

number of the local problems turns out to be constant, i.e. independent of the separation

distance.

Figure 15.4 shows the reconstruction of the Beethoven relief using two different partition

techniques in the partition-of-unity approach. The reconstruction on the left uses a regularly

sized local covering. Owing to this local nature, the interpolant becomes zero when the

Fig. 15.4 Partition-of-unity reconstruction of the Beethoven relief.

280 Numerical methods

evaluation points are far from the initial data. The second reconstruction (middle picture)

is based on an overlapping kd-tree decomposition of the bounding box (see the right-

hand picture of Figure 15.4). This reconstruction is better adapted to nonuniform data. For

example, it was used for the reconstructions of the glacier and the dragon in the ¬rst chapter.

However, the background in the middle picture of Figure 15.4 indicates that problems may

come up when different parts of the surface are merged together.

15.5 Multilevel methods

All the numerical methods discussed so far are for globally supported radial basis functions.

In particular all methods except the partition-of-unity technique need a far-¬eld expansion.

Even though compactly supported basis functions can be used in the same way as globally

supported ones, this does not take the local character of these functions into account. Hence,

it is time to discuss ideas for the ef¬cient use of compactly supported functions.

As mentioned earlier, one possibility is to adjust the support radius as a function of the

data density. Instead of using the function , it is better to employ the function δ = (·/δ)

and to choose δ as a function of the ¬ll distance h X, . The choice δ = ch X, would always

lead to a sparse matrix and we know that its condition number would be independent of the

number N of centers and of h X, . Unfortunately, we also know, by the trade-off principle,

that we cannot then expect convergence for h X, ’ 0. Hence, the right choice of the support

radius is a delicate problem and we now want to describe a possible solution.

The idea of a multilevel method is again based an a decomposition of the set of centers

X , but this time in a nested sequence of subsets,

X 1 ⊆ X 2 ⊆ · · · ⊆ X k = X. (15.10)

If X is quasi-uniform then the subsets X j should also be quasi-uniform. Moreover, they

should satisfy q X j+1 ≈ ca q X j and h X j+1 , ≈ ca h X j , with a ¬xed constant ca . A good choice

for ca would be 1/2.

Now, the multilevel method is simply one cycle of the alternating-projection method

discussed earlier, but this time we use compactly supported basis functions with a different

support radius at each level. We could even use different basis functions at the different

levels. Hence, we need to formulate the algorithm more generally. For every 1 ¤ j ¤ k we

choose a basis function j . As in Section 15.3 we denote the interpolation operator by

Pj f = ’ x j ),

cx j ( f ) j (·

x j ∈X j

but using now the basis function j at level j. We will take j as (·/δ j ) with a compactly

supported basis function and scaling parameter δ j proportional to h X j , . A more thorough

discussion will follow.

The idea behind this algorithm is that one starts with a very thin, widely spread, set of

points and uses a smooth basis function to recover the global behavior of the function f . In

15.5 Multilevel methods 281

Fig. 15.5 Multilevel reconstruction of the Beethoven relief.

the next level a ¬ner set of points is used and a less smooth function possibly with a smaller

support is employed to resolve more details, and so on.

As was said before, the algorithm performs one cycle of the alternating-projection algo-

rithm. This means that we set f 0 = f and s0 = 0 and compute

s j = s j’1 + P j f j’1 ,

f j = f j’1 ’ P j f j’1

for 1 ¤ j ¤ k.

Even though the multilevel algorithm resembles the alternating-projection algorithm, the

idea behind it is completely different. The most obvious differences are that we use different

basis functions at each level and that we perform only one cycle. The latter is reasonable

since any further cycle would not change our approximant because the data sets are nested.

Figure 15.5 demonstrates the multilevel algorithm in the case of the Beethoven-relief data

set. We set up ¬ve levels and the support radius on each level was chosen so that on each

level the interpolation matrices had a bandwidth of roughly 70“80 points. The interpolation

equations were solved using an unpreconditioned conjugate gradient method. The ¬rst row

of Figure 15.5 shows the accumulated interpolants s j while the second row shows the level

interpolants P j f j’1 .

Lemma 15.22 The function sk interpolates f on X .

Proof We ¬rst remark that f j |X j ≡ 0 for 1 ¤ j ¤ k by de¬nition. Next, we have again

the relation f ’ s j = f j for 1 ¤ j ¤ k between the interpolants and the residuals. Hence,

sk |X = sk |X k = f |X k ’ f k |X k = f |X k = f |X .

282 Numerical methods

Since f k is zero on X k it is also zero on X 1 ⊆ X k , and thus the interpolant to f k on X 1

is zero and would add nothing to the new approximant. A consequence of this is that we

cannot use the technique of Section 15.3 to prove convergence. Another reason why we

cannot use those ideas is that we are now using different basis functions, meaning that the

interpolation operators are orthogonal projections in different Hilbert spaces with different

norms. When scaling a ¬xed basis function the latter problem can be avoided if the scaling

factor is absorbed into the function that is interpolated. But the ¬rst problem still remains.

Let us have a closer look at the idea of scaling a compactly supported function. We

assume that the decomposition (15.10) of X is such that

h X j+1 , ≈ ca h X j , qX j ≈ h X j ,

and (15.11)

j’1

for all 1 ¤ j ¤ k. This gives h X j , ≈ ca h X 1 , . Hence, if we choose as the support radius

j’1

δ j = ca then the classical volume trick shows that the number of centers in the support of

δ j is bounded by a constant. Thus the interpolation matrix associated with the interpolation

operator on level j has only a constant number of nonzero entries in each row, the constant

being independent of the level j and of the current number of centers. This means that

matrix“vector multiplication can be done in linear time. Moreover, we can interpret inter-

polation on X with δ as interpolation on X/δ = {x1 /δ, . . . , x N /δ} with . The separation

distance of X/δ is obviously given by q X/δ = q X /δ, showing that the condition number of

the interpolation matrices is uniformly bounded and independent of the level. Hence, on

each level a conjugate gradient method would need only a constant number of iterations to

converge.

Proposition 15.23 Suppose that the decomposition (15.10) of the set X = {x1 , . . . , x N }

satis¬es (15.11) for 1 ¤ j ¤ k. If we choose a compactly supported radial basis function

j’1

and set the support on level j as δ j = c a then the multilevel algorithm needs O(k N ) time

to determine the interpolant sk .

Practical tests show that the ¬rst step is crucial for a successful interpolant. They also

show that the number of levels k can in general be chosen much smaller than the number

of data sites N , making the multilevel algorithm very ef¬cient. Nonetheless, for proving

convergence theoretically, k has to tend to in¬nity. By setting Q j := I ’ P j again, we see

that the error has the representation

f ’ sk = f k = Q k Q k’1 · · · Q 1 f.

Hence, the method converges if for example all the Q j have norms less than a constant

c < 1. Unfortunately, convergence has not yet been proven in the general case and we ¬nish

this section by discussing some ideas on how to prove it. But before we do this we want

to point out that an improvement can be achieved if a preconditioning technique is used,

meaning that Q j is replaced by Q j := R j (I ’ P j ) with a smoothing operator R j . But so

far it is not really clear which smoothing operator is a good choice.