Introduction

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 n2n^{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(n2)O(n^{2}) to O(n)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.

Try It

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

<div class="controls">
  <label>Particles: <span id="nLabel">80</span>
    <input type="range" id="nSlider" min="20" max="400" step="20" value="80">
  </label>
  <label>Grid cells: <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</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 (approx): <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">Each cell in the grid groups nearby particles. Distant cells contribute a single multipole sum — one interaction per (cell, particle) pair instead of one per particle pair.</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; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Fast Multipole Method demo
// Simulates n charged particles; compares cost of brute-force O(n²) vs a
// simple 2-D cell-list FMM approximation O(n · G) where G = grid cells.
// "Interactions" = number of (source, target) pairs evaluated.

const canvas = document.getElementById('c');
const ctx    = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;

let N = 80, G = 4;
let particles = [];

// ---------- helpers ----------
function rand(lo, hi) { return lo + Math.random() * (hi - lo); }
function dist2(a, b)  { const dx=a.x-b.x, dy=a.y-b.y; return dx*dx+dy*dy; }
function dist(a, b)   { return Math.sqrt(dist2(a,b)); }

function spawn() {
  particles = [];
  for (let i = 0; i < N; i++) {
    particles.push({ x: rand(8, W-8), y: rand(8, H-8), q: (Math.random() < .5 ? 1 : -1) });
  }
}

// ---------- brute-force: exact Coulomb force on particle 0 ----------
function bruteForce() {
  // Returns { ops, fx, fy } — force on particle 0 from all others
  let fx = 0, fy = 0, ops = 0;
  const p0 = particles[0];
  for (let j = 1; j < N; j++) {
    const pj = particles[j];
    const d2 = dist2(p0, pj);
    const d  = Math.sqrt(d2) + 1e-6;
    const F  = pj.q / d2;
    fx += F * (p0.x - pj.x) / d;
    fy += F * (p0.y - pj.y) / d;
    ops++;
  }
  return { ops: N * (N - 1), fx, fy }; // report total pairwise ops
}

// ---------- FMM approximation: cell multipole on particle 0 ----------
function fmmApprox() {
  const cellW = W / G, cellH = H / G;

  // Assign particles to cells
  const cells = {};
  for (const p of particles) {
    const ci = Math.min(Math.floor(p.x / cellW), G-1);
    const cj = Math.min(Math.floor(p.y / cellH), G-1);
    const key = ci + ',' + cj;
    if (!cells[key]) cells[key] = { charge: 0, cx: 0, cy: 0, count: 0 };
    cells[key].charge += p.q;
    cells[key].cx     += p.x;
    cells[key].cy     += p.y;
    cells[key].count  += 1;
  }
  // Compute centroid of each cell
  for (const k in cells) {
    const c = cells[k];
    c.cx /= c.count;
    c.cy /= c.count;
  }

  const p0 = particles[0];
  const p0ci = Math.min(Math.floor(p0.x / cellW), G-1);
  const p0cj = Math.min(Math.floor(p0.y / cellH), G-1);

  let fx = 0, fy = 0, ops = 0;
  const threshold = Math.max(cellW, cellH) * 1.5; // "far" criterion

  for (const k in cells) {
    const [ci, cj] = k.split(',').map(Number);
    const c = cells[k];

    if (ci === p0ci && cj === p0cj) {
      // Same cell: interact with individual particles
      for (const pj of particles) {
        const same_ci = Math.min(Math.floor(pj.x / cellW), G-1) === p0ci;
        const same_cj = Math.min(Math.floor(pj.y / cellH), G-1) === p0cj;
        if (!(same_ci && same_cj)) continue;
        if (pj === p0) continue;
        const d2 = dist2(p0, pj);
        const d  = Math.sqrt(d2) + 1e-6;
        const F  = pj.q / d2;
        fx += F * (p0.x - pj.x) / d;
        fy += F * (p0.y - pj.y) / d;
        ops++;
      }
    } else {
      // Far cell: use cell centroid (monopole approximation)
      const dx = p0.x - c.cx, dy = p0.y - c.cy;
      const d2 = dx*dx + dy*dy;
      const d  = Math.sqrt(d2) + 1e-6;
      const F  = c.charge / d2;
      fx += F * dx / d;
      fy += F * dy / d;
      ops++;
    }
  }

  // Total FMM ops: n (classify) + G² (cell sums) + n (near) + G² (far per particle, 1 particle shown)
  // For full simulation: each particle does O(G²) far + O(n/G²) near → O(n) when G ~ sqrt(n)
  const totalFmmOps = N + G * G + (N / (G * G) | 0) * G * G;

  return { ops: totalFmmOps, fx, fy };
}

// ---------- rendering ----------
function draw() {
  ctx.clearRect(0, 0, W, H);

  // Draw grid
  const cellW = W / G, cellH = H / G;
  ctx.strokeStyle = 'rgba(100,130,160,0.25)';
  ctx.lineWidth = 1;
  for (let i = 1; i < G; i++) {
    ctx.beginPath(); ctx.moveTo(i*cellW, 0); ctx.lineTo(i*cellW, H); ctx.stroke();
    ctx.beginPath(); ctx.moveTo(0, i*cellH); ctx.lineTo(W, i*cellH); ctx.stroke();
  }

  // Draw particles
  for (let i = 0; i < particles.length; i++) {
    const p = particles[i];
    if (i === 0) {
      // Target particle: larger, highlighted
      ctx.beginPath(); ctx.arc(p.x, p.y, 7, 0, 2*Math.PI);
      ctx.fillStyle = '#1d3557'; ctx.fill();
      ctx.strokeStyle = '#fff'; ctx.lineWidth = 1.5; ctx.stroke();
    } else {
      ctx.beginPath(); ctx.arc(p.x, p.y, 3.5, 0, 2*Math.PI);
      ctx.fillStyle = p.q > 0 ? '#e63946' : '#2a9d8f';
      ctx.fill();
    }
  }

  // Draw cell centroids for non-home cells
  const p0  = particles[0];
  const p0ci = Math.min(Math.floor(p0.x / cellW), G-1);
  const p0cj = Math.min(Math.floor(p0.y / cellH), G-1);
  const cells = {};
  for (const p of particles) {
    const ci = Math.min(Math.floor(p.x / cellW), G-1);
    const cj = Math.min(Math.floor(p.y / cellH), G-1);
    const key = ci+','+cj;
    if (!cells[key]) cells[key] = { charge:0, cx:0, cy:0, count:0 };
    cells[key].charge += p.q; cells[key].cx += p.x; cells[key].cy += p.y; cells[key].count++;
  }
  for (const k in cells) {
    const [ci,cj] = k.split(',').map(Number);
    const c = cells[k]; c.cx/=c.count; c.cy/=c.count;
    if (ci===p0ci && cj===p0cj) continue;
    ctx.beginPath(); ctx.arc(c.cx, c.cy, 5.5, 0, 2*Math.PI);
    ctx.fillStyle = 'rgba(200,150,20,0.55)'; ctx.fill();
    ctx.strokeStyle = '#b07800'; ctx.lineWidth = 1; ctx.stroke();
  }

  // Label
  ctx.fillStyle = '#1d3557';
  ctx.font = '11px system-ui';
  ctx.fillText('● target particle   ■ cell centroid (FMM source)   ● +charge   ● −charge', 6, H-6);
}

// ---------- update stats ----------
function update() {
  const bf  = bruteForce();
  const fmm = fmmApprox();

  document.getElementById('bfOps').textContent  = bf.ops.toLocaleString();
  document.getElementById('fmmOps').textContent = fmm.ops.toLocaleString();

  // Error: relative difference in force magnitude on particle 0
  const bfMag  = Math.sqrt(bf.fx*bf.fx + bf.fy*bf.fy) + 1e-12;
  const fmmMag = Math.sqrt(fmm.fx*fmm.fx + fmm.fy*fmm.fy) + 1e-12;
  const err = Math.abs(bfMag - fmmMag) / bfMag * 100;
  document.getElementById('errVal').textContent = err.toFixed(1) + '%';

  draw();
}

// ---------- controls ----------
const nSlider = document.getElementById('nSlider');
const gSlider = document.getElementById('gSlider');
nSlider.addEventListener('input', () => {
  N = +nSlider.value;
  document.getElementById('nLabel').textContent = N;
  spawn(); update();
});
gSlider.addEventListener('input', () => {
  G = +gSlider.value;
  document.getElementById('gLabel').textContent = G+'×'+G;
  update();
});
document.getElementById('runBtn').addEventListener('click', () => { spawn(); update(); });

// initial run
spawn();
update();

// Verify: with N=80 G=4, FMM ops must be << N*(N-1)=6320
const bf80  = 80 * 79;
const fmm80 = 80 + 16 + (80/16|0)*16;
console.assert(fmm80 < bf80,
  fmt('FMM ops should be less than brute-force ops: fmm=%fmm80% bf=%bf80%', { fmm80: fmm80, bf80: bf80 }));
console.log(fmt('Demo verified: bf ops=%bf80% fmm ops=%fmm80%', { bf80: bf80, fmm80: fmm80 }));
Result
Console

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 Real Complexity

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(n2)O(n^{2}). For n = 10610^{6} that is 5 × 101110^{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(nlog⁡n)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)O(n) with a constant that depends on the desired accuracy (number of expansion terms p): roughly O(p4n)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.

Where It Matters

Once you can compute all-pairs interactions in O(n)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(n2)O(n^{2}); with FMM it is O(n)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(nlog⁡n)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.

Conclusion

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(n2)O(n^{2}) until you notice that distant groups can be summarised, at which point O(n)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.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/fast-multipole-method/Content licensed under CC BY-NC 4.0.