BallTreeBoruvkaAlgorithm returns non-minimal spanning trees even with approx_min_span_tree=False

[matteodellamico]

Disclaimer: I found this bug while studying the Borůvka code for a different
purpose, and it was an LLM (Claude Fable 5) that discovered it and provided a
first analysis. This is the result of me interacting with it and verifying
everything to the best of my ability, but I do not have a deep knowledge of
the code base or the cited papers. Please correct any misunderstandings.

Summary

BallTreeBoruvkaAlgorithm can return a spanning tree that is not a minimum
spanning tree of the mutual reachability graph, even when
approx_min_span_tree=False. The cause is the bound_min tightening term in
the upward bound propagation, which is not a valid bound in the per-component
(Borůvka) setting. KDTreeBoruvkaAlgorithm, whose propagation uses only the
max of the children's bounds, did not fail in any of my tests.

This may also explain the residual ~0.5% discrepancies left unexplained at the
end of #404: my failing datasets contain no duplicate points.

Reproducer

The multiset of edge weights of an MST is unique for a given graph, so a
mismatch in sorted edge weights proves the returned tree is not minimal.

import numpy as np
from scipy.sparse.csgraph import minimum_spanning_tree
from scipy.spatial.distance import pdist, squareform
from sklearn.neighbors import BallTree
from hdbscan._hdbscan_boruvka import BallTreeBoruvkaAlgorithm

MIN_SAMPLES = 5

def exact_mst_weights(X):
    D = squareform(pdist(X))
    # same core-distance convention as the Boruvka code:
    # distance to the MIN_SAMPLES-th neighbor, excluding the point itself
    core = np.sort(D, axis=1)[:, MIN_SAMPLES]
    mr = np.maximum(D, np.maximum(core[:, None], core[None, :]))
    np.fill_diagonal(mr, 0.0)
    return np.sort(minimum_spanning_tree(mr).data)

def boruvka_mst_weights(X, leaf_size):
    alg = BallTreeBoruvkaAlgorithm(
        BallTree(X, leaf_size=leaf_size),
        min_samples=MIN_SAMPLES, metric='euclidean',
        leaf_size=max(3, leaf_size // 3),  # same as hdbscan_.py wrapper
        approx_min_span_tree=False, n_jobs=1)
    return np.sort(alg.spanning_tree()[:, 2])

for seed in range(50):
    X = np.random.default_rng(seed).uniform(size=(500, 3))
    ref = exact_mst_weights(X)
    got = boruvka_mst_weights(X, leaf_size=5)
    bad = np.where(~np.isclose(got, ref, rtol=1e-9))[0]
    if len(bad):
        print(f"seed={seed}: {len(bad)} edge weights differ, "
              f"total weight {got.sum():.6f} vs true MST {ref.sum():.6f}")

Output on my machine:

seed=13: 8 edge weights differ, total weight 71.176836 vs true MST 71.175762
seed=19: 16 edge weights differ, total weight 71.428521 vs true MST 71.425838

The returned edges are all genuine mutual-reachability distances between their
endpoints, and the tree is a valid spanning tree — it is just not minimal:
the differing Borůvka edge weights are always ≥ the exact ones.

leaf_size=5 is used only to make the reproducer compact: smaller leaves
mean a deeper tree with more bound-propagation steps, so failures are dense
enough that 50 seeds suffice. The default leaf_size=40 is affected too,
just more rarely (see table below). The attached zip contains two concrete
failing datasets for the default settings (mixed-density Gaussian blobs,
n=1500, d=8, min_samples=15, leaf_size=40) plus a verify_failing_dataset.py
script: on one of them the Borůvka tree contains an edge of weight 13.7295
where the true MST edge weight at that rank is 12.1116 — a 13% error on a
single merge height, so the errors are not confined to near-ties.

How often

A sweep over 468 configurations per (variant, leaf_size): uniform, Gaussian
blobs with mixed densities, and filament-shaped data; n ∈ {500, 1500},
d ∈ {3, 8}, min_samples ∈ {5, 15}, 26 seeds; all with
approx_min_span_tree=False:

variant	leaf_size	non-minimal trees / runs	worst relative weight inflation
ball	5	18 / 468	5.3e-4
ball	10	18 / 468	5.3e-4
ball	20	2 / 468	4.3e-4
ball	40 (default)	2 / 468	4.3e-4
kd	5–40	0 / 1872	—

The weight error is small, but it perturbs merge heights in the single-linkage
hierarchy, and on one dataset (uniform, n=500, min_samples=5,
min_cluster_size=15, leaf_size=5 — data on which cluster selection is
inherently marginal) the final labels from algorithm='boruvka_balltree'
vs algorithm='generic' differed with an adjusted Rand index of 0.23.

Root cause

In BallTreeBoruvkaAlgorithm.dual_tree_traversal, bounds are propagated up
the tree as

bound_max = max(self.bounds_ptr[left], self.bounds_ptr[right])
bound_min = min(self.bounds_ptr[left] + 2 * (parent_info.radius - left_info.radius),
                self.bounds_ptr[right] + 2 * (parent_info.radius - right_info.radius))
if bound_min > 0:
    new_bound = min(bound_max, bound_min)
else:
    new_bound = bound_max

The bound_min term is the recursive child term of the tighter
k-nearest-neighbor pruning bound of Curtin et al. 2013 ("Tree-Independent
Dual-Tree Algorithms", one of the two papers cited in this file's header):

B(N_q) = min{ …, min_{N_c ∈ C_q} (B(N_c) + 2(λ(N_q) − λ(N_c))), … }

(and the leaf-level new_lower_bound + 2 * radius matches their B₂ term).
That bound is justified in kNN because each candidate distance D_p is
anchored at a query point p inside the node: p's candidates are
physically near p, so nearby query points can inherit them through the
triangle inequality. In Borůvka, candidate_distance belongs to a
component, and a component's current candidate pair can be anywhere in
the dataset — there is nothing inside the child node to anchor that
inequality to. (Curtin et al. §7 do suggest reusing B(N_q) "with k = 1"
for dual-tree Borůvka; as far as I can tell the transfer argument does not
survive the per-component anchoring.)

Concretely: suppose the left child holds points of component C, which has
just found a tiny candidate somewhere else in the dataset
(bounds[left] small), while the right child holds points of component D
whose candidate was reset to infinity by a merge (bounds[right] = ∞).
The propagated bound bounds[left] + 2 * (R_parent - r_left) then asserts
that no point in the parent — including D's points — needs to look beyond
that radius. But D's true nearest outgoing edge can be farther than that,
and with mutual-reachability weights max(d(p,q), core(p), core(q))
arbitrarily farther, since core distances are not bounded by any node's
radius — consistent with the default-leaf-size failures occurring on data
with widely varying densities. D then misses its true minimum outgoing
edge under this node and instead accepts a heavier edge found under a node
with a looser bound: the result is a valid but non-minimal spanning tree,
silently.

A short case analysis shows how little room there is for anything tighter
in this setting. If a node's points span at least two components, every
point has a foreign point inside the node, so
max(2 * R_node, c_max(node)) — with c_max the largest core distance in
the node, needed because the weights are mutual reachabilities — is a
valid bound. If the node's points all belong to one component, that
component's nearest outgoing edge is a global property of the point set
that no local geometry can bound, and the max of the children's bounds is
all there is. Either way, no sound bound is ever a function of a sibling's
candidate distances, which is exactly what bound_min is.

Notably, March, Ram & Gray (KDD 2010) — the original dual-tree Borůvka
paper, whose Theorem 4.1 is the correctness proof for this algorithm —
use the pure max bound only: d(Q) = max_{q∈Q} d(C_q) at leaves and
d(Q) = max{d(Q.left), d(Q.right)} for internal nodes, and the proof
relies on exactly that.

Ablation confirms this in both directions, on a 10-seed sweep (840 runs)
that produces 6 ball-tree failures with pristine code:

• replacing the propagation with new_bound = bound_max (as the kd-tree
  variant already does) → 0 failures;
• keeping bound_min in the propagation but removing the leaf-level
  min(new_upper_bound, new_lower_bound + 2 * radius) term → the same
  6 failures return identically (so the leaf-level term appears benign;
  only the propagation term is implicated). This asymmetry matches the
  case analysis above: the leaf-level term satisfies
  new_lower_bound + 2r ≥ 2r, so it never undercuts the diameter bound
  and is sound in plain distances (only mutual-reachability cores could
  break it, and since core distances are 1-Lipschitz they vary by at most
  2r within a leaf, capping any potential damage at about the leaf
  diameter), whereas the propagation term undercuts the diameter bound
  2 * R_parent already in plain distances whenever
  bounds[child] < 2 * r_child. A targeted adversarial search against
  the leaf-level term alone (~9,600 exactness checks, including
  hand-crafted near-counterexample geometry and duplicate-heavy data)
  found no failure.

KDTreeBoruvkaAlgorithm propagates with plain
max(bounds[left], bounds[right]) — March et al.'s proven bound — which is
consistent with it never failing in these tests.

References: Curtin et al. 2013, arXiv:1304.4327;
March, Ram & Gray, KDD 2010.

Suggested fix

Two changes, one per term:

1. Remove the bound_min propagation term unconditionally
   (new_bound = bound_max) — i.e. revert to the bound March et al.
   actually proved correct, which is also what the kd-tree variant
   already does. It is demonstrably unsound and benchmarks as buying
   nothing: timing BallTreeBoruvkaAlgorithm.spanning_tree()
   (leaf_size=40 tree / 13 internal, min_samples=5, best of 3 runs) on
   uniform and mixed-density blob data at n ∈ {20000, 50000},
   d ∈ {3, 10}, in both approx modes, the patched build was within noise
   of pristine on every configuration (worst case +3.5%: 1.42 s → 1.47 s
   on uniform n=50000 d=3, exact mode). Even in approx mode there is no
   reason to keep it.

2. Gate the leaf-level min-term on approx_min_span_tree=True.
   This term is different: removing it as well costs a consistent ~5–15%
   on the same benchmark (e.g. uniform n=50000 d=3 exact:
   1.42 s → 1.65 s), so it genuinely buys pruning — but it has no
   soundness proof, only the empirical evidence and damage cap above.
   approx_min_span_tree=False is an explicit request for a guaranteed
   minimum spanning tree, and a guarantee mode should not include a
   bound that is merely believed safe; a user who selects the
   non-default exact mode should be fine paying ~10% for the guarantee.
   Under approx_min_span_tree=True the term is entirely in the spirit
   of the flag's documented meaning ("take shortcuts and only
   approximate the min spanning tree"). The same gating logically
   applies to the kd-tree variant's leaf term (also unproven, and
   computed in rdist space where the triangle-inequality reasoning is
   even less clear), though no failure has ever been observed there.

If tighter provably sound pruning is wanted in exact mode, the case
analysis above gives a replacement for both sites: for nodes whose
points span more than one component, use
min(bound_max, max(2 * R_node, c_max(node))) — per-node
single-component status is already tracked in component_of_node, and
c_max (the largest core distance among a node's points) is fixed once
core distances are computed, so it can be precomputed per node.

Environment

• hdbscan master (692c2a5, v0.8.44), built from source
• python 3.14.4, numpy 2.5.0, scipy 1.18.0, scikit-learn 1.9.0, Linux x86-64`

boruvka-nonminimal-mst-repro.zip

`

[lmcinnes]

I'm on vacation right now, so it might be a while before I get to this.

[Abubakar-rashid]

i want to work on this issue, please assign this to me