Algorithmic redistricting: practical or not?
So I got myself entangled in redistricting back in November, after
writing a post here about how solving a related problem (that I thought
was the same as districting) could be done through linear programming,
if the centers were given.
Some months and lots of code later, I thought I'd post about what I
found out.
The summary, first:
- Bad news: the linear program I gave doesn't ensure contiguity or
convexity, and my attempts to forcing it were too expensive.
- Good news: the linear program isn't useless and can be used as a
component in a branch-and-bound-like approach. Once it has found a good
choice of centers, a local search can be used to find a convex solution.
- All of this is much harder than just using proportional
representation, so it's mostly an intellectual exercise.
To recap, the linear program is:
For centers i = 1...n, with cneterx_i, centery_i being the (x,y)
coordinates of center i
minimize the sum over centers i, and points x, y inside the state:
(1) assign[i, x, y] * (pop(x, y) + eps) * dist^2(centerx_i,
centery_i, x, y)
subject to
(2) for all i, x, y: 0 <= assign[i, x, y] <= 1
(3) for all i: (sum over x, y: assign[i, x, y] * pop(x, y)) - tpop/n = 0
(4) for all x, y: (sum over i: assign[i, x, y]) = 1
(1) defines the problem: minimize expected squared distance from a
random inhabitant to their assigned center. Assign[i, x, y] is a free
variable that is one if the point (x,y) is assigned to center i, zero
otherwise, and pop(x, y) is the number of people at that reside closer
to (x,y) than any other point. The epsilon value eps ensures that a
point with no population is assigned to its closest center.[1]
(2) defines the range of the assign values.
(3) implements the population constraint that every district has equal
population.
(4) requires that every point be assigned to some district.
The linear program does solve the problem as specified. However, the
resulting districts aren't contiguous.
You'd think they were because suppose district A picks a point that's
inside district B, then shouldn't we be able to improve the solution for
both by swapping this assignment with points closer to A? The problem is
the population constraint. Sometimes it's easier to skip points because
otherwise, a bunch of other points would have to change to keep equal
populations, and the sum of these changes could be worse.
So instead you need additional constraints, and these are expensive.[2]
That's the bad news.
The good news is that it's actually not as bad as you'd think, because
the linear program I gave provides a good lower bound for a contiguous
solution. Furthermore, with integer programming, it's possible to add a
large number of candidate points and have the solver choose only some of
these as centers. For instance, we may pick 150 candidate centers with
something like kmeans++; and then use an integer program to say it
should choose 8 of them, like this:
minimize the sum over centers i, and points x, y inside the state:
assign[i, x, y] * (pop(x, y) + eps) * dist^2(centerx_i, centery_i, x, y)
subject to
(2b) for all i, x, y: 0 <= assign[i, x, y] <= chosen[i]
for all i: (sum over x, y: assign[i, x, y] * pop(x, y)) - tpop/k *
chosen[i] = 0
for all x, y: (sum over i: assign[i, x, y]) = 1
(5) sum over i: chosen[i] = k
free variables: assign[i, x, y] continuous, chosen[i] binary
(2b) says: if chosen[i] is zero, then nobody can be assigned to center
i. Otherwise, points may be assigned to center i.
(5) forces exactly k districts to be chosen. The population constraint
(3) has also been updated accordingly: if chosen[i] is zero, then
population must be zero.
(Thanks to Joshua Boehme for this idea.)
When run, the program will find the best k centers and give, as its
objective value, a lower bound on the sum of squared distances that
we're actually interested in. If this is worse than the current record,
we just discard the solution as it can't be better anyway; otherwise, we
check it more closely with contiguity constraints, which will return the
true objective value.
This seems to work reasonably well in theory, but the datrasets are much
too large to run contiguity-constrained programs on. Instead I use a
method from Janáček and Gábrišová[2], which uses the property that,
suppose we start with some given centers,
(centerx_i, centery_i) for i = 1..k
Then there exist weights w_i so that the "ordinary" k-means clustering
(Voronoi clustering) objective
minimize sum over i, x, y: assign[i, x, y] * (pop(x, y) + eps) *
(dist^2(centerx_i, centery_i, x, y) + w_i)
will ensure that the districts have equal populations. (This makes
intuitive sense: if w_i is very large, district i is penalized and
shrinks relative to the other districts.) And this "uncapacitated" or
unconstrained problem can be easily solved by checking each point's
value for each district and picking the best.
So I ended up with a two-step procedure: use integer programming to find
reasonably good centers from a randomly selected large enough candidate
pool, and then use a local search to update the weights w_i until the
districts had nearly equal populations, thus enforcing contiguity and
convexity.[3][4]
In practice, there were a bunch of additional problems due to census
district shapes, non-Euclidean geometry, and how to handle large amounts
of data quickly. But I won't go into those here; perhaps later, as this
post is long enough.
-km
[1] I made a mistake back when I stated this back at the end of October:
the objective value should be weighted by population, not just sum of
squared distances. Otherwise every district wants to be exactly the same
size as the others, which leads them all to carve inwards to cut shares
of cities. See the unweighted objective picture here:
https://github.com/kristomu/redistricter/tree/master/undesired_outcomes
[2] See JANÁČEK, Jaroslav; GÁBRIŠOVÁ, Lýdia. A two‐phase method for the
capacitated facility problem of compact customer sub‐sets. Transport,
2009, 24.4: 274-282.
[3] As far as I understand, it's possible to use a postprocessing phase
to force exact equality. The Janáček and Gábrišová paper gives more detail.
In my tests for Colorado, the local search converged to around 0.02%
difference in population between the most populated district and the
least, so I didn't implement it.
[4] There may be better ways to find the weights than local search, see
e.g. SPANN, Andrew; KANE, Daniel; GULOTTA, Dan. Electoral redistricting
with moment of inertia and diminishing halves models. The UMAP Journal,
2007, https://rangevoting.org/SpannGulottaKaneMCM2007.pdf, particularly
pages 9 and 10.
Kristofer,
Are you using actual geography and population data from somewhere?
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.
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 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.
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.
I have tried a couple other things with the cost function:
-
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.
-
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.
Or maybe I just need to try using the 6000 census blocks and see what happens.
~ Andy Jennings
On Tuesday, February 4th, 2025 at 2:39 PM, Kristofer Munsterhjelm km-elmet@munsterhjelm.no wrote:
So I got myself entangled in redistricting back in November, after
writing a post here about how solving a related problem (that I thought
was the same as districting) could be done through linear programming,
if the centers were given.
Some months and lots of code later, I thought I'd post about what I
found out.
The summary, first:
- Bad news: the linear program I gave doesn't ensure contiguity or
convexity, and my attempts to forcing it were too expensive. - Good news: the linear program isn't useless and can be used as a
component in a branch-and-bound-like approach. Once it has found a good
choice of centers, a local search can be used to find a convex solution. - All of this is much harder than just using proportional
representation, so it's mostly an intellectual exercise.
To recap, the linear program is:
For centers i = 1...n, with cneterx_i, centery_i being the (x,y)
coordinates of center i
minimize the sum over centers i, and points x, y inside the state:
(1) assign[i, x, y] * (pop(x, y) + eps) * dist^2(centerx_i,
centery_i, x, y)
subject to
(2) for all i, x, y: 0 <= assign[i, x, y] <= 1
(3) for all i: (sum over x, y: assign[i, x, y] * pop(x, y)) - tpop/n = 0
(4) for all x, y: (sum over i: assign[i, x, y]) = 1
(1) defines the problem: minimize expected squared distance from a
random inhabitant to their assigned center. Assign[i, x, y] is a free
variable that is one if the point (x,y) is assigned to center i, zero
otherwise, and pop(x, y) is the number of people at that reside closer
to (x,y) than any other point. The epsilon value eps ensures that a
point with no population is assigned to its closest center.[1]
(2) defines the range of the assign values.
(3) implements the population constraint that every district has equal
population.
(4) requires that every point be assigned to some district.
The linear program does solve the problem as specified. However, the
resulting districts aren't contiguous.
You'd think they were because suppose district A picks a point that's
inside district B, then shouldn't we be able to improve the solution for
both by swapping this assignment with points closer to A? The problem is
the population constraint. Sometimes it's easier to skip points because
otherwise, a bunch of other points would have to change to keep equal
populations, and the sum of these changes could be worse.
So instead you need additional constraints, and these are expensive.[2]
That's the bad news.
The good news is that it's actually not as bad as you'd think, because
the linear program I gave provides a good lower bound for a contiguous
solution. Furthermore, with integer programming, it's possible to add a
large number of candidate points and have the solver choose only some of
these as centers. For instance, we may pick 150 candidate centers with
something like kmeans++; and then use an integer program to say it
should choose 8 of them, like this:
minimize the sum over centers i, and points x, y inside the state:
assign[i, x, y] * (pop(x, y) + eps) * dist^2(centerx_i, centery_i, x, y)
subject to
(2b) for all i, x, y: 0 <= assign[i, x, y] <= chosen[i]
for all i: (sum over x, y: assign[i, x, y] * pop(x, y)) - tpop/k *
chosen[i] = 0
for all x, y: (sum over i: assign[i, x, y]) = 1
(5) sum over i: chosen[i] = k
free variables: assign[i, x, y] continuous, chosen[i] binary
(2b) says: if chosen[i] is zero, then nobody can be assigned to center
i. Otherwise, points may be assigned to center i.
(5) forces exactly k districts to be chosen. The population constraint
(3) has also been updated accordingly: if chosen[i] is zero, then
population must be zero.
(Thanks to Joshua Boehme for this idea.)
When run, the program will find the best k centers and give, as its
objective value, a lower bound on the sum of squared distances that
we're actually interested in. If this is worse than the current record,
we just discard the solution as it can't be better anyway; otherwise, we
check it more closely with contiguity constraints, which will return the
true objective value.
This seems to work reasonably well in theory, but the datrasets are much
too large to run contiguity-constrained programs on. Instead I use a
method from Janáček and Gábrišová[2], which uses the property that,
suppose we start with some given centers,
(centerx_i, centery_i) for i = 1..k
Then there exist weights w_i so that the "ordinary" k-means clustering
(Voronoi clustering) objective
minimize sum over i, x, y: assign[i, x, y] * (pop(x, y) + eps) *
(dist^2(centerx_i, centery_i, x, y) + w_i)
will ensure that the districts have equal populations. (This makes
intuitive sense: if w_i is very large, district i is penalized and
shrinks relative to the other districts.) And this "uncapacitated" or
unconstrained problem can be easily solved by checking each point's
value for each district and picking the best.
So I ended up with a two-step procedure: use integer programming to find
reasonably good centers from a randomly selected large enough candidate
pool, and then use a local search to update the weights w_i until the
districts had nearly equal populations, thus enforcing contiguity and
convexity.[3][4]
In practice, there were a bunch of additional problems due to census
district shapes, non-Euclidean geometry, and how to handle large amounts
of data quickly. But I won't go into those here; perhaps later, as this
post is long enough.
-km
[1] I made a mistake back when I stated this back at the end of October:
the objective value should be weighted by population, not just sum of
squared distances. Otherwise every district wants to be exactly the same
size as the others, which leads them all to carve inwards to cut shares
of cities. See the unweighted objective picture here:
https://github.com/kristomu/redistricter/tree/master/undesired_outcomes
[2] See JANÁČEK, Jaroslav; GÁBRIŠOVÁ, Lýdia. A two‐phase method for the
capacitated facility problem of compact customer sub‐sets. Transport,
2009, 24.4: 274-282.
[3] As far as I understand, it's possible to use a postprocessing phase
to force exact equality. The Janáček and Gábrišová paper gives more detail.
In my tests for Colorado, the local search converged to around 0.02%
difference in population between the most populated district and the
least, so I didn't implement it.
[4] There may be better ways to find the weights than local search, see
e.g. SPANN, Andrew; KANE, Daniel; GULOTTA, Dan. Electoral redistricting
with moment of inertia and diminishing halves models. The UMAP Journal,
2007, https://rangevoting.org/SpannGulottaKaneMCM2007.pdf, particularly
pages 9 and 10.
Election-Methods mailing list - see https://electorama.com/em for list info
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:
-
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. -
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.