force-directed at scale: barnes-hut and multilevelforce-directed at scale: barnes-hut and multilevel

part 1 ended with a problem. every step, each node gets pushed by every other node. that's 10,000 pairs at 100 nodes, a million at 1,000, and a hundred million at 10,000. multiply by a few hundred steps and a big graph takes minutes. too slow to use.

this post covers the two tricks that fix it. the first makes each step cheaper. the second makes the graph smaller.

the idea: far away, things blur togetherthe idea: far away, things blur together

think about a node on the left edge of the drawing, and a tight group of 20 nodes on the far right. does it matter exactly where each of those 20 sits? barely. from far away, they push almost like one big node of weight 20, sitting at their center.

so instead of 20 calculations, do one. that's the whole trick. the hard part is organizing nodes so you can find "tight far-away groups" quickly. for that, you use a quadtree.

the quadtreethe quadtree

a quadtree chops the drawing into boxes:

1. start with one box around everything.
2. if a box has more than one node, split it into four equal boxes.
3. repeat until every box has at most one node.

each box also remembers how many nodes it holds and where their center is. that's all Barnes-Hut needs.

six nodes in a 100 by 100 area:

 +--------------------------------+
 |         2                      |
 |                                |
 |   1                  5    6    |
 |                                |
 |        3                       |
 |                    4           |
 +--------------------------------+

split once. top-left has three nodes, top-right has two, bottom-right has one, bottom-left is empty:

 +---------------+---------------+
 |  1, 2, 3      |  5, 6         |
 |  (3 nodes)    |  (2 nodes)    |
 +---------------+---------------+
 |  (empty)      |  4            |
 |               |               |
 +---------------+---------------+

top-left and top-right still have more than one node, so they split again, until every node has its own box.

the barnes-hut walkthe barnes-hut walk

to find the push on one node, start at the biggest box and ask one question: is this box far enough away to treat as a single node?

• if yes: one calculation, using the box's total weight and center. skip everything inside it.
• if no: open the box and ask the same question about its four smaller boxes.

"far enough" is measured by comparing the box's width s with its distance d from the node. a small box far away has a small s/d. the rule:

if s / d < θ   treat the box as one node
otherwise      look inside it

θ (theta) is your knob.

an example. the node is at (10, 10), and a box of 20 nodes spans (70, 70) to (90, 90). width s = 20, distance to its center d ≈ 99. so s/d ≈ 0.2, which is under θ = 0.5: one calculation instead of 20, with about 4% error.

move the node to (65, 65), right next to the box. now d ≈ 21 and s/d ≈ 0.95. too close, so the walk looks inside.

in code:

function push(node, box) {
  if (box is empty) return 0;
  if (box holds only this node) return 0;

  const d = distance(node, box.center);
  if (box.width / d < theta || box is a single node) {
    return pushFrom(box.center, box.weight, node);   // treat as one
  }
  return box.children.map(child => push(node, child)).sum();  // look inside
}

the θ knobthe θ knob

• θ = 0: nothing is ever far enough, so every node is computed one by one. that's the slow version again.
• θ around 0.5 to 0.9: the sweet spot. usually 10 to 100 times faster, and you can't see the error. d3-force uses 0.9.
• θ very large: even nearby boxes get treated as one node. fast, but the forces point the wrong way near groups, and nodes slide along edges of clusters they should be pushed away from.

try it on 150 nodes. tap any node, then drag θ up from 0:

how much faster?how much faster?

building the quadtree and walking it both take about n log n work per step instead of n². at 10,000 nodes that's roughly 130,000 operations instead of 100,000,000.

the tree is rebuilt every step, since nodes move. that cost is included in the n log n.

two cases where it helps less:

• very dense clusters. inside a packed cluster, nearby boxes are never "far enough", so you're back to near one-by-one work there. Barnes-Hut shines on spread-out graphs.
• nodes almost on top of each other. splitting until two nearly identical points separate can take dozens of levels. d3 handles exactly-identical points with a special case, but near-identical ones still make deep trees.

reading the source: d3-quadtreereading the source: d3-quadtree

d3-quadtree (at commit d01b12a0) is the cleanest real implementation. here's how a point gets added:

function add(tree, x, y, d) {
  if (isNaN(x) || isNaN(y)) return tree; // ignore invalid points

  var parent,
      node = tree._root,
      leaf = {data: d},
      x0 = tree._x0,
      y0 = tree._y0,
      x1 = tree._x1,
      y1 = tree._y1,
      xm, ym, xp, yp, right, bottom, i, j;

  // If the tree is empty, initialize the root as a leaf.
  if (!node) return tree._root = leaf, tree;

  // Find the existing leaf for the new point, or add it.
  while (node.length) {
    if (right = x >= (xm = (x0 + x1) / 2)) x0 = xm; else x1 = xm;
    if (bottom = y >= (ym = (y0 + y1) / 2)) y0 = ym; else y1 = ym;
    if (parent = node, !(node = node[i = bottom << 1 | right])) return parent[i] = leaf, tree;
  }

  // Is the new point exactly coincident with the existing point?
  xp = +tree._x.call(null, node.data);
  yp = +tree._y.call(null, node.data);
  if (x === xp && y === yp) return leaf.next = node, parent ? parent[i] = leaf : tree._root = leaf, tree;

  // Otherwise, split the leaf node until the old and new point are separated.
  do {
    parent = parent ? parent[i] = new Array(4) : tree._root = new Array(4);
    if (right = x >= (xm = (x0 + x1) / 2)) x0 = xm; else x1 = xm;
    if (bottom = y >= (ym = (y0 + y1) / 2)) y0 = ym; else y1 = ym;
  } while ((i = bottom << 1 | right) === (j = (yp >= ym) << 1 | (xp >= xm)));
  return parent[j] = node, parent[i] = leaf, tree;
}

three things to notice:

• a box with children is just an array of 4. the child index is bottom << 1 | right, so 0 is top-left and 3 is bottom-right.
• the while (node.length) loop walks down the tree to find the right box for the new point.
• two points at exactly the same spot aren't split forever. they're chained together in one box with leaf.next.

this is the walk Barnes-Hut uses:

export default function(callback) {
  var quads = [], q, node = this._root, child, x0, y0, x1, y1;
  if (node) quads.push(new Quad(node, this._x0, this._y0, this._x1, this._y1));
  while (q = quads.pop()) {
    if (!callback(node = q.node, x0 = q.x0, y0 = q.y0, x1 = q.x1, y1 = q.y1) && node.length) {
      var xm = (x0 + x1) / 2, ym = (y0 + y1) / 2;
      if (child = node[3]) quads.push(new Quad(child, xm, ym, x1, y1));
      if (child = node[2]) quads.push(new Quad(child, x0, ym, xm, y1));
      if (child = node[1]) quads.push(new Quad(child, xm, y0, x1, ym));
      if (child = node[0]) quads.push(new Quad(child, x0, y0, xm, ym));
    }
  }
  return this;
}

the callback decides whether to look inside a box. return true and the walk skips that box's children: that's "far enough, treat as one". return false and it opens the box.

d3-force passes exactly that kind of callback when computing the push:

function apply(quad, x1, _, x2) {
  if (!quad.value) return true;

  var x = quad.x - node.x,
      y = quad.y - node.y,
      w = x2 - x1,
      l = x * x + y * y;

  // Apply the Barnes-Hut approximation if possible.
  if (w * w / theta2 < l) {
    if (l < distanceMax2) {
      if (x === 0) x = jiggle(random), l += x * x;
      if (y === 0) y = jiggle(random), l += y * y;
      if (l < distanceMin2) l = Math.sqrt(distanceMin2 * l);
      node.vx += x * quad.value * alpha / l;
      node.vy += y * quad.value * alpha / l;
    }
    return true;  // prune: don't descend into this cell's children
  }

  // Otherwise, process points directly.
  else if (quad.length || l >= distanceMax2) return;

  if (quad.data !== node || quad.next) {
    if (x === 0) x = jiggle(random), l += x * x;
    if (y === 0) y = jiggle(random), l += y * y;
    if (l < distanceMin2) l = Math.sqrt(distanceMin2 * l);
  }

  do if (quad.data !== node) {
    w = strengths[quad.data.index] * alpha / l;
    node.vx += x * w;
    node.vy += y * w;
  } while (quad = quad.next);
}

the key line is w * w / theta2 < l. it's the same s/d < θ rule, just squared so it can skip a square root. before this runs, d3 does one pass over the tree to compute each box's total weight and center.

coarsening: solve a small version firstcoarsening: solve a small version first

Barnes-Hut makes each step cheaper. but at a million nodes, even cheap steps add up. the second trick attacks the problem from the other side: make the graph smaller.

1. coarsen. merge pairs of connected nodes into one "super-node". repeat until only a handful are left.
2. lay out the small graph. with a few nodes, force-directed finishes almost instantly.
3. refine. split each super-node back into its parts, placed near where the super-node was. run a few quick steps to tidy up. repeat until you're back to the full graph.

ten nodes, shrunk twice:

original (10 nodes):
  1 - 2 - 3 - 4 - 5
  |           |
  6 - 7   8 - 9 - 10

after one merge (5 super-nodes):
  [12] - [34] - [59]
          |      |
         [67]  [810]

after another (3 super-nodes):
  [1234] - [59810]
               |
             [67]

here it is step by step, all the way down and back up:

which pairs to merge? two common choices: merge each node with its most strongly connected neighbor (Walshaw 2003), or pick a spread-out set of "leader" nodes and merge every other node into a neighboring leader (Hu 2005). Hu's version tends to give more even super-nodes.

why it workswhy it works

a big layout has two jobs: put each cluster in the right general area, and arrange the nodes inside each cluster.

on the full graph, moving a whole cluster across the drawing is painfully slow. every node can only move a little per step, and each one is mostly held by its neighbors. on the shrunken graph, that same cluster is one or two super-nodes, and it moves in a single step.

so the small graph does the big moves cheaply, and the refinement steps handle local detail. it's the same idea as multigrid solvers in numerical math.

in practice, Graphviz's sfdp and OGDF's FM³ both use this plus Barnes-Hut, and lay out a million nodes in seconds. d3-force only does Barnes-Hut. past about 10,000 nodes in the browser, compute positions with sfdp or FM³ and let d3 just draw them.

papers to read nextpapers to read next

Barnes and Hut (1986), "A hierarchical O(N log N) force-calculation algorithm." four pages, written for simulating stars, not graphs. the math is the same because both forces weaken with distance.

Walshaw (2003), "A Multilevel Algorithm for Force-Directed Graph Drawing." the paper that brought multilevel methods to graph layout.

Hu (2005), "Efficient and High Quality Force-Directed Graph Drawing." the most readable of the three. free PDF. if you read one, read this.

code: sfdp is part of Graphviz (gitlab.com/graphviz/graphviz). FM³ is part of OGDF (github.com/ogdf/ogdf). reading sfdp next to Hu 2005 is the fastest way to build your own.

next: part 3 covers the Sugiyama framework, for graphs where everything flows one way.