# 1953: Fermi’s “little discovery”
and the birth of the numerical experiment

Stefano Lepri  Roberto Livi  Stefano Ruffo  Consiglio Nazionale delle
Ricerche, Istituto dei Sistemi Complessi, Via Madonna del Piano 10 I-50019
Sesto Fiorentino, Italy  Dipartimento di Fisica e Astronomia, Università di
Firenze, Via G. Sansone 1, I-50019
Sesto Fiorentino, Italy  SISSA Via Bonomea 265, I-34136 Trieste, Italy

###### Abstract

The year 1953 is pivotal for computational physics: the first application of
the Monte-Carlo method is published and calculations of the so-called Fermi-
Pasta-Ulam-Tsingou experiment are started. It is the beginning of the massive
use in the physical sciences of numerical methods implemented on electronic
computers and a decisive step in the development of modern nonlinear dynamics.
This will lead to an unpredictable development during the following 70 years.
We briefly review the unfolding of these events and present some recent
results that show how the issues raised are still relevant today.

††journal: Giornale di Fisica

##  1 A step back in time

In 1953, the Los Alamos National Laboratory was a complex of buildings in the
New Mexico desert. The facility had been established ten years earlier as a
top secret research center for the development of nuclear weapons during World
War II, and was still heavily involved in that task. Let’s put ourselves in
the shoes of a hypothetical visitor, describing his or her experience:

> As I walked through the dusty streets of the lab’s main campus, I saw
> scientists and engineers hurrying back and forth, carrying stacks of papers
> and computer printouts. Many of them wore thick glasses and lab coats, and
> all of them seemed focused and serious, as if they had the weight of the
> world on their shoulders.
>
> Everywhere I looked, there were signs of the lab’s cutting-edge research. In
> one building, I glimpsed a group of scientists peering through microscopes,
> studying the structure of atoms. In another, I saw a group of engineers
> huddled around a massive machine, making adjustments to its gears and wires.
>
> But perhaps the most impressive sight was the lab’s newest creation: the
> Mathematical Analyzer, Numerical Integrator, and Computer, or MANIAC for
> short. This massive machine, which took up an entire room, was one of the
> fastest and most powerful computers in the world. Its thousands of vacuum
> tubes and punch-card programming made it an essential tool for the nuclear
> weapons lab’s work.
>
> Despite the seriousness of the lab’s work, there was also a palpable feeling
> of excitement and possibility in the air. Los Alamos was at the forefront of
> a new era of science and technology, and the scientists and engineers
> working there knew that their work had the potential to change the world in
> profound ways.
>
> As I left the lab and returned to the desert, I could not help but feel a
> sense of wonder at the incredible work being done behind those walls. Los
> Alamos was a place of innovation, determination and possibility, and I knew
> that its work would continue to shape the world for years to come.

Very suggestive. But even more suggestive to know that the “testimony” was
generated by asking the now-famous ChatGPT to narratively describe Los Alamos
in 1953. Apart from some naivete (such as studying the structure of atoms by
peering into a microscope of the time) what better way to emphasize the
upheaval that the advent of computers and software has brought in (only) 70
years? The world has really changed as profoundly as presumably the scientists
envisioned by artificial intelligence expected.

Even limiting ourselves to (our) specific field of Theoretical Physics the
progress is truly impressive: as every researcher knows well, the possibility
of simulating a model of given physical system is a very powerful tool. And we
are not only referring to the possibility to deal with an enormous amount of
data, but above all to devise Gedankenexperiments (ideal experiments) aimed at
understanding a phenomenon and perhaps, as each of us hopes in our daily work,
to discover a new one.

Well, all this originated precisely in the dusty streets of Los Alamos 70
years ago when, for the first time, people thought of using the impressive
computational tools as real conceptual laboratories. It is therefore with good
reason that one considers what took place in those years as the birth moment
of computational physics or of the experimental mathematics [1, 2]. And it is
in the year 1953 that two crucial steps take place: the publication of the
first article where the use of the Monte-Carlo technique is introduced and
demonstrated and the first calculations of the famous numerical experiment
conceived by Fermi, Pasta and Ulam.

On the occasion of this anniversary, we aim here to provide a brief historical
account of the experiment by Fermi and collaborators and an equally brief
update of the ideas that continue to flow from it. The relevance of the
problem and approach has been emphasized in so many publications, including
those directed to a wider audience [3, 2, 4]. We will therefore be concise,
referring to the literature for more specific issues. We begin by introducing
the protagonist of these scientific enterprises.

Figure 1: John von Neumann in front of MANIAC in 1952. The calculator was over
20 tons heavy and took up an entire room. It consisted of thousands of vacuum
tubes, which were used to perform logical operations, and was programmed using
punched cards.

##  2 MANIAC

John von Neumann, inspired by the work of Alan Turing, is considered to be the
mind behind the machines of the MANIAC family (fig. 1). According to
Metropolis [5] they were both convinced that

> modern computers could have profound effects in nonlinear mathematics by
> virtue of observing the unfolding of mathematical solutions that could
> stimulate thoughts about models and structures in those very refractory
> areas.

The internal architecture of today’s computers is, in fact, that conceived by
von Neumann. His intention was that MANIAC would initiate a computer
revolution, enabling the solution of scientific problems that had been deemed
impossible. To maximize the impact of such a program, he avoided applying for
patents and published detailed reports on its progress – a spirit that today
we would call of open science. The development of the first machine at the
Institute for Advanced Study in Princeton had already begun in 1946. Soon,
seventeen computers were built according to its specifications, including the
MANIAC in Los Angeles, identically named to that of Los Alamos, and IBM’s
first machine. Since it was a copy of the one in Princeton, the implementation
of the Los Alamos MANIAC I proceeded apace and, after a development phase, the
first run took place in March 1952 (fig. 2) [5].

Figure 2: MANIAC I in 1952 in color!

Despite its size and complexity, MANIAC was surprisingly fast, capable of
performing thousands of calculations per second. Obviously, by today’s
standards this is laughable performance. You have to think that the entire
computer memory (five kilobytes) is equivalent to an amount of memory less
than that occupied by a single icon on our laptops or cell phones. Also no one
had yet invented a modern programming language: the development of the first
FORTRAN compiler, the number crunching language par excellence, began shortly
thereafter, in 1954 [6].

Among its many uses, the machine seemed particularly suitable for solving of
differential equations that described novel problems, often (unfortunately)
motivated by wartime research. As Nicholas Metropolis, one of the leading
figures of that pioneering season, recalls [5]

> But the war was to change all that by demanding solutions to all kinds of
> nonlinear differential equations. At Los Alamos, from the very beginning in
> April 1943, it was clear that nonlinearities would come into their own — a
> new experience for the theorists who had learned to touch only linear
> phenomena and had a built-in shyness for anything beyond.

In this respect, the MANIAC manual states that
111www.bitsavers.org/pdf/lanl/LA-1725_The_MANIAC_Jul54.pdf

> The differential equations are of such complexity that analytical methods
> are not known for obtaining their solutions. The only recourse is to
> numerical procedures; therefore these problems are ideally suited for the
> computer. The first step in the solution of the problem is to replace the
> differential equations by a set of finite-difference equations. We do not
> discuss here the stability or convergence of such methods, but only mention
> them as necessary considerations in writing the difference equations.

So the new possibilities were truly remarkable, and the theorists soon
overcame their shyness toward nonlinear problems. One example is the Monte-
Carlo method which, although already known to statisticians, started with the
possibility of generating sets of pseudo-random numbers. And indeed, 1953 is
also the year of publication of the famous article ”Equation of State
Calculations by Very Fast Computing Machines” di N. Metropolis, A. W.
Rosenbluth e M. N. Rosenbluth, e M. Teller e E. Teller [7]. Its abstract
states that Results for the two‐dimensional rigid‐sphere system have been
obtained on the Los Alamos MANIAC and are presented here. Such work marked the
beginning of the use of the Monte-Carlo method for problem solving in the
physical sciences. The method described in this publication later became known
as the Metropolis algorithm, undoubtedly the most famous and most widely used
Monte Carlo algorithm, as evidenced by the impressive number of citations that
approaches 50000 to date. Interestingly, none of the authors of the paper made
subsequent use of the algorithm, and they remained unknown to the simulation
physics community that has developed since this publication. Their role in the
development of the algorithm became a subject of mystery and legend, as
recounted in [8]. The impact of Monte-Carlo methods on the physical sciences
was undoubtedly enormous, finding applications in all fields, from particle
physics to condensed matter.

Later, in 1956, MANIAC I had a further moment of glory to the general public
as well, becoming the first computer to beat a human being in a game of chess.
It was actually a simplified version of the game (a 6-by-6 board, without
bishops) known as Los Alamos Chess. To that result contributed also Stanislaw
(Stan) Ulam [9], one of the protagonists of the story we will now discuss.

##  3 The Los Alamos Report LA-1940

The results of the numerical experiment are reported in a now famous report
[10], which has never been published in a scientific journal. The work was
conceived in the summer of 1952 and the calculations were carried out in the
following summer. The model studied, which from the name of the authors,
Enrico Fermi, John Pasta and Stan Ulam (fig. 3) has long been referred to by
the acronym FPU, is presented as the discretization of a one-dimensional
continuum with the extremes fixed containing nonlinear terms. When the number
of equations tends to infinity, the system tends to the well-known wave
equation 222In this respect, it is curious to recall that in 1918 Fermi passed
brilliantly the entrance examination to the Scuola Normale in Pisa precisely
with a dissertation entitled Caratteri distintivi dei suoni where he writes si
trova l’equazione differenziale della verga vibrante, risolta per mezzo dello
sviluppo in serie di Fourier, dopo aver trovato le autofunzioni e gli
autovalori del problema [the differential equation for a vibrating bar is
found and solved by Fourier series, after having determined the eigenfunctions
and eigenvalues of the problem] [11].  with nonlinear terms of a complicated
nature. In the original notation, the equations thereby considered are,
denoting by $x_{n}$ the displacement

$\ddot{x}_{n}=x_{n+1}-2x_{n}+x_{n-1}+\alpha[(x_{n+1}-x_{n})^{2}-(x_{n}-x_{n-1})^{2}]$
(1)

where $n=1,\ldots,N$ ($N=64$ in the report) and

$\ddot{x}_{n}=x_{n+1}-2x_{n}+x_{n-1}+\beta[(x_{n+1}-x_{n})^{3}-(x_{n}-x_{n-1})^{3}],$
(2)

that are respectively called “model $\alpha$” and “model $\beta$”. Fixed
boundary are enforced letting $x_{0}=x_{N+1}=0$. The total mechanical energy
is conserved and the associated Hamiltonian has the general form

$H=\sum_{n=1}^{N}\frac{p_{n}^{2}}{2}+\sum_{n=0}^{N}U(x_{n+1}-x_{n})\,,$ (3)

where

$U(x)\equiv\frac{1}{2}x^{2}+\frac{\alpha}{3}x^{3}+\frac{\beta}{4}x^{4}.$ (4)

This formulation includes both model $\alpha$ ($\alpha\neq 0,\beta=0$) and
$\beta$ ($\alpha=0,\beta\neq 0$), but also the so-called “$\alpha\beta$ model”
($\alpha\neq 0,\beta\neq 0$), where both nonlinear terms are present.

Considering the case of fixed boundary conditions, the normal modes  of the
problem are defined by the transformations that diagonalize the quadratic part
of the Hamiltonian,

$Q_{k}=\sqrt{\frac{2}{N+1}}\sum_{i}x_{i}\sin\left(\frac{\pi
ki}{N+1}\right)\quad
P_{k}=\sqrt{\frac{2}{N+1}}\sum_{i}\dot{x}_{i}\sin\left(\frac{\pi
ki}{N+1}\right),$ (5)

and the main observables are the normal modes energies ($k=1,\ldots,N$)

$E_{k}=\frac{1}{2}P_{k}^{2}+\frac{\omega_{k}^{2}}{2}Q_{k}^{2},\quad\omega_{k}\equiv
2\sin\left(\frac{\pi k}{2(N+1)}\right).$ (6)

The numerical experiment consists in choosing initial conditions strongly out
of equilibrium, in which all the energy is assigned to one or few normal modes
of large wavelength, and in monitoring the energy distribution of the modes
over time. Fermi’s idea was to test on what time scale energy was transferred
from the macroscopic to the microscopic scale, reaching equilibrium. More
precisely, the state of equilibrium corresponds to the equipartition of energy
among the normal modes, that is, to the situation in which each mode, on
average, shares the same fraction of the total energy, initially attributed to
the initial condition. The question then is whether the time-averaged $E_{k}$
converges to the equilibrium value prescribed by the energy equipartition
principle and, if so, on what time scale.

Figure 3: Enrico Fermi, Jim Pasta, Stan Ulam and Mary Tsingou Mentzel.

###  3.1 “…features which were, from the beginning, surprising to us”

To the Authors’ surprise, on the time scale accessible to MANIAC, the
equilibrium state was not reached. The result is stated in the summary, which
concludes by summarizing the results concisely and effectively:

> The results show very little, if any, tendency toward equipartition of
> energy among the degrees of freedom.

Rather, a (quasi)steady state is observed, in which only some normal modes are
excited. In addition, the phenomenon of recurrence is observed. that is, the
energy of a mode returns to a value very close to its initial value, as
visible in the famous figure reproduced in fig.4. This apparent paradox has
been the source of many studies and we will discuss some recent developments
in a moment. But first we must talk about another protagonist in the story.

Figure 4: Figure 1 of the original report showing the time evolution of the
energies of the first 5 normal modes starting from a single sine-wave initial
condition.

###  3.2 “We thank miss Mary Tsingou…”

Another excerpt from the MANIAC user manual:

> There are several phases to the solution of a problem by an electronic
> computer. First, there is the formulation of the problem itself by the
> mathematician or theoretical physicist. Second, this is followed by the
> detailed preparation of the problem by the programmer for the specific
> computer.

It was therefore necessary for someone to translate the problem into a
flowchart and implement the design into a code usable by the machine. Fermi
himself had shown keen interest in programming already back in the 40s when,
with his characteristic approach, he self-taught the necessary techniques [5].
Indeed, as testified by Ulam himself [10], Fermi learned very quickly to
program MANIAC and worked personally on the flowcharts and the preparation of
the code. A footnote to the report thanks Mary Tsingou for efficient coding of
the problems and for running the computations on the Los Alamos MANIAC
machine. In the era when scientific journals demand to specify each author’s
contributions to an article and it is not considered rightly permissible to
exclude contributors, this seems unconceivable. Even more so considering that
we are talking about a young woman. 333As a matter of fact, several female
researchers were employed by the staff of MANIAC I, as listed by Table I in
[5].

In recent times, Mary Tsingou’s role has been widely recognized [12, 13], and
it has become common to replace the acronym FPU with FPUT. Although belated, a
due recognition to those who have spent a great deal of time and energy to
achieve the results, laboriously turning flowcharts of the code into punch
cards, with all the difficulties that anyone who had to deal with the
procedures of debugging can figure out.

In this regard, female researchers played a major role. On the occasion of
this anniversary it is appropriate to mention among them Arianna W.
Rosenbluth, considered to be the author of the implementation of the
Metropolis algorithm and the first person to implement the Monte-Carlo method
[7].

###  3.3 ”After the untimely death of Professor E. Fermi ..”

Let us come back to the report: at the end of the summary we learn that after
the Fermi’s passing in 1954, computational activity continued. What would
Enrico Fermi, the Papa as he was nicknamed in Via Panisperna, have thought,
from the height of his boundless knowledge of physics, of the developments of
his litte discovery? The profound influence of the ”never published article”
[3] cannot be underestimated. Just to give a rough idea of its significance
let us mention the following items:

  * 1.

As mentioned at the beginning, it demonstrates an innovative use of the
computer as a tool for investigation through what today we would call
numerical experiments.

  * 2.

It may be considered one of the starting points of the molecular dynamics, the
approach that allows the calculation of dynamical and statistical properties
based on numerical integration of the microscopic equations of motion. This
approach has numerous applications in physics, and is complementary to Monte-
Carlo methods, which originated in precisely the same context. 444For a
detailed and comprehensive historical approach to the birth and impact of this
important discipline on theoretical physics, see the recent volume [14].

  * 3.

Has motivated the study of deterministic chaos and its role in securing one of
the foundations of statistical mechanics namely the ergodic hypothesis,
bringing this problem back to the attention of physicists.

  * 4.

The attempt to explain the recurrence paradox led to the modern theory of
nonlinear wave equations (Boussinesq, KdV…) and of integrability for nonlinear
systems with many degrees of freedom.

We refer the reader to [15, 16, 17, 18] for an account of the subsequent
developments. In what follows we summarize some more recent ones, without
claiming to be complete, with the aim of showing that this is a still an
active field of research .

##  4 The shadow of Toda

From the very beginning, thanks to the contribution of Zabusky and Kruskal, it
had been assumed that the phenomenon of recurrences originated from the quasi-
integrability of the model. More recently, a new point of view has emerged
[19]. The new viewpoint is actually based on an old knowledge, namely the
Hamiltonian of the form (3) with the exponential potential

$U_{Toda}=U_{0}\left(e^{\lambda x}-1-\lambda x\right).$ (7)

Morikazu Toda formulated it in 1967 (and has since then beared his name) and
determined its analytical solutions. Subsequently [20, 21] it was proven to
have a special property: it is one of the rare examples of a nonlinear
integrable system, which have as many constants of motion as there are degrees
of freedom. As is known from analytical mechanics, by means of a canonical
transformation to action-angle variables $I,\varphi$, the corresponding
Hamiltonian can be rewritten as a function of actions only, $H(I)$. Thus, the
dynamics turns out to be quasi-periodic, remaining confined on a multi-
dimensional torus.

The main point can be summarized simply: Toda’s model represents an integrable
approximation of FPUT more accurate than the harmonic chain. Thus, the slow
chaotic motion leading to equipartition is not so much due to the fact that
FPUT is a discretization of an integrable wave equation, but rather to the
fact that the dynamics is (at least at low enough energies) essentially
indistinguishable from Toda’s on very long time scales (fig. 5). This
viewpoint has been established only recently and represents a significant
paradigm shift in explaining the problem, even in quantitative terms through
scaling laws [22]. Furthermore, the metastable state can be seen as a
generalized Gibbsian state, the generalization to the integrable case of the
canonical state of usual statistical mechanics [23].

Figure 5: The motion of an integrable system, such as the Toda model, occurs
on an invariant torus in which the actions are constant and the angles grow
linearly with time. The presence of terms that break integrability induces a
slow motion in the transverse direction, but mantain the trajectories close to
the torus. On this time scale the dynamics is practically indistinguishable
from that of the underlying integrable system: This results in a metastable,
non-ergodic state (Courtesy of G. Benettin and A. Ponno).

##  5 Further developments

The amount of research problems generated by the FPUT numerical experiment is
so vast that it cannot be completely accunted for in this brief article. In
this section we report a selection of some of these issues, referring readers
to the cited literature.

###  5.1 The Anti-FPUT problem

The original numerical experiment considered initial conditions focused on the
normal modes of large wavelength. What happens if instead we concentrate the
initial energy on those at the other end of the spectrum, namely the modes of
wavelengths comparable to the lattice spacing? In this case, through a
mechanism known as modulational instability, localized excitations are
spontaneously formed that correspond to so-called discrete breathers or, in
the limit of small amplitudes, envelope solitons [24, 25]. The thermalization
process is determined by the interaction of this ”gas” of breathers that
collide and exchange energy in a peculiar way: the excitation tends to
concentrate over time in a single breather that is eventually destroyed by the
fluctuations. The system relaxes to equipartition, on time scales that
increase as energy density decreases.

###  5.2 The mathematics of FPUT

In 1954, a year before the Los Alamos report, Kolmogorov announced in
Amsterdam one of the most important results in twentieth-century mathematics.
That theorem, now known as KAM (from the names of the three mathematicians
credited with its authorship: Kolmogorov himself, Arnold and Moser), provides
sufficient conditions for the ”integrable” behavior of a Hamiltonian system to
survive for all times even when the system is perturbed. However, the
applicability of the KAM Theorem to the FPUT model has only recently been
demonstrated. The first positive result, valid only for specific chain lengths
and only for the $\beta$ model, was due to Nishida in 1971. It took another 30
years before the applicability of this result was proved first for the $\beta$
model [26] and then also for the other models [27] . As a consequence of
applying the KAM theorem, we now know that, when the number of particles in
the chain is fixed, in a low-energy regime the system will never reach
thermodynamic equilibrium.

Following the route opened by Zabusky and Kruskal regarding the link between
the Korteweg de Vries (KdV) equation and FPUT, the first mathematically
rigorous results date back to the early 2000s. In [28] it was shown that in an
appropriate continuous limit and after a ”normal form” step, the FPUT model is
described by a pair of KdV equations describing waves traveling in opposite
directions. As for the connection with Toda’s model, the first to realize it
was Manakov in 1974, followed by a very important work by Flashka, Ferguson
and McLaughlin [29]. But to arrive to two theorems linking the length of
equipartition times in FPUT and Toda dynamics we have to wait for the work
[30]. Here, the authors show that if the energy of the system is low enough,
the Toda actions are very close to the normal modes of the chain. Thus the
quasi-conservation of Toda actions implies lack of equilibration for the FPUT
model on the timescale for which the approximation holds, in agreement with
the numerical results. Only very recently in [31] and in [32] the authors
manage to show that, if the temperature is low enough, Toda actions are
adiabatic invariants for the FPUT system even in the thermodynamic limit.

Finally, we mention another approach to the FPUT problem based on the wave-
kinetic (or wave turbulence) method that has been recently proposed, see [33]
for an account.

###  5.3 Almost-integrability and transport

The quasi-periodic recurrence of initial states discovered by FPUT not only
challenged spontaneous evolution toward thermodynamic equilibrium, but was
also incompatible with the presence of any diffusive mechanism for energy
transport.

One might therefore wonder what happens if the FPUT model is subjected to an
external thermodynamic force, such as a temperature gradient. In a normal
transport regime, we would expect it to to behave like an ordinary thermal
conductor for which the well-known Fourier’s law states that the heat flow
passing through it is proportional to the temperature gradient itself.
Instead, it turns out that this relationship is not verified and that the
system behaves as a ”heat superconductor”, in the sense that the thermal
conductivity diverges in the thermodynamic limit [34]. This phenomenon of so-
called anomalous transport, although not directly related to the slow
thermalization described above, represents another surprise of the model and
has been studied in detail, even experimentally [35].

But the problem becomes even more intriguing in the limit of low temperatures
where FPUT dynamics is, as described, well approximated by the Toda chain. In
the latter, solitons propagate freely and give rise to ballistic transport,
characteristic of independent quasi-particle. From the point of view described
in section 4, that is, considering the FPUT model as a weak perturbation of
Toda, near the integrable limit, the quasi-particles interact and acquire a
finite but large mean free path $\ell$. Therefore, if the chain size $L$ is
such that $L<\ell$ purely ballistic behavior is still observed, whereas if
$L>\ell$ the solitons diffuse giving rise to normal transport. Thus Toda’s
shadow manifests again, restoring Fourier’s law at least in a regime of
intermediate sizes. But this does not last forever: at even higher $L$ sizes,
one recovers the anomalous behavior mentioned above [36], and the effects of
quasi-integrability are lost in the chaos of hydrodynamic fluctuations.

##  6 Dèjà vu in the lab

There is also some experimental evidence of the physics of FPUT, particularly
of the phenomenon of recurrence [37, 38] which, suggestively, has been defined
as the analogue of the Déjà vu in the inanimate world [39]. Nonlinear optical
systems are good examples. For example, the propagation of the electric field
in an optical fiber is accurately described by the nonlinear Schrödinger
equation.

$i\psi_{z}+\psi_{xx}+2\psi|\psi|^{2}=0$ (8)

ruling the evolution of the envelope $\psi(x,t)$ of a quasi-monochromatic
field in a medium with weak optical nonlinearity ($z$ represents the direction
of propagation). The equation is fully integrable and also has solitonic
solutions. In addition, the plane wave-like solutions are modulationally
unstable, and if the perturbation consists, for example, of a single mode,
after growing to a certain maximum amplitude, it reduces in size and
eventually disappears. From this point of view, the recurrence of FPUT can be
seen as is a periodic solution of the wave equation (8) [39]. On the basis of
these considerations, several experimental observations of the recurrences in
various optical systems have been demonstrated [40, 41, 42].

##  7 Conclusion

A first conclusion we would like to emphasize is how models, however
mathematically simple, can lead to new and surprising discoveries and, as in
the case of FPUT, even open up new fields of research, even interdisciplinary
ones.

Moreover, before concluding, it seems to us that this is the appropriate place
to emphasize the contribution of the Italian physics community to the FPUT.
Even from the brief bibliography provided here, it is clear that many of the
works on the topic, since its inception, are signed by Italian researchers
(Rome, Milan, Padua, Florence) who continue to be passionate about the problem
and provide new advances, both theoretical and experimental.

Finally, as we have seen, the many current developments in the problem
originate from the creative use of the computer as a tool for investigation.
We began this article with the help of the most modern and powerful
algorithmic resource of today. Indeed, we are witnessing a new and impressive
technological breakthrough that will surely impact on the methods of
scientific research, as well as on the rest of the society:

> […] there is now a new and very real possibility for a novel synthesis of
> advances in experimental technique, high-performance computing, and theory:
> automatically building theories from data. That is,to the extent we
> understand pattern, we can use machines to find emergent organization in our
> vast data sets. And from there, the machines can build our theories, most
> likely with guidance from a new generation of theorists. [43]

Understanding the potential and methodological implications offered by such an
evolution, in a creative and original approach as happened at Los Alamos 70
years ago, could pave the way for new and unexpected scientific advances.

## References

  * [1] T. Dauxois, M. Peyrard, S. Ruffo, The Fermi–Pasta–Ulam ‘numerical experiment’: history and pedagogical perspectives, European Journal of Physics 26 (5) (2005) S3.
  * [2] M. A. Porter, N. J. Zabusky, B. Hu, D. K. Campbell, Fermi, Pasta, Ulam and the birth of experimental mathematics: a numerical experiment that Enrico Fermi, John Pasta, and Stanislaw Ulam reported 54 years ago continues to inspire discovery, American Scientist 97 (3) (2009) 214–221.
  * [3] M. Falcioni, A. Vulpiani, Il contributo di Enrico Fermi ai sistemi non lineari: L’influenza di un articolo mai pubblicato, Conoscere Fermi nel centenario della nascita 29 (2001).
  * [4] J. G. De Tullio, Il modello Fermi-Pasta-Ulam: La nascita della simulazione numerica, Lettera Matematica Pristem 96 (1) (2016) 34–40.
  * [5] N. Metropolis, The Los Alamos experience, 1943–1954, in: A history of scientific computing, 1990, pp. 237–250.
  * [6] W. B. Maynard, Daybreak of the digital age, Princeton Alumni Weekly 112 (10) (2012) 28–33.
  * [7] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, Equation of state calculations by fast computing machines, The journal of chemical physics 21 (6) (1953) 1087–1092.
  * [8] J. E. Gubernatis, Marshall Rosenbluth and the Metropolis algorithm, Physics of plasmas 12 (5) (2005) 057303.
  * [9] J. Kister, P. Stein, S. Ulam, W. Walden, M. Wells, Experiments in chess, Journal of the ACM (JACM) 4 (2) (1957) 174–177.
  * [10] E. Fermi, J. Pasta, S. Ulam, Studies of nonlinear problems, Los Alamos Report LA-1940 (1955) 978.
  * [11] E. Segrè, Enrico Fermi, fisico. Una biografia scientifica, Plunkett Lake Press, 2020.
  * [12] T. Dauxois, Fermi, Pasta, Ulam and a mysterious lady, Physics Today 61 (1) (2008) 55–57.
  * [13] T. Dauxois, La (dis-) simulation FPU: une Femme Physicienne Underground-Mary Tsingou et l’histoire de la première expérience numérique, Reflets de la physique (67) (2020) 26–30.
  * [14] G. Battimelli, G. Ciccotti, P. Greco, Computer Meets Theoretical Physics: The New Frontier of Molecular Simulation, Springer Nature, 2020.
  * [15] J. Ford, The Fermi-Pasta-Ulam problem: paradox turns discovery, Physics Reports 213 (5) (1992) 271–310.
  * [16] A. Carati, L. Galgani, A. Giorgilli, A. Ponno, The Fermi-Pasta-Ulam problem, Il Nuovo Cimento B 117 (2002) 1017–1025.
  * [17] G. Gallavotti, The Fermi-Pasta-Ulam problem: a status report, Vol. 728, Springer, 2007.
  * [18] D. K. Campbell, P. Rosenau, G. M. Zaslavsky, Focus issue: the Fermi–Pasta–Ulam problem—the first fifty years, Chaos: An Interdisciplinary Journal of Nonlinear Science 15 (1) (2005) 015101.
  * [19] G. Benettin, H. Christodoulidi, A. Ponno, The Fermi-Pasta-Ulam problem and its underlying integrable dynamics, Journal of Statistical Physics 152 (2) (2013) 195–212.
  * [20] M. Henón, Integrals of the Toda lattice, Physical Review B 9 (4) (1974) 1921\.
  * [21] H. Flaschka, The Toda lattice. II. Existence of integrals, Physical Review B 9 (4) (1974) 1924.
  * [22] G. Benettin, A. Ponno, Understanding the FPU state in FPU-like models., Mathematics in Engineering 3 (3) (2020) NA–NA.
  * [23] T. Goldfriend, J. Kurchan, Equilibration of quasi-integrable systems, Physical Review E 99 (2) (2019) 022146.
  * [24] Y. A. Kosevich, S. Lepri, Modulational instability and energy localization in anharmonic lattices at finite energy density, Physical Review B 61 (1) (2000) 299\.
  * [25] T. Dauxois, R. Khomeriki, F. Piazza, S. Ruffo, The anti-FPU problem, Chaos: An Interdisciplinary Journal of Nonlinear Science 15 (1) (2005) 015110.
  * [26] B. Rink, Proof of Nishida’s conjecture on anharmonic lattices, Communications in mathematical physics 261 (3) (2006) 613–627.
  * [27] A. Henrici, T. Kappeler, Results on normal forms for FPU chains, Communications in mathematical physics 278 (1) (2008) 145–177.
  * [28] D. Bambusi, A. Ponno, On metastability in FPU, Communications in mathematical physics 264 (2006) 539–561.
  * [29] W. Ferguson Jr, H. Flaschka, D. McLaughlin, Nonlinear normal modes for the Toda chain, Journal of computational physics 45 (2) (1982) 157–209.
  * [30] D. Bambusi, A. Maspero, Birkhoff coordinates for the Toda lattice in the limit of infinitely many particles with an application to FPU, Journal of Functional Analysis 270 (5) (2016) 1818–1887.
  * [31] A. Maiocchi, D. Bambusi, A. Carati, An averaging theorem for FPU in the thermodynamic limit, Journal of Statistical Physics 155 (2) (2014) 300–322.
  * [32] T. Grava, A. Maspero, G. Mazzuca, A. Ponno, Adiabatic invariants for the FPUT and Toda chain in the thermodynamic limit, Communications in Mathematical Physics 380 (2) (2020) 811–851.
  * [33] M. Onorato, Y. V. Lvov, G. Dematteis, S. Chibbaro, Wave turbulence and thermalization in one-dimensional chains, arXiv preprint arXiv:2305.18215 (2023).
  * [34] S. Lepri, R. Livi, A. Politi, Heat conduction in chains of nonlinear oscillators, Physical review letters 78 (10) (1997) 1896.
  * [35] G. Benenti, D. Donadio, S. Lepri, R. Livi, Non-Fourier heat transport in nanosystems, La Rivista del Nuovo Cimento 46 (3) (2023) 105–161.
  * [36] S. Lepri, R. Livi, A. Politi, Too close to integrable: Crossover from normal to anomalous heat diffusion, Phys. Rev. Lett. 125 (2020) 040604.
  * [37] G. Van Simaeys, P. Emplit, M. Haelterman, Experimental demonstration of the Fermi-Pasta-Ulam recurrence in a modulationally unstable optical wave, Physical review letters 87 (3) (2001) 033902.
  * [38] M. Wu, C. E. Patton, Experimental observation of Fermi-Pasta-Ulam recurrence in a nonlinear feedback ring system, Physical review letters 98 (4) (2007) 047202\.
  * [39] N. N. Akhmediev, Déjà vu in optics, Nature 413 (6853) (2001) 267–268.
  * [40] D. Pierangeli, M. Flammini, L. Zhang, G. Marcucci, A. Agranat, P. Grinevich, P. Santini, C. Conti, E. DelRe, Observation of Fermi-Pasta-Ulam-Tsingou recurrence and its exact dynamics, Physical Review X 8 (4) (2018) 041017.
  * [41] J.-W. Goossens, H. Hafermann, Y. Jaouën, Experimental realization of Fermi-Pasta-Ulam-Tsingou recurrence in a long-haul optical fiber transmission system, Scientific reports 9 (1) (2019) 18467.
  * [42] G. Vanderhaegen, P. Szriftgiser, A. Kudlinski, M. Conforti, S. Trillo, M. Droques, A. Mussot, Observation of four Fermi-Pasta-Ulam-Tsingou recurrences in an ultra-low-loss optical fiber, Optics Express 28 (12) (2020) 17773–17781.
  * [43] J. P. Crutchfield, The dreams of theory, Wiley Interdisciplinary Reviews: Computational Statistics 6 (2) (2014) 75–79.


2022

[1]Jared Grove

1]Mathematics, University of Iowa, 2 West Washington St, Iowa City, 52242,
Iowa, United States

2]Computer Science, University of Iowa, 2 West Washington St, Iowa City,
52242, IA, United States

3]Metric Geometry and Gerrymandering Group, Tufts University, 10 Upper Campus
Rd, Medford, 02155, MA, United States

# A $k$-medoids Approach to Exploring Districting Plans

<EMAIL_ADDRESS>Dr. Suely P. Oliveira<EMAIL_ADDRESS>Anthony
Pizzimenti<EMAIL_ADDRESS>Dr. David Stewart<EMAIL_ADDRESS>[  [
[

###### keywords:

Redistricting, Gerrymandering, Algorithms, Computer Districting

##  1 Introduction

###  1.1 Districting Plans

Every $10$ years the United States goes through a reapportionment process,
where states may gain or lose seats in the House of Representatives based on
population changes. States with more than one representative may need to
redistrict, or redraw the lines for their congressional districts based on how
the populations have shifted within the state or because the number of
representatives for that state changed.

When a state redistricts, lines must be drawn to divide the state into $k$
districts, one for each Congressional seat. This can be represented
mathematically by using a graph or network structure, where the state is a
graph $G=(V,E)$ where $V$ is the set of vertices and $E$ is the set of edges.
Each state is composed of many voter tabulation districts (VTDs) that need to
be assigned to districts, and some states require that larger existing
political structures remain intact. For example, Iowa requires that counties
are not split apart. We will refer to the basic building blocks of a
districting plan as tabulation blocks (TBs), which can represent counties,
VTDs, precincts, or any combination of them. If there are $N$ TBs, they will
be represented in $G$ as vertices $\left\\{\,v_{i}\mid
i=1,\ldots,N\,\right\\}=V$. The edges in $G$,
$\left\\{v_{i},\,v_{j}\right\\}\in E$, will represent a shared geographical
border between vertices $v_{i}$ and $v_{j}$. As an example of county level
TBs, the state map and graphical representation of Iowa can be seen in Figure
1.

(a) (b)

Figure 1: County map of Iowa given in (1(a)) and graphical representation of
Iowa given in (1(b))

To build a districting plan, each vertex $v_{j}$ needs to be assigned to a
district $v_{j}\in D_{i}$ in the districting plan $\left\\{\,D_{i}\mid
i=1,\ldots,k\,\right\\}$ which is a partition of the TBs. The TBs must be
assigned to districts such that the districts fit defined criteria for
contiguity, compactness, and population equality. The criteria for contiguity
requires a district to be a single connected region without any holes. The
criteria for compactness is less strictly defined, in that a district is
considered compact when it covers the smallest possible geographic region.
This is often interpreted to mean that a district should be as close to square
or circular as possible. A variety of compactness measures have been used
previously, including Polsby-Popper Polsby , Reock Reock , Schwartzberg
Schwartz , Convex Hull RestOfComp , X-Symmetry RestOfComp , and Length-Width
Ratio LWR . All listed measures return a compactness value between $0$ and $1$
with $1$ being considered compact and $0$ being not compact. Non-compact
districts are often long and winding with a bizarre shape Bizarre , and are
often rejected because non-compactness is a hallmark of gerrymandered
districts. In this work we will be using the Convex Hull as our compactness
measure as it has been used in many court cases chen_rodden_2015  and, as
shown in Article4-Jumble , is less affected by choice of map projection. The
Convex Hull measure is defined by the ratio of district area to the area of
the minimum convex polygon bounding the district. We will report the
compactness value as the mean of the compactness values of each district in
the districting plan. Lastly, population equality requires that each district
has approximately equal populations. Thus an ideal district would have a
population equal to the total population of the state divided by the number of
districts $k$:

$\mathit{Pop}^{*}=\frac{1}{k}\sum_{i=1}^{N}\mathit{Pop}(v_{i}),$ (1)

where $\mathit{Pop}(v_{i})$ is the population of the TB represented by vertex
$v_{i}$. In order to evaluate a districting plan we will compare the
population of each district, $\mathit{Pop}(D_{i})=\sum_{j:v_{j}\in
D_{i}}\mathit{Pop}(v_{j})$ to $\mathit{Pop}^{*}$. The deviation of a
districting plan will be the maximum absolute percent deviation of the
districts:

$\mathit{Dev}=\max_{i\in\\{1,\ldots,k\\}}\left(\frac{\lvert\mathit{Pop}(D_{i})-\mathit{Pop}^{*}\rvert}{\mathit{Pop}^{*}}\cdot
100\%\right)$ (2)

We aim to create an algorithm to generate districting plans with
$\mathit{Dev}<1$, which is equivalent to a districting plan that has maximum
population deviation of less than $1\%$. We will do this only using the
geographic data and population information.

###  1.2 Gerrymandering

While redistricting is a necessary process that can allow the people to be
properly represented in Congress, it can also be manipulated by the agents
that draw the lines. This is known as gerrymandering and can be done to retain
political power, take power away from opponents, or to further various other
political goals. In the worst cases, political officials are essentially able
to select their voters, rather than voter selecting their officials.

The most obvious gerrymanders are often identified by the bizarreness of their
shape Bizarre , for example long narrow districts that wind throughout the
state. These types of districts are often easily identifiable using
compactness measures. However, some gerrymanders are much less obvious to the
eye and can be much harder to identify. To avoid this problem various
algorithms have been proposed to create districting plans, effectively
removing politicians and their agents from the redistricting process. However,
algorithm based districting does not necessarily create politically neutral
districts as the geographic and population data can correlate with political
affiliation.

Generating districting plans using methods that reduce political and other
bias can also help identify gerrymandering EnsambleOutliers . A gerrymandered
districting plan would have political outcomes (that is, election outcomes)
well outside the range of election outcomes for districting plans created
without political bias. Generating large numbers of districting plans without
political bias, and then determining the expected election outcomes for these
generated plans shows the range of outcomes to be expected where there is no
gerrymander. Expected election outcomes for a given districting plan well
outside this range then becomes strong evidence of a gerrymander

###  1.3 Other Approaches

The majority of computer-based districting approaches are either sampling or
optimization based methods.

This is because many of the formulations for this type of problem have been
shown to be NP-Complete (Subset Sum computersAndIntractibility ).

Sampling based methods are generally concerned with generating many plans
quickly. Often times these are flood fill or multi-kernel methods that start
by randomly selecting several initial centers to iteratively build districts
from. Typically these will completely build one district, then the second, and
so on until there are the proper number of districts. This process has been
discussed and implemented many times vickrey_1961-Sampling ; GROUP-Sampling ;
SouthCarolina-Sampling ; thoreson_liittschwager_1967-Sampling . A potential
issue with this method is that it can create enclaves, groups of TBs that
cannot be added to any district, and implementations have to monitor this and
will restart if enclaves appear. A further variation of the basic sampling
method is mentioned in Article11-Algorithms , growing the districts in
simultaneously instead of sequentially. This should keep the shape of the
districts more consistent compared to the sequential version. A similar method
chen_rodden_2015 , considers all TBs as their own district and combines them
until the correct number of districts has been achieved.

The second family of approaches is based on optimization, and seeks to find
the best plan under some objective function. These can be roughly subdivided
into two more categories; clustering and heuristics.

Clustering algorithms typically follow a $k$-means or warehouse location
approach to assign TBs to cluster centers while minimizing the sum of squares
distances from the TBs to the assigned cluster center. Several approaches
using these techniques can be found, weaver_hess_1963-KMEANS ;
guest_kanayet_love_2019-KMEANS ; cohen-addad_klein_young_2018-KMEANS ;
Article1-Heuristic . When following a clustering approach special care needs
to be taken to ensure that the dividing lines do not break up the TBs that
make up the state; some of these legal technicalities are ignored by certain
implementations of these methods (see, for example,
guest_kanayet_love_2019-KMEANS ; cohen-addad_klein_young_2018-KMEANS ).

Heuristic algorithms often involve using a mixed integer linear program to
model districting and find an optimal solution according to some criteria,
often minimizing population deviation. They can be solved with various local
searches nagel_1965-LS ; kaiser_1966-LS ; ricca_simeone_2008-LS , tabu search
BOZKAYA200312-TABU ; ricca_simeone_2008-LS , genetic algorithms PEAR ,
simulated annealing ricca_simeone_2008-LS ; fifield_higgins_imai_tarr_2020-SA
, and many more heuristic approaches. While these can do a great job of
finding plans according to their objective function, there is always the issue
of agreeing on what the plans should optimize for: should they be focused on
minimizing cut edges in the graph, getting the populations exactly equal, or
some other criteria?

A third type of algorithm that seems to have less representation in the
literature uses random walks over the set of districting plans
fifield_higgins_imai_tarr_2020-SA ; Gerrychain . Random walks move from one
districting plan to another in an attempt to understand what the solution
space of districting plans looks like. These types of methods need many good
quality initial seeds in order to effective sample the solution space.
Combining them with either the sampling or optimization methods can provide
these random walk methods with the starting points they need to explore the
space Gerrychain .

###  1.4 Our Approach

We use a $k$-medoids algorithm to approach this problem. It is a sampling
method that borrows ideas from some of the optimization methods. We begin from
the graph representation of the state and follow the multikernel procedure
mentioned in Section 1.3. However, instead of following the sequential
building pattern used by most implementation we loosely follow an alternating
building pattern mentioned in Article11-Algorithms . The basics of the
$k$-medoids approach is to start with $k$ TBs as initial district centers,
then alternate between assigning TBs to the nearest centers and computing the
centers of districts. It is similar to the $k$-means method for clustering
weaver_hess_1963-KMEANS ; guest_kanayet_love_2019-KMEANS ; cohen-
addad_klein_young_2018-KMEANS ; Article1-Heuristic .

However, instead of using the physical distances between units, we use graph-
based distances. Furthermore, in $k$-medoids, the centers, or medoids, are
required to be vertices on the graph as opposed to any point in space.
Following the $k$-medoids process we move into a local search phase. We use
the traditional local search strategies, flip and swap neighborhoods, along
with a tabu criteria to fine-tune solutions. In states with many TBs we also
implement a coarsening strategy to work with fewer TBs, then follow an
uncoarsening schedule to allow for more fine-tuning as better solutions are
found. To the authors’ knowledge, the first publication of this
coarsening/uncoarsening strategy for districting plans is Magkeby and
Mosessoon magleby_mosesson_2018-Coarsening . We explore the implemented
algorithms in more detail in Section 2.

##  2 Algorithms

The base of the algorithm described here is a $k$-medoids algorithm, which
alternates between assigning points to clusters based on proximity to a medoid
and computing the medoids of the clusters. However, in the most basic form,
$k$-medoids algorithms fail to create districts with equal populations.
Several additional steps are added to the basic $k$-medoids framework to
account for this issue and improve overall performance. Before describing the
algorithm as a whole, it is first useful to walk through these additional
steps. These include various forms of local search as well as graph coarsening
and uncoarsening.

###  2.1 Local Search

With the $k$-medoids process used here, deviation does not strictly decrease,
so finding a districting plan with a low deviation is not guaranteed.

To increase the likelihood of finding a good solution we look to fine tune the
output of the $k$-medoids algorithm by making small changes to the resultant
districting plans. These small changes can be thought of as finding
neighboring districting plans. Methods of identifying neighboring districting
plans are known as neighborhood functions.

We consider three types of neighborhood functions - Flip, Swap, and a
Combination Search (CMB). The Flip neighborhood function generates all
districting plans that have exactly one TB assigned to a different district.
The Swap neighborhood function generates all districting plans that have
exactly two TBs that have mutually changed districts. The CMB neighborhood
function generates all districting plans that are in either the Flip
neighborhood or Swap neighborhood. In general, we will abbreviate an arbitrary
neighborhood function as $\mathit{NF}$.

An issue that arises with these neighborhood functions is that the number of
neighboring districting plans generated can be very large, making it
infeasible to enumerate them and find the neighboring districting plan that
would have the lowest deviation.

To get around this, we only generate neighboring districting plans according
to specific criteria. From the initial districting plan we consider all pairs
of neighboring districts $\left\\{D_{m},D_{\ell}\right\\}$ with
$\mathit{Pop}(D_{m})>\mathit{Pop}(D_{\ell})$. From these we find the pair of
districts with the largest population disparity. We use these districts
differently for each method:

  * Flip:

Find the set of all vertices $v_{i}\in D_{m}$ that share a geographical border
with any vertex in $D_{\ell}$ and that are not articulation points of the
graph of $D_{m}$; articulation points are vertices that when removed from a
graph increase the number of connected components. This will be stored as
$\mathit{FLIP}=\left\\{\,(v_{i},v_{-1})\,\right\\}$ where $v_{-1}$ is a dummy
vertex with $\mathit{Pop}(v_{-1})=0$.

  * Swap:

Find the set of all pairs of vertices
$\mathit{SWAP}=\left\\{\,(v_{i},v_{j})\mid v_{i}\in D_{m},\,v_{j}\in
D_{\ell}\,\right\\}$ such that neither $v_{i}$ nor $v_{j}$ are articulation
points of their respective districts, $v_{i}$ shares a geographical border
with a vertex in $D_{\ell}\backslash\\{v_{j}\\}$, and $v_{j}$ shares a
geographical border with a vertex in $D_{m}\backslash\\{v_{i}\\}$.

  * CMB:

$\mathit{CMB}=\mathit{FLIP}\,\cup\,\mathit{SWAP}$. The set of pairs of
vertices that satisfy either Flip or Swap.

The sets of pairs of vertices derived from the above rules represent the
districting plans generated according to the corresponding neighborhood
function, which we will denote $\mathit{NH}$. If no districting plans were
generated, we select the pair of neighboring districts with the next highest
population disparity until there is at least one districting plan generated.

Once a districting plan is generated, the local search considers all generated
districting plans and selects the one with the minimum population disparity.
This is done by finding

$\displaystyle\begin{split}(v^{*}_{i},v^{*}_{j})=\operatorname*{arg\,min}_{(v_{i},v_{j})\in\mathit{NH}}\max\left(\,\left\lvert\,\frac{[\mathit{Pop}(D_{m})-\mathit{Pop}(v_{j})+\mathit{Pop}(v_{j})]-\mathit{Pop}^{*}}{\mathit{Pop}^{*}}\,\right\rvert,\right.\\\
\left.\left\lvert\,\frac{[\mathit{Pop}(D_{\ell})+\mathit{Pop}(v_{i})-\mathit{Pop}(v_{j})]-\mathit{Pop}^{*}}{\mathit{Pop}^{*}}\,\right\rvert\,\right)\end{split}$
(3)

This pair of vertices represents the districting plan with the minimum
population disparity.

To construct the districting plan we move $v_{i}$ to $D_{\ell}$ and $v_{j}$ to
$D_{m}$, if $j$ is not the dummy node $-1$. This can be seen in Figure 2 for
Flip and Figure 3 for Swap.

(a) (b) (c) (d)

Figure 2: Example of Flip neighborhood function. First select Districts that
Flip will be performed on, $D_{m}$:  and $D_{\ell}$:  (2(a)), then find
candidate vertices in $D_{m}$:  that can be flipped (2(b)). Next find
candidate vertex that minimizes Equation (3) (2(c)) and Flip minimizing vertex
from $D_{m}$ to $D_{\ell}$ (2(d))

(a) (b) (c) (d)

Figure 3: Example of Swap neighborhood function. First select Districts that
Swap will be performed on, $D_{m}$:  and $D_{\ell}$:  (3(a)), then find pairs
of candidate vertices in $D_{m}$:  and $D_{\ell}$:  that can be swapped
(3(b)). Next find candidate vertex pair that minimizes Equation (3) (3(c)) and
Swap minimizing vertex pair from (3(d))

Starting from this new districting plan, we repeat the process of finding the
pair of neighboring districts with the largest disparity, generating
neighboring districting plans, finding ($v_{i}^{*}$, $v_{j}^{*}$) according to
Equation (3), and changing districting plans until a districting plan with
$\mathit{Dev}<1$ is found or until $\mathit{LI}$ iterations of local search
are reached, where $\mathit{LI}$ is a user defined maximum number of allowed
local search iterations.

We noticed that these neighborhood functions could cycle one TB between a few
districts, getting stuck in a loop. To get around this we also included a tabu
list in the search process, with tabu tenure set as $0.1\times\mathit{LI}$.
This means that after each successful neighborhood move is made, we store the
move in the tabu list, and cannot repeat it for $0.1\times\mathit{LI}$
iterations. The entire Local Search process is presented in Algorithm 1.

Algorithm 1 Local Search

1:Choose $\mathit{LI}$ and $\mathit{NF}$

2:while $\mathit{Dev}>1$ and Iterations $<\mathit{LI}$ do

3: Store pairs of neighboring districts and population disparity between them

4: while a valid move has not been found do

5: Set $\\{D_{m},D_{\ell}\\}$ as pair of districts with largest population
disparity

6: Find the set of valid, non-tabu moves in $\mathit{NH}$ according to
$\mathit{NF}$

7: if there is a valid move then

8: Make move that minimizes Equation (3)

9: Make move tabu for $0.1\times\mathit{LI}$ iterations

10: else

11: Remove $\\{D_{m},D_{\ell}\\}$ from the pairs of neighboring districts

12: end if

13: end while

14: Increment Iteration and compute $\mathit{Dev}$

15:end while

###  2.2 Coarsening and Uncoarsening

There are $27$ states with counties whose population is larger than
$\mathit{Pop}^{*}$ for that state. Thus, these states cannot be districted at
the county level and instead must be broken down further to the VTD level.
This causes the number of TBs to grow very large, and can greatly slow down
our method. For example, with Florida there are $67$ counties and $9435$ VTDs.
The original graph has $N$ vertices to represent the TBs. The user can provide
a parameter $uc_{0}\in\left(0,1\right]$ in order to coarsen $G$ to a graph
$G_{0}$ which has $\lfloor\,uc_{0}N\,\rfloor$ vertices. To reduce the number
of vertices, the coarsening process randomly selects two neighboring vertices.
If the sum of the populations of these two vertices is less than
$\mathit{Pop}^{*}$, the two neighboring (parent) vertices are combined,
resulting in a new (child) vertex.

The child vertex will inherit edges from the union of the edges of the parent
vertices so that the child vertex has the same neighbors as the parent
vertices. The population of the child vertex is the sum of the parents’
populations. As we coarsen the graph, we store the parent vertices, their
populations, edges (or neighbors), and the order they are combined in, to
ensure we can reconstruct the original graph, $G$. This coarsening process is
repeated until there are only $\lfloor\,uc_{0}N\,\rfloor$ vertices in the
graph. This process can be seen in Figure 4.

(a) (b) (c) (d) (e)

Figure 4: Coarsening Example with $uc_{0}=0.6$. Beginning with graph $G$ which
has $N=10$ vertices (4(a)) and $uc_{0}=0.6$, coarsen to $G_{0}$ with
$\lfloor\,uc_{0}N\,\rfloor=6$ vertices. Randomly select pairs of neighboring
vertices and combine them until there are only $6$ vertices remaining. First
vertices $v_{4}$ and $v_{5}$ are combined to form child vertex $c_{1}$ (4(b)),
then vertices $v_{4}$ and $v_{5}$ are combined to form child vertex $c_{1}$
(4(c)), then vertices $c_{1}$ and $v_{9}$ are combined to form child vertex
$c_{3}$ (4(d)), lastly $v_{7}$ and $v_{10}$ are combined to form child vertex
$c_{4}$ and $G_{0}$ is reached in (4(e)). Colors and hatching are used to show
coarsened vertices, with darker colors representing more parent vertices,
northeast lines representing vertices with two parent vertices, and northwest
lines representing vertices with three parent vertices

In order to rebuild the original graph $G$ from a coarsened graph $G_{0}$ we
use the following uncoarsening process. First we find the most recently
created child vertex, remove any edges in $G_{0}$ that connect to the child
vertex and then remove the child vertex from $G_{0}$. We split the child
vertex into it’s parent vertices and add both of parent vertices, along with
adding any edges that the parent vertices were a part of, into $G_{0}$. We
ensure that the parent vertices retain their original populations and both
inherit the district assignment of the child vertex. We repeat this until the
graph has $N$ nodes again, this will be the original graph $G$.

While the coarsening process can help speed up the $k$-medoids process by
working with fewer vertices, it has the trade-off of working with larger
blocks and does not allow very precise fine tuning, especially in the local
search process. To allow for more fine tuning, instead of a single value
$uc_{0}$, the user can provide an uncoarsening schedule with $q$ uncoarsening
steps

$UC:0<uc_{0}<uc_{1}<\ldots<uc_{q-1}<uc_{q}=1.$

This will first coarsen $G$ to $G_{0}$ with $\lfloor\,uc_{0}N\,\rfloor$
vertices. Then, instead of uncoarsening to $G$, uncoarsen to $G_{1}$, a graph
with $\lfloor\,uc_{1}N\,\rfloor$ vertices. At these partially coarsened steps
the local search procedures can be applied to allow for more fine tuning of
the solutions. Continue this process of uncoarsening graph $G_{i}$ to
$G_{i+1}$ for $i\in\\{0,\ldots,q-1\\}$ and doing a local search until $G_{q}$
is reached. Since $G_{q}$ is the same as $G$, the original graph will have
been reconstructed. One final local search will be run and the uncoarsening
process is complete. This process can be seen in Figure 5. Since coarsening
may not be needed in every situation, the uncoarsening schedule

$UC0:1$

can be used. This has $uc_{0}=uc_{q}=1$ and the graph will be coarsened to
$\lfloor 1\times N\rfloor=N$ vertices, which is equivalent to not coarsening.

(a) (b) (c) (d) (e) (f)

Figure 5: Uncoarsening Example with $UC:0.6,0.8,1$ on $G_{0}$ from Figure 4
(5(a)) with $q=2$ uncoarsening steps. From $G_{0}$ assign each vertex to a
district (5(b)). First uncoarsening step removes most recently created child
vertex, $c_{4}$, from graph and add the parent vertices $v_{7}$ and $v_{10}$
along with their respective edges (5(c)). Repeat until there are
$\lfloor\,uc_{1}N\,\rfloor=\lfloor\,0.8\times 10\,\rfloor=8$ vertices and
$G_{1}$ has been found (5(d)), then reassign vertices to districts through
local search (5(e)). Repeat until $G_{2}$ has been found and the graph is
uncoarsened

###  2.3 k-medoids

The final updated $k$-medoids algorithm begins with the selection of an
uncoarsening schedule ($UC$), a maximum number of $k$-medoids iterations
($\mathit{MI}$), a neighborhood function ($\mathit{NF}$), and a maximum number
of local search iterations ($\mathit{LI}$). Then $G$ is then coarsened to
$G_{0}$ and $k$ vertices are randomly selected from $G_{0}$ to be the intial
medoids. Each of these new medoids are considered a district, while the
remaining vertices are considered unassigned. To assign vertices to districts,
the algorithm first determines the district with the smallest population. Then
performs one iteration of a Breadth First Search
cormen_leiserson_rivest__1990-BFS  from each of the vertices on the border of
this minimum population district to find adjacent vertices. These adjacent
vertices are assigned to the district as they are found, provided they are not
already assigned to a district and their addition will not cause the district
population to exceed $\mathit{Pop}^{*}$. If there are no vertices that can be
added to the district, the algorithm then moves on to the next smallest
district. If no vertices can be added to the adjacent district without
exceeding $\mathit{Pop}^{*}$ then each unassigned vertex is added to the
adjacent district with the smallest population.

When all vertices are assigned to districts, the algorithm then determines the
medoid, $m$, of each district. This new medoid is selected by first using
Broder’s Algorithm to randomly cut the edges of a given district until it
forms a tree $T$. We let $P_{T,m}$ be the set of paths in $T$ that originate
at $m$. For each path $p\in P_{T,m}$ we compute the length, $C_{p}$, by
counting the number of edges along the path. Then we find the path $p^{*}$
such that $C_{p^{*}}>\sum_{p\in P_{T,m}\backslash p^{*}}C_{p}$, if it exists.
Next, we find the vertex neighboring $m$ along $p^{*}$ and make this vertex
the new medoid by reassigning $m$. We repeat until there is no path $p^{*}$,
in that case the current medoid is the new medoid. This new medoid is
sensitive to which edges are cut in Broder’s Algorithm and since the algorithm
determines the new medoid on $T$, the medoid may not appear to be visually
centered within the district. An example of the selection process for a new
medoids is diagrammed in Figure 6.

(a) (b) (c) (d) (e) (f)

Figure 6: Selecting a new medoid example with Iowa. Begin with medoid $v_{2}$:
in district:  (6(a)). Use graphical representation (6(b)) and apply Broder’s
Algorithm to randomly select edges to remove to change graph into a tree
(6(c)). Dashed lines are used to show the selected edges. Compute length,
$C_{p}$, for each path $p$ beginning at medoid $v_{2}$ (6(d)). Since
$C_{p3}>C_{p1}$ assign first vertex along path $p3$ as new medoid. Compute
$C_{p}$ for each path starting at medoid $v_{3}$ (6(e)). Since there is no
path where $C_{p^{*}}>\sum_{p\in P_{T,m}\\\ p^{*}}C_{p}$ the medoid does not
change and the process terminates. Geographical view of new medoid location
(6(f))

Once a new medoid has been found for each district, compute $\mathit{Dev}$ for
the districting plan. If $\mathit{Dev}$ is below $1$ or the algorithm has
completed $\mathit{MI}$ iterations, the $k$-medoids process is terminated.
Otherwise, all vertices are unassigned from their districts, the iteration
variable is incremented, and the district assignment process begins again.
Since this process can cause the deviation to increase between iterations, the
districting plan with the lowest deviation across all iterations is stored.
Once the algorithm has terminated, $\mathit{LI}$ iterations of local search
are performed on the districting plan with the lowest deviation. The user
defined uncoarsening schedule is then followed to alternate between
uncoarsening and performing a local search until a local search is performed
on $G_{p}$. The complete algorithm is presented in Algorithm 2. In addition to
the main goal of minimizing the population deviation, this method ensures the
generated districting plans are contiguous. In the $k$-medoids process the
Breadth First Search construction forces the initial plan to be contiguous as
TBs can only be added to adjacent districts. As each new set of medoids is
found the Breadth First Search is repeated and the districting plans remain
contiguous. In the local search process both Flip and Swap moves maintain
contiguity. A Flip will not move a TB if it is an articulation point, thus
keeping the districts contiguous. A Swap move will not move a pair
$(v_{i},v_{j})$ if either TBs are articulation points of their respective
districts or if they do not share a border with any other TB in the other
district, thus keeping the districts contiguous. Since the CMB method is just
the best between Flip and Swap moves, it also maintains contiguity. Lastly the
uncoarsening process also maintains contiguity, as when a child node is split
the parent nodes are assigned to the same district as the child node. The
final districting plan criteria is compactness, but this is not something we
directly optimize for in our method. Growing the districts from medoids in the
$k$-medoids process should keep the districts relatively compact, but we only
measure the compactness at the very end of the algorithm after the plan with
the lowest deviation has been found. This lesser focus on compactness comes
from the Indiana Citizens Redistricting Commission’s 2021 public mapping
competition where compactness was considered less important than keeping
cities and counties intact ReducedCompactness .

Algorithm 2 $k$-Medoids

1:Choose $UC$, $\mathit{MI}$, $\mathit{LI}$, and $\mathit{NF}$

2:Coarsen the graph to $G_{0}$

3:Randomly select $k$ vertices to be the initial medoids

4:while $\mathit{Dev}>1$ and Iterations $<\mathit{MI}$ do

5: while there are vertices not assigned to a district do

6: WorkingDistrict $\leftarrow$ district with minimum population

7: U $\leftarrow$ Set of all unassigned vertices adjacent to WorkingDistrict

8: Add vertices of U to WorkingDistrict provided district population does not
exceed $\mathit{Pop}^{*}$

9: if no more vertices can be added to any district then

10: Add each unassigned vertex to adjacent district with minimum population

11: end if

12: end while

13: Compute medoid of each district, these will be medoids in next iteration

14:end while

15:Perform $\mathit{LI}$ iterations of Local Search with $\mathit{NF}$ on
districting plan with smallest deviation

16:for i = 1, …, q do

17: Uncoarsen local search output to $G_{i}$

18: Perform $\mathit{LI}$ iterations of local search with $\mathit{NF}$ on
$G_{i}$

19:end for

##  3 Tests and Results

###  3.1 Tests

For testing we have chosen to use both Iowa and Florida to evaluate our
method. Iowa because it is one of the simplest nontrivial cases, having only
$99$ TBs and $k=4$ districts. We chose Florida for a more challenging case
with $9435$ VTDs and $k=27$ districts. Furthermore, Florida has been used
previously in the literature chen_rodden_2015  and can allow us to more
directly compare to the method of Chen and Rodden. To maintain existing state
boundaries we added a preprocessing step that combines all VTDs in a county,
with total county population less than $\mathit{Pop}^{*}$, to a single TB.
With this preprocessing step we reduced Florida to $4700$ TBs while keeping
$60$ of $67$ counties intact. The population data we used comes from the 2010
census by the United State Census Bureau PopulationInfoIowa ,
PopulationInfoFlorida  and the shape files we used to make the graphs also
come from the United States Census Bureau ShapeFiles . As a further note, we
will only consider a pair of TBs neighbors if they share a border with nonzero
length; TBs that only share a corner point are not considered neighbors. The
tests are all run on a single CPU (AMD 2950X at 3.5 GHz) with the algorithms
implemented in Python 3.8.5. First, we test parameters to see which allows our
model to perform the best. Next, we look at the results with the best
parameters. Lastly, we compare to the method used by Chen and Rodden.
chen_rodden_2015 .

####  3.1.1 Parameter Selection

In building the algorithm we evaluated $90$ possible parameter sets where each
set consisted of a neighborhood function, a maximum number of local search
iterations, and an uncoarsening schedule. We tested three neighborhood
functions ($\mathit{NF}\in\\{\text{Flip},\text{Swap},\text{CMB}\\})$, five
maximum local search iteration values
($\mathit{LI}\in\\{100,250,500,750,1000\\}$), and six uncoarsening schedules
($UC\in\\{UC1,UC2,UC3,UC4,UC5,UC6\\}$) defined here:

UC1: | 0.3 | 0.7 | 1 |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---
UC2: | 0.25 | 0.5 | 0.75 | 1 |  |  |  |  |
UC3: | 0.3 | 0.5 | 0.7 | 0.9 | 1 |  |  |  |
UC4: | 0.1 | 0.3 | 0.5 | 0.7 | 0.9 | 1 |  |  |
UC5: | 0.3 | 0.5 | 0.7 | 0.9 | 0.95 | 1 |  |  |
UC6: | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.8 | 0.9 | 1.

Using each of these parameter sets we ran the resulting algorithm $18$ times
on Florida, each time with a fixed number of $k$-medoids iterations
($\mathit{MI}=100$). The summaries of results for these tests can be seen in
Tables 1,2,3. The best found parameters are used in the full tests, if
multiple parameters perform well we used both.

Examining the results for the uncoarsening schedules in Table 1, we can see
that $UC6$ outperforms the other uncoarsening schedules with a minimum
deviation of $4.28\%$, whereas no other uncoarsening schedules produce
districting plans with under $10\%$ deviation. While UC6 is slower than the
other uncoarsening schedules, we feel that the lower average deviation and
minimum deviation values more than make up for the extra $60$-$120$ seconds of
runtime. Even though UC3 is modestly faster than UC5, we claim that UC5 is the
second best uncoarsening schedule because it has the next lowest mean
deviation after UC6 and approximately matches the minimum deviation of UC3.
With this in mind our full tests have been run with uncoarsening schedules UC5
and UC6.

Table 1: Uncoarsening Schedule Tests

Uncoarsening | Mean | Min | Mean | Min
---|---|---|---|---
Schedule | $\mathit{Dev}$ | $\mathit{Dev}$ | Runtime | Runtime
(UC) | ($\%$) | ($\%$) | (s) | (s)
UC1 | 97.206 | 13.686 | 173.903 | 97.630
UC2 | 92.494 | 14.207 | 192.596 | 107.234
UC3 | 91.716 | 10.318 | 200.839 | 101.113
UC4 | 93.272 | 10.153 | 258.924 | 144.469
UC5 | 86.719 | 10.897 | 225.245 | 107.370
UC6 | 85.920 | 4.285 | 298.919 | 153.151
|  |  |  |

00footnotetext: Uncoarsening Schedule: Uncoarsening Schedule being tested;
Mean $\mathit{Dev}$: mean deviation across all tests; Min $\mathit{Dev}$:
minimum deviation across all tests; Mean Runtime: mean time to coarsen to
$G_{0}$ and run $k$-medoids algorithm - including all uncoarsening and local
searches; Min Runtime: minimum time to coarsen to $G_{0}$ and run $k$-medoids
algorithm - including all uncoarsening and local searches

Examining the neighborhood function results in Table 2, we can see that there
is very little difference in the runtimes between the three neighborhood
functions. However, there is a clear difference in terms of deviation, with
CMB attaining lower deviation in the mean and minimum cases. While we expected
CMB to perform better, as it picks between the better of Flip and Swap, the
difference between Flip and Swap was surprising. With this in mind our full
tests have been run with CMB as the neighborhood function.

Table 2: Neighborhood Function Tests

Neighborhood | Mean | Min | Mean | Min
---|---|---|---|---
Function | $\mathit{Dev}$ | $\mathit{Dev}$ | Runtime | Runtime
($\mathit{NF}$) | ($\%$) | ($\%$) | (s) | (s)
Flip | 89.375 | 34.173 | 213.650 | 97.630
Swap | 135.545 | 49.255 | 226.173 | 98.584
CMB | 48.744 | 4.285 | 235.389 | 98.182
|  |  |  |

00footnotetext: Neighborhood Function: Neighborhood Function being tested;
Mean $\mathit{Dev}$: mean deviation across all tests; Min $\mathit{Dev}$:
minimum deviation across all tests; Mean Runtime: mean time to coarsen to
$G_{0}$ and run $k$-medoids algorithm - including all uncoarsening and local
searches; Min Runtime: minimum time to coarsen to $G_{0}$ and run $k$-medoids
algorithm - including all uncoarsening and local searches

Examining the results for maximum number of local search iterations in Table
3, we see that the $750$ and $1000$ iterations cases outperform the other
options. The $750$ and $1000$ iteration cases attain the lowest deviation
values of those tested and are only about $30-40$ seconds slower than the
fastest cases. While the $500$ iterations case achieves similar deviation
values to the $750$ and $1000$ iteration cases, the runtimes are much slower.
With this in mind our full tests have been run with $750$ and $1000$ maximum
local search iterations.

Table 3: Local Search Iterations Tests

Local Search | Mean | Min | Mean | Min
---|---|---|---|---
Iterations | $\mathit{Dev}$ | $\mathit{Dev}$ | Runtime | Runtime
($\mathit{LI}$) | ($\%$) | ($\%$) | (s) | (s)
100 | 101.981 | 19.446 | 173.893 | 97.630
250 | 91.352 | 15.368 | 229.778 | 116.340
500 | 89.099 | 9.995 | 308.498 | 163.629
750 | 86.525 | 4.285 | 189.575 | 125.957
1000 | 85.396 | 7.446 | 216.219 | 138.463
|  |  |  |

00footnotetext: Local Search Iterations: Maximal number of local search
iterations being tested; Mean $\mathit{Dev}$: mean deviation across all tests;
Min $\mathit{Dev}$: minimum deviation across all tests; Mean Runtime: mean
time to coarsen to $G_{0}$ and run $k$-medoids algorithm - including all
uncoarsening and local searches; Min Runtime: minimum time to coarsen to
$G_{0}$ and run $k$-medoids algorithm - including all uncoarsening and local
searches

####  3.1.2 Results

We tested the algorithm on Iowa and Florida with slightly different parameter
sets. For both states we used $100$ iterations of $k$-medoids
($\mathit{MI}=100$), CMB as the neighborhood function ($\mathit{NF}=CMB$), and
two different maximum local search iteration values
($\mathit{LI}\in\\{750,1000\\}$). Since Iowa has so few TBs, the algorithm
performs well without any uncoarsening schedule so we chose to not coarsen
Iowa at all ($UC=UC0$). Since Florida has many more TBs than Iowa, it produces
a much more complicated graph. To handle the complexity of the resultant
graph, we used two different uncoarsening schedules ($UC\in\\{UC5,UC6\\}$).
For each parameter set for each state we ran the resulting algorithm $100$
times - the summary of results for these tests are available in Table 4. We
split the results into two parts, the first being the results from performing
the $k$-medoids part of the algorithm on the fully coarsened graph with no
local search, the second being the results obtained after performing all of
the uncoarsening and local searches.

Table 4: Results for $k$-medoids algorithm on Iowa and Florida

|  |  |  $k$-medoids | Local Search | Algorithm | Total
---|---|---|---|---|---|---
State | $\mathit{LI}$ | UC | Mean | Mean | Mean | Min | Mean | Runtime | Runtime
|  |  |  $\mathit{Dev}$($\%$) | Comp |  $\mathit{Dev}$($\%$) |  $\mathit{Dev}$($\%$) | Comp | (s) | (s)
IA | 750 | UC0 | 1.126 | 0.773 | 0.693 | 0.181 | 0.763 | 0.718 | 3.478
IA | 1000 | UC0 | 1.171 | 0.769 | 0.699 | 0.271 | 0.760 | 0.690 | 3.450
FL | 750 | UC5 | 152.221 | 0.707 | 44.059 | 13.915 | 0.626 | 276.853 | 469.880
FL | 750 | UC6 | 163.727 | 0.676 | 40.769 | 12.511 | 0.612 | 328.559 | 512.829
FL | 1000 | UC5 | 153.110 | 0.707 | 42.043 | 12.061 | 0.616 | 321.880 | 523.083
FL | 1000 | UC6 | 169.677 | 0.675 | 36.707 | 10.876 | 0.607 | 391.102 | 579.940
|  |  |  |  |  |  |  |  |

00footnotetext: State: state the tests were done on; $\mathit{LI}$: number of
local search iterations; UC: uncoarsening schedule followed; $k$-medoids Mean
$\mathit{Dev}$: mean deviation of the algorithm without any local search;
$k$-medoids Mean Comp: mean compactness score of the algorithm without any
local search; Local Search Mean $\mathit{Dev}$: mean deviation after
$\mathit{LI}$ iterations of local search; Local Search Mean Comp: mean
compactness after $\mathit{LI}$ iterations of local search; Local Search Min
$\mathit{Dev}$: minimum deviation across all tests with given initial
conditions; Algorithm Runtime: mean time to coarsen to $G_{0}$ and run
$k$-medoids algorithm - including all uncoarsening and local searches; Total
Runtime: Algorithm Runtime plus mean time to read in data and summarize
results

From Table 4 we can see that our method works very well for Iowa. Without
local search the $k$-medoids algorithm is able to attain a mean deviation of
about $1.1\%$, just over the legally allowed $1\%$ deviation. Including local
search decreases the mean deviation to about $0.69\%$, just below the legal
limit. The local search also slightly decreases in the mean compactness value
(only a difference of $\sim 0.01$). There is hardly any difference between the
two parameter settings for Iowa, with the $750$ iterations case finding
slightly lower deviation values, slightly higher compactness scores, and
running a bit slower than the $1000$ iterations case.

Our method is less consistent for the more complicated case of Florida. While
including local search decreases the mean deviation value by more than
$100\%$, the minimum deviation values found were still higher than $1\%$. As
with the Iowa case, we can see that adding local search decreases the mean
compactness score. Comparing the parameter settings we can see that the final
districting plans produced using UC6 tend to have lower deviations than final
plans produced using UC5. We can also see that districting plans produced with
$1000$ iterations of local search tend to have lower deviations than plans
produced with $750$ iterations of local search. Here we also note the random
nature of this algorithm, as the minimum deviation in Table 4 is $10.876\%$
for UC6 with $1000$ iterations of local search, while the minimum deviation
found in the Section 3.1.1 parameter setting tests was $4.285\%$ for UC6 with
$750$ iterations of local search.

With our method consistently finding districting plans with deviation below
$1\%$ on Iowa, we added in an extra step to find more good districting plans
over the course of a single run. Anytime the algorithm finds a districting
plan that has a deviation less than $5\%$ during the $k$-medoids phase, we
save it. Then we perform the local search step on all of these additional
districting plans as well. The results of these tests are in Table 5. With
this additional step we are able to find another nine districting plans that
have a mean deviation of about $0.6\%$ each time we run our algorithm.
Moreover, these extra districting plans come at a cost of about $2.77$ seconds
each, slightly lower than the approximately $3.46$ seconds for the first
districting plan. Overall this leads to finding about ten districting plans
below $1\%$ deviation in just under thirty seconds.

Table 5: Additional District Results

|  |  | Mean | Mean | Mean | Mean | Total
---|---|---|---|---|---|---|---
State | $\mathit{LI}$ | UC | Additional | $\mathit{Dev}$ | Comp | Runtime | Runtime
|  |  | Districts | Per ($\%$) | Per | Per (s) | (s)
IA | 750 | UC0 | 9.290 | 0.667 | 0.755 | 2.775 | 28.577
IA | 1000 | UC0 | 8.830 | 0.671 | 0.775 | 2.775 | 28.642
|  |  |  |  |  |  |

00footnotetext: State: state tests were done on; $\mathit{LI}$: number of
local search iterations; UC: uncoarsening schedule followed; Mean Additional
Districts: mean number of districting plans with deviation below $5\%$ found
during $k$-medoids process; Mean $\mathit{Dev}$ Per: mean deviation of each
additional districting plan after local search is performed; Mean Comp Per:
mean of compactness value for each additional districting plan after local
search is performed; Mean Runtime Per: mean time to perform a local search on
each additional districting plan and summarize results; Total Runtime: mean
time to read in data, run $k$-medoids algorithm, perform local searches on
additional districting plans and summarize all results

###  3.2 Comparison to Alternative Method

Next we wanted a direct comparison to a different, but similar, method. We
implemented our own version of the algorithm described by Chen and Rodden in
chen_rodden_2015 . We ran tests on this method to directly compare to the best
versions of $k$-medoids. Since coarsening allows the local search to occur
many times, directly comparing local search iterations between the two methods
would be an uneven number of searches. To make up for this difference, we
allowed Chen and Roddens’ method take $\mathit{LI}\times\lvert UC\rvert$
iterations. This way they each use the same number of total local search
iterations. We compare the methods using the same tests from Section 3.1.2.

Table 6: Comparison of $k$-medoids Algorithm with Chen and Rodden method

|  |  |  | Local Search | Algorithm | Total
---|---|---|---|---|---|---
Method | State | $\mathit{LI}$ | UC | Mean | Min | Mean | Runtime | Runtime
|  |  |  |  $\mathit{Dev}$($\%$) |  $\mathit{Dev}$($\%$) | Comp | (s) | (s)
KMED | IA | 750 | UC0 | 0.693 | 0.181 | 0.763 | 0.718 | 3.478
CR | IA | 750 | UC0 | 7.162 | 0.426 | 0.782 | 28.155 | 29.663
KMED | IA | 1000 | UC0 | 0.699 | 0.271 | 0.760 | 0.690 | 3.450
CR | IA | 1000 | UC0 | 6.967 | 0.104 | 0.784 | 34.914 | 36.388
KMED | FL | 750 | UC5 | 44.059 | 13.915 | 0.626 | 276.853 | 469.880
CR | FL | 4500 | UC5 | 143.884 | 65.195 | 0.692 | 706.139 | 728.052
KMED | FL | 750 | UC6 | 40.769 | 12.511 | 0.612 | 328.559 | 512.829
CR | FL | 7500 | UC6 | 147.769 | 41.759 | 0.697 | 1080.773 | 1102.027
KMED | FL | 1000 | UC5 | 42.043 | 12.061 | 0.616 | 321.880 | 523.083
CR | FL | 6000 | UC5 | 147.096 | 50.079 | 0.694 | 774.962 | 797.021
KMED | FL | 1000 | UC6 | 36.707 | 10.876 | 0.607 | 391.102 | 579.940
CR | FL | 10000 | UC6 | 143.654 | 51.167 | 0.699 | 1201.172 | 1223.036
|  |  |  |  |  |  |  |

00footnotetext: Method: method used $k$-medoids algorithm (KMEDS) or the Chen
and Rodden method (CR); State: state the tests were done on; $\mathit{LI}$:
number of local search iterations; UC: uncoarsening schedule followed; Local
Search Mean $\mathit{Dev}$: mean deviation after $\mathit{LI}$ iterations of
local search; Local Search Mean Comp: mean compactness after $\mathit{LI}$
iterations of local search; Local Search Min $\mathit{Dev}$: minimum deviation
across all tests with given initial conditions; Algorithm Runtime: average
time to run the algorithm (including uncoarsening and local search time if
used); Total Runtime: Algorithm Runtime plus average time to read in data and
summarize results

From Table 6 we can see that our $k$-medoids method produces districts with
lower deviation than the other method in the same number of iterations. For
example, in the Iowa tests our method has a mean deviation below $1\%$, while
the other method has a mean deviation of about $7\%$. In the Florida tests our
method has a mean deviation of about $40\%$, while the other method has a mean
deviation of about $145\%$. Since our method runs faster than the Chen and
Rodden method ($\sim 20-30$ seconds for Iowa and $\sim 300-500$ seconds for
Florida), our method is able to find more districting plans faster. With our
method also tending to have lower deviation values, our method is more likely
to find a good districting plan in less time. Furthermore, if we account for
the extra districting plans found with by local searching on any plan found
with $\mathit{Dev}<5\%$, our method is able to find about ten districting
plans with $\mathit{Dev}<1\%$ in the approximately the same time it takes the
other method to finish. However, the Chen and Rodden method does produce
districts that are slightly more compact than ours, with the mean compactness
values being $0.02$ to $0.09$ higher. This indicates to us that our method
outperforms the Chen and Rodden method.

###  3.3 Current Districting Plans

For a baseline comparison we look at what the official districting plans for
Iowa and Florida are from the $2010$ Census.

The official plan Iowa used in 2021 can be seen in Figure 7(a). It has a
maximum population deviation of $0.005\%$ and compactness of $0.78$. Iowa law
requires all districting plans to have deviation less than $1\%$ and ours are
able to achieve the legal requirement, but with deviation slightly higher than
the official Iowa districting plan. Additionally, our districting plans have
very comparable compactness values, averaging about $0.76$.

The official plan for Florida used in 2021 can be seen in Figure 7(b). It has
a maximum population deviation of $0.0007\%$ and compactness of $0.767$, which
equates to districts with populations that differed by at most $1$ person.
However, this was done using a finer grain TB than we used, breaking the state
into Census Blocks. For Florida there are $488553$ Census Blocks, allowing for
much more fine tuning at the cost of less intuitive borders. This plan has a
deviation at least $4.2847$ lower than the best plan our method has created
and has a compactness value about $0.1$ higher than the plans produced by our
method.

(a) (b)

Figure 7: Districting plans used by Iowa (7(a)) and Florida (7(b)) in 2021

##  4 Conclusions and Future Work

Our method works very well on the simple case of Iowa, and performs reasonably
on the more challenging case of Florida. In both cases it outperforms the
method used by Chen and Rodden in terms of both time and deviation achieved.
Our method is able to combine many of the smaller ideas used throughout the
years into one method. We also introduce a seemingly novel method to build
initial districts by allowing the smallest district to expand next at each
step.

Further work needs to be done working on the remaining states and with the
updated census data. Furthermore, we could add additional methods to
potentially improve our algorithm. Examples would be the addition of multiple
object moves in the local search step from jin_2017 , the implementation of
geographs for faster local search Article3-GeoGraphs ;
Article5-GeographPractice , forcing the local search to only make moves that
strictly decrease the deviation, different selection criteria for which pair
of districts to use in the local search phase.

We could also add in a specific piece to create minority-majority districts in
states where they are required. An approach for this would be to merge
neighboring TBs with high minority populations before any clustering is done.
This would force the minority communities of interest to stay intact and could
lead to the minority-majority districts.

## Declarations

### Competing Interests

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

### Funding

We wish to acknowledge the support of the National Science Foundation for
Anthony Pizzimenti through their grant NSF-1407216.

### Availability of data and materials

The census data analyzed during the current study are available from the
census bureau’s website. The Florida data is at
https://data.census.gov/cedsci/table?text=P1&g=0400000US12%24700000
&y=2010&tid=DECENNIALPL2010.P1 and the Iowa data is at
https://data.census.gov/cedsci/table?text=P1&g=0400000US19%24700000
&y=2010&tid=DECENNIALPL2010.P1. Alternatively, the data can be reached from
https://data.census.gov, searching for table P1, selecting the year 2010,
selecting Geography tab $\rightarrow$ Voting District $\rightarrow$ (Iowa or
Florida) $\rightarrow$ All Voting Districts (VTD).

The Shape files used during the current study are available from
https://www2.census.gov/geo/tiger/TIGER2012/VTD/ .

The code and results data that support the findings of this study are
available from the corresponding author upon request.

## References

  *   * (1) Polsby, D.D., Popper, R.D.: The third criterion: Compactness as a procedural safeguard against partisan gerrymandering. Yale Law & Policy Review 9, 301–353 (1991)
  * (2) Reock, E.C.: A note: Measuring compactness as a requirement of legislative apportionment. Midwest Journal of Political Science 1, 70–74 (1961). https://doi.org/10.2307/2109043
  * (3) Schwartszberg, J.E.: Reapportionment, gerrymanders, and the notion of compactness. Minnesota Law Review 50 (1965)
  * (4) McGlone, D.: Measuring district compactness in postgis. Azavea (2016)
  * (5) Harris, C.C.: A scientific method of districting. Behavioral Science 3, 219–225 (1964). https://doi.org/10.1002/bs.3830090303
  * (6) Chambers, C.P., Miller, A.D.: A measure of bizarreness. Quarterly Journal of Political Science 5(1), 27–44 (2010). https://doi.org/10.1561/100.00009022
  * (7) Chen, J., Rodden, J.: Cutting through the thicket: Redistricting simulations and the detection of partisan gerrymanders. Election Law Journal: Rules, Politics, and Policy 14(4), 331–345 (2015). https://doi.org/10.1089/elj.2015.0317
  * (8) Bar-Natan, A., Najt, L., Schutzman, Z.: The gerrymandering jumble: Map projections permute districts’ compactness scores. Cartography and Geographic Information Science 47(4), 321–335 (2020). https://doi.org/10.1080/15230406.2020.1737575
  * (9) DeFord, D., Duchin, M., Solomon, J.: Recombination: A family of markov chains for redistricting. Harvard Data Science Review (2021). https://doi.org/10.1162/99608f92.eb30390f. https://hdsr.mitpress.mit.edu/pub/1ds8ptxu
  * (10) Garey, M.R., Johnson, D.S.: Computers and Intractability; A Guide to the Theory of NP-Completeness. W. H. Freeman & Co., USA (1990)
  * (11) Vickrey, W.: On the prevention of gerrymandering. Political Science Quarterly 76(1), 105–110 (1961). https://doi.org/10.2307/2145973
  * (12) Rossiter, D.J., Johnston, R.J.: Program group: The identification of all possible solutions to a constituency-delimitation problem. Environment and Planning A: Economy and Space 13(2), 231–238 (1981). https://doi.org/10.1068/a130231
  * (13) Cirincione, C., Darling, T.A., Orourke, T.G.: Assessing south carolinas 1990s congressional districting. Political Geography 19(2), 189–211 (2000). https://doi.org/10.1016/s0962-6298(99)00047-5
  * (14) Thoreson, J.D., Liittschwager, J.M.: Computers in behavioral science. legislative districting by computer simulation. Behavioral Science 12(3), 237–247 (1967). https://doi.org/10.1002/bs.3830120309
  * (15) Becker, A., Solomon, J.: Redistricting Algorithms. arXiv preprint 2011.09504 cs.DS (2020)
  * (16) Weaver, J.B., Hess, S.W.: A procedure for nonpartisan districting: Development of computer techniques. The Yale Law Journal 73(2), 288–308 (1963). https://doi.org/10.2307/794769
  * (17) Guest, O., Kanayet, F.J., Love, B.C.: Gerrymandering and computational redistricting. Journal of Computational Social Science 2(2), 119–131 (2019). https://doi.org/%****␣ms.bbl␣Line␣275␣****10.1007/s42001-019-00053-9
  * (18) Cohen-Addad, V., Klein, P.N., Young, N.E.: Balanced centroidal power diagrams for redistricting. Proceedings of the 26th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems (2018). https://doi.org/10.1145/3274895.3274979
  * (19) Mehrotra, A., Johnson, E.L., Nemhauser, G.L.: An optimization based heuristic for political districting. Management Science 44(8), 1100–1114 (1998). https://doi.org/10.1287/mnsc.44.8.1100
  * (20) Nagel, S.S.: Simplified bipartisan computer redistricting. Stanford Law Review 17(5), 863–899 (1965). https://doi.org/10.2307/1226994
  * (21) Kaiser, H.F.: An objective method for establishing legislative districts. Midwest Journal of Political Science 10(2), 200–213 (1966). https://doi.org/10.2307/2109148
  * (22) Ricca, F., Simeone, B.: Local search algorithms for political districting. European Journal of Operational Research 189(3), 1409–1426 (2008). https://doi.org/10.1016/j.ejor.2006.08.065
  * (23) Bozkaya, B., Erkut, E., Laporte, G.: A tabu search heuristic and adaptive memory procedure for political districting. European Journal of Operational Research 144(1), 12–26 (2003). https://doi.org/10.1016/S0377-2217(01)00380-0
  * (24) Liu, Y., Cho, W., Wang, S.: Pear: a massively parallel evolutionary computation approach for political redistricting optimization and analysis. Swarm and Evolutionary Computation 30, 78–92 (2016). https://doi.org/10.1016/j.swevo.2016.04.004
  * (25) Fifield, B., Higgins, M., Imai, K., Tarr, A.: Automated redistricting simulation using markov chain monte carlo. Journal of Computational and Graphical Statistics 29(4), 715–728 (2020). https://doi.org/10.1080/10618600.2020.1739532
  * (26) Metric Geometry and Gerrymandering Group: Gerrychain. https://github.com/mggg/GerryChain (2013)
  * (27) Magleby, D.B., Mosesson, D.B.: A new approach for developing neutral redistricting plans. Political Analysis 26(2), 147–167 (2018). https://doi.org/10.1017/pan.2017.37
  * (28) Cormen, T.H., Leiserson, C.E., Rivest, R.L.: Introduction to Algorithms. The MIT Press, Cambridge, Massachusetts (1990)
  * (29) Indiana Citizens Redistricting Commission: Mapping Competition Rules. https://portal.indiana-mapping.org/ (2021)
  * (30) United States Census Bureau: 2010: DEC Redistricting Data (PL 94-171). https://data.census.gov/cedsci/table?text=P1&g=0400000US19%24700000
&y=2010&tid=DECENNIALPL2010.P1 (08/2018) (2010)

  * (31) United States Census Bureau: 2010: DEC Redistricting Data (PL 94-171). https://data.census.gov/cedsci/table?text=P1&g=0400000US12%24700000
&y=2010&tid=DECENNIALPL2010.P1 (01/2019) (2010)

  * (32) United States Census Bureau: TIGER/Line Shapefiles. https://www2.census.gov/geo/tiger/TIGER2012/VTD/ (07/2018) (2012)
  * (33) Jin, H.: Spatial optimization methods and system for redistricting problems. PhD thesis, University of South Carolina (2017)
  * (34) King, D.M., Jacobson, S.H., Sewell, E.C., Cho, W.K.T.: Geo-graphs: An efficient model for enforcing contiguity and hole constraints in planar graph partitioning. OPERATIONS RESEARCH 60(5), 1213–1228 (2012). https://doi.org/10.1287/opre.1120.1083
  * (35) King, D.M., Jacobson, S.H., Sewell, E.C.: The geo-graph in practice: creating united states congressional districts from census blocks. Computational Optimization and Applications 69(1), 25–49 (2018). https://doi.org/10.1007/s10589-017-9936-3