[BUG] kmeans::fit is not run-to-run reproducible; KmeansFitBatchedTestF.Result/4 flakes ~50% on GB10 and should not gate CI

Related: #2100 (still open), #2256 (changed Result/4's tol 1e-4 -> 1e-5)

Summary

KmeansFitBatchedTests/KmeansFitBatchedTestF.Result/4 fails about half the time
on an NVIDIA GB10 (SM 12.1), on both CUDA 13.3 and 13.4, even with the #2256
change applied. A standalone reproducer (below, no gtest, no batching, no host
data, just the in-memory cuvs::cluster::kmeans::fit) shows the cause is not
batched-vs-in-memory ordering: the in-memory fit alone is not bit-reproducible
between identical runs, and for this test's shape that noise decides whether the
convergence check stops at iteration 11 or iteration 14. The two stopping points
give centroids that differ by ~3%, so the test's 1e-2 centroid comparison fails.

Because the outcome is a coin flip on a legitimate numerical property (both results
are equally good clusterings), we think this test should not be part of a gating
run as written.

Environment

• GPU: NVIDIA GB10 (DGX Spark, aarch64), compute capability 12.1, driver 580.173.02
• CUDA toolkit: 13.3 and 13.4 (both reproduce)
• cuvs: e0f8a4eb (upstream main, 2026-09-17); raft f197acbc
• Built single-arch sm_121a, -O3

Failure rate (--gtest_repeat=50, filter ...KmeansFitBatchedTestF.Result/4)

toolkit	failed / 50
CUDA 13.3	22 (44%)
CUDA 13.4	26 (52%)

(Upstream measured 9/50 on an RTX PRO 6000 with CUDA 13.2 before #2256; #2256 did
not remove the flake on GB10.)

Every failure is at the same coordinate and shows the same two values:

Value of: centroids_match
  Actual: false (actual=0.5148894 != expected=0.49995017 @1,1)

Across all 48 recorded failures the pair is always {0.4999, 0.5148} at @1,1, with
the two values swapping roles between the batched result and the in-memory
reference from run to run. So the in-memory reference itself takes both values.

Root cause (from the reproducer)

Same inputs, same seed, same starting centroids (init = Array), in-memory
kmeans::fit, 40 fits in one process, shape 10000 x 16, k = 10, tol = 1e-5,
max_iter = 20:

run  0: centroid[1,1]=0.49995014 inertia=582834.562500 n_iter=14
run  1: centroid[1,1]=0.49995008 inertia=582834.937500 n_iter=14
run  2: centroid[1,1]=0.51488936 inertia=582841.937500 n_iter=11
run  3: centroid[1,1]=0.51488930 inertia=582841.750000 n_iter=11
...
--- distinct centroid[1,1] (x1000, truncated) over 40 runs:
  499 : 15 runs   (n_iter=14, lower inertia)
  514 : 25 runs   (n_iter=11)

• Two outcomes, split by n_iter (11 vs 14). The final inertia differs by only
  ~1.3e-5 relative (582834.6 vs 582842.0): these are equally good clusterings.
• inertia also varies in its low digits within each outcome, i.e. the
  reduction is not bit-reproducible. This shows up at every size we tried (30
  identical fits per process; ranges over repeated process runs): 1000 x 2, k=3
  gives 1-2 distinct inertia values (4 of 6 runs gave 2); 5000 x 16, k=10 gives
  6-8; 100000 x 2, k=3 gives 10-13.
• Most shapes are stable anyway, because the noise is far too small to change
  the iteration at which convergence trips. 10000 x 16, k=10 happens to sit exactly
  on the tolerance threshold, so the noise flips the stopping iteration (e.g. 11 x 17
  / 14 x 13; split varies per run). Nearby shapes (9000, 8000, 12000 rows, d=8/12, k=8/12) are stable.
  10000 x 2, k=10 also flips, less often (8 x 27 / 12 x 3).

So this is the convergence check (change in a noisy, atomically-reduced inertia
compared with tol) meeting a non-deterministic reduction, not a
GB10-specific miscompute. The test bit-matches centroids after an early stop,
which turns any iteration-count flip into a failure. Lowering tol in #2256 moved the
threshold but cannot make the stopping iteration deterministic.

(Our earlier hypothesis, that the batched path's different summation order versus
the in-memory path is the cause, is not needed: the in-memory path alone flips.)

Reproducer

Self-contained, ~35 lines, no gtest. Usage: ./kmeans_repro N D K [tol].
./kmeans_repro 10000 16 10 shows the iteration-count flip;
./kmeans_repro 100000 2 3 shows bitwise-different inertia between identical fits
(10-13 distinct values out of 30).

#include <cuvs/cluster/kmeans.hpp>
#include <raft/core/device_mdarray.hpp>
#include <raft/core/resource/cuda_stream.hpp>
#include <raft/core/resources.hpp>
#include <raft/random/make_blobs.cuh>
#include <raft/random/rng.cuh>
#include <cstdio>
#include <cstdlib>
#include <map>
int main(int argc, char** argv)
{
  using T = float;
  int n = atoi(argv[1]), d = atoi(argv[2]), k = atoi(argv[3]), reps = 30;
  float tol = argc > 4 ? atof(argv[4]) : 1e-5f;
  raft::resources handle;
  auto s = raft::resource::get_cuda_stream(handle);
  auto X = raft::make_device_matrix<T, int>(handle, n, d);
  auto L = raft::make_device_vector<int, int>(handle, n);
  raft::random::make_blobs<T, int>(X.data_handle(), L.data_handle(), n, d, k, s.get(), true,
                                   nullptr, nullptr, T(1.0), false, T(-10), T(10), 1234ULL);
  auto init = raft::make_device_matrix<T, int>(handle, k, d);
  raft::random::RngState rng(1);
  raft::random::uniform(handle, rng, init.data_handle(), k * d, T(-1), T(1));
  cuvs::cluster::kmeans::params p;
  p.n_clusters = k; p.tol = tol; p.rng_state.seed = 1; p.oversampling_factor = 0;
  p.init = cuvs::cluster::kmeans::params::Array; p.max_iter = 20;
  std::map<int, int> iters; std::map<float, int> inertias;
  for (int r = 0; r < reps; ++r) {
    auto c = raft::make_device_matrix<T, int>(handle, k, d);  // identical start every time
    raft::copy(c.data_handle(), init.data_handle(), k * d, s);
    T inertia = 0; int n_iter = 0;
    cuvs::cluster::kmeans::fit(handle, p, raft::make_const_mdspan(X.view()), std::nullopt,
                               c.view(), raft::make_host_scalar_view<T>(&inertia),
                               raft::make_host_scalar_view<int>(&n_iter));
    iters[n_iter]++; inertias[inertia]++;
  }
  printf("n=%d d=%d k=%d tol=%g  n_iter histogram:", n, d, k, tol);
  for (auto& kv : iters) printf(" %d(x%d)", kv.first, kv.second);
  printf("  distinct inertia values: %zu\n", inertias.size());
}

Example output on GB10 / CUDA 13.4 (30 identical fits each; exact counts vary
from run to run):

n=10000  d=16 k=10 tol=1e-05  n_iter histogram: 11(x24) 14(x6)  distinct inertia values: 11
n=1000   d=2  k=3  tol=1e-05  n_iter histogram: 8(x30)          distinct inertia values: 2
n=100000 d=2  k=3  tol=1e-05  n_iter histogram: 7(x30)          distinct inertia values: 12

Suggested directions

1. Test (smallest change): don't require centroids to match to 1e-2 after a
   tolerance-based early stop. Both outcomes above are equivalent (inertia within
   ~1e-5 relative). Compare inertia / ARI within a tolerance instead, or fix the
   iteration count for the comparison (note tol <= 0 is rejected by fit, so
   this needs a very small positive tol with a matching max_iter).
2. Until then: exclude or quarantine Result/4 from gating runs. It has been
   open as #2100 since May, and retrying once is not enough at a ~50% failure rate
   (a retry loses 1 time in 4).
3. Library (optional): if run-to-run reproducibility of kmeans::fit matters,
   the atomic reductions on this path would need a fixed summation order. #2256
   attributes the numerical differences to raft::reduce_rows_by_key; we did not
   isolate which reduction is responsible.

We have not tested any of these against upstream CI, only on GB10.

---

Assisted-by: Claude

Start with KmeansFitBatchedTests/KmeansFitBatchedTestF.Result/4 and build the provided kmeans_repro program for the 10000 x 16, k=10 case. Inspect how the test compares centroids after tolerance-based early stopping, then determine whether the test should use an equivalent-result criterion or be excluded from gating; done means repeated runs no longer produce a CI-blocking false failure.