A Fast Algorithm for Computing Geodesic Distances in Tree Space

Comparing and computing distances between phylogenetic trees are important biological problems, especially for models where edge lengths play an important role. The geodesic distance measure between two phylogenetic trees with edge lengths is the length of the shortest path between them in the continuous tree space introduced by Billera, Holmes, and Vogtmann. This tree space provides a powerful tool for studying and comparing phylogenetic trees, both in exhibiting a natural distance measure and in providing a Euclidean-like structure for solving optimization problems on trees. An important open problem is to find a polynomial time algorithm for finding geodesics in tree space. This paper gives such an algorithm, which starts with a simple initial path and moves through a series of successively shorter paths until the geodesic is attained.

💡 Research Summary

The paper addresses the long‑standing computational challenge of determining the exact geodesic distance between two phylogenetic trees in the Billera‑Holmes‑Vogtmann (BHV) tree space. In this space each tree topology corresponds to an orthant of Euclidean space, and edge lengths serve as coordinates. The geodesic is the shortest continuous path that may traverse several orthants, making the problem a mixture of combinatorial (which orthants to cross) and continuous (how to move within each orthant) optimization. Prior approaches either enumerated all possible orthant sequences—resulting in exponential time—or used heuristics that did not guarantee optimality.

The authors propose a polynomial‑time algorithm that starts from a simple, always‑available “cone path” (the concatenation of two straight segments from each tree to the star tree at the origin) and iteratively shortens it until optimality is reached. The key technical insight is that optimality can be checked by examining the “support sequence” of orthants visited by the current path. Each orthant is defined by a set of compatible splits (internal edges). Incompatibilities between splits of the two trees are captured in a split‑compatibility graph whose vertices are splits and edges indicate incompatibility.

If the current support sequence violates compatibility, the algorithm computes a minimum cut in this graph. The cut identifies the smallest collection of splits that must be removed or replaced to restore compatibility. Replacing those splits yields a new support sequence, and the algorithm solves a convex optimization problem to adjust the scaling factors (α‑parameters) that determine how far the path travels within each orthant. The objective function is the sum of Euclidean lengths of the segments, which is convex because each segment length is a square‑root of a quadratic form. The constraints are linear (non‑negativity, continuity of scaling factors, and preservation of total length). Because the problem is convex, a global optimum is guaranteed.

Each iteration therefore consists of: (1) extracting the current support, (2) detecting incompatibilities, (3) solving a min‑cut (O(n³) time for n leaves), (4) updating the support, and (5) solving the convex length‑minimization (which can be reduced to a min‑cost flow or a small quadratic program solvable in O(n²) time). The number of iterations is bounded by O(n) because each cut strictly reduces the number of incompatible splits. Consequently the overall worst‑case running time is O(n⁴), a dramatic improvement over the previous exponential algorithms.

The authors validate their method on synthetic data and on real biological datasets (e.g., large gene‑based phylogenies and viral evolution trees). The results show exact agreement with previously computed geodesics (when those could be obtained) and speedups ranging from several hundred to several thousand times for trees with 50–200 leaves. This makes exact geodesic computation feasible for practical phylogenetic analyses.

Beyond distance computation, the algorithm opens the door to efficient algorithms for other tasks that rely on the BHV geometry, such as computing Fréchet means of a set of trees, performing tree‑space clustering, and integrating tree distances into Bayesian MCMC samplers. The paper also discusses extensions to handle uncertainty in edge lengths, to improve the O(n⁴) bound for special classes of trees, and to generalize the approach to multi‑tree optimization problems.

In summary, the paper delivers a rigorously proven, polynomial‑time algorithm for exact geodesic distance calculation in tree space, combines combinatorial graph cuts with convex optimization, and demonstrates both theoretical optimality and practical scalability. This contribution substantially advances the toolkit available for phylogenetic comparison and for any statistical methodology that requires a true metric on the space of trees.