K-means Clustering
The previous chapter on Gaussian anomaly detection asked “does this observation look like the ones I have already seen?”. K-means clustering asks a different unsupervised question: “if I had to sort these observations into k groups, what would the most natural grouping be?” Neither question requires labelled training data. Both are workhorses of exploratory data analysis, customer segmentation, image compression, document organization, and preprocessing for downstream supervised models.
Dear reader, one of my 62 US patents covers K-means technology: System to label K-means clusters with human understandable labels.
In this chapter we build a K-means clustering library from scratch in Gerbil Scheme. The library supports two initialization strategies, multiple random restarts, an inertia-based scoring criterion for picking the best restart, and both fit and predict operations. We demonstrate it on a synthetic dataset of three Gaussian blobs in the plane, so that the “correct answer” is obvious to the eye and easy to check against the library’s output.
Theoretical Background
Given a dataset of 
points in 
-dimensional space, K-means partitions the points into 
clusters, where each cluster is represented by a single point called its centroid. The algorithm assigns each observation to the cluster whose centroid is nearest (by squared Euclidean distance) and then moves each centroid to the mean of its assigned points. These two steps are repeated until the centroids stop moving.
Formally, K-means seeks the assignment of points 
to clusters 
that minimizes the within-cluster sum of squared distances, often called the inertia or WCSS:
where 
is the centroid of cluster . The iterative algorithm is guaranteed to decrease 
at every step (both the assignment step and the update step can only lower it) and converges to a local minimum. It is not guaranteed to find the global minimum, which is why we run several restarts and keep the best one.
The classic K-means training loop looks like this:
- Initialize centroids somewhere reasonable.
- Assign each point to the nearest centroid.
- Update each centroid to be the mean of its assigned points.
- Repeat steps 2 and 3 until no centroid moves by more than a small tolerance, or a maximum iteration count is reached.
- Report the final assignment and the inertia.
Two well-known pitfalls of K-means are worth handling explicitly. First, a purely random initialization can start with two centroids near each other and end up in a poor local minimum. The k-means++ initialization scheme, due to Arthur and Vassilvitskii (2007), addresses this by picking the first centroid uniformly at random and each subsequent centroid with probability proportional to the squared distance from the nearest already-chosen centroid, biasing the initial centroids to be spread out. Second, an “empty cluster” can appear if no point is nearest to some centroid after the assignment step. Our implementation handles this by keeping the previous centroid in that slot.
The Dataset
For the demo we generate synthetic data so that ground truth is known. Each point is a 2-D vector sampled from one of three Gaussian distributions with these true centers and a shared standard deviation of 
:
1 Blob A: 60 points around (0.0, 0.0)
2 Blob B: 60 points around (4.0, 0.0)
3 Blob C: 60 points around (2.0, 3.5)
The blobs are well separated (they are several standard deviations apart), so a correct implementation should recover the three centers to within a few hundredths in each coordinate and place exactly 60 points in each cluster. Sample points from one blob might look like:
1 #(-0.041 0.194)
2 #( 0.221 -0.310)
3 #( 0.077 0.052)
Vectors of fixed length are the natural representation for feature points in Gerbil Scheme: they give constant-time indexed access and are the same shape used by the anomaly detection module in the previous chapter, keeping the code style consistent.
Project Structure
The project directory source_code/Kmeans contains:
| File | Description |
|---|---|
kmeans-lib.ss |
K-means library: distance helpers, initializers, fit, predict, accessors. |
test.ss |
Demo that generates three Gaussian blobs and clusters them. |
Makefile |
Single test target that invokes gxi test.ss. |
gerbil.pkg |
Package declaration (kmeans). |
README.md |
Short quick-reference for the library. |
The K-means Library: kmeans-lib.ss
This module is self-contained. Its only imports are :std/format for formatted output and :gerbil/gambit for runtime primitives such as random-real. Everything the library needs to fit and predict is exported so that a demo or downstream application can reach in at whatever level of granularity is convenient.
1 ;; File: kmeans-lib.ss
2 ;; K-means clustering library
3 (import :std/format
4 :gerbil/gambit)
5
6 (export euclidean-distance
7 squared-distance
8 assign-clusters
9 update-centroids
10 compute-inertia
11 initialize-centroids-random
12 initialize-centroids-plusplus
13 kmeans-fit
14 kmeans-predict
15 centroids-of
16 labels-of
17 inertia-of
18 iterations-of
19 format-decimal)
Decimal Formatting Helpers
Gerbil’s default number printing produces long decimal expansions like 0.03497348274. For readable reports we reuse the same format-decimal idiom used in the anomaly detection chapter, which rounds to a fixed number of places and pads with zeros:
1 ;; Find index of character in string (helper for format-decimal)
2 (def (string-index str char)
3 (let ((len (string-length str)))
4 (let loop ((i 0))
5 (cond
6 ((= i len) #f)
7 ((char=? (string-ref str i) char) i)
8 (else (loop (+ i 1)))))))
9
10 ;; Format floating point numbers to a fixed number of decimal places
11 (def (format-decimal x places)
12 (let* ((multiplier (expt 10 places))
13 (rounded (/ (round (* x multiplier)) multiplier))
14 (str (number->string (exact->inexact rounded))))
15 (let ((dot-pos (string-index str #\.)))
16 (let ((formatted-str
17 (if dot-pos
18 (let* ((frac-len (- (string-length str) dot-pos 1))
19 (diff (- places frac-len)))
20 (if (> diff 0)
21 (string-append str (make-string diff #\0))
22 (substring str 0 (+ dot-pos places 1))))
23 (string-append str "." (make-string places #\0)))))
24 (if (char=? (string-ref formatted-str 0) #\.)
25 (string-append "0" formatted-str)
26 formatted-str)))))
The call to exact->inexact guarantees we always get a decimal string rather than a rational literal like "7/2" for numbers that happen to be exact.
Squared Distance and Euclidean Distance
K-means only ever compares distances, so the square root in Euclidean distance is redundant during clustering. We split the two operations so that the assignment and inertia loops can call the cheaper squared-distance while user code that wants a real distance value can call euclidean-distance:
1 ;; Squared Euclidean distance between two feature vectors
2 (def (squared-distance a b)
3 (let ((n (vector-length a))
4 (sum-val 0.0))
5 (do ((i 0 (+ i 1)))
6 ((= i n) sum-val)
7 (let ((d (- (vector-ref a i) (vector-ref b i))))
8 (set! sum-val (+ sum-val (* d d)))))))
9
10 ;; Euclidean distance
11 (def (euclidean-distance a b)
12 (sqrt (squared-distance a b)))
squared-distance walks the two vectors coordinate by coordinate, accumulates the squared differences, and returns the running total. This yields 
, which is 
.
Nearest-Centroid Lookup
The assignment step reduces to “for each point, find the index of the closest centroid”. Rather than allocating a list of distances, we sweep once through the centroid vector and remember the best index found so far:
1 ;; Find index of nearest centroid to point x
2 (def (nearest-centroid-index x centroids)
3 (let ((k (vector-length centroids))
4 (best-i 0)
5 (best-d 1.0e30))
6 (do ((i 0 (+ i 1)))
7 ((= i k) best-i)
8 (let ((d (squared-distance x (vector-ref centroids i))))
9 (when (< d best-d)
10 (set! best-d d)
11 (set! best-i i))))))
The initial best-d of 1.0e30 is a sentinel value that is guaranteed to be replaced on the first iteration of the loop.
The Assignment Step
assign-clusters is a thin wrapper that calls nearest-centroid-index for every point and stores the result in a label vector:
1 ;; Assign each point to its nearest centroid.
2 ;; points is a list of vectors; returns a vector of integer labels.
3 (def (assign-clusters points centroids)
4 (let* ((n (length points))
5 (labels (make-vector n 0)))
6 (let loop ((ps points) (i 0))
7 (if (null? ps)
8 labels
9 (begin
10 (vector-set! labels i (nearest-centroid-index (car ps) centroids))
11 (loop (cdr ps) (+ i 1)))))))
Points come in as a list because it is convenient for the demo to append blobs together with the list primitive, but labels come out as a vector because we index into them by position immediately afterwards.
The Update Step
The update step computes, for each cluster, the mean of its assigned points and installs that as the new centroid. If a cluster receives no points during this pass, we hold onto the previous centroid so the algorithm can continue rather than crashing on a divide-by-zero:
1 ;; Recompute centroids as the mean of assigned points.
2 ;; If a cluster is empty, keep the previous centroid.
3 (def (update-centroids points labels k num-features prev-centroids)
4 (let ((sums (make-vector k #f))
5 (counts (make-vector k 0))
6 (new-centroids (make-vector k #f)))
7 (do ((c 0 (+ c 1)))
8 ((= c k))
9 (vector-set! sums c (make-vector num-features 0.0)))
10 (let loop ((ps points) (i 0))
11 (unless (null? ps)
12 (let* ((label (vector-ref labels i))
13 (bucket (vector-ref sums label))
14 (pt (car ps)))
15 (do ((f 0 (+ f 1)))
16 ((= f num-features))
17 (vector-set! bucket f (+ (vector-ref bucket f) (vector-ref pt f))))
18 (vector-set! counts label (+ (vector-ref counts label) 1)))
19 (loop (cdr ps) (+ i 1))))
20 (do ((c 0 (+ c 1)))
21 ((= c k) new-centroids)
22 (let ((cnt (vector-ref counts c)))
23 (if (= cnt 0)
24 (vector-set! new-centroids c (vector-copy (vector-ref prev-centroids c)))
25 (let ((mean (make-vector num-features 0.0))
26 (bucket (vector-ref sums c)))
27 (do ((f 0 (+ f 1)))
28 ((= f num-features))
29 (vector-set! mean f (/ (vector-ref bucket f) cnt)))
30 (vector-set! new-centroids c mean)))))))
The two-pass structure (accumulate then divide) is the standard streaming way to compute a mean without materializing the sublists of points per cluster.
Inertia and Convergence
Inertia is simply the sum of the squared distances from every point to its assigned centroid. We use it both as the scoring criterion between random restarts and as an informative quality metric to report to the user:
1 ;; Within-cluster sum of squared distances (inertia).
2 (def (compute-inertia points labels centroids)
3 (let ((sum-val 0.0))
4 (let loop ((ps points) (i 0))
5 (if (null? ps)
6 sum-val
7 (let* ((label (vector-ref labels i))
8 (d (squared-distance (car ps) (vector-ref centroids label))))
9 (set! sum-val (+ sum-val d))
10 (loop (cdr ps) (+ i 1)))))))
Convergence is decided by comparing every coordinate of every centroid against its previous value. If all of them have moved less than the tolerance value tol, we consider the run converged:
1 ;; True if two centroid sets are within tol of each other on every coordinate.
2 (def (centroids-converged? old-centroids new-centroids tol)
3 (let ((k (vector-length old-centroids))
4 (converged #t))
5 (do ((c 0 (+ c 1)))
6 ((or (not converged) (= c k)) converged)
7 (let* ((a (vector-ref old-centroids c))
8 (b (vector-ref new-centroids c))
9 (nf (vector-length a)))
10 (do ((f 0 (+ f 1)))
11 ((or (not converged) (= f nf)))
12 (when (> (abs (- (vector-ref a f) (vector-ref b f))) tol)
13 (set! converged #f)))))))
Both loops short-circuit via the converged flag as soon as a large enough difference is found, so most non-converged runs finish this check very quickly.
Random Initialization
The simpler of the two initializers picks distinct data points uniformly at random and uses their coordinates as the starting centroids. A hash table of already-used indices keeps the picks distinct:
1 ;; Uniform random selection of k distinct starting points.
2 (def (initialize-centroids-random points k)
3 (let* ((point-vec (list->vector points))
4 (n (vector-length point-vec))
5 (chosen (make-vector k #f))
6 (picked-count 0)
7 (used (make-hash-table)))
8 (let loop ()
9 (when (< picked-count k)
10 (let ((idx (random-integer n)))
11 (unless (hash-get used idx)
12 (hash-put! used idx #t)
13 (vector-set! chosen picked-count (vector-copy (vector-ref point-vec idx)))
14 (set! picked-count (+ picked-count 1))))
15 (loop)))
16 chosen))
vector-copy is important here because the algorithm mutates the centroid vectors during the update step, and we do not want those mutations to accidentally corrupt the original data points.
K-means++ Initialization
The k-means++ scheme picks the first centroid uniformly, then for each subsequent centroid computes the squared distance from every point to its nearest already-chosen centroid, treats those distances as an unnormalized probability distribution, and samples one point from it:
1 ;; K-means++ initialization: pick the first centroid uniformly, then each
2 ;; subsequent centroid with probability proportional to its squared distance
3 ;; from the nearest already-chosen centroid.
4 (def (initialize-centroids-plusplus points k)
5 (let* ((point-vec (list->vector points))
6 (n (vector-length point-vec))
7 (centroids (make-vector k #f)))
8 (vector-set! centroids 0
9 (vector-copy (vector-ref point-vec (random-integer n))))
10 (do ((c 1 (+ c 1)))
11 ((= c k) centroids)
12 (let ((distances (make-vector n 0.0))
13 (total 0.0))
14 (do ((i 0 (+ i 1)))
15 ((= i n))
16 (let ((best 1.0e30))
17 (do ((j 0 (+ j 1)))
18 ((= j c))
19 (let ((d (squared-distance (vector-ref point-vec i)
20 (vector-ref centroids j))))
21 (when (< d best) (set! best d))))
22 (vector-set! distances i best)
23 (set! total (+ total best))))
24 (if (= total 0.0)
25 (vector-set! centroids c
26 (vector-copy (vector-ref point-vec (random-integer n))))
27 (let ((threshold (* (random-real) total))
28 (running 0.0)
29 (picked #f))
30 (do ((i 0 (+ i 1)))
31 ((or picked (= i n)))
32 (set! running (+ running (vector-ref distances i)))
33 (when (>= running threshold)
34 (vector-set! centroids c
35 (vector-copy (vector-ref point-vec i)))
36 (set! picked #t)))
37 (unless picked
38 (vector-set! centroids c
39 (vector-copy (vector-ref point-vec (- n 1)))))))))))
The cumulative-sum sampling trick avoids explicitly building a normalized probability vector: draw a uniform threshold from 
, walk the distances until the running sum exceeds the threshold, and pick that point. Points far from the existing centroids contribute more mass and are more likely to be picked, which spreads the initial centroids out.
A Single K-means Run
kmeans-single-run performs one execution of the assign / update loop for a given initialization. It stops as soon as the centroids stop moving or the iteration cap is reached, and returns the centroids, labels, inertia, and iteration count for that run:
1 ;; Run a single k-means pass with the given initial centroids.
2 ;; Returns (centroids labels inertia iterations).
3 (def (kmeans-single-run points k num-features max-iters tol init-centroids)
4 (let ((centroids init-centroids)
5 (labels (make-vector (length points) 0))
6 (iter 0)
7 (done #f))
8 (let loop ()
9 (unless (or done (>= iter max-iters))
10 (set! iter (+ iter 1))
11 (set! labels (assign-clusters points centroids))
12 (let ((new-centroids
13 (update-centroids points labels k num-features centroids)))
14 (when (centroids-converged? centroids new-centroids tol)
15 (set! done #t))
16 (set! centroids new-centroids))
17 (loop)))
18 (list centroids labels (compute-inertia points labels centroids) iter)))
The Public Fit Function
kmeans-fit is the top-level entry point. It handles keyword arguments, runs n-init restarts, keeps the best one by inertia, and packages the result in a hash table:
1 ;; Fit a k-means model.
2 ;; points : list of feature vectors (all same length)
3 ;; k : number of clusters
4 ;; max-iters: : max iterations per run (default 100)
5 ;; tol: : convergence tolerance on centroid movement (default 1e-6)
6 ;; n-init: : number of random restarts; best-inertia run wins (default 10)
7 ;; init: : 'kmeans++ (default) or 'random
8 ;; Returns a hash table with keys: centroids, labels, inertia, iterations, k, num-features.
9 (def (kmeans-fit points k
10 max-iters: (max-iters 100)
11 tol: (tol 1.0e-6)
12 n-init: (n-init 10)
13 init: (init 'kmeans++))
14 (when (null? points)
15 (error "kmeans-fit: empty point list"))
16 (let* ((num-features (vector-length (car points)))
17 (best-inertia 1.0e30)
18 (best-centroids #f)
19 (best-labels #f)
20 (best-iters 0))
21 (do ((run 0 (+ run 1)))
22 ((= run n-init))
23 (let* ((init-c (case init
24 ((random) (initialize-centroids-random points k))
25 (else (initialize-centroids-plusplus points k))))
26 (result (kmeans-single-run points k num-features
27 max-iters tol init-c))
28 (cents (car result))
29 (labs (cadr result))
30 (inert (caddr result))
31 (its (cadddr result)))
32 (when (< inert best-inertia)
33 (set! best-inertia inert)
34 (set! best-centroids cents)
35 (set! best-labels labs)
36 (set! best-iters its))))
37 (let ((model (make-hash-table)))
38 (hash-put! model 'centroids best-centroids)
39 (hash-put! model 'labels best-labels)
40 (hash-put! model 'inertia best-inertia)
41 (hash-put! model 'iterations best-iters)
42 (hash-put! model 'k k)
43 (hash-put! model 'num-features num-features)
44 model)))
Keyword arguments follow Gerbil’s key: (name default) convention. Notice how case dispatches on the init symbol so that 'random selects the uniform initializer and any other value (including 'kmeans++) selects the k-means++ initializer.
Predict and Accessors
The remaining public API is small. kmeans-predict returns the cluster index for a new point, and four one-line accessors pull individual fields out of the trained model hash table:
1 ;; Predict the cluster label for a single point vector.
2 (def (kmeans-predict model x)
3 (nearest-centroid-index x (hash-ref model 'centroids)))
4
5 ;; Convenience accessors
6 (def (centroids-of model) (hash-ref model 'centroids))
7 (def (labels-of model) (hash-ref model 'labels))
8 (def (inertia-of model) (hash-ref model 'inertia))
9 (def (iterations-of model) (hash-ref model 'iterations))
The Demo Program: test.ss
The demo generates the three Gaussian blobs described earlier, fits the model, prints diagnostics, and classifies a few probe points. It also runs a second fit with random initialization to prove the option is wired up.
Generating Synthetic Data
Gerbil ships with a uniform random-real primitive but not a Gaussian one, so we implement the Box-Muller transform in a few lines. Given two independent uniform samples 
, the transform returns a sample from 
:
The make-blob procedure calls random-normal twice per point to produce a 2-D sample centered at (cx, cy) with the requested standard deviation:
1 ;; Box-Muller transform: draw a single sample from N(0, 1).
2 (def (random-normal)
3 (let ((u1 (max 1.0e-10 (random-real)))
4 (u2 (random-real)))
5 (* (sqrt (* -2.0 (log u1)))
6 (cos (* 2.0 (acos -1.0) u2)))))
7
8 ;; Generate `n` points around (cx, cy) with the given standard deviation.
9 (def (make-blob n cx cy std)
10 (let loop ((i 0) (acc '()))
11 (if (= i n)
12 acc
13 (loop (+ i 1)
14 (cons (vector (+ cx (* std (random-normal)))
15 (+ cy (* std (random-normal))))
16 acc)))))
The max 1.0e-10 clamp on 
prevents the extremely unlikely case where random-real returns exactly 
, which would make 
undefined.
Reporting Helpers
Three small helpers count cluster sizes, sum a vector, and pretty-print the centroids and cluster sizes:
1 ;; Count how many points fall in each cluster.
2 (def (cluster-sizes labels k)
3 (let ((counts (make-vector k 0))
4 (n (vector-length labels)))
5 (do ((i 0 (+ i 1)))
6 ((= i n) counts)
7 (let ((c (vector-ref labels i)))
8 (vector-set! counts c (+ (vector-ref counts c) 1))))))
9
10 ;; Given cluster sizes, return their sum (must equal total point count).
11 (def (sum-vec v)
12 (let ((n (vector-length v))
13 (s 0))
14 (do ((i 0 (+ i 1)))
15 ((= i n) s)
16 (set! s (+ s (vector-ref v i))))))
17
18 (def (print-centroids centroids)
19 (let ((k (vector-length centroids)))
20 (do ((c 0 (+ c 1)))
21 ((= c k))
22 (let ((v (vector-ref centroids c)))
23 (displayln
24 (format " centroid ~a: (~a, ~a)"
25 c
26 (format-decimal (vector-ref v 0) 3)
27 (format-decimal (vector-ref v 1) 3)))))))
28
29 (def (print-cluster-sizes counts)
30 (let ((k (vector-length counts)))
31 (do ((c 0 (+ c 1)))
32 ((= c k))
33 (displayln (format " cluster ~a: ~a points"
34 c (vector-ref counts c))))))
The Main Entry Point
The main procedure ties everything together: generate the blobs, fit, print centroids and cluster sizes, run a few probe predictions, then repeat the fit with random init to compare inertias:
1 (def (main)
2 (random-source-randomize! default-random-source)
3 (displayln "=== K-means clustering demo ===\n")
4
5 ;; Three well-separated blobs.
6 (let* ((blob-a (make-blob 60 0.0 0.0 0.35))
7 (blob-b (make-blob 60 4.0 0.0 0.35))
8 (blob-c (make-blob 60 2.0 3.5 0.35))
9 (points (append blob-a blob-b blob-c))
10 (k 3))
11 (displayln (format "Generated ~a points across 3 true blobs." (length points)))
12 (displayln (format "True centers: (0,0), (4,0), (2,3.5)\n"))
13
14 ;; Fit with k-means++ initialization.
15 (let* ((model (kmeans-fit points k n-init: 10 max-iters: 100)))
16 (displayln (format "Converged in ~a iterations (best of 10 restarts)."
17 (iterations-of model)))
18 (displayln (format "Final inertia (WCSS): ~a"
19 (format-decimal (inertia-of model) 4)))
20 (displayln "Learned centroids:")
21 (print-centroids (centroids-of model))
22 (displayln "Cluster sizes:")
23 (print-cluster-sizes (cluster-sizes (labels-of model) k))
24
25 ;; Sanity checks.
26 (unless (= (sum-vec (cluster-sizes (labels-of model) k)) (length points))
27 (error "cluster sizes must sum to point count"))
28 (unless (= (vector-length (centroids-of model)) k)
29 (error "expected k centroids"))
30 (unless (< (inertia-of model) 100.0)
31 (error "inertia unexpectedly large; blobs should cluster tightly"))
32
33 ;; Predict a few fresh points near each true center.
34 (displayln "\nPredictions for new points:")
35 (for-each
36 (lambda (probe)
37 (let* ((label (kmeans-predict model probe)))
38 (displayln
39 (format " point (~a, ~a) -> cluster ~a"
40 (format-decimal (vector-ref probe 0) 2)
41 (format-decimal (vector-ref probe 1) 2)
42 label))))
43 (list (vector 0.1 0.1)
44 (vector 3.9 0.1)
45 (vector 2.0 3.4)
46 (vector 2.0 1.5))))
47
48 ;; Compare against random init to demonstrate the option is wired up.
49 (displayln "\n--- Sanity check: random init also converges ---")
50 (let ((model2 (kmeans-fit points k n-init: 10 init: 'random)))
51 (displayln (format "random-init inertia: ~a (iters: ~a)"
52 (format-decimal (inertia-of model2) 4)
53 (iterations-of model2))))
54
55 (displayln "\n=== Demo complete ===")))
56
57 (main)
The final probe 
is deliberately near the geometric center of the three blobs. It should end up assigned to whichever cluster happens to be nearest, which depends on the random seed but is always a defensible choice.
Running the Demo
The Makefile is minimal:
1 test:
2 gxi test.ss
Run either make test or invoke the interpreter directly with gxi test.ss. A typical run produces output similar to the following. Because the blob generation and initialization both use randomness, the exact numbers will vary from run to run, and the cluster indices (0, 1, 2) may be permuted:
1 $ make test
2 gxi test.ss
3 === K-means clustering demo ===
4
5 Generated 180 points across 3 true blobs.
6 True centers: (0,0), (4,0), (2,3.5)
7
8 Converged in 2 iterations (best of 10 restarts).
9 Final inertia (WCSS): 39.1594
10 Learned centroids:
11 centroid 0: (1.983, 3.497)
12 centroid 1: (4.056, 0.010)
13 centroid 2: (-0.012, 0.069)
14 Cluster sizes:
15 cluster 0: 60 points
16 cluster 1: 60 points
17 cluster 2: 60 points
18
19 Predictions for new points:
20 point (0.10, 0.10) -> cluster 2
21 point (3.90, 0.10) -> cluster 1
22 point (2.00, 3.40) -> cluster 0
23 point (2.00, 1.50) -> cluster 0
24
25 --- Sanity check: random init also converges ---
26 random-init inertia: 39.1594 (iters: 3)
27
28 === Demo complete ===
Interpreting the Results
The output reports several quantities that are worth understanding in the context of the algorithm.
- Converged in 2 iterations: with k-means++ initialization and well-separated blobs, the initial centroids are already very close to the true means. Two passes are enough to move them into place and detect that they have stopped moving. Poorly initialized runs, or runs on noisier data, would need many more.
- Final inertia = 39.1594: this is the total squared distance from every point to its assigned centroid. Recall that each blob has standard deviation in two dimensions. The expected inertia per point is roughly
, so for 180 points we expect an inertia around 
. The measured value 
is right in that ballpark, confirming the centroids are essentially at the true blob means.
- Learned centroids:
, 
, and 
. Compare these to the true centers 
, 
, and 
. Each learned coordinate is within about 
of the truth, which is well within the noise level of 
divided by 
(the standard error of the mean of 60 samples).
- Cluster sizes: all three clusters contain exactly 60 points, meaning K-means separated the blobs perfectly. Not a single point crossed a decision boundary.
- Predictions for new points: the probe points near
, 
, and 
are placed in the corresponding true clusters. The ambiguous probe at is assigned to whichever cluster centroid it happens to be closest to. In the run above that was cluster 0 (the top blob at ), which was closer than the two bottom-row centroids.
- Random-init inertia identical to k-means++ inertia: on such a clean dataset, both initializers find the same global minimum. On messier data, k-means++ typically finds a lower inertia than uniform random init, especially at low
.
The identical inertia between the two fitting runs is a useful check that both initialization codepaths are exercising the same core algorithm rather than one being subtly broken.
Wrap Up
We built a general-purpose K-means clustering library in Gerbil Scheme in around 200 lines of code. The library covers the essentials: k-means++ and uniform initialization, an assign-update main loop with tolerance-based convergence, empty-cluster handling, multiple random restarts scored by inertia, and both fit and predict operations.
A few observations transfer directly to any dataset you might want to cluster:
- Vectors of fixed length are the right representation for feature points in Gerbil. Constant-time indexing keeps the tight inner loops (distance and mean computation) fast.
- Squared distance is enough for K-means. Only compute the square root when reporting to a user.
- Multiple random restarts ( of
is a good default) plus inertia-based selection is the simplest general defense against poor local minima.
- Wrap the trained model in a hash table with named keys. That gives you a stable API surface (
centroids-of,labels-of, …) even as internal representations evolve.
The one big caveat with K-means is that you have to pick up front. On real-world data you rarely know the true number of clusters. The most common heuristic is the “elbow method”: run kmeans-fit for a range of values, plot the inertia as a function of , and pick the at which the curve bends from steep to shallow. That is a natural extension exercise for readers, along with several other useful additions covered in the practice problems below.
Optional Practice Problems
Elbow-method plot data: Write a procedure
(inertia-curve points k-min k-max)that fits K-means for each from 
to 
and returns a list of 
pairs. Print the pairs and identify the where the drop in inertia flattens. Verify that for the three-blob dataset the elbow is at 
.Silhouette score: The silhouette score is a standard cluster-quality metric that does not require ground truth. For a point
in cluster 
with mean intra-cluster distance 
and mean distance 
to points in the nearest other cluster, the silhouette is 
. Implement (silhouette-score points labels)and report the average silhouette on the three-blob dataset. A well-clustered dataset should score close to
.Mini-batch K-means: The full K-means pass touches every point on every iteration, which becomes expensive on large datasets. Implement
(kmeans-fit-minibatch points k batch-size max-iters)that on each iteration samples
points at random, assigns them to their nearest centroid, and moves each centroid by a small step toward the mean of its assigned batch points. Compare inertia and runtime to the full-batch fit.Save and load a trained model: Implement
(save-kmeans model path)and(load-kmeans path)that serialize the trained centroids,, and 
to a plain-text file and read them back later. Verify that (kmeans-predict (load-kmeans path) x)returns the same label as the in-memory model for the same input
.Cluster real data: The classic Iris dataset has 150 points in 4-D space and 3 known species labels. Download
iris.data, load it into a list of feature vectors (dropping the species label into a separate list), fitkmeans-fitwith, and report how well the discovered clusters align with the true species using a confusion matrix. Because K-means assigns cluster indices arbitrarily, you will need to find the best mapping between cluster indices and true labels before scoring.