The Mixing of Markov Chains on Linear Extensions in Practice

1 Introduction

Computing the number of implicitly given objects-or more generally a weighted sum-is a central task in both the application and theory of artificial intelligence. Not only various counting problems, like probabilistic inference in graphical models, need be solved within intelligent systems, but hard counting problems also call for algorithm design paradigms that go beyond the usual worst-case view of theoretical computer science. Leveraging modern satisfiability solvers in approximate counting is a recent example of research in this direction [Gomes et al., 2006; Chakraborty et al., 2015].

When the interest is in good probabilistic guarantees for an approximate count, the theoretically most efficient algorithms known to date often rely on sampling objects along a rapidly mixing Markov chain [Sinclair and Jerrum, 1989]. The practical value of these algorithms is, however, unclear: typical algorithms always simulate a Markov chain the same number of steps given by a worst-case upper bound of the mixing time; the bound may only depend on basic parameters, like instance size, and be impractically large, even if polynomial. It is largely an open question, whether such algorithms can be

Mikko Koivisto

Dept. of Computer Science University of Helsinki mikko.koivisto@helsinki.fi

made to run fast in practice, that is, on moderate-size "typical" instances, for which analytic bounds may be pessimistic.

The present work is motivated by the following basic observation. Suppose a Markov chain mixes fast on typical instances. Then there is hope of discovering better, instancesensitive upper bounds or termination rules that are efficient even if good analytic upper bounds were unknown. In making such analytic or algorithmic discoveries, which can be technically very difficult, it could greatly help if one knows in advance which Markov chains exhibit varying mixing times and how the variance depends on the instance structure. To this end, empirical results on mixing times should be valuable.

In this paper, we investigate the question of practical Markov chain algorithms, focusing on the classical problem of counting the linear extensions of a given partial order. The problem is equivalent to counting the topological sorts of a given digraph and has applications for instance in sequence analysis [Mannila and Meek, 2000], preference reasoning [Lukasiewicz et al., 2014], partial-order plans [Muise et al., 2016], and Bayesian network learning [Niinim¨aki et al., 2016]. The problem is #P-complete [Brightwell and Winkler, 1991] but admits a fully polynomial time randomized approximation scheme [Dyer et al., 1991]. The worst-case running times of the fastest known schemes, for n-element instances and relative error within ϵ > 0, scale as ϵ-2n5 , ignoring factors polylogarithmic in n and 1/ϵ [Bubley and Dyer, 1999; Banks et al., 2010]. The former scheme relies on the frequently studied Karzanov-Khachiyan chain [Karzanov and Khachiyan, 1991; Bubley and Dyer, 1999; Wilson, 2004], whose worst-case mixing time is Θ(n3 log n). For the scheme to be practical, we would need a chain that mixes much faster on real-world instances. Is this possible?

We make three main contributions. Our first challenge is to estimate the mixing time of a fixed Markov chain on a per-instance basis. This is generally very hard both computationally [Bhatnagar et al., 2011] and information theoretically [Hsu et al., 2015], requiring work at least linear in the size of the state space. As a remedy, we measure the mixing time not in relation to the uniform distribution on linear extensions, but to the induced distribution on ordered pairs of elements only; this is feasible (on moderate-size instances) using a recent exact algorithm for counting linear extensions [Kangas et al., 2016]. We show that a small relative mixing time suffices for fast approximate counting and give a sampling based estimator of relative mixing time. With modifications our method could also apply outside the context of linear extensions.

Second, we explore empirically, using various classes of instances, the gap between actual mixing times and the known upper bounds. Our study includes the Karzanov-Khachiyan chain and a less studied "insertion chain." In addition, we introduce a novel "shuffle chain," which trades the time requirement of a single transition step for better mixing time.

Finally, we study two algorithms for exact sampling of linear extensions [Huber, 2006; 2014]. The algorithms are based on so-called perfect simulation of variants of the KarzanovKhachiyan chain and the modified insertion chain. As perfect simulation results in exact sampling, we may expect the time requirements be significantly larger than for approximate sampling. On the other hand, the perfect simulation algorithms are, by nature, adaptive; they run faster on easier instances for which the chain mixes faster. As far as we know, the performance of these algorithms has not been studied empirically prior to this work.

2 Preliminaries

Consider a Markov chain (Xt)∞t=0 on a finite state space Ω with stationary distribution π. Denote by pts the distribution of Xt when the initial state is X0 = s. The mixing time of the chain is the smallest t such that the total variation distance between the distributions pts and π, given by maxA⊂Ω |pts (A) - π(A)|, is less than 1/4 for all s ∈ Ω.

Throughout this paper the state space Ω will be the set of linear extensions of some partial order P and the distribution π is the uniform distribution. For convenience, we assume the ground set of P is [n] := {1, 2, . . . , n}. Another partial order P ′ on [n] is an extension of P if it contains P and a linear extension if it additionally is a total order. We write Ω(P) for the set of linear extensions of P, or Ω for short. We denote order relations by capital letters in a somewhat nonstandard manner. For (a, b) ∈ P we write ab ∈ P for short (but not aPb). Often it will be handy to view a total order L on [n] as a sequence x = x1x2 . . . xn satisfying xixj ∈ L if i ≤ j.

To express how close a random variable Z is to a target value v > 0, we say Z is an (ϵ, δ)-approximation of v if (1 + ϵ)-1 v ≤ Z ≤ (1 + ϵ)v with probability at least 1 - δ. A draw X from a finite set S is a variable that is uniformly distributed on S, that is, Pr[X = s] = 1/|S| for all s ∈ S.

3 Relative Mixing Time

Brightwell and Winkler [1991] give the following reduction from approximate counting of linear extensions to sampling. Construct a sequence of partial orders P0, P1, . . . , Pk where P0 = P, Pk ∈ Ω(P), and Pi+1 is obtained from Pi by imposing an ordering on a pair of elements (ai, bi) that are incomparable in Pi. Let Li be a draw from Ωi := Ω(Pi). We have

To estimate |Ω(P)|, it thus suffices to draw almost uniformly from each Ωi, which can be accomplished by simulating a Markov chain for a number of steps given by the mixing time.

In fact, as we show next, it suffices to draw from a distribution that is close to the uniform distribution only in relation to the probabilities of the events aibi ∈ Li. To formalize this idea, we introduce the notion of relative mixing time:

Definition 1 (Relative Mixing Time). Let M be a Markov chain (Xt)∞t=0 on Ω with stationary distribution π. The mixing time of M in relation to a set of events A ⊂ P(Ω), initial state X0 = s, and distance ϵ > 0 is the smallest t′ such that Dts (A) := maxA∈A |pts (A) - π(A)| ≤ ϵ for all t ≥ t′ .

The relative mixing time, for ϵ = 1/4, is a lower bound for the standard mixing time, and the two coincide if we put A = P(Ω) and take the maximum over all initial states. An advantage of relative mixing time is that it enables more efficient estimation of convergence, since we only need to consider a limited set of events A instead of the full set of states required for detecting convergence in the total variation distance; see the next subsection. Rabinovich et al. [2016] introduce a similar notion, the function-specific mixing time.

For approximate counting of linear extensions it suffices to consider mixing times in relation to event sets of the form

we call the all-pairs event set of the partial order P.

Theorem 1. Let P be a partial order on [n] and 0 < ϵ ≤ 1/2. Suppose that for every extension P ′ of P there is a Markov chain on Ω(P ′ ) with uniform stationary distribution and fixed initial state such that simulating one transition takes at most time I and the mixing time in relation to A(P ′ ) and distance ϵ′ = ϵ/(12n log2 n) is at most a known bound T . Then there exists an algorithm that computes an (ϵ, 1/4)-approximation of |Ω(P)| in time O(ϵ-2 ITn2 log2 n).

The idea of the proof is similar to that of Brightwell and Winkler [1991, Theorem 5]. We improve the time requirement by adopting a technique of Jerrum et al. [2004].

Proof (sketch). We construct the algorithm in two phases. First, we form a sequence of partial orders P0, P1, . . . , Pk as in the decomposition (1) such that k ≤ n log2 n and ri := Pr[aibi ∈ XiT ] ≥ 1/4 for all i = 0, 1, . . . , k - 1, where (Xi t ) is the given Markov chain on Ω(Pi). We will show how to succeed in this with probability at least 7/8. Second, we compute an (ϵ/2, 1/8)-approximation ˆR of R := r0r1 · · · rk-1 and output 1/ ˆR. As both phases succeed with probability at least 7/8, the whole algorithm succeeds with probability at least 3/4. By the definition of relative mixing time, we have |ri - Pr[aibi ∈ Li]| ≤ ϵ′ for each i; hence 1/R is within a factor (1 + 4ϵ′ )k ≤ (1 + 2ϵ/5) of |Ω|.

The first phase starts from P0 = P and simulates the merge sort algorithm on the elements in [n]. For each comparison between two elements a, b ∈ [n] made by the algorithm we check whether a and b are comparable in the current partial order Pi. If not, we estimate Pr[ab ∈ XiT ] by running the chain s1 times. We obtain the next partial order Pi+1 from Pi by adding the restriction aibi (and what follows by transitivity), where aibi equals ab if the estimate is at least 1/2 and ba otherwise. Because merge sort makes at most n log2 n comparisons and each comparison adds at most one new partial order into the sequence, we have that k ≤ n log2 n. Using Hoeffding's inequality [1963] and a union bound we get that all the estimates are within 1/4 absolute error with probability at least 1 - 2k exp(-s1/8), so we can set s1 = O(log n) and get that ri ≥ 1/4 for all i with probability at least 7/8.

For the second phase, we estimate R by ˆR = ˆ r0 ˆ r1 · · · ˆ rk-1, where the estimate ˆ ri for each ri is obtained using s2 samples. As E[ ˆR] = R, Chebyshev's inequality gives us Pr[| ˆR-R| < ϵR/2] ≥ 1 - 4 Var[ ˆR]/(ϵR)2 . Adopting an idea of Jerrum et al. [2004, p. 692], we get that, if ri ≥ 1/4 for all i, then

Thus we can set s2 = O(ϵ-2n logn) and get that ˆR is an (ϵ/2, 1/8)-approximation of R.

Using the mixing time bound O(n3 log n log(1/ϵ′ )) of the Karzanov-Khachiyan chain as an upper bound for the relative mixing time, Theorem 1 yields a randomized approximation sheme running in time O(ϵ-2n5 log3 n log(n/ϵ)).

3.1 An Estimator of Relative Mixing Time

Consider the estimation of the relative mixing time of a given Markov chain for a fixed set A, initial state s, and distance ϵ. We break the estimation task into two subtasks: the estimation of the distance Dts (A) = maxA∈A |pts (A) - π(A)| for a given t; and the estimation of the smallest t after which the distance stays below ϵ.

We solve the first task by simulating the Markov chain a sufficiently large number K of times independently to obtain an accurate estimate of pts (A) for each A ∈ A. Assuming the exact values π(A) are available, we then obtain a straighforward estimate of Dts (A). In the application to linear extensions, each event A of interest is {L ∈ Ω(P) : ab ∈ L} for some a, b ∈ [n], and the value π(A) can be computed exactly using the algorithms of Kangas et al. [2016].

For the second task, a brute-force search would estimate the distance Dts (A) for each and every t = 1, 2, . . . , U, where U is some known upper bound for the relative mixing time (or standard mixing time). This approach can however be very inefficient, as good upper bound may not be available and the multitude of the test points t requires strong control of error in every point. To overcome these challenges, we use binary search: first double t until Dts (A) decreases below a threshold and then use binary search to refine the result. To mitigate the effect of sampling error, we do not output a point estimate but a confidence interval. The lower bound of the interval corresponds to distance ϵ+ϵ′ and the upper bound to distance ϵ - ϵ′ for some small ϵ′ > 0. For the binary search to work correctly, we have to assume that t ↦→ Dts (A) is an essentially non-increasing function when its value is close to the threshold ϵ. More precisely, we require that it is ϵbounded after value ϵ - ϵ′ /2 in the following sense.

Definition 2. Let α < β. A function f : N→R is β-bounded after value α if f(t) ≤ α implies that f(t′ ) ≤ β for all t′ ≥ t.

This property holds for the total variation distance t ↦→ Dts (P(Ω)) [Mitzenmacher and Upfal, 2005]. Our empirical results (not shown) suggest that it also holds for the particular function t↦→Dts (A) of our interest for small positive ϵ and ϵ′ .

If T is an upper bound of the found interval, we evaluate Dts (A) for at most 4⌈log2 T⌉ values of t. By using Hoeffding's inequality [1963] to bound the probability of sampling error being at most ϵ′ /2, we get the following result:

Theorem2. If t ↦→Dts (A) is ϵ-bounded after value ϵ-ϵ′ /2, then the method outputs an interval that contains the mixing time relative to the set of events A and distance ϵ with probability at least 1 - 8⌈log2 T⌉|A| exp(-Kϵ′2 /2), where T is the upper bound of the output interval.

The method requires time O(K(IT + E|A|) + D|A|), where I is the time requirement of simulating a single step of the chain, E is the time requirement of testing membership of a state in an event of A, and D is the time requirement of computing the exact probability of an event in A.

The method does not require knowledge of any upper bound for the relative mixing time, and the running time depends only on the measured relative mixing time. However, the method may need to be rerun with a larger parameter K if the relative mixing time is larger than anticipated and the probability of success becomes small. We do not have any guarantees on the size of the confidence interval, but in Section 5 we will see that they tend to be very narrow in practice at least in the case of sampling linear extensions.

4 Markov Chains

We proceed to describe concrete Markov chains, which enable sampling linear extensions either almost uniformly or exactly uniformly. Empirical results are reported in Section 5.

4.1 Approximate Sampling

We describe the transitions of the well known KarzanovKhachiyan chain, a less studied insertion chain, and a novel shuffle chain. In each of these chains, a transition "does nothing" with probability 1/2, e.g., to ensure aperiodicity; below we describe what the chain does with the remaining probability. Let x = x1x2 . . . xn be the current state of the chain.

The Karzanov-Khachiyan Chain

The chain attributed to Karzanov and Khachiyan [1991] swaps adjacent elements as follows:

Draw i from {1, 2, . . . , n - 1} and swap the positions of the elements xi and xi+1 if xixi+1 ̸∈ P.

A single transition of the chain can be implemented to run in constant time. Bubley and Dyer [1999] originally proved the bound O(n3 log n) for the mixing time of a weighted variant that draws i proportionally to i(n-i). Wilson [2004] showed that the same bound holds for the unmodified chain.

The Insertion Chain

Bubley and Dyer [1999] mention the insertion chain that removes a random element from the ordering and reinserts it to a random position as follows:

Draw positions i and j from {1, 2, . . . , n}. Let y1y2 . . . yn-1 = x1 . . . xi-1xi+1 . . . xn. Move to state y1 . . . yj-1xiyj . . . yn-1 if it is a linear extension of P.

A straightforward implementation of the transition takes time Θ(n) in the worst case, but for dense partial orders it is often fast in practice. Bubley and Dyer give a lower bound Ω(n2 ) and an upper bound O(n4 log n log |Ω|) for the worstcase mixing time of the insertion chain.

We also consider a version of the chain that draws j so that the resulting state is always a linear extension of P. This tends to expedite mixing on dense P .

The Shuffle Chain

Let G be the comparability graph of P , that is, the undirected graph with the vertices [n] and an edge between vertices that are comparable in P. If G is disconnected, we can draw from Ω by independently sampling a linear extension of P restricted to each component of G and interleaving these total orders randomly. The shuffle chain uses this idea:

Choose an interval {i, i + 1, . . . , i + l - 1} ⊂ [n] by first drawing the length l from a predefined distribution q and then the starting position 1 ≤ i ≤ n - l + 1. Find the components of the subgraph of G induced by the vertices xi, xi+1, . . . , xi+l-1. Choose an interleaving of the components uniformly at random such that the elements in the same component stay in the same order, and reorder the vertices accordingly.

Since the transition matrix of the chain is symmetric, the stationary distribution of the chain is the uniform distribution.

The most time consuming part of the transition is finding the components of the subgraph. For a subgraph of l vertices, a depth-first search can take time O(l2 ) since there can be O(l2 ) edges. However, by using bit vectors to store the adjacency lists in G and bitwise logic operators for set union and intersection, this can be optimized to O(l(⌈n/w⌉ + logn)) where w is the number of bits in a word.

The choice of the distribution q affects the expected complexity of the transition, because a transition is more demanding for longer intervals. We let q(l) = n/[l(l - 1)(n - 1)] for all l ∈ {2, 3, . . . , n}. With this choice, the expected interval length is O(log n) and the expected time requirement of a transition is O((⌈n/w⌉ + log n) log n).

The shuffle chain mixes in polynomial time:

Proposition 1. If the distribution q is as above, then the mixing time of the shuffle chain is O(n4 log2 n).

Proof (sketch). With probability q(2) > 1/2, the transition is the same as in the Karzanov-Khachiyan chain. Combining the O(n3 log n) mixing time bound of that chain with a result of Dyer and Greenhill [2000, Theorem 5.3] yields a bound O(n3 log n(log |Ω| + log ϵ-1 )).

4.2 Exact Sampling

Sometimes one can turn an approximate Markov chain sampler into an exact sampler by using a technique called coupling from the past (CFTP) [Propp and Wilson, 1996]. In general this requires that either the Markov chain is monotonic or we can construct a so called bounding chain for the Markov chain. A bounding chain is a Markov chain itself that "bounds" in which states the original Markov chain can be. In the beginning of the simulation the bounding chain allows all states. Once the bounding chain has converged to bound only one state, we know it must be the current state of the original Markov chain, regardless of the starting state. Discovering a bounding chain that converges quickly can be challenging.

Two bounding chains have been proposed for the exact sampling of linear extensions.

Exact Karzanov-Khachiyan Chain

Huber [2006] gives a bounding chain for the KarzanovKhachiyan chain. The idea is to keep track of a right bound for the position of each element in the order. Like in the Karzanov-Khachiyan chain, each step takes a constant time.

The expected running time of the resulting CFTP algorithm is O(n3 log n). If implemented as described in the original article, the algorithm always takes the same number of steps for a fixed number of elements n (and a fixed output sequence of the random number generator). A modification that initializes the bounding chain adaptively according to the input partial order can terminate after a smaller number of steps (Mark L. Huber, personal communication, August 12, 2016).

Exact Relocation Chain

Rcently Huber [2014] gave another type of bounding chain, which is related to the insertion chain but instead operates in a continuous state space, the unit n-cube [0, 1]n . The idea is that any point u = (u1, . . . un) (except a zero-measure subset) corresponds to a unique total order x satisfying ux1 < ux2 < · · · < uxn . Moreover, ua ≤ ub for all ab ∈ P if and only if the corresponding x is a linear extension of P. Drawing a point from {u ∈ [0, 1]n : ua ≤ ub for all ab ∈ P} then corresponds to drawing a linear extension of P .

A straightforward Gibbs sampler makes a transition by drawing an element a from [n] and then resampling the coordinate ua. This roughly corresponds to reinserting an element to a randomly selected location as in the insertion chain. It turns out that such a chain is monotonic and yields an efficient CFTP algorithm. As in the case of the insertion chain, the straightforward implementation of a single iteration runs in time O(n). Huber shows that the expected running time of the algorithm is O(∆2n logn) on height-2 partial orders where every element is comparable with at most ∆ other elements, but he does not give any bound for the general case.

5 Experimental Results

We have studied empirically the performance of the Markov chain algorithms of Section 4 using partial orders of different types and sizes. For the Karzanov-Khachiyan chain we here only show the results for the original version as there was not much difference to the weighted version. For the insertion chain we only show the results of the modified version as it seemed to mix consistently faster than the basic variant. Because all the algorithms are iterative, we compare both the average running time of a single iteration (i.e., transition) and the number of required iterations.

We consider the following classes of DAGs of sizes 10 ≤ n ≤ 50 generating the partial orders.

AVGDEG(k) is a random DAG with expected average degree k, generated by adding an arc (a, b) with probability k/(n - 1) for all 1 ≤ a < b ≤ n.

MAXINDEG(k) is a random DAG of maximum indegree at most k. The DAG is generated by first drawing a number nb from {0, 1, . . . , min{k, b-1}} for each vertex b, and then adding an arc (a, b) for nb vertices a drawn from {1, 2, . . . , b - 1} without replacement.

BIPARTITE(p) is a randomDAG generated by adding an arc (a, b) with probability p for all a ≤ n/2 < b.

GRIDTREE(k) is a randomDAG generated as follows: Start with a directed k by k grid. Until the DAG contains at least n vertices, repeatedly join a new directed k by k grid along a randomly selected side of a grid already in the tree, so that the directions of the shared arcs agree. Remove potential extra vertices from the last added grid.

TWOPATHS is the DAG consisting of the two length-n/2 paths (1, 2, . . . , n/2) and (n/2 + 1, n/2 + 2, . . . , n).

LADDER is obtained from TWOPATHS by adding the arcs (1, n/2 + 2), (2, n/2 + 3), . . . , (n/2 - 1, n).

The first four classes are adopted from Kangas et al. [2016]. The last two classes are derived from the hard instances we found for the insertion and shuffle chains using a local search with the objective of maximizing the (estimated) mixing time.

5.1 Performance of a Single Iteration

To measure how time consuming a single iteration in each of the chains is in practice, we implemented the simulation algorithms in C++, also paying attention to code optimization. We measured the performance over all the test instances. The results are summarized in Figure 1.

As expected, the constant time transitions of the KarzanovKhachiyan chain are the fastest. For the other chains the time requirements are substantially larger, however, growing only slowly in n. While the shuffle chain is expected to behave logarithmically, the large constant factors involved make it slower than the linearly behaving insertion and relocation chains in the range of n in question.

5.2 Number of Iterations

We estimate the mixing times of the Markov chains in relation to the all-pairs event set A(P) of the input partial order P, initial state s = {ab : 1 ≤ a ≤ b ≤ n}, and distance

ϵ = 1/4, using the method of Section 3.1. Putting ϵ′ = 0.01 and K = 5 · 105 , the confidence of each returned interval is at least 0.999 (cf. Theorem 2). For the exact samplers, we simply report the average number of iterations over 105 runs. The results are shown in Figure 2.

The Karzanov-Khachiyan chain mixes typically an order of magnitude faster than suggested by the known worst case bound. Yet, the chain appears to require a nearly cubic number of steps on all instance types. The insertion chain and the shuffle chain not only scale quadratically in most cases but also mix faster than the Karzanov-Khachiyan chain already on very small instances. On the TWOPATHS and LADDER instances, where the insertion chain seems cubic, the shuffle chain remains about quadratic. Despite the advantage of the Karzanov-Khachiyan chain in transition performance, the difference of relative mixing times becomes more significant in larger instances. For example, when n = 50, the median total performance (the product of the running time of a single iteration and the relative mixing time) is 92 µs for the Karzanov-Khachiyan chain and 57 µs for the shuffle chain.

For the exact Karzanov-Khachiyan sampler and the exact relocation sampler the numbers of iterations seem to follow the relative mixing times of the Karzanov-Khachiyan chain and the insertion chain, respectively, with a constant-factor overhead between 10 and 100. While the comparison of exact and approximate samplers is not completely fair, it gives some indication of the overhead of these exact samplers.

6 Concluding Remarks

We introduced a notion of relative mixing time. Like standard mixing time, it is a sufficient in applications to approximate counting; but unlike standard mixing time, it renders estimation computationally feasible, provided that a related counting problem is feasible. We demonstrated this methodology in the context of linear extensions of moderate-size partial orders. The approach could be applicable in other domains as well. While we utilized an exact algorithm for counting linear extensions [Kangas et al., 2016], extending the method to tolerate approximate counts is conceptually straightforward.

Our experimental results reveal that the known worst-case bounds poorly characterize the actual mixing times of several Markov chains on linear extensions. In particular, the insertion chain and the shuffle chain exhibit much smaller mixing times than suggested by the worst-case upper bounds. These two chains appear to be the most efficient chains in practice; while the transitions have nonconstant complexity, they are amenable to tools from algorithm engineering. It remains an open question, whether the promise of the shuffle chain can be realized via better analytic upper bounds or variants that enable CFTP (coupling from the past). Regarding the latter, our findings about exact relocation, a recent CFTP algorithm [Huber, 2014], are encouraging: it has just a constant-factor overhead to the related insertion chain, and it is often superior to the worst-case bound of the Karzanov-Khachiyan chain.

Acknowledgements

This work was supported in part by the Academy of Finland, under Grant 276864 and Grant 303816.