006 — blob structure, centroid count, and five implementations

Same budget for every full fit: at most 12 Lloyd iterations, or stop when the center shift falls under 1e-8 times the mean feature variance. Float32. The host is the L4 box (4 CPUs, 15 GiB RAM). The array-API code is eager PyTorch, CPU or CUDA, with no torch.compile. The request listed PyTorch CPU under the GPU heading; the GPU array-API run is PyTorch CUDA eager.

Scikit-learn 1.7.2 is the stock build. The optimized build is scikit-learn 1.10.dev0 from cakedev0/scikit-learn#23 at 7e26a5c, which also contains the chunked Lloyd kernel from cakedev0/scikit-learn#26. cuML is 26.10.00. cuML has no public k-means++ function, so that row is KMeans(max_iter=0, init="k-means++"), which leaves the initial centers in place and then assigns labels. That assignment is included in the cuML init time.

Seeds are not a shared stream. The array-API draws use NumPy Generator(0). Both scikit-learn builds use random_state=0. cuML uses random_state=0 as well. On the smaller matrices the CPU and CUDA array-API draws stay identical. On the two largest they diverge, because a float rounding difference changes a weighted sample and every later center follows.

Peak memory is an increment. PyTorch reports allocator bytes above the input. Scikit-learn reports the rise in VmHWM on the first call. cuML reports nvidia-smi, whose step is 1 MiB; a 2.0 MiB increment is that step, not a measured workspace.

Cases

Blob data places k random centers, scales them by separation, and adds Gaussian noise of scale std. The ratio separation / std is the blob-iness. noise is standard normal. Lloyd timing always starts from the first k rows, so that duration is not mixed with the initializer.

Quality

Mean squared distance to the assigned center. For well separated blobs this lands on features * std² once every blob owns a center (32 for the mild family, 512 for speed-m).

Init, mean squared distance:

Fit, mean squared distance and iterations actually run:

Where the table shows one inertia repeated across a row, the relative gap is under 0.1%. That happens for tight, mild, noise, tight-d128, noise-d128, k4, tall, and wide. On those, every implementation uses the same small iteration count when the blobs are separated, and the full budget of 12 when the data is noise.

Disclaimer: several gaps are the initializer, and they are not all in one direction. Extending the array-API run on loose from 12 to 40 iterations only moved the inertia from 73.65 to 73.64, so the extra iterations are not going to catch scikit-learn's 53.99. That fit started from a worse k-means++ draw (129.60 against 98.05). The same pattern, with the opposite winner:

• k16: cuML's init is 96.51 against about 62, and its fit stays at 52.76 while the others reach 31.94 in 2 iterations.
• k64 and k128: cuML's init is better and its fit reaches the blob floor near 32, while the array API is still near 38 after 12 iterations.
• mixed: scikit-learn's init is 189 against about 127, and its fit stays at 94.62 while the array API and cuML reach 63.99 in 2 iterations.
• speed-m: cuML's init is 1017 against about 1334, and its fit reaches 511 (the 512 * std² floor) in 2 iterations. The others are still near 670 after 7 to 12 iterations.
• overlap: cuML's fit is 84.63 against 81.77 for the other four, about 3.5% higher, after a full budget of 12 iterations.
• speed-l: the array API and cuML land together (1199 and 1189). Both scikit-learn builds stay near 1360, about 14% higher, after a worse init (about 2700 against 2380).
• speed-xl: cuML reaches 2044, next to the 2048 * std² floor. The array-API CPU fit is 2214 and the CUDA fit is 2398; those two draws diverged. Both scikit-learn builds finish near 2382. The spread from the floor to the CUDA fit is about 17%.

Empty clusters are not the cause of the loose gap: the array-API labeling at iteration 12 used every center. The implementation still does not relocate an empty center, which scikit-learn does. That difference did not appear on the rows where the inertias match.

Duration and memory

Calls below a few tens of milliseconds are real second-call times, but cuML's one-iteration measurement sits on a 40 ms floor on every small shape, so it is not resolved under that floor. The speed-m fit is the first row where several implementations take longer than a second.

k-means++, seconds and peak increment:

One Lloyd iteration from the first k rows:

Full fit:

On speed-xl (131072 samples, 2048 features, 256 clusters, a 1.00 GiB matrix) one Lloyd step is 0.026 s on CUDA, 0.285 s in cuML, 1.855 s for the array API on CPU, and 2.03 s for both scikit-learn builds. The scikit-learn step peaks 2.00 GiB above the input, two extra copies of the matrix. The array-API CPU step peaks at 164 MiB and the CUDA step at 131 MiB. cuML peaks about 1.02 GiB above the memory already in use.

The full fit on that shape is 1.52 s for array-API CUDA, 4.00 s for cuML, 30.3 s for the patched scikit-learn, 49.6 s for array-API CPU, and 123 s for stock scikit-learn. cuML's fit is the one that sits on the blob floor. The patched initializer is the fastest of the CPU inits (24.1 s against 41.8 s for the array API and 118 s for stock scikit-learn) and its draw is the worst of the five.

A PyTorch CPU matmul of this width uses the MKL build with 4 threads. The per-center product is a short N dimension, and the CPU initializer remains slower than the patched scikit-learn kernel. No further rewrite beat that gap inside the array API.