[0,1/2) [1/4, 1/2)

[0,1/8) [1/8, 1/4) [1/4, 3/8) [3/8, 1/2) [1/2, 5/8) [5/8, 3/4) [3/4, 7/8) [7/8, 1)

Fig. 15.2 Hierarchical decomposition of the panel [0, 1] into smaller panels. Four levels of the hier-

archy are shown.

consider it as constant and we need therefore O(N ) time to compute the coef¬cients {βk }

and constant time for each evaluation of s.

One could say that we have averaged the information given at the sources x j to one single

information given at the center t0 of the panel. In this sense, (15.2) is a unipole expansion

of . Note that in the case of translation-invariant kernels it is easy to obtain the far-¬eld

expansion of around any t0 if the far-¬eld expansion around zero is known. To see this,

suppose that (15.2) is the far ¬eld expansion of (x, t) = (x ’ t) around zero. Then we

can write

p

(x ’ t) = ((x ’ t0 ) ’ (t ’ t0 )) = φk (x ’ t0 )ψk (t ’ t0 ) + R(x ’ t0 , t ’ t0 )

k=1

to derive the far-¬eld expansion around t0 .

In general, the evaluation points are close to at least a few of the centers. In the next step

we will re¬ne our approach to ¬t to this situation. The solution is an hierarchical subdivision

of the region of interest into panels or cells of sources. This is not new to us, since we

encountered this concept in the last chapter. The reader should be reminded in particular

of kd-trees and bd-trees. Now, how can we use this concept here? Let us explain the idea

by an example. Suppose that all data points are contained in the interval [0, 1] and that we

divide this interval hierarchically, into “panels” as shown in Figure 15.2. With every source

panel T we associate the part of the function s that corresponds to the sources in that panel,

by setting

sT = ± j (·, x j ).

x j ∈T

Moreover, we also assign the far-¬eld expansion sT of sT to the panel T .

Then the approximation s to the function s at an evaluation point x is computed by adding

the associated functions sT for the panel that contains x itself and all neighboring panels.

Those panels that are well separated from x contribute only by their far-¬eld expansion sT

to s. Here, we say that a point x is well separated from a panel T if it has distance at least

diam(T ) from T . A panel U is called well separated from a panel T if all points in U are

256 Numerical methods

well separated from T . Since we want to avoid ¬‚oating-point operations, we always use the

largest possible source panel for the far-¬eld series.

Let us return to the example given in Figure 15.2. If we want to approximate s(x) at a

point x in the panel [0, 1/8) we form

s(x) = s[0,1/8) (x) + s[1/8,1/4) (x) + s[1/4,3/8) (x) + s[3/8,1/2) (x) + s[1/2,3/4) (x) + s[3/4,1] (x).

Note, that we use the two level-2 approximants s[1/2,3/4) and s[3/4,1] instead of the four

level-3 approximants s[1/2,5/8) , . . . , s[7/8,1] . This halves the computational complexity in this

case. We can do this because the panels [1/2, 3/4) and [3/4, 1] are well separated from

[0, 1/4), the panel containing x on the same level as them. However, we could not use the

approximant s[1/2,1] because its panel [1/2, 1] is not well separated from the panel [0, 1/2],

which is the panel containing x on the same level.

Similarly, to evaluate s(x) approximately in the panel [3/8, 1/2) we would use

s(x) = s[3/8,1/2) (x) + s[1/4,3/8) (x) + s[1/2,5/8) (x)

+ s[0,1/8) (x) + s[1/8,1/4) (x) + s[5/8,3/4) (x) + s[3/4,1] (x).

Having considered this example, we now return to the general situation and describe how

to set up the data structure and how to evaluate the approximant.

Since we discussed several kinds of tree in the previous chapter, we do not need to go

into the details of setting up the data structure. If we are using kd-trees or bd-trees to

decompose space then we have to use a bucket size greater than 1, which is nonetheless

small when compared to N . The same is true if quadtrees and their higher-dimensional

relatives are used. Another approach would be to ¬x the maximum depth of the resulting

tree. Obviously, for computational reasons this should be O(log N ) if N is the number of

source points, which is always satis¬ed for certain splitting rules in the case of kd-trees.

Now, to set up the data structure we ¬rst choose the accuracy and from this the number

p of necessary terms in the far-¬eld expansion. Then we build the tree as usual, assigning

the points to the panels. For the leaf panels we also store the coef¬cients of the interpolant.

After this, we compute the far-¬eld expansions bottom-up. The reason for this is that we

can use the results from lower levels to compute the expansions for higher levels. For each

panel we have to compute the coef¬cients βk , 1 ¤ k ¤ p. If we do this bottom-up then every

source x j is considered only once so that all coef¬cients can be computed in O(N ) time.

But even if we do it top-down then the complexity is bounded by O(N log N ) and this is

the time we need in any case to build the tree. Hence, building the structure takes O(N log N )

time.

The evaluation of the approximation s to the radial basis sum s can be done recursively,

as will be described in Algorithm 13, and we are going to estimate the computational cost

of a single evaluation of s. The far-¬eld part of the work involves the evaluation of only

a constant number of far-¬eld expansions on every level of the tree, each consisting of p

coef¬cients. Hence the time needed to sum up the far-¬eld expansion is O( p log N ). The

15.1 Fast multipole methods 257

direct part consists of a ¬nite number (this number, as well as the ¬nite number mentioned

in the case of the far-¬eld expansion, depends only upon the space dimension) of panels;

thus, this takes only O(1) time, so that the time needed for one evaluation is O(log N ).

Algorithm 13 Recursive evaluation of a multipole method

Input: Panel T and evaluation point x.

Output: Approximate value s(x).

if T is well separated from x then

Return sT (x).

else

if T is a leaf then

Return sT (x).

else

Apply this algorithm to x and the children of T .

Return the sum of their results.

Lemma 15.1 The data structure for a fast multipole algorithm can be built in O(N log N )

time. One evaluation takes O(log N ) time.

In several cases the evaluation time can be further reduced to O(1), but this will not

bother us any further.

The rest of this section is devoted to giving far-¬eld expansions for some of the most

interesting kernels. Unfortunately, for optimal results this has to be done explicitly for each

kernel and even for each space dimension. In many cases these far-¬eld expansions are

truncated Taylor or Laurent expansions. We will give an example for the case of thin-plate

splines in two dimensions. After this we will show how the fast Gauss transform gives

not only a multipole expansion for the Gaussian but also a general framework for deriving

far-¬eld expansions for other conditionally positive de¬nite functions.

The idea of expanding thin-plate splines in two dimensions is based upon identifying R2

with C. Hence, the Euclidean norm in R2 is the usual norm in C and for the moment we will

denote it by | · |. Let us use the notation φt (z) = |z ’ t|2 log |z ’ t| for z, t ∈ C. The idea

of the far-¬eld expansion is based upon the simple equation φt (z) = |z ’ t|2 [log(z ’ t)]

and the Taylor expansion for log(1 ’ z). This gives the Laurent-type expansion

∞

1 tk

t t

log(z ’ t) = log z 1 ’ = log z + log 1 ’ = log z ’ ,

k zk

z z k=1

which holds whenever |t| < |z|. This proves the ¬rst part of

258 Numerical methods

Proposition 15.2 Let z, t ∈ C and φt (z) = |z ’ t|2 log |z ’ t|. Then for |z| > |t| we

have

∞

1 tk

φt (z) = |z ’ t|2 log z ’

k zk

k=1

∞

(ak z + bk )z ’k ,

= |z| ’ 2 (t z) + |t| log |z| +

2 2

k=0

where bk = ’tak , k ≥ 0, and a0 = ’t, ak = t k+1 /[k(k + 1)], k ≥ 1. Moreover, if the ap-