Fast multipole method

The Barnes–Hut algorithm makes the N-body problem cheap, as long as one percent of error is good enough. A graph layout can live with that. Electrostatics in a protein, or an integral equation in engineering, needs six to twelve correct digits. Barnes–Hut can only reach them by opening more and more cells, until it is slower than the direct sum.

The figure computes the same 4000 potentials three ways. Barnes–Hut at θ = 0.1 reaches an error of $10^{-5}$ and takes 1.6 s, longer than the 1.3 s of the exact direct sum. The fast multipole method (FMM) reaches $6\times10^{-12}$ in 82 ms. This note shows how.

Open the figure full-size.

The idea: let boxes talk to boxes

In 1986, Barnes–Hut had two limits. A far cell acted only through its total mass, so accuracy came only from opening cells. Each body also walked the tree on its own, which kept a $\log N$ factor.

Greengard and Rokhlin removed both limits in 1987 (Greengard and Rokhlin 1987). Their first move keeps more information about a far box: not one number, but a series of $p$ terms. Each extra term buys a fixed number of digits. Their second move lets a far box act on a whole target box at once, through a second series centred on the target. Then a box costs the same, however many bodies it holds, and the total work becomes $O(N)$.

The method needs four steps:

1. describe every box by a multipole expansion of its charges;
2. convert far multipoles into a local expansion around each target box;
3. decide which boxes translate to which, with the interaction list;
4. organise all of it as an upward and a downward pass over the tree.

This note works in 2D with the potential $\phi(z)=\sum_j q_j\log|z-z_j|$, where points are complex numbers $z=x+iy$. That is the setting of the 1987 paper. In 3D the series are spherical harmonics, but the four steps stay the same.

Step 1: one box as one series

The charges $q_i$ at $z_i$ in a box with centre $z_c$ produce, at any point further away than all of them, $$\phi(z)=\mathrm{Re}\Big[\,Q\log(z-z_c)+\sum_{k=1}^{p}\frac{a_k}{(z-z_c)^k}\Big],$$ $$Q=\sum_i q_i,\qquad a_k=-\sum_i\frac{q_i\,(z_i-z_c)^k}{k}.$$ The coefficients $a_k$ are the multipole moments of the box. Computing them from the particles is the operator P2M.

Greengard and Rokhlin bound the error after $p$ terms by $$|\text{error}|\le\frac{A}{c-1}\left(\frac1c\right)^{p},\qquad A=\sum_i|q_i|,\quad c=\frac{|z-z_c|}{r},$$ where $r$ is the radius of a circle around $z_c$ that holds all the charges. So the error falls geometrically with $p$, and it falls faster further away.

In the figure, the dashed circle is that radius. With p = 0 only the total charge is left, which is what Barnes–Hut uses: the error is $10^{-2}$ even four box widths away. With p = 4 it is $10^{-5}$ at two box widths, and with p = 12 it is $8\times10^{-10}$. Inside the circle the series diverges at every $p$: no number of terms helps there.

• Why does the error pattern stay round, although the box is square?
• How many terms does one extra digit cost at two box widths, and at four?

Step 2: collect the far field in a local series

A multipole series still has to be evaluated once for every target body. The second idea removes that cost. Around the centre $z_t$ of a target box, the field of a far box is smooth, so it has a Taylor series: $$\phi(z)\approx\mathrm{Re}\sum_{l=0}^{p} b_l\,(z-z_t)^l .$$ The operator M2L computes the $b_l$ directly from the far box’s multipole moments $a_k$, without touching its particles. With $z_0$ the source centre relative to the target centre, $$b_0=Q\log(-z_0)+\sum_{k=1}^{p}\frac{a_k}{z_0^k}(-1)^k,$$ $$b_l=-\frac{Q}{l\,z_0^l}+\frac{1}{z_0^l}\sum_{k=1}^{p}\frac{a_k}{z_0^k}\binom{l+k-1}{k-1}(-1)^k .$$ Local series from many far boxes simply add up. One set of $p$ coefficients then carries the whole far field of a target box, and every body in the box reads its potential from it.

In the figure the two boxes are five box widths apart. With p = 2 the worst error inside the target box is $4\times10^{-4}$. With p = 8 it is $10^{-9}$, and with p = 16 it reaches the rounding error of double precision. The good region is a disc around the target. It grows with $p$ but never reaches the source box.

Step 3: who talks to whom

The FMM uses a quadtree of uniform depth $L$: level $l$ has $4^l$ boxes. Two boxes on the same level are well separated when they do not touch. For a target box, the space splits into three parts:

• its 8 touching neighbours, plus itself, are the near field, summed particle by particle at the finest level;
• its interaction list holds the children of its parent’s neighbours that do not touch it, at most 27 boxes, and each sends one M2L;
• everything further away was already well separated from the parent, so the parent’s local series carries it.

At depth 3 one leaf needs 7 + 27 translations and 9 near boxes, out of 64 leaves. At depth 4 it needs 7 + 27 + 27 out of 256, and at depth 6 7 + 4 × 27 out of 4096. The work grows by 27 per level, while the number of leaves grows by a factor of 4. Click other leaves: boxes at the edge have shorter lists.

Step 4: the two passes

The four ideas become one algorithm with five operators.

The upward pass starts with P2M at the leaves, then shifts the multipoles up one level at a time with M2M. The shift is exact up to the truncation at $p$ terms, so no particle is visited twice. The downward pass starts at level 2 with the M2L translations. On each finer level, a box takes its parent’s series through L2L and adds the translations from its own list. At the leaves, L2P evaluates the series at the particles and adds the direct near-field sum.

Why the cost is O(N)

• Let each leaf hold about $s$ particles, so there are about $N/s$ leaves and $4N/(3s)$ boxes in all.
• Every box does at most 27 M2L translations of cost $p^2$, one M2M and one L2L. That is $O(N p^2/s)$ in total.
• The near field costs $9s$ terms per particle, so $9Ns$ in total. P2M and L2P cost $p$ per particle.
• The total, $O(N(p^2/s+s))$, is linear in $N$ for a fixed $p$ and leaf size. The best leaf size grows like $p$, and then the cost is $O(Np)$.
• The precision sets $p$. The closest M2L partners are two box widths apart, so in the worst case each term gains about $\log_{10}(1/0.55)\approx0.26$ digits. In the benchmark below the measured gain is about 0.45 digits per term: 7 digits for 16 more terms.

Performance analysis

All three methods run in the same JavaScript. The charges are uniform random unit charges, and the error is the largest potential error, relative to the largest potential. The leaf holds about 16 to 31 particles. The button measures again on your machine, up to N = 8000.

• At N = 16,000 the FMM with 8 terms is 92 times faster than the direct sum. It is also 3 times faster than Barnes–Hut, while its error is 500 times smaller.
• Going from 8 to 16 terms costs little here, because the near field dominates the time at these sizes. The digits come almost for free, as the accuracy wall shows.
• At N = 500 the FMM is not faster than the direct sum (27 ms against 22 ms). The tree and the translations have a fixed cost, so small problems should stay direct.
• The time grows in steps, not smoothly: the tree depth is a whole number, and each new level changes the leaf size by a factor of 4.

1992: Anderson’s rings

• Greengard and Rokhlin derived the five operators for one kernel, the log potential. Every new kernel needs a new derivation, and in 3D the algebra of spherical harmonics is heavy.
• Anderson removed the series algebra (Anderson 1992). A potential that is harmonic outside a circle is fixed by its values on that circle. So each box stores the potential at $K = 2M+1$ points on a ring instead of series coefficients.
• An outer ring rebuilds the potential outside the circle and plays the part of a multipole expansion. An inner ring rebuilds it inside and plays the part of a local expansion. The only kernel the method ever evaluates is the potential itself.
• This is the variant of the Scientific Computing 2 course (worksheet 5). anderson_rings.py checks both rings for 50 unit charges, with the target about 4.7 box widths away:

• In 2000 the FMM was listed among the ten algorithms of the 20th century with the greatest influence on science and engineering (Dongarra and Sullivan 2000).

Code

• fmm-figure.html: the figure of this note. It contains a complete 2D FMM in JavaScript, with Barnes–Hut and the direct sum for comparison, and produced every number above.
• fmm_2d.py: the same algorithm in Python, with the five operators as separate methods.
• anderson_rings.py: Anderson’s rings, cleaned from the course solution. The course weighted the ring sum by $h/2\pi = a/K$, which is correct only for ring radius $a = 1$; the vault code uses $1/K$.

Further reading

• Greengard and Rokhlin’s paper (Greengard and Rokhlin 1987) is short and readable. Its Lemmas 2.1 to 2.5 are the five operators of this note.
• Jeffrey Heer, The Barnes-Hut Approximation (Heer 2017) ends where this note begins, and points to the FMM for higher precision.

Links:

• Barnes–Hut algorithm — the predecessor: particle–cell interactions and only the total-charge term.
• Multigrid method — the other $O(N)$ hierarchical solver. Multigrid works on a differential equation and its grid levels; the FMM works on the integral form and its box levels.
• Molecular dynamics — electrostatics is the long-range force that MD cannot cut off.
• Algorithms and data structures — the hall of fame there lists the FMM.

References:

• Scientific Computing 2, TUM (2020), lecture 5 and worksheet 5.