Voronoi tessellations and Lloyd's algorithm

Place $k$ points $p_1, \ldots, p_k$ in the plane (called generators). The Voronoi cell of $p_i$ is

V_i \;=\; \{x \in \mathbb{R}^2 : \|x - p_i\| \le \|x - p_j\| \text{ for all } j \ne i\}.

Each cell is a convex polygon, and the cells together partition the plane. Random generators give irregular cells of varying sizes and shapes. Carefully placed generators give regular tilings. There is a procedure that takes any starting configuration and pulls it toward a maximally regular one.

Lloyd’s algorithm

At each step:

Compute the Voronoi cells $V_1, \ldots, V_k$ for the current generators.
Replace each generator $p_i$ with the centroid of its cell: $p_i' \;=\; \frac{1}{|V_i|} \int_{V_i} x \, dx.$
Repeat.

This is Lloyd’s algorithm (Lloyd, 1957/1982). It decreases the total squared distance from each point in the plane to its nearest generator at every step (a Lyapunov / energy function), so it converges to a local minimum.

The fixed points are configurations where every generator is already at the centroid of its own cell. These are called centroidal Voronoi tessellations (CVTs). On flat domains, the cells of a CVT are nearly hexagonal: hexagons are the polygon shape that minimizes the second moment of area per unit area, so the algorithm seeks them out (see The bees for more details on this minimization).

Click on the canvas to add generators. Drag a generator to move it. Press step to run one Lloyd iteration, or run to step until convergence. Watch the cells round out and the configuration regularize.

Press random 24 generators for a fresh starting configuration, then run 30 steps to watch Lloyd converge. Press step if you want to see the iteration in single-step mode.

After enough iterations on a rectangular domain, the cells become approximately hexagonal in the interior and adapt to the boundary near the edges. The convergence is monotone in the energy

F(p_1, \ldots, p_k) \;=\; \sum_{i=1}^{k} \int_{V_i} \|x - p_i\|^2 \, dx,

the total squared distance to the nearest generator. Each step strictly decreases $F$ unless the configuration is already a CVT.

Where this connects

k-means clustering. Lloyd’s algorithm is the algorithm typically called “k-means”. Given a finite point cloud and $k$ cluster centers, k-means alternates “assign each point to the nearest center” (the discrete analog of computing Voronoi cells) and “move each center to the centroid of its cluster” (the discrete analog of step 2). Continuous Lloyd’s is the limiting case where the point cloud is replaced by a uniform measure on a domain. k-means inherits its convergence properties from the same energy-decrease argument.
Vector quantization. Designing optimal codebooks for compressing continuous signals. The optimal $k$ -point codebook for squared error is exactly a CVT of the source distribution.
Mesh generation. CVT meshes are popular for finite-element methods because they have well-conditioned cells (no slivers, near-uniform sizes).
Stippling and natural patterns. Lloyd’s algorithm produces visually pleasing “blue noise” point distributions used in computer graphics; biological tissues and bee compound eyes settle into similar hexagonal arrangements through different mechanisms but for similar reasons.

The bees

Bees build approximately hexagonal honeycombs, and they are often described as running Lloyd’s algorithm in nature. Not quite: there are two distinct mathematical optimality statements that both have the hexagonal tiling as their answer in 2D, and bees and Lloyd’s algorithm solve different ones.

Hales’s Honeycomb Theorem (Hales, 2001). Among all tilings of the plane into regions of equal area, the regular hexagonal tiling minimizes the total perimeter. This is the optimization that minimizes wax usage for a given total storage volume. The conjecture goes back at least to Pappus of Alexandria in the fourth century AD; the proof is modern.

Lloyd’s CVT in the flat-domain, uniform-density limit minimizes a different functional, the second moment of area per unit area:

\frac{1}{|V|} \int_V \|x - p\|^2 \, dx,

where $p$ is the centroid of $V$ .

The integrand $\|x - p\|^2$ is the squared distance from a point of the cell to its centroid. Its integral $\int_V \|x - p\|^2 \, dx$ is the second moment of area of $V$ about its centroid: a measure of how spread out the cell is, small for compact round cells and large for elongated or branching ones. Dividing by $|V|$ gives the average value of $\|x - p\|^2$ over a uniformly random point of $V$ .

For shapes of fixed area, this average depends only on the shape of $V$ , not its size. The minimum over all shapes of a given area is achieved by the disk (round cells have the lowest spread), but disks do not tile the plane. The relevant question for Lloyd’s algorithm is therefore: among shapes that do tile the plane, which one is closest to a disk in this sense? With cells of unit area, the comparison among the three regular polygons that tile, plus the disk for reference, is

disk: $\dfrac{1}{2\pi} \approx 0.1592$ (the absolute minimum over all shapes),
regular hexagon: $\dfrac{5}{18\sqrt{3}} \approx 0.1604$ ,
square: $\dfrac{1}{6} \approx 0.1667$ ,
equilateral triangle: $\dfrac{\sqrt{3}}{9} \approx 0.1925$ .

The hexagon comes within $0.7\%$ of the disk lower bound; the square is about $5\%$ over, the triangle about $21\%$ over. Among polygons that tile the plane (not only the regular ones), the regular hexagon minimizes this functional (L. Fejes Tóth, 1959). The analogous statement in $d \ge 3$ , that the optimal tiling for moment of inertia is by truncated octahedra, is Gersho’s conjecture, still open in general dimensions.

So the perimeter minimum and the moment-of-inertia minimum are mathematically distinct objectives, but they happen to be solved by the same shape in 2D. That is the reason bees and k-means converge on hexagons through completely unrelated processes.

The bee mechanism is not direct geometry. Bees do not measure angles. The accepted picture, going back to D’Arcy Thompson (1917) and experimentally confirmed by Karihaloo, Zhang, and Wang (2013), is that bees first build approximately circular cells out of soft wax, and the wax then relaxes under the heat of the bees’ bodies and the equal mutual pressure of neighboring cells. Surface tension drives the boundaries to a configuration that minimizes wax for fixed cell area, which by the Honeycomb Theorem is the hexagonal tiling.

In this sense, bees do calculus the same way soap films do calculus: the physical system minimizes the same variational integral that the mathematical statement names. Lloyd’s algorithm solves a different but adjacent variational problem by direct iteration, and arrives at the same answer because the answer is hexagons either way.

References

Stuart P. Lloyd. Least squares quantization in PCM. Bell Laboratories Technical Note, 1957; published in IEEE Transactions on Information Theory, 28(2):129-137, 1982.
Qiang Du, Vance Faber, Max Gunzburger. Centroidal Voronoi tessellations: Applications and algorithms. SIAM Review, 41(4):637-676, 1999.
Atsuyuki Okabe, Barry Boots, Kokichi Sugihara, Sung Nok Chiu. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. Wiley, 2000.
Thomas C. Hales. The honeycomb conjecture. Discrete and Computational Geometry, 25(1):1-22, 2001.
L. Fejes Tóth. What the bees know and what they do not know. Bulletin of the American Mathematical Society, 70(4):468-481, 1964.
B. L. Karihaloo, K. Zhang, J. Wang. Honeybee combs: how the circular cells transform into rounded hexagons. Journal of the Royal Society Interface, 10(86):20130299, 2013.
D’Arcy Wentworth Thompson. On Growth and Form. Cambridge University Press, 1917.
Allen Gersho. Asymptotically optimal block quantization. IEEE Transactions on Information Theory, 25(4):373-380, 1979.