Quadratic Dynamical Systems

and

Theoretical Aspects of Genetic Algorithms

Thesis Submitted for the Degree

“Doctor of Philosophy” by

Yuri Rabinovich

Submitted to the Senate of the Hebrew University of Jerusalem in December 1993 This work was carried out under the supervision of Professors Nathan Linial and Avi Wigderson

Acknowledgements

I am deeply grateful to my advisors, Prof. Nati Linial and Prof. Avi Wigderson, for their care and for the best of scientific guidance. They shall ever provide for me a model for a scientist and a teacher.

Many thanks to Dr. Alistair Sinclair, with whom Avi and I had a pleasure to work with at Princeton, and ICSI, Berkeley, for his many contributions to the development of ideas and proofs presented in this thesis.

Many thanks to Prof. Shafi Goldwasser for her support over the years, and for her interest in my work.

Many thanks also to my friends, Yossi Azar, Yossi Ben Asher, Aviad Cohen, Renata Elitzur, Igal Galperin, Aviad Kobrin, Eran London, Ilan Newman, Assaf Schuster, Mario Szegedy, Yuval Yonai, for all their help.

I am grateful to Profs. Hillel Furstenberg, Yossi Rinott and Benjamin Weiss for enlightening discussions and concrete effective suggestions.

Finally, the most important, I wish to thank my parents and my family for their love and care.

Abstract

The main purpose of this thesis is to develop general techniques for the study of computational aspects, primarily the convergence rate, of nonlinear dynamical systems. A class of interesting quadratic systems is indentified, then studied in the most general setting as well as in the particular cases of special interest. Albeit the inherent hardness of the analysis, strict time bounds are obtained in a number of cases.

From the point of view of Computer Science, an attempt is made to introduce nonlinear dynamical systems as a usefull tool in designing novel randomised algorithms for problems of combinatorial optimisation and, possibly, random generation. In view of the recent success of algorithms based on the finite Markov Chains, this seems to be a promising approach. However, only first few steps in this direction are accomplished here.

Another goal pursued in this thesis is a systematic employment of ideas of the Genetic Algorithms in a purely theoretical combinatorial setting. The resulting systems display interesting properties, reflecting the structure of the underlying combinatorial objects. A rigorous analysis of some relatively simple (but not trivial) systems of this kind is performed, yielding explicit time bounds on convergence and the grow rate of some parameters. To our best knowledge, no such results for Genetic Algorithms appeared in the literature so far.

Contents

1 General Introduction 1 1.1 Dynamical Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Genetic Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Symmetric quadratic systems . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 The summary of results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.5 Some definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2 Optimising a Linear Function by a Genetic Algorithm: An Analysis 9 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 QDS with the Bitwise Independent Mating Operator . . . . . . . . . . . . . 11 2.4 GA: Fitness + Bitwise Independent Mating . . . . . . . . . . . . . . . . . . 16

2.4.1 Symmetric case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.2 General case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.5 From Scheme to a Randomised Algorithm: a sketch . . . . . . . . . . . . . . 20 2.6 Remark on the Crossover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3 Analysis of the Binomial and the Uniform Systems on Integers 25 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2 Some preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3 The Binomial System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.3.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3.2 Mating Operator and the e-sequence of a distribution . . . . . . . . . 28 3.3.3 A Lower Bound on av (p(k)) . . . . . . . . . . . . . . . . . . . . . . 29 3.3.4 An Upper Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.3.5 Stationary distributions . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.4 The Uniform System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.4.1 Description of the system . . . . . . . . . . . . . . . . . . . . . . . . 33 3.4.2 Some facts related to the Mating Operator . . . . . . . . . . . . . . . 33 3.4.3 Bounds on av (p(n)) . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.4.4 Stationary distributions . . . . . . . . . . . . . . . . . . . . . . . . . 36

4 Symmetric Quadratic Dynamical Systems 37 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2 The general theory: Entropy, Invariants and Attractors (limit points) …. 39 4.3 The Boundary Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.3.2 The Difference Vector . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.3.3 More about the Universal Families . . . . . . . . . . . . . . . . . . . 46

5 Matchings in Graphs 49 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.2 Matchings in general graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.3 Matchings in trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

5.3.1 More about the Stationary Populations . . . . . . . . . . . . . . . . . 55 5.3.2 The rate of convergence . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.3.3 Remark on the Projection Property . . . . . . . . . . . . . . . . . . . 65

Chapter 1 General Introduction

1.1 Dynamical Systems

Many natural phenomena can be described by dynamical systems in Euclidean space. In such a description there is a fixed set of types, one for each dimension, and a point in the system (often called a state, a probability distribution on types or a population) specifies how many elements of each type exist at a given time instant. The system evolves under a fixed map F, which attempts to capture the underlying phenomenon. We mention below a few examples; for a detailed treatment see, e.g., [8].

In Physics, the system might describe the behaviour of gas molecules in a container. Here the types correspond to, say, velocity values, and the state p specifies how many molecules of each type there are at a certain point in time. The map F includes the Newtonian laws so as to produce from each state p a new state F (p) in the next time step, under some assumption about the spatial distribution of the molecules [14]. (One can also easily define a continuous-time analogue of this system). This view of dynamics is quite general and can in principle be used to describe a wide range of physical phenomena.

In Biology, the types may be the genotypes of some species (specified by the values of their genes). A population p is then simply the number of individuals of each type. The map F determines the population in the next generation according to a fixed set of rules that includes the genetic outcome of mating, the survival capacity of different types, random mutations etc. [2]. A completely different biological example describes the ecological coexistence of different species. Here the types are species, p specifies the number of creatures of each type, and F describes the outcome of their interaction: for example, one species may be the food supply of the other, or compete with it for the same resource [8].

The basic object of study in any such system is the trajectory of a point p in the state space, i.e., the sequence of points p, F (p), F (F (p)),… There are many questions one can ask about these objects, each of which has given rise to a major research area. For example, there are geometric and topological studies of the structure of the limit points or sets (attractors) of trajectories and of the initial points converging to them; ergodic and measure-theoretic questions on the recurrence of visits of trajectories to subsets of the state space; “chaos”-type questions on the sensitivity of the trajectory and its limit to small perturbations in the initial point; and computational questions regarding the time it takes for transient behaviour to die out and the system to reach its limiting behaviour. This last type of question is of central importance both to understanding the phenomenon captured by the dynamical system and in computer simulations of the system, which are often performed as a means of realizing complex limiting behaviour. Our main interest is in the last type of question, though when studying it one often needs to tackle the other questions as well.

The simplest dynamical systems arise from linear maps F. The ubiquitous example of a linear system is a Markov chain, viewed as acting on the simplex of probability distributions over types. (When F preserves the “size” of the population, it is often more convenient to take the state space to be just these probability distributions rather than all of Euclidean space; we shall adopt this view throughout.) It is well known that, if F is irreducible and aperiodic, there is a unique limit point (or stationary distribution) which attracts all initial points. Moreover, there is a beautiful spectral theory (see e.g. [22]) that explains how the second largest eigenvalue of the matrix representing F controls the rate of convergence of any trajectory to the limit. Although originally motivated by the study of natural phenomena, this theory has proved extremely useful in many mathematical disciplines. In particular, there has been a recent upsurge of interest in algorithmic applications of artificial Markov chains in Computer Science, spurred by the development of new tools for analysing the rate of convergence. Examples include probabilistic algorithms for approximate counting [7, 11], pseudo-random number generators [1] and “simulated annealing” type algorithms for combinatorial optimisation [10, 11, 21].

For nonlinear systems, in stark contrast, not only there are no global convergence rate bounds available, but even the question of whether and where the system converges can be extremely complicated (as in chaotic systems). Accordingly, such systems appeared very rarely in the Theoretical Computer Science so far.

In this thesis we introduce a special class of nonlinear systems that we call symmetric quadratic systems. Arguably, they capture several natural and artificial phenomena for which obtaining convergence rate information is computationally important. On the other hand, they seem to have a rich theory, and one may hope to obtain a quantative description of their behaviour. Study of such systems is the core of this work.

1.2 Genetic Algorithms

Genetic algorithms provide a heuristic paradigm for solving combinatorial optimisation problems based on a biological analogy. They were introduced by Holland in [9], and ever since are continuously gaining popularity as an appealing and efficient tool for finding satisfactory solutions for hard practical optimisation problems. In what follows, we slightly deviate from the standard definition, replacing the binary strings by arbitrary structures, and the crossover operator (cut and glue in a cross-over manner) by a mating operator peculiar and tailor-made for these structures.

Suppose one is searching for an optimal (or near-optimal) solution amongst a set of feasible solutions; as a concrete example, consider the problem of finding a matching of maximum cardinality amongst all matchings in a graph. The feasible solutions play the role of “creatures” that “evolve” in a space-homogeneous world. Each feasible solution defines a type. The algorithm starts with an initial distribution (or, using a more lively term, population) over simple types that is easily constructed; in our example, this might be the uniform distribution over all single-edged matchings.

Then the quality of a solution is defined in such way that “high quality” corresponds to “nearly optimal”. Thus, the quality is precisely the parameter we seek to optimise. For example, it is natural to define a quality of a matching as its size.

The population evolves in discrete steps under the application of two operators, called the mating operator (denoted p × p, or M(p)) and the fitness operator (denoted W (p)). Usually a third operator, known as mutation, is added, but for the sake of simplicity we shall ignore it here.

The mating operator is designed to increase variety in the population. It acts by forming random (according to the current distribution) pairs of solutions, and replacing them by their random (according to definition of the mating operator) children. Vaguely, all it does is creating new legal combinations of the parents’ genes. In all the cases we discuss, the mating operator preserves the expected quality of the population.

The fitness operator W is the one responsible for making the average quality increase. It acts by discarding or reducing the numbers of the low quality solutions.

In the overall genetic algorithm, the population quality constantly increases (due to W), while the possibility of appearance of an optimal solution in the population is, hopefully, never lost (due to the productive power of M).

Despite the apparent success of this approach in many experimental applications (see e.g. the proceedings of the annual International Conference on Genetic Algorithms, [17], [12] and many many others), few rigorous theoretical results exist to explain and motivate the experimental results or guide the design of such algorithms. This situation is, perhaps, not so surprising for complex systems, however there is no reason why it should be so for the simpler ones. In particular, though it is a common knowledge that Genetic Algorithms perform well in optimising a linear function on the set of binary strings (as it has short and independent schemata; see [9], [19]), even for this highly artificial and apparently easy problem no rigorous results existed. In Chapter 2 we close this gap and provide a theoretical foundation for this good performance (extending and generalizing [16]).

It is our feeling that when the underlying combinatorial structures admit a “natural” definition of the mating operator and thus the whole scheme of a Genetic Algorithm (creating a “combinatorial evolution”), the properties of the resulting system may teach us something new and important about the structures themselves. Perhaps, a good example is β^{match}, the mating operator for matchings in a graph, discussed below. Not only it provides a scheme for a completely new randomised algorithm for finding a maximum matching (which, unfortunately, we are not able at the present to analyse for general graphs), it turns out to have a surprisingly rich and mathematically interesting properties related to matching polytope, certain systems of high degree equations, etc. (see Chapter 5).

Many well known structures (like cliques, matroids, paths etc.) do in fact have a natural mating operator. Though the behaviour of the resulting systems is far from being understood, some of their general properties are established in this thesis, together with basic tools needed for their analysis.

One lesson learned from the results of Chapters 2 and 3 is that the interaction between the mating and the fitness operators is much easier to understand if M is applied repeatedly sufficiently many times, (so the system is close to stationarity), before W is applied. The reason for it is clear: while the effect on the population of a single step of the mating operator is often rather hard to quantify, a sufficiently long sequence of mating operations bring the population close to a stationary state, which is well understood. This is a reasonable strategy with mating operators which, like β^{match} (the mating operator on matchings in a graph; discussed in Chapter 5), preserve or do not substantially reduce the average quality of solutions. Interestingly, we have recently learned that this idea has led to improved performance of genetic algorithms in some experiments [23]. Given that the effect of the fitness operator is usually relatively straightforward to describe, we contend that the principal challenge in analysing the performance of genetic algorithms of the above kind lies in estimating the rate of convergence to a stationary population under the mating operator. Since the mating operator is quadratic and almost always symmetric, this is precisely the central question addressed in Chapter 4.

1.3 Symmetric quadratic systems

Informally, let us illustrate such systems by their oldest example, namely the (space homogeneous) kinetic gas model of Maxwell. Gas molecules interact by means of elastic collisions that change the distribution of their velocities. Thus the types are velocity vectors in ^{3}, and when two types i, j collide they disappear and are replaced by a random pair k, l which preserves momentum and energy. (In classical mechanics the resulting pair would be uniquely determined, but statistical mechanics makes the outcome random by postulating a uniform distribution on all possible relative positions of the balls at collision.) It is clear that the map M describing the state of the system (distribution p over the velocity vectors) after a short period of time is quadratic in terms of the old state, (i.e. M(p) is a quadratic function of p), since outcomes depend only on pairwise collisions. This map is the classical Boltzmann equation.

Note that in this system the event that the pair i, j produces k, l has the same probability as its time reversal, in which k, l produces i, j. This is the property of symmetricity which is central in our study. Systems and operators having this property will be called symmetric.

Physics is not the only branch of science in which such systems arise. In Genetics, the feature of pairwise interaction (of parents) to produce the next generation is almost universal. Here symmetry is less common, but exists in simple genetic models where, for example, a child randomly picks each gene (or each interval of genes) from either the father’s or the mother’s DNA. Here it is convenient to introduce a “complementary” child who makes the opposite choices.

The importance of understanding the evolution of natural systems of this kind, and in particular the computational issue of the rate at which they approach limiting behaviour, is well known. We were surprised to discover, however, how much common symmetric quadratic operators are in genetic-like systems of combinatorial origin. When we were looking for fitness-preserving operators, usually symmetric operators were the most natural ones for the problem. Later, when we studied the properties of such systems (in [15]), it turned out that the general framework of Quadratic Dynamical Systems clarifies and makes uniform the remarkable properties of particular systems, with little or no loss for the obtained level of generality.

The study of these systems and, in particular, their long term behaviour, is of central importance for Genetic Algorithms − when the mating operator is symmetric . Furthermore, since even simple non-symmetric quadratic operators can give rise to systems with chaotic behaviour, the GA designer should prefer, when available, symmetric mating operator.

1.4 The summary of results

In Chapters 2 and 3 we obtain rigorous results on behaviour of several systems related to Genetic Algorithms. It is also shown that essentially the same results hold when the system is simulated in a bounded (small) space. The importance of these results lies not in understanding those concrete system, but in the better understanding of how Genetic Algorithms work in general. We could find no results of this kind in the literature.

Realizing the importance of symmetric mating operators (and, in fact, proposing a new scheme, where a mating is replaced by a “supermating” - sufficiently many matings, so that convergence occurs), we dedicate Chapter 4 to the study of symmetric dynamical systems. The general results about such systems include the convergence theorem (with some unresolved exclusions), the complete description of the stationary distributions, the complete description of the invariants of the system. The theorem showing that entropy always increases in such systems (a discrete generalization of Bolzmann’s H-Theorem) plays an important role in obtaining these results. We also discuss some approaches to the open problems.

In Chapter 5 we address an important special case of such systems, involving matchings in a graph. Having obtained some interesting general results, we address the rate of convergence of the system when the underlying graph is a tree. The system is shown to have the property of rapid mixing, e.g. it converges polynomially fast, even though it has exponential size. This is apparently the first rapid mixing result for a nonlinear system.

1.5 Some definitions

In this section we introduce some terms which will be used throughout the thesis.

Let N be a set of some structures, to be called types. A probability distributions p = (p_{i}) on N shall be often refered to as a population. A fitness function f on N attaches to each l ∈ N a nonnegative value, which can be viewed as the quality or fitness of l. The average or expected fitness of a distribution p on N is

av (p) = ∑_{l∈N} f(l) p_{l}

The fitness operator W acting on distributions on N (called also populations) is completely defined by the fitness function:

W (p)_{l} = \frac{f(l) p_{l}}{av (p)}

Note that it increases the probability of each string proportionally to its fitness. Thus, the proportion of the above-arevage strings in the population grows, and so does the average fitness.

The mating operator M is somewhat more complicated. For each quadruple (i, j, k, l) ∈ N^{4}, let β_{ijkl} denote the probability that, in a mating between the two “parent” types i and j, the two types k and l are produced as “offspring”. Since any mating must have a definite outcome, we always have ∑_{k,l∈N} β_{ijkl} = 1 for all i, j ∈ N. We will adopt the view that the order of parents and offspring is not significant, so that β is symmetric in i and j, and in k and l, i.e., β_{ijkl} = β_{jikl} = β_{ijlk}. If, furthermore, β_{ijkl} = β_{klij} (meaning that the mating operation is locally time reversible), β is called symmetric.

Given two populations p and q, we define × as

(p × q)_{l} = ∑_{i,j,k∈N} p_{i}q_{j}β_{ijkl}  ∀l ∈ N .  (1.1)

The mating operator M is defined as M(p) = p × p. It should be clear that the definition of M describes a mating process in which two parents are selected independently at random from the population p and mated to form random offspring according to the distribution specified by β. The operator M (or, equivalently, β) defines a (deterministic) quadratic mapping from the simplex of probability distributions over N to itself. In what follows, M will always preserve the average fitness.

Given such N , β, f and the initial population p(0), a Genetic Algorithm (GA) in each (descrete) moment of time t = 1, 2, .. produces a new population p(t) from the old one:

p(t) = W M p(t − 1)  (1.2)

Sometimes M will be replaced by a “supermating” M̃, which is obtained by applying M a specified number of times. Thus, unlike in most application, we define Genetic Algorithms in an abstract setting, and for the most part of the thesis will not be conserned with simulations having finite population size. (See, however, Section 2.5).

A Quadratic Dinamical System (QDS) requires only N , β and p(0), and is defined by

p(t) = M p(t − 1)  (1.3)

Next we discuss three important mating operators on binary strings of a fixed length n.

The first is called the bitwise independent mating operator. The simplest way to describe it is as follows: given a pair of strings {s, r} the i-th bit of the first offspring is determined by choosing randomly and independently, the i-th bit of s or the i-th bit of r. The i-th bit of the second offspring is the one rejected by the first offspring. Though this description is clearer than a direct definiton of β, we bring the later as well: β_{ijkl} is 2^{−d(i,j)} if k ∪ l = i ∪ j (viewed as multisets corresponding to strings), and is 0 otherwise. The metric used is the Hamming distance on {0, 1}^{n}.

Another interesting mating operator is the crossover. This was the mating operator originally used by Holland [9]. Given a pair of strings {s, r}, they are cut at the same random location 0 ≤ k ≤ n, then the head of s is glued together with the tail of r, and vise versa, producing two offsprings.

The third operator is called a generalized crossover, and incudes the first two as a special cases. Let U = {u_{1}, .., u_{m}} ⊆ _{2}^{n} a set of binary strings. Each u_{i} has an associated probability 0 < p_{i} < 1, and ∑_{i=1}^{m} p_{i} = 1. The mating operator is defined as follows: Given strings s and t, choose a random vector u in U according to the probability distribution defined by p_{i}-s. The two offsprings are produced by exchanging the values of s and t in positions corresponging to 1-es of u. The system is somewhat easier when span(U) = _{2}^{n}; in this case any exchange is in principle possible.

It is readilly checked that all three mating operators defined above are locally reversible, and thus symmetric.

The last notion we introduce is the variation distance. It appears frequently in the Probability Theory, and is defined as follows: For distributions p, q on a finite set S,

∥p − q∥ = \frac{1}{2} ∑_{i∈S} |p_{i} − q_{i}| = max_{A⊆S} |p(A) − q(A)|.

Sometimes the notion Δ(p, q) will be used as well.

Chapter 2 Optimising a Linear Function by a Genetic Algorithm: An Analysis

2.1 Introduction

The optimisation of a linear function on the set of all binary strings of length n is a standard example of a problem that is well solved by genetic algorithms (see e.g. the discussion on the short and nonconflicting schemata in [9], [19]). Computer simulations and rather convincing heuristic arguments show that this is indeed the case. However, no rigorous analysis of the long-term behaviour of this system (e.g. the rate of convergence, the structure of distributions stationary under the mating operator) was so far performed.

In fact, to our knowledge, no nontrivial system involving a genetic algorithm was so far provided with rigorous analysis. Since most applications of genetic algorithms involve problems of a very complicated nature where no other methods perform well, it is sometimes argued that development of good general heuristics is much more important than the quantitative analysis of (relatively) simple cases, which is not likely to extend to systems of real interest. But, important as they are, heuristics cannot replace mathematical foundations. If one really believes that Genetic Algorithms have a future as a powerful tool for solving complicated problems (and we do), their mathematics cannot be neglected. This applies even more to this thesis, one of the main themes of which is that Genetic Algorithms can also be used for classic, much theoretically studied problems of combinatorial optimisation. And how can one hope to understand a complex problem, if not by starting with simple instances?

Linear fitness functions appear whenever the quality of solution depends on a number of factors, such that the contribution of each factor is independent from the others. Such fitness functions, although at the first glance not too powerful, were successfully used in many complicated applications, e.g. in the famous classical program by Samuel for the game of checkers [20].

Linear functions are a good starting point for a number of reasons. First, many fundamental problems of combinatorial optimisation can be seen as optimisation of a linear function (often just the number of 1’s) on some subset of binary strings of a given length (for example, those corresponding to independent sets of vertices of a graph). Before making any attempt to understand systems of this sort, it seems to be unavoidable to realize first what is going on when this subset consists of all strings of a given length.

Second, such systems turn out to be conceptually (as well as mathematically) interesting. In particular, they display some fundamental properties shared by a much wider class of systems, and thus can be used for forming new conjectures and checking new proof techniques for this class of systems. Some of the results in Chapters 4,5 were in fact obtained this way.

The realization of the special features of the mating operator is, perhaps, the most important lesson of this chapter. It turns out that a consecutive application of a mating operator (whenever it is the bitwise independent operator or the crossover) results in a rapid convergence to a distribution where each bit is totally independent from the others, and has the same expectation as in p(0). On the other hand, the action of W on a bitwise independent distribution can be easily described and evaluated. The whole approach of this chapter rests on these findings.

We shall also sketch briefly how the theoretical results of this chapter can be implemented to yield a randomised Monte Carlo algorithm. This algorithm is obtained by simulating the theoretical (exact) GA in bounded (small) space. Since statistical errors will inevitably occur, this system will deviate from the original one, and a separate analysis is required.

2.2 Description

The types are binary strings of length n, i.e r = r_{1}...r_{n} ∈ {0, 1}^{n}. Sometimes we call them the n-strings. For the most part of this chapter, the bitwise independent mating operator is used. In the last section it is shown that using the usual crossover would not significantly alter the results.

The fitness of a binary string s is given by a linear function f,

f(s) = f(s_{1}, ..., s_{n}) = A + ∑_{i=1}^{n} a_{i}s_{i}

Since we want f to be always nonnegative on the domain (so that the genetic algorithm scheme applies), the coefficients are assumed to satisfy A + ∑_{a_{i}<0} a_{i} ≥ 0.

2.3 QDS with the Bitwise Independent Mating Operator

In this section we temporarily forget about W, and concentrate on × alone. It is natural to begin the study of our mating operator by looking for a convenient representation or description of its “global” action. The original definition of × is “local” in nature (i.e. gives a good picture of what happens when two strings meet, but is rather unsuitable for the entire distribution).

The following claim describes the invariants of our system; it is an easy consequence of the “local” definition:

Claim 1 p and p × p have the same bit expectations and, consequently, the same average fintess.

We start by introducing an important notion of a normal distribution.

Given a distribution p over binary strings of length n, let s denote a random string according to p, and s_{i} - its i-th bit. A distribution is called normal if {s_{i}}_{i=1}^{n} are entirely independent. We shall use N(\overline{e}) = N(e_{1}, ..., e_{n}), for a normal distribution with bit expectations E(s_{i}) = P r(s_{i} = 1) = e_{i}.

Normal distributions play in our system a role analogous to that of eigenvectors in the linear systems, and will be used for a sort of “spectral decomposition”. The following proposition establishes he main property of such distributions:

Proposition 2

N (e_{1}, ..., e_{n}) × N (d_{1}, ..., d_{n}) = N (\frac{e_{1} + d_{1}}{2} , ..., \frac{e_{n} + d_{n}}{2} )

in particular, N(e_{1}, .., e_{n}) × N(e_{1}, .., e_{n}) = N(e_{1}, .., e_{n}).

Proof: Operator × acts independently on each bit, and thus preserves the bitwise independence. Also, it is easily verified that for any pair of distributions p, q the expectation of the i-th bit in p × q is the average of those in p and q.

In what follows, we construct a many-to-one mapping from distributions on {0, 1}^{n} to distributions over the continuous n-dimensional unit cube. This mapping will prove to be very helpful, and is constructed using the normal distributions.

Observe first that for r ∈ {0, 1}^{n} the normal distribution N(r) has all its weight on the string r. Therefore it is possible to represent any distribution p on {0, 1}^{n} as a convex combination of normal distributions:

p = ∑_{r∈{0,1}^{n}} p_{r} N(r)

In general there is no uniqueness, and infinitely many such representations are possible.

Now, for an arbitrary convex combination of normal distributions (which is itself a distribution on n-strings)

q = ∑_{v∈[0,1]^{n}} w_{v}N(v)

define a random variable X on n-dimensional unit cube [0, 1]^{n} ⊂ ^{n} by

P r(X = v) = w_{v}

X is said to be associated with q.

Thus, every distribution on n-strings can be associated with a random variable on ncube (generally, in many ways), and every such random variable can be associated (uniquely!) with a distribution on n-strings. The usefulness of this association is demonstrated by the following theorem:

Theorem 3 Let p, q be distributions on binary n-strings, and let X, Y be the associated random variables. Then p × q is associated with \frac{X+Y}{2} .

Proof: Immediate from proposition 2.

We have now rather convenient way of looking at our system; we proceed to obtain yet another representation, and then address the rate of convergence. Let p be a distribution on n-strings, X an associated random variable and B ⊆ {1, 2, .., n} a set of indices. Also, let s_{i} be (as before) the i-th bit of a random (according to p) string s. Define e_{B} = e_{B}(p) as

e_{B} = P r ( ⋀_{i∈B}(s_{i} = 1)) = E (∏_{i∈B} s_{i})

By convention, e_{}is 1. It is not hard to see that e_{B} equals also to

e_{B} = E (∏_{i∈B} X_{i})

where X_{i} is the i-th coordinate of X, an associated random variable. What is the value of e_{B}(p × q)? Using the previous equality we obtain

e_{B} (p × q) = E (∏_{i∈B} (\frac{X_{i} + Y_{i}}{2} )) = 2^{−|B|} ∑_{K⊆B} e_{K}(p) e_{B−K}(q)  (2.1)

Observe that given a distribution p over the n-strings, e_{B} -s are obtained from p_{s} by a linear transformation (not depending on p), and, furthermore, this transformation is invertible (as can be shown either by ordering lexicographically the columns and rows of the transition matrix and noticing it is upper triangular, or by explicitly expressing p_{s} by e_{B} -s using the inclusion-exclusion formula). Therefore, e_{B} -s provide an alternative basis for the space of distributions over the n-strings. This basis is particularly convenient in view of a handy equation (2.1).

Before passing to the next issue, let us deduce a nice and useful corollary about the covariances:

Corollary 4 Let p be a distribution on n-strings, s_{i}, s_{j} are random variables corresponding to the i-th and the j -th bits, and let si^{′}, sj^{′} be the corresponding bits in p′ = p × p. Then

Cov (si^{′}, sj^{′}) = \frac{1}{2}Cov (s_{i}, s_{j})

Proof: By equation 2.1,

Cov (si^{′}, sj^{′}) = ei,j^{′} − ei^{′} ej^{′} = \frac{1}{2}(e_{i,j} + e_{i} e_{j}) − e_{i} e_{j} = \frac{1}{2}Cov (s_{i}, s_{j})

The rest of this section deals with convergence of the system.

We use the following notation: Δ[∗; ∗] is the variationdistance between two distributions. Given a distribution p over binary n-strings, M^{k}p is obtained by k successive applications of M. As before, M(p) = p × p. Usually, e_{i} will be used for the expectation of the i-th bit in p, and define \overline{e} = (e_{1}, e_{2}, ..., e_{n}).

The main result of this section is:

Theorem 5 Given an initial distribution p, M^{k}p converges to N(\overline{e}). Moreover, the distance between M^{k}p and N(\overline{e}) is bounded by

Δ[M^{k}p; N(\overline{e})] ≤ n\sqrt{k} 2^{− \frac{k}{2}}

Proof: One approach to proving the theorem would be as follows: first, generalize the statement of corollary 4 to the k-th moments (provided all lower moments satisfy e_{B} = ∏_{i∈B} e_{i} ). Then, inductively assuming that the i-th moments already converged up to a negligibly small error, compute the time needed for the i + 1-th moments to converge. Finally, sum up these times. This line of reasoning is described in detail in Chapter 5, for much more sophisticated system (see the section on the rate of convergence). Here we choose a different route.

Observe first that the variationdistance between two normal distributions can be bounded as follows:

Claim 6 For u = (u_{1}, .., u_{n}), v = (v_{1}, .., v_{n}) ∈ [0, 1]^{n},

∥u − v∥_{∞} ≤ Δ[N(u); N(v)] ≤ ∥u − v∥_{1}

Proof: By induction on n. We start with the second, the more important, inequality.

Δ[N(u_{1}, .., u_{n−1}, u_{n}); N(v_{1}, .., v_{n−1}, v_{n})] ≤
Δ[N(u_{1}, ., u_{n−1}, u_{n}); N(u_{1}, ., u_{n−1}, v_{n})] + Δ[N(u_{1}, ., u_{n−1}, v_{n}); N(v_{1}, ., v_{n−1}, v_{n})] =
|u_{n} − v_{n}| + Δ[N(u_{1}, .., u_{n−1}); N(v_{1}, .., v_{n−1})]

Applying inductively the same argument, the statement follows.

To prove the first inequality, assume w.l.o.g. that ∥u − v∥_{∞} = max_{i} |u_{i} − v_{i}| is achieved at i = 1. Clearly, since N(u_{1}, ., u_{n−1}, u_{n}) refines N(u_{1}, ., u_{n−1}),

Δ[N(u_{1}, ., u_{n−1}, u_{n}); N(v_{1}, ., v_{n−1}, v_{n})] ≥ Δ[N(u_{1}, ., u_{n−1}); N(v_{1}, ., v_{n−1},)]

Therefore,

Δ[N(u); N(v)] ≥ Δ[N(u_{1}); N(v_{1})] = |u_{1} − v_{1}| = ∥u − v∥_{∞}

Next, let us bound the distance between a distribution q and the normal distribution N(\overline{e}), where \overline{e} represents the bit probabilities of q, in terms of the expected l_{1} distance between the associated random variables. Note that N(\overline{e}) is associated with the random variable having all its weight on \overline{e}.

Claim 7 Let Y be a random variable on [0, 1]^{n} associated with q, \overline{e} ∈ [0, 1]^{n}. Then

Δ[q; N(\overline{e})] ≤ E(∥Y − \overline{e}∥_{1})

Proof: By claim 6,

E(∥Y − \overline{e}∥_{1}) = ∑_{v∈[0,1]^{n}} P r(Y = v) ∥v − \overline{e}∥_{1} ≥ ∑_{v∈[0,1]^{n}} P r(Y = v) Δ[N(v) ; N(\overline{e})]

Since ∑_{v∈[0,1]^{n}} P r(Y = v) N(v) = q by definition of an associated random variable, the triangle inequality for the variation norm yields

E(∥Y − \overline{e}∥_{1}) ≥ Δ[ ∑_{v∈[0,1]^{n}} P r(Y = v) N(v); N(\overline{e})] = Δ[q; N(\overline{e})]

As claimed.

Let X be the random variable associated with the initial distribution p. Let S_{2^{k}} denote the sum of 2^{k} independent random variables equidistributed with X. Then, according to theorem (3), 2^{−k} S_{2^{k}} is associated with M^{k}p. Call this random variable Y . By claim 7,

Δ[M^{k}p; N(\overline{e})] ≤ E(∥Y − \overline{e}∥_{1}) = ∑_{i=1}^{n} E(|Y_{i} − e_{i}|)

and it suffices to bound the sum ∑_{i=1}^{n} E(|Y_{i} − p_{i}|). Notice that Y_{i}, the i-th coordinate of Y , is a normalized sum of X_{i} (as Y is of X -es), is supported on interval [0, 1], and has mean e_{i}. By Hoeffding Inequality,

P r (|Y_{i} − e_{i}| > t) ≤ 2e^{−2nt^{2}}

Therefore, for any parameter t ∈ (0, 1),

E(|Y_{i} − e_{i}|) ≤ t + P r (|Y_{i} − p_{i}| > t) ≤ t + 2 e^{−2^{k+1} t^{2}}

Choosing t = 2^{−k/2} \sqrt{k} (1 − ϵ), we obtain E(|Y_{i} − p_{i}|) ≤ 2^{−k/2} \sqrt{k}. Finally,

Δ[M^{k}p; N(\overline{e})] ≤ E (∥Y − \overline{e}∥_{1}) = ∑_{i=1}^{n} E(|Y_{i} − p_{i}|) ≤ n \sqrt{k} 2^{− \frac{k}{2}}

Corollary 8 Δ[M^{k}p; N(\overline{e})] ≤ ϵ, provided k ≥ 2 log_{2} (\frac{n}{ϵ}) + (1 + o(1)) log_{2} log_{2} (\frac{n}{ϵ}) .

2.4 GA: Fitness + Bitwise Independent Mating

Having gained some understanding of how the mating operator works, let us proceed with more complex systems − those involving both the mating and the fitness operators.

The distribution p(t) (t = 1, 2, 3, ...) will obtained from p(t − 1) by first applying the mating operator until the distribution becomes “sufficiently normal”, then applying the fitness operator. As we shall shortly see, “sufficiently normal” in our case can be replaced by the weaker “sufficiently pairwise independent”.

For a distribution p, let s denote a random string (according to p). Then av (p) = E(f(s)) = E(A + ∑ a_{j}s_{j}) = A + ∑ a_{j}e_{j}. Let ei^{′} denote the expectation of the i-th bit in W (p).

ei^{′} = \frac{∑_{r,r_{i}=1} f(r) p_{r}}{∑_{r} f(r) p_{r}} = \frac{∑_{r} f(r)r_{i} p_{r}}{∑_{r} f(r) p_{r}}

where the sum is taken over all binary strings r of length n. Thus

ei^{′} = \frac{E( f(s) · s_{i})}{E( f(s))} = \frac{av (p) · e_{i} + Cov( f(s), s_{i})}{av (p)};

ei^{′} = e_{i} + \frac{a_{i}}{av (p)}(1 − e_{i})e_{i} + \frac{1}{av (p)} ∑_{j≠i} a_{j}Cov (s_{j}, s_{i})  (2.2)

For av (W (p)) we have a general expression (which holds for any fitness function):

Proposition 9

av (W (p)) = \frac{∑_{r} f(r)^{2}p_{r}}{∑_{r} f(r) p_{r}} = \frac{E( f(s)^{2})}{E( f(s))} = E( f(s)) + \frac{Var( f(s))}{E( f(s))}

which in our case implies

av (W (p)) = av (p) + \frac{1}{av (p)} (∑_{i} a_{i}^{2}Var(s_{i}) + ∑_{i≠j} a_{i}a_{j}Cov (s_{i}, s_{j}))  (2.3)

As we have seen (corollary 4), M acting on p preserves all e_{i}-s while reducing all the covariances Cov(s_{i}, s_{j}) by a factor of 2. Combining this with equation 2.3 we obtain

av (W M^{k}(p)) = av (p) + \frac{1}{av (p)} (∑_{i} a_{i}^{2}Var(s_{i}) + ∑_{i≠j} a_{i}a_{j} \frac{1}{2^{k}} Cov (s_{i}, s_{j}))  (2.4)

2.4.1 Symmetric case

We shall first analyse a special case of a (GA) system involving both mating and fitness, with f(x_{1}, .., x_{n}) = ∑ x_{i}, and with a symmetric initial distribution. The analysis in this case comes out to be both clear and accurate.

The symmetry of p(0) means that for any permutation π ∈ S_{n} and any binary n-string r, the probabilities of π(r) and of r are equal. Since the symmetry is preserved under the mating and the fitness operators, all p(t)-s will be symmetric.

Let p(t + 1) = W M p(t) (e.g. a single application of M followed by W). Let av (t) denote (for convenience) av (p(t)). Assuming p(0) does not have all its weight on the allzeroes string, all av (t) are positive.

By (2.4),

av (t + 1) = av (t) + \frac{1}{av (t)} (∑_{i} Var(s_{i}) + \frac{1}{2} ∑_{i≠j} Cov (s_{i}, s_{j})) =

av (t) + \frac{1}{2} \frac{∑_{i} Var(s_{i})}{av (t)} + \frac{1}{2} \frac{Var(∑_{i} s_{i})}{av (t)} ≥ av (t) + \frac{1}{2} \frac{∑_{i} Var(s_{i})}{av (t)} = av (t) + \frac{1}{2} (1 − \frac{av (t)}{n})

(2.5)

since Var(∑_{i} s_{i}) ≥ 0, Var(s_{i}) = (1 − e_{i})e_{i}, and e_{i} = \frac{av (p(t))}{n} for any i, t. This implies

Proposition 10 n − av (t) ≤ (n − av (0)) (1 − \frac{1}{2n})^{t}. In other words, n − av (t) ≤ ϵ provided that t ≥ 2n log \frac{n−av (0)}{ϵ}.

How sharp is the bound of proposition 10? The e_{B} -s defined in the previous section prove to be very useful in answering this question. Define (as before, but because of the symmetry, only the size |B| = k matters) e_{(k)} = E(∏_{i=1}^{k} s_{i}). We call e_{(k)} the k-th moment of p. It is readily checked that the k-th moment of W (p) (expresses by the moments of p) equals to \frac{n−k}{n} \frac{e_{(k+1)}}{e_{(1)}} + \frac{k}{n} \frac{e_{(k)}}{e_{(1)}} , while the k-th moment of M(p) is 2^{−k} ∑_{i=0}^{k} \binom{k}{i}e_{(i)}e_{(k−i)}. (Compare with the formula 2.1).

The inequality 2.5 becomes loose when the covariance Cov(s_{i}, s_{j}) is large, or even just positive. Note that Cov(s_{i}, s_{j}) = e_{(2)} − e_{(1)}^{2}. How will this quantity behave in the course of iterative applications of M and W? It turns out that if the initial distribution has a certain property, this covariance is never positive, implying a lower bound.

Apparently, one needs to consider the entire sequence E = {e_{(0)} = 1, e_{(1)}, e_{(2)}, ..., e_{(n)}} and not only the first three terms.

Definition 1 A distribution p is called balanced if sequence E is log-concave, e.g. e_{(k)}^{2} ≥ e_{(k−1)}e_{(k+1)} for all natural k.

Remark A somewhat mysterious appearance of this definition was suggested by computer simulations, showing that after sufficiently many iterations any initial distribution eventually becomes balanced. Since the support of E is always a continuous block of integers, it is also what is sometimes called a P F_{2} sequence (see [3]). The real meaning of balancedness, as well as the reason for its apparent usefulness, remains obscure to us.

The simplest example of a balanced distribution is a symmetric distribution on singletons. The following lemma explains the importance of such distributions:

Lemma 11 The property of being balanced is preserved under both M and W.

The proof is quite simple for W; it uses only the expressions for new e_{(k)}-s, the assumption that the original distribution is balanced, and elementary inequalities. For M, on the other hand, the lemma is surprisingly difficult. By a lucky coincidence, it is proven (for general P F_{2} sequences) in Walkup [25].

If p(0) is balanced, e_{(2)} ≤ e_{(1)}^{2} for all p(t)-s. Equivalently, Cov(s_{i}, s_{j}) is always negative. Therefore, inequality 2.5 has the following counterpart for balanced distributions:

av (t + 1) ≤ av (t) + \frac{∑ Var(s_{i})}{av (t)} = av (t) + (1 − \frac{av (t)}{n})  (2.6)

and thus the bound of proposition 10 is tight:

Proposition 12 For a balanced p(0), n − av (t) ≥ (n − av (0)) (1 − \frac{1}{n})^{t} . Equivalently, n − av (t) ≥ ϵ for a t = (n − 1) log \frac{n−av (0)}{ϵ}.

2.4.2 General case

We proceed to analyse the case of a general linear function (positive on [0, 1]^{n}), and an arbitrary initial distribution.

Observe first that without loss of generality all a_{i}-s can be assumed to be positive. Indeed, suppose that a_{1} is negative. f can be viewed as f(x_{1}, x_{2}, .., x_{n}) = (A + a_{i}) − a_{1}(1 − x_{1}) + ∑_{i=2}^{n} a_{i}x_{i}. Let g(y_{1}, x_{2}, .., x_{n}) = (A + a_{i}) − a_{1}y_{1} + ∑_{i=2}^{n} a_{i}x_{i}. Consider a new system obtained from the original one: start with p̃(0) obtained from p(0) by flipping the first bit; each time the fitness operator was applied, use g instead of f. Clearly, the two systems are isomorphic, and, in particular, have the same average fitness in every moment of time. But g has less negative coefficients than f! In this way all negative coefficients can be got read of, and we may assume that all a_{i}-s are positive.

Before defining formally the GA for the task (the question is how many matings should be performed in a row), let us have a closer look on what we already know.

By equation 2.4, k successive applications of M followed by a single W, change e_{i} to

ei^{′} = e_{i} + \frac{a_{i}}{av (p)}(1 − e_{i})e_{i} + 2^{−k}R_{i}  (2.7)

where

R_{i} = \frac{1}{av (p)} ∑_{j≠i} a_{j}Cov (s_{j}, s_{i}) = Cov (\frac{∑_{j≠i} a_{j}s_{j}}{av (p)} , s_{i}) ≥ −e_{i}

Despite the fact that the average fitness increases, as long as it is far from the optimum, many s_{i} can be negatively correlated, and R_{i} may have large negative value, contrary to what we want. We seek to reduce the possible negative influence of R_{i} by sufficiently many applications of M.

Assume for a moment that R_{i} is negligible. Define L = A + ∑_{i} a_{i}, the maximum of f on \{0, 1\}^{n}. Then ei^{′} = h(e_{i}), where h(x) = x + α(1 − x)x with α = \frac{a_{i}}{L}. h can be described as follows:

• For x ∈ [0, 0.25], h(x) ≥ (1 + \frac{1}{2}α)x

• For x ∈ [0.25, 0.75], h(x) ≥ x + \frac{3}{16}α

• For x ∈ [0.75, 1], (1 − h(x)) ≤ (1 − \frac{1}{2}α)(1 − x)

Therefore, if we start with e_{i} = δ and apply h repeatedly t times, where t ≥ \frac{1}{α}(2 log \frac{1}{4δ} + \frac{8}{3} + 2 log \frac{1}{4ϵ}), we are within an ϵ distance from 1. The first term estimates time needed to get from δ to 0.25, the second - from 0.25 to 0.75, the third - from 0.75 to 1 − ϵ. Note that e_{i} increases proportionally to a_{i}, and for sufficiently small ϵ the rate of convergence of av (p(t)) is governed by the smallest a_{i}!

Now let us try to perform a similar analysis without neglecting R_{i}. Assume that in the initial distribution p(0), all e_{i} ≥ δ. Define m = max_{i} \frac{L}{a_{i}}. The GA we propose is defined as follows:

Definition 2 p(t + 1) is obtained from p(t) by r = 4 + log_{2} m successive applications of M, followed by a single W.

Theorem 13 In the above GA, L − av (t) ≤ ϵ, provided that t ≥ 2m (log \frac{1}{δ} + 2 log \frac{L}{ϵ} + O(1)).

Proof: Two distinct parts of the process will be considered separately: first, from the beginning until each e_{i} is at least 0.5, second, from this point to the end.

In the first part we proceed essentially along the same lines as in the case of negligible R_{i}. Equation 2.7 yields for our scheme, for all i

e_{i}(t + 1) ≥ e_{i}(t) + \frac{1}{m}(1 − e_{i}(t)) e_{i}(t) − \frac{1}{2^{4}m}e_{i}(t)

The same type of argument as used in the investigation of h(x), shows that the number of steps needed to get from δ to 0.5 is bounded by 2m(log \frac{1}{δ} + O(1)).

An important observation is that e_{i} will strictly increase at least until e_{i} < \frac{15}{16} , and once it passed \frac{14}{16} it can never go down below this value. (Though it can oscillate in the interval [\frac{14}{16}, 1].) Therefore, in the second part of our analysis, we can assume that all e_{i} are always ≥ 0.5.

The analysis of the second part resembles that of the symmetric case. In fact, we do not need r mating operators anymore, a single mating will do (though we can only gain if there are more). As before, in equation (2.5), we have

av (t + 1) ≥ av (t) + \frac{1}{2} \frac{∑_{i} Var (a_{i}s_{i}(t))}{av (t)} = av (t) + \frac{1}{2} \frac{∑ a_{i}^{2}(1 − e_{i}(t)) e_{i}(t)}{av (t)}

By our assumptions, e_{i} ≥ \frac{1}{2} , and \frac{a_{i}e_{i}(t)}{av (t)} ≥ \frac{1}{2} \frac{a_{i}}{L} ≥ \frac{1}{2m} for all i.

av (t + 1) ≥ av (t) + \frac{1}{4m} ∑ a_{i}(1 − e_{i}(t)) = av (t) + \frac{1}{4m}(L − av (t))

The standard (by now) argument shows that 4m(log \frac{L}{ϵ} + O(1)) iterations suffice to get from \frac{L}{2} to L − ϵ.

We conclude with remark that when f has very small coefficients, and it is required only that av (t) ≥ L − ϵ for some given ϵ, the above scheme can be improved.

More concretely, if m > \frac{100n}{ϵ}, define r to be r = 4 + log \frac{100n}{ϵ} rather than 4 + log m, as before. It is a purely technical matter to check that an analogue of the above proof works in this case as well. The only difference is that the coefficients smaller than ϵ\frac{L}{n} should be ignored. The result is

Claim 14 In the scheme with r = 4 + log_{2} \frac{100n}{ϵ} (the number of M -s between two W -s), an ϵ-approximation will be obtained after O(\frac{n}{ϵ} log \frac{Ln}{δϵ}) steps.

2.5 From Scheme to a Randomised Algorithm: a sketch

In this section we sketch how to design and analise a genetic algorithm that terminates in polynomial time and finds with high probability the maximum of a linear function on {0, 1}^{n}.

So far we have discussed dynamic systems associated with distributions on binary strings, with operators defined on these distributions. Here we shall be interested in finite-size simulations of such systems, and thus we have to review our definitions.

The basis object is now a multiset containing N binary strings of length n. The new mating operator (denoted as M̃) acting on this object, is defined as follows: given a multiset of N strings, choose randomly and independently (with replacement) N items. Pair the first item with the second, the third with the fourth, etc. Each pair produces two offsprings; consider one such pair. For each position 1 ≤ i ≤ n, independently from the other, the i-th bit of the first offspring is taken from the random parent, the other offspring takes the remaining bit. The implementation of this operator requires time polynomial in N and n.

The new W̃ acts by choosing randomly and independently (with replacement) N items according to their relative frequency times the fitness. This operator can also be implemented in time polynomial in N and n (and the length of description of the coefficients of f).

Note that previously M and W were deterministic, and now they are probabilistic.

It is convenient to associate with a multiset of N binary strings a distribution on binary strings, induced by the relative frequencies of strings in the multiset. W̃ and M̃ can be viewed as probabilistic operators defined on such distributions.

Now we turn to the description of genetic algorithm: Suppose that we want to optimise a function f(x_{1}, .., x_{n}) = ∑ x_{i} (for simplicity of discussion, but in fact without loosing generality). To this end, we choose sufficiently large N (to be specified later), and start with a multiset containing singletons, each with the same multiplicity. At each step, r = 2 log_{2} n (new) mating operators are applied, followed by a single W̃. Our aim is to obtain a multiset with average fitness > n − 1, so that there is at least one maximum fitness string. The main questions we ask is how large N and t should be so that after t steps this aim is achieved with high probability.

Equivalently, we define p̃(t + 1), a random distribution associated with the multiset after t steps, as W̃ M̃^{r}p̃(t).

The main problem in analysing the present system is that unlike in linear systems, the errors (deviations from the expected behaviour) may (apriori) accumulate exponentially fast. Thus, instead of bounding the distance from the expected distribution after t steps, we perform a separate analysis, along the same lines as before, but in the presence of random noise. The following lemma establishes some stability properties of our system:

Lemma 15 Let p̃ be the distribution induced by some multiset of N strings. Comparing the expectations of the i-th bit, and the covariances of the i-th and j -th bits, in W (p̃) and W̃ (p̃), M(p̃) and M̃(p̃), it always holds (for any of these four quantities, and any i, j)

P r (|q − q̃| ≥ ϵ) ≤ 2e^{−Nϵ^{2}}

where q is the one obtained by using M, W, and q̃ - by using M̃, W̃.

The proof follows from the Hoeffding Inequality.

Therefore, since error in e_{i} accumulates linearly, and Cov(s_{i}, s_{j}) is still reduced by a factor of two, plus a small constant error term under M̃ (with respect to M), from equation 2.7 we obtain (as in the proof of theorem 13):

Proposition 16 For all i, and r = the number of mating operators applied in a row,

ẽ_{i}(t + 1) ≥ ẽ_{i}(t) + \frac{1}{n}(1 − ẽ_{i}(t)) ẽ_{i}(t) − \frac{1}{2^{r}} ẽ_{i}(t) − O(rϵ) with prob. ≥ Ω ((1 − 2e^{−Nϵ^{2}})^{r})

Recall that r = 2 log_{2} n, and therefore the coefficient of the third term is \frac{1}{n^{2}} . Now let us choose ϵ = \frac{1}{10n^{2} log^{2} n}. Analysing the behaviour of ẽ_{i}(p) as in the first part of the proof of theorem (13), we obtain

Proposition 17 Let t = 2n log n. Then, for all i,

ẽ_{i}(2n log n) ≥ 1 − \frac{1}{2n} − O(rtϵ) < 1 − \frac{1}{n}

with probability ≥ Ω ((1 − e^{−Nϵ^{2}})^{rt}).

Choosing N properly, we arrive at

Theorem 18 Applying the above scheme with N = Ω(n^{4}log^{5}n), after t = O(n log n) iterations (where each iteration involves 2 log n + 1 single operators) the average fitness > n − 1 with arbitrary high constant probability (depending on the constants in the above O-s and Ω).

2.6 Remark on the Crossover

The goal of this section is to show that the crossover mating operator, very common in the field of genetic algorithms, is also provably good for the task of optimising a linear funtion.

Crossover mating operator, like the bitwise independent operator, preserves the bit probabilities, and is symmetric. Therefore, the general theory presented in Chapter 4, implies that if the initial distribution p on the binary strings on length n has bit expectations e_{1}, ..., e_{n}, then M^{k}p converges to N(\overline{e}), as before. However, the crossover system is much less “round” than the bitwise independent one, and its analysis is harder. (See however the section on the generalized crossover in Chapter 5.)

In spite of those difficulties, optimisation of a linear function using the crossover can still be performed (providing a fast randomised algorithm) and analysed. The key point is that the only properties of M required for the previous analysis are the preservation of bit expectations, and the fact that after small number k of iterations M^{k}p becomes almost pairwise independent. But the crossover has these properties as well! We shall see it using essentially the same method as before: the n-bit system is projected on two bits.

Let p be a distribution on n-strings, s_{i}, s_{j} the i-th and j -th bits of a random string. Denote also the corresponding bits of p × p as si^{′}, sj^{′}. Then

Proposition 19

Cov (si^{′}, sj^{′}) = (1 − c) Cov (s_{i}, s_{j})

where c = c_{ij} = \frac{|i−j|}{n+1} .

Proof: Consider the entire system from the point of view of two bits: the i-th and the j -th. If the crossover does not hit between them (it happens with probability 1 − c), nothing happens. If it does, some interchange occurs:

(∗0 ∗ 0∗) × (∗1 ∗ 1∗) → {(∗1 ∗ 0∗) ; (∗0 ∗ 1∗)}
(∗1 ∗ 0∗) × (∗0 ∗ 1∗) → {(∗0 ∗ 0∗) ; (∗1 ∗ 1∗)}

Using notation p_{xy} for P r ((s_{i} = x) ∧ (s_{j} = y)), we obtain

p11^{′} = p_{11} + c (p_{01}p_{10} − p_{11}p_{00})

Expressing this in terms of e_{B} -s (defined as before) and keeping in mind the invariance of bits’ expectations, we have

eij^{′} = e_{ij} + c (e_{i} − e_{ij})(e_{j} − e_{ij}) − c e_{ij}(1 − e_{i} − e_{j} + e_{ij}) = (1 − c) e_{ij} + ce_{i}e_{j}

Finally,

Cov (si^{′}, sj^{′}) = eij^{′} − ei^{′}ej^{′} = (1 − c) (e_{ij} − e_{i}e_{j}) = (1 − c) Cov (s_{i}, s_{j}).

Since c ≥ \frac{1}{n+1} , n + 1 successive applications of the crossover operator result in an even better pairwise independence than the one caused by a single bitwise independent mating. This is precisely the sort of property we needed.

24

Chapter 3 Analysis of the Binomial and the Uniform Systems on Integers

Based on a joint work with Avi Wigderson

The two systems investigated in this chapter, though not directly related to concrete optimisation problems, belong to the class of quadratic evolutionary systems studied in this thesys. As it is usually the case with systems in this class, there is no obvious way to predict their long-term behaviour, and no standard tools for that need exist. Our aim here is to understand as well as we can their various aspects, in hope of creating a better basis for study of similar, but more involved systems.

Abandoning the traditional realm of genetic algorithms, the binary strings, let our types be the natural numbers, and let the fitness of type “i” be simply i. As always, we are interested mostly in the fitness (in this case, expectation) - preserving mating operators. Note that for such an operator the transition probabilities β_{ijkl} (probability that the parents {i, j} will have offsprings {k, l}) are 0 whenever i + j ≠ k + l. Letting the random variable defined by the first offspring of {i, j} be binomially/uniformly distributed on the interval [0, i + j], we obtain, correspondingly, the first and the second systems studied in this chapter. Note that, unlike before or afterwards, the number of types is infinite. This, however, presents no special difficulties.

As usual, p(t + 1) is defined by p(t + 1) = W M p(t). The two main objects of study in this chapter, are, as before, the behaviour of the expected fitness of p(t), av (p(t)), and the behaviour of M^{t}(p), in particular the rate of convergence and the structure of the limit stationary distribution.

In spite of the seeming similarity of definitions, it turns out that the two systems behave quite differently: av (p(t)) grows linearly for the binomial system, and exponentially for the uniform. The stationary distributions are, correspondingly, Poissonian and Geometrical.

The binomial system is the only example of a non-time reversible system discussed in this thesis. That is, the equation β_{ijkl} = β_{klij} does not hold for it. The general theory developed in the next chapter does not apply to such systems, and thus its analysis is of a special interest. On the other hand, being non-time reversible, it is a limit (in a certain sense) of such systems! Indeed, consider the bitwise independent mating operator on the binary strings of length n of the previous chapter. For a symmetric p (i.e. the probability of a string depends solely on the number of 1-es it has; see discussion on the symmetric case in the previous chapter), let p_{i} be the total weight of all strings with i 1-es. It is not hard to see that for a large n, M(p)_{i} and W (p)_{i}, after an appropriate normalization, behave essentially in the same way as in the binomial system. This analogy (though never employed explicitly) provides a possibility to transfer the methods of the previous chapter (the symmetric case) to the new setting, obtaining a meaningful and elegant proofs, as well as tight bounds. The possibility to use properties of a better understood time reversible system for the analysis of a non-reversible one is surprising, and raises a question how general this phenomenon is.

The uniform system is time reversible, and this is reflected in the structure of its stationary distributions. To our regret, not always could we find a conceptual approach for its analysis, and some of the proofs are done by brute force. The importance of this system lies in the main result − the exponential rate of growth of av (p(t)), and in the future use for designing/checking general methods for estimating the rate of convergence of M^{t}(p). In particular, it would be of interest to understand whether the method of [4] can be employed in this case.

3.2 Some preliminaries

The infinity of our systems raises a problem we did not have before: it may happen that the expected fitness of a distribution is ∞, in which case all the questions we ask are meaningless. In what follows, we assume that the initial distribution p(0) has a finite support (e.g. only a finite number of types have positive probability). Since the property of finite support is preserved under × and W, all p(k) will have a finite support.

The types are the nonnegative integers. The fitness of type ”i” is i. In this chapter we shall often use generating functions, obtained by associating the distribution p with polynomial P (x) = ∑_{i} p_{i}x^{i}. Operators on polynomials corresponding to av (p) and W (p) are given by

Claim 20 Let p be a distribution, and P (x) the associated polynomial. Define

av (P (x)) = \frac{P′(1)}{P (1)} ; W (P (x)) = \frac{xP′(x)}{P′(1)}

Then av (P (x)) = av (p), and the polynomial associated with W (p) is W (P (x)).

Observe also that for two polynomials P and Q, av (P Q) = av (P) + av (Q), and for a constant c, av (cP) = av (P).

The following nice inequality will repeatedly prove useful:

Claim 21 Let α, a_{i}, b_{i}, c_{i} be nonnegative reals (i = 1, 2, ...) and suppose that for each i,

\frac{a_{i}}{b_{i}} ≥ α\frac{b_{i}}{c_{i}}

Then

\frac{∑_{i=1}^{n} a_{i}}{∑_{i=1}^{n} b_{i}} ≥ α \frac{∑_{i=1}^{n} b_{i}}{∑_{i=1}^{n} c_{i}} .

Proof: Define x_{i} = \sqrt{\frac{a_{i}}{α}} , y_{i} = \sqrt{\frac{c_{i}}{α}}. By the Cauchy - Schwarz inequality, (∑ x_{i}y_{i})^{2} ≤ (∑ x_{i}^{2})(∑ y_{i}^{2}). Substituting the expressions for x_{i}, y_{i}, we deduce the statement of the claim.

3.3 The Binomial System

3.3.1 Description

The “types” N are all the natural numbers and zero. The transition probabilities β_{ijkl} are

β_{ijkl} = \begin{cases}\binom{i+j}{k}2^{−(i+j)} , & i + j = k + l \\ 0, & otherwise \end{cases}

(E.g. when types i and j meet, the type of the first offspring is distributed binomially on the interval [0, i + j] with p = \frac{1}{2} ). As we have already mentioned, this × operator is not time reversible. By definition, the average fitness is preserved under ×.

3.3.2 Mating Operator and the e-sequence of a distribution

The following claim shows that employment of generating functions in our case is indeed very appropriate:

Claim 22 Let p, q be distributions, and P (x), Q(x) the associated polynomials. Then the polynomial associated with p × q is (P · Q)(\frac{x+1}{2}).

Proof: Since the × operator is bilinear, it suffices to check the statement for distributions concentrated on a single point. Let p be concentrated on some type i, and q - on j. By definition of ×, in this case (p × q)_{k} = \binom{i+j}{k}2^{−(i+j)}. But this is exactly the coefficient of x^{k} in (\frac{x+1}{2})^{j}(\frac{x+1}{2})^{i}.

Definition 3 Given a distribution p, define its e-sequence (e_{0}, e_{1}, e_{2}, ...) by generating relation

P (x) = ∑_{i=0}^{∞} p_{i}x^{i} = ∑_{k=0}^{∞} \frac{e_{k}}{k!} (x − 1)^{k} .

In particular, e_{0} = 1, e_{1} = av (p). Sometimes we shall call e_{k} a (generalised) k-th moment. Observe that all e_{i}-s are nonnegative, and if e_{i} = 0, so are all e_{k}, k ≥ i (because p_{i}-s are nonnegative).

Claim 23 Let p, q be two distributions with corresponding e-sequences {c_{i}}_{i=0}^{∞} and {d_{i}}_{i=0}^{∞} respectfully. Then

• The k-th moment of p × q is

∑_{i=0}^{k} \frac{\binom{k}{i}}{2^{k}}c_{i}d_{k−i}

• The k-th moment of W (p) is

\frac{e_{k+1}}{e_{1}} + k \frac{e_{k}}{e_{1}}

In particular, e1^{′}, the first moment of W (p), equals \frac{e_{2}}{e_{1}} + 1 .

Proof: Immediate from claims (22), (20) and the definition of the e-sequence.

3.3.3 A Lower Bound on av (p(k))

Proposition 24 The generalised moments of p (denoted as e_{i},) and those of W (p) (denoted as ei^{′}) are related as follows:

\frac{ek^{′}}{ek−1^{′}} ≥ \frac{e_{k}}{e_{k−1}} for all k with e_{k} > 0.

Proof: Expanding the ek^{′} -s in terms of e_{k} -s (from Claim 23), we obtain an equivalent statement

\frac{e_{k+1} + ke_{k}}{e_{k} + (k − 1)e_{k−1}} ≥ \frac{e_{k}}{e_{k−1}} .

or

\frac{e_{k+1}}{e_{k}} + 1 ≥ \frac{e_{k}}{e_{k−1}} .

Now expressing the e_{k} -s in the terms of the p_{i}-s (by Definition 3, e_{k} = P^{(k)}(1)), we obtain for the left side

\frac{∑_{i=0}^{∞} i(i − 1)...(i − k + 1)(i − k)p_{i}}{∑_{i=0}^{∞} i(i − 1)...(i − k + 1)p_{i}} + 1 = \frac{∑_{i=0}^{∞} i(i − 1)...(i − k + 1)(i − k + 1)p_{i}}{∑_{i=0}^{∞} i(i − 1)...(i − k + 1)p_{i}}

while the right side becomes

\frac{∑_{i=0}^{∞} i(i − 1)...(i − k + 2)(i − k + 1)P_{i}}{∑_{i=0}^{∞} i(i − 1)...(i − k + 2)P_{i}} .

Now the left and the right sides can be represented (respectively) in form

\frac{∑_{i=k}^{∞} a_{i}}{∑_{i=k}^{∞} b_{i}} and \frac{∑_{i=k}^{∞} b_{i}}{∑_{i=k}^{∞} c_{i}} .

where \frac{a_{i}}{b_{i}} = \frac{b_{i}}{c_{i}} for all i ≥ k. The assumpions of Claim (21) hold, and therefore the left side is greater or equal to the right side, as claimed.

The parameter \frac{e_{2}}{e_{1}} is important, as the average fitness of the next stage distribution depends solely on it (e1^{′} = \frac{e_{2}}{e_{1}} + 1). The following lemma plays the key role in obtaining a lower bound on av (p(n)):

Lemma 25 Let p be a distribution with p_{0} < 1. Let {e_{i}}, {ei^{′}} and {e_{i}^{∗}} denote the esequences of p, W (p) and W (p) × W (p) respectively. Then

\frac{e_{2}^{∗}}{e_{1}^{∗}} ≥ \frac{e_{2}}{e_{1}} + \frac{1}{2} .

Proof:

\frac{e_{2}^{∗}}{e_{1}^{∗}}  =^{(1)}  \frac{e2^{′} + e1^{′}^{2}}{2e1^{′}} = \frac{1}{2}e1^{′} + \frac{1}{2} \frac{e2^{′}}{e1^{′}}  ≥^{(2)}
\frac{1}{2}e1^{′} + \frac{1}{2} \frac{e_{2}}{e_{1}}  =^{(3)}  \frac{1}{2} (\frac{e_{2}}{e_{1}} + 1) + \frac{1}{2} \frac{e_{2}}{e_{1}} = \frac{e_{2}}{e_{1}} + \frac{1}{2} .

(1) the first two moments of MW (p) expressed in terms of those of W (p) ;

(2) by proposition 24;

(3) e1^{′} of W (p) expressed in terms of e_{1}, e_{2} of p.

Lemma 26 Let p be some distribution with p_{0} = 0. Then

av (W M(p)) ≥ av (p) + \frac{1}{2} .

Proof: Since av (p) equals to e_{1}, it is sufficient to prove the corresponding inequality for e_{1}-s.

For every distribution p with p_{0} = 0, there exists a distribution q such that W (q) = p.

Let {e_{i}}, {ei^{′}} {e_{i}^{∗}} and {e_{i}^{+}} denote the e-sequences of q, W (q) = p, MW (q), W MW (q) =

W M(p) respectively. Then

e_{1}^{+} =^{(1)} \frac{e_{2}^{∗}}{e_{1}^{∗}} + 1 ≤^{(2)} \frac{e_{2}}{e_{1}} + \frac{1}{2} + 1 =^{(3)} e1^{′} + \frac{1}{2} .

Where

(1) e_{1}^{+} of W MW (q) expressed in terms of e_{1}^{∗}, e_{2}^{∗} of MW (q);

(2) by the previous lemma;

(3) e1^{′} of W (q) expressed in terms of e_{1}, e_{2} of q.

Now we can establish a lower bound:

Theorem 27

av (p(n + 1)) ≥ av (p(1)) + \frac{n}{2} .

Proof: The proof follows by induction from lemma (26). Indeed, since all p(k), k > 0, assign zero probability to the type “0” (as p(k) = W M(p(k − 1)) ),

av (p(n + 1)) ≥ av (p(n)) + \frac{1}{2} ≥ ... ≥ av (p(1)) + \frac{n}{2} .

3.3.4 An Upper Bound

Now we would like to show that the obtained lower bound is tight, at least for certain kinds of distributions.

Definition 4 A distribution p will be called balanced if its e-sequence is log-concave, e.g. e_{k}^{2} ≥ e_{k−1}e_{k+1} for all natural k.

We have already encountered with balanced distributions in the previous chapter. Here they are used in a very similar fashion. The simplest example of a balanced distribution is that with p_{1} = 1. The corresponding e-sequence is (1, 1, 0, 0, ..).

Balanced distribution have (as before) the following important property:

Lemma 28 If distributions p, q are balanced, so are W (p) and p × q.

The proof (and the statement) are essentially the same as in Lemma 11; the second part is proved by Walkup in [25]. Lemmas (25), (26) and theorem (27) all have reversed counterparts for balanced distributions.

Lemma 29 Let P be a balanced distribution. Using the notions of the lemma 25, we claim that

\frac{e_{2}^{∗}}{e_{1}^{∗}} ≤ \frac{e_{2}}{e_{1}} + 1 .

Proof:

\frac{e_{2}^{∗}}{e_{1}^{∗}}  =^{(1)}  \frac{e2^{′} + e1^{′}^{2}}{2e1^{′}} = \frac{1}{2}e1^{′} + \frac{1}{2} \frac{e2^{′}}{e1^{′}}  ≤^{(2)}
\frac{1}{2}e1^{′} + \frac{1}{2}e1^{′} = e1^{′}  =^{(3)}  \frac{e_{2}}{e_{1}} + 1 .

(1) the first two moments of MW (p) expressed in terms of those of W (p) ;

(2) by preservation of balancedness under W;

(3) e1^{′} of W (p) expressed in terms of e_{1}, e_{2} of p.

Lemma 30 Let p(0) be some nonzero balanced distribution. Then all p(i) and W (p(i)) are balanced, and for all natural k it holds

av (p(k + 1)) ≤ av (p(k)) + 1 .

Proof: Let q = p(k − 1) × p(k − 1). Let {e_{i}}, {ei^{′}} {e_{i}^{∗}} and {e_{i}^{+}} denote the e-sequences of q, W (q) = p(k), MW (q) and W MW (q) = p(k + 1) respectively. Then

e_{1}^{+} =^{(1)} \frac{e_{2}^{∗}}{e_{1}^{∗}} + 1 ≤^{(2)} \frac{e_{2}}{e_{1}} + 1 + 1 =^{(3)} e1^{′} + 1 .

Where

(1) e_{1}^{+} of W MW (q) expressed in terms of e_{1}^{∗}, e_{2}^{∗} of MW (q);

(2) by the preservation of balancedness and the previous lemma. (3) e1^{′} of W (q) expressed in terms of e_{1}, e_{2} of q.

Thus, the bound of theorem (27) is tight:

Theorem 31 For a balanced initial distribution p(0), av (p(n + 1)) ≤ av (p(1)) + n

3.3.5 Stationary distributions

As we have mentioned, the present system is not time reversible. Therefore the structure of its stationary distributions, and even the fact of their existence, cannot be inferred from the general theory of Chapter 4. However, keeping in mind the analogy with the symmetric systems of Chapter 2, a natural guess would be that the stationary distributions are Poissonian. As we shall see, this is indeed the case.

Definition 5 Let N(α) denote the distribution with P r(i) = \frac{α^{i}}{i!} e^{−α} for i = 0, 1, 2, ...

The generating function of N(α) is e^{α(x−1)}, and its e-sequence is (1, α, α^{2}, ...). Observe that av (N(α)) = e_{1} = α.

Theorem 32 For any α ≥ 0 there exists a unique stationary distribution with expected fitness α, and that is N(α). Moreover, if we start with an arbitrary p with e_{1} = α, and apply M repeatedly, the distribution M^{n}p converges to N(α) in the variation norm.

Proof: Using generating functions, we see that N(α) × N(β) = N(\frac{α+β}{2} ) . Thus our special distributions are stationary. Moreover, it it not hard to infer, using the e-sequences and Claim (23), that these are the only stationary distributions. But e_{1} = α is preserved under ×, and thus the only possible stationary point in the closure of {M^{n}p}_{n=0}^{∞} is N(α).

The convergence and its rate can now be obtained using the following observation: If the polynomial associated with p is P (x), then, by Claim (22), the polynomial corresponding to M^{n}(p) is

[P (1 + \frac{x − 1}{2^{n}})]^{2^{n}}

Now a tedious but straight-forward elementary calculus argument shows that the i-th derivative of the above polynomial at x = 0 (that is, p_{i}) tends (with explicit quantitative bounds on the rate) to α^{i}e^{−α} , where α = P′(1) = av (p), as claimed. To get an estimate on the rate of convergence to N(α), one needs only to consider a sufficiently long “head” of the distribution, since the “tail” tends very fast to zero.

3.4 The Uniform System

3.4.1 Description of the system

The types are (as previously) the natural numbers and zero. The transition probabilities β_{ijkl} are this time

β_{ijkl} = \begin{cases}\frac{1}{i+j+1}, & i + j = k + l \\ 0, & otherwise \end{cases}

(E.g. when types i and j meet, the type of the first offspring is uniformly distributed on the interval [0, i + j].) This system is time reversible.

The average fitness is preserved under ×.

3.4.2 Some facts related to the Mating Operator

The generating functions for the present system suit less than before, still we prefer to retain them.

Claim 33 Let P (X), Q(x) be polynomials associated with p, q. Then the polynomial associated with p × q is

\frac{∫_{1}^{x} P (t)Q(t)dt}{x − 1}

Proof: Immediate from the definition of ×.

The following two simple facts will prove useful:

Claim 34 Let F (x) be some function analytic in the neighbourhood of 1, and f(x) =

\frac{∫_{1}^{x} F (t) dt}{x−1} . Then f^{(n)}(1) = \frac{F^{(n)}(1)}{n+1} for every natural n.

Proof: It suffices to notice that for F (x) = ∑_{i=0}^{∞} a_{i}(x − 1)^{i} , f(x) has expansion f(x) = ∑_{i=0}^{∞} \frac{a_{i}}{i+1}(x − 1)^{i}.

Claim 35 Let Q(x) = ∑_{i=0}^{n} a_{i}x^{i} , P (x) = \frac{Q(x)−Q(1)}{x−1}. Then P (x) = ∑_{i=0}^{n−1} A_{i}x^{i} , with A_{k} = ∑_{i=k+1}^{n} a_{i} .

Proof:

(x − 1) ∑_{i=0}^{n−1} A_{i}x^{i} = ∑_{i=1}^{n}(A_{i−1} − A_{i})x^{i} − A_{0} = ∑_{i=1}^{n} a_{i}x^{i} − ∑_{i=1}^{n} a_{i}

= ∑_{i=0}^{n} a_{i}x^{i} − ∑_{i=0}^{n} a_{i} = Q(x) − Q(1) .

Consequences :

(1) If the coefficients of Q(x) are nonnegative, the coefficients of P (x) will also be nonnegative; the sequence {A_{i}}_{i=0}^{n−1} will be monotone nonincreasing.

(2) If a_{i} = 0 in Q(x) , then A_{i−1} = A_{i} in P (x).

3.4.3 Bounds on av (p(n))

The upper bound is easy. Suppose that p(0) is a concentrated on the first d types. Then, obviously, p(n) is concentrated on the first 2^{n}d types, implying av (p(n)) ≤ 2^{n}d = O(2^{n}).

In order to establish the lower bound, we need two lemmas, both rather technical. Let P_{n}(x) denote the polynomial associated with p(n).

Lemma 36

av (P_{n}(x)) = 1 + \frac{2}{3}av (P_{n−1}(x)) + \frac{2}{3}av (Pn−1^{′}(x)) .

Proof:

av (P_{n}(x)) = av ( W M( P_{n−1}(x))) = av (x[M(P_{n−1}(x))]′) =^{(1)} 1 + av ([M(P_{n−1}(x))]′) =

1 + \frac{[M(P_{n−1}(x))]′′|_{1}}{[M(P_{n−1}(x))]′|_{1}} =^{(2)} 1 + \frac{\frac{1}{3}[ 2P_{n−1}(1)Pn−1^{′′}(1) + 2Pn−1^{′}(1)Pn−1^{′}(1) ]}{\frac{1}{2}[2Pn−1^{′}(1)Pn−1^{′}(1) ]}

= 1 + \frac{2}{3}av (P_{n−1}(x)) + \frac{2}{3}av (Pn−1^{′}(x))

where

(1) since av (P · Q) = av (P) + av (Q)

(2) expressing the derivatives M(P_{n−1}(x)) in terms of those ofP_{n−1}(x), according to Claim (34).

In what follow we assume that the initial distribution p(0) is not concentrated on 0.

Lemma 37 av (Pn^{′}(x)) ≥ \frac{2}{3}av (P_{n}(x)) for all n ≥ 2.

Proof: First, note that the coefficients of M(P_{n}(x)) do not increase with degree, and that the coefficients of x^{0}, x^{1}, x^{2} are equal. This is a consequence of Claim (35), since

M(P_{n}(x)) = ∫_{1}^{x}[tM(Pn−1^{′}(x))]^{2}dt / (x − 1)

Now, we represent P_{n}(x) as P_{n}(x) = S(x) + T (x), where T (x) is the tail of P_{n}(x), and consists of terms of degree ≥ 5, and S(x) is the head of P_{n}(x), consisting of terms of degree ≤ 4. By claim (21), in order to show av (Pn^{′}(x)) ≥ \frac{2}{3}av (P_{n}(x)), it suffices to establish av (S′(x)) ≥ \frac{2}{3}av (S(x)) and av (T′(x)) ≥ \frac{2}{3}av (T (x)) .

Start with T. Let T (x) = ∑_{i>4} b_{i}x^{i} . Then

\frac{4}{5} av (T (x)) = \frac{4}{5} \frac{∑_{i>4} ib_{i}}{∑_{i>4} b_{i}} ≤ \frac{∑_{i>4} i(i − 1)b_{i}}{∑_{i>4} ib_{i}} = av (T′(x))

The inequality follows from Claim (21), since for all i > 4 the corresponding terms satisfy

\frac{i(i − 1)b_{i}}{ib_{i}} = i − 1 ≥ \frac{4}{5}i = \frac{4}{5} \frac{ib_{i}}{b_{i}}

The proof for S requires more work. It turns out that, up to a multiplicative factor, S(x) can be explicitly computed. Consider the chain

M(P_{n−2}(x)) →^{W} P_{n−1}(x) → P_{n−1}^{2}(x) → ∫_{1}^{x} P_{n−1}^{2}(t)dt → \frac{∫_{1}^{x} P_{n−1}^{2}(t)dt}{x − 1} = M(P_{n−1}(x)) →^{W} P_{n}(x)

What happens to the first 5 terms of the generating function along this chain? The head of M(P_{n−2}(x)) is proportional to ax^{3} + x^{2} + x + 1, with 1 ≥ a ≥ 0 (by the second consequence of Claim (35)). Using, in that order, the expressions for W (P (x)) (Claim 20), for squaring, for integrating, for dividing by x − 1 (Claim 35), and, finally, again for W (P (x)), we obtain:

ax^{3} + x^{2} + x + 1  ⟶  3ax^{3} + 2x^{2} + x  ⟶  (6a + 4)x^{4} + 4x^{3} + x^{2}  ⟶

\frac{6a+4}{5}x^{5} + x^{4} + \frac{1}{3}x^{3}, which is proportional to (18a + 12)x^{5} + 15x^{4} + 5x^{3}, ⟶

Ax^{5} + (A + 18a + 12)x^{4} + (A + 18a + 27)x^{3} + (A + 18a + 32)(x^{2} + x + 1)

for some A ≥ 0. Substituting B for A + 18a, we proceed with

(B + 12)x^{4} + (B + 27)x^{3} + (B + 32)(x^{2} + x + 1) ⟶
(4B + 48)x^{4} + (3B + 81)x^{3} + (2B + 64)x^{2} + (B + 32)x .

Thus, up to a multiplicative factor, S(1) = 10B + 225; S′(1) = 30B + 595; S′′(1) = 70B + 1090 for some B ≥ 0. Applying once more the Claim (21) we conclude that

av (S′(x)) = \frac{S′′(1)}{S′(1)} ≥ \frac{2}{3} \frac{S′(1)}{S(1)} = \frac{2}{3} av(S(x))

Now we are in the position to prove the main theorem of this section:

Theorem 38 For all n ≥ 3, av (p(n)) > \frac{10}{9}av (p(n − 1)) .

Proof: Using the generating functions, our aim is to show that av (P_{n}(x)) > \frac{10}{9}av (P_{n−1}(x)) . By lemma (36), av (P_{n}(x)) > \frac{2}{3}av (P_{n−1}(x)) + \frac{2}{3}av (Pn−1^{′}(x)). By lemma (37), av (Pn−1^{′}(x)) ≥ \frac{2}{3}av (P_{n−1}(x)). Combining the two results, we obtain

av (P_{n}(x)) > \frac{2}{3}av (P_{n−1}(x)) + \frac{4}{9}av (P_{n−1}(x) = \frac{10}{9}av (P_{n−1}(x))

Corollary 39 av (p(n + 1)) > (\frac{10}{9})^{n−3} av (p(3)) ≥ (\frac{10}{9})^{n−3} av (p(0)) .

3.4.4 Stationary distributions

Theorem 40 For any α ≥ 0 there exists a unique stationary distribution with expected fitness α, and that is the geometrical distribution with parameter \frac{α}{1+α}.

Proof: Since the present system is time reversible, it falls into the general framework of Chapter 4. Therefore a stationary distribution π must satisfy

π_{i}π_{j} = π_{k}π_{l} for i + j = k + l

and consequently is geometric, with parameter determined by the average fitness of p(0) (which remains invariant). This fact can also be verified directly: the generating function of a stationary distribution π must satisfy the integral equation

\frac{∫_{1}^{x} f^{2}(t) dt}{x − 1} = f(x)

and, therefore, the differential equation

f^{2}(x) = (x − 1)f′(x) + f(x)

with boundary conditions f(1) = 1, f′(1) = α. The solution of this equation is:

f(x) = \frac{1}{(1 + α) − αx} = ∑_{i=0}^{∞} \frac{1}{1 + α}(\frac{α}{1 + α})^{i} x^{i} .

as expected.

In fact, a stronger version of the above theorem also holds, namely, if we start with an arbitrary p with av (p) = α, and apply M repetedly, the distribution M^{n}p converges to the appropriate geometrical distribution. The question of convergence of M^{n}p is more delicate than for the binomial system. The occurrence of convergence can be inferred from the facts that the entropy of the system rises with each iteration, and that for each α there is a unique stationary distribution π(α) with average fitness α (see the next chapter for explanations). However, the rate of convergence remains an open problem, which we would like to address in the general framework of time reversible systems.

A new approach for this kind of problems (successfully applied to two classical systems and yielding new proofs) was recently developed in [4]. We believe it can be used for the present system as well, but it is left for the future study.

Chapter 4 Symmetric Quadratic Dynamical Systems

Based on a joint work with Alistair Sinclair and Avi Wigderson

This chapter is the core of the entire work, and the discussion in the general introduction leaves little to be added. Here we give a (slightly) more detailed description of the results of this chapter, along with some additional remarks. The exposition, in general lines, follows that of [15]. Section 4.3 contains some new results.

Recall that in QDS, mating (often denoted here as β) is the only operator acting on the population. In a symmetric system, for any i, j, k, l (types), the transition probabilities satisfy β_{ijkl} = β_{klij} (i.e. β is locally time-reversible). By definition of QDS, β is assumed to satisfy also β_{ijkl} = β_{jikl} = β_{ijlk}. Finally, we call β aperiodic if β_{ijij} > 0 for all i, j. Aperiodicity is an essentially technical requirement, again satisfied in most applications, that simplifies the statement of our results. As indicated in the general introduction, the operators encountered in many applications of interest are symmetric. Moreover, as we shall see presently, symmetric operators have much useful structure which we will be able to exploit.

We would like to remark that although we are focusing on systems described by pairwise interactions (as these capture the majority of natural examples and exhibit the essential structure), it should be clear that the results can be generalized in a straightforward manner to the case of m-ary interactions (in which m parents produce m offspring) for any fixed m, provided the notion of symmetry is suitably extended.

The repeated application of β defines a the sequence of points p(0), p(1), p(2),… in the simplex of probability distributions over types, p(t + 1) = p(t) × p(t) for all t ∈ . This sequence will be called the trajectory of the point p(0).

We consider here only discrete-time, finite-dimensional systems. For any symmetric quadratic system of this form, we prove that every trajectory that does not approach the boundary of the simplex converges to a fixed point; this is enough to imply global convergence of many naturally occurring systems. The conjecture is this is always the case; the last section is dedicated to discussion of this fundamental open problem. It will be refered to as the “boundary problem”. This strong property is not shared by general, asymmetric systems, where both cycling and chaotic behaviour may occur.

Next, we give a characterization of all fixed points, and describe all the linear invariants of the system. Our main analytical tool is the fact that entropy is strictly increasing along any non-trivial trajectory. This is a discrete analogue of Boltzmann’s famous H-Theorem, (see, e.g., [14]), which he proved when studying the kinetic gas model described above. Although such results, and sometimes much more, are known for certain special cases studied in Physics, our setting is more general and focuses on the combinatorial structure of the system.

The part dealing with concrete combinatorial systems, where the results of this chapter are implemented (along with developing additional theory), is left for Chapter 5.

It is our hope that some of the techniques and results presented here might prove helpful in the future study of the quadratic dynamical systems.

4.2 The general theory: Entropy, Invariants and Attractors (limit points)

As explained earlier, it is our main aim here to study the evolution of p(t) with time. The obvious first question to ask is whether p(t) tends to some limit point π as t → ∞. Clearly, since β is continuous any such population π must be a fixed point, i.e., π × π = π; we call such a population stationary. In order to investigate this and other properties of our systems, we study the entropy of the population. For a probability distribution p on any set S, define the entropy

H(p) = − ∑_{i∈S} p_{i} lg p_{i}.

Theorem 41 For a symmetric, aperiodic operator β, the entropy of any non-stationary population is strictly increased by mating; i.e., if p is not stationary then H(p × p) > H(p).

Proof: It is helpful to work with probability distributions on the set N^{2} of ordered pairs of types. For any population p on N, define the probability distribution p̃ on N^{2} by p̃_{ij} = p_{i}p_{j}. The relation between the entropies of p and p̃ is easily seen to be

H(p̃) = 2H(p).  (4.1)

We observe first that a single application of the quadratic operator β can be decomposed into three stages: creating the pair distribution p̃ from p, the applying a linear operation on p̃, followed by a recomputation of the singleton probabilities. The first stage is defined by a linear operator B as follows:

(Bp̃)_{kl} = ∑_{ij∈N^{2}} β_{ijkl}p̃_{ij} ∀ kl ∈ N^{2}.

In the second stage, we compute (p × p)_{l} for each l ∈ N by

(p × p)_{l} = ∑_{k∈N}(Bp̃)_{kl}.

It is immediate from the definition of the operator × (see equation 1.1), that these three stages do indeed correctly compute p × p. We now proceed to show that, if p is not stationary, each of the two stages causes the entropy of the pair distribution to increase; by (4.1) this will imply the theorem. I.e., we will show that

2H(p) = H(p̃) < H(Bp̃) ≤ H(p̃ × p) = 2H(p × p).

Consider the first stage. Note that B = (β_{ij,kl}) is the transition matrix of a Markov chain on N^{2}, since ∑_{k,l∈N} β_{ijkl} = 1 for all i, j ∈ N. By the assumption on the symmetry and aperiodicity of β, B is symmetric and aperiodic. Now partition the state space N^{2} into classes of “communicating states”, i.e., equivalence classes under the equivalence relation ∼ defined by ij ∼ kl iff B_{ij,kl}^{n} > 0 for some natural number n, where B^{n} is the nth power of B. Each of these classes behaves under B as an independent Markov chain, which is symmetric and aperiodic and thus has the uniform distribution as its unique stationary distribution. It is well known (see, e.g., [18, Section 5.6]) that such a Markov chain causes the entropy (or, in fact, any strictly convex function) of a non-uniform probability distribution to increase strictly on every step. This in turn implies that the entropy of the entire distribution over N^{2} increases under B, i.e., H(Bp̃) > H(p̃), unless p̃ is uniform on all the classes.

We turn now to the second stage. The effect of this stage on the pair distribution is to transform it from Bp̃ into p̃ × p. Note that this operation preserves the new marginals q_{l} = ∑_{k}(Bp̃)_{kl} but makes the components k, l independent, i.e., (p̃ × p)_{kl} = q_{k}q_{l}. It is a routine matter to check that this latter distribution has maximum entropy among all distributions on N^{2} with the given marginals q_{l}, so H(p̃ × p) ≥ H(Bp̃) as claimed.

Finally, we have seen that the increase in entropy in the first stage is strict unless p̃ is uniform on each communicating class. But in this case we have Bp̃ = p̃, which implies by the definition of the second stage that p × p = p, so p is stationary. Thus we get a strict increase in entropy unless p is stationary.

Remark: As mentioned in the introduction, Theorem 41 is a discrete version of a classical result in statistical mechanics known as the Boltzmann’s H-Theorem, which asserts that the entropy of an ideal gas increases monotonically with time (see, e.g., [14]). Here the interaction probabilities β_{ijkl} have values that describe the collisions of pairs of gas molecules. The H-Theorem is generally proved for continuous time systems and assumed to hold for discrete time without rigorous justification. Theorem 41 gives a rigorous proof in the discrete-time setting and makes clear the essential combinatorial properties required of β. The proof is different from the one in the continuous case.

Theorem 41 is a powerful analytical tool. We use it first to prove the following fundamental result, which states that every trajectory that does not approach the boundary of the simplex converges to some fixed point (stationary population). To make this statement precise, we shall employ the variation distance. We will write p(t) → π to denote the fact that ∥p(t) − π∥ → 0 as t → ∞. We say that a trajectory p(0), p(1),…, p(t),… does not approach the boundary if there exists some δ > 0 (which may depend on p(0)) such that p_{i}(t) ≥ δ for all i ∈ N and all sufficiently large t.

Theorem 42 For any symmetric, aperiodic operator β and any initial population p(0) such that the trajectory from p(0) does not approach the boundary, there exists a stationary population π such that p(t) → π.

Proof: Consider the trajectory p(0),…, p(t),…, and define the function h(t) = H(p(t)). By Theorem 41, h(t) increases monotonically with t, and since it is evidently bounded above (by lg N) it must tend to some finite limit, h say, as t → ∞. Note that this is not sufficient to ensure that p(t) itself converges. However, since the simplex is a compact set, there exists a convergent subsequence {p(t_{j}) : j = 0, 1,…} with some limit π. By our assumption that the trajectory does not approach the boundary, π must have full support (i.e., π_{i} ≥ δ > 0 for all i ∈ N). Moreover, it is easy to verify using continuity of the functions H and × that π must be stationary and that H(π) = h.

We will show that p(t) → π. To do this, we consider the information divergence, defined for distributions p, q on N by D(p||q) = ∑_{i∈N} p_{i} lg(p_{i}/q_{i}); note that D(p||q) ≥ 0, with equality iff p = q. By the well-known relationship between variation distance and information divergence (see, e.g., [5, page 58])

D(p||q) ≥ 2 ∥p − q∥^{2}

and thus it suffices to show that D(p(t)||π) tends to zero as t → ∞.

Now for any t ∈ we may write

D(p(t)||π) = ∑_{i∈N} p_{i}(t) lg p_{i}(t) − ∑_{i∈N}(lg π_{i})p_{i}(t).  (4.2)

The first term in (4.2) is just −h(t), while the second term is constant in time: this can be seen by noting that it is just the invariant function ι_{π}(p) ≡ ∑_{i∈N} (lg π_{i})p_{i}, defined in Theorem 44, evaluated at the point p(t). But since π is the limit of the subsequence p(t_{j}), by continuity the value of the invariant ι_{π}(p(t)) must be ι_{π}(π) = −H(π) = −h. Thus (4.2) yields

D(p(t)||π) = −h(t) + h → 0,

as required.

Remarks: (a) Theorem 42 is a very special property that depends crucially on the symmetry of the operator β. It is well known that asymmetric quadratic systems need not converge to a point but may exhibit more complex behaviour such as cycling or chaos [8].

(b) The above proof actually tells us rather more: it suffices to make the weaker assumption that some convergent subsequence has a limit with full support. Thus a trajectory can only fail to converge to a limit point if it approaches a connected set of equi-entropic stationary points on the boundary of the simplex.

(c) A major open question is whether the extra condition in the above theorem, that the trajectory of any initial population p(0) with full support cannot approach the boundary, is really necessary. Our conjecture is that it is not, and can in fact be proved to hold. A detailed discussion of one approach to this problem appears in the following section. However, even if it is not so, this condition is satisfied for many naturally occurring quadratic systems, and hence the trajectory converges to a limit point (with full support). This holds in particular for the quadratic system on matchings in graphs which is the subject of Chapter 5.

A further corollary of Theorem 41 is the following elegant characterization of the stationary populations.

Theorem 43 For a symmetric, aperiodic operator β, a population p is stationary iff it satisfies

p_{i}p_{j} = p_{k}p_{l}  ∀ijkl with β_{ijkl} > 0.  (4.3)

Proof: Recall from the proof of Theorem 41 that (4.3) is precisely the condition that the pair distribution p̃_{ij} = p_{i}p_{j} is a stationary distribution for the Markov chain B defined in the second stage. Therefore if it does not hold, the entropy will increase, and therefore p cannot be stationary. On the other hand, if it does hold, it implies the stationarity of p for the quadratic process, since p itself changes if and only if p̃ changes in the second stage. For if the second stage has no effect (that is, p̃ is unchanged), the third stage is performed directly after the first, and they cancel each other.

Note that the characterization of Theorem 43 is purely structural, i.e., it depends only on which of the β_{ijkl} are non-zero (more precisely, what are the equivalence classes defined by β on N^{2}) and not on the numerical values.

Theorem 43 confirms that in general there will be a continuum of stationary populations, corresponding to the set of solutions of the system of equations (4.3). Each stationary population will have a certain domain of attraction in the simplex. The question then arises of how to determine, for a given initial population satisfying the condition of Theorem 42, the particular stationary population which attracts it. The answer is to look for invariants of the operator β, i.e., functions ι on populations that satisfy ι(p × p) = ι(p) for all p.

The following family of linear invariants plays a crucial role in our systems. Recall that a population p has full support if p_{i} > 0 for all i ∈ N.

Theorem 44 Let π be any stationary population with full support. Then the function ι_{π}(p) = ∑_{l∈N} (lg π_{l}) p_{l} is an invariant.

Proof: For each l ∈ N define the function Δp_{l}(t) = p_{l}(t + 1) − p_{l}(t). We may expand Δp_{l} as

Δp_{l}(t) = ∑_{ijk}(p_{i}(t)p_{j}(t) − p_{k}(t)p_{l}(t))β_{ijkl},  (4.4)

where we have used the equation 1.1 in Section 1.5) to rewrite p_{l}(t + 1), and symmetry of β to rewrite p_{l}(t).

Now let α_{l} = lg π_{l}, and note from Theorem 43, since π is stationary, that

α_{i} + α_{j} = α_{k} + α_{l}  ∀ijkl with β_{ijkl} > 0.  (4.5)

Now using (4.4) we may write

ι_{π}(p × p) − ι_{π}(p) = ∑_{ijkl} α_{l}(p_{i}p_{j} − p_{k}p_{l})β_{ijkl} = \frac{1}{4} ∑_{ijkl}(α_{k} + α_{l} − α_{i} − α_{j})(p_{i}p_{j} − p_{k}p_{l})β_{ijkl}

where we have used the symmetry properties of β. But (4.5) implies that every term in this sum is zero, so ι_{π} is indeed invariant.

Remarks:

(a) Equation (4.4) above can be easily modified to work in continuous time: simply replace Δp_{l}(t) on the left-hand side by the derivative dp_{l}/dt. Thus, the continuous time sytems have the same invariants.

(b) Since the invariants constitute a linear space, it suffices to look at a basis of this space. It is an interesting question what can be said about the dimension of this space. Clearly, it is finite and bounded from above by |N |, but, because of dependencies as in Theorem 4.5, it can be much smaller.

(c) Though in this section we never discuss finite size implementations of dynamical systems, we still would like to bring an important property of the invariants which plays a key role in the analysis of such implementations. Suppose for a moment that p(t + 1) is obtained from p(t) by β^{match} in the presence of random noise. Then, since the systems under discussion are quadratic, the error in p(t) may (apriori) accumulate exponentially fast in t. However, it can be shown that the error in values of the invariants grows only linearly.

The significance of the above invariants is expressed in our next theorem, which states that they completely determine the stationary population with full support to which a given initial population is attracted (assuming such a limit point exists). That is, they constitute a complete systems of invariants for full support populations.

Theorem 45 Suppose p(t) → p^{∗} and q(t) → q^{∗}, where p^{∗} and q^{∗} have full support. Then ι_{π}(p(0)) = ι_{π}(q(0)) for all invariants ι_{π} as above, iff p^{∗} = q^{∗}.

Proof: We know in particular that ι_{p^{∗}}(p(0)) = ι_{p^{∗}}(q(0)), and hence by continuity of the operator × that ι_{p^{∗}}(p^{∗}) = ι_{p^{∗}}(q^{∗}). Similarly, ι_{q^{∗}}(p^{∗}) = ι_{q^{∗}}(q^{∗}). Writing out these two equalities yields

∑_{l}(lg p_{l}^{∗})p_{l}^{∗} = ∑_{l}(lg p_{l}^{∗})q_{l}^{∗}
∑_{l}(lg q_{l}^{∗})p_{l}^{∗} = ∑_{l}(lg q_{l}^{∗})q_{l}^{∗}

which can be subtracted to give

D(p^{∗}||q^{∗}) = ∑_{l}(lg p_{l}^{∗}/q_{l}^{∗})p_{l}^{∗} = ∑_{l}(lg p_{l}^{∗}/q_{l}^{∗})q_{l}^{∗} = −D(q^{∗}||p^{∗}).

By the nonnegativity of informational divergence (or just by Jensen’s inequality applied to the convex function x lg x), the left-hand side of this equation is always positive and the right-hand side always negative. Hence both sides must in fact equal zero, which happens iff p^{∗} = q^{∗}.

Remark: Theorem 45 can be proven differently, considering the set U of of all populations with some given values of invariants. Thus, U is an intersection of the simplex of probability distributions with some affine space. On one hand, since the entropy function is strictly convex, there exists a unique population in U that maximizes the entropy, and this population is necessarily stationary. On the other hand, one can show using Lagrange multipliers that any full support stationary distribution in U is a local maximum. Since the entropy function on U has no local maxima but the global one, there exists a unique stationary full support population in U.

Theorem 45 gives us a finite system of equations whose unique solution is the limit point of the trajectory of p(0). The equations consist of (4.3) together with a finite basis for the set of linear equations ι_{π}(p) = c_{π}, where the values c_{π} are determined by p(0), i.e., c_{π} = ι_{π}(p(0)).

4.3 The Boundary Problem

4.3.1 Introduction

In this section we present some approaches to the boundary problem, and it contains as much conjectures as theorems.

Recall that the fundamental Theorem 42 works only under condition that the trajectory does not approach the boundary. Being able to show that this condition does hold in certain cases, we do not have at the moment a proof (or a refutation) showing that the theorem holds without any extra assumption but the full support of p(0). The major conjecture is that the above condition follows from this assumption as well.

Here we would like to discuss a certain direction towards the proof of his conjecture. The idea is to use a special kind of a potential function, which will prevent the trajectory of an inner point to converge to the boundary. The Entropy, being a potential function, fails to do so. (Notice, however, that since for the given values of invariants the Entropy is maximized in an inner point π, the points in some open domain containing π cannot leave this domain).

Given a symmetric quadratic dynamical system Q on n types, we identify distributions with points of the n − 1-dimensional simplex S_{n}. Suppose we could construct a family of n − 1-dimensional bodies {B_{t}}_{t∈[0,1]} such that

1. B_{0} is the unit simplex S_{n}, B_{1} is its center (a single point);

2. B_{t} ⊂ interior(S_{n}) for all t > 0; lim_{t→0} B_{t} = S_{n}.

3. For all P ∈ B_{t}, t ∈ [0, 1], M(P) ∈ B_{t} as well.

The non-convergence to the boundary would be established quite easily: Since any interior point of S is included in some B_{t} with t > 0, by definition its trajectory lies in B_{t} as well, and thus cannot attain the boundary.

In what follows we attempt to construct (or, rather, characterize) such a family of bodies for a sufficiently large class of our quadratic systems. Systems in this class have two reasonable restrictions on β (besides the symmetry):

• For all i, β_{iiii} = 1;

• For all i and j, β_{ijij} + β_{ijji} ≥ 1 − \frac{1}{n^{2}} .

First restriction is introduced mainly for our convenience, and can probably be removed. The reasons for the second restriction will be brought later. Note, however, that any system can be altered to satisfy this restriction − on expense of n^{2} slowdown.

Note also that we don’t have to bother anymore about the invariants and the problems they introduce.

A very good candidate for an universal family would be

E_{t} = {p ∈ S_{n} | log(n)^{−1}H(p) ≥ t}  (4.6)

where H(p) is the entropy function. Unfortunately, there is a problem with this definition: condition (2) above does not hold for t ≤ \frac{log(n−1)}{log(n)} , and E_{t} ⊄ interior(S_{n}).

4.3.2 The Difference Vector

Given p, what can be said about the geometric locus of next stage distributions p × p (for all possible β -s)? The difference vector p × p − p turns out to have some interesting properties. But first some definitions.

Let V_{ijkl} ∈ ^{n} be a vector with i-th and j -th entries equal −1, k-th and l-th entries 1, and the rest 0.

Given a distribution p, define Δ_{ijkl} = sign(p_{i}p_{j} − p_{k}p_{l})V_{ijkl}.

The difference vector p × p − p is a nonnegative linear combination of Δ-s:

p × p − p = ∑_{i,j,k,l} β_{ijkl} p_{i}p_{j}V_{ijkl} = 4 ∑_{{{i,j},{k,l}}} β_{ijkl} |p_{i}p_{j} − p_{k}p_{l}| Δ_{ijkl}  (4.7)

The cone of x_{1}, .., x_{k} ∈ ^{n} is a set of all vectors y representable as y = ∑ α_{i}x_{i}, where all α-s are nonnegative. As we have seen, the cone of Δ-s includes, and, in fact, coincides with all possible directions of the difference vector. Our first observation is

Proposition 46 For each p ∈ S_{n}, the cone of vectors {Δ_{ijkl}} lies on one side of some hyperplane.

Proof: (Using the ideas of the original Bolzmann’s proof): For an interior point p of S_{n}, let Log(p) ∈ ^{n} be an n-dimensional vector whose i-th coordinate equals log(p_{i}). We show that for any i, j, k, l it holds ⟨Log(p), Δ_{ijkl}⟩ ≤ 0, and therefore the whole cone lies on the negative side of the hyperplane defined by ⟨Log(P), x⟩ = 0.

Indeed, without loss of generality assume that p_{i}p_{j} ≥ p_{k}p_{l}. Then

⟨Log(p), Δ_{ijkl}⟩ = − log p_{i} − log p_{j} + log p_{k} + log p_{j} = log (\frac{p_{k}p_{l}}{p_{i}p_{j}} ) ≤ 0

A closer look at the proof reveals that more can be said. For p ∈ S_{n} define a cell of p as the set of all (interior) points q in S_{n} with a property that for all i, j, k, l q_{i}q_{j} − q_{k}q_{l} and p_{i}p_{j} − p_{k}p_{l} have always the same sign (to zero corresponds zero). Now the proposition (46) can be strengthened as follows:

Proposition 47 Let C_{p} denote the cone corresponding to p, and let H_{q} be the half-space {x | ⟨Log(q), x⟩ ≤ 0}. Then C_{p} is the closure of intersection of H_{q} -s, where Q runs over all interior points in the cell of p.

An interesting property of this cone is that C_{p} with apex placed at p (e.g. p + C_{p}, containing in particular all possible p × p) contains the center of the simplex S_{n}.

Proposition 48 For every p ∈ S_{n}, the center of S_{n} belongs to the cone p + C_{p}.

Proof: Let U = (\frac{1}{n}, \frac{1}{n}, .., \frac{1}{n}) be the center of S_{n}. Using proposition (47), the statement of proposition reduces to

∑ p_{i} log q_{i} ≥ \frac{1}{n} ∑ log q_{i}  (4.8)

for any two p and q, interior points in the same cell. Define p̃ (and, similarly, q̃) as distribution on the pairs of types, with p̃_{ij} = p_{i}p_{j}. By a well known property of products of probability spaces,

∑_{i} p_{i} log q_{i} = \frac{1}{2} ∑_{i,j} p̃_{ij} log q̃_{ij}  (4.9)

Now the assumption that p and q are in the same cell implies that q̃_{ij} and p̃_{ij} are ordered (as real numbers) in exactly the same order. Performing small order-preserving perturbations on p̃_{ij}, we see that the quantity on the right hand side of equation (4.9) is minimized when all p̃_{ij} -s are equal. Thus

∑_{i} p_{i} log q_{i} = \frac{1}{2} ∑_{i,j} p̃_{ij} log q̃_{ij} ≥ \frac{1}{2} \frac{1}{n^{2}} ∑_{i,j} log q̃_{ij} = \frac{1}{n} ∑ q_{i}

Before finishing this subsection, we bring one simple but useful observation:

Observation The vector β_{ijkl} |p_{i}p_{j} − p_{k}p_{l}| Δ_{ijkl}, appearing in (4.7) is always directed (from P) towards the manifold {X | x_{i}x_{j} = x_{k}x_{l}}, but never crosses it. Also, the line p + αΔ_{ijkl} intersects this manifold at a single point given by α = \frac{|p_{i}p_{j}−p_{k}p_{l}|}{p_{i}+p_{j}+p_{k}+p_{l}}.

4.3.3 More about the Universal Families

Let us restrict ourselves to universal famillies of a special kind: the bodies will be convex and symmetrical (e.g. closed under the permutations of variables). B_{t}-s must now satisfy some local and global restrictions:

• Local Requirements: For a p ∈ boundary(B_{t}), the cone p + C_{p} must lie on the correct side of any supporting hyperplane of B_{t} at p. Since B_{t} is convex, it is enough to require that any Δ_{ijkl} is directed into B_{t} from p (e.g. only 2\binom{n}{4} conditions). Note that when the cone is strictly (locally) interior to B_{t}, some quantitative bounds on the change of t can be obtained. For systems with continuous time, the local requirements are both necessary and sufficient.

• Global Requirements In the case of discrete time, the danger that for a p ∈ B_{t}, p × p might jump out remains, even if the local requirements are satisfied. The following theorem shows that in our setting it cannot happen:

Proposition 49 Assuming that B_{t}-s are convex, and β -s satisfy the second restriction (here it comes in!), the local requirements imply the global ones.

Proof: Since the sum of all nontrivial β -s is ≤ 1, and B_{t} is convex, it is enough to require that for p ∈ B_{t}, p + |p_{i}p_{j} − p_{k}p_{l}| Δ_{ijkl} belongs to B_{t} for all i, j, k, l (as in the local case). Now, suppose that one of these vectors jumped out of B_{t}, crossing the boundary of the body at the point q. We claim that the local conditions are not satisfied at q. Indeed, keeping in mind the last observation of the previous subsection, one sees that the vector q + |q_{i}q_{j} − q_{k}q_{l}| Δ_{ijkl} coincides (locally) with p + |p_{i}p_{j} − p_{k}p_{l}| Δ_{ijkl}, and has the same direction. Thus, it must leave the body.

The following simple observation provides an alternative (equivalent) way of checking the local/global requirements. Let Proj_{(ijkl)} be a linear projection from ^{n} to ^{n−1} with kernel spanned by V_{ijkl}. Then, for all i, j, k, l it must hold

Proj_{(ijkl)} (B_{t} ∩ {X | x_{i}x_{j} = x_{k}x_{l}}) = Proj_{(ijkl)}(B_{t})  (4.10)

In view of the global requirements above, the following conjecture is rather interesting: Given p, an elementary operation on p is p → p + α|p_{i}p_{j} − p_{k}p_{l}| Δ_{ijkl} for some i, j, k, l and 0 ≤ α ≤ 1.

Conjecture 1 Starting from some interior point p ∈ S_{n}, and performing only elementary operations, it is impossible to get arbitrary close to the boundary of S_{n}.

It is our impression that conjecture (1) is more tractable than the original one, and it is possible either to prove it or to produce a counterexample (killing the whole approach of universal families). Also, it seems that once it is resolved, the solution of the original problem must be close at hand.

Perhaps a modification of the method used in the proof of following proposition could work in the general case:

Proposition 50 Let n = 5, and p any point with H(p) > 1. Then the conjecture (1) holds.

Proof: Let ϵ be very small. The proof is by contradiction, and goes in two stages:

First, we claim that at the moment ϵ was created, there already existed two coordinates with values ≤ ϵ^{.5}. The claim is verified by simple case analysis.

Second, we claim that at the moment that this situation (two coordinates with sum ≤ 2ϵ^{.5}) first occurred, there already existed three coordinates that summed up to ≤ 3ϵ^{.25}. Again, a simple case analysis is required.

Since elementary operations increase the entropy, we arrive at contradiction.

Remark: This method can also be applied directly to any symmetric quadratic system, showing that some set positive measure (those with entropy close to the maximal possible value) cannot approach the boundary. Another approach would be a direct construction of an universal family by an intelligent (and exceptionally lucky) guess. The following demonstration shows that this tactics sometimes works:

Proposition 51 For n = 4, and the following family is universal:

B_{t} = Symm (\frac{1}{(1 + t)^{2}} ; \frac{t^{2}}{(1 + t)^{2}} ; \frac{t}{(1 + t)^{2}} ; \frac{t}{(1 + t)^{2}} )  (4.11)

where Symm(P) is the set of points obtainable from p by permutations of coordinates.

Proof: The picture of B_{t} would be very instructional; but, since we don’t really need it for the proof, we bring here only a verbal description. It is a sort of truncated simplex S, with four (triangular) facets corresponding to vertices of S, four (triangular) facet corresponding to facets of S, and six (quadralitteral) facets corresponding to the edges of S. Each quaralitteral facet is spanned by two Δ-s. The manifold {X | x_{i}x_{j} = x_{k}x_{l}} cuts the boundary of B_{t} at four quadralitteral facets, by straight lines. It desects the body into two parts.

To obtain the proof, one has only to check that the equations (4.10) hold. Which is readily done.

Chapter 5 Matchings in Graphs

Based on a joint work with Alistair Sinclair and Avi Wigderson

Extending the corresponding part of [15], we define and study in this chapter a natural quadratic system on probability distributions over the matchings in a graph. The motivation comes both from its inherent combinatorial interest and its application to algorithms for generating matchings at random from certain probability distributions. We show that for this system the trajectory of a full support p(0) converges to an inner stationary point. Unfortunately, the most important ingredient needed to design of effective algorithms, the good bounds on convergence rate of such systems for general graphs, is still missing.

However, we give a complete analysis of this system when the underlying graph is a tree. The analysis is in two parts. The first part gives an explicit closed-form formula for the limit distribution in terms of simple parameters of the initial distribution, and seems to be of an independent combinatorial interest. The second, and more significant part is a proof that the system always converges to its limit in time that depends only polynomially on the size of the tree. This is a strong result in view of the fact that the number of matchings in the tree is in general exponentially large. Moreover, it appears to be the first bound of its kind for a complex nonlinear system defined on combinatorial structures.

(We note here that elegant global convergence rate bounds have recently been obtained, by completely different techniques, for the Boltzmann equation [4]; these are considerably weaker than our bound, but apply to a physically realistic system.)

In the last section of this chapter we discuss an interesting property of a generalized crossover, which, in principle, makes it liable to a quantitative analysis along essentially the same lines as for the matchings in trees. Here is, perhaps, the place to express our regret that many interesting combinatorial systems (e.g. on subtrees, paths and cliques in graphs) remain outside of this thesis. Consider, for instance, the following operator on subtrees of a given graph: Let T_{1}, T_{2} be a pair of such trees. Their offspring are obtained by first deciding (randomly) who is going to be the “father”, and who − the “mother”. Then the father chooses a random rib (of his), and adds it to the mother. The mother chooses a random rib from the cycle created, such that this rib will restore the father’s connectivity. It is not hard to see that such a rib always exists, and (under some assumptions) is never unique. The described operator can be shown to be symmetric, and gives rise to a rather interesting system. We hope that in the future the time will permit us to study this and others such systems in detail.

5.2 Matchings in general graphs

We define below a symmetric quadratic operator that acts on probability distributions over matchings in graphs. This operator is a natural one to choose as part of an algorithm for generating a random matching in a graph from a certain probability distribution. Applications of such algorithm would include approximate counting of matchings and computation of the partition function of monomer-dimer systems [11], as well as a genetic algorithm for finding a large matching. However, this algorithm seems at the moment to be out of reach. Let G = (V, E) be an undirected graph and M(G) the set of matchings in G. Given matchings i, j ∈ M(G), construct a random pair of offspring k and l according to the following scheme. Let G′(i, j) = (V, E′) be the subgraph of G whose edge set E′ is the symmetric difference of the edge sets of i and j. Thus each connected component of G′(i, j) is a path or a cycle whose edges belong alternately to i and to j. For each such component independently, flip a fair coin; if the coin comes up heads, make k agree with i and l agree with j on the component; if it comes up tails, reverse the roles of i and j. Also, include in both k and l all edges common to both i and j, and no other edges. (Thus the unions of the edge sets of i and j and of k and l, viewed as multisets, are equal.) It is a simple matter to check that the above scheme specifies a symmetric, aperiodic operator as defined in the previous section. We shall refer to this operator as β^{match}.

First let us consider the stationary populations under β^{match}. Using Theorem 43 of the previous chapter, we can obtain a simple description of these. We call a population p on M(G) normal if it is of the form

p_{i} = \frac{1}{Z} ∏_{e∈i} λ_{e},  (5.1)

where λ_{e} ∈ ^{+} for each edge e in G, and Z = ∑_{i∈M(G)} ∏_{e∈i} λ_{e} is a normalizing factor. In other words, in a normal population the probability of a matching is proportional to the product of the weights λ_{e} of edges e in the matching. Observe a remarkable property of normal populations: for every fixed matching M in G, the probabilities assigned to the submatchings of M (up to normalization) are as if we flipped independent (biased) coins for every edge in M. The reader might find this consistent with the fact that stationary populations maximize entropy. Actually, a suitable generalization of this fact is true in the general setting as well; we omit the details.

Theorem 52 Under the operator β^{match}, a population with full support is stationary iff it is normal.

Proof: It is easy to check that any normal population is stationary. For the other direction, let p be a stationary population on M(G), and for each edge e of G define λ_{e} = p_{{e}}/p_{0}, where 0 ∈ M(G) denotes the empty matching and {e} ∈ M(G) denotes the matching consisting of the single edge e. Since we are assuming that p has full support, each λ_{e} is positive and finite. Now we claim that, for each matching i, p_{i} = p_{0} ∏_{e∈i} λ_{e}; note that this is enough to establish the theorem. We proceed by induction on the size of i. The claim is trivially true when the size of i is 0 or 1. Suppose i = j ∪ {f}, where the edge f does not belong to matching j. Then it is easy to check that β_{i0j{f}}^{match} > 0, so Theorem 43 of the previous chapter implies that p_{i}p_{0} = p_{j}p_{{f}} = λ_{f}p_{0}p_{j}. The claim now follows by applying the inductive hypothesis to j. This concludes the proof of the claim and the theorem.

Following the approach of the previous chapter, we identify a family of invariants of β^{match}. Let e be any edge of G, and define q_{e} = ∑_{i∋e} p_{i}, the probability that a matching in population p contains e. It is easy to check that, for every e ∈ E, q_{e} is an invariant of β^{match}. Moreover, q_{e}-s form a complete basis for the linear space of invariants. Indeed, by Theorem 44 of the Chapter 4, and and the characterization of the previous theorem, any invariant linear form looks as follows:

ι_{π}(p) = ∑_{l∈N}(lg π_{l}) p_{l} = ∑_{l∈M(G)} ∑_{e∈l} lg(λ_{e})p_{l} = ∑_{e∈G} lg(λ_{e}) ∑_{l,e∈l} p_{l} = ∑_{e∈G} lg(λ_{e})q_{e}(p)

What about stationary populations with deficient support? The following surprising theorem asserts that the support of such distribution has a most interesting structure.

Given a distribution p, let S = support (p). The mating operator × is extended to supports of distributions in the obvious manner (and becomes purely combinatorial). For an edge e ∈ E define S_{e}^{+} and S_{e}^{−} the matchings in S that contain e and do not contain e, respectively. Clearly, S is the disjoint union of these two sets.

Theorem 53 If S is stationary, then at least one of the following four conditions holds:

1. S = M(G).

2. Some vertex in V is covered by all matchings in S.

3. For some (odd) subset U of V , every matching in S has exactly \frac{|U|−1}{2} inner edges (e.g. with both endpoints) in U.

4. Some egde in E is not covered by any matching in S.

Remark: Note the similarity of the last three conditions to the three types of facets defining the matching polytope. Thus, Theorem 53 can be equivalently restated as S is either all the vertices of the matching polytope, or is included in some facet. In fact the proof of our theorem is based on the proof of the famous theorem of Edmonds characterizing the facets of the Matching Polytope (as appearing in [13]).

Proof: First notice that without loss of generality, G can be assumed to be connected. Otherwise S is clearly (by its stationarity) the Cartesian Product of its projections on the connected components, and it suffices to consider each component separately.

Assume for contradiction that the conditions (1),(2),(3) and (4) do not hold. We first show that there is an edge e such that some matching in S_{e}^{−} misses both its endpoints. As (3) is false for U = V in particular, and (2) is false, there must be a pair of vertices u, v which are missed by some matching M, and choose the closest pair for which this holds. If (u, v) ∈ E we are done. Otherwise let w be a vertex on a shortest path connecting them (recall that G is connected). Clearly M covers w by our choice. Let N be any matching which does not cover w (there must be such a matching, since (2) is false). Let A be the alternating path in M ∪ N starting at w. Being a path, and as M misses both u and v, it cannot cover both of them. So consider the matching K = (N ∩ A) ∪ (M − A).

As S is stationary, K ∈ S. It misses either the pair u, w or w, v, in contradiction to our choice of closest pair.

So we conclude that for some edge e ∈ E and matching K we have M ∈ S_{e}^{−} and M misses both endpoints of e. But notice that by stability of S this implies that for every matching N ∈ S_{e}^{+}, the matching (N − e) belongs to S_{e}^{−}. It is so since in K ∪ N, e constitutes a separate component; therefore the edges of N without e form a matching in K × N. Conversely, for every K ∈ S_{e}^{−} that misses both endpoints of e, we have (K + e) ∈ S_{e}^{+}. Observe that, by our assumption, (4) does not hold, and thus S_{r}^{+} is not empty. We deduce that both S_{e}^{−} and S_{e}^{+} are stationary.

Now we can use induction on the number of edges in G. By induction S_{e}^{−} (on G − e) satisfies one of the conditions (1), (2), (3) or (4), and from the previous paragraph it follows that S will satisfy exactly the same condition (with same vertex in (2), odd set in (3) and edge in (4)).

Observe that there exists 1 − 1 correspondence between the possible values of invariants (q_{e}-s, and, subsequently all the linear invariants) and points of the Matching Polytope of G. The following theorem shows that this correspondence has additional properties. Let and initial p(0) be an initial distribution on M, and let S(0) and S(∞) denote the support of p(0) and of a limit point of the trajectory of p(0).

Theorem 54 S(∞) = M if and only if S(0) does not belong to a facet of the Matching Polytope. If it does, S(∞) belongs to the same facet. The values of invariants completely determine which of the two cases takes place.

Proof: Since S(∞) is stationary, it is either the whole M, or satisfies one of the conditions (2),(3),(4). Equivalently, it belongs to some facet of the matching polytope. It remains to be shown that the conditions (2),(3),(4) can be deduced from the values of the invariants. Recall that q_{e}-s, the probabilities that some edge appears, are invariants of our system. Let r_{v}, v ∈ V (G), denote the probabilities that some vertex is not covered. Since r_{v} = 1 − ∑_{e,v∈e} q_{e}, r_{v} -s are also invariants.

Condition (2) is equivalent to r_{v} = 0. Condition (4) is equivalent to q_{e} = 0. Condition (3) implies that ∑_{e∈U} q_{e} = \frac{|U|−1}{2} , where the sum is taken over all inner edges of U. On the other hand, this is is the maximum possible value of this sum, and is achieved if and only if (3) holds. Thus (3) can also be deduced from the values of invariants.

Remark: In view of this proof, it is very tempting to define for an arbitrary quadratic symmetric system its polytope in the following manner: for some basis of invariants {ι_{k}}_{k=1}^{n} define a polytope ⊆ ^{n} as the set of all projections (ι_{1}(p), .., ι_{n}(p)), where p runs over all distributions. The study of properties of this polytope looks really promising. Using the new terms, the above theorem states that a point in the simplex of distributions over matchings is projected to the facet of the matching polytope if and only if this point lies in a face S of the simplex, and that S = S × S (as support).

Remark: In would be tempting also to go one step further, and to relate the stationary distributions inside a facet as to its subfaces of codimention 1. This, however, does not work. The counterexample is a set of three disjoint matchings of size 3 in K_{3,3}.

We say that p(0) generates M(G) if some point p(t) on the trajectory from p(0) has full support (i.e., if the trajectory from p(0) enters the interior of the simplex). In particular, it is not hard to see that a sufficient condition for this is that p(0) assigns positive weight to the empty matching in G and, for each edge e of G, to some matching containing e.

An immediate consequence of the above theorem is the following important corollary:

Corollary 55 For any full support generating initial population on M, its trajectory cannot approach the boundary, and therefore converges to a full support stationary population.

Proof: By Theorem 54, the trajectory of p(0) cannot approach the boundary of the simplex, for if it did, its projection to the matching polytope (the values of q_{e}-s) must have been lying on the facet, contrary to our assumption. Applying Theorem 52 we complete the proof.

We turn now to the problem of determining the normal population to which the trajectory of a given initial population p(0) converges. The family of invariants {q_{e}}_{e∈E} provides us with a family of polynomial equations satisfied by the parameters λ_{e} of our normal population. Define the matching polynomial of a graph G = (V, E) over variables \overline{x} = (x_{e})_{e∈E} by g(G, \overline{x}) = ∑_{m∈M(G)} ∏_{e∈m} x_{e}. This was exactly the normalizing factor in equation 5.1. We may write the system of equations as follows:

q_{e} = \frac{λ_{e}g(G/e, \overline{λ})}{g(G, \overline{λ})}  ∀e ∈ E,  (5.2)

where G/e denotes the graph obtained from G by removing the endpoints of e, and the values q_{e} are determined by the initial population p(0).

The invariants {q_{e}}_{e∈E} form a basis for the family of linear invariants ι_{π} (defined in Theorem 44 of Chapter 4). Therefore, by Theorem 45 they uniquely determine the limit stationary population. We summarize this discussion as follows:

Theorem 56 Let p(0) be any initial population that generates M(G). Then the equations (5.2) have a unique solution {λ_{e}}_{e∈E} which determines the normal population π with the same values of invariants as p(0). Moreover, p(t) → π as t → ∞.

For general graphs (or even for C_{5}, the 5-cycle), the solution to these equations is not a rational function of the quantities q_{e}. However, by simulating the system (if it is possible), one may obtain an approximate solution to these equations. Thus it is of interest to ask what sort of systems of algebraic equations can be reduced to equations 5.2.

However, as we shall see in the next section, in the case that the underlying graph G is a tree these equations do have a simple closed-form solution. More significantly, we shall also show in this case that the quadratic system defined by β^{match} converges to a normal population very rapidly, in time bounded by a low-degree polynomial in the size of the tree.

5.3 Matchings in trees

In this section, we present a complete analysis of the mating operator β^{match} in the case where the underlying graph G is a tree T. For convenience, in what follows it will be assumed that the initial population p(0) generates M(T), which implies Theorem 56. However, a closer examination of the results of this section reveals that many of them (and in particular the convergence) can be extended to systems without full support.

5.3.1 More about the Stationary Populations

First, we give an explicit formula for the stationary population corresponding to any initial population which generates M(T). This formula will be in terms of the invariants q_{e}-s and r_{v} -s, where, as usual, q_{e} = ∑_{i∈M(G),e∈i} p_{i}, the probability of appearance of edge e, and r_{v} = 1 − ∑_{e=(v,u)} q_{e}, the probability that v is not covered by a matching.

Theorem 57 Let T be a tree and p(0) an initial population that generates M(T). Then p(t) converges to the normal population π defined by

λ_{e} = \frac{q_{e}(1 − q_{e})}{r_{v}r_{u}},

for each edge e = (v, u) of T, where the values of q_{e} and r_{v} are determined by p(0). The population π satisfies

π_{i} = \frac{∏_{e∈i} q_{e}∏_{v ∉i} r_{v}}{∏_{e∉i}(1 − q_{e})}  ∀i ∈ M(T).

(The notation v ∉ i indicates that vertex v is not covered by the matching i.)

Proof: We already know that p(t) converges to a normal population having the same values of invariants q_{e}. Since a normal population is, by definition (see equation 5.1), of a form

π_{i} = \frac{∏_{e∈i} λ_{e}}{g(T, \overline{λ})},

it is enough to obtain expressions for λ_{e}, and for the matching polynomial with these λ_{e}-s, in terms of q_{e}-s r_{v} -s. We start with the former.

Claim 58 Let T be a tree and π(\overline{λ}) a normal distribution on M(T). Then

λ_{e} = \frac{q_{e}(1 − q_{e})}{r_{v}r_{u}},  (5.3)

for each edge e = (v, u) of T.

Observe first that

q_{e} = \frac{λ_{e}g(T/e, \overline{λ})}{g(T, \overline{λ})}; 1 − q_{e} = \frac{g(T − e, \overline{λ})}{g(T, \overline{λ})}  (5.4)

where (T/e) is a subgraph of T obtained by deleting the vertices u, v, and (T − e) is obtained by deleting the edge e. Similarly

r_{v} = \frac{g(T − v, \overline{λ})}{g(T, \overline{λ})}; r_{u} = \frac{g(T − u, \overline{λ})}{g(T, \overline{λ})}  (5.5)

Recall the multiplicativity of a matching polynomial: if G_{1}, G_{2} are disjoint components of a graph G, G = G_{1} ∪ G_{2}, then

g(G, \overline{x}) = g(G_{1}, \overline{x}) g(G_{2}, \overline{x})  (5.6)

Let T_{v} be the maximal subtree of T rooted at v and not including e. In a similar way, let T_{u} be maximal subtree rooted at u and not including e. Let also S_{v} = T_{v} − v, S_{u} = T_{u} − u.

Now, using (5.4) , (5.5) and the multiplicativity of a matching polynomial, one has:

\frac{q_{e}(1 − q_{e})}{r_{v}r_{u}} = \frac{λ_{e}g(S_{v}, \overline{λ})g(S_{u}, \overline{λ}) g(T_{v}, \overline{λ})g(T_{u}, \overline{λ})}{g(S_{v}, \overline{λ})g(T_{u}, \overline{λ}) g(S_{u}, \overline{λ})g(T_{v}, \overline{λ})} = λ_{e}

At this the first goal (expressing λ-s) is completed, and we proceed to obtain an expression for the matching polynomial.

Claim 59 Let T be a tree and π(\overline{λ}) a normal distribution on M(T). Then

g(T, \overline{λ}) = \frac{∏_{e}(1 − q_{e})}{∏_{v} r_{v}}  (5.7)

Call a subtree S of T proper if it is maximal not containing a certain edge e. Let S denote the set of all proper subtrees of T. Using (5.4) and the multiplicativity of the matching polynomial, one gets the following expression for the numinator of (5.7):

∏_{e}(1 − q_{e}) = g(T, \overline{λ})^{−|E|} ∏_{S∈S} g(S, \overline{λ})

On the other hand, with the use of (5.5), the denominator can be expressed as:

∏_{v} r_{v} = g(T, \overline{λ})^{−|V |} ∏_{S∈S} g(S, \overline{λ})

Since in a tree |E| = |V | − 1, what remains after cancelations is

\frac{∏_{e}(1 − q_{e})}{∏_{v} r_{v}} = g(T, \overline{λ})

as claimed.

So far we were concerned with expressing λ-s and p_{m}-s in the terms of q_{e}-s. But how to compute q_{e}-s given the λ-s? How to compute the matching polynomial given the λ-s? Of course it is possible to add the appropriate terms, but there can be exponentially many of them. It turns out that it can be done in a quite nice and effective manner; while postponing (in the meanwhile) this problem, we bring here a theorem that will play a key role in the related discussion:

Theorem 60 Let T be a tree and π(\overline{λ}) a normal distribution on M(T). Regard T as a rooted tree, defining a father - son relation on the vertices of T. Let u be vertex, v is one of its sons, and {v_{i}}_{i=1,...,k} a complete list of v’s sons. Then

\frac{r_{u}}{q_{(v,u)}} = λ_{(v,u)}^{−1} (1 + ∑_{i=1}^{k} \frac{q_{(v_{i},v)}}{r_{v}} )

Proof: Denote by T_{w} a subtree of T containing the vertex w and all its descendants. The first observation is

\frac{r_{u}}{q_{(v,u)}} = λ_{(v,u)}^{−1} \frac{g(T_{v}, \overline{λ})}{g(T_{v} − v, \overline{λ})}  (5.8)

(the arguments needed here are similar to those used at the proof of the theorem 5.3). On the other hand, g(T_{v}, \overline{λ}) can be expanded in the following manner:

g(T_{v}, \overline{λ}) = g(T_{v} − v, \overline{λ}) + ∑_{i=1}^{k} λ_{(v_{i},v)} g(T_{v_{i}} − v_{i}, \overline{λ}) ∏_{j≠i} g(T_{v_{j}}, \overline{λ})  (5.9)

The first term in the above expansion stands for matchings in (T_{v} − v); the i-th term in the sum stands for matchings in T_{v} containing the edge (v_{i}, v). From the multiplicativity of the matching polynomial it follows

g(T_{v} − v, \overline{λ}) = ∏_{i=1}^{k} g(T_{v_{i}}, \overline{λ})

Therefore, dividing both sides of (5.9) by g(T_{v} − v, \overline{λ}), one gets

\frac{g(T_{v}, \overline{λ})}{g(T_{v} − v, \overline{λ})} = 1 + ∑_{i=1}^{k} λ_{(v_{i},v)} \frac{g(T_{v_{i}}, \overline{λ})}{g(T_{v_{i}} − v_{i}, \overline{λ})}

To finish the proof, it remains only to interpret each term of the above sum as in (5.8):

λ_{(v_{i},v)} \frac{g(T_{v_{i}}, \overline{λ})}{g(T_{v_{i}} − v_{i}, \overline{λ})} = \frac{q_{(v_{i},v)}}{r_{v_{i}}} .

5.3.2 The rate of convergence

In order to analyse the rate of convergence of β^{match} on M(T), we split the process conceptually into a sequence of stages. Each stage is associated with a partition of M(T) into equivalence classes, with the partition becoming successively finer. During each stage, we analyse the time for the probability distribution over its classes to converge. The crucial observation is that this process on equivalence classes can be viewed as a linear system, which can then be tackled using established techniques plus some additional insight.

Let v be a leaf of T, and e = (u, v) the edge adjacent to v. Let T′ be the subtree obtained by removing e and v from T. Now partition the set M(T) into equivalence classes by placing two matchings in the same class iff they agree on T′; thus each class contains at most two matchings. Note that there is a 1-1 correspondence between equivalence classes and matchings in T′. To fix notation, let [i] denote the equivalence class of matching i ∈ M(T′), and if i contains two elements denote them i^{+} and i^{−} according to whether or not they contain e. Also, for any population p on M(T) define the population p′ on M(T′) by p′_{i} = ∑_{j∈[i]} p_{j}.

Our first claim is that the operator β^{match} acting on the equivalence classes of M(T) induces a similar mating process on M(T′).

Lemma 61 For any population p on M(T), we have (p × p)′ = p′ × p′.

Proof (sketch): Since the mating operator is bilinear, and the operation p ↦ p′ is linear, it suffices to show that, for any two matchings m_{1}, m_{2} ∈ M(T), if p_{1} and p_{2} are populations concentrated at m_{1}, m_{2} respectively, then (p_{1} × p_{2})′ = p′_{1} × p′_{2}, This fact is straightforward to check by a simple case analysis depending on whether one, both or neither of m_{1}, m_{2} contains e.

In fact, this lemma is true in a more general setting. Let T_{i} be a subtree of T, and define p′′ to be a population on M(T_{i}) induced by projection p′′_{m} = ∑_{j, m=j∩T_{i}} p_{j}.

Lemma 62 For any population p on M(T), (p × p)′′ = p′′ × p′′.

Proof: Let T_{i} ⊆ T_{i+1} ⊆ ... ⊆ T_{n} = T be a chain of trees where T is the original tree, and T_{j} is obtained from T_{j+1} by plucking a leaf. Observe that the projection p ↦ p′′ can also be viewed as a composition of n − i simple projections

M(T) → M(T_{n−1}) → ... → M(T_{i})

The proof can now be completed by n − i inductive applications of Lemma 61.

An important property of the induced system on M(T_{i}) is that for all edges f ∈ T_{i}, the induced q_{f} -s are equal to the original ones.

In the light of Lemma 61, we can understand the evolution of the mating process on M(T) inductively as follows. First we view the process as working on equivalence classes: by Lemma 61, this is equivalent to a mating process on matchings in the smaller tree T′, so we may argue inductively about its behaviour. Then, assuming that the process on equivalence classes has reached stationarity, we may reason about the effects of mating within the classes. As we shall see in a moment, this latter process is much more tractable than the original quadratic system.

Let us call a population p on M(T) pseudo-stationary if the induced population p′ on M(T′) is stationary. Our next observation is that, if we start with a pseudo-stationary population, the mating operator on M(T) may be viewed as a linear system. In fact, it defines a Markov chain on the space M(T)^{+} of those matchings in T that include the edge e.

Lemma 63 Let p be any pseudo-stationary population on M(T). Then the probability distribution on matchings in M(T)^{+} behaves under β^{match} as a time-homogeneous Markov chain with transition probabilities

α_{i^{+}l^{+}} = ∑_{jk} p′_{j}β′_{ijkl}  ∀i^{+}, l^{+} ∈ M(T)^{+}  (5.10)

where β′_{ijkl} defines the mating operator in T′.

Proof: Let l^{+} ∈ M(T). By definition of the system (and ignoring for a moment the possibility that some matching appearing in the sum may not be defined), we may write (p × p)_{l} as

∑_{[i][j][k]}(p_{i^{+}}p_{j^{+}}β_{i^{+}j^{+}k^{+}l^{+}} + p_{i^{−}}p_{j^{+}}β_{i^{−}j^{+}k^{−}l^{+}} + p_{i^{+}}p_{j^{−}}β_{i^{+}j^{−}k^{−}l^{+}})  (5.11)

Consider first the case where both [i] and [j] consist of two elements. It is not hard to check that in this case β_{i^{+}j^{+}k^{+}l^{+}} = 2β_{i^{−}j^{+}k^{−}l^{+}} = 2β_{i^{+}j^{−}k^{−}l^{+}} = β′_{ijkl}. Using these relationships and the fact that p_{i^{−}} = p′_{i} − p_{i^{+}} , we may rewrite the corresponding term in the sum (5.11) as \frac{1}{2}(p′_{i} β′_{ijkl})p_{j^{+}} + \frac{1}{2}(p′_{j} β′_{ijkl})p_{i^{+}} . The case where [i] consists of two elements and [j] of one contributes only a single term to (5.11), which may be written as p_{i^{+}}p_{j^{−}}β_{i^{+}j^{−}k^{−}l^{+}} = (p′_{j} β′_{ijkl})p_{i^{+}} . The remaining case ([i] consists of one element and [j] of two) is symmetrical.

Bringing all these terms together we obtain

(p × p)_{l^{+}} = ∑_{i^{+}}(∑_{jk} p′_{j} β′_{ijkl})p_{i^{+}} = ∑_{i^{+}} α_{i^{+}l^{+}}p_{i^{+}}.  (5.12)

Since q_{e} is invariant we have ∑_{l^{+}}(p × p)_{l^{+}} = ∑_{i^{+}} p_{i^{+}} , so (5.12), suitably normalized, does indeed define a Markov chain as claimed.

Suppose now that p is a pseudo-stationary population on M(T). Then clearly p′_{i} = ∑_{j∈[i]} π_{j}, where π is the normal population corresponding to p. Therefore, by normality of π_{j}, the ratio \frac{π_{j^{+}}}{p′_{j}} does not depend on j, and equals \frac{λ_{e}}{1+λ_{e}}. Now it is easy to verify the following:

Theorem 64 The Markov chain defined in Lemma 63 is reversible and converges to the stationary distribution π̂_{l^{+}} = π_{l^{+}}/q_{e} over M(T)^{+}.

Proof: We need to show that π̂_{i^{+}}α_{i^{+}l^{+}} = π̂_{l^{+}}α_{l^{+}i^{+}} for any i, l. Since π̂_{i^{+}} and π̂_{l^{+}} are proportional, respectively, to p′_{i} and p′_{l}, this is equivalent to establishing

p′_{i}(∑_{jk} p′_{j} β′_{ijkl}) = p′_{l}(∑_{kj} p′_{k} β′_{lkji})

We claim that the (j, k)-th term in the first sum equals the (k, j)-th term in the second one. Indeed, β′_{ijkl} = β′_{lkij} by symmetricity of β′, and p′_{i}p′_{j} = p′_{l}p′_{k} whenever β′_{ijkl} ≠ 0, since p′ is normal.

We now address the question of the rate of convergence of the above Markov chain. Let Δ(t) denote the variation distance of the distribution at time t from π̂ maximized over initial states, and for ϵ ∈ (0, 1) define

τ(ϵ) = min {t : Δ(t′) ≤ ϵ for all t′ ≥ t} .  (5.13)

A bound on τ(ϵ) may be obtained by viewing the chain as a network whose vertices are the states and which has an oriented edge of capacity C(i^{+}, l^{+}) = π̂_{i^{+}}α_{i^{+}l^{+}} between vertices i^{+} and l^{+}. Suppose it is possible to route π̂_{a^{+}}π̂_{b^{+}} units of a commodity (a^{+}, b^{+}) from a^{+} to b^{+}, for each pair of distinct states a^{+} and b^{+} simultaneously, in such a way that no oriented edge (i^{+}, l^{+}) carries more than ρC(i^{+}, l^{+}) units of flow in total and the length of any flow-carrying path is at most ℓ. Then Corollary 6′ of [24] provides the following bound (see also [6, 11]):

Theorem 65 [24] The quantity τ(ϵ) satisfies

τ(ϵ) ≤ ρℓ (ln π̂_{min}^{−1} + ln ϵ^{−1}) ,

where π̂_{min} = min_{i^{+}} π̂_{i^{+}} .

The key to rapid convergence is the construction of a good flow in our Markov chain:

Theorem 66 For the Markov chain of Lemma 63, there exists a flow with ρ ≤ 2q_{e}^{−1} and ℓ = 2.

Proof: Note from (5.10) that we may view each edge (i^{+}, l^{+}) of the Markov chain as a set of parallel edges, one for each j ∈ M(T′): the edge associated with j has capacity π̂_{i^{+}}p′_{j}β′_{ijkl} and corresponds to mating i^{+} with j^{+} and producing l^{+} as an offspring. We refer to this edge as the j -edge from i^{+} to l^{+}. We will describe our flow in terms of these edges.

Let a^{+}, b^{+} be any pair of distinct states in M(T)^{+}, and consider the set P_{ab} of paths of length 2 from a^{+} to b^{+} of the following form: the first edge of the path is any b-edge leaving a^{+}, to c^{+} say, and the second edge is the unique \overline{c}-edge from c^{+} to b^{+}, where \overline{c} is defined by c ∪ \overline{c} = a ∪ b (viewed as multisets). We distribute the π̂_{a^{+}}π̂_{b^{+}} units of flow from a^{+} to b^{+} evenly over all 2^{r} paths in P_{ab}, where r is the number of non-trivial connected components in a^{+} ∪ b^{+}.

Now consider any j -edge e from i^{+} to l^{+}. We may assume that the equivalence class [j] contains two elements since in the above flow we have only used j -edges for which j^{+} ∈ M(T)^{+}. Let us compute the flow through e. Clearly only two flow-carrying paths use e; for if e is the first edge of a path then this path must belong to P_{ij}, while if it is the second edge of a path then this path must belong to P_{ml}, where m is the complement of l in the multiset i ∪ j. Suppose e is the first edge of a path from i^{+} to j^{+}. The flow along this path is precisely 2^{−r}π̂_{i^{+}}π̂_{j^{+}} , where r is the number of non-trivial connected components in i^{+} ∪ j^{+}. Dividing this by the capacity of e we get

\frac{2^{−r}π̂_{i^{+}}π̂_{j^{+}}}{π̂_{i^{+}}p′_{j}β′_{ijkl}} = \frac{2^{−r}π_{j^{+}}/q_{e}}{(π_{j^{+}} + π_{j^{−}})2^{−r}} ≤ \frac{1}{q_{e}} ,

where we have used the fact that β′_{ijkl} = 2^{−r}. By a symmetrical argument, the ratio of flow to capacity in the case where e is the second edge on a path can be similarly bounded. Hence ρ ≤ 2q_{e}^{−1}. Evidently all flow-carrying paths have length 2 so ℓ = 2.

Finally, we may combine the above results to derive a bound on the rate of convergence of the original quadratic mating process. Define the quantity τ^{match}(ϵ) in analogous fashion to τ(ϵ) in (5.13) but using the variation distance of p(t) from the corresponding stationary population π defined in Theorem 57.

Theorem 67 Let T be a tree with n vertices, and p(0) any population that generates M(T). Then τ^{match}(ϵ) satisfies

τ^{match}(ϵ) = O((n lg n + n lg π_{min}^{−1} + lg ϵ^{−1})q_{min}^{−1}),

where q_{min} = min_{e} q_{e} and π_{min} = min_{i∈M(T)} π_{i} .

Notice that, by Theorem 57, π_{min} ≤ qmin^{\frac{n}{2}} r_{min}^{n}.

Proof: The proof is somewhat tedious. First we provide a sketch: The chain of subtrees T_{1} ⊆ T_{2} ⊆ ... ⊆ T_{n} = T, induces, by Lemmas 61,62, a chain of induced systems on these subtrees. The induced system on T_{1} is stationary; assume by inductive hypothesis that after sufficiently many steps, the induced system on M(T_{i−1}), is close to stationarity. By Lemma 63, it implies that the induced system on M(T_{i})^{+} is, roughly, a Markov chain. The transition probabilities appearing in (5.10) are perturbed by a small time-varying amount. However, a pipelining effect occurs, and the errors do not accumulate. Applying Theorem 65 together with the flow bound of Theorem 66, we obtain a bound on the rate of convergence of this Markov chain. This in turn bounds the distance from the stationary state of the induced system on M(T_{i}). (Since the values of M(T_{i})^{+} uniquely define the values of M(T_{i}), by the structure of the equivalence classes of M(T_{i})). Continuing in a similar fashion, we arrive eventually to i = n, and bound the global convergence time.

The norm used throughout this proof is the variation norm, so that for p, q distributions ∥p − q∥ = Δ(p; q), and for an arbitrary vector v, ∥v∥ = \frac{1}{2}∥v∥_{1}. Distributions will often be considered as real valued vectors of a certain length.

In what follows, we use π_{i} and p_{i}(t) (without risk of confusion with the probability of matching i in p(t)) for the stationary distribution and the induced population at time t on M(T_{i}). Instead of considering the p_{i}(t), we shall usually consider δ_{i}(t) = p_{i}(t) − π_{i}. V_{i} is defined as a vector space on M(T_{i}) with coordinates summing up to 0. Clearly, δ_{i}(t) ∈ V_{i}.

The linear projection s_{i} : V_{i} ↦ V_{i−1} is defined by summing over each equivalence classes of M(T_{i}). That is, for v ∈ V_{i} and l ∈ M(T_{i−1}), the l-th coordinate of s_{i}(v) is

s_{i}(v)_{l} = ∑_{m∈[l]} v_{m}

As we have seen (Theorems 61,63), if the induced system on M(T_{i−1}) is in a stationary state, the induced system on M(T_{i})^{+} is a Markov chain. Consequently, the induced system on M(T_{i}) can be described by a linear operator. This operator will be called η̄_{i}. Clearly η̄_{i}(π_{i}) = π_{i}, and η̄_{i} is defined over ker(s_{i}) ⊆ V_{i}. Later we shall extend it to all V_{i}, while preserving its essential properties.

Claim 68 Let v ∈ ker (s_{i}), a nonzero vector in V_{i}. Then

\frac{∥η̄_{i}^{t}(v)∥}{∥v∥} ≤ \begin{cases}1 & always \\ e^{(−_{\frac{t−c}{4q{−1 \atop min}}} )} & if t > c \end{cases}

where c = 4q_{min}^{−1} ln(π_{min}^{−1}).

Proof: Recall that by Theorem 63, the induced system on M(T_{i})^{+} (normalized by its total weight, q_{e}, e = T_{i} − T_{i−1}) is a rapidly converging Markov Chain. Therefore, the whole system on M(T_{i}) is also a rapidly converging linear system. The time needed to η̄_{i} to approach π_{i}, the stationary distribution of M(T_{i}), within distance ϵ, is bounded (by Theorem 65 with the flow bound as in Theorem 66) by 4q_{e}^{−1}(ln(\frac{q_{e}}{π_{min}} ) + ln(\frac{q_{e}}{2ϵ}) ). Restating this fact and replacing q_{e} by q_{min}, the second part of the inequality follows. The first part is simply a consequence of the convexity of the variation norm.

Next we extend η̄_{i}, which is defined on ker(s_{i}), to a linear operator η_{i} : V_{i} ↦ V_{i}. The Claim 68 will hold for the extension as well. Let us first define a linear mapping r_{i} from V_{i} to ker(s_{i}). Consider a vector v in V_{i}. For all coordinates k such that the corresponding equivalence class in M(T_{i−1}) consists of this coordinate alone, the k-th coordinate of r_{i}(u) will be zero. The remaining coordinates are paired according to their equivalence class. Let {j, k} be such a pair of coordinates. The corresponding coordinates of r_{i}(u) will be defined as \frac{u_{j}−u_{k}}{2} and \frac{u_{k}−u_{j}}{2} , respectively. An easy check shows that r_{i} does indeed map V_{i} to ker(s_{i}). The η_{i} is now defined as

η_{i}(v) = η̄_{i}r_{i}(v)

Claim 69

• η_{i} coincides with η̄_{i} on ker(s_{i}).

• Claim 68 continues to hold when η̄_{i} is replaced by η_{i} and v is any vector in V_{i}.

Proof: It is enough to check that r_{i} is an identity on ker(s_{i}), and that

∥r_{i}(v)∥ ≤ ∥v∥  (5.14)

for any v ∈ V_{i}.

We start now with the main part of the proof. Let, as before, p_{i}(t) be the induced distribution on M(T_{i}) after t steps, π_{i} the limit stationary distribution on M(T_{i}), and δ_{i}(t) = p_{i}(t) − π_{i}, an error term. Our main goal is to bound ∥δ_{n}(t)∥.

Consider δ_{i}(t + 1): where does this error term originate? If δ_{i−1}(t) = 0 (i.e. p_{i−1}(t) is already stationary), then δ_{i}(t + 1) = η̄_{i}δ_{i}(t). Otherwise there will be an additional term reflecting the presence of “noise” δ_{i−1}(t). This term will be expressed using an “operator” γ = γ_{i,t}:

γ_{i,t+1}δ_{i−1}(t) ≡ δ_{i}(t + 1) − η_{i}δ_{i}(t)

Actually, γ is but a shorthand, and is defined only within our framework. Its main property is

Claim 70 ∥γδ∥ ≤ 4∥δ∥

Proof: First, we estimate the Lipschitz constant of β^{match}; in fact the statement holds for any QDS. Let p, q be distributions on types of some QDS, and × the mating operator in this system. The following inequality provides a bound for sensitivity of × to small perturbations:

∥q × q − p × p∥ ≤ 3∥q − p∥  (5.15)

Here is a proof: The operator β^{match} can be extended to a bilinear operator from V × V to V, where V is a real vector space with coordinates labeled by types of the QDS. Since vectors corresponding to the probability distributions form a basis for V, such extension is possible and unique. The same × will be used to denote it.

For any two vectors v, u ∈ V the properties of × on the probability distributions imply

∑_{i∈M(T^{∗})} (v × u)_{i} = ( ∑_{i∈M(T^{∗})} v_{i}) ( ∑_{i∈M(T^{∗})} u_{i})

This, together with the observation that all the coefficients β_{ij→k} of × are nonnegative, implies:

∥v × u∥ ≤ 2∥v∥ ∥u∥

Let p, q be distributions, and r = q − p. Clearly ∥r∥ ≤ 1, ∥p∥ = \frac{1}{2} . We have

q × q − p × p = (p + r) × (p + r) − p × p = 2(p × r) + r × r;
∥q × q − p × p∥ ≤ 4 ∥r∥ ∥p∥ + ∥r∥^{2} = 2 ∥r∥ + ∥r∥^{2} ≤ 3 ∥r∥

as claimed.

To establish the Claim 70 we need one more definition: Let p̂_{i}(t) be a distribution on M(T_{i}) given by

p̂_{i}(t) = π_{i} + r_{i} (p_{i}(t) − π_{i})

where r_{i} is the linear function from V_{i} to ker(s_{i}) defined above. From the definition of r_{i} it follows

∥p̂_{i}(t) − p_{i}(t)∥ = ∥π_{i−1} − p_{i−1}(t)∥ = ∥δ_{i−1}(t)∥

Now,

δ_{i}(t + 1) = p_{i}(t + 1) − π_{i} = p_{i}(t) × p_{i}(t) − π_{i} =
[p̂_{i}(t) × p̂_{i}(t) − π_{i}] + [p_{i}(t) × p_{i}(t) − p̂_{i}(t) × p̂_{i}(t)]

Call the second term E_{2}; we proceed with the first term. Since, by definition, p̂_{i}(t) − π_{i} belongs to ker(s_{i}),

p̂_{i}(t) × p̂_{i}(t) − π_{i} = η_{i} [p̂_{i}(t) − π_{i}] = η_{i} [p_{i}(t) − π_{i}] + η_{i} [p̂_{i}(t) − p_{i}(t)] =
η_{i}δ_{i}(t) + η_{i} [p̂_{i}(t) − p_{i}(t)]

Call the second term of the last expression E_{1}. The above equalities imply

γδ_{i−1}(t) = δ_{i}(t + 1) − η_{i}δ_{i}(t) = E_{1} + E_{2}

But, keeping in mind that ∥p̂_{i}(t) − π_{i}(t)∥ = ∥δ_{i−1}(t)∥ , one has

∥E_{2}∥ ≤ 3 ∥δ_{i−1}(t)∥ by (5.15)  and  ∥E_{1}∥ ≤ ∥δ_{i−1}(t)∥ by Claim 68

establishing the statement of the claim.

We continue with the proof of the theorem. Using the convention γ, one can express ∥δ_{n}(t)∥ as follows:

∥δ_{n}(t)∥ = ∥η_{n}δ_{n}(t − 1) + γδ_{n−1}(t − 1)∥ ≤ ∥η_{n}δ_{n}(t − 1)∥ + 4 ∥δ_{n−1}(t − 1)∥;

Expressing once more the δ-s in the same form, and using bounds of Claims 68,70 and the linearity of η, one has

∥δ_{n}(t)∥ ≤ ∥η_{n}η_{n}δ_{n}(t − 2)∥ + ∥η_{n}γδ_{n−1}(t − 2)∥ + 4 ∥η_{n−1}δ_{n−1}(t − 2)∥ + 4 ∥γδ_{n−2}(t − 2)∥

≤ ∥η_{n}η_{n}δ_{n}(t − 2)∥ + 4 ∥δ_{n−1}(t − 2)∥ + 4 ∥η_{n−1}δ_{n−1}(t − 2)∥ + 4^{2} ∥δ_{n−2}(t − 2)∥

It is easier to understand what happens after t stages of this expansion by associating each term with a string of length t of characters η and γ, which encode the “path” that created the term. Consider a string containing k ≤ n − 1 characters γ:

ϕ = η_{n}^{(s_{n})}γη_{n−1}^{(s_{n−1})}γ...η_{n−k}^{(s_{n−k})}

where all s_{j} ≥ 0 and ∑_{j} s_{j} = t − k. Each long block of η-s results, by Claim 68, in contraction of the norm; each γ increased the norm, by Claim 70, by at most 4. Therefore, the norm of term corresponding to this string is bounded by

4^{k} ∏_{j} e^{−(_{\frac{s_{j} −c}{4q{−1 \atop min}}} )}∥δ_{n−k}(0)∥

where the product is taken over all s_{j} of length exceeding c.

Let Φ_{k,t−k} be the set of all strings with k γ -s and t − k η-s. Using the same line of reasoning, we obtain

∥δ_{n}(t)∥ ≤ ∑_{k=0}^{n−1} ∑_{ϕ∈Φ_{k,t−k}} ∥ term corresponding to ϕ ∥  (5.16)

All δ_{j}(0) have norm ≤ 1. Taking the worst case (where all blocks of η-s have equal size) and substituting in (5.16), we get

∥δ_{n}(t)∥ ≤ ∑_{k=0}^{n−1} \binom{t}{k}4^{k}e^{−(_{\frac{t−k(c+1)}{4q{−1 \atop min}}})}∥δ_{n−k}∥ ≤ (\frac{c_{1}t}{n})^{n} e^{−b_{1}t}

where c_{1} = 5 e^{\frac{c+1}{b_{1}}} , b_{1} = \frac{1}{4}q_{min} (keeping in mind that t ≫ n).

Now, to obtain a lower bound on τ^{match}(ϵ), it suffices to find t so that

(\frac{c_{1}t}{n})^{n} e^{−b_{1}t} ≤ ϵ

Solving this inequality for t, we obtain the desired lower bound. The method of solution uses a standard technical approach; the details are omitted.

Theorem 67 is a surprisingly strong result. It says that the time required for the mating process to approach stationarity depends on the size of the tree only as a low degree polynomial, despite the fact that the number of types |M(T)| is in general exponentially large in n. This is apparently the first analysis of its kind for a combinatorially non-trivial problem. Moreover, although this is a somewhat simplified setting, we conjecture that the analytical tools introduced in this chapter are applicable to other systems.

5.3.3 Remark on the Projection Property

The technique of proof used in the last section can be exploited sometimes in analysis of different systems as well, yielding the rapid mixing. One such system is the generalized crossover (see Section 1.5 for a definition). This system was studied in the context of Populational Genetics, and (as we have recently learned, see [2]) many of the general results of the previous chapter are known for it. However, the rate of convergence, to our best knowledge, was never addressed nor solved.

The reason for certain similarity between the generalized crossover and the matchings in trees is that the two systems share the same projection property. In both cases it is possible to construct a chain of systems S_{n} → S_{n−1} → ... → S_{1} such that

• S_{n} is the original system and S_{1} is trivial.

• The dynamics of S_{i} is completely determined by that of S_{i+1}.

• Provided that S_{i} reached stationarity, S_{i+1} is linear (a Markov chain).

Remark: Observe a certain formal similarity between the definitions of the projection property in dynamical systems and solvability in groups. We do not know for the present whether it is purely coincidental, or reflects some deeper analogy.

The existence of such a chain for system on matchings of a tree is a contents of Lemmas 61,63. For the generalized crossover, such chain is obtained by considering the distributions induced by projecting on the first i bits:

{1, 2, .., n} → {1, 2, .., n − 1} → ... → {1}.

We omit the proof that indeed this is a valid chain of systems. Constructing an appropriate flow (as in Theorem 66) and using the same arguments, a nontrivial bound on the rate of convergence can be established for this system.

References

[1] M. Ajtai, J. Komlós and E. Szemerédi. Deterministic simulation in LOGSPACE. Proceedings of 19th ACM STOC (1987), p. 132.

[2] E. Akin. The geometry of population genetics. Lecture Notes in Biomathematics, Vol. 31, Springer-Verlag, Berlin, 1979.

[3] F. Brenti. Unimodal, log-concave and Polya sequences in combinatorics. AMS, Providence, R.I., 1989.

[4] E. Carlen and M. Carvalho. Strict entropy production bounds and stability of the rate of convergence to equilibrium for the Boltzmann equation. Preprint, 1992.

[5] I. Csiszár and J. Körner. Information Theory: Coding Theorems for Discrete Memoryless Systems. Academic Press, New York, 1981.

[6] P. Diaconis and D. Stroock. Geometric bounds for eigenvalues of Markov chains. Annals of Applied Probability 1 (1991), pp. 36–61.

[7] M. Dyer, A. Frieze and R. Kannan. A random polynomial time algorithm for approximating the volume of convex bodies. Proceedings of 21st ACM STOC (1989), pp. 375–381.

[8] M. Hirsch and S. Smale. Differential equations, dynamical systems, and linear algebra. Academic Press, New York, 1974.

[9] J. Holland. Adaptation in natural and artificial systems. MIT Press, 1992.

[10] M. Jerrum. The elusiveness of large cliques in a random graph. Report CSR-9-90, Dept. of Computer Science, University of Edinburgh, 1990.

[11] M. Jerrum and A. Sinclair. Approximating the permanent. SIAM Journal on Computing 18 (1989), pp. 1149–1178.

[12] J.R. Koza. Genetic Programming. MIT Press, Cambridge, Mass, 1992.

[13] L. Lovász and M.D. Plummer. Matching Theory. North-Holland, Amsterdam, 1986.

[14] O. Penrose. Foundations of statistical mechanics. Pergamon Press, Oxford, 1970.

[15] Y. Rabinovich, A. Sinclair and A. Wigderson. Quadratic dynamical systems. Proceedings of the 33rd IEEE FOCS (1992), pp. 304-313.

[16] Y. Rabinovich and A. Wigderson. Analysis of a simple genetic algorithm. Proceedings of the 4th International Conference on Genetic Algorithms 1991, pp. 215–221.

[17] G.J.E. Raulings, ed. Foundation of Genetic Algorithms. Morgan Kaufmann Publ., Sant Mateo, 1991.

[18] A. Rényi. Foundations of probability. Holden-Day, San Francisco, 1970.

[19] R.L. Riolo. Survival of the fittest bits. Scientific American 267, July 1992.

[20] A.L. Samuel. Some studies in Machine Learning using the game of chechers. In: Computers and Thought, McGraw-Hill, New York, 1963.

[21] G. Sasaki and B. Hajek. The time complexity of maximum matching by simulated annealing. Journal of the ACM 35 (1988), pp. 387–403.

[22] E. Seneta. Nonnegative matrices and Markov chains. Springer-Verlag, New York, 1981.

[23] H. Siegelmann. Personal communication.

[24] A. Sinclair. Improved bounds for mixing rates of Markov chains and multicommodity flow.Report ECS-LFCS-91-178, Dept. of Computer Science, University of Edinburgh, 1991. To appear in Combinatorics, Probability & Computing.

[25] D.W. Walkup. Polya sequences, binomial convolution and Unions of Random Sets. Journal of Applied Probability 13 (1976), pp. 76–85.