( j) 2 2

( j) ( j)

s f,X ’ sk = s f,X ’ sk’1 ’ θk u k

N() N()

2

( j)

s f,X ’ sk’1 , u k N()

( j) 2

= s f,X ’ ’ .

sk’1

|u k |2

N ()

15.2 Approximation of Lagrange functions 269

( j)

This means for every 1 ¤ k ¤ N ’ q that s f,X ’ sk’1 , u k N ( ) tends to zero as j ’ ∞.

Hence for k = 1 we have s f,X ’ s ( j) , u 1 ’ 0 as j ’ ∞. Next, for k = 2, we ¬nd

( j)

0 = lim s f,X ’ s1 , u 2 N()

j’∞

s f,X ’ s ( j) , u 1 N()

= lim s f,X ’ s ’ u1, u2 .

( j)

|u 1 |2 ( )

j’∞

N N()

Since we already know that s f,X ’ s ( j) , u 1 N ( ) tends to zero, this shows that

s f,X ’ s ( j) , u 2 N ( ) also tends to zero for j ’ ∞. Proceeding in the same way, we can

deduce that s f,X ’ s ( j) , u k tends to zero as j tends to in¬nity for every 1 ¤ k ¤ N ’ q.

However, these inner products can be written as

s f,X ’ s ( j) , u k = »ki s f,X (xi ) ’ s ( j) (xi ) ,

N()

i∈Lk

so that we now have

»ki f i ’ s ( j) (xi ) = 0, 1 ¤ k ¤ N ’ q.

lim (15.5)

j’∞

i∈Lk

Finally we consider the equations (15.5) in reverse order, remembering that s ( j) (xi ) = f i ,

N ’ q + 1 ¤ i ¤ N and that »kk > 0. By induction on k, we conclude that s f,X (xk ) ’

s ( j) (xk ) tends to zero as j tends to in¬nity for k = N ’ q, N ’ q ’ 1, . . . , 1.

We will end this section by analyzing the computational complexity of this method and

giving some hints about choosing the index sets Lk . To start with the latter, we ¬rst sort our

points in such a way that the last Q = dim P points are P-unisolvent. For example, in the

case of thin-plate splines in R2 , we would need three noncollinear points. Here, we could

choose the two points having the greatest and least ¬rst coordinate. The third point would

then be the point whose perpendicular distance from the in¬nite straight line through the

other two points is greatest.

Every set Lk should contain these last Q points. Apart from these Q points, it will

contain the data point xk and the q ’ Q ’ 1 data points in {x j : j > k} nearest to xk . With

one of the data structures from the last chapter, we can expect that this will be done in

O(N (q + log N )) time.

To analyze the complexity of the entire algorithm, we ¬rst have to say a few words about

its realization. The main idea is to store and update only the coef¬cients of the current

approximant

N

s(x) = »k (x, xk ) + p(x) (15.6)

k=1

270 Numerical methods

and the current residual vector r . Moreover, we solve the local problems for the local

cardinal functions in advance. This can be done in constant time for each cardinal function

and hence in linear time for all of them.

Since we have an iterative algorithm we restrict ourselves here to estimating the com-

plexity for one iteration. Each iteration consists of O(N ) steps. In each step we have to

update the coef¬cients »k , which can be done in O(1) time since at most Q coef¬cients are

affected. More complicated is the update of the residual vector. An update of each entry

would cost O(N ) time, so that one iteration sums up to O(N 2 ) time. This is comparable to

the time needed for one iteration in the C G-algorithm. However, a combination with the

multipole expansions of the last section improves the results. If we present the approxi-

mant not in the form (15.6) but using a multipole data structure based upon a tree of depth

O(log N ) then updating the coef¬cients will cost us O(log N ) since the constant number

of points is contained in a constant number of leaf panels and thus only a constant number

of coef¬cients in a constant number of panels on each level of the tree have to be updated.

This gives us O(log N ) complexity, which is slightly worse than the complexity in the direct

representation. But the ef¬ciency in updating the residual vector is much greater. The ar-

gument just given shows that the new iteration is concentrated on O(log N ) panels. Hence,

its value is zero apart from these panels. This shows that only O(log N ) residuals have to

be updated. Hence, the time necessary for one complete iteration reduces to O(N log N ).

Finally, we have to mention that an additional O(N log N ) time is necessary for initializing

the data structure, but since this is a one-off problem, it is more than compensated by what

we gain in the updating phases.

15.3 Alternating projections

In the last section we investigated one possible way of using the far-¬eld expansion theory

to derive an ef¬cient method for solving the interpolation equations iteratively. Now we

want to discuss another approach, which again exploits the fact that a radial basis function

interpolant is also an orthogonal projection.

In the following, we will use only positive de¬nite kernels. This is no restriction, since

every conditionally positive de¬nite kernel has an associated positive de¬nite kernel by

Theorem 10.20.

Remember that we have assigned to a set X ⊆ of centers the subspace

VX := span{ (·, x) : x ∈ X }

of the native space N ( ). By Theorem 13.1 we know that the orthogonal projection

P : N ( ) ’ VX is given by the interpolation process, i.e. P f = s f,X . In addition, we

need the following simple result on the orthogonal complement of VX .

Lemma 15.10 The orthogonal complement of VX in N ( ) is given by

⊥

VX = {s ∈ N ( ) : s(x) = 0 for all x ∈ X }.

15.3 Alternating projections 271

Proof An arbitrary element of VX has the form g = N ± j (·, x j ). Using the

j=1

reproducing-kernel property, the inner product can be computed via

N

= ± j s(x j ).

(s, g)N ()

j=1

Hence, on the one hand, if s vanishes on X then it is orthogonal to VX . On the other hand,

if s is orthogonal to VX then we can choose g = (·, x j ), to see that s(x j ) = 0.

The idea of the algorithm starts with a decomposition of the set X of data sites into subsets

X 1 , . . . , X k . These subsets need not be disjoint but their union must be X . The algorithm

starts to interpolate on the ¬rst set X 1 , forms the residual, interpolates this on X 2 , and so on.

After k steps, one cycle of the algorithm is complete and it starts over again. A more formal

description is given in Algorithm 15. The algorithm terminates if a prescribed error bound

> 0 is achieved. Again, we use the notation P j f = s f,X j to emphasize the projectional

character of the interpolant. We denote the sequence of residuals by f j and the sequence of

interpolants by s j .

Algorithm 15 Alternating projection

> 0, right-hand side f |X .

Input: Set of centers X , accuracy

Output: Approximate interpolant.

Set f 0 = f , s0 = 0.

for n = 0, 1, 2, . . . do

for r = 1, . . . , k do

f nk+r = f nk+r ’1 ’ Pr f nk+r ’1

snk+r = snk+r ’1 + Pr f nk+r ’1

if f (n+1)k L ∞ (X ) < then

stop.