`label()` is O(n²) — `UnionFind.fast_find`'s path compression never touches the queried node

[matteodellamico]

File: hdbscan/_hdbscan_linkage.pyx — UnionFind.fast_find, used by label()
Type: performance (results are correct; time is quadratic)
Related: the same path-compression defect (bug A below) is also present in
scikit-learn's own sklearn/cluster/_hierarchical_fast.pyx::UnionFind.fast_find. Bug report: scikit-learn/scikit-learn#34626
LLM usage: bug found by using Claude, reviewed manually

Summary

The single-linkage labelling step label() runs in O(n²) time. The root cause is bug (A): the path-compression loop in fast_find never actually shortens the queried node's path, so compression is effectively disabled and fast_find degenerates to a plain (uncompressed) find. On single-linkage inputs — which tend to chain, producing deep merge trees — repeated uncompressed finds sum to O(n²). A second, independent defect (B) makes each individual access slower (Python ndarray indexing instead of the cached C pointer) but is not what causes the quadratic blowup — (A) alone is sufficient for that, as confirmed by the sklearn UnionFind (which doesn't have defect B, see "Related" above) exhibiting the same asymptotic behavior.

Measured on a standalone reproduction of this UnionFind + label() pair, label() alone:

n	label() time
50 000	2 593 ms
200 000	39 565 ms

4× the points → ~15× the time.

The code (upstream _hdbscan_linkage.pyx)

cdef class UnionFind (object):
    cdef np.ndarray parent_arr
    cdef np.ndarray size_arr
    cdef np.intp_t next_label
    cdef np.intp_t *parent        # C pointer into parent_arr.data
    cdef np.intp_t *size

    def __init__(self, N):
        self.parent_arr = -1 * np.ones(2 * N - 1, dtype=np.intp, order='C')
        ...
        self.parent = (<np.intp_t *> self.parent_arr.data)   # cached, and used by union()
        ...

    cdef void union(self, np.intp_t m, np.intp_t n):
        ...
        self.parent[m] = self.next_label     # union() correctly uses the C pointer
        self.parent[n] = self.next_label
        ...

    cdef np.intp_t fast_find(self, np.intp_t n):
        cdef np.intp_t p
        p = n
        while self.parent_arr[n] != -1:      # (B) indexes the ndarray, not self.parent
            n = self.parent_arr[n]
        # label up to the root
        while self.parent_arr[p] != n:       # (A) "path compression"
            p, self.parent_arr[p] = self.parent_arr[p], n
        return n

Root cause

(A) — the path-compression loop never compresses the queried node [primary defect, causes the O(n²)]

The line

p, self.parent_arr[p] = self.parent_arr[p], n

is a chained tuple assignment target1, target2 = rhs1, rhs2. Both Python and Cython evaluate the full right-hand side first — as a tuple of values, using the current value of p — and only then perform the assignments, left to right. So the sequence of events is:

1. Evaluate rhs = (self.parent_arr[p], n) using the old p.
2. Assign target 1: p = rhs[0] — p is rebound now, to its old parent.
3. Assign target 2: self.parent_arr[p] = rhs[1] — but this indexes with the already-updated p from step 2, not the node the loop meant to compress. The write lands on p's old parent, one hop up from where it was intended.

The node originally passed to fast_find(n) is therefore never repointed at the root; only its ancestors shift by one hop per call. Reproducing step 2–3 directly:

>>> parent = [1, 2, 3, -1]          # chain 0 -> 1 -> 2 -> 3(root)
>>> p, n = 0, 3
>>> p, parent[p] = parent[p], n     # one iteration of the "compression" loop
>>> p, parent
(1, [1, 3, 3, -1])                  # node 0 STILL points at 1 — untouched

Every subsequent fast_find(0) re-walks 0 → 1 → 3. Because label() builds a dendrogram where union() allocates a new parent node per merge, single-linkage chaining yields O(n)-deep trees; with compression effectively disabled, the repeated fast_find walks sum to O(n²).

This same tuple-assignment pattern, with the same defect, is used in scikit-learn's own UnionFind.fast_find (sklearn/cluster/_hierarchical_fast.pyx, shared by AgglomerativeClustering's linkage labelling and HDBSCAN's make_single_linkage). sklearn's version stores parent as a typed memoryview rather than a boxed ndarray, so it doesn't have defect (B) below and each access is cheaper — but the compression is just as broken, so it is just as asymptotically quadratic on chaining inputs. See the sklearn bug
report linked at the top of this report.

(B) — indexes the Python ndarray instead of the cached C pointer [secondary, makes each access slower]

fast_find reads self.parent_arr[n] (and [p]), a full ndarray.__getitem__ — bounds-checked, boxes a np.intp scalar, compared as
a Python int — instead of self.parent[n], the C pointer that's already cached in __init__ and is what union() correctly uses. This defect was simply missed when fast_find was written. It multiplies the cost of every access by a large constant factor, but by itself would still be linear-ish per call; it's defect (A) that turns the algorithm quadratic.

Impact

• label() — and therefore any single-linkage path that calls it — is quadratic in the number of points. It is invisible on small inputs but dominates at scale (≈40 s at n=200k in the measurement above; the actual union-find work is
  milliseconds).
• Results are unaffected. fast_find still returns the correct root, so the emitted single-linkage tree is identical; this is purely a performance bug.

Suggested fix

Use the cached C pointer and compress every node on the path to the root, walking to the root first and then compressing in a second pass (this avoids the tuple-assignment aliasing trap entirely):

cdef np.intp_t fast_find(self, np.intp_t n):
    cdef np.intp_t root, nxt
    root = n
    while self.parent[root] != -1:
        root = self.parent[root]
    while self.parent[n] != -1 and self.parent[n] != root:
        nxt = self.parent[n]
        self.parent[n] = root
        n = nxt
    return root

Roots are unchanged, so label()'s output is bit-identical. Verified in a local vendored copy against an independent oracle (hdbscan.hdbscan(..., metric='precomputed') on the same MST) at ARI 1.0, and the whole downstream test
suite. After the fix:

n	label() before	label() after
50 000	2 593 ms	~3 ms
200 000	39 565 ms	~30–180 ms

Notes

• The sibling TreeUnionFind in _hdbscan_tree.pyx is not affected — it stores state in a typed memoryview (C-level indexing) and uses a correct compression loop.
• Reproduction: any single-linkage label() call at n ≳ 10⁴, most visibly on chaining data.