PERFECT SAMPLING FROM THE LIMIT OF DETERMINISTIC PRODUCTS OF STOCHASTIC MATRICES

We illustrate how a technique from the theory of random iterations of functions can be used within the theory of products of matrices. Using this technique we give a simple proof of a basic theorem about the asymptotic behavior of (deterministic) ``backwards products'' of row-stochastic matrices and present an algorithm for perfect sampling from the limiting common row-vector (interpreted as a probability-distribution).


Preliminaries
An N × N matrix P = (p i j ) is called stochastic if j p i j = 1, and p i j ≥ 0 for all i and j.A stochastic matrix P is called stochastic-indecomposable-aperiodic (SIA) if exists and has all rows equal.Such matrices are of interest in the theory of Markov chains since if {X n } is a Markov chain with a transition matrix, P, which is SIA then in the long run X n will be distributed according to the common row vector of Q independent of the value of X 0 .Basic theory of Markov chains gives soft conditions ensuring a stochastic matrix to be SIA.Essentially "periodicity" has to be ruled out, and one fixed closed irreducible set should eventually be reached by the chain from any given starting point.To have something in mind, the matrix      0.5 0.5 0 0 0.5 0.5 0 0 0 0.5 0 0.5 0 0 0.5 0.5 are not SIA since decomposability and periodicity respectively are present in the latter two.
Wolfowitz [11] proved that if A 1 , ..., A n 0 are N × N matrices such that any finite product of them is SIA, then for any ε > 0 there is an integer ν(ε) such that m ≥ ν(ε) implies that every product of m A i 's has rows which differ by at most ε.In particular, if ω = ω 1 ω 2 . .., with ω j ∈ {1, ..., n 0 } for any j, then the limit exists with all rows of L ω being identical.Such products are known as infinite "backwards products" in the literature, see e.g.Seneta [8].In contrast, "forwards products", defined by multiplying the matrices in reversed order will typically not converge.These notions have their roots in the theory of time-inhomogeneous Markov chains where "forwards/backwards products" corresponds to running the Markov chain forwards/backwards in "time".It is not sufficient to require merely that A 1 , ..., A n 0 are SIA in Wolfowitz' theorem since it is e.g.easy to construct two SIA matrices A 1 and A 2 such that their product A 1 A 2 is decomposable and (A 1 A 2 ) n will therefore not have approximately identical rows for any n, see e.g.Hajnal [4], Equation 4(b).The conditions of Wolfowitz' theorem may therefore be hard to verify.Some rather weak sufficient conditions for Wolfowitz' condition however follow from Sarymsakov [7].In addition, Thomasian [9] proved a result which implies an algorithm, requiring a bounded number of arithmetical operations, for determining whether any finite product (of whatever length) of SIA-matrices is SIA.See Paz [6] for a sharp bound.
In the present paper we will extend Wolfowitz' theorem with an algorithm for simulating from the common row vector of L ω .In the particular case when n 0 = 1 our algorithm reduces to a "Coupling From The Past (CFTP)"-algorithm for "perfect sampling" from the stationary distribution of a homogeneous Markov chain with a SIA transition matrix.
Let us present the simple key idea in the simplest "clearest" case when n 0 = 1; Any probability matrix, P = (p i j ), can be represented by an Iterated Function System (IFS) with probabilities, i.e. there exists a finite family of functions f k : {1, ..., N } → {1, ..., N }, and associated probabilities p k , k = 1, 2, .., M , for some constant M , (p k > 0 and Intuitively we may thus "run" a Markov chain with transition probability matrix P by random iterations of the functions in a representing IFS choosing functions f k with the associated probabilities p k .The concept of IFSs was coined by Barnsley and Demko [3] and has mainly been used on general state-spaces as a tool to construct fractals.If P has at least one column with only positive elements then we can choose the IFS with the property that (at least) one of the maps f k , say f k 0 is a constant map (i.e.f k 0 (x) = c, for some constant c), see Athreya and Stenflo [2] for details in a more general setup.If {I n } ∞ n=1 is a sequence of independent random variables with P(I n = k) = p k for any n and k then it follows that i j denotes the element on row i and column j in the matrix P n .If the IFS contains a constant map, then it follows that the random variable of composed maps will not depend on i, (i.e. will be constant) for all n ≥ T , where T = inf{ j : I j = k 0 } is the first time when the constant map is applied (a geometrically distributed random variable).Thus is a random variable having the same distribution as the common row vector of lim n→∞ P n .See Athreya and Stenflo [2] for further details of the above method and Wilson [10] for further literature related to CFTP.

Results
Let Ω be an arbitrary set, and let for each ω ∈ Ω, be a row-stochastic matrix, i.e. a matrix with j p ω i j = 1, and p ω i j ≥ 0 for all i and j.Our basic theorem is the following: Theorem 1. Suppose all matrices P ω satisfy a Doeblin condition uniformly in ω i.e. suppose there exists a constant c > 0 such that for any ω ∈ Ω, where c ω jmin = min{p ω 1 j , p ω 2 j , . . ., p ω N j } denotes the minimum value of the elements in the j:th column of P ω .Let ω = ω 1 ω 2 . .., be an arbitrary sequence of elements of Ω i.e. ω j ∈ Ω, j ≥ 1. a) Convergence towards the limit with an exponential rate: exists, and M ω is a matrix with identical rows, i.e.

b) Perfect sampling:
Let U 1 , U 2 , U 3 ,... be independent random variables uniformly distributed on the unit interval, and let T be a geometrically distributed random variable with and g The result in b) suggest an obvious algorithm for generating perfect samples from the common row-distribution of M ω by simulating sets of a geometrically distributed size of uniformly distributed random numbers.See Athreya and Stenflo [2] for generalizations of the above theorem in the case when Ω contains one point.
Proof The matrices P ω may be expressed as where P ω ⋆ = (p ⋆ω i j ), with p ⋆ω i j = c ω jmin /( N j=1 c ω jmin ) and P ω ⋄ = (p ⋄ω i j ), with for each given ω ∈ Ω, and s ∈ (0, 1).The function g ⋄ω is an IFS representation of P ω ⋄ in the sense that if U is a uniformly distributed random variable on the unit interval, then P(g ⋄ω U (i) = j) = p ⋄ω i j .Similarly the matrix P ω ⋆ can be represented by an IFS containing only the constant maps If for s = (s 1 , s 2 ) ∈ (0, 1) 2 , where χ denotes the indicator function, then if U = (U, Û) is a uniformly distributed random variable on (0, 1) 2 , then and we have thus created an explicit IFS representation of P ω .
If {U n = (U n , Ûn )} is a sequence of independent uniformly distributed random variables on (0, 1) 2 , and we define then if T = inf{n : Ûn > 1−c}, then it follows that T is a geometrically distributed random variable with parameter c, independent of the random variables U n , and for all n, ω, i and j.In order to see this we use induction; First note by construction that P(X ω Therefore (8) follows by induction.Define µ ω j = P(X ω = j).By construction This completes the proof of the theorem.

Remark 2. Let d denote the Dobrushin ergodic coefficient defined for a stochastic matrix A by
for any two stochastic matrices A and B. By using the right hand side of (9) we see that assumption (4) implies that d(P ω ) ≤ 1 − c, for any ω and thus Now fix j and ω and note that min i M ω n (i, j) is increasing in n.Also note, similarly, that max i M ω n (i, j) is decreasing in n.By (10), and therefore there exists a limit µ ω j such that We can thus prove that the exponential convergence rate of Theorem 1a) holds without involving IFSs.A good reason to involve IFSs in the proof above is however that it simultaneously allows us to prove b).
Remark 3. Note that we have no restrictions on the number of matrices in our theorem above.The theorem corresponds to backwards limits for Markov chains in random environments or timeinhomogeneous Markov chains with state space S = {1, ..., N }.It is straightforward to extend the theorem to the case when S is an arbitrary complete separable metric space.Transition matrices P ω = (p ω i j ) are then replaced by transition kernels P ω (x, B) (interpreted as the transition probability of moving from state x into the set B) satisfying the Doeblin condition P ω (x, B) ≥ cν ω (B), for all sets B for some probability measure ν ω depending on ω ∈ Ω.In the proof we then use the fact that any Markov chain on a complete separable metric space can be represented by an IFS, see Kifer [5] or Athreya and Stenflo [2].
As a consequence of Theorem 1 we can now, as our main result, extend the theorem by Wolfowitz [11] with a perfect sampling algorithm; Suppose A 1 , ..., A n 0 are N × N matrices such that any finite product of them is SIA.If ω = ω 1 ω 2 . .., with ω j ∈ {1, ..., n 0 } for any j, then it follows from [11] that for any c ∈ (0, 1) there exists an integer n c such that the matrices N j } denotes the minimum value of the elements in the j:th column of P (k) , k ≥ 1.By applying Theorem 1 we obtain; Corollary 1. Suppose A 1 , ..., A n 0 are N × N matrices such that any finite product of them is SIA.Then for any ω = ω 1 ω 2 . .., with ω j ∈ {1, ..., n 0 } for any j, the limit exists with all rows of L ω being identical and the convergence rates of the row distributions towards the limit is exponential (in total-variation norm) uniformly in the row-index and ω.Remark 4. The exponential convergence rate in Corollary 1 is a known consequence of the theorem in [11], see e.g.Anthonisse and Tijms [1].See also Seneta [8].
An advantage with our method based on IFSs is that it gives the perfect sampling algorithm "for free".