PROBABILITY AND ALGORITHMICS: A FOCUS ON SOME RECENT DEVELOPMENTS

This article presents different recent theoretical results illustrating the interactions between probability and algorithmics. These contributions deal with various topics: cellular automata and calculability, variable length Markov chains and persistent random walks, perfect sampling via coupling from the past. All of them involve discrete dynamics on complex random structures. Résumé. Cet article présente différents résultats récents de nature théorique illustrant les interactions entre probabilités et algorithmique. Ces contributions traitent de sujets variés : automates cellulaires et calculabilité, chaînes de Markov à mémoire variable et marches aléatoires persistantes, échantillonage parfait par la méthode de couplage par le passé. Leur point commun est de faire intervenir des dynamiques discrètes sur des structures aléatoires complexes. Probability and algorithmics are two domains that mutually enrich each other. On the one hand, the analysis of algorithms naturally involves the study of random discrete objects. On the other hand, in order to study and generate complex combinatorial structures, one often needs to develop sophisticated algorithms. In that context, a topical question is also the classification of probability distributions that can be described with the help of different models of computation. This article presents some recent theoretical results illustrating the interactions between probability and algorithmics. It is the reflect of the session of the Journées MAS 2016 entitled “Probability and Algorithmics”. The three sections are based on the talks given in Grenoble on this occasion by Mathieu Sablik, Peggy Cénac, Christelle Rovetta, and Rémi Varloot. Irène Marcovici, the organizer of that session of the Journées MAS 2016, has written the introduction and coordinated the present article. Section 1 is an introduction to the advances of Mathieu Sablik on the possible asymptotic measures that can be reached by a cellular automaton.a Cellular automata are used to model all kinds of phenomena with local interactions that appear in biology, physics, or theoretical computer science. S. Wolfram [52] was the first to propose a classification of cellular automata, based on the observation of space-time diagrams. His empirical classification was then followed by many other attempts, the aim being to understand the complex phenomena that can appear. Here, we focus on the typical asymptotic configurations observed from a random initial configuration. They are described by the set of limit measures. From a modelisation point of view, this 1 Institut de Mathématiques de Bourgogne, UMR 5584, Université de Bourgogne, CNRS, F-21078 Dijon Cedex, France 2 Institut Élie Cartan de Lorraine, UMR 7502, Université de Lorraine, CNRS, F-54506 Vandoeuvre-lès-Nancy Cedex, France 3 Inria Paris and LIP6, UMR 7606, Université Pierre et Marie Curie, CNRS, F-75252 Paris Cedex 05, France 4 Institut de Mathématiques de Toulouse, UMR 5219, Université de Toulouse, CNRS, Université Paul Sabatier, F-31062 Toulouse Cedex 9, France 5 Microsoft Research Inria Joint Centre, F-91120 Palaiseau, France aThis work was partially supported by the ANR project QuasiCool ANR-12-JS02-011-01 c © EDP Sciences, SMAI 2017 Article published online by EDP Sciences and available at http://www.esaim-proc.org or https://doi.org/10.1051/proc/201760203 204 ESAIM: PROCEEDINGS AND SURVEYS set corresponds to the set of measures that are “physically” relevant, starting from a simple measure such as the uniform Bernoulli measure. From a computer science point of view, cellular automata can be seen as a model of computation with high parallelism and it is natural to study the possible complexities of the asymptotic behavior, and its dependence on the initial distribution. Section 2, presented by Peggy Cénac, provides a short zoology of variable length Markov chains, with necessary and sufficient conditions for existence and uniqueness of invariant measure. Such conditions can be related to recurrence/transience properties of random walks built from Markov chains with unbounded memory. We will follow some examples illustrating the unstable nature of the random walk behavior under slightly disturbed settings. The works described in this section are extracted from two reference articles [10,11]. Section 3 is based on works of Christelle Rovetta and Rémi Varloot on perfect sampling. In [42], Propp and Wilson have introduced the coupling from the past method to generate random variables according to the stationary distribution of a given Markov chain. The advantage over the classical Monte Carlo approach is that it provides a perfect sampling algorithm: it gives an actual algorithm that terminates in finite expected time, as opposed to a converging algorithm, and the output distribution is unbiased, rather than arbitrarily close to the target distribution. This initially comes at an important cost in complexity, which becomes proportional to the size of the state space. Different strategies have been introduced so as to remove this dependency and achieve acceptable complexity. The basic building blocks for such strategies are presented, alongside two recent applications stemming from collaborations with Anne Bouillard and Ana Bušić. The first is a novel means of sampling from the stationary distribution of closed queueing networks in an efficient way [4]. The second is an approach to speed up the simulation of “very lazy” random walks, i.e. random walks for which the probability of changing states during any given transition is small [47], presented here in the context of sampling random independent sets from a graph. 1. Limit measures of cellular automata: a computational approach In this section, we present (in a simplified version) some new advances on the study of typical asymptotic measures obtained when iterating a cellular automaton from some initial measure. There are two possible approaches for the study of the sets of limit measures of cellular automata: ‚ taking a class of cellular automata having empirically the same typical asymptotic behavior, and trying to determine their sets of limit measures; ‚ characterizing sets of measures that can be reached as sets of limit measures for some cellular automata. In Subsection 1.1, we recall the main definitions and present some results about the asymptotic measures for some classes of cellular automata. Generally, one observes a behavior on space-time diagrams, and then tries to determine the set of limit measures for this class. This type of study can be difficult and often requires to develop new tools for each class that is considered. In Subsection 1.2, we are interested in characterizing the limit measures that can be reached asymptotically when iterating a cellular automaton from some initial measure. A first step is to determine some obstructions, that is, some necessary conditions that these measures must satisfy. As cellular automata are model of computation, in addition to the classical topological obstructions, some necessary computational obstructions appear. The main problem is to prove a reverse statement: given a measure satisfying the computational obstructions, we want to construct a cellular automaton which, starting from any simple initial measure, reaches this measure asymptotically. Similar computational obstructions appear in topological dynamics, when characterizing possible properties of subshifts of finite type or cellular automata: possible entropies [29], possible growth-type invariants [38], possible sub-actions [1,28]... However, the construction is quite different here, since starting from a random configuration, the construction requires an ability to self-organize the space. The results presented here are part of a collaboration with Benjamin Hellouin de Menibus [25]. In the last subsection we give some perspectives and open questions. ESAIM: PROCEEDINGS AND SURVEYS 205 1.1. Definitions and study of some classes of cellular automata 1.1.1. Problematic A cellular automaton is a complex system defined by a local rule which acts synchronously and uniformly on the configuration space AZ, where A is a finite alphabet. Formally, it is a function F : AZ Ñ AZ defined by a local rule F : Ar ́r,rs ÝÑ A, where r P N is called the radius, such that: F pxqi “ F pxri ́r,i`rsq for all i P Z. Equivalently, cellular automata are exactly continuous functions which commute with the shift map σ : AZ Ñ AZ, defined by σpxqi “ xi`1 for all x P AZ and i P Z. Example 1.1. Let A “ t0, 1u and consider the cellular automaton defined by F pxqi “ xi ́1 ` xi ` xi`1; [mod 2] for all i P Z. Given an initial configuration, a representation of the behavior of this cellular automaton is given by the superposition of the different configurations obtained by iteration of the dynamics. Figure 1.1 shows an example of such a space-time diagram, starting from a random initial configuration (represented at the bottom). Figure 1.1. Space time diagram of the cellular automaton defined by F pxqi “ xi ́1 ` xi ` xi`1 [mod 2] for all i P Z. Cellular automata are simple models but have a wide variety of different dynamical behaviors. We are interested in their typical asymptotic behavior starting from a random configuration, as observed on simulations. We consider that the space of configurations is homogeneous. This is rendered by an initial configuration chosen randomly according to a σ-invariant measure μ, meaning that μpBq “ μpσ ́1pBqq for any Borel set B. Denote by MσpAq the set of σ-invariant measures. This space is compact and metrizable for the weak ̊ topology. We consider here the following metric:

set corresponds to the set of measures that are "physically" relevant, starting from a simple measure such as the uniform Bernoulli measure. From a computer science point of view, cellular automata can be seen as a model of computation with high parallelism and it is natural to study the possible complexities of the asymptotic behavior, and its dependence on the initial distribution. Section 2, presented by Peggy Cénac, provides a short zoology of variable length Markov chains, with necessary and sufficient conditions for existence and uniqueness of invariant measure. Such conditions can be related to recurrence/transience properties of random walks built from Markov chains with unbounded memory. We will follow some examples illustrating the unstable nature of the random walk behavior under slightly disturbed settings. The works described in this section are extracted from two reference articles [10,11]. Section 3 is based on works of Christelle Rovetta and Rémi Varloot on perfect sampling. In [42], Propp and Wilson have introduced the coupling from the past method to generate random variables according to the stationary distribution of a given Markov chain. The advantage over the classical Monte Carlo approach is that it provides a perfect sampling algorithm: it gives an actual algorithm that terminates in finite expected time, as opposed to a converging algorithm, and the output distribution is unbiased, rather than arbitrarily close to the target distribution. This initially comes at an important cost in complexity, which becomes proportional to the size of the state space. Different strategies have been introduced so as to remove this dependency and achieve acceptable complexity. The basic building blocks for such strategies are presented, alongside two recent applications stemming from collaborations with Anne Bouillard and Ana Bušić. The first is a novel means of sampling from the stationary distribution of closed queueing networks in an efficient way [4]. The second is an approach to speed up the simulation of "very lazy" random walks, i.e. random walks for which the probability of changing states during any given transition is small [47], presented here in the context of sampling random independent sets from a graph.

Limit measures of cellular automata: a computational approach
In this section, we present (in a simplified version) some new advances on the study of typical asymptotic measures obtained when iterating a cellular automaton from some initial measure. There are two possible approaches for the study of the sets of limit measures of cellular automata: ‚ taking a class of cellular automata having empirically the same typical asymptotic behavior, and trying to determine their sets of limit measures; ‚ characterizing sets of measures that can be reached as sets of limit measures for some cellular automata.
In Subsection 1.1, we recall the main definitions and present some results about the asymptotic measures for some classes of cellular automata. Generally, one observes a behavior on space-time diagrams, and then tries to determine the set of limit measures for this class. This type of study can be difficult and often requires to develop new tools for each class that is considered.
In Subsection 1.2, we are interested in characterizing the limit measures that can be reached asymptotically when iterating a cellular automaton from some initial measure. A first step is to determine some obstructions, that is, some necessary conditions that these measures must satisfy. As cellular automata are model of computation, in addition to the classical topological obstructions, some necessary computational obstructions appear. The main problem is to prove a reverse statement: given a measure satisfying the computational obstructions, we want to construct a cellular automaton which, starting from any simple initial measure, reaches this measure asymptotically. Similar computational obstructions appear in topological dynamics, when characterizing possible properties of subshifts of finite type or cellular automata: possible entropies [29], possible growth-type invariants [38], possible sub-actions [1,28]... However, the construction is quite different here, since starting from a random configuration, the construction requires an ability to self-organize the space. The results presented here are part of a collaboration with Benjamin Hellouin de Menibus [25].
In the last subsection we give some perspectives and open questions.

Definitions and study of some classes of cellular automata
A cellular automaton is a complex system defined by a local rule which acts synchronously and uniformly on the configuration space A Z , where A is a finite alphabet. Formally, it is a function F : A Z Ñ A Z defined by a local rule F : A r´r,rs ÝÑ A, where r P N is called the radius, such that: Equivalently, cellular automata are exactly continuous functions which commute with the shift map σ : A Z Ñ A Z , defined by σpxq i " x i`1 for all x P A Z and i P Z. Example 1.1. Let A " t0, 1u and consider the cellular automaton defined by Given an initial configuration, a representation of the behavior of this cellular automaton is given by the superposition of the different configurations obtained by iteration of the dynamics.
Cellular automata are simple models but have a wide variety of different dynamical behaviors. We are interested in their typical asymptotic behavior starting from a random configuration, as observed on simulations. We consider that the space of configurations is homogeneous. This is rendered by an initial configuration chosen randomly according to a σ-invariant measure µ, meaning that µpBq " µpσ´1pBqq for any Borel set B. Denote by M σ pA Z q the set of σ-invariant measures. This space is compact and metrizable for the weak˚topology. We consider here the following metric: where rus " x P A Z : x r0,n´1s " u ( is the cylinder centered on the word u P A n . For simulations, the initial configuration is generally chosen randomly according to a Bernoulli measure: the probability that a cell has the value a P A is p a P r0, 1s, independently for different cells. In other words, the Bernoulli measure associated to pp a q aPA with ř aPA p a is defined by λ ppaq aPA prusq " p u0 p u1 . . . p un´1 for all u P A n . For the uniform Bernoulli measure, all the parameters are chosen to have the same value p a " 1 CardA , we denote it λ A Z . It is also possible to define measures with finite ranges of spatial dependences, thanks to Markov processes.
Another important class of measures consists in the σ-invariant measures supported by periodic orbits: for w P A n , we denote by 8 w 8 the σ-periodic configuration associated to w, and we define: Let us recall that A˚denotes the set of words on the alphabet A. Even if the measure x δ w is algorithmically simple, the set Most of the results we obtain are not restricted to Bernoulli initial measures. A measure is σ-ergodic if σ-invariant sets have measure 0 or 1. We denote by M σ´erg pA Z q the set of σ-ergodic measures. By Birkhoff's theorem, the frequency of a word u P A n in a random configuration chosen according to a σ-ergodic measure is equal to the measure of the cylinder rus (see [46] for details). Sometimes, we also need to evaluate the independence of cells in a random configuration. A measure is full support if all cylinders are charged, and we denote by M full σ pA Z q the set of σ-invariant measures of full support. Since a cellular automaton F commutes with the shift, it acts naturally on M σ pA Z q by: The measure F t µ represents the measure obtained at time t when iterating the cellular automaton F from the initial measure µ. For an initial measure µ P M σ pA Z q, we are interested in the asymptotic behavior of the sequence pF t µq tPN , and of its Cesàro mean, defined by In particular, we focus on the following sets of limit points, which are non-empty, by compactness of M σ pA Z q: ‚ the µ-limit measures set VpF, µq, defined as the set of limit points of the sequence pF t µq tPN ; ‚ the Cesàro mean µ-limit measures set V 1 pF, µq, defined as the set of limit points of pϕ F t µq tPN . The union of the support of the µ-limit measures set corresponds to the µ-limit sets introduced in [34] in a different way. Very complex µ-limit sets can be constructed [6,7], and the construction of Subsection 1.2 is partly inspired from these works.

Equicontinuous dynamics
The notions of equicontinuity and equicontinuous points in some directions, introduced in [45], allow to detect some rigidity in a particular direction. In that case, it is easy to obtain informations on the µ-limit measures sets, which are σ-invariant (see Figure 1.2 for some examples). Given a direction α P R, we say that ‚ F is equicontinous if for all e P N there exists d P N such that x r´d,ds " y r´d,ds ùñ F n pσ tnαu pxqq r´e,es " F n pσ tnαu pyqq r´e,es for all n P N, ‚ x P A Z is an equicontinous point of F if for all e P N, there exists d P N such that x r´d,ds " y r´d,ds ùñ F n pσ tnαu pxqq r´e,es " F n pσ tnαu pyqq r´e,es for all n P N. In this case, the word x r´d,ds is said a blocking word.

Nilpotence.
A cellular automaton F is nilpotent if there exists k P N and a P A such that F k pxq " 8 a 8 for all x P A Z . This is equivalent to having two different directions of equicontinuity. In this case, VpF, µq " for all µ P M σ pA Z q.
Equicontinuity. If a cellular automaton F is equicontinuous of direction α, then it is ultimately periodic in this direction. Thus, there exist k, p P N such that Equicontinuous points in two directions. If a cellular automaton F has equicontinuous points in two directions, there exists a P A such that VpF, µq " for all µ P M full σ´erg pA Z q (see Theorem 4.8 of [45]). Equicontinuous points in one direction. If a cellular automaton F has equicontinuous points in one direction, there exists a word u, called a blocking word, such that for any word v, all the configurations in the cylinder ruvus have the same ultimate periodic behavior in direction α. In this case, V 1 pF˚, µq is a singleton when µ P M full σ´erg pA Z q, and we have some information on the limit measure [3,21].
Remark 1.1. In the last two classes, the condition of full support can be weakened by considering an initial measure charging a blocking word. Moreover, the notion of equicontinuous points can be relaxed by the notion of µ-equicontinous points, which takes into account the initial measure [45].

Cellular automata with particles
For some cellular automata, starting from a random configuration, a system of interacting particles emerges. Particles correspond to interfaces between homogenous zone in the space time diagram. Generally, in this case, there is convergence towards a "simple" measure (see Figure 1.

for some examples).
An interesting example is given by the following gliders cellular automaton, defined on the alphabet t0, L, Ru. Under the action of the gliders cellular automaton, the letter L is shifted to the left, the letter R is shifted to the right, and two letters disappear if they intersect. It is shown in [2] that for µ P M σ´erg pA Z q, one has VpF, µq " tνu, with: Another interesting example is the n-states cyclic cellular automaton on the alphabet t0, . . . , n´1u, defined by F n pxq i " x i`1 rmod ns if x i´1 or x i`1 is equal to x i`1 rmod ns, and F n pxq i " x i otherwise. For the 3states cyclic cellular automaton, if µ is a Bernoulli measure, one has VpF 3 , µq " tµpr2sq p δ 0`µ pr0sq p δ 1`µ pr1sq p δ 2 u (see [24] Section 2.3 for a proof). For the 4-states cyclic cellular automaton, one has VpF 4 , µq " tα 0 p δ 0`α1 p δ 1ὰ Gliders CA Rule 184 or traffic rule One-sided captive Another one-sided captive 3-states cyclic CA 4-states cyclic CA 5-states cyclic CA Random walk CA Examples of cellular automata with emergence of particles (the typical asymptotic behavior is described in [26]).

Algebraic cellular automata
When A is a finite group, if the cellular automaton F is a group morphism for the product group A Z , we say that F is linear (the cellular automaton of Example 1.1 was such a linear cellular automaton). In this case, a randomisation phenomenon appears: as soon as the initial measure is in a large class which contains Markov measures, the Cesàro mean sequence converges to the uniform Bernoulli measure [18,36,37,41]. In other words, V 1 pF, µq " tλ A Z u. In some simple cases, it is possible to describe VpF, µq with the help of the convolution measure of µ (see [24], Chapter III).

Possible convergence towards a unique measure
The first natural question is to determine what are the measures that can be reached asymptotically by the iteration of a cellular automaton from an initial measure. Given a measure ν, this measure is clearly reached asymptotically if we iterate the identity cellular automaton from this measure. In fact, the question is not precise enough: we want to determine the measures that are reached asymptotically starting from a large class of measures, for example, starting from all Bernoulli measures. Computability obstruction. As cellular automata can be seen as models of computation, computability obstructions naturally appear. A function f : A˚Ñ B˚is computable if there exists an algorithm which, for a given word u over the alphabet A given as input, returns the word f puq over the alphabet B. This definition can be generalized to N, Q or A˚ˆN, by the introduction of a natural encodage. We distinguish two notions of computable measure, depending on whether we know or not the speed of convergence of the algorithm approximating the measure.
The set of computable (resp. limit-computable) measures is very large, even if it is countable, since the number of programs is countable. Let us present some examples: ‚ any measure supported by a periodic orbit is computable; ‚ any Bernoulli measure or Markov measure with computable (resp. limit-computable) parameters is computable (resp. limit-computable); ‚ if an effective subshift (i.e. a subset of A Z for which the set patterns which does not appear in a configuration can be enumerated by a Turing machine) has a unique σ-ergodic measure µ, then µ is computable. The image of a computable measure by a cellular automaton F of radius r is computable. Indeed, given a word u P A n , one can compute the set of words in A n`2r which are pre-images of u, and for each of these words, one can compute an approximation of the measure of the corresponding cylinder. One thus obtains an approximation of F˚µprusq, with a bound on the error. As a consequence, we obtain the following obstruction for the limit.
Theorem of realization. In fact, there is no other obstruction: given a limit-computable measure ν, it is always possible to "program" a cellular automaton forcing the convergence to the limit ν starting from any σ-ergodic measure of full support. This is given by the following realization theorem. Theorem 1.1 (see [25]). Let ν P M σ pA Z q be a limit-computable measure, there exists a cellular automaton The property for a measure ν P M σ pA Z q to be limit-computable is equivalent to the existence of a sequence pw n q nPN of words, enumerated by an algorithm, such that y δ wn ÝÑ nPN ν. The idea of the proof is thus to construct a cellular automaton enumerating the words pw n q nPN and copying them successively on the space. In order to do so, we extend the initial alphabet with auxiliary states. The dynamics of this cellular automaton, described in [25], is very simple. It consists in three different phases, which coexist: Formatting: Computation takes place on segments, defined as regions of Z located between two symbols W , that play the role of walls, and persist in time except under special circumstances. Since we have no control over the initial content of each segment, we first want to format them, that is, to erase uninitialized content of auxiliary states which can perturb the computation. A segment is said to be initialized if the walls come from a special state I , which can only appear in the initial configuration, and send a signal to the right. This signal deletes all auxiliary states until finding another initialized wall W . The difficult task is to distinguish uninitialized walls from initialized walls. To do that, the right signal sent by I keeps track of its age using a binary counter. Meanwhile, each initialized wall also keeps track of its age under the form of a binary counter to its left, incrementing at each step. Counters already present in the initial configuration (uninitialized) have a nonnegative value at time 0, whereas those created by a symbol I (initialized) have value 0 at time 1, and they increment at the same rate. Thus, uninitialized walls have older time counters, and by comparing time counters and formatting counters as they cross, we can erase older counters and uninitialized walls. This version with binary counter is detailed in [25], there exists also a geometric version of this process of formatting, described in [5,6,15] Computation and copy: A Turing machine is simulated in a small part of the segment. Precisely, if the segment has size k, the Turing machine needs a length opkq of space for the computation. In [25], the space is delimited by the time counter, whereas in [5,6], the size of the segment is measured by a signal, and a small part is used to simulate the Turing machine, whose evolution can be seen as a particular cellular automaton. The machine considered computes successively each w n on the part allowed for the computation. When a word w n is produced, concatenated copies are written progressively on all the segment. By the process of copying, the cellular automaton produces a sample of y δ wn on the segment associated. Merging: In order to have larger and larger domains of computation with a low density of auxiliary states, a segments of size k merges with a neighboring segment, when we need more than opkq cells for auxiliary states. Different processes of merging can be proposed, either unsynchronized as in [6], or synchronized as in [25]. When walls are destroyed in the merging process, only the computation of destroyed walls are stopped. The sample written on the newly merged segment can be only deleted by a copy of a new sample computed with more space. It remains to prove the convergence towards ν. In short, for a given and a word u P A l , there exists K such that | y δ w k prusqq´νprusq| ă for all k ě K. By the merging process, for a time n large enough, with probability larger than 1´ , the pattern F n pxq r0,l´1s is in a sufficiently large segment which allows to have enough space to compute w K , and where the density of auxiliary states is less than (this is possible since auxiliary states occur in a negligible part of the segment compared to its size). Thus, with probability 1´ , the configuration F n pxq r0,l´1s contains auxiliary states, or is a part of a sample composed of successions of w k for some k ě K. We deduce that for sufficiently large n, one has max kěK |F n µprusq´y δ w k prusqq ă 2 . Thus, |F n µprusq´νprusq| ÝÑ Generally, the sequence pF n µq nPN does not converge. So, the second step in the problematic is to determine the set of limit points, denoted by VpF, µq. By compactness of M σ pA Z q, it is a non-empty set. Some necessary conditions can be described using the setting of computable analysis on metric spaces (the formalism is introduced in [48]). In [25], the authors give a partial characterization when the set VpF, µq is connected.
Open question 1.1. What are the sets of measures that can be reached asymptotically iterating a measure by a cellular automaton?

What happens if we do not allow additional letters?
The theorem of realization of Section 1.2 uses additional states to encode the different processes. In particular, we need a special symbol, denoted by I , that only appears in the initial configuration. What can we do if we do not allow additional states in the construction? Given a limit-computable measure ν, we want to construct a cellular automaton defined on A Z , which realizes ν as limit measure. If ν does not charge a word u, it is possible to code auxiliary states using this word and we obtain the following corollary. Corollary 1.1 (see [25]). Let ν P M σ pA Z q be such that there exists a word u which is not charged by ν. Then, there exists a cellular automaton F : If µ is full support, the cellular automaton constructed must be surjective, since every word has pre-images. In fact, there is no hope having a result similar to Corollary 1.1: indeed, the uniform Bernoulli measure is invariant by surjective cellular automata, so it is the only limit measure that can be reached starting from the uniform Bernoulli measure. The following question remains.

Open question 1.2. What are the limit measures that can be reached by surjective cellular automata?
The difficulty is that we cannot initialize the computation by a special pattern that would appear only in the initial configuration. In fact, there are some other obstructions that are not well understood yet. For example, the action of a surjective cellular automaton on random configurations preserves the quantity of information. Formally, this means that the metric entropy of the shift according to µ P M σ pA Z q, denoted by h µ pσq, is equal to the metric entropy according to F˚µ (see [32]). By the semi-continuity of the entropy, the limit measures can only have a greater entropy than the initial measure. The following obstruction follows.

Proposition 1.2. Let F : A Z Ñ A Z be a surjective cellular automaton and let
We remark that this obstruction implies the obstruction given by the uniform Brenoulli measure since it is the unique measure of maximal entropy.

What is the speed of convergence?
In Theorem 1.1, we can estimate the speed of the convergence to ν. The merging process used in [25] is very slow and might be improved for some classes of measures. A natural question is to determine what are the classes of measures that can be quickly reached, and we can wonder if this property is related to the algorithmic complexity of the measures. Let us recall briefly the probabilistic presentation of Variable Length Markov Chains (VLMC), following [10] (cf [44] and [14, pp. 117-134]). Introduce the set L " A´N of left-infinite words on the alphabet A :" td, uu or A :" t0, 1u and consider a complete tree on this alphabet, i.e. a tree such that each node has 0 or 2 children, whose leaves C are (possibly infinite) words on A. The set C is supposed finite or countable. To each leaf c P C, called a context, is attached a probability distribution q c on A. Endowed with this probabilistic structure, such a tree is named a probabilized context tree. The related VLMC is defined as the Markov Chain on L whose transitions are given by PpU n`1 " U n |U n q " q ÐÝ pref pUnq p q, where ÐÝ pref pwq P C is defined as the shortest prefix of w P L, when w is read from right to left, appearing as a leaf of the context tree. In the example of Figure 2.4, on the left side, the context tree is finite of height 4 and, for instance, PpU n`1 " U n d|U n "¨¨¨duduuudq " q duu pdq because ÐÝ pref p¨¨¨duduuudq " duu (read the word¨¨¨duduuud right-to-left and stop when finding a context). The rightmost letter of the sequence U n P L will be denoted by X n so that @n ě 0, U n`1 " U n X n`1 .
The final letter process pX n q ně0 is not Markov as soon as the context tree has at least one infinite context. When the tree is finite, pX n q ně0 is a Markov chain whose order is the height of the tree, i.e. the length of its longest branch. The vocable VLMC is somehow confusing but commonly used.
Consider the probabilized context tree given on the right side of Figure 2.4, called Infinite Comb probabilized context tree. For this model, [13] established a sufficient condition for the existence of a stationary probability measure. In this case, there is one infinite leaf 0 8 and countably many finite leaves 0 n 1, n P N. The data of a corresponding VLMC consists thus in probability measures on A " t0, 1u: q 0 8 and q 0 n 1 , n P N. converges. Assuming the series ř c n converges, the stationary probability measure π on L is unique.

Towards the bestiary
The following propositions give some examples of VLMC admitting a stationary probability measure without necessary and sufficient condition of series convergence. [19]). Let pU n q ně0 be a VLMC. Assume there is a finite number of infinite contexts. Moreover, assume there exists ε ą 0 such that @c P C, @α P A, ε ă q c pαq ă 1´ε.

Proposition 2.3 (Meyn and Tweedies
If 1 P C and for some integer a ě 1, 0 a P C (or 0, 1 a are in C), then the Markov process pU n q ně0 admits a unique invariant measure.

Persistent random walks
Classical random walks are usually defined from a sequence of independent and identically distributed (i.i.d.) increments pX n q ně1 by Assume that @c P C, @α P A, q c pαq ‰ 0, 1. For all these context trees, the Markov process pU n q ně0 admits a unique invariant measure without any condition. For the first and the second tree, the existence and unicity is a consequence of Proposition 2.3. For the last tree, called the bamboo blossom, an explicit formula for the invariant measure can be found in [10]. The existence and unicity can also be viewed as a corollary of Proposition 2.3. For the third tree, a pedestrian and quite long calculation leads to an explicit formula as for the bamboo.
When the jumps are defined as a finite-order Markov chain, a short memory in the dynamics of the stochastic paths is introduced and the random walk pS n q ně0 itself is no longer Markovian. Such a process is called in the literature a Persistent Random Walk (PRW), a Goldstein-Kac random walk or also a correlated random walk. Concerning the genesis of the theory, we allude to [16,23,31,43,49,50] as regards the discrete-time situation but also its connections with the continuous-time telegraph process and we refer to [27,35] concerning recurrence and transience features. Here, we aim at investigating the asymptotic behavior of one-dimensional PRW for which the increments are driven by a VLMC (an infinite-order Markov chain) built from a probabilized context tree. This construction furnishes an extented model for the dependence of the increments of PRWs which can be easily adapted to various situations. Let pU n q ně0 be a VLMC defined on the alphabet A " td, uu. The n th increment X n of the corresponding PRW is given as the latest letter of U n , with the one-to-one correspondance d "´1 (for a descent) and u " 1 (for a rise). X n :"`1 if U n " U n´1 u whereas X n :"´1 if U n " U n´1 d. Different probabilized context trees lead to different probabilistic impacts on the asymptotic behavior of resulting PRWs. Besides, the characterization of the recurrent versus transient behavior, the so-called type problem, is difficult in general. Here, we state exhaustive criteria for PRWs defined from a double-infinite comb context tree introduced in [9]. Roughly speaking, the leaves -coding for the memory -are the words on td, uu » t´1, 1u of the form d n u and u n d. In other words, the probability to invert the current direction depends only on its length. In addition, we derive sufficient conditions for the type of PRWs built from a larger class of context trees obtained by some grafting of the original double-infinite comb model. Foremost, we refer to Figure 2.7 that illustrates our notations and assumptions by a realization of S, for the so called double-infinite comb PRW. In order to avoid trivial cases, we assume that S cannot be frozen in one of the two directions with a positive probability -it makes infinitely many U-turns almost surely. Besides, for the sake of simplicity, we deal throughout this paper with the conditional probability with respect to the event pX 0 , X 1 q " pu, dq -the initial time is suppose to be an up-to-down turn. Obviously, there is no loss of generality supposing this and the long time behavior of S is not affected as well. Therefore, we assume the following.   Let τ u n and τ d n be respectively the length of the n th rise and of the n th descent, that is Then, by the renewal property, pτ d n q and pτ u n q are independent sequences of i.i.d. random variables. In order to deal with a more tractable random walk built with possibly unbounded but i.i.d. increments, we introduce the underlying skeleton random walk pM n q ně0 associated with the even U-turns -the orginal walk observed at the random times of up-to-down turns. Note that the expectation d M of an increment Y n :" τ u n´τ d n is meaningful whenever one of the persistence times (at least) is integrable. In this situation, we can set (extended by continuity whenever one of the persistence times is not integrable) The case when the characteristics J u|d and J d|u are both finite do not appear in the latter theorem since in that situation it follows from [17] that the persistence times are both integrable and hence it reduces to the well-defined drift case in Proposition 2.4.
Consider a double-infinite comb and attach to each finite leaf c another context tree T c (possibly trivial) as in Figure 2.8. The leaves of the related graft are denoted by C c and this one is endowed with Bernoulli distributions tq : P C c u on tu, du. Note that any probabilized context tree on tu, du can be constructed this way. We denote by S g the corresponding PRW. In that case, the random walk is particularly persistent, in the sense that the rises and descents are no longer independent. A renewal property may persist but it is heavier As a consequence, if S and S are of the same recurrent or transient type, then S is recurrent or transient accordingly. We can extend our results to some probabilized context trees satisfying Assumption 2.1, mainly when the persistence times of S or S are both non-integrable.

Perfect sampling
Perfect sampling was introduced in 1996 by Propp and Wilson [42]. This technique allows one to generate unbiased samples from the stationary distribution of an ergodic Markov chain over a finite state space. The key feature is the use of coupled Markov chains originating from all possible states.

Markov automata
For the purpose of simulations, a Markov chain pX n q nPN with transition matrix P is represented by means of a Markov automaton A (not to be confused with the cellular automata discussed in Section 1). In addition to the state space S, A comprises a finite event space A, a distribution ρ and a transition function¨: SˆA Ñ S. The Markov chain pX n q nPN is obtained by generating an i.i.d. sequence pA n q nPN according to ρ, and recursively setting X n`1 " X n¨An`1 . The sequence pX n q nPN is clearly a Markov chain, and if the Markov automaton is chosen such that Ppi¨A 1 " jq " P i,j for all i, j P S, then its transition matrix is P .

Coupling From The Past
The algorithm described here is based on the following intuition: a Markov chain that has been running "infinitely long" is now distributed according to its stationary distribution π. Starting this Markov chain arbitrarily far back in time and determining its state at time 0 gives a sample distributed according to π. The proof underlying this method can be found in [42].
The objective is to determine this final state using only a finite amount of information. The key idea is to consider only the last k 1 steps of the Markov chain, generated by the sequence of events A´p k1´1q , . . . , A 0 " ρ bk1 . As we have no means of determining the value of X´k 1 , we take into account all possibilities, i.e. we compute X 0 for all possible starting states at time´k 1 . This gives us a set of plausible end states If S 0 is a singleton, then regardless of X´k 1 , X 0 is the unique element of S 0 , which is therefore distributed according to π. Otherwise, it simply means we have not considered a large enough portion of the past; we therefore take a k 2 ą k 1 , generate new events A´p k2´1q , . . . , A´k 1 " ρ bk2´k1 independently from the previous ones, and compute S 2 " ts¨A´p k2´1q¨. . .¨A 0 , s P Su. By reiterating for increasing values of k i until S i is a singleton, we have an algorithm which returns a state distributed according to π. The exact choice of the pk i q iPN does not change the distribution of the output, but can impact complexity. When taking k i " i deterministically, for example, one can maintain a unique vector s in memory, updating it at each iteration by setting s i`1 pxq " s i px¨A´iq. We then have that S i " ts i pxq|x P Su, i.e. the algorithm terminates when all the entries in s are equal. In the case of monotonicity (described below), it is no longer possible to work with this unique vector: the entire trajectories must be computed once more. Doubling the value of k between successive iterations, however, greatly reduces the execution time of the algorithm [42].
Furthermore, call τ the first k i such that S i is a singleton; regardless of the choice of pk i q iPN , so long as one can construct a finite sequence A 1 , . . . , A l with positive probability such that ts¨A 1¨. . .¨A l , s P Su is a singleton, then τ has finite expectation, i.e. the algorithm 1 terminates in finite expected time. For bounds on this running time, see [42].

Reducing complexity through envelopes
The above algorithm can be refined in many ways, such as that described in [51], which does not require storing the pA´nq nPN . In general, sampling techniques based on Markov chains, often referred to as Monte Carlo methods, are used when sampling from a large, multi-dimensional set S for which computing the exact distribution is not an option. Coupling From The Past is interesting in that it generates unbiased samples, but in its raw form, its complexity is linear in the size of S. Techniques have therefore been introduced to overcome this limitation, the most common of which use some form of envelope.
Informally, an envelope E implicitly or explicitly defines a subset of S, comprising the states it encompasses. The idea is to define a low-dimension set of envelopes which allow us to keep track of the evolution of our Markov chains at a reduced cost: the transition function¨is generalized in order to apply to these envelopes, after which we run the simulation over the envelopes (S i becomes S 1 i " S 1 0¨A´pki´1q¨. . .¨A 0 , where S 1 0 encompasses all of S). The algorithm now terminates when S 1 i comprises only one state. Many envelope-based strategies exist [22,30,33,40]. The most common examples rely on an ordering of the set space; this is the case for both the monotonous [42] and non-monotonous [8] examples described below.

Monotonous chains
Suppose that one can define a partial order ĺ over S such that there exists an upper bound s max and a lower bound s min for ĺ, and such that the following holds: @i, j P S, @A P A, s i ĺ s j ùñ s i¨A ĺ s j¨A .
Such a chain is said to be monotonous for ĺ.
Consider the set of intervals tpm, M q, m ĺ M u, where rm, M s " ts P S | m ĺ s ĺ M u. Taking S 1 0 " ps min , s max q and pm, M q¨A " pm¨A, M¨Aq, one can verify that, for all i P N, S 1 i contains all the states in S i [42]. A direct consequence of this is that, if ever S 1 i " pm i , M i q is such that m i " M i , then S i is a singleton, and its unique element, m i , is distributed according to π.

Non-monotonous chains
Suppose now that the state space is still partially ordered, but that the dynamic is no longer monotonous. We still use envelopes of the form pm, M q, with m 0 " s min and M 0 " s max , but the generalization of the transition function is altered: Like before, if this dynamic over the envelopes leads to a point where m i " M i , then m i is distributed according to the stationary distribution π [8].
Notice that with this technique, pm i q and pM i q are not trajectories of the initial dynamic. Furthermore, one must verify that the algorithm converges; the fact that there exists a sequence A 1 , . . . , A l such that ts¨A 1¨. . .Ä l , s P Su is a singleton does not imply that there exists a sequence A 1 1 , . . . , A 1 k such that ps min , s max q¨A 1 1¨. . .¨A 1 k " ps, sq for some s P S.

Closed queuing networks model
Consider a closed queuing network of K queues and M customers in total. Each queue k has a finite capacity C k ď M and a single server with independent exponential services of rate µ k . After a customer has been served in queue i, it is routed to queue j with probability p i,j , independently of the current state and past evolution of the network. If queue j is full, then the customer remains in queue i and waits to be served again (a new destination will then be drawn, independently from the previous one). A configuration is represented by a vector x P S where x i represents the number of customers in queue i. The state space S being defined as If all queues have an infinite capacity (i.e. if for each k P t1, . . . , Ku, C k " M ) then |S| "`K`M´1 K´1˘. Thus for M " K the complexity of the state space is OpM K q.
The dynamics of the system defines a Markov chain pX n q nPN , represented by the Markov automaton composed of the event space tt i,j |1 ď i, j ď Ku, the distribution ρ such that ρpt i,j q " µipi,j ř K k"1 µ k , and the transition function t i,j which describe routing from queue i to queue j is defined by The aim is to sample the stationary distribution with perfect sampling techniques. The chain is nonmonotonous and the cardinality of the state space is OpM K q, however Bouillard et al. [4] developed a new envelope-based strategy to improve the perfect sampling algorithm for this model. The idea is to exploit the constraint of the total number of customers in the network to represent any set of states as a directed graph, called diagram.

Set of states as a diagram
In a diagram, a state is encoded by a path of length K starting at node p0, 0q and ending at node pK, M q, such that the arc ppk´1, mq, pk, m 1 qq belongs to the path if and only if there are m 1´m customers in queue k. A given diagram envelopes all configurations for which the encoding is present. A diagram is said to be complete if it encodes the state space S. A complete diagram contains at most OpKM 2 q arcs. In Figures 3.12 and 3.13, we illustrate the encoding of the configuration p0, 0, 2, 0, 1q, as well as the complete diagram for a network with K " 5 queues, M " 3 customers and a capacity vector C " p2, 1, 3, 1, 2q.

Transition in a diagram
The transition induced by event t i,j is applied directly on the diagram D. To compute this transition, we need to know what paths of D will remain unchanged, either because the number of customers in queue i is equal to 0 (set of arcs Stay), or because queue j is full (set of arcs Full). From there, we deduce which arcs will be modified (set of arcs T ransit). In accordance with t i,j , these arcs are modified as follows: ‚ In column i, the slopes of the arcs are decreased by one. ‚ In columns k P tminpi, jq`1, . . . , maxpi, jq´1u, the arcs are shifted up or down, depending whether i ă j or i ą j. ‚ In column j, the slopes of the arcs are increased by one. Figure 3.14 shows how the diagram D is modified by event t 4,2 . First, the sets Stay, Full and T ransit are determined according the values of arcs in columns 4 and 2. As an example, arcs in Stay are along paths such that arcs in column 4 represent 0 customers. Afterward, arcs in T ransit are updated in accordance with t 4,2 , that corresponds to the set T ransit 1 . Finally, T 4,2 pDq is the diagram that contains arcs in StayYFullYT ransit 1 .
Each transition has complexity OpKM 2 q, that corresponds to the number of arcs in a diagram.

Perfect simulation
To use diagrams with Coupling From the Past algorithm (Algorithm 1), we start by replacing the set S with the complete diagram, and apply the transitions t i,j on this diagram. The algorithm terminates when the final diagram contains only one path. This path corresponds to a state. Bouillard et al. [4] proved that there exists a finite sequence of events pt i1,j1 , . . . , t i N ,j N q that reduces the complete diagram into a diagram that contains only one path, thereby proving the following theorem: Theorem 3.1. The coupling from the past algorithm using diagrams terminates in finite expected time and produces an exact sample from the stationary distribution.

The hard-core gas model
For a simple graph G " pV, Eq, define an independent set as a subset I of V , such that no two vertices of I are connected by an edge in E. Denote I the set of independent sets.
The hard-core gas model, originating from statistical physics [20], looks at the feasible configurations for a gas. In this model, particles occupy sites, represented by vertices, and incompatible sites are connected by an edge. A configuration is considered feasible if no two incompatible sites are occupied, i.e. if it forms an independent set. The likelihood of each configuration is of the form The positive parameter λ is called the fugacity of the system, and Z λ is a normalizing constant. We are interested in generating random samples according to the distribution π. This cannot be done in a straightforward fashion, as that would require computing Z λ , which is 7P -complete a . We therefore rely of Monte Carlo methods, and namely on Coupling From the Past.
A simple dynamic for defining a Markov chain with stationary distribution π is the following: let A " Vˆt`,´u be the set of events, and ρ be the distribution over A such that ρpu`q " 1 |V | λ 1`λ and ρpu´q " 1 |V | 1 1`λ (u is a vertex chosen uniformly at random, after which we pick`or´independently with respective probabilities # I Y tuu if I Y tuu P I, I otherwise. The resulting Markov automaton defines a Markov chain whose stationary distribution is π. An enhanced variant of this dynamic can be found in [30].
For coupling from the past, we build an envelope based upon the partial ordering induced by inclusion [30]: we maintain two sets m, M Ď V such that, regardless of the initial state of a trajectory X, we have m i Ď X i Ď M i at each time step i. Note that m and M need not be independent sets themselves; indeed, to envelope trajectories emanating from all states, it is necessary to take m 0 " ∅ and M 0 " V . The transitions are then applied to the envelope as followed: pm, M q¨u´" pmztuu, M¨ztuuq and pm, M q¨a`" if m Y tuu P I, pm, M q otherwise. This envelope can be used to run the Coupling From the Past algorithm for the hard-core gas model. Bounds on its complexity can be found in [30].

Event skipping
The overall Monte Carlo scheme for the hard-core gas model suffers from an inherent weakness: it can easily be verified that the probability of having I¨a " I for a " ρ is at least 1 2 , and can be arbitrarily close to 1 for eccentric values of λ (either very large or very small with respect to 1). For typical Monte Carlo, this is not necessarily an issue, as one can sample a conditioned to having I¨a ‰ I. However, one must pay special attention to the time at which the simulation is stopped, in order to make sure that the output is not biased. Since the conditional distribution of a could be obtained through rejection sampling (though this method would not reduce the complexity, since "skipped" events are still sampled), a simple coupling argument suffices to show that the number of omitted transitions one should account for at each step of the algorithm is geometrically distributed. Counting these transitions allows one to stop the algorithm in such a way as to obtain an unbiased output.
For Coupling From the Past, this strategy can be adapted, but requires even more care. Indeed, this simple conditioning cannot be applied, as a given event A´i does not always apply to the same state: it is first applied a 7P is the complexity class of counting problems associated to decision problems in NP.
to pm 0 , M 0 q¨A´p k0´1q¨. . .¨A´p i`1q , but then as we start the simulation from further back, it is also applied to pm 0 , M 0 q¨A´p k1´1q¨. . .¨A´p i`1q , and so on (cf. Figure 3.15). The sequence of events cannot simply be re-sampled independently from the preceding ones every time we re-run the simulation from further back in time; it is therefore necessary to introduce an update scheme for the sequence of events. This scheme must both remove events which do not modify the "current" state of the envelope, and insert new events according to a precise distribution, so as to compensate any bias that might result from removals in previous states. Note that we cannot simply keep track of which events have been removed: this presents the same caveat as rejection sampling, where the overall complexity is unchanged, since all the events are still sampled. The exact details of this scheme, as well as the proof of its correctness, can be found in [47].