To Kill A Gerrymander

The Method

How Convex Optimal
Districting Works

A plain-language explanation of the algorithm — and the mathematics that proves it works.

The Problem with Gerrymandering

Gerrymandering doesn't just produce ugly maps. It reduces the number of genuinely contested elections, pushing representatives toward the extremes and away from the voters in the middle. The diagram below shows how different redistricting strategies affect the number of competitive districts — and therefore the health of democracy itself.

Diagram showing how gerrymandering reduces contested districts and increases political polarization

How different redistricting strategies affect contested elections. Gerrymandering — whether partisan or bipartisan — reduces competition and accelerates polarization.

The Convex Optimal Solution

The key insight is simple: if the algorithm doesn't know how people vote, it cannot gerrymander. Convex optimal districting takes only population geography as its input — where people live, nothing else. It then finds the unique partition of the state into districts that minimizes the average distance from each person to their district's center, using a notion of "average" that guarantees every district is convex.

Because the solution is mathematically unique, there is no room for manipulation. The map draws itself.

The Algorithm

The method works iteratively, refining districts until they are both equal in population and geometrically stable:

Start with random initial district centers placed across the region.

Assign each population block to the center minimizing f(distance to center) − bias, where f(x) = x² for flat geometry and f(x) = ln(sec²x) for spherical geometry. The bias controls population balance.

Measure population imbalance. Increase the bias for underpopulated districts (making them more attractive, so they grow) and decrease it for overpopulated ones (making them less attractive, so they shrink). Keep doing this until equal populations are achieved in each district.

Move each center to the population-weighted point that minimizes the average f(distance) to all assigned population blocks.

Repeat steps 2–4 until districts are equal in population and boundaries stop changing.

The result is a set of districts that are compact, convex, and equal in population — determined entirely by where people live.

The Earth, Partitioned

To demonstrate the algorithm at global scale, Andrew Bray applied it to the entire world's population, dividing Earth into equal-population regions. The maps below show the result for 5 to 60 regions. Each colored region contains the same number of people. Although the boundaries appear curved on these maps, that is purely an illusion of the Mercator projection. On a globe, every boundary is a portion of a great circle — like the equator or a meridian — the straightest possible line on a curved surface. Every region is mathematically convex.
World divided into 5 equal-population convex regions
5 regions. Five convex regions of equal population. All boundaries are great circles.
World divided into 7 equal-population convex regions
7 regions. Two additional regions appear as population centers separate. All regions remain convex.
World divided into 10 equal-population convex regions
10 regions. The algorithm continues to produce convex regions with great-circle boundaries, regardless of how complex the population distribution becomes.
World divided into 13 equal-population convex regions
13 regions. Dense population zones in Asia require smaller regions; sparse zones remain large. All convex throughout.
World divided into 16 equal-population convex regions
16 regions. The bias terms automatically adjust region sizes to maintain equal population while preserving convexity.
World divided into 20 equal-population convex regions
20 regions. As the number of regions grows, boundaries multiply — but every region remains a convex polygon on the sphere.
World divided into 40 equal-population convex regions
40 regions. A fine-grained partition. The apparent curvature of boundaries is entirely due to the Mercator projection — on a globe, all boundaries are straight.
World divided into 60 equal-population convex regions
60 regions. The most detailed partition shown. Every region is convex, every region contains the same population, and every boundary is a great circle.

For the Mathematically Curious

The following section presents the mathematical foundations of convex optimal districting — why squared distance produces straight boundaries on a flat map, and why a modified distance function produces great-circle boundaries on a sphere.

Why the Boundaries Are Straight

The diagram below shows the key insight: the same algebraic structure that makes d² produce straight-line boundaries on a flat surface has a spherical counterpart — using ln(sec²d) — that produces geodesic (great-circle) boundaries on a sphere. This means the algorithm works correctly on the curved Earth, not just on flat maps.

Plane vs spherical geometry comparison showing the Pythagorean structure behind straight boundaries

Left: flat geometry where f(x) = x². Right: spherical geometry where f(x) = ln(sec²x). In both cases the same algebraic cancellation produces straight boundaries.

The Energy Function and Equal Population

Shows why the bias terms enforce equal population — and then vanish entirely from the solution.

Key insight: The bias terms are the mechanism that enforces equal population during optimization. But at the solution, they cancel out exactly — leaving a pure minimization of total f(distance) subject to equal population.
The Energy Function

The algorithm minimizes (over centers and assignments) and maximizes (over biases) the following energy function. For concreteness, consider 8 regions and 24 population points:

E =   8∑k=1 3bk  +  24∑i=1 ( f(distance(Xi, Cκ(i))) − bκ(i) )

where Xi is the i-th population point, Ck is the center of region k, bk is the bias for region k, and κ(i) is the assignment of point i to a region. The function f is f(x) = x² for flat geometry and f(x) = ln(sec²x) for spherical geometry.

Equal Population from Maximizing Over Biases

Taking the partial derivative of E with respect to bk and setting it to zero:

∂E/∂bk = 3 − (number of points assigned to region k) = 0

This forces each region to contain exactly 3 points — equal population. The equal population constraint is not imposed externally; it emerges naturally from maximizing E over the biases.

The cancellation: At the min-max solution, each region contains exactly 3 points. Substituting back into E: the first summation contributes 3·(b₁ + b₂ + ⋯ + b₈), and the bias terms in the second summation contribute −3·(b₁ + b₂ + ⋯ + b₈). These cancel exactly, regardless of the values of the biases, leaving:
E = 24∑i=1 f(distance(Xi, Cκ(i)))
Conclusion: At the solution, the biases vanish. We are genuinely minimizing the total f(distance) — and therefore the average f(distance) — subject to equal population. Minimizing the average f(distance) is equivalent to minimizing f⁻¹ of the average f(distance), which is what we define as the average distance with respect to f.

What is an average distance with respect to f? Consider a simple example: what is the average of 1 and 7? The usual answer is (1 + 7)/2 = 4. But with respect to f(x) = x², it is \(\sqrt{\frac{1^2 + 7^2}{2}} = \sqrt{\frac{50}{2}} = \sqrt{25} = 5\). This f-average gives extra weight to larger values because squaring amplifies them.

Why use this unusual average? As shown in the proofs above and below, only by minimizing the average distance with respect to f(x) = x² on a flat surface — or f(x) = ln(sec²x) on a sphere — do the boundary equations reduce to straight lines. Convexity is not an added constraint. It is a mathematical consequence of the choice of f.
Local vs. Global Minimum

The min-max argument guarantees a local minimum, not necessarily a global one. Different starting configurations may converge to different solutions. To find the best map, the algorithm is run many times from random initial district centers, and the solution with the smallest value of E — the smallest total f(distance) — is selected. This is the convex optimal solution.

Proof: Straight-Line Boundaries on a Flat Surface

Minimizing squared Euclidean distance produces linear boundaries between any two regions.

Claim: If each point is assigned to the center minimizing (distance² − bias), the boundary between any two regions is a straight line.
Setup

Let C₁ = (x₁, y₁) and C₂ = (x₂, y₂) be two centers with biases b₁ and b₂. The boundary is the set of points where both centers are equally attractive:

distance₁² − b₁ = distance₂² − b₂
Step 1 — Pythagorean Theorem

Expand using the standard Euclidean formula:

distance₁² = (x − x₁)² + (y − y₁)² distance₂² = (x − x₂)² + (y − y₂)²
Step 2 — Expand and subtract
(x² − 2xx₁ + x₁² + y² − 2yy₁ + y₁²) − (x² − 2xx₂ + x₂² + y² − 2yy₂ + y₂²) = b₁ − b₂
Critical observation: The x² and y² terms appear identically on both sides and cancel completely — even with bias terms present. This cancellation is unique to squared Euclidean distance.
x² − x² = 0     y² − y² = 0
Step 3 — What remains is a line

Only linear terms survive, matching the standard form Ax + By + C = 0:

x(2x₂ − 2x₁) + y(2y₂ − 2y₁) + [(x₁² + y₁²) − (x₂² + y₂²) − (b₁ − b₂)] = 0

The biases shift the line's position but introduce no curvature.

Conclusion: The boundary is a straight line. Using unsquared Euclidean distance would leave cross-terms that produce a hyperbola instead.

Proof: Geodesic Boundaries on a Curved Surface

Replacing d² with ln(sec²d) recovers straight boundaries — great circles — on a sphere.

Claim: On a sphere, if each point is assigned to the center minimizing ln(sec²d) − bias, the boundaries are great circles (geodesics).
Setup

Let C₁ and C₂ be two centers on the equator with biases b₁ and b₂. The boundary satisfies:

F(d₁) − b₁ = F(d₂) − b₂    where F(d) = ln(sec²d)
Step 1 — Spherical Pythagorean Theorem

On a sphere, d² = p² + h² becomes the Spherical Pythagorean Theorem for a right triangle:

cos(d) = cos(p) · cos(h)
Step 2 — Transform to additive form

Apply the secant identity and take the natural log. Multiplication becomes addition:

sec²(d) = sec²(p) · sec²(h) ln(sec²d) = ln(sec²p) + ln(sec²h)

This is the Curved Pythagorean Theorem in additive form — structurally identical to the flat case.

Step 3 — Substitute and apply biases

Where h is the shared off-axis height above the equator:

[ln(sec²p₁) + ln(sec²h)] − b₁ = [ln(sec²p₂) + ln(sec²h)] − b₂
Critical observation: ln(sec²h) appears identically on both sides and cancels completely — just as x² cancelled in the flat proof.
ln(sec²h) − ln(sec²h) = 0
Step 4 — What remains is a great circle

In this setup, with both centers on the equator, the boundary depends only on p (position along the equator), not on h (height above it). Since the solution is independent of h, the boundary does not curve as you move away from the equator — it forms a meridian, perpendicular to the equator and to the line connecting the two centers. A meridian is a great circle.

In the general case, where the two centers are anywhere on the sphere, the same algebraic cancellation occurs along whatever great circle connects them. The boundary is always a great circle arc — the straightest possible line on a curved surface — perpendicular to the great circle connecting the two centers.

Conclusion: Minimizing F(d) = ln(sec²d) produces great-circle boundaries on the curved surface of the Earth — regardless of where the district centers are located. This is the spherical counterpart of straight-line boundaries on a flat map, and it is why the regions in the world maps are convex on the globe, even though they appear to have curved boundaries under the Mercator projection.

References

Cohen-Addad, V., Klein, P. N., & Young, N. E. (2018). Balanced power diagrams for redistricting. arXiv preprint arXiv:1710.03358v3.

NASA GPW v4 (2020) Population Density Data. Global grid of population estimates.

U.S. Census Bureau (2020). Redistricting Data (PL 94-171). Official block-level population data.

Convex optimal maps and proofs by Andrew Bray. World partition maps computed in MATLAB.