As rabbit holes usually are, you never expect falling into them. I watched a YouTube video, took a Graph Algorithms college course, and realized I was working on a long-term graph visualization engine.

Why visualize graphs?

Many ideas in graph theory hurt to think about and explain, but are highly intuitive when you see a picture. Think k-core, community detection, and chordal graphs. To add to this, there is something oddly satisfying about creating a picture any non-technical person can see and understand from such an abstract data structure. There is much to be gained from the ability to spot anomalies in a graph at a glance, see what the World Wide Web looks like, or just watch spring simulations converge.

Different classes of graphs call for differently shaped drawings. And different algorithms can produce highly different drawings for the same input graph. Each algorithm usually defines a set of aesthetics it achieves and guarantees. This means the embedding engine should provide many options for layout algorithms, such that the user can correctly decide on which one fits their graph data best.

Force Directed Framework

A graph in itself is a pure mathematical structure with no idea of geometry. Having to go from this structure to a drawing can be a daunting experience through the abstract. Such layout algorithms are deterministic and work on combinatorial ideas in graph theory.

However, it turns out that we can elegantly model they graph as an n-body physics simulation. This can take out all the complicated graph theory operations and replace them with two primatives that define the final layout. Repulsive forces, applied globally on all vertices, and attractive forces, applied only between vertices with an edge between them.

The two functions which define the force magnitudes in turn define the whole layout. And while anything can be used for them, through the years of research a few became standards. Most of the good ones are usually modeled after real-world force functions that are in similar n-body problems.

Noack proposed a general energy model for these equations. In general, force models can be defined by a tuple (a, r), where a and r are the exponents of the distance variable in the attractive and repulsive force functions.

Depending on how the two exponents relate to eachother in terms of growth, leads to very different output layouts. You can get seperated clumps of vertices, or add constants to the force equations to control the edge length.

For example, the Fruchterman-Reingold force equation has the tuple (2, -1), and defines an ideal edge length k.

$$ f_a(d) = \frac{d^2}{k}, \qquad f_r(d) = \frac{k^2}{d} $$

Solving for equilebrium, we see that that is achieved once the distance is equal to the ideal edge length. And in cases with the distance is smaller, repulsion starts increasing quickly, with attraction decreasing. And the opposite happens when the vertices are too far apart.

The layout on top is based on LinLog, a force model designed to seperate hubs. And we can see that happening as the groups of vertices that are close together are tightly connected, and seperated from other groups. Additionally, we can see how in the Fruchterman Reingold layout edge lengths are mostly similar. With an exception to the areas where heavy connections occur, which leads to the attractive forces overpowering.

While the different force models produce vastly different layouts, one shared drawback is present among all of them. That is the O(|V|^2) time complexity that comes with pairwise force calculations. So, I'll show you some significant optimization algorithms I've ran into which help mediate this time complexity.

Barnes-Hut

The repulsive forces are the most expensive part of the force calculations, as attractive forces are only between vertices with an edge between them, while repulsive forces are globally pairwise. Barnes-Hut is an optimization that can be added to force-directed models that takes the time complexity down to O(|V|log(|V|)), while keeping the main structure and loop the same.

Barnes and Hut achieve this by maintaining a spatial index data structure of the plane, commonly a quadtree, wherein distant vertices are approximated as one body when calculating repulsive forces. We construct the quadtree such that leaf nodes only include 1 vertex. For a given vertex V, we calculate its repulsive forces like so, with variables theta a constant, node initially set to the quadtree root, and vPos V's position.

1. If node is a leaf, we calculate the repulsive force normally.
2. We calculate a ratio = node.width / dist, where dist is between the center of mass of the node and vPos
3. If the ratio < theta, we treat the quadtree node as one vertex and calculate repulsive forces against its center of mass.
4. Otherwise, we repeat recursively on each child node.

void ForceAtlas::AccumulateBHRepulsion(const QuadTree::Node *node,
                                       size_t selfIdx, const double *vPos,
                                       double *acc) const {
  // forces are accumlated in acc

  if (!node || node->Mass() == 0.0) // empty
    return;

  if (node->IsLeaf()) {
    // only one vertex in the quadtree node, by construction
      size_t idx = node->PointAt(0); 
      if (idx == selfIdx)
        return;
      const double *uPos = positions_.data() + idx * 2;
      model_->Repulsive(vPos, uPos, mass_[selfIdx], mass_[idx], acc);
    return;
  }

  double comX, comY;
  node->CenterOfMass(&comX, &comY);
  double dx = comX - vPos[0];
  double dy = comY - vPos[1];
  double dist = std::sqrt(dx * dx + dy * dy);
  double ratio = (2.0 * node->HalfSize()) / dist;

  if (ratio < theta_) {
    double com[2] = {comX, comY};
    model_->Repulsive(vPos, com, mass_[selfIdx], node->Mass(), acc);
    return;
  }

  for (size_t q = 0; q < QuadTree::kQuadrantCount; q++)
    AccumulateBHRepulsion(node->Child(static_cast<QuadTree::Quadrant>(q)),
                           selfIdx, vPos, acc);
}

This clever modification to how the repulsive forces are calculated takes algorithms from being unrunnable on large graphs to running smoothly. A 10,000 vertex graph takes 23ms to run a single force iteration without Barnes-Hut, and 4ms with! Additionally, the fact that Barnes-Hut is agnostic of the force model being used makes it essentially a free on/off switch for optimization!

Graph Drawing with Intelligent Placement

In contrast to Barnes-Hut, GRIP is not a plug-and-play optimization to an existing force model. Instead, it defines its own process and produces the graph layout in a multilevel fashion that is designed to accommodate extremely large graphs. GRIP, and other multilevel algorithms, work by dissecting the input graph G into multiple subgraphs G1 -> G2 -> ... -> Gn, working on the layout from the smallest subgraph, where operations are cheap, and propagating information up the chain, repeating until a final layout is achieved.

GRIP uses a vertex subset filtration, where each layer is a subset of the previous, and contains fewer vertices that are also further apart. The core idea is a repeated BFS: pick an unvisited vertex, keep it, and mark everything within some radius as unpickable; this ensures the distance constraint at each layer. The radius shrinks as we approach the full graph; this means later layers are smaller in size but run a deeper BFS, keeping search cost balanced across layers.

Given the nested subset structure, each layer is just a contiguous slice of a single array (MisFiltration), with MisBorders marking where each slice starts. We use the following iterative process to build it, stopping once a layer has vertex count equal to the layout dimension plus one.

bool GRIP::IterMISFiltration(size_t i, BitSet &prevSubset) {
  size_t nvertices = Structure().VertexCapacity();
  BitSet newSubset(nvertices);
  BitSet marked(nvertices);

  size_t radius = size_t{1} << (i - 1);

  for (size_t curr : prevSubset) {
    if (marked.Test(curr)) {
      continue;
    }

    newSubset.Set(curr);
    marked.Set(curr);

    // VerticesWithinRadius performs a BFS from curr with maxDepth 
    // radius, and sets the reached vertices in the marked BitSet.
    VerticesWithinRadius(curr, radius, marked);
  }

  size_t writePos = MisBorderAt(i - 1) - 1;
  for (size_t curr : prevSubset) {
    if (!newSubset.Test(curr)) {
      misFiltration_[writePos--] = curr;
    }
  }

  size_t border = newSubset.Popcount();

  // There are edge cases here, but handling them isn't complex
  bool cont = border == Dim(); 

  misBorder_.push_back(border);
  prevSubset = newSubset;

  return cont;
}

With our MIS filtration created, we can place the last layer, the one with Dim() + 1 vertices, in a simplex. A triangle for 2D, tetrahedron for 3D, and so on. From there, we iterate over the layers in reverse order and, for each layer Vi, we take the following two phase process:

1. This phase is the namesake of the algorithm, intelligent placement. For each vertex V in Vi but not Vi+1, we place V at the barycenter of its K-nearest neighbors in Vi+1, in terms of graph distance.
2. The second phase is local refinement. For ALL the vertices in Vi, we calculate attractive and repulsive forces with its K-nearest neighbors, also in terms of graph distance. We can run the refinement for multiple rounds.

That's it! After running these two phases on all layers, the final layout is produced. The two simple phases, which all perform local updates with no knowledge of global state, give us the following guarantees:

1. The initial position of each vertex will be near its position in the final layout.
2. All global structures and symmetries of the graph will appear in the final layout, even though all updates were local!

This O(|V|) algorithm never stopped blowing my mind, and can process and layout awesome fractal graphs.

Trees

Visualizing trees is an extremely common place graph visualization show up. As apposed to layout of general graphs, tree layouts are mathematically solved and optimized to use the least amount of space possible. While you can use normal force directed methods, there are mathematical visualizers which are deterministic and have concrete gaurantees. In fact, you have probably used such tree visualizers regularly.

The tree command can visualize a subtree of your file system and draw it to the terminal. While no fancy rendering or physics is being done, it gives useful and readable detail about the hierarchical structure of your directory.

% tree -L 2
.
├── CLAUDE.md
├── CMakeLists.txt
├── compile_commands.json
├── include
│   ├── BitSet.hpp
│   ├── BreadthFirst.hpp
│   ├── ConnectedComponents.hpp
│   ├── DepthFirst.hpp
│   ├── EmbeddedGraph.hpp
    ...
├── README.md
├── src
│   ├── BreadthFirst.cpp
│   ├── ConnectedComponents.cpp
│   ├── DepthFirst.cpp
│   ├── EmbeddedGraph.cpp
│   ├── ForceAtlas.cpp
│   ├── Graph.cpp
    ...
└── third-party
    ├── boyerMyrvold
    ├── msf_gif.h
    └── unity

Force-Directed Tree Layouts

The tree command, and other tree layout algorithms, give extremely normalized drawings regardless of the tree shape. This is in part due to the fact that they rely on deterministic mathematical processes to find vertex offsets. However, we can still use our force-directed methods to produce intereseting layouts.

LinLog defines forces that fit Noack's framework as the tuple (0, -1), where 0 means we scale with the logarithm of the distance. It has properties that seperate hubs, vertices with many edges, from eachother in the layout. As well as keeping the smaller, non-hub, vertices close to the hubs they connect to. We can see these properties reflected in the tree layout. Where large branching points having comfortable distance between them, and their leaves being tucked close.

The Fruchterman-Reingold (FR) layout, which uses (2, -1) forces, is very tree-like, moreso than LinLog. This is primarily because the quadratic attraction forces overpower the linear repulsion forces. This leads to much more uniform edge lengths, which gives its branching this tree-like look.

Reingold-Tilford Algorithm

When visualizing user-facing layouts for small trees, for example when showing Kubernetes deployments in a tool like ArgoCD, the layout should be deterministic and similar across trees. We can extend the recursive methods in the tree program to implement an extended version of the Reingold-Tilford (RT) algorithm (the original RT algorithm only supports binary trees) for tree layout. Algorithms like this give us strong gauruntees for the shape of the output, Reingold and Tilford define the following aesthetics:

1. Vertices at the same depth level of the tree should lie along a straight line, and the straight lines defining the levels should be parallel.
2. The children of a vertex should appear in their adjacency-list order, left to right.
3. A parent should be centered over its children.
4. A tree and its mirror image should produce drawings that are reflections of one another; moreover, a subtree should be drawn the same way regardless of where it occurs in the tree.

The last two of which are not achieved by tree. Reingold and Tilford propose an O(|V|) algorithm to achieve all of them, which has become a standard. The y-axis is trivial. It is simply the depth of each node multiplied by some scaling factor. To find x-axis positions, the algorithm calculates the relative x-offsets in each subtree independently. This independence is what ensures aesthetic (4) is achieved, as no global state is required to find what any one subtree will look like.

Visualizing another 500 vertex tree with RT produces an extremely normal and optimally packed layout. However, Since one axis is scaling with depth, and the other is scaling with branching. Our layout is extremely wide, as vertex count scales exponentially with depth.

The algorithm follows a divide-and-conquer approach through a post-order traversal of the tree (children left-to-right then parent), and places each node as a subtree root.

void ReingoldTilford::CalculateOffsets(size_t root, size_t level) {
  size_t degree = graph_.Degree(root);

  // Divide
  for (size_t i = 0; i < degree; i++)
    CalculateOffsets(graph_.Neighbor(root, i), level + 1);

  // Initialization. Includes base case (leaf node)
  InitializeRTSubtreeRoot(root, level);

  // Conquer
  for (size_t i = 0; i < degree; i++)
    CombineSubtreeLeft(root, i);
}

As is usual, the base case (a node with no children) is trivial. It is simply placed with no child offsets to calculate. After that, when non-leaves are reached in the traversal, they are combined based on a contour-walk.

CombineSubtreeLeft(root, child) takes the a subtree and its parent, it calculates the minimum x-coordinate distance between subtree root and its left sibling such the subtree does not overlap with any siblings on its left. It repeats this in a loop until all children are merged. This is the conquer step, and I think the process is best explained via a physical analogy.

First, with your child subtrees each drawn on a piece of paper (which happens in the previous recursive calls that returned, via the same process), cut along the left and right boundaries of each subtree. Each silhouette of a subtree is made up of the left and right contours. Where the contours of a subtree are defined as the leftmost/rightmost vertex at each of its depth layers.

When merging a subtree with its left siblings, begin by overlapping its root with the sibling directly to its left. Then, you walk the right contour of the left subtree and left contour of the right subtree in lockstep, and calculate the x-coordinte difference between the contour vertices at each step. If that difference is less than a constant minimum distance, you add to the seperation the difference. You repeat this process until you reach the bottom of either contour.

The tricky part is that most trees will have unequal depths, so the two contours being compared will usually have one side reach a leaf before the other is done walking its contour. The fix is threading. When a contour walk hits a leaf on the shallow side while the other side still has levels left, it creates a thread, a temporary fake edge, from that leaf down to the next vertex on the deeper contour.

These threads effectievely merge the subtree with its siblings, so that later walks that pass through the leaf will now follow the thread. This will make the shallow subtree's contour behave as if it extended all the way down. These threads are purely a bookkeeping tool, and are usually stored via gross pointer repurposing.

Once all children are merged this way, the parent is simply placed centered based on the calculated relative offsets of its children. These offsets are still relative to each node's own parent though, so a second pass, this time top-down from the root, walks the tree once more accumulating each node's offset with its parent's to produce final absolute coordinates.