Euclidean MST, in depth: why the tree hides inside the Delaunay triangulation, dense Prim and dual-tree Borůvka in higher dimensions, single-linkage clustering and real-data failure modes

Why the tree lives inside the Delaunay triangulation

Two facts about any minimum spanning tree do all the work. The cut property: the lightest edge crossing any partition of the vertices belongs to some MST. The cycle property: the heaviest edge on any cycle belongs to no MST (with distinct weights). Geometry lets us apply the cycle property without enumerating cycles.

Take an EMST edge uv and draw the circle that has uv as its diameter. Suppose a third point w sat inside it. The angle uwv would then exceed 90 degrees, so both |uw| and |vw| are shorter than |uv|. Deleting uv splits the tree into two parts and w lies in one of them, so either uw or vw reconnects the parts more cheaply, which contradicts minimality. Therefore every EMST edge has an empty diametral circle: it is a Gabriel edge. A slightly stronger version of the same argument shows that no w can be closer to both endpoints than they are to each other, which makes it a relative neighbourhood graph edge. And an empty diametral circle is a special case of the empty circumcircle that defines Delaunay edges.

So, for points in general position, EMST ⊆ RNG ⊆ Gabriel ⊆ Delaunay. In the plane the Delaunay triangulation has at most 3n-6 edges and is computed in O(n log n) time, so the recipe is: triangulate, then run Kruskal or Prim on those O(n) edges. The whole EMST costs O(n log n), which is optimal, because sorting reduces to it. The Delaunay triangulation article covers how the triangulation itself is built; here we only consume it.

One more consequence is useful for checking output. Two EMST edges meeting at a point must make an angle of at least 60 degrees, or the third side of the triangle would be shorter than the longer of the two. So in the plane no vertex has degree above 6, and with points in general position some EMST has maximum degree 5. A vertex of degree 9 in your output means a bug.

Choosing an algorithm by size and dimension

The Delaunay shortcut stops working as the dimension grows. In d dimensions a Delaunay triangulation can contain on the order of n^⌈d/2⌉ simplices, and even the typical case is costly beyond about 3 or 4 dimensions. Embedding vectors, sensor features and image descriptors live in tens or hundreds of dimensions, so a different strategy is needed there. Pick the algorithm from n and d:

The second row deserves more respect than it gets. Prim's algorithm on the implicit complete graph never stores the distance matrix: it keeps one array holding each outside point's distance to the growing tree and updates it with one vectorised distance row per step. That is n² distance evaluations but only O(n) memory, it is exact in any dimension, and for 20,000 points it finishes in seconds. The fourth row gives a spanning forest that may not be the exact EMST: if the true MST edge between two clusters is not among either endpoint's k nearest neighbours, it is missing. Check the number of connected components before you trust it.

Tested code: Delaunay filtering and dense Prim

Here are the two exact methods, using NumPy and SciPy. The planar version deduplicates points, triangulates, builds a sparse graph from the Delaunay neighbour lists, and lets scipy.sparse.csgraph.minimum_spanning_tree do Kruskal. The duplicate and collinear handling is there for reasons explained under failure modes.

import numpy as np
from scipy.spatial import Delaunay, QhullError
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import minimum_spanning_tree, connected_components


def emst_2d(points):
    """Euclidean MST of 2-D points via Delaunay filtering. Returns (i, j, length) edges."""
    pts = np.asarray(points, dtype=float)
    uniq, first, inverse = np.unique(pts, axis=0, return_index=True, return_inverse=True)
    inverse = inverse.ravel()
    if len(uniq) < 2:
        return [(int(first[0]), k, 0.0) for k in range(len(pts)) if k != first[0]]
    if len(uniq) == 2:
        edges = [(0, 1)]
    else:
        try:
            tri = Delaunay(uniq)
            indptr, nbrs = tri.vertex_neighbor_vertices
            edges = [(i, j) for i in range(len(uniq))
                     for j in nbrs[indptr[i]:indptr[i + 1]] if i < j]
        except QhullError:                       # every point on one line
            order = np.lexsort(uniq.T[::-1])
            edges = list(zip(order[:-1], order[1:]))
    rows, cols = map(np.array, zip(*edges))
    w = np.linalg.norm(uniq[rows] - uniq[cols], axis=1)
    g = coo_matrix((w, (rows, cols)), shape=(len(uniq), len(uniq))).tocsr()
    t = minimum_spanning_tree(g).tocoo()
    assert connected_components(t, directed=False)[0] == 1, "forest, not tree"
    out = [(int(first[i]), int(first[j]), float(d)) for i, j, d in zip(t.row, t.col, t.data)]
    # each duplicate hangs off its first copy at length 0
    out += [(int(first[inverse[k]]), k, 0.0) for k in range(len(pts)) if first[inverse[k]] != k]
    return out


def emst_dense(X):
    """Exact EMST in any dimension: Prim on the implicit complete graph, O(n) memory."""
    X = np.asarray(X, dtype=float)
    n = len(X)
    in_tree = np.zeros(n, dtype=bool)
    best = np.full(n, np.inf)
    parent = np.full(n, -1)
    best[0] = 0.0
    edges = []
    for _ in range(n):
        u = int(np.argmin(np.where(in_tree, np.inf, best)))
        in_tree[u] = True
        if parent[u] >= 0:
            edges.append((int(parent[u]), u, float(best[u])))
        d = np.linalg.norm(X - X[u], axis=1)     # one row of the distance matrix
        closer = (~in_tree) & (d < best)
        best[closer] = d[closer]
        parent[closer] = u
    return edges

Use the two against each other as a test oracle: on a few thousand random points with a planted duplicate, the total lengths should agree to floating-point tolerance and both should return exactly n-1 edges. Keep that test in CI; it catches the forest bug below.

Dual-tree Borůvka for large point sets

For large n in moderate dimension, the standard exact method is dual-tree Borůvka, described by March, Ram and Gray at KDD 2010 and shipped in mlpack as its EMST tool. It combines two ideas. Borůvka's algorithm proceeds in rounds: every current component finds its single cheapest edge to any other component, all those edges are added at once (the cut property guarantees each is safe), and the number of components at least halves, so there are at most log₂ n rounds. The expensive part, "nearest point outside my component", is a nearest-neighbour query, and a space-partitioning tree answers many of those together.

The dual-tree search walks a query node and a reference node of the same kd-tree simultaneously and prunes whole pairs of boxes:

def boruvka_round(Q, R):                    # Q, R: nodes of the same kd-tree
    if Q.component is not None and Q.component == R.component:
        return                              # every pair already connected: prune
    if box_distance(Q, R) > Q.bound:
        return                              # cannot beat any query's current best
    if Q.is_leaf and R.is_leaf:
        for q in Q.points:
            for r in R.points:
                if comp[q] != comp[r] and dist(q, r) < best[comp[q]].d:
                    best[comp[q]] = Edge(q, r, dist(q, r))
        Q.bound = max(best[comp[q]].d for q in Q.points)
        return
    for Qc, Rc in children_pairs(Q, R, closest_first=True):
        boruvka_round(Qc, Rc)
    Q.bound = max(child.bound for child in Q.children)

while number_of_components > 1:
    best = {c: Edge(None, None, inf) for c in components}
    reset_bounds(tree)
    boruvka_round(tree.root, tree.root)
    for e in best.values():
        union_find.union(e.u, e.v, record=e)    # union-find drops duplicate edges
    relabel(tree)    # a node whose points share one component gets that label

The component labels on tree nodes are what make it fast: late in the run, large regions are a single component and the first test prunes them instantly. Ties need care, because two components can pick the same edge, or two equal-length edges that would close a cycle; routing every addition through a union-find structure and breaking ties by point index keeps the result a tree. Like every kd-tree method it degrades as dimension rises; the kd-tree deep dive explains why pruning stops working past roughly 20 dimensions, and that is the point at which dense Prim or an approximate kNN graph becomes the better choice. In Python, the hdbscan package exposes this family through algorithm='boruvka_kdtree'.

Worked example and single-linkage clustering

Take the six points in the diagram: A(0,0), B(2,0), C(2,1), D(5,1), E(5,3), F(1,4). Sorted, the shortest distances are BC = 1, AB = 2, DE = 2, AC = √5 ≈ 2.24, CD = 3, then BD and CF, both √10 ≈ 3.16. Kruskal walks the list with a union-find:

1. BC = 1: join B and C.
2. AB = 2: join A to {B, C}.
3. DE = 2: join D and E.
4. AC = 2.24: A and C are already connected, so skip it. This is the cycle property: AC is the heaviest edge of triangle ABC.
5. CD = 3: joins {A, B, C} with {D, E}.
6. BD = 3.16: B and D are already connected, so skip it. CF = 3.16: F joins. That is five edges for six points, so stop.

The total length is 1 + 2 + 2 + 3 + 3.16 = 11.16. The tree also gives single-linkage clusterings for free. Removing the longest edge, CF, gives two clusters, {F} and the rest; removing CD as well gives three, {A, B, C}, {D, E} and {F}. That is exactly what agglomerative clustering with single linkage would produce, at the cost of one sort. HDBSCAN uses the same idea with a twist: it builds the MST under mutual reachability distance, max(core_k(a), core_k(b), d(a, b)), so that sparse noise points cannot form cheap bridges between dense regions, then condenses the resulting hierarchy. Compare it with the density-based approach in DBSCAN.

Failure modes on real data

The first three were reproduced by running the code above with SciPy; the rest come from general practice with point data.

• Duplicate points produce a forest. Qhull does not triangulate a repeated point; SciPy lists it in Delaunay.coplanar and it gets no neighbours. Feed the neighbour lists straight to the MST and that point is isolated, so you get n-2 edges and two components. GPS fixes and rounded coordinates make duplicates common. Deduplicate first, as the code does, and assert one component at the end.
• Zero-length edges vanish. minimum_spanning_tree treats an explicit zero in a sparse matrix as a missing edge, so two coincident points joined by a 0.0 weight come back disconnected. Deduplicate, or add a tiny epsilon to every weight.
• Collinear input crashes the triangulation. If every point lies on one line, Delaunay raises QhullError. The EMST is then just the points in sorted order, which is the fallback branch above. Nearly collinear data (a road, a scan line) is triangulated but with slivers; the result is still correct.
• Ties make the tree non-unique. Grids and integer coordinates create many equal distances, so two correct implementations can return different trees with the same total. Compare totals in tests, never edge lists, and do not build downstream logic on which of two equal edges was chosen.
• Unscaled features. In more than two dimensions "Euclidean" is only meaningful if the axes share units. A feature measured in metres beside one in millimetres decides the whole tree. Standardise first, or the clustering you read off the MST describes your unit choices.
• Quadratic memory by accident. Calling scipy.spatial.distance.pdist to build the complete graph needs 8 bytes per pair: about 40 GB for 100,000 points. Dense Prim does the same work in O(n) memory.
• Approximate kNN graphs silently disconnect. Two well-separated clusters may share no k-nearest-neighbour edge. Count components and, if there are several, join them with exact nearest-pair queries between components.

Trade-offs

If the tree feeds clustering, exactness of the longest edges matters most, because those are the ones you cut. An approximate method that perturbs short intra-cluster edges is harmless; one that misses a bridge edge changes the answer.

What to do next

1. Run both functions above on your own data and assert they agree on total length and return n-1 edges.
2. Deduplicate points before any triangulation, and log how many duplicates you removed.
3. Pick the method from the table: Delaunay for 2-D, dense Prim below about 50k points, dual-tree Borůvka or a repaired kNN graph above that.
4. Standardise features before computing distances in more than two dimensions.
5. For clustering, plot the sorted MST edge lengths; a sharp jump marks a natural cut.
6. Review Kruskal in depth if the cycle-property reasoning above felt fast.