[EM] Algorithmic redistricting: practical or not?

[Kristofer Munsterhjelm]

On 2025-02-11 05:13, Andrew B Jennings (elections) wrote:
> Kristofer,
> 
> Are you using actual geography and population data from somewhere?

Yes, I'm using the TIGER2023 census block dataset. My code, though quite 
messy, can be found here: https://github.com/kristomu/redistricter

> I got involved in a similar project recently. I met with some people
> who want to propose computer-based redistricting for the city I live in
> (Mesa, Arizona, USA).

> I do like the argument of Brian Olson (bdistricting.com) to use a
> metric-based approach instead of mandating a specific algorithm. I'm
> not especially concerned whether an optimal solution is tractable. I
> convinced those I was talking with that we could define a
> mathematical compactness metric and then allow public submissions (in
> addition to a city-proposed one) and choose the best one according to
> the metric. Now I just need to make sure I suggest to them the right
> metric! Ideally, we would set the quality metric that had zero
> tunable parameters, but I'm not sure that's possible.

An approach that mirrors mine would be:

- Each participant submits the latitude and longitude of each proposed 
district's center, as well as a weight associated with that district.
- The method assigns each census block with center coordinates c to 
district x that minimizes (d(c, x_c)^2 + w_x), where x_c is the provided 
district center and w_x is the provided weight for that district.
- If the difference between the most populated and least populated 
district exceeds a threshold, then the solution is disqualified.
- Otherwise, the solution is recorded with an objective (penalty) value 
equal to the sum of populated-weighted squared distances between the 
census block centers and the assigned district centers.
- Lowest objective value wins.

Unless you're dealing with very strange (en-/exclave style) census 
blocks, or inlets count,[1] this will give you contiguity due to the 
form of the objective function: if you move a census block assigned to a 
district closer to that district, it will still be assigned to that 
district. But the block boundaries will most likely spoil any chance of 
convexity.[2]

The only tunable here is the equipopulation threshold (and possibly what 
metric d you use: Euclidean, Earth surface distance, geoid, or something 
that takes land features into account).

You could then supply a local optimizer which, when given center 
coordinates, will find a reasonably good set of weights that enforces 
the population constraint.[3] Janáček and Gábrišová provide an algorithm 
for this. I implemented a similar one in kmeans.py, though my version 
would sometimes fail to converge, which I'm assuming is due to a bug. 
J&G also suggest a postprocessing step to get very close to population 
equality without losing compactness, but I haven't implemented that one.

It might also be possible to get a good initial set of weights by 
looking at the duals of the MIP, but I haven't investigated that in detail.

> Our city touches 6000 populated census blocks or 333 populated block
> groups or 137 populated census tracts.
> 
> I feel like the census tracts are too coarse to allow a good 
> solution. The blocks seem too small, so my first attempt has been 
> with the block groups. (Though I note that the current districts 
> split block groups in a couple places.) I don't feel like the 
> districts have to have  equal population to within 1 person. The 
> group I was talking with  thought 3% was reasonable. That was before
> they knew that the current districts have a population discrepancy of
> 8% between the largest and the smallest. (That is from my analysis
> based on the 2020 census and the maps. There may be some more
> official calculation where the discrepancy is smaller.)

I used Colorado as the standard problem when writing my code. It has 
around 200k census blocks. In practice I would sample around 1000 of 
them using population-weighted k-means++ for the initial step of finding 
reasonable centers. For the weight determination, I would use all 200k 
blocks since optimizing to pass the population constraint has to be done 
fully disaggregated.

This worked, though the second part was pretty slow. I don't think 6000 
blocks would be a problem.

> I wasn't planning on mandating convexity, though I hope the right
> compactness function would favor it as much as possible. I was also
> hoping that contiguity would come for free. It's probably not acceptable
> to have a non-contiguous map. Maybe we will need to say that we'll
> accept that best map (according to the compactness metric) that also
> preserves contiguity.

The k-means approach that I use (and that Brian Olson also seems to use) 
would give you convexity if you were dealing with point sets. But even 
small census blocks can have borders that don't fit. E.g. suppose a 
sloping district made out of square census blocks; you'll get a 
stairstep boundary which isn't convex.

I looked a bit into assuming a uniform distribution of people within 
each block and then making boundaries that cut across them. This *is* 
convex, but dealing with point in polygon queries is a real hassle. I 
don't know if legislation allows it either.

Some example images can be found at 
https://github.com/kristomu/redistricter/tree/master/examples.

> Olson is promoting as a metric "lowest average distance over all
> people to the geographic center of their district". I think he's using
> census blocks. Not sure if he's mandating contiguity, or what level of
> population discrepancy is tolerated.
> 
> Like you, I am leaning toward using the square of distance. I know I
> said I don't care about tractability, but it's just so much easier to
> compute without the square roots. It is also what is minimized when
> finding any geometric centroid, so it seems defensible.
> 
> Also the sum of (population weighted) squares of distances to the
> centroid is basically equivalent to the sum of squares of distances
> between every pair, which feels elegant. And it can be converted into
> "for every pair of block groups, there is a cost whenever they are in
> the same district" max-cut graph problem. Alternatively, it makes naive
> brute-force searching much faster.
> 
> Taking that metric (sum of squares of distances) and finding local
> minima gives good looking districts, but they vary in population by
> 200% or more. I doubt that brute force or MILP optimizer can get it
> down to 3-8%, but I'm going to keep trying things.
One thing I haven't tried that you could is to use a hard-capacitated 
(population-constrained, but not necessarily contiguous) solution as a 
Lloyd's algorithm step:
	1. Start with some initial centers.
	2. Calculate the hard-cap solution at a convenient granularity level 
(e.g. 1000 population-weighted k-means++ samples).
	3. Calculate the centroid points for each district. Mixed points 
(divided between more than one district) contribute proportionally to 
their centroids.
	4. Set the district centers to these centroids and go to 2.

Alternatively (at the cost of having to do a coarser sampling) you could 
use J&G's swap compactness constraint (theorem 4) to get contiguity in 
step 2.

If you're using or decide to use k-means++, this paper might also be 
useful: https://proceedings.mlr.press/v119/choo20a/choo20a.pdf

> 
> I have tried a couple other things with the cost function:
> 
> 1. Instead of using the distance between the census-given "interior
> point" of each block, use the distance between the closest points for
> each pair of blocks.
> 
> 2. A penalty for uneven district sizes (add to the cost function a
> large constant times the population of both block groups in the pair)
> 
> I can tune the constant to get population discrepancies down to 3-8%
> and districts that are almost contiguous, although I hate having a
> constant that can be tuned by politicians. Maybe I can find some way to
> formulate it based on city population, city geography, and the number of
> districts, or something.

I agree, having tunable parameters is a problem. If the method were to 
be used by a neutral commission, then that commission could decide what 
values would provide the best trade-off, but not so much if it's partisan.

One thought I've had as far as uneven district sizes is to limit maximum 
distance to the center relative to minimum distance. This can 
unfortunately only be done as a constraint, not as a penalty. For MIP, 
something like:

for all i, x,y:
	min_dist[i] <= assign_binary[i, x, y] * d(centeri_x, centeri_y, x, y) + 
(1 - assign_binary[i, x, y]) * M
	max_dist[i] >= assign_binary[i, x, y] * d(centeri_x, centeri_y, x, y)

for all i:
	max_dist[i] <= K * min_dist[i]

for a large value M and for some constraint parameter K, and 
assign_binary[i, x, y] being 1 if assign[i, x, y] is greater than some 
threshold value (say 1e-4), zero otherwise.

While this might be interesting, in practice, plain old center of 
gravity optimization doesn't seem too bad, see eg. ch 8 and 9 of 
https://rangevoting.org/SpannGulottaKaneMCM2007.pdf.

> Or maybe I just need to try using the 6000 census blocks and see what happens.
> 
> ~ Andy Jennings

[1] By "inlets" I'm thinking of a situation where some feature (inlet, 
estuary, etc), or concavity of the state boundary itself, divides a 
region of land, and there's one block on each side. Then it's 
theoretically possible that the two blocks are assigned to the same 
district, but every inland connection between the two passes through 
some other district.

[2] A better option might be to assign a block to a district A instead 
of B if more points in that block are closer to A than to B. That should 
fix exclave district problems, if I'm not mistaken, but area 
calculations for arbitrary (nonconvex) polygons can again be a pain.

[3] I have a hunch that you could get rid of the weights by calculating 
the centroids (or Frechet mean) of the assigned districts to get a new 
set of centers, and that the simple k-means allocation based on these 
centers would pass the population constraint. But I haven't tested or 
proved this.