Every site in the middle of its own cell

Scatter some points at random and draw the diagram of who is nearest to what. The cells come out wildly uneven — a few enormous ones where the scatter left a hole, a huddle of slivers where three points landed almost on top of each other. That is what randomness looks like, and it looks like a mistake.

Now apply one rule. Every point moves to the centre of mass of its own cell. Then redraw the diagram, because the cells have changed, and do it again.

The rule is due to Stuart Lloyd, who wrote it down in 1957 inside Bell Labs and did not publish it for twenty-five years. It has been rediscovered repeatedly under other names, which is usually a sign that a rule is the obvious thing to try. What is not obvious is that it works, that it can be proved to work in two lines, and that the thing it converges to is a honeycomb.

The rule is a descent, and the quantity it descends #

Attach a number to any arrangement of sites: for every point of the square, take the squared distance to its nearest site, and add all of those up. Write it as

$$E = \sum_i \int_{C_i} |x - s_i|^2 \, dA,$$

where $C_i$ is the cell of site $s_i$. It is a measure of how badly the sites cover the square — small when every point of the square is close to some site, large when there is a region far from all of them.

Lloyd’s rule lowers $E$, and the argument splits it into two halves that each lower it for a different reason.

The first half holds the cells fixed and moves the sites. For a fixed region, which point minimises the mean squared distance to it? The centroid, and the reason is an identity a page of algebra shorter than it looks. For any point $p$ and any region $C$ with area $A$ and centroid $c$,

$$\int_C |x - p|^2 dA = \int_C |x - c|^2 dA + A|p - c|^2.$$

The cross term vanishes because $c$ is the centroid — that is what being the centroid means. So the integral is a constant plus something proportional to the squared distance from $p$ to $c$, and it is smallest when $p$ is $c$.

That identity is the parallel-axis theorem, which arrives in mechanics as a statement about moments of inertia and here as a statement about a cost. It is the same theorem, and the reason it is the same is that both quantities are $\int |x - p|^2$ over a region — mechanics calls it the moment about $p$ and this problem calls it the cost of representing the region by $p$. The mechanical reading gives the result an intuition worth having: a body spins most easily about its centre of mass, and a region is best represented by its centre of mass, because the two sentences are one calculation.

The second half holds the sites fixed and redraws the cells. After the sites have moved, some points of the square are no longer nearest to the site whose cell they are in. Reassigning them to their actual nearest site can only reduce the contribution of those points, because the contribution is the distance to the site they are assigned to and the new assignment is by definition the smallest such distance.

So each round is two moves, each of which lowers $E$ or leaves it alone, and the whole rule is a descent on a quantity bounded below by zero. That is the entire proof that it cannot run away.

Where the fall actually happens #

Descent arguments are cheap and they are usually uninformative, because knowing a quantity falls says nothing about how fast. Here the shape of the fall is the practically important fact and it is not what the proof suggests.

Almost all of the improvement is bought in the first two rounds. What happens afterwards is a slow, geometric crawl: each round multiplies the remaining distance to the fixed point by a factor slightly under one, so the arrows shorten steadily and never quite reach zero in exact arithmetic.

The reason for the split is that the two rounds are doing different jobs. The first round destroys the clumps, which are where nearly all of the excess cost lives — a pair of sites almost on top of each other contributes two slivers and one enormous neighbouring cell, and one round pulls them apart. What is left afterwards is an arrangement that is already roughly even and is being refined, and refining an even arrangement is a matter of adjusting each cell against its neighbours, which is a local negotiation that propagates one cell per round.

That shape has a practical consequence and it is the reason the rule is used rather than admired. In every application below, the first two rounds are the ones anyone runs, and running two hundred rounds instead of five buys a fraction of a percent. A theorem guaranteeing eventual convergence is not what makes an algorithm useful; the guarantee here is a licence to stop early.

The fixed point, and what it looks like #

The rule stops moving exactly when every site is already at the centroid of its own cell. A diagram with that property is called centroidal, and it is not an arrangement anybody would think to write down.

What is startling is what those arrangements look like. From any starting scatter, the settled cells are close to equal in area and overwhelmingly six-sided.

Nothing in the rule mentions six. The rule says centre of mass, which is a statement about integration, and the diagram is defined by nearest, which is a statement about distance. A hexagon appears because it is the answer to a question neither of them asked: among tilings of the plane by equal-area cells, which one has the smallest mean squared distance from a cell’s points to its centre? The conjectured answer is the regular hexagon, which is Gersho’s conjecture — proved in two dimensions by Fejes Tóth, open in three and above.

So the iteration is not producing hexagons because hexagons are built into it. It is descending a cost whose minimiser happens to be a hexagonal arrangement, and the picture is a piece of evidence about that minimiser rather than a demonstration of it. The same shape turns up for the same reason wherever equal regions are asked to be compact — it is the isoperimetric question asked of a tiling rather than of a single curve, and a circle is what each cell would like to be, with the hexagon as the closest a tiling can come.

The whole thing, on a line #

Everything above happens in one dimension too, and in one dimension it can be carried out by hand, which is worth doing once because it shows exactly where the difficulty in higher dimensions comes from.

Put three sites on the interval from nought to one. The cells are the intervals cut by the midpoints between consecutive sites, and the centroid of an interval is its midpoint. So a round of the rule is: compute the two boundaries as midpoints of the sites, then replace each site by the midpoint of its own interval.

Start at $0.1$, $0.2$, $0.9$. The boundaries are $0.15$ and $0.55$, so the cells are $[0, 0.15]$, $[0.15, 0.55]$ and $[0.55, 1]$, and the new sites are $0.075$, $0.35$ and $0.775$. One round has taken a badly clumped arrangement to a nearly even one. Another round gives boundaries $0.2125$ and $0.5625$ and sites $0.106$, $0.388$, $0.781$, and the arrangement is creeping towards the fixed point at $\frac{1}{6}$, $\frac{1}{2}$, $\frac{5}{6}$ — the arrangement that cuts the interval into three equal pieces, each site at its own centre.

Two things are visible in that computation and both survive into two dimensions. The first is the speed: the first round removes most of the error and the rest is a geometric crawl, with the ratio between successive errors settling to a constant. The second is that the fixed point is the even arrangement, which in one dimension is unique up to relabelling — and that uniqueness is exactly what fails in the plane.

In one dimension the cost has one minimum; in two it has many. Nothing about the rule changes. What changes is the shape of the landscape it descends, and that is a fact about the plane rather than about the algorithm.

Where it fails #

The descent argument is airtight and proves less than it appears to. It shows that the cost cannot rise, which means the iteration converges to a critical point of the cost — and critical points are not minima.

The cost function is genuinely not convex in the sites, and it has many local minima. Two nearby sites can settle into sharing a region that a single site would serve at a similar cost, and no local move separates them. Different starting scatters converge to different arrangements, all centroidal, with different costs. Nothing in the iteration escapes a local minimum, and nothing detects that it is in one.

That is the same shape as every other descent method on this site. Newton’s method converges to a root, and which root depends on where it starts in a way no local reasoning predicts. The honest summary is that the rule finds a good arrangement quickly and an optimal one never.

The failure has a concrete shape worth being able to recognise. Suppose the region is a long thin rectangle and there are four sites. The good arrangement puts them in a row. An arrangement with two of them side by side across the width, and the other two in a row beyond, is also centroidal — each of the four sits at the centre of its own cell, the arrows are zero, the iteration has stopped — and it costs more, because the two side-by-side cells are wider than they are tall and a cell’s cost grows with its longest dimension. No single site can improve matters by moving, because moving one alone makes its own cell worse. Escaping needs two sites to move at once, and the rule moves them one at a time by construction.

That is what a local minimum is in this setting: not a subtle numerical phenomenon, but a shape the arrangement has got itself into, which is locally stable and globally wrong.

There is a second failure and it is the more common one in practice. The square frame is not a neutral container: a hexagonal packing cannot meet a straight boundary cleanly, so the cells within two or three rows of the edge are pentagons and squares no matter how long the iteration runs. That is why the count above excludes the cells that touch the frame — a cell cut by the boundary has the number of sides the boundary gave it, and including them would report a fact about the frame as a fact about the rule.

The same rule under four other names #

The rule is worth following into its applications because in each of them it was invented independently.

Quantisation, which is where Lloyd found it. A signal takes continuous values and has to be stored as one of $k$ levels; the cost of the approximation is the mean squared error, which is exactly $E$ in one dimension. The optimal levels are the centroids of the intervals they represent, and the optimal interval boundaries are the midpoints between levels — a centroidal diagram on a line.

$k$-means clustering, where the sites are cluster centres and the square is replaced by a finite set of data points. Lloyd’s rule is the standard algorithm, its two half-steps are called assign and update, and the descent argument above is the standard proof that it terminates. It also inherits the local minima, which is why the algorithm is run several times from different starts.

Stippling, where an image is redrawn as a field of dots with the local density matching the image’s darkness. Weighting the integral by the image’s brightness makes the cells shrink where the image is dark, and the settled dots are the stipple. The visual quality of the result is a direct consequence of the hexagonal tendency: dots at random look blotchy, and dots on a settled diagram look even without looking gridded.

Mesh generation, where a region has to be filled with triangles that are not too thin. Sites settled by Lloyd’s rule have a Delaunay triangulation close to equilateral, which is precisely the property a numerical method wants of its mesh. This is the application that most often runs the iteration to convergence rather than stopping at two rounds.

Four subjects, one rule, and none of the four cites the others. That is normal for a rule this short: the cost of rediscovering it is smaller than the cost of finding out that somebody else has.

Lloyd’s own paper illustrates the point. It was written in 1957 as a Bell Labs internal memorandum about quantising a signal for transmission, and it was published in the open literature in 1982, by which time the same rule had appeared as the $k$-means algorithm — under that name from 1967 — and in half a dozen engineering settings. The 1982 paper is now the standard citation for something almost nobody learned from it.

The rediscoveries are not carelessness. Written in the vocabulary of quantisation the rule is “put each level at the mean of the values it represents”; written in the vocabulary of clustering it is “move each centre to the mean of its cluster”; written geometrically it is “move each site to its cell’s centroid”. Those three sentences do not share a word, and there is no index under which a person looking for one would find another. The same object arriving in two vocabularies is the recurring theme of this ladder, and here it produced four separate literatures.

What the pictures cannot show #

The cost is an integral and the pictures are polygons. Every number printed under a panel comes from an exact calculation over the drawn cells — the moment of a polygon about a point has a closed form, and the figure uses it rather than sampling — but that calculation happens inside the generator and no part of it is visible. A reader can see the cells even out; the claim that a particular quantity fell by a particular amount rests on the assertion, not on the drawing.

Convergence is asymptotic and the last panel is a round number. The arrows in the settled view are not zero, they are short, and no finite picture can distinguish a diagram at its fixed point from one very near it. The figure prints the longest remaining arrow for exactly that reason.

And the local minima are invisible. Every picture here shows a run that went somewhere sensible. A run that settled into a poor arrangement would look almost identical — even cells, mostly hexagons, arrows at zero — and would cost more. Drawing the failure would require drawing the better arrangement beside it, and finding that better arrangement is the problem the iteration cannot solve.

Where the ladder goes next #

This rung lets the sites move and keeps them equal. The next one keeps them still and makes them unequal: give each site a weight and watch the boundaries slide, which produces the diagram that describes a foam and turns out to be the version the lifting map really wanted. After that, the triangulation is read as a graph, and a spanning tree is found already inside it.

Sideways, the descent-to-a-critical-point argument is the same argument as gradient descent and has the same gap between what it proves and what is wanted, and the settled hexagons are the two-dimensional case of the sphere-packing question asked about cells rather than about spheres.

What is worth carrying away #

A rule that provably improves something is not the same as a rule that finds the best arrangement, and the difference is where all the interesting behaviour lives.

Lloyd’s iteration has a two-line proof of convergence, a first two rounds that do almost all of the work, a fixed point nobody designed that turns out to be a honeycomb, and a set of local minima it cannot see. Every one of those is a separate fact, and only the first is what the descent argument establishes.

The habit worth taking is to ask what a monotone quantity is bounded below by, and then to ask what its critical points look like. The first question is usually easy and gives convergence. The second is the one that says whether convergence is worth anything, and here it is the question whose answer is a hexagon.