The Fast Multipole Method

Picture a million stars, each pulling on every other. The gravitational force on star i depends on the position of every other star j. Compute that naively and you need $n^{2}$ calculations — a trillion operations for a million bodies. Even at a billion operations per second, that is a thousand seconds per time step. Simulation of a galaxy becomes impossible.

Now ask yourself: does a star in the Milky Way really need to know the exact position of every star in Andromeda? Or would it be enough to know that Andromeda is roughly over there, with a certain total mass and a shape that can be summarised by a handful of numbers?

That intuition is the heart of the Fast Multipole Method (FMM), introduced by Leslie Greengard and Vladimir Rokhlin in 1987. By grouping distant particles into clusters and representing each cluster with a compact multipole expansion, the FMM reduces all-pairs force computation from $O(n^{2})$ to $O(n)$. It was named one of the top ten algorithms of the 20th century by Computing in Science & Engineering in 2000.

The key insight is hierarchical: nearby particles must be handled individually, but for particles that are far enough away, the whole group can be summarised — and the error of that summary can be made arbitrarily small by including more terms.

Below is a live simulation of n charged particles. The Brute Force method computes every one of the $n^{2}$ pairs. The FMM approximation groups distant particles into cells and sums their contributions with a multipole expansion, touching only $O(n)$ interactions.

<div class="controls">
  <label>{{particles_label}}: <span id="nLabel">80</span>
    <input type="range" id="nSlider" min="20" max="400" step="20" value="80">
  </label>
  <label>{{grid_cells_label}}: <span id="gLabel">4×4</span>
    <input type="range" id="gSlider" min="2" max="8" step="1" value="4">
  </label>
  <button id="runBtn" type="button">{{randomise_btn}}</button>
</div>
<canvas id="c" width="460" height="220"></canvas>
<div class="stats">
  <div class="stat brute"><span class="dot bf"></span>{{brute_force_ops}} <b id="bfOps">—</b></div>
  <div class="stat fmm"><span class="dot fm"></span>{{fmm_ops}} <b id="fmmOps">—</b></div>
  <div class="stat err"><span class="dot er"></span>{{force_error}} <b id="errVal">—</b></div>
</div>
<p class="note">{{note_text}}</p>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; padding: 14px; }
.controls { display: flex; flex-wrap: wrap; gap: .6rem 1.2rem; align-items: center; margin-bottom: .5rem; }
label { font-size: .88rem; display: flex; align-items: center; gap: .4rem; }
input[type=range] { width: 110px; accent-color: #1d3557; }
button { font: 600 13px system-ui; padding: .35rem .8rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 6px; cursor: pointer; }
canvas { border: 1px solid #cdd9e3; border-radius: 8px; width: 100%; max-width: 460px;
         display: block; background: #f7fafc; }
.stats { display: flex; flex-wrap: wrap; gap: .5rem 1.4rem; margin-top: .55rem; font-size: .85rem; }
.stat { display: flex; align-items: center; gap: .35rem; }
.dot { width: 10px; height: 10px; border-radius: 50%; display: inline-block; }
.dot.bf { background: #e63946; }
.dot.fm { background: #2a9d8f; }
.dot.er { background: #e9a11a; }
b { font-weight: 700; }
.note { font-size: .82rem; color: #555; margin: .55rem 0 0; line-height: 1.45; }

// Code not found

Drag the Particles slider upward and watch the brute-force operation count soar quadratically while the FMM count grows almost linearly. At 200 particles the gap is already visible; at 1 000 it is enormous. The error bar shows how close the approximate total force is to the exact result — typically within 1 % even with a coarse expansion.

The FMM is a solved, proven-correct algorithm — not NP-hard, not open. Its complexity analysis is rigorous:

• Brute force requires exactly n(n−1)/2 pairwise evaluations — $O(n^{2})$. For n = $10^{6}$ that is 5 × $10^{11}$ operations, far beyond real-time.
• Barnes-Hut (1986) was the first tree method: build an octree, then for each particle walk the tree and use a cell's centre-of-mass whenever the cell is far enough away (controlled by the opening angle θ). Cost: $O(n \log n)$, error O(θ²). A huge improvement, but still super-linear.
• The FMM (Greengard & Rokhlin, 1987) goes further. Instead of re-computing a cell's effect for every query particle, the FMM pre-translates the multipole expansion of a far cell into a local expansion at the query particle's cell. This translation is done once per cell pair at each level of the tree, not once per particle. The result is $O(n)$ with a constant that depends on the desired accuracy (number of expansion terms p): roughly $O(p^{4} n)$ in 3-D, where p ≈ 10 gives single-precision accuracy.
• Error control is a key feature: the approximation error decays exponentially with p. You choose p to match the required accuracy, and the algorithm delivers it with a rigorous bound.

The FMM also generalises beyond gravity and electrostatics — any interaction kernel that satisfies Laplace's equation (or more generally is smooth at distance) admits multipole expansions. Variants handle the Helmholtz equation for wave scattering, making the method applicable to acoustics and electromagnetism at a range of frequencies.

Learn more about the related problem of simulating n interacting bodies in our article on n-body simulation, or see how the Barnes-Hut tree relates to closest pair search.

Once you can compute all-pairs interactions in $O(n)$, a vast territory of formerly prohibitive problems becomes routine:

• Astrophysics and cosmology: simulating the formation of galaxies and large-scale structure of the universe requires tracking millions to billions of particles over cosmic time. FMM-based codes such as PKDGRAV3 run on supercomputers and produce the universe-scale maps used to study dark matter and dark energy.
• Molecular dynamics: the dominant cost of simulating proteins, membranes and drug-receptor binding is computing electrostatic forces among millions of charged atoms. AMBER and NAMD use FMM variants to make millisecond-scale protein simulations possible.
• Integral equations: many boundary-value problems in engineering (fluid flow around an aircraft, acoustic scattering from a submarine hull) are solved by discretising the boundary into panels and summing interactions. Without FMM the matrix-vector multiply is $O(n^{2})$; with FMM it is $O(n)$, enabling panels in the millions.
• Kernel methods in machine learning: Gaussian process regression, support vector machines and kernel density estimation all require summing a kernel function over all training points. The FMM can accelerate these sums and is behind fast Gaussian transform implementations used in large-scale statistical computing.
• Radiosity and rendering: computing how light from every patch of a scene reaches every other patch is an all-pairs problem. FMM-inspired hierarchical radiosity methods bring it into $O(n \log n)$ or better.

The FMM is an example of how a mathematical insight — smooth functions can be approximated by a small number of terms — can change the practical complexity class of an entire field of computation.

The Fast Multipole Method is a reminder that algorithmic complexity is not fixed by the problem — it is fixed by how you look at it. The all-pairs force problem looks $O(n^{2})$ until you notice that distant groups can be summarised, at which point $O(n)$ becomes reachable.

Greengard and Rokhlin did not change the laws of physics; they changed the representation. By encoding the effect of a cluster of particles in a compact multipole expansion and translating that expansion hierarchically through a tree, they turned a trillion operations into a million — and made computational science of the cosmos, the cell, and the circuit board possible at modern scales.

The algorithm is solved, its error bounds are proven, and it has been in production in the world's largest supercomputer codes for decades. It stands as a model for what happens when the right mathematical idea meets the right computational structure: not an incremental speedup, but a phase transition in what is computable.