ROTOR-ROUTER AUTOMATA (Spring, 2002)
I've done some work recently that combines my interests in probability,
cellular automata, number theory, and strongly convergent games. The
following is a status report on this work.
In the one-dimensional setting, the newest results are findings of my
student Lionel Levine, whose boundedness theorem for the one-dimensional
version of the "rotor-router" model implies that these automata are able
to generate good approximations to quadratic irrationals. (See section
4.)
In the two-dimensional setting, I give both experimental and theoretical
evidence for the proposition that the asymptotic growth shape for the
growth of the occupied region is a perfect circular disk. (See section
6.)
Non-probabilists may wish to skip sections 1 and 2 and the first two
paragraphs of section 3; this material is only included for motivation.
Number theorists will probably find section 4 the most interesting
section.
CONTENTS:
1. Diffusion-limited aggregation: external and internal
2. The Diaconis-Fulton model
3. The rotor router model
4. The one-dimensional case
5. The directed case
6. A more general circularity conjecture
7. The cellular automaton version
1. Diffusion-limited aggregation: external and internal
Diffusion-limited aggregation (DLA for short) is a (probabilistic) growth
process that creates a sequence of finite subsets of the lattice Z^d,
which I'll denote by S_1, S_2, ... . Here S_1 is the singleton consisting
of the origin, and for all n>1, S_n is equal to S_{n-1} with a single
lattice-point x_n adjoined, where x_n is adjacent to some element of
S_{n-1}. To choose x_n, one imagines a particle that wanders in "from
infinity", and lets x_n be the first site that the particle visits that
is adjacent to S_{n-1}. See "Diffusion-Limited Aggregation: A Model for
Pattern Formation", http://www.aip.org/pt/vol-53/iss-11/p36.html (Physics
Today On-Line).
Internal DLA is a variant of DLA in which particles are added at the origin
rather than at infinity. Unlike external DLA, which gives rise to feathery
dendritic shapes of a conjecturally fractal nature that have so far mostly
defied rigorous analysis, internal DLA is mathematically tractable: Lawler,
Bramson, and Griffeath ("Internal diffusion limited aggregation", Ann. Prob.
14 (1992), 2117-2140) proved that the set S_n, rescaled by sqrt(n), converges
in probability to a disk.
2. The Diaconis-Fulton model
Diaconis and Fulton, in an unjustly neglected article ("A growth model, a
game, an algebra, Lagrange inversion, and characteristic classes", Commutative
algebra and algebraic geometry, II (Italian) (Turin, 1990). Rend. Sem. Mat.
Univ. Politec. Torino 49 (1991), no. 1, 95--119 (1993)), define an amusing
way of "adding" two finite subsets of Z^d. If the sets S and T are disjoint,
their sum is just their union. If the intersection of S and T consists of a
single point x, the sum of S and T is S union T together with a random point
y obtained by letting a particle execute a random walk in Z^d starting from
x until it first lands on a point y not in S union T. More generally, if S
and T overlap in a set {x_1,...,x_k}, the sum of S and T is a random set of
cardinality |S|+|T|, constructed by letting k particles do random walk
starting from the respective points x_1,x_2,...,x_k, terminating in distinct
points
y_1 not in S union T,
y_2 not in S union T union {y_1},
y_3 not in S union T union {y_1,y_2},
...
y_k not in S union T union {y_1,y_2,...,y_{k-1}}
and adjoining these k new points.
This might not appear to be a well-defined definition of addition, since one
might suppose that the probability law governing the resulting set
S union T union {y_1,...,y_k}
would depend on the choice of ordering of the points x_i in the intersection
of S and T. However, it turns out that the ordering of the x_i's does not
matter. A good way to see this is to use a "stacks" picture, in the manner
of David Wilson's work on loop-erased random walk and random spanning trees
(see J. Propp and D. Wilson, "How to get a perfectly random sample from a
generic Markov chain and generate a random spanning tree of a directed graph",
Journal of Algorithms 27 (1998), 170-217).
In the case of the Diaconis-Fulton operation, one associates a stack with
each vertex v, and the stack at v consists of an infinite sequences of
vertices, each of which is a random neighbor of v; all the elements of all
the stacks are independent of one another. Then, to do a random walk from
a vertex v, one pops the top element of the stack of v (call it v') and walks
to v', then one pops the top element of _that_ stack (call it v") and walks
to v", etc., until the stopping-condition is satisfied (that is, the particle
arrives at a vertex that is not already occupied). Relative to the randomness
embodied in the entries of the stacks, the Diaconis-Fulton operation is
deterministic, and one need only prove that the outcome does not depend on
the order in which one resolves conflicts (i.e. situations in which two
particles start out attempting to occupy the same site). Proving this claim
can be reduced to the simpler claim that (S + {x}) + {y} is the same set as
(S + {y}) + {x}.
In the Diaconis-Fulton framework, internal DLA can be described very
succinctly as the result of adding the singleton {0} to itself repeatedly.
3. The rotor-router model
Let us focus on the case of the graph Z^2 (with the usual way of viewing the
four cardinal directions as North, South, East, and West), and replace the
random stacks by deterministic stacks, so that the stack at v consists of
the Eastward neighbor of v,
the Northward neighbor of v,
the Westward neighbor of v,
the Southward neighbor of v,
the Eastward neighbor of v,
the Northward neighbor of v,
the Westward neighbor of v,
the Southward neighbor of v,
...
etc. This turns internal DLA into a deterministic game satisfying the "strong
convergence property": if you add a finite number of particles to the system
(at different locations), letting each particle find an empty site before the
next particle is added, the resulting state does not depend on the order in
which the particles were added. (Note that this is an example of a strongly
convergent game, of the sort studied by Kimmo Erickson; see e.g. "Strong
convergence and the polygon property of 1-player games", Discrete Math. 153
(1996), 105-122.)
Here's a description of this evolution rule that doesn't depend on your
having read anything that appears above:
Particles are added to sites in the grid. Each site can be either
empty or occupied. When a particle is added to an empty site, the site
becomes occupied. When a particle is added to a site that is already
occupied, the new particle is ejected to one of the four neighboring
sites. The first time a specific site ejects a particle, ejection
occurs to the EAST; the second time, ejection occurs to the NORTH; then
to the WEST; then to the SOUTH; then to the EAST again, and so on,
rotating ad infinitum. As long as there are only finitely many
particles on the grid, the new particle will eventually be routed to an
empty site it can call its own (this requires proof but is not hard).
Note that each site has its own mod-4 counter, telling how many
particles (mod 4) it has already ejected, thereby determining the
direction into which the next visiting particle will be routed.
We could call these mod-4 numbers "states" rather than "counters",
and loosen the restriction that their initial values must be 0.
I call this the "rotor-router" model (though these are not the same
sort of rotors as have been studied by Serge Troubetzkoy and others).
I conjecture that:
If one adds n particles at the origin (allowing the system
to equilibrate after the addition of each new particle),
the set of occupied sites, when rescaled by the square root
of n, converges to a disk as n goes to infinity.
This guess is suggested by three things:
1) The analogy with internal DLA (defined above). Recall that in the
original (stochastic) version, a particle that arrives at an occupied
site is ejected in a random direction, with each of the four possible
directions being equally likely. Here we rotate among the four possible
directions in cyclic fashion. Recall that Lawler et al. proved that
the set of occupied sites, rescaled by sqrt(n), converges (in probability)
to a disk.
2) The proven result that in one dimension, the deterministic version
of internal DLA has the same limiting behavior as the stochastic version
(see section 4).
3) The simulation I've done with n going from 1 to 100,000. The
blob one sees at the 100,000th step appears to be extremely round!
What do I mean by "very round"? Well, one way to measure whether
a set of lattice points is round is to measure its radius in a
bunch of directions from a putative "center" and see how closely
the radii agree. I did this with the blob in sixteen different
directions, and the radii were all 179 plus or minus 1. More
specifically:
Direction of radius Length of Radius
------------------- -----------------
(+1, 0) 179
(+4,+3) 179
(+1,+1) 180
(+3,+4) 179
( 0,+1) 178
(-3,+4) 178
(-1,+1) 179
(-4,+3) 178
(-1, 0) 178
(-4,-3) 178
(-1,-1) 179
(-3,-4) 178
( 0,-1) 179
(+3,-4) 179
(+1,-1) 180
(+4,-3) 179
4. The one-dimensional case
It is natural to consider an analogous system in one dimension,
where the routers flip between pointing east and pointing west.
Say we start with all rotors pointing east, with one particle at each
site in the interval [-x(0),y(0)] with -x(0) <= 0 <= y(0), and we add
one particle at a time, letting the system equilibrate before a new
particle is added; let [-x(n),y(n)] denote the set of occupied sites
after n new particles have been added to the system. Then, as Harvard
undergraduate Lionel Levine proved (as part of a much more general
result),
y(n)^2 - x(n)^2
is bounded below by a constant, while
(y(n)-1)^2 - (x(n)+1)^2
is bounded above by a constant.
Why is this good news? Because it corresponds to the continuum limit
of IDLA!
I'll explain...
Imagine we add the particles one at a time, at a rate that is much slower
than the rate at which the particles find unoccupied sites. If we let -x
denote the location of the leftmost particle, and +y denote the location
of the rightmost particle (with both x and y functions of time), then for
IDLA we have dx:dy::y:x. This is an easy consequence of the fact that if
you do Brownian motion in the interval [a,c] with absorption at a and c,
and you start from b in (a,c), then
P(absorption at a):P(absorption at c)::(c-b):(b-a).
Solving x dx = y dy we get y^2 - x^2 = constant.
So the near-constancy of (y(n))^2 - (x(n))^2 that we see in the discrete
case mirror the exact constancy we see in the continuous case.
Here's an amusing aside on the discrete case, also proved by Levine: let
A,B,C be any Pythagorean triple with A^2+B^2=C^2. If we begin with x(0)=0
and y(0)=A, then after n = B+C-A additional particles have been added, we
get x(n)=B and y(n)=C, _exactly_.
Levine's work has other interesting consequences. Suppose we change the
rules so that when a particle that we've added to the system at the origin
ends up at a hitherto unoccupied site to the left (resp. right) of the
other occupied sites. we add r (resp. s) occupied sites at the left (resp.
right) of the growing interval. (Thus the rules described above correspond
to the case r=s=1.) If we let -x and y denote the positions of the leftmost
and rightmost occupied sites, the heuristics for the continuum limit predict
that x/y will converge to sqrt(r/s). My experiments in the case r=2, s=1
strongly suggested that in fact x - y sqrt(r/s) is bounded, and Levine has
proved this. I suspect that there is a link between what the system does
and the continued fraction for sqrt(r/s), though neither Levine nor I has
explored this. (Any of you who wish to look into these matters are very
welcome to do so, provided you let us know what you find!)
5. The directed case
I also ran a 2-dimensional directed IDLA in which each site can eject
particles only in the (+1,0) and (0,+1) directions, and in which you
put down lots of particles at (0,0) and see where they go. At large
times, the set of occupied sites seems to be the intersection of the
first quadrant with an increasingly eccentric ellipse that passes
through the origin and whose major axis has slope +1.
It is not surprising that the occupied sites stay so close to the
line y=x, since directed random walk stays close to this line.
However, I see no reason why the set of occupied sites should be
elliptical. (Perhaps a look at the corresponding stochastic IDLA
model would be enlightening.)
6. A more general circularity conjecture
One way to prove the circularity conjecture stated above (for the rotor-
router model in two dimensions) would be to find some regularities for
the behavior of the model, analogous to (but more complex than) the ones
found by Levine for one dimension. However, simulations have not been
encouraging.
A different approach would be to try to prove something much more general
and robust, in which what matters is not the fact that the ejections from
each vertex v occur in the order East, North, West, South, etc., but
the fact that the ejections from each vertex v occur in each of the four
directions in nearly equal numbers.
Here is a broadened version of the circularity conjecture: Suppose we
replace the rotating routers by some more general routing mechanism,
with the property that at any stage, and for any vertex v, the number
of particles that have been ejected from v at that stage in the process
is the same in each of the four directions, to within an error of plus
or minus one (or some other fixed bound). Then I claim that the shape
of the occupied region (after the system has equilibrated), rescaled by
the square root of the number of particles, converges to a disk.
A (related) suggestion is that one could prove the original circularity
conjecture by focussing on the function
N(v,w) = # of particles ejected from vertex v to vertex w
and using the fact that N(v) = sum_w N(v,w) is approximately harmonic
away from the origin. One might develop a theory of "approximately
harmonic" Z-valued functions and use it to show that circularity holds,
at least to within an error of o(sqrt(n)). I suspect (based on empirical
evidence) that the actual deviations from circularity are much smaller,
but even a crude circularity theorem would be worth having as a first
step.
7. The cellular automaton version
One can convert the rotor-router into a more standard CA with 24 states.
With 24=6x4 states we can keep track of the occupancy level of a site
(which can be anything between 0 and 5), as well as the direction of
the rotor (EAST, NORTH, WEST, or SOUTH). At each tick of the clock,
all the sites that are occupied by more than 1 particle eject all the
excess particles (and at the same time may receive some particles from
neighboring sites). If at time t no site has more than 5 particles,
the same will be true at time t+1, since ejection of particles brings
a site's occupancy down to 1 and at most 4 new particles can be received
during that same time-step. The direction of the rotor determines which
neighbors of a site receive the ejected particles. The 24-state model
does not allow for new particles to be added to the system, but we can
modify the model to allow for particle creation. My approach is to make
the origin a special site, with the property that whenever its occupancy
is 1, four new particles get created there (which will be sent to the
site's neighbors at the next time-step). The set of occupied sites
becomes more and more disk-like as time passes. This growth can be
modeled by an ordinary nearest-neighbor-interaction cellular automaton
(using the 24 states mentioned above along with extra states for the
origin).
I've run the CA version of the rotor-router model for just under
five million steps, and the blob formed by the occupied sites is
extremely well approximated by a disk.
Specifically, after the 4,706,517th global update, when the millionth
particle has been created at the origin, the blob of occupied sites has
radius 548-plus-or-minus-1 in each of sixteen different directions.
(I remind you that the blob grows asymmetrically, so the radii in
these sixteen different directions are all independent of one another.)
Are there other CA rules that do as good a job of making circles?
The directions in which I measured the radius were N,S,E,W,NE,NW,SE,SW,
and the eight directions (+-.6,+-.8), (+-.8,+-.6). Starting with the
E-direction, and proceeding in counterclockwise order, the radii were
(+1.0, 0.0) direction: 549
(+0.8,+0.6) direction: 549
(+0.7,+0.7) direction: 549
(+0.6,+0.8) direction: 549
( 0.0,+1.0) direction: 548
(-0.6,+0.8) direction: 549
(-0.7,+0.7) direction: 549
(-0.8,+0.6) direction: 549
(-1.0, 0.0) direction: 548
(-0.8,-0.6) direction: 547
(-0.7,-0.7) direction: 548
(-0.6,-0.8) direction: 547
( 0.0,-1.0) direction: 547
(+0.6,-0.8) direction: 547
(+0.7,-0.7) direction: 548
(+0.8,-0.6) direction: 547
The deviations from circularity, with the biggest radii in the (+,+)
quadrant and the smallest radii in the (+,-) quadrant (more or less),
reproduce on a macroscopic scale the rules of the CA (where particles
get ejected from occupied sites in the E, N, W, and S directions,
cyclically, with the first ejection occurring in the E direction).
I'd be interested to know whether there are other CA models with
simple update rules and a simple initial state that are either
rigorously proved or plausibly conjectured to have a circular
asymptotic shape. (Norm Margolus has sent me some references to
articles in the physics literature that describe models of wave
propagation and models of diffusion that have this property, such
as
http://xxx.lanl.gov/abs/nlin.CG/0106038,
http://xxx.lanl.gov/abs/comp-gas/9811002, and
http://xxx.lanl.gov/abs/comp-gas/9403004,
but I haven't had a chance to look at them yet.)
Jim Propp
There's been a new outgrowth of my work on the rotor-router model,
involving the use of rotors to simulate not random aggregation but
random walk.
Define an "over-random walk" on a graph G as an infinite path on G with
the property that there exists a universal bound B such that if v is
any vertex of G and w_1 and w_2 are any two of its neighbors, then for
every N, the number of steps the walk takes from v to w_1 between time 0
and time N is within B of the number of steps the walk takes from v to w_2
between time 0 and time N.
That is to say, a walk (viewed as a sequence of directed edges) is
over-random if the directed edges emanating from a given vertex are
traversed nearly equally often (where "nearly" means O(1)).
Note that true random walk on a finite graph is (with probability 1)
NOT over-random in my sense, since for random walk one expects the
quantities in question to differ by an unbounded amount as N gets large.
(So "over-random" seems like a more apt term than "quasi-random", which
I used when I posted an earlier version of this message to the domino
forum.)
In the case where G is finite, it's easy to show that over-random walk
visits each site its fair share of the time asymptotically (where "fair"
means "in accordance with the steady-state distribution for random walk).
Are there any theorems that imply that over-random walk (in the presence
of additional assumptions) must satisfy suitable laws of large numbers
when G is _infinite_?
Note that additional assumptions are clearly needed, since a walk in an
infinite graph that never visits any vertex twice counts as over-random
under my definition.
As a special case of an over-random walk, one could consider a particle
travelling in Z2 (starting from (0,0)) where the particle, upon its
first visit to a vertex v, chooses a random neighbor of v as the vertex
to visit next, and upon each subsequent visit to v, always moves to
whichever neighbor of v is "next" in cylic (say counterclockwise) order.
Has anyone studied this kind of walk?
I have some startling experimental evidence involving a closely related
random walk, much like the one described in the preceding paragraph,
but completely deterministic: on its first visit to a vertex v, the
particle moves to the neighbor of v that is closest to the origin
(with a simple rule for breaking ties). The other difference is that
when the particle arrives at (1,1), it is placed back at (0,0) before
continuing its walk. Let A(n) be the number of times that the particle
visits (1,1) before its nth visit to (0,0). Then I find that A(n)/n
tends quite rapidly towards p = pi / 8 (which is the probability that a
random walker starting at (0,0) visits (1,1) before returning to (0,0)).
Specifically, the discrepancy A(n) minus pn stays stays in the interval
[-.22, +3.39] for all n up to 10,000 !
This is much smaller than the discrepancy one would see for true random
walk, which would be on the order of the square root of n. (E.g., for
n=10,000, the root-mean-square deviation would be about 50: an order of
magnitude larger than 3.39.)
The interval [-.22, +3.39] seems even more surprisingly narrow when one
considers that for n between 1 and 5, the discrepancy A(n) minus pn ranges
between 0 and 2.963495. So the discrepancy interval doesn't grow much
during the particle's next 9,995 visits to (0,0). It seems that A(n) is
trying really hard to stay close to pn, with error on the order of log(n)
or maybe even less.
Has anyone looked at things like this before? I could try to prove some
theorems, but first I'd like to make sure I'm not just rediscovering
things people already know. (E.g., do numerical analysts know things
about approximately harmonic functions on the grid?)
Related questions could be formulated about weaker forms of over-randomness,
such as requiring discrepancies to be O(sqrt(n)) rather than O(1), or maybe
just o(n). In fact, I suspect that I would find results analogous to the
ones I'm looking for in the o(n) case if I looked in the early literature
on the foundations of probability (specifically, theorems that say that
if a sequence is "von-Mises random" then it satisfies various sorts of
laws of large numbers).
But the O(1) case seems so silly that maybe nobody has studied it, since
(e.g.) an "over-random" sequence of coin flips would have the property
that the cumulative numbers of heads and tails stay within a constant
bound of each other. This is the sort of axiomatization of randomness
that might be devised by a frequentialist at a casino when his dwindling
supply of cash and his fifth free drink lead him to revert to primitive
belief in a "law of averages" that surely must rescue him at any moment.
Any thoughts you may have will be appreciated!
Thanks,
Jim