The Fast Multipole Method: It

1
The Fast Multipole Method Its All about Adding
Functions

• Alan Edelman
• MIT Dept of Mathematics,
• Lab for Computer Science
• FOCM 2002
• Saturday, August 10

2
Outline

• Multipole What are we adding?
• The exclusive add without the multipole
• Representing functions on a computer
• Adding functions and the Pascal Matrix
• Research Opportunities

3
The matrix-vector viewpoint
f V(y, x) q
potentials aty1 , ..., yN
charges atx1 , ..., xN
Vij r(yi - xj)
qj yi-xj
fiS
Example
j
Direct evaluation O(N2) work Greengard, Rokhlin
FMM (1987) (N)
4
Some applications

• N-body Problems
• Potential Evaluation
• Future Fast Fourier Transform
• Divide and Conquer Eigenvalue Solvers
• Polynomial Roots
• More waiting to be found .

5
Outline

• Multipole What are we adding?
• The exclusive add without the multipole
• Representing functions on a computer
• Adding functions and the Pascal Matrix
• Research Opportunities

6
Its all about adding functions
7
What might these functions look like?
f(y)S fj(y)
j
f1(y)
f2(y)
8
Outline

• Multipole What are we adding?
• The exclusive add without the multipole
• Representing Functions on a computer
• Adding functions and the Pascal Matrix
• Research Opportunities

9
The Exclusive Add
yi ?j?i xj

• Why not sum and subtract each element?
• sum(x)-x
• What if NaNs or numbers on different scales are
  present?

10
The Exclusive Add
yi ?j?i xj

• A Good Algorithm is worth seeing five ways!
• 1. Pseudocode
• 2. MATLAB
• 3. On a tree
• 4. Matrix Notation
• 5. Kronecker Product Notation

11
Exclusive Add (pseudocode)
yi ?j?i xj
12
Exclusive Add in MATLAB
yi ?j?i xj
function yexclude(v) nlength(v) if n1,
y0v else wv(12n)v(22n)
Pairwise adds
wexclude(w)
Recur y(12n)wv(22n)
y(22n)wv(12n) Update Adds end
13
Exclusive Add on a Tree
yi ?j?i xj
36 10
26 3 7 11 15 1 2
3 4 5 6 7 8
0 26
10 33 29 25 21 35 34
33 32 3130 29 28
Pairwise sums up the tree
Modify evens/odds down
14
With Matrices
T
0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0


15
Kronecker product Approach
T
0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0


A8 (I4?( )) A4 (I4?( ))T I4 ? ( )
11
11
0 1 1 0
1)Pairwise Sums 2)Recursive Excl Add 3)Update
evens odds
16
Kronecker product Approach
T
0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0


A8 (I4?( )) A4 (I4?( ))T I4 ? ( )
11
11
0 1 1 0
1)Pairwise Sums 2)Recursive Excl Add 3)Update
evens odds
17
Kronecker product Approach
T
0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0


A8 (I4?( )) A4 (I4?( ))T I4 ? ( )
11
11
0 1 1 0
1)Pairwise Sums 2)Recursive Excl Add 3)Update
evens odds
18
Kronecker product Approach
T
0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0


A8 (I4?( )) A4 (I4?( ))T I4 ? ( )
11
11
0 1 1 0
1)Pairwise Sums 2)Recursive Excl Add 3)Update
evens odds
19
The Family
All Algorithms 1)Pairwise Sums 2)Recursive
foo 3)Update
20
Multipole is a doubly exclusive add

• We add functions with poles not numbers
• We exclude ourselves and our nearest neighbors
• For reasons analogous to the floating point
  example
• The function can not be represented accurately
  near the singularities.
• This is 1d, higher dimensions analogous
• Algorithm Pairwise Add, Recursion, Update
  Missing Pieces

21
Outline

• Multipole What are we adding?
• The exclusive add without the multipole
• Representing Functions on a computer
• Adding functions and the Pascal Matrix
• Research Opportunities

22
Rounding functions to p coefficients
Best you can do 1)Pick an interval 2)Pick a
function space inside or outside the
interval. 3)Find the best approximation on that
function space.
23
Rounding functions to p coefficients

• Ideal polynomial Approximation
• Truncated Taylor Series
• f(x) ? ak(x-c)k
• Theory Converges in a disk from c to nearest
  singularity
• Practical Accurate where smooth, far from the
  singularity
• Alternative Polynomial Interpolate rather than
  Taylor

24
Rounding functions to p coefficients

• Ideal Multipole Approximation
• Truncated Multipole Series
• f(x) ? ak/(x-c)k1
• Theory Converges outside a disk from c to
  nearest singularity
• Practical Accurate where smooth, far from the
  singularity
• Alternative Multipole Interpolate

25
Rounding functions to p coefficients

• Examples
• R multipole/Taylor Gu
• R Chebyshev interpolation Dutt, Gu, Rokhlin
• S1 Chebyshev interpolation Dutt, Rokhlin
• R singular funs of int ops Yarvin, Rokhlin
• C multipole/Taylor Greengard, Rokhlin
• C discretized Poisson form Anderson
• C singular funs of int ops Hrycak, Rokhlin
• R3 multipole/Taylor Greengard
• R3 discretized Poisson Anderson
• R3 singular funs of int ops Greengard, Rokhlin
• Any Virtual Charges E, McCorquodale

26
Finite Precision Arithmetic

• Idea adding real numbers
• We all know what this means
• x11
• xe?
• On a computer
• Storage must round to d-bit precision
• Arithmetic need a finite precision algorithm

27
Finite Precision Arithmetic

• Idea adding functions
• We all know what this means
• f(x)sin2(x)cos2(x) 1
• f(x)log(x)log(x) log(x2)
• On a computer
• Storage must round to d-bit precision
• Arithmetic need a finite precision algorithm

28
Functions to add

• Slowly growing singularity
• f(x)q/(x-c)
• f(x)qlog(x-c)
• f(x)qcot(x-c) but not f(x)exp(-(x-c)2)

29
Outline

• Multipole What are we adding?
• The exclusive add without the multipole
• Representing Functions on a computer
• Adding functions and the Pascal Matrix
• Research Opportunities

30
Adding Functions

• Issues with adding polynomials
• common center
• accuracy
• Same issues appear with multipole

31
The Pascal Matrix

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1
• Gil Strangs latest favorite matrix

ij i
32
The Pascal Matrix

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1

ij i
33
The Pascal Matrix

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1

ij i
34
The Pascal Matrix

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1

ij i
35
The Pascal Matrix

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1

ij i
36
The Pascal Matrix

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1

ij i
37
The Pascal Matrix

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1

ij i
38
The Pascal Matrix

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1

ij i
39
The Pascal Matrix

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1

ij i
40
The Pascal Matrix

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1

ij i
41
Pascal and Cholesky(Pascal)

• Lchol(pascal(5))
• 1 0 0 0 0
• 1 1 0 0 0
• 1 2 1 0 0
• 1 3 3 1 0
• 1 4 6 4 1
• Lij( ) i,j 0,,n-1

• Ppascal(5)
• 1 1 1 1 1
• 1 2 3 4 5
• 1 3 6 10 15
• 1 4 10 20 35
• 1 5 15 35 70
• Pij( ) i,j 0,,n-1

i j
ij i
ij j
PLL ( ) ? ( )( )
i k
j j-k
42
Binomial Identities
j i
j

• (1x) j ?i0 ( ) xi ? L

i j
?

• (x-1)-(j1) ?ij ( ) x-(i1) ? L

43
Flipping (multipole?taylor)
w F v
Fij (c-d)-i Pij (d-c)-(j1)
44
Shifting (multipole?multipole)
w S v
? wi /(x-d) i1 ? vj / (x-c) j1
j0
i0
Sij (c-d) (i1) Lij (c-d)-(j1)
45
Shifting (taylor?taylor)
w S v
? wi (x-d)i ? vj(x-c)j
j0
i0
T
Sij (d-c)-i Lij (d-c) j
46
Outline

• Multipole What are we adding?
• The exclusive add without the multipole
• Representing Functions on a computer
• Adding functions and the Pascal Matrix
• Research Opportunities

47
Group Representations
Let V be a vector space of functions of x e.g.
polynomials, rational functions, etc. Let G be a
group acting on x, e.g. rotations,
non-singular matrices. Clearly the map ? from
f(x) to h(x)f(g-1x) is linear. Clearly ?(gh)
?(g) ?(h). We say that ? is a representation of
G.
48
Group Representations Research
Generalize Fast Multipole to Arbitrary
Representations! Algebraic Multipole Theory???