k-means clustering over any numeric dataset, parallelised with OpenMP, deterministic under a seed, and tested.
If you came to use it, start at Build. If you came to learn why parallel code is hard, start at The bug — this repository contains a textbook data race that survived six years in public, and the evidence of what it did.
300 identical runs of the 2019 standardisation loop on eight threads. The correct answer is one
number; the racy loop returns a different one almost every time, and usually one that is orders of
magnitude too small. Regenerate with python3 scripts/make_figures.py.
Written in 2019 for the graduate course Matrices Distribuidas during the MSc in Data Science at Universidad Ricardo Palma, Lima, by Alexander Castro, Octavio Palomino and John Edisson Tapias Zarrazola.
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build
ctest --test-dir buildOr without CMake:
gcc -std=c17 -O2 -fopenmp -Iinclude src/kmeans.c src/main.c -lm -o kmeansOpenMP is optional: without it the library builds and runs correctly, serially.
OMP_NUM_THREADS=8 ./kmeans -k 6 -i 50 -s 42 -l data/credit-52918.tsv-k CLUSTERS number of clusters (required)
-i ITERATIONS iteration budget (default 100)
-s SEED seed for k-means++ initialisation (default 1)
-l the first field of each row is a label, not a coordinate
-r use the data as given, without standardising columns
-q print the summary only, without the per-row assignments
The summary reports fit_seconds, timed around the fit alone. Timing the whole process instead
would measure the cost of printing 52 919 lines, which parallelises very differently from k-means.
Input is whitespace-separated numeric rows. Any width, any number of columns.
The 2019 version standardised each column like this:
#pragma omp parallel for
for (i = 0; i < observaciones; i++) {
desviaciones[0] += pow(clientes[i].diasMora - medias[0], 2);
desviaciones[1] += pow(clientes[i].lineaCredito - medias[1], 2);
...
}Four shared accumulators, written by every thread, with no reduction and no atomic. This
is a data race: x += y is a read, an add and a write, and two threads interleaving those steps
lose updates.
Here is what it actually did — same input, same parameters, the standard deviation of one column:
| Threads | desviaciones[1] |
|---|---|
| 1 | 786.540344 ← the correct value |
| 4 | 763.643555 |
| 8 | 763.996704 |
| 8 | 637.565979 |
| 8 | 779.569397 |
Three runs at eight threads, three different answers, none of them right. And because standardisation happens before clustering, every distance in every subsequent iteration was computed on corrupted data. Every result that program ever produced on more than one thread was wrong, and nothing in its output said so.
That table is what the 2019 program printed. The figure at the top of this README is the same bug
isolated in bench/race_demo.c, which runs both versions of the loop side by
side in one process, on the same data, with the same threads:
cc -std=c17 -O2 -fopenmp -Iinclude src/kmeans.c bench/race_demo.c -lm -o race_demo
OMP_NUM_THREADS=8 ./race_demo data/credit-52918.tsv 300 1On the machine these figures were generated on, the racy loop understates the true standard deviation of every column: by roughly 2.5× on three of them, and by up to two orders of magnitude on the fourth. The exact factor changes between runs, which is why no fixed number is quoted here — your machine will print different ones, and that is precisely the finding.
A warning about what you may see. The race is not guaranteed to appear. On a different machine, a smaller dataset or a colder cache, the interleaving may never happen and the racy loop may agree with the correct one. That is the danger, not a reassurance: undefined behaviour that usually works is still undefined behaviour, and it will pick its moment.
The fix is one clause:
#pragma omp parallel for reduction(+ : variance)
for (size_t i = 0; i < n; i++) { ... variance += delta * delta; }reduction gives each thread a private accumulator and merges them safely at the end. It costs
nothing and it is the entire difference between a correct program and a plausible one.
#pragma omp parallel for
for (i = 0; i < observaciones; i++)
clientes[i].grupo = (rand() % numCentroides) + 1;rand() is not thread-safe: it mutates hidden global state. Called from parallel threads it
is undefined behaviour, and in practice threads serialise on a lock inside libc — so this loop
was slower and wrong.
This version uses an explicit PRNG state, called only from serial sections, which is also what makes a seeded run reproducible.
The program could cluster exactly one dataset in the world. The observation type was a hard-coded credit record:
struct Cliente { char id[9]; float diasMora, lineaCredito, anioCastigo, anioUltimoPago; int grupo; };
#define NUM_ATRIBUTOS 4Four columns, fixed names, one domain. It is now an n × d matrix of doubles read from any
whitespace-separated file — which is the difference between coursework and a tool.
A real defect in the shipped data, found by refusing to guess. Row 49 569 of
data/credit-52918.tsv is:
3 N105382 228.00 22132.00 2005 2005
The identifier contains a space. The 2019 reader used fscanf("%s") and silently read 3 as the
id and N105382 as the first numeric field — strtod would have given 0 — misaligning every
column of that row. The row had been quietly corrupt for six years.
The rule here is explicit and documented: with -l, the trailing dimension fields are the
coordinates and everything before them is the label. Identifiers with spaces are ordinary in
real data; guessing is not.
k-means++ instead of random assignment. The original assigned each point to a random cluster.
This version seeds centres with probability proportional to squared distance from the nearest
chosen centre — which carries a proven O(log k) approximation guarantee that random assignment
has never had, at the cost of one extra pass per centre.
A convergence signal. result.converged distinguishes "no assignment changed" from
"the iteration budget ran out". The original could not tell you which had happened.
Not everywhere, and knowing the difference is the point of the course this came from.
| Step | Cost | Parallel? |
|---|---|---|
| Assignment: each point to its nearest centre | O(n · k · d) |
yes — points are independent; this is where the time goes |
| Update: recompute centres | O(n · d) |
yes, with private accumulators merged in a critical section |
| Column statistics | O(n · d) |
yes, with reduction |
Looping over d columns |
O(d) |
no — the original parallelised a four-iteration loop, which costs more in thread setup than it saves |
That last row is worth sitting with: #pragma omp parallel for over a loop of four iterations is
pure overhead. Parallelism is not free, and applying it uniformly is how programs get slower.
Speedup measured against this program's own single-threaded time — the only comparison that means anything. All 52 919 records, 100 iterations, the fastest of five runs at each point.
Twelve threads buy roughly 3×, not 12×, and raising k twentyfold does not improve it. The
reason is visible in the cost table above: with d = 4, the assignment step moves far more memory
than it does arithmetic, so it is bandwidth-bound, and extra threads contend for the same memory
bus rather than adding throughput. The dip at eight threads is reproducible; the machine reports
twelve hardware threads, which is not the same as twelve cores.
None of this is a defect. It is what correct parallel code looks like when you measure it instead of assuming it, and it is the difference between a speedup and a speedup claim.
Running k = 6 on the shipped credit records produces cluster sizes
[4723, 12057, 12576, 992, 1, 22570]. Two of those six clusters are not describing borrowers:
- 992 records give the last payment year as 1900. It is a sentinel meaning never paid, and because it is stored as a number, k-means treats it as a date roughly seven standard deviations in the past and isolates it. The band is clearly visible in the right-hand panel.
- One cluster holds exactly one record — the credit-line outlier, more than 200 standard deviations from the mean. A single row is consuming a sixth of the model's capacity.
The lesson generalises well past this dataset: k-means cannot distinguish a data-entry convention from a phenomenon. Sentinel values and un-winsorised outliers do not announce themselves; they quietly consume clusters, and the only reason they were caught here is that somebody plotted the result. Deciding what to do about them — drop, impute, or model separately — is a domain question the algorithm has no way to answer.
ctest --test-dir build --output-on-failureTen tests, asserting properties rather than outputs: that standardised columns really have zero
mean and unit variance, that a zero-variance column does not become NaN, that two obvious
clusters are recovered, that the same seed gives bit-identical assignments, that more clusters
never fit worse, and that ragged or non-numeric input is rejected instead of guessed.
The determinism test is the one that would have caught the original bug.
pip install matplotlib numpy
python3 scripts/make_figures.pyEvery number in every figure above is measured when you run this, on your machine — nothing is stored. The script also prints the tables behind the plots, so the claims can be checked without reading a picture.
The original is preserved as an annotated git tag rather than a directory, so the repository reads as a tool while the evidence stays one command away:
git checkout coursework-2019| 2019 | Now | Why |
|---|---|---|
shared accumulators in a parallel for |
reduction(+ : ...) |
data race — different answers on every run |
rand() inside a parallel for |
explicit PRNG state, serial | rand() is not thread-safe |
srand(time(NULL)) |
-s SEED |
runs could not be reproduced |
| hard-coded 4-field credit record | n × d matrix of doubles |
the program could cluster one dataset in the world |
fscanf("%s") into char[9] |
tokenised line, documented label rule | misparsed a real row; also a buffer overflow waiting to happen |
| random cluster assignment | k-means++ | a proven approximation guarantee instead of none |
| no convergence signal | converged |
budget exhaustion looked like success |
unchecked malloc |
status codes on every allocation | a failed allocation dereferenced NULL |
float throughout |
double |
52 918 rows accumulate visible error in float |
| parallelised 4-iteration loops | left serial | thread setup cost more than the loop |
Original coursework: Alexander Castro, Octavio Palomino and John Edisson Tapias Zarrazola, Universidad Ricardo Palma, Lima, 2019. First commit 2019-12-07; the source header dates the assignment 2019-11-18.
- Arthur & Vassilvitskii (2007), k-means++: The Advantages of Careful Seeding — the seeding method used here, and its guarantee.
- OpenMP Application Programming Interface, the
reductionclause — the two paragraphs that would have prevented this repository's central bug. - Lloyd (1982), Least squares quantization in PCM — the algorithm itself.
See CITATION.cff, or:
Castro, A., Palomino, O. & Tapias Zarrazola, J. E. parallel-kmeans: k-means clustering in C with OpenMP. Version 1.0.0, 2026. https://github.com/Kemquiros/parallel-kmeans
MIT — see LICENSE.