A Scalable Scheme for Counting Linear Extensions

1 Introduction

The generic problem of computing the number of objects in an implicitly given (multi)set plays a fundamental role in artificial intelligence research. On the one hand, various instantiations of the problem appear in applications, particularly in probabilistic inference. On the other hand, the computational hardness of such problems calls for heuristics that are efficient in practice and able to beat the known, perhaps pessimistic, complexity upper bounds. In particular, researchers have sought alternatives to the popular Markov chain Monte Carlo (MCMC) method, which often does not provide practical means for controlling the accuracy and time requirement. For example, there has been recent interest in approximate counting methods that leverage the power of modern SAT solvers [Gomes et al., 2006; Chakraborty et al., 2015]. Generally, it remains largely an open question, what sort of approximate counting schemes are the most efficient in practice.

This paper focuses on the classic problem of counting linear extensions: given a partial order ≺ (i.e., an irreflexive and transitive binary relation) on n elements, compute the number of those total orders on the elements that are supersets of ≺. The problem has applications, for example, in sequence analysis [Mannila and Meek, 2000], sorting [Peczarski, 2004], preference reasoning [Lukasiewicz et al., 2014], convex rank tests [Morton et al., 2009], partial order plans [Muise et al., 2016], and learning graphical models [Wallace et al., 1996;

Niinim¨aki et al., 2016]. The problem is #P-hard [Brightwell and Winkler, 1991] but admits a fully polynomial time approximation scheme based on MCMC [Dyer et al., 1991]. The best known asymptotic bounds for the expected running time are O(ϵ-2n3 log2 ℓ log n) [Banks et al., 2017] and O(ϵ-2n5 log2 n) [Talvitie et al., 2018], where ℓ is the number of linear extensions and ϵ the allowed relative error.

The large degree of the polynomial bounds raises the question of the practical value of the schemes. This question was recently addressed by Talvitie et al. [2018] who studied variants of the schemes empirically on benchmark instances. Although their polynomial-time scheme, enhanced by a number of heuristic tricks, was superior to prior schemes, it was inferior to an exponential-time scheme they coin adaptive relaxation Monte Carlo (ARMC), even when allowing hours of running time. Their finding suggests that the MCMC method, though asymptotically faster, may not lend itself to implementations that are competitive in practice.

In this paper, we refute that conclusion. We present a novel MCMC scheme that stems from modifying and combining ideas scattered in previous works. Like Banks et al. [2017] we apply the generic Tootsie Pop algorithm (TPA) of Huber and Schott [2010], however, with two major differences: First, we employ a different continuous embedding, which also requires us to design a novel sampler-our perfect sampler is the main technical contribution of this work. Second, unlike all previous efficient reductions from counting to sampling, which add constraints one after another to the input partial order until it becomes a linear extension, our reduction constructs intermediate problems between the input partial order and its relaxation we obtain by removing some constraints; this is inspired by the ARMC scheme. We demonstrate empirically using the collection of benchmark instances of Talvitie et al. [2018] that our scheme, which we call relaxation Tootsie Pop, is superior to all previous schemes. Moreover, its running time appears to grow polynomially, even if we currently have no proof for a polynomial bound.

In addition, we make an attempt to apply state-of-the-art SAT solver based schemes for counting linear extensions. Specifically, we instantiate the algorithms of Chakraborty et al. [2016], Thurley [2006], and Lagniez and Marquis [2017] using two alternative, natural problem encodings.

Section 2 begins by formulating the problem of counting linear extensions and reviewing existing approaches; we also give our instantiations of the SAT solver based method. Section 3 presents the relaxation Tootsie Pop scheme, and Section 4 the required perfect sampler. The experiments are summarized in Section 5. We conclude in Section 6 by discussing the significance of our findings for counting linear extensions in particular and for approximate counting in general.

2 The Problem and Existing Techniques

Let A be a finite set of n elements and let ≺ be an irreflexive and transitive binary relation on A. We say that ≺ is a partial order on A and call the pair P = (A, ≺) a partially ordered set (poset). The elements (x, y) ∈ ≺ are a called order constraints of P and denoted simply x ≺ y. We depict P with its cover graph, a directed acyclic graph where the vertex set is A and the arc set is the transitive reduction of ≺ (Fig. 1). If for all distinct x, y ∈ A we have either x ≺ y or y ≺ x, we say that P is a linear order.

Let ≺′ be another partial order on the set A. If ≺ ⊆ ≺′ , we say that P ′ = (A, ≺′ ) is an extension of P or, conversely, that P is a relaxation of P ′ . If P ′ is an extension of P and also a linear order, we say that P ′ is a linear extension of P (Fig. 1). We denote by L(P) the set of linear extension of P and by ℓ(P) := |L(P)| the number of its linear extensions.

For a quantity q, an estimate ¯ q, and parameters ϵ, δ > 0, we say that ¯ q is an (ϵ, δ)-approximation of q if

We consider the problem of giving an (ϵ, δ)-approximation of ℓ(P) for a given poset P and parameters ϵ, δ.

We begin with an overview of known schemes for approximate counting.

2.1 Telescopic Product

Brightwell and Winkler [1991] presented an approximate linear extension counting scheme based on uniform sampling of linear extensions. The method forms a sequence of posets (P0, P1, . . . , Pk) starting from the original poset P0 = P and iteratively adding more constraints until reaching a linear order Pk. Thus L(P0) ⊃ L(P1) ⊃ · · · ⊃ L(Pk) and we can write ℓ(P)-1 as a telescopic product ∏ki=1 ℓ(Pi)/ℓ(Pi-1). Each factor ℓ(Pi)/ℓ(Pi-1) is estimated by drawing samples uniformly from L(Pi-1) and taking the fraction of the samples that are also elements of L(Pi).

A comparison sorting based approach results in a sequence of posets of length k = O(n logn), and drawing in total O(n2 log2 nϵ-2 log δ-1 ) linear extension samples is sufficient for obtaining an (ϵ, δ)-approximation [Talvitie et al., 2017]. Using Huber's [2006] O(n3 log n)-expected-time algorithm-asymptotically the fastest know sampler-we obtain an O(n5 log3 nϵ-2 log δ-1 ) expected total running time. Talvitie et al. [2018] improved this bound by a factor of log n by exploiting the special structure of the posets in the sequence. They also achieved a significant performance improvements in practice by switching to a Gibbs sampler of Huber [2014].

2.2 Tootsie Pop Algorithm

The Tootsie Pop algorithm (TPA) [Huber and Schott, 2010] is a scheme for computing the ratio of the measures of two nested sets A′ ⊃ A′′ in a continuous space under certain conditions. Specifically, there must exist a family of intermediate sets {Aβ | β ∈ R} where A′ = Aβ′ and A′′ = Aβ′′ , the family must be nested in the sense that Aβ1 ⊂ Aβ2 for any β1 < β2, and the measure µ(Aβ) must be a continuous function of β. The scheme works by iteratively constructing a random sequence (βi)∞i=0 that starts from β0 = β′ , and each subsequent element βi+1 is computed by drawing Xi uniformly at random from Aβi and setting βi+1 = inf{β ∈ R|Xi ∈ Aβ} (whether it is easy to draw Xi and compute βi+1 depends on the set family). Now the number of iterations min{i ≥ 0 | βi+1 ≤ β′′ } required to reach A′′ is Poisson distributed with expected value r = ln(µ(A′ )/µ(A′′ )), and we can estimate the log-ratio by averaging over multiple independent runs. In total, O(r2 ϵ-2 log δ-1 ) calls to the sampling subroutine are sufficient for an (ϵ, δ)-approximation.

Banks et al. [2017] apply TPA to linear extension counting by transforming the discrete counting problem into a measure computation problem in a continuous space. This is achieved by affixing a continuous dimension to the set of linear extensions: Aβ = ⋃L∈L(P) {L} × [0,wL(β)) where β ↦→ wL(β) is a continuous function with values in range [0, 1]. The weights wL(β) are chosen carefully such that uniform sampling from Aβ is possible and A′ = ⋃L∈L(P ) {L}× [0, 1) and A′′ = {Lhome} × [0, 1), where Lhome is an arbitrary linear extension of P. Then ℓ(P) = µ(A′ )/µ(A′′ ) can be estimated by TPA. The sampling algorithm runs in O(n3 log n) time, which means that the (ϵ, δ)-approximation algorithm for counting the number of linear extensions runs in O(n3 log2 ℓ(P ) log nϵ-2 log δ-1 ) expected time.

2.3 Adaptive RelaxationMonte Carlo

The adaptive relaxation Monte Carlo scheme (ARMC) [Talvitie et al., 2018] utilizes an exact exponential-time algorithm for counting linear extensions [De Loof et al., 2006; Kangas et al., 2016], which in practice beats approximation schemes on small posets. The ARMC scheme starts by choosing a relaxation R for the given poset P, in such a way that each connected component in the cover graph of R has a size bounded by a parameter k. With this division into small components the exact algorithm may compute ℓ(R) in O(2k k2 ) time and also enables fast uniform sampling from L(R). The ratio µ = ℓ(P)/ℓ(R) (and by proxy, ℓ(P)) is then estimated by sampling, similarly to the telescopic product scheme.

The expected number of samples from L(R) required for an (ϵ, δ)-approximation is roughly O(µ-1 ϵ-2 log δ-1 ). If R is not too different from P, then µ is small and sampling is fast. On the other hand, letting R deviate from P allows a smaller k and thus makes the exact computation of ℓ(R) faster. An adaptive procedure is used to choose R so as to obtain a good tradeoff between these two phases.

2.4 BooleanModel Counting

The Boolean model counting problem, #SAT, asks for the number of satisfying truth assignments to a given Boolean formula. Given a poset P = (A,≺), we encode the problem of computing ℓ(P) in #SAT by producing a Boolean formula ϕ (in conjunctive normal form) whose satisfying assignments are bijectively mapped to L(P). We investigate two natural encodings and evaluate them by employing existing model counting software, specifically, implementations of two exact model counters [Thurley, 2006; Lagniez and Marquis, 2017] and an algorithm by Chakraborty et al. [2016], which outputs an (ϵ, δ)-approximation of the number of satisfying assignments.

In the first encoding, SAT #1, we let ϕ contain a variable vxy = ¬vyx for each pair of distinct elements x, y ∈ A to represent the order between x and y: if vxy is true, then x precedes y in the linear extension, otherwise y precedes x (n(n - 1)/2 variables). Let ϕ further contain the clauses

where clauses (i) enforce transitivity and clauses (ii) guarantee that the corresponding linear order is an extension of P. The total number of clauses is n(n - 1)(n - 2) + c, where c is the number of edges in the cover graph of P .

In the second encoding, SAT #2, we let ϕ contain a variable vxp for each x ∈ A and p ∈ [n] with the interpretation that vxp is true if and only if x is in position p in the linear extension (n2 variables). Let ϕ contain the clauses

1. ∨np=1 vxp for all x ∈ A (each element has a position),
2. ∨x∈A vxp for all p ∈ [n] (each position has an element),
3. ¬vxp ∨¬vxq for all x ∈ A and p, q ∈ [n] such that p ̸= q (no element is in two positions),
4. ¬vxp∨¬vyp for all p ∈ [n] and x, y ∈ A such that x ̸= y (no position has two elements).

where clauses (i-iv) guarantee that a satisfying assignment maps each element to exactly one position and vice versa, and clauses (v) enforce compliance with P . The total number of clauses is 2n + 2n2 (n - 1) + cn(n - 1)/2.

3 Relaxation Tootsie Pop

In this section, we present a new linear extension counting method based on TPA [Huber and Schott, 2010]. In order to translate the counting problem into the continuous space, we use the same embedding as the linear extension Gibbs sampler due to Huber [2014]: For a poset P = (A,≺), each point p ∈ [0, 1]A specifies a position px ∈ [0, 1] for every element x ∈ A, and every order constraint x ≺ y in P maps to a numerical constraint px ≤ py. The domain A ⊆ [0, 1]A

is defined as the set of points satisfying the numerical constraints. The linear order f(p) = (A,≺′p ) corresponding to a point p ∈ Ω is obtained from the ordering of the positions: x ≺′p y if and only if px < py (we can choose arbitrarily the cases where px = py for x ̸= y since their measure in A is zero). Consequently, if we sample p ∈ A uniformly at random, the distribution of f(p) in L(P) is also uniform.

Similarly to the ARMC method [Talvitie et al., 2018], we count the linear extensions of P = (A,≺) by estimating ℓ(R)/ℓ(P) where R = (A,≺′ ) is a relaxation of P for which the number of linear extensions ℓ(R) can be computed quickly. However, instead of direct sampling, we use TPA for estimating the ratio, as illustrated in Figure 2. This has the benefit that the required number of samples for TPA is only polylogarithmic in the ratio ℓ(R)/ℓ(P), compared to the linear growth in ARMC.

We define the family of intermediate sets required by TPA as follows:

As a result, the ratio ℓ(R)/ℓ(P) is equal to the corresponding ratio of measures µ(A1)/µ(A0). In Section 4 we present an algorithm for uniform sampling from [0, 1]A in the presence of generalized order constraints, of which the sets Aβ are a special case. For now, we assume that we can sample from Aβ and thus the requirements of TPA introduced in Section 2.2 are satisfied, yielding the following result:

Theorem 1. Given a relaxation R for P with ℓ(R) known, the Relaxation Tootsie Pop algorithm computes an (ϵ, δ)- approximation for ℓ(P) in O(log2 (ℓ(R)/ℓ(P))ϵ-2 log δ-1 ) calls to a generalized order constraint sampling subroutine.

3.1 Relaxation Finding

Choosing the empty relaxation R = (A, ∅) with ℓ(R) = |A|! is always valid, but finding a relaxation where ℓ(R) is closer to ℓ(P) will make the algorithm run faster. Even though Talvitie et al. [2018] investigated relaxation finding heuristics for the ARMC method, their heuristics rely heavily on exponential exact counting algorithms. Investing that much time to relaxation finding does not pay off in our case: the gains are much smaller due to the polylogarithmic dependence of the running time on the ratio ℓ(R)/ℓ(P). For that reason, we will use faster heuristics to find reasonable relaxations.

Our heuristics stem from the observation that for many special classes of posets P = (A, ≺), the number of linear extensions is easy to compute exactly or the problem can at least be reduced to counting the linear extensions of some induced subposets, i.e., posets P ′ = (A′ ,≺) where A′ ⊂ A. First we check if the poset is already a member of such a special class:

1. If A can be partitioned into sets A1 and A2 that are unrelated in the sense that for all a1 ∈ A1 and a2 ∈ A2, a1 ̸≺ a2 and a2 ̸≺ a1, then we recursively find relaxations R1 = (A1,≺1) and R2 = (A2,≺2) for the respective subposets, and combine them into relaxation R = (A, ≺1 ∪ ≺2) where ℓ(R) = ( |A| |A1|)ℓ(R1)ℓ(R2).
2. If A can be partitioned into sets A1 and A2 that are ordered in the sense that for all a1 ∈ A1 and a2 ∈ A2, a1 ≺ a2, then similarly to the previous part we combine the relaxations for the subposets into relaxation R = (A, ≺1 ∪ ≺2) where ℓ(R) = ℓ(R1)ℓ(R2).
3. If the cover graph of the poset P is a tree, then we use P as the relaxation and count its linear extensions using the O(n2 ) algorithm due to Atkinson [1990].
4. If any exact linear extension counter finishes computing ℓ(P) within a short time limit, we can use P as the relaxation. We run the counter due to Kangas et al. [2016] until its dynamic programming table has 10,000 items.

If none of these attempts succeeds, we try the next two methods and select the relaxation with fewer linear extensions.

1. If A can be partitioned into three sets A1, A2, and A3 such that A1 and A2 are ordered (a1 ≺ a2 for all a1 ∈ A1 and a2 ∈ A2), then we recursively find relaxations Rk = (Ak, ≺k) for k ∈ {1, 2, 3} and combine them into relaxation R = (A,≺′ ) where we keep the constraints between A1 and A2 but remove all constraints between A3 and A1 ∪ A2, i.e., ≺′ = (⋃3k=1 ≺k ) ∪ (A1 × A2). Now ℓ(R) = ( |A| |A3|)ℓ(R1)ℓ(R2)ℓ(R3). Our partition finding attempts to maximize the heuristic (|A1|+|A2| |A1| ), which is inversely correlated with the number of linear extensions if we assume that the relaxation for all the parts Ak is the empty partial order. For every x ∈ A, we use greedy hill climbing, starting from A1 = {x} and adding more vertices to A1 while maintaining A2 as the maximal set such that A1 and A2 are ordered.
2. We relax P by taking a spanning tree of its cover graph and then apply the O(n2 ) time algorithm. We sample

6 minimum-weight spanning trees, with the arcs of the cover graph sorted randomly by weight, and choose the one with the smallest number of linear extensions.

4 The Sampler

As the last remaining piece of relaxation Tootsie Pop algorithm, we consider the problem of sampling from sets of type Aβ defined in equation (1). More exactly, we consider the more general problem of sampling uniformly from the set

4.1 Continuous Relocation Chain

To sample from Ω, we define the continuous relocation Markov chain (Xk)∞k=0 in Ω in which each state transition is given by the rule Xk+1 = φ(Xk, ik, tk), where dimension ik ∈ [n] and position tk ∈ [0, 1] are sampled uniformly at random, and the transition function φ : Ω × [n] × [0, 1] → Ω maps (x, i, t) to y as follows:

• Construct the proposition state y′ by setting y′i = t and y′j = xj for all j ∈ [n] \ {i}.
• Move to state y = y′ if y′ ∈ Ω; otherwise stay in y = x.
• See Figure 3 for an example. From the symmetry of the tran-

sition distributions we get the following result:

Theorem 2. The uniform distribution on Ω is a stationary distribution for the continuous relocation chain.

4.2 Perfect Sampling

Although by Theorem 2 the stationary distribution of the chain is uniform, we do not have a bound for its mixing time, and so we do not know how many iterations are required for the desired precision. Nevertheless, we next see that one can sample exactly from the uniform distribution using monotonic coupling from the past [Propp and Wilson, 1996]. For this we need the chain to be monotonic in relation to the update function φ and some partial order ⊏ on Ω:

Lemma 1. Let x ⊏ y if xi ≤ yi for all i ∈ [n] and x ̸= y.

• (a) The relation ⊏ is a partial order on Ω;
• (b) there exist 0, 1 ∈ Ω such that 0 ⊑ x ⊑ 1 for all x ∈ Ω;
• (c) if x, y ∈ Ω such that x ⊑ y, then φ(x, i, t) ⊑ φ(y, i, t) for all i ∈ [n] and t ∈ [0, 1].

Proof. It is easy to verify that ⊏ is a partial order and that 0 = (0, 0, . . . , 0) and 1 = (1, 1, . . . , 1) satisfy the conditions in (b). Assume that x, y ∈ Ω and x ⊑ y. To prove that φ(x, i, t) ⊑ φ(y, i, t), it suffices to prove that φ(x, i, t)i ≤ φ(y, i, t)i because for all j ∈ [n] \ {i} we have that φ(x, i, t)j ≤ φ(y, i, t)j . For each z ∈ {x, y} let

By the definition of Ω we have that zi ∈ [Lz,Rz] and if t ∈ [Lz,Rz], then φ(z, i, t)i = t and else φ(z, i, t)i = zi. Because x ⊑ y we also have that Lx ≤ Ly and Rx ≤ Ry.

We prove that φ(x, i, t)i ≤ φ(y, i, t)i in cases: First, if t ∈ [Ly,Rx], then t ∈ [Lx,Rx] ∩ [Ly,Ry] and thus φ(x, i, t)i = t = φ(y, i, t)i. Second, if t < Ly (case t > Rx is symmetric), then φ(x, i, t)i is either xi or t. It holds that xi ≤ yi because x ⊑ y and t < Ly ≤ yi because yi ∈ [Ly,Ry]. Combined, these yield that φ(x, i, t)i ≤ yi = φ(y, i, t)i.

The method of Propp and Wilson [1996] now yields the following algorithm: We run two instances (Xk)0k=-T and (Zk)0k=-T of the continuous relocation chain for T > 0 iterations starting from time -T with initial states X-T = 0 and Z-T = 1. We use the same sequences of parameters ik and tk for both chains. If X0 = Z0, then X0 is the output sample, and otherwise we rerun the whole algorithm with the value of T doubled (i.e. starting from further in the past) still using the same parameters ik and tk. After the algorithm has terminated, we know that if we started the chain (Yk)0k=-T from any initial state and used the same parameters ik and tk, the output Y0 would match the output X0 of the algorithm because by Lemma 1 we have that X-T ⊑ Y-T ⊑ Z-T and inductively Xk ⊑ Yk ⊑ Zk for all -T + 1 ≤ k ≤ 0, yielding X0 ⊑ Y0 ⊑ X0. By considering a chain (Yk)∞k=-T where the initial state Y-T is sampled uniformly at random from Ω, we see that the distribution of the output X0 = Y0 is also uniform on Ω because it is a stationary distribution of the chain.

While we do not know a good upper bound for the expected running time of the algorithm, we can show that it is finite.

Theorem 3. The monotonic coupling-from-the-past algorithm for the continuous relocation chain outputs a sample from the uniform distribution on Ω in finite expected time.

Proof. Let c = 2n - 1 and define the event Al by (i-lc, i-lc+1, . . . , i-lc+c-1) = (1, 2, . . . , n, n - 1, . . . , 1) and tk ∈ [(ik - 1)/n, ik/n] for all -lc ≤ k < -(l - 1)c. We see that if on some iteration T ≥ lc an event Al occurs, then X-(l-1)c = Z-(l-1)c, which means that the algorithm terminates after that iteration. All the events A1,A2, . . . are mutually exclusive and have the same probability p > 0.

Let L be the smallest l such that the event Al occurs. Because T is doubled on every iteration, at least half of the time is spent in the last iteration, and thus the running time is O(L) for constant n. Because L is geometrically distributed with parameter p, its expected value is 1/p <∞.

Our algorithm is similar to Huber's [2014] sampler. The main difference is that whereas our method resamples xi along a randomly chosen dimension from the whole set [0, 1], rejecting the moves that fall outside, that method uses Gibbs sampling, meaning that it samples uniformly at random from only the set of allowed values. Even though this is beneficial in the sense that the chain stalls less in the same state, the downside is that the sampling never terminates (i.e., reaches X0 = Z0), and therefore additional tricks, which do not work with generalized order constraints, are required.

5 Experiments

We empirically evaluated the following schemes:

Telescopic Product: The scheme of Brightwell and Winkler [1991], with optimizations [Talvitie et al., 2018].

Extension Tootsie Pop (ETP): The scheme based on TPA due to Banks et al. [2017].

Adaptive Relaxation Monte Carlo (ARMC): The Monte Carlo scheme based on exponential-time counting by Talvitie et al. [2018].

SAT #1 and #2: Reducing the problem into SAT model counting using encodings given in Section 2.4 and solving it using sharpSAT [Thurley, 2006] or D4 [Lagniez and Marquis, 2017] exact model counter or ApproxMC2 approximate model counter [Chakraborty et al., 2016].

(Trivial) Relaxation Tootsie Pop (RTP): The TPA based scheme given in Section 3. For comparison, we also consider a version with the trivial relaxation (empty poset).

We adopted the experimental setup and benchmark instances of Talvitie et al. [2018]: The instances are randomly generated and include posets of average indegree k ∈ {3, 5}, bipartite posets of density p ∈ {0.2, 0.5}, and posets extracted from subgraphs of the networks Andes, Diabetes, Link, Munin, and Pigs from the Bayesian Network Repository (www.cs.huji.ac.il/˜galel/Repository). For every instance class and size between 8 and 512 elements, we have five different posets. We ran each algorithm for every instance for at most 24 hours of CPU time and memory limited to 8 gigabytes. All the algorithms were instantiated to produce a (1, 1/4)-approximation. We use the implementations due to Talvitie et al. [2018] for the telescopic product, extension Tootsie Pop and adaptive relaxation Monte Carlo methods. We will make our implementation for relaxation Tootsie Pop and the SAT encoding generators publicly available. 1

Figure 4 shows the results of the experiments. The growth of the running time of RTP clearly resembles polynomial growth, as opposed to the exponential growth of the running time of ARMC. While ARMC is faster in small cases, when the number of elements exceeds roughly 200 elements, RTP becomes the fastest scheme in most instance classes. The version of RTP that uses the trivial relaxation is already faster than the telescopic product scheme and ETP, and using proper relaxation instead of the trivial relaxation improves the running time roughly by an order of magnitude. The performance of the SAT based schemes is poor: they take hours to run already in the smallest cases.

1 github.com/ttalvitie/le-counting-practice

6 Discussion

The relaxation Tootsie Pop scheme is a novel application of TPA [Huber and Schott, 2010] for approximate counting of linear extensions. Our empirical results suggest that the scheme is the first that scales to large instances in practice. In contrast to the exponential-time ARMC scheme, relaxation Tootsie Pop requires only little memory and enables efficient large-scale parallelization (not demonstrated here). In applications like Bayesian network structure learning [Niinim¨aki et al., 2016], the new scheme should enable the handling of partial orders (i.e., networks) with several hundreds of elements, which are beyond the reach of ARMC; for smaller partial orders ARMC is superior, though.

Our key technical contribution was the perfect sampler in the constrained continuous space. We observed that, even if the sampler does not have-or we do not know how to prove-polynomial worst-case time complexity bounds, the sampler runs fast on all tested benchmark instances. This success should encourage working on heuristic designs of perfect samplers also for other counting and sampling problems. Compared to theoretical works that are concerned with provable uppers bounds, the heuristic approach gives access also to samplers whose complexity is difficult to analyze.

For the poor performance of the SAT solver based schemes we offer three possible explanations. First, our encodings of the problem were straightforward and could potentially be beaten by more sophisticated encodings. Second, while the tested algorithms represent the state of the art, other related schemes [Gomes et al., 2006; Zhao et al., 2016] could be more efficient in the present application. Third, the idea of reducing a counting problem to the NP-hard CNF-SAT problem may be ill-suited in general-even if powerful solvers exist-when the counting problem is "easy" in the sense that it admits a fully polynomial-time approximation scheme.

Acknowledgments

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