274 Numerical methods

that is not contained in X j . Let us consider an element of VY j of the form

s= c j (·, x j ) + (·, z).

x j ∈X j

⊥

This element is in VX j if

c j (x, x j ) = ’ (x, z) for all x ∈ X j ,

x j ∈X j

and this system has a unique solution because is positive de¬nite. Hence s is a nontrivial

⊥ ⊥

element in VY j © VX j . The result for VY j+1 © VY j can be established in a similar way.

With these preparatory steps we are able to show that the projection algorithm described

earlier converges at least linearly.

Theorem 15.15 Let f ∈ N ( ) be given. Suppose that X 1 , . . . , X k are weakly distinct

subsets of ⊆ Rd . Set Y j = ∪i= j X i , 1 ¤ j ¤ k. Denote by s ( j) the approximant after j

k

completed cycles. Then there exists a constant c ∈ (0, 1) such that

s f,Y1 ’ s (n) ¤ cn f N ( ).

N()

Proof We start our iteration with f 0 = s f,Y1 instead of f 0 = f and set f (n) = f nk . In the

⊥

notation of Theorem 15.12 we then ¬nd that U j = VX j and

k

⊥ ⊥

U= VX j = VY1

j=1

⊥

or, more generally, A j = VY j+1 . This means in particular that the orthogonal projection of

⊥

s f,Y1 ∈ VY1 to U = VY1 is given by Qs f,Y1 = 0. Hence we have

s f,Y1 ’ s (n) = f (n) ’ Qs f,Y1

N() N()

¤c s f,Y1 ’ Qs f,Y1

2n

N()

¤c 2n

f N()

by Theorem 15.12, taking s f,Y1 N ( ) ¤ f N ( ) (cf. Theorem 10.26) into account also.

It remains to show that c is strictly less than one or in other words that all angles ± j are

⊥ ⊥

different from zero. ± j is the angle between U j = VX j and A j = VY j+1 and from its de¬nition

we have to consider the intersection space

B j := (U j © A j )⊥ = (VX j © VY j+1 )⊥ = VY j .

⊥ ⊥

A zero angle means that cos ± j = 1 or

⊥ ⊥

1 = sup (u, v)N : u ∈ VX j © VY j , v ∈ VY j+1 © VY j ,

()

N ( ), v ¤1 .

u N()

Since we are working entirely in the space VY1 we can use compactness to ¬nd two elements

⊥ ⊥

u — ∈ VX j © VY j and v — ∈ VY j+1 © VY j both with norm one and (u — , v — )N ( ) = 1. Equality in

15.4 Partition of unity 275

Fig. 15.3 Alternating-projection reconstruction of a relief of Beethoven.

the Cauchy“Schwarz inequality is only achieved if the two elements are linearly dependent,

⊥ ⊥ ⊥

in fact meaning here that u — = v — . Thus u — ∈ VX j © VY j+1 = VY j . However, u — is also in

VY j . Hence it has to be zero, which contradicts the fact that it has norm unity. This means

that all the angles have to be different from zero and therefore that c < 1.

The analysis of the complexity of this algorithm is very similar to the analysis for the

method of approximated Lagrange functions. We will leave the details to the reader.

Figure 15.3 shows some results of the alternating-projection method when applied to a

data set representing a Beethoven relief. This moderate data set consists of about 30 000

regularly distributed points and will serve us as a model in the rest of this chapter. Here we

have subdivided the data set into equal-sized slightly overlapping boxes, so that each box

contains at most 200 centers. The ¬rst picture in Figure 15.3 shows the reconstruction after

a half cycle. The next shows the reconstruction after one complete cycle. Here, the sum of

residuals has already dropped from 57 708 to 1654.63 and the reconstruction already looks

authentic. However, a side view, which is presented in the middle picture, reveals that it is

still far from a satisfactory reconstruction. The last two pictures show the ¬nal result after

20 cycles. The error has dropped to below 10’6 .

15.4 Partition of unity

The idea of the partition-of-unity method is the following. We start with a mildly overlapping

covering { j } M of the region ; we will make the term “mildly overlapping” more

j=1

precise in just a moment. For this covering we choose a partition of unity, i.e. a family of

compactly supported, nonnegative, continuous functions {w j } with M w j = 1 on and

j=1

supp(w j ) ⊆ j . Moreover, we choose for every cell j an approximation space V j . Then

a function f is approximated on each cell by a local approximant s j ∈ V j , and the local

approximants are joined by forming

M

sf = sjwj. (15.7)

j=1

To be more precise, we make the following de¬nition.

276 Numerical methods

De¬nition 15.16 Let ⊆ Rd be a bounded set. Let { j } M be an open and bounded

j=1

covering of . This means that all j are open and bounded and that is contained

in their union. Set δ j = diam( j ) = supx,y∈ j x ’ y 2 . We call a family of nonnegative

functions {w j } M with w j ∈ C k (Rd ) a k-stable partition of unity with respect to the covering

j=1

{ j } if

(1) supp(w j ) ⊆ j,

M

w j ≡ 1 on

(2) ,

j=1

(3) for every ± ∈ Nd with |±| ¤ k there exists a constant C± > 0 such that

0

¤ C± /δ |±|

D± w j L∞( j ) j

for all 1 ¤ j ¤ M.

So far, we have not made further assumptions about the covering , but for ef¬ciency it

is necessary that

I (x) := #{ j : x ∈ j} (15.8)

is uniformly bounded on . Nonetheless, even without this assumption we can give a ¬rst

convergence result.

Theorem 15.17 Let ⊆ Rd be bounded. Suppose that { j } M is an open and bounded

j=1

covering of and {w j } M is a k-stable partition of unity. Let f ∈ C k ( ) be the function to

j=1

be approximated. Let V j ⊆ C k ( j ) be given. Assume that the local approximation spaces V j

have the following approximation property. In each region j © , f can be approximated

by a function s j ∈ V j such that

D± f ’ D± s j ¤ µ j (±).

L∞( © j)