Utilizing Network Structure to Accelerate Markov Chain Monte Carlo Algorithms

We consider the problem of estimating the measure of subsets in very large networks. A prime tool for this purpose is the Markov Chain Monte Carlo (MCMC) algorithm. This algorithm, while extremely useful in many cases, still often suffers from the drawback of very slow convergence. We show that in a special, but important case, it is possible to obtain significantly better bounds on the convergence rate. This special case is when the huge state space can be aggregated into a smaller number of clusters, in which the states behave approximately the same way (but their behavior still may not be identical). A Markov chain with this structure is called quasi-lumpable. This property allows the aggregation of states (nodes) into clusters. Our main contribution is a rigorously proved bound on the rate at which the aggregated state distribution approaches its limit in quasi-lumpable Markov chains. We also demonstrate numerically that in certain cases this can indeed lead to a significantly accelerated way of estimating the measure of subsets. The result can be a useful tool in the analysis of complex networks, whenever they have a clustering that aggregates nodes with similar (but not necessarily identical) behavior.


Introduction
The Markov Chain Monte Carlo (MCMC) method is one of the most frequently used algorithms to solve hard counting, sampling and optimization problems.This is relevant for many areas, including complex networks, physics, communication systems, computational biology, optimization, data mining, big data analysis, forecast problems, prediction tasks, and innumerable others.The success and influence of the method is shown by the fact that it has been selected as one of the top 10 of all algorithms in the 20th century, see [1].
The MCMC algorithm also plays an important role in large, complex networks.In this paper, building on our earlier conference presentations [2,3], we consider the following regularly occurring application of the MCMC method: Consider a very large graph G, with node set S, and let A ⊆ S be a subset of the nodes.We would like to estimate the relative size of A, that is, the goal is to obtain a good estimate of the value More generally, if a random walk is considered on the graph, with stationary distribution π, then we would like to estimate π(A), the stationary probability of being in A. In the special case when π is the uniform distribution, we get back the formula (1).
If we can take random samples from S, according to the stationary distribution, then an obvious estimate with good properties is the relative frequency of the event that the sample falls in A.
Unfortunately, in most nontrivial cases of interest, this sampling task is not feasible.The reason is that often the large set S is defined implicitly.Examples are the set of all cliques in a graph, or the set of all feasible solutions to an optimization problem, and many others.No efficient general method is known to sample uniformly at random from such complex sets.
An important application in telecommunication networks is to estimate blocking probabilities, see [4,5].More generally, if we have a large system, with an enormous state space, we may want to estimate that the actual state falls in a specific subset.For example, if the state space consists of all possible load values of the network links, which leads to a state space of astronomical size, we may want to know what the probability is that at most k links are overloaded, for some value of k.
At this point, the MCMC does a very good service.If we define a Markov chain in which the states are the elements of S and the transitions are based on simple local operations, then we can very often obtain a Markov chain with uniform, or some other simple stationary distribution over S.Then, if we run this chain long enough so that it gets close to the stationary distribution, then the state where we stop the chain will be a good approximation of a random sample over S, distributed according to the stationary distribution.Then by repeating the process sufficiently many times, and by counting the relative frequency that the random sample falls in A, we can get a good estimate of the probability measure of A.
The key difficulty is, however, that we should run the chain long enough to get sufficiently close to the stationary distribution.This time is often referred to as mixing time [6].If the mixing time grows only polynomially with the size of the problem, e.g. with the size of the graph, then we say that the chain is rapidly mixing.Unfortunately, in many cases of interest the mixing time grows exponentially with the problem parameters, so in many important cases the Markov chain is mixing very slowly.
What we are interested in is whether it is possible to speed up the running time.It is clear that if we want to estimate the size of any possible subset, then we really need to get close to the stationary distribution, since only this distribution can guarantee that the probability of the random state falling in the set is really the relative size of the set.On the other hand, if we only want to estimate the relative size of a specific subset A, then it is enough for us if we reach a distribution in which the measure of A is close to the stationary measure, but this does not have to hold for every other set.In other words, if π t denotes the state distribution after t steps and π is the stationary distribution, then we want to choose t such that |π t (A) − π(A)| is small, but the same does not have to hold for all other sets.This makes it possible to reduce the required value of t, that is, to speed up the algorithm.In this paper we investigate under what conditions it is possible to obtain such a speed-up.
The main result is that the structure of the chain, that is, the network structure, can significantly help, if it has some special properties.Specifically, if the Markov chain is close to a so called lumpable chain, then remarkable speedup is possible.In other words, in this case we can indeed capitalize on the particular network structure to accelerate the method.Below we informally explain what the concept of lumpability means.The formal definition will follow in the next section.
The concept of lumpability stems from the following observation: it is very useful if the state space can be partitioned such that the states belonging to the same partition class "behave the same way," in the sense defined formally in the next section.This is the concept of lumpability [7].Informally speaking, it means that some sets of states can be lumped together (aggregated) and replaced by a single state, thus obtaining a Markov chain which has a much smaller state space, but its essential behavior is the same as the original.
In some cases, the lumpability of the Markov chain can have a very significant effect on the efficiency of the model.A practical example is discussed in [8,9], where the authors present a fast algorithm to compute the PageRank vector, which is an important part of search engine algorithms in the World Wide Web.The PageRank vector can be interpreted as the stationary distribution of a Markov chain.This chain has a huge state space, yielding excessive computation times.This Markov chain, however, is lumpable.Making use of the lumpability, the computation time can be reduced to as low as 20% of the original, according to the experiments presented in [8].
Unfortunately, it happens relatively rarely that the Markov chain satisfies the definition of lumpability exactly.This motivates the concept of quasi-lumpability [10,11].Informally, a Markov chain is quasi-lumpable if its transition matrix is obtainable by a small perturbation from a matrix that exactly satisfies the lumpability condition (see the formal definition in the next section).
In this paper we are interested in the following problem, which is often encountered in applications: how long do we have to run the Markov chain if we want to get close to the stationary distribution within a prescribed error?While the general question is widely discussed in the literature (see, e.g., [6,12]), we focus here on a less researched special case: how much gain can the convergence speed enjoy, if we can capitalize on the special structure of quasi-lumpability.

Aggregation in Markov Chains
We assume the reader is familiar with the basic concepts of Markov chains.We adopt the notation that a Markov chain M is given by a set S of states and by a transition probability matrix P, so we write M = (S, P).This notation does not include the initial distribution, because it is assumed arbitrary.
Let us first define the concept lumpability of a Markov chain.Informally, as mentioned in the Introduction, a chain is lumpable if its states can be aggregated into larger subsets of S, such that the aggregated (lumped) chain remains a Markov chain with respect to the set-transition probabilities (i.e., it preserves the property that the future depends on the past only through the present).Note that this is generally not preserved by any partition of the state space.Let us introduce now the formal definition.

Definition 1. (Lumpability of Markov chain)
holds for any t, k, j, i 1 , . . ., i k , whenever these conditional probabilities are defined (i.e., the conditions occur with positive probability).If the chain is lumpable, then the state set of the lumped chain is Q and its state transition probabilities are defined by Checking whether a Markov chain is lumbable would be hard to do directly from the definition.That is why it is useful to have the following characterization, which is fundamental result on the lumpability of Markov chains, see [7].For simple description, we use the notation p(x, A) to denote the probability that the chain moves into a set A ⊆ S in the next step, given that currently it is in the state x ∈ S. Note that x itself may or may not be in A.

Theorem 1. (Necessary and sufficient condition for lumpability, see [7])
A Markov chain M = (S, P) is lumpable with respect to a partition Q = {A 1 , . . ., A m } of S if and only if for any i, j the value of p(x, A j ) is the same for every x ∈ A i .These common values define the transition probabilities p(A i , A j ) for the lumped chain, which is a Markov chain with state set Q and state transition probabilities where x is any state in A i .
Informally, the condition means that a move from a set A i ∈ Q to another set A j ∈ Q happens with probability p(x, A j ), no matter which x ∈ A i is chosen.That is, any x ∈ A i has the property that the probability of moving from this x to the set A j in the next step is the same for every x ∈ A i .The sets A i , A j are partition classes of Q.We also allow i = j, so they may coincide.
Whenever our Markov chain is lumpable, we can reduce the number of states by the above aggregation, and it is usually advantageous for faster convergence (a specific bound will be proven in Section 3).
It is worth noting that lumpability is a rather special property, and one has to be quite lucky to encounter a real-life Markov chain that actually has this property.Sometimes it happens (see, e.g., the example in the Introduction about PageRank computation), but it is not very common.Therefore, let us now relax the concept of lumpability to broaden the family of the considered Markov chains.The extended condition, as explained below, is called quasi-lumbability.
Informally, a Markov chain is called quasi-lumpable or -quasi-lumpable or simply -lumpable, if it may not be perfectly lumpable, but it is "not too far" from that.This " -closeness" is defined in [10,11] in a way that the transition matrix can be decomposed as P = P − + P .Here P − is a component-wise non-negative lower bound for the original transition matrix P, such that P − satisfies the necessary and sufficient condition of Theorem 1.The other matrix, P , represents a perturbation.It is an arbitrary non-negative matrix in which each entry is bounded by .This definition is not very easy to visualize, therefore, we use the following simpler but equivalent definition.
Note that if we take = 0, then we get back the ordinary concept of lumpability.Thus, quasi-lumpability is indeed a relaxation of the original concept.It can also be interpreted in the following way.If > 0, then the original definition of lumpability may not hold.This means, the aggregated process may not remain Markov.i.e., it does not satisfy (2).On the other hand, if is small, then the aggregated process will be, in a sense, "close" to being Markov, that is, to satisfying (2).
What we are interested in is the convergence analysis of quasi-lumpable Markov chains, typically for a small value of .For the analysis we need to introduce another definition.
Note that it always holds (component-wise) that L ≤ U.If the chain is lumpable, then these matrices coincide, so then L = U = P, where P is the transition matrix of the lumped chain.If the chain is -lumpable, then L and U differ at most by in each entry.
Generally, L and U are not necessarily stochastic matrices (A vector is called stochastic if each coordinate is non-negative and their sum is 1.A matrix is called stochastic if each row vector of it is stochastic.), as their rows may not sum up to 1.

Convergence Analysis
An important concept in Markov chain convergence analysis is the ergodic coefficient, see, e.g., [12].It is also called coefficient of ergodicity.

Definition 4. (Ergodic coefficient)
Let P = [p ij ] be an n × n matrix.Its ergodic coefficient is defined as The ergodic coefficient is essentially the largest L 1 distance that occurs between different row vectors of the matrix P. That is, in a sense, it captures how diverse are the row vectors of the matrix.The 1/2 factor is only for normalization purposes.For stochastic matrices two important properties of the ergodic coefficient are the following [12]: The importance of the ergodic coefficient lies in its relationship to the convergence rate of the Markov chain.It is well known that the convergence rate is determined by the second largest eigenvalue of the transition matrix (that is, the eigenvalue which has the largest absolute value less than 1), see, e.g., [6].If this eigenvalue is denoted by λ 1 , then the convergence to the stationary distribution happens at a rate of O(λ t 1 ), where t is the number of steps, see [12].It is also known [12] that the ergodic coefficient is always an upper bound on this eigenvalue, it satisfies λ 1 ≤ ρ(P) ≤ 1.Therefore, the distance to the stationary distribution is also bounded by O(ρ(P) t ).Thus, the smaller is the ergodic coefficient, the faster convergence we can expect.Of course it only provides any useful bound if ρ(P) < 1.If ρ(P) = 1 happens to be the case, then it does not directly provide a useful bound on the convergence rate, since then ρ(P) t remains 1.In this situation a possible way out is considering the k-step transition matrix P k for some constant integer k.If k is large enough, then we can certainly achieve ρ(P k ) < 1, since it is known [12] that lim k→∞ ρ(P k ) = 0. Now we are ready to present our main result, which is a bound on how fast will an -lumpable Markov chain converge to its stationary distribution on the sets that are in the partition, which is used in defining the -lumpability of the chain.We are going to discuss the applicability of the result in the next section.Theorem 2. Let ≥ 0 and M = (S, P) be an irreducible, aperiodic Markov chain with stationary distribution π.Assume the chain is -lumpable with respect to a partition Q = {A 1 , . . ., A m } of S. Let ρ be any upper bound on the ergodic coefficient of the lower transition matrix L (Definition 3), that is, ρ(L) ≤ ρ.Let π 0 be any initial probability distribution on S, such that P(X t ∈ A i ) > 0 holds for any i, and t = 0, 1, 2, . .., whenever the chain starts from π 0 .Then for every t ≥ 1 the following estimation holds: Remark: Recall that the parameter quantifies how much the Markov chain deviates from the ideal lumpable case, see Definition 2. In the extreme case, when = 1, every Markov chain satisfies the definition.This places an"upward pressure" on : the larger it is, the broader is the class of Markov chains to which -lumpability applies.On the other hand, a downward pressure is put on by Theorem 2: the convergence bound is only meaningful, if ρ + m/2 < 1 holds.This inequality can be checked for any particular , since it is assumed that ρ and m are known parameters.Furthermore, the smaller is , the faster is the convergence.Therefore, the best value of is the smallest value which still satisfies Definition 2 for the considered state partition.
Proof of Theorem 2. Let π 0 be an initial state distribution of the Markov chain M, let π t be the corresponding distribution after t steps and π = lim t→∞ π t be the (unique) stationary distribution of M. For a set A ⊆ S of states the usual notations π t (A) = P(X t ∈ A), π(A) = lim t→∞ π t (A) are adopted.
Using the sets A 1 , . . ., A m of the partition Q, let us define the stochastic vectors for t = 0, 1, 2, . . .and the m × m stochastic matrices for t = 1, 2, . ... Let us call them aggregated state distribution vectors and aggregated transition matrices, respectively.Note that although the entries in (4) involve only events of the form {X t ∈ A k }, they may also depend on the detailed state distribution within these sets, which is in turn determined by the initial distribution π 0 .In other words, if two different initial distributions give rise to the same probabilities for the events {X t ∈ A k } for some t, they may still result in different conditional probabilities of the form P(X t+1 ∈ A j | X t ∈ A i ), since the chain is not assumed lumpable in the ordinary sense.This is why the notations Pt (π 0 ), p (π 0 ) t (i, j) are used.Also note that the conditional probabilities are well defined for any initial distribution allowed by the assumptions of the lemma, since then P(X t ∈ A i ) > 0.
For any fixed t the events {X t ∈ A i }, i = 1, . . ., m, are mutually exclusive with total probability 1, therefore, by the law of total probability, follows.
We next show that for any t = 1, 2, . . . the matrix Pt (π 0 ) falls between the lower and upper transition matrices, i.e., L ≤ Pt (π 0 ) ≤ M holds.Let us use short notations for certain events: for any i = 1, . . ., m and for a fixed t ≥ 1 set Applying the definition of conditional probability and the law of total probability, noting that P(H i ) > 0 is provided by the assumptions of the lemma, we get Therefore, it is enough to take the summation over A i , instead of the entire S. For x ∈ A i , however, , so they are non-negative and sum up to 1. Further, must hold, since l ij , u ij are defined as the minimum and maximum values, respectively, of p(x, A j ) = P(X t+1 ∈ A j | X t = x) over x ∈ A i .Since the weighted average must fall between the minimum and the maximum, therefore, we have that is, for any t ≥ 1 and for any initial distribution π 0 allowed by the conditions of the Theorem.
Let us now start the chain from an initial distribution π 0 that satisfies the conditions of the Theorem.We are going to compare the arising aggregated state distribution vectors (3) with the ones resulting from starting the chain from the stationary distribution π.Note that, due to the assumed irreducibility of the original chain, π(x) > 0 for all x ∈ S, so π is also a possible initial distribution that satisfies the conditions P(X t ∈ A i ) > 0.
If the chain happens to be exactly lumpable, then we get a "cleaner" result.Let πt be the state distribution of the lumped chain after t steps and let π be its stationary distribution.For concise description let us apply a frequently used distance concept among probability distributions.If p, q are two discrete probability distributions on the same set S, then their total variation distance D TV (p, q) is defined as It is well known that 0 ≤ D TV (p, q) ≤ 1 holds for any two probability distributions.It is also clear from the definition of the ergodic coefficient that it is the same as the maximum total variation distance occurring between any two row vectors of the transition matrix.
Note that exact lumpability is the special case of -lumpability with = 0. Therefore, we immediately obtain the following corollary.

Corollary 1.
If the Markov chain in Theorem 2 is exactly lumpable, then in the lumped chain for any t = 0, 1, 2, . . . the following holds: D TV ( πt , π) ≤ ρ t Furthermore, we can take ρ = 1 − p 0 − q 0 − as an upper bound on the ergodic coefficient of L. Then we obtain from Theorem 2, expressing the estimation in terms of the total variation distance: where the distributions πt , π are over the sets of the partition (A, Ā), not on the original state space.Note that in our case we actually have Therefore, we obtain the estimation directly for the set A: If p 0 + q 0 is not extremely small, then the term (1 − p 0 − q 0 ) t will quickly vanish, as it approaches 0 at an exponential rate.Therefore, after a reasonably small number of steps, we reach a distribution π t from any initial state, such that approximately the following bound is satisfied: It is quite interesting to note that neither the precise estimation (11), nor its approximate version (12) depend on the size of the state space.Now we demonstrate via numerical results that the obtained bounds indeed hold.Moreover, they are achievable after a small number of Markov chain steps, that is, with fast convergence.We simulated the example with the following parameters: n = 100 states, p 0 = q 0 = 0.25, and = 0.1.The set A was a randomly chosen subset of 50 states.The transition probabilities were also chosen randomly, with the restriction that together with the other parameters they had to satisfy conditions (A), (B), (C).
Figure 1 shows the relative frequency of visiting A, as function of the number of Markov chain steps.It is well detectable that the chain converges quite slowly.Even after many iterations the deviation from the stationary probability π(A) does not visibly tend to 0. On the other hand, it indeed stays within our error bound: |π t (A) − π(A)| ≤ p 0 + q 0 = 0.1 0.25 + 0.25 = 2×0.1 as promised.Having observed this, it is natural to ask, how soon can we reach this region, that is, how many steps are needed to satisfy the bound?This is shown in Figure 2. We can see that after only 10 iterations, the error bound is already satisfied.Note that this is very fast convergence, since the number of steps to get within the bound was as little as 10% of the number of states.

Conclusions
We have analyzed the convergence rate of quasi-lumpable Markov Chains.This represents the case when the large state space can be aggregated into a smaller number of clusters, in which the states behave approximately the same way.Our main contribution is a bound on the rate at which the aggregated state distribution approaches its limit in such chains.We have also demonstrated that in certain cases this can lead to a significantly accelerated convergence to an approximate estimation of the measure of subsets.The result can serve as a useful tool in the analysis of complex networks, when they have a clustering that approximately satisfies the conditions lumpability.

Definition 3 .
(Lower and upper transition matrices) Let M = (S, P) be a Markov chain which is -lumpable with respect to a partition Q = {A 1 , . . ., A m }.The lower and upper transition matrices L = [l ij ] and U = [u ij ] are defined as m × m matrices with entries

Figure 1 .
Figure 1.Deviation from the stationary measure for many iterations.

Figure 2 .
Figure 2. Very fast convergence to satisfy the error bound.