A Sequential Importance Sampling Algorithm for Estimating Linear Extensions

TL;DR

This paper introduces a practical sequential importance sampling algorithm to estimate the number of linear extensions of a partially ordered set, improving upon traditional methods by incorporating importance functions for better efficiency.

Contribution

The paper presents a novel SIS algorithm for counting linear extensions, enhancing the practicality of probabilistic sampling methods in this domain.

Findings

Implemented two importance functions and evaluated their performance.

Achieved more efficient estimation compared to traditional algorithms.

Demonstrated the method's applicability to complex counting problems.

Abstract

In recent decades, a number of profound theorems concerning approximation of hard counting problems have appeared. These include estimation of the permanent, estimating the volume of a convex polyhedron, and counting (approximately) the number of linear extensions of a partially ordered set. All of these results have been achieved using probabilistic sampling methods, specifically Monte Carlo Markov Chain (MCMC) techniques. In each case, a rapidly mixing Markov chain is defined that is guaranteed to produce, with high probability, an accurate result after only a polynomial number of operations. Although of polynomial complexity, none of these results lead to a practical computational technique, nor do they claim to. The polynomials are of high degree and a non-trivial amount of computing is required to get even a single sample. Our aim in this paper is to present practical Monte Carlo methods for one of these problems, counting linear extensions. Like related work on estimating the coefficients of the reliability polynomial, our technique is based on improving the so-called Knuth counting algorithm by incorporating an importance function into the node selection technique giving a sequential importance sampling (SIS) method. We define and report performance on two importance functions.

Equations

$$\left(\frac{r(a)+r(b)}{r(b)}\right)\left(\frac{r(a)+r(d)}{r(a)}\right)\left(\frac{r(c)+r(d)}{r(d)}\right)\left(\frac{r(c)}{r(c)}\right)$$

$$\mathrm{RV}(P)={⟨\frac{f_P(\lambda )}{L}⟩}_u-1$$

$$\mathbb{E}[f_P(\lambda {)}^2]=\sum _{\lambda \in \Lambda (P)}{f}_P(\lambda {)}^2\frac{1}{f_P(\lambda )}=\sum _{\lambda \in \Lambda (P)}{f}_P(\lambda )$$

$$\begin{array}{rl}\mathrm{RV}(P)& =\frac{\mathbb{E}[f_P(\lambda {)}^2]-\mathbb{E}[f_P(\lambda ){]}^2}{\mathbb{E}[f_P(\lambda ){]}^2}\\ & =\frac{\mathbb{E}[f_P(\lambda {)}^2]}{L^2}-1\\ & =\frac{1}{L^2}\left(\sum _{\lambda \in \Lambda (P)}{f}_P(\lambda )\right)-1\\ & =\frac{1}{L}\left(\sum _{\lambda \in \Lambda (P)}\frac{f_P(\lambda )}{L}\right)-1\\ & ={⟨\frac{f_P(\lambda )}{L}⟩}_u-1\end{array}$$

$$\mathrm{RV}(P)+1=\frac{r(M)}{L^{2}}\sum_{m\in M}{\frac{L_{m}^{2}}{r(m)}\big{(}\mathrm{RV}(P\setminus m)+1\big{)}}$$

$$f_P(\lambda )=\frac{r(M)}{r(m)}{f}_{P\setminus m}({\lambda }^{\prime })$$

$$\mathrm{RV}(P)+1=\frac{1}{L^2}\sum _{\lambda \in \Lambda (P)}{f}_P(\lambda )$$

$$\begin{array}{rl}\mathrm{RV}(P)+1& =\frac{1}{L^2}\sum _{m\in M}\sum _{\lambda \in \Lambda (P\setminus m)}\frac{r(M)}{r(m)}{f}_{P\setminus m}(\lambda )\\ & =\frac{r(M)}{L^2}\sum _{m\in M}\frac{1}{r(m)}\sum _{\lambda \in \Lambda (P\setminus m)}{f}_{P\setminus m}(\lambda )\\ & =\frac{r(M)}{L^2}\sum _{m\in M}\frac{L_m^2}{r(m)}\frac{1}{L_m^2}\sum _{\lambda \in \Lambda (P\setminus m)}{f}_{P\setminus m}(\lambda )\end{array}$$

$$\mathrm{RV}(P\setminus m)+1=\frac{1}{L_m^2}\sum _{\lambda \in \Lambda (P\setminus m)}{f}_{P\setminus m}(\lambda )$$

$$\mathrm{RV}(P)+1=\frac{r(M)}{L^{2}}\sum_{m\in M}{\frac{L_{m}^{2}}{r(m)}\big{(}\mathrm{RV}(P\setminus m)+1\big{)}}$$

$$A_i:=\underset{P\in {\mathcal{P}}_i}{\max }\left(\underset{m\in M}{\max }\frac{r(M)}{r(m)}\frac{L_m}{L}\right)$$

$$\mathrm{RV}(P)\le A_1{A}_2\cdots A_n-1$$

$$\mathrm{RV}(P\setminus m)+1\le A_1{A}_2\cdots A_{n-1}$$

$$\begin{split}\mathrm{RV}(P)+1&=\frac{r(M)}{L^{2}}\sum_{m\in M}{\frac{L_{m}^{2}}{r(m)}\big{(}\mathrm{RV}(P\setminus m)+1\big{)}}\\ &\leq\frac{r(M)}{L^{2}}\sum_{m\in M}{\frac{L_{m}^{2}}{r(m)}A_{1}A_{2}\cdots A_{n-1}}\\ &=\sum_{m\in M}{\frac{L_{m}}{L}\frac{r(M)}{r(m)}\frac{L_{m}}{L}A_{1}A_{2}\cdots A_{n-1}}\\ &\leq\sum_{m\in M}{\frac{L_{m}}{L}\max_{P\in\mathcal{P}_{i}}\left(\max_{m\in M}\frac{r(M)}{r(m)}\frac{L_{m}}{L}\right)A_{1}A_{2}\cdots A_{n-1}}\\ &=\sum_{m\in M}{\frac{L_{m}}{L}A_{1}A_{2}\cdots A_{n}}\\ &=A_{1}A_{2}\cdots A_{n}\sum_{m\in M}{\frac{L_{m}}{L}}\\ \end{split}$$

$$\sum _{m\in M}\frac{L_m}{L}=1$$

$$\mathrm{RV}(P)+1\le A_1{A}_2\cdots A_n$$

$$|\Lambda (P)|\ge \frac{n!}{d(v_1)d(v_2)\cdots d(v_n)}$$

$$|\Lambda (F)|=\frac{n!}{d_F(v_1)d_F(v_2)\cdots d_F(v_n)}$$

$$r(v)=\frac{i+d(v)-1}{i-d(v)+1}$$

$$|\Lambda (D)|=\frac{|D|!}{|D_1|!|D_2|!\cdots |D_k|!}|\Lambda (D_1)|\cdot |\Lambda (D_2)|\cdots |\Lambda (D_k)|$$

Full text

Although of polynomial complexity, none of these results lead to a practical computational technique, nor do they claim to. The polynomials are of high degree and a non-trivial amount of computing is required to get even a single sample. Our aim in this paper is to present practical Monte Carlo methods for one of these problems, counting linear extensions. Like related work on estimating the coefficients of the reliability polynomial, our technique is based on improving the so-called Knuth counting algorithm by incorporating an importance function into the node selection technique giving a sequential importance sampling (SIS) method. We define and report performance on two importance functions.

1 Introduction

For a partially ordered set, any total order which respects the partial order is known as a linear extension. Determining the number of linear extensions of a given partially ordered set (poset) is a fundamental problem in the study of ordering with many applications.

The most common applications are in scheduling problems, where the number of linear extensions gives the size of the solution space. For example, if we are trying to find a schedule that is optimal with respect to some cost function, the size of the solution space can be valuable in deciding how to carry out the search. If the solution space is small, then we may be able to perform an exhaustive search for the best order, but if the solution space is large, then we may wish to settle for an approximate solution.

Computing the number of linear extensions exactly is $\#$P-complete [3], but there are fully polynomial-time randomized approximation schemes to estimate the number, such as the Markov Chain Monte Carlo (MCMC) method presented in [3]. These algorithms, while of great theoretical importance, are not practical because the polynomials are of high degree, even for the improved versions due to Bubbly-Dyer [4] and Wilson [8].

Our aim is a practical method based on sequential importance sampling (SIS) in the spirit of the applications discussed in [2] and [5]. Our method may have exponential rather than polynomial complexity in some cases, but, as is typical of SIS, a modest amount of computation is often sufficient to learn at least something plausible about the number of linear extensions for any poset of interest. This is not the case for MCMC methods where significant computation may be required to obtain even a single sample chosen from a probability distribution close to the limit distribution.

In Section 2, we give definitions and notation necessary to describe our algorithm. In Section 3, we give the algorithm and explain how it works. The algorithm takes an importance function as an input, and in Section 5 we discuss the relative merits of several importance functions. In Section 4, we give an expression for the variance of the estimates produced by the algorithm and discuss methods for bounding the variance. In Section 6, we discuss possible improvements to the algorithm, and in Section 7 we give results of our numerical experiments with randomly generated posets.

2 Preliminaries

For a partially ordered set (poset), $(P,\leq)$, we use the notation $x\leq y$ to mean that $x$ precedes $y$ in the partial order or that $x=y$. If $x$ and $y$ are not equal but $x\leq y$, we write $x<y$. If $x<y$ and there is no $z$ such that $x<z<y$, we say that $y$ covers $x$.

By the descendants of an element $v\in P$ we mean all $w\in P$ such that $w\leq v$. We will refer to the number of descendants of $v\in P$ as $d(v)$. By the ancestors of $v\in P$ we mean all $w\in P$ such that $v\leq w$.

A linear extension of $P$ is any permutation $\sigma:P\to\{1,2,\dots,|P|\}$ such that for all $x,y\in P$, $x\leq y$ implies that $\sigma(x)\leq\sigma(y)$. We denote the set of all linear extensions of a poset $P$ by $\Lambda(P)$, and the number of linear extensions by $|\Lambda(P)|$. If $P$ has $n$ elements, then $|\Lambda(P)|$ is some integer between $1$ and $n!$.

In practice, posets are often arise from directed acyclic graphs (DAGs), and in this context, a linear extension is usually called a topological ordering or topological sort. Each DAG has a unique poset that corresponds to its reachability relation; that is, $w\leq v$ in the poset if and only if there is a directed walk from $v$ to $w$ in the DAG.

Determining the reachability poset of a DAG amounts to taking the transitive reduction or the transitive closure of the DAG. This can be accomplished by boolean matrix multiplication and hence has the same complexity as that procedure [1].

In general, different DAGs can have the same reachability relation, but we can create a unique DAG from a poset by directing an edge from $v$ to $w$ if and only if $v$ covers $w$ in the poset.

3 The Algorithm

Generating a single linear extension of a poset is well-understood and can be done in time $O(n)$ for a poset with $n$ elements. For example, any depth-first or breadth-first search of a DAG will generate a linear extension as a side-effect. For a breadth-first search, we can order the elements by their first visit times in ascending order, and for a depth-first search, we can order the elements by their last visit times in descending order. A greedy search will also generate a linear extension.

The standard method of obtaining a linear extension, which was first described in a 1962 paper by Kahn [6], corresponds to a breadth-first search. The procedure is as follows. We choose some minimal element of the poset and delete it, possibly causing some new elements to become minimal. We choose another minimal element of the poset and delete it. We repeat this until there are no elements left in the poset. The order in which the elements were deleted is then a linear extension of the poset.

We can also perform the same procedure with maximal elements instead of minimal, and in this case the reverse order of the deletions gives the linear extension. In this paper, we will use maximal elements instead of minimal.

Since all linear extensions can be obtained by sequential deletion of maximal elements, $|\Lambda(P)|$ is equal to the number of ways to execute this procedure. We can think of the choices available to us at each step as forming a decision tree, where each branching corresponds to a choice of a maximal element, and each leaf corresponds to one complete linear extension. Hence the number of leaves is $|\Lambda(P)|$.

In the classic Knuth algorithm [7] for estimating the number of leaves of a tree, we walk a path from the root to a leaf of the decision tree by choosing which branch to take at each level uniformly at random from the branches available. We then take the product of the number of branches available at each step along the path to be an estimate of the number of leaves in the tree. Knuth showed that this estimate is unbiased, and so the mean of the estimates over many samples will approach the number of leaves.

In our case, this means that to produce an estimate of $|\Lambda(P)|$, we begin our estimate with the number of maximal elements in the poset, then delete one of the maximal elements uniformly at random. We then multiply our estimate by the new number of maximal elements, and again delete one of the maximal elements uniformly at random. We continue in this fashion until the poset is empty.

One of the chief difficulties with Knuth’s algorithm is that if the decision tree is far from uniform, the variance of the estimates may be large, and an exponential number of samples will be required for an accurate estimate. Unfortunately, the degree of non-uniformity is unknown in general. One way to mitigate the problem of non-uniform decision trees is to choose paths through the tree non-uniformly, according to an importance function. This is known as sequential importance sampling (SIS).

Consider a function $r$ from the elements of a poset $P$ to the positive real numbers. We will refer to $r$ as the importance function and to $r(v)$ as the importance of $v\in P$. By the importance of a set $S\subseteq P$, denoted $r(S)$, we mean the sum over $S$ of the importance of its elements.

To obtain an estimate of $|\Lambda(P)|$ using SIS, we again as before traverse a single path from the root to a leaf of the decision tree. However, at each decision point, instead of choosing from among the maximal elements uniformly at random, we instead choose a maximal element with probability proportional to its importance.

To be more precise, if $v$ is a maximal element and $S$ is the set of all currently maximal elements, we choose $v$ from $S$ with probability $p(v)=r(v)/r(S)$. The factor multiplied into our estimate at that branch point is then $1/p(v)$ rather than $|S|$, and our estimate of $|\Lambda(P)|$ is the product over each $v$ chosen of $1/p(v)=r(S)/r(v)$.

For example, in the decision tree shown in figure 1, the highlighted path gives an estimate of

[TABLE]

Because the probability of obtaining any given linear extension is the reciprocal of the estimate associated with that linear extension, when we take the expected value of the estimates, each linear extension contributes exactly 1 to the sum, and so the expected value of the estimates is $|\Lambda(P)|$. Hence SIS gives an unbiased estimate for $|\Lambda(P)|$.

Pseudocode for finding a single estimate using this method appears in Algorithm 1. For a poset with $n$ elements, the outer loop runs exactly $n-1$ times, while the inner loops run at most $n$ times. Finding maximal elements can also be accomplished in linear time. Hence the complexity of the algorithm is $O(n^{2})$.

4 Variance

Here we calculate the relative variance of our estimates. Note that by “relative variance” we mean the ratio of the variance to the square of the mean. This is also the square of the coefficient of variation. We denote this quantity $\mathrm{RV}(P)$ for a poset $P$.

The relative variance of our algorithm is given explicitly by the following theorem.

Theorem 4.1**.**

For a poset $P$ with $L$ linear extensions, the relative variance of the estimates given by Algorithm 1 is

[TABLE]

where $f_{P}(\lambda)$ is the estimate associated with linear extension $\lambda$ and $\langle\cdot\rangle_{u}$ denotes the mean of the uniform distribution over all linear extensions.

Proof.

Let $\Lambda(P)$ be the set of all linear extensions of poset $P$, let $f_{P}(\lambda)$ be the SIS estimate associated with extension $\lambda\in\Lambda(P)$, and let $L=|\Lambda(P)|$. Note that $L$ is also the expected value of $f_{P}(\lambda)$.

Recall that the probability of selecting extension $\lambda$ is precisely $1/f_{P}(\lambda)$. Hence

[TABLE]

Then, using our definition of relative variance, we have

[TABLE]

∎

We can also express the relative variance as a recursive function of the relative variance of certain sub-posets, as detailed in the following theorem.

Theorem 4.2**.**

For a poset $P$, let $M$ be the set of maximal elements of $P$, where for $m\in M$, $L_{m}$ is the number of linear extensions of $P$ that begin with $m$. Then the relative variance of the estimates produced by Algorithm 1 using importance function $r$ is given by

[TABLE]

Proof.

First, observe that the linear extensions of $P$ that begin with a particular $m\in M$ are equivalent to the extensions in $\Lambda(P\setminus m)$, but with $m$ appended to the beginning of each extension in $\Lambda(P\setminus m)$.

Furthermore, $\lambda$ is an extension in $\Lambda(P)$ that begins with $m$, then the estimate produced by Algorithm 1 for $\lambda$ in $P$ is related to the estimate for the equivalent extension, call it $\lambda^{\prime}$, in $\Lambda(P\setminus m)$, by the following formula.

[TABLE]

Note that every extension in $\Lambda(P)$ begins with some element $m\in M$, and so for each $\lambda\in\Lambda(P)$, there is $m\in M$ such that equation (1) holds.

In the proof of Theorem 4.1 we showed that

[TABLE]

Then, by substituting equation (1) into equation (2) and manipulating the summation expressions, we get

[TABLE]

Theorem 4.1 tells us that

[TABLE]

Hence, from the last line of equation (3), we have

[TABLE]

∎

This recursive form of the relative variance suggests a way to bound it, which is described in our next theorem.

Theorem 4.3**.**

Let $\mathcal{P}_{i}$ be the set of all posets of size $i$. Let $A_{i}$ be given by

[TABLE]

where $M$, $L$, and $L_{m}$ are defined as before with respect to $P$. Then for all posets $P$ of size $n$, the relative variance of the estimates produced by Algorithm 1 using importance function $r$ is bounded by

[TABLE]

Proof.

The proof is by induction on the size of the poset $P$.

For the basis step, consider a poset of size 1. The single element, call it $m$, is also the only maximum, and so $r(m)=r(M)$ and $L_{m}=L$. Thus $A_{1}=1$. Since this poset has only one possible linear extension, the algorithm always produces the same estimate, making the relative variance zero. Thus the theorem statement is true for $n=1$.

Assume the theorem statement is true for all posets of size $n-1$ and less. Then for a poset $P$ of size $n$, the theorem statement holds for the sub-posets $P\setminus m$, where $m\in M$. That means

[TABLE]

Then, beginning with Theorem 4.2, we have

[TABLE]

Since every extension of $P$ must start with some $m\in M$,

[TABLE]

Hence,

[TABLE]

∎

From Theorem 4.3, we can see that the ideal importance function is $r(m)=L_{m}$, that is, the importance function which maps each maximal element to the number of extensions that begin with that element. This importance function would give us $A_{i}=1$ for all $i\in\mathbb{N}$, making $\mathrm{RV}(P)=0$ for all posets $P$.

Of course, this ideal importance function is impossible to achieve, since if we knew its values, we would not need an estimation algorithm. In Section 5, we discuss more practical importance functions.

5 The Importance Function

If the importance function $r$ is uniform across all elements at each decision point, then our algorithm is simply Knuth’s algorithm again. On the other hand, if a good importance function is available, one that reflects the actual structure of the decision tree, excellent results can sometimes be obtained, as in [2], [5]. Here we present two importance functions, each with its own strengths and weaknesses, which we will discuss.

5.1 The Number of Descendants

For the first importance function, we define $r(v)$ to be $d(v)$, the number of descendants of $v$, including $v$ itself. An example is shown in Figure 2.

The idea for this importance function came from the observation that each vertex must appear in a linear extension before all of its descendants, which themselves have fewer descendants. In fact, if we order the vertices by decreasing number of descendants, we obtain a linear extension.

We also obtain a useful lower bound from this importance function. The sum over the maximal elements of the number of their descendants must always be at least the number of vertices which have not been chosen yet, hence a lower bound for our samples and so for $|\Lambda(P)|$ itself is

[TABLE]

Here the denominator is the product over all $v_{i}\in P$ and is simply the observation that every $v_{i}$ is processed exactly once when generating a sample.

If the poset can be represented by a DAG that is a forest, then no maximal elements share descendants, and so every sample is equal to this lower bound, which is therefore exact. In fact, these observations serve as an alternative proof of the formula for the number of linear extensions of a forest. We also note in passing that this gives a way to sample from linear extensions of a forest uniformly at random.

In addition to a lower bound, this importance function also gives us an upper bound for each sample. Consider that when constructing a linear extension, each element of the poset becomes maximal exactly once. An element $v$ which was not maximal in the original poset $P$ becomes maximal when its last remaining ancestor and the ancestor’s edge to $v$ are deleted. If we collect this edge for each element of the poset, we have a spanning forest, $F$, of the poset $P$.

The number of linear extensions of this forest is clearly an upper bound for $|\Lambda(P)|$, and can be calculated exactly as

[TABLE]

where $d_{F}(v)$ is the number of descendants of $v$ in $F$. Hence this upper bound depends on which forest we construct.

5.2 The Available Spaces Quotient

The second importance function we will discuss is

[TABLE]

where $i$ is the number of elements that have not yet been added to the linear extension.

Unlike the first importance function, this importance function changes depending on how far along in constructing the linear extension we are. This clearly adds to the running time of the algorithm, but it does not increase the order of its complexity.

The motivation for this function came from the intuition that in addition to favoring elements with more descendants, we should also favor elements with fewer available spaces left. At any point during the construction of a linear extension, if there are $i$ elements left to be added to the extension, then the number of available spaces for element $v$ is $i-d(v)+1$, since $v$ must be chosen before all of its proper descendants.

By placing the formula for the number of available spaces in the denominator, the importance of an element increases dramatically as the number of spaces available for that element decreases. The particular formula for the numerator, on the other hand, was chosen as a result of numerical experimentation.

This importance function, unlike the first, does not give exact estimates for any known set of posets, nor has it given us any bounds on the number of linear extensions. Its sole virtue is that it significantly reduces the observed variance of estimates in numerical testing. See Section 7 for numerical results.

6 Improvements on the Algorithm

Here we discuss a method of reducing the variance of our estimates by modifying the way in which we construct linear extensions.

Consider a DAG $D$ that, when undirected, consists of several connected components $D_{1},D_{2},...,D_{k}$. Each linear extension of $D$ can be obtained by first partitioning the positions $1,2,...,n$ into sets of size $|D_{1}|,|D_{2}|,...,|D_{n}|$, and then filling the positions selected for each component with a linear extension of that component. Hence

[TABLE]

This suggests a recursive algorithm whose pseudocode appears in Algorithm 2.

The order of the time complexity for Algorithm 2 is, at worst, $O(n^{3})$. This is because we have added the additional step of searching for connected components, which executes $n-1$ times and in the worst case takes quadratic time in the size of the poset being searched. See Section 7 for numerical results on the variance of Algorithm 2.

7 Numerical Results

Numerical tests of Algorithm 1 were implemented in C++ using a sparse representation of the posets. Given the poset elements $v_{1},v_{2},\dots,v_{n}$, for each pair of elements $v_{i}$ and $v_{j}$ with $i<j$, the relation $v_{i}>v_{j}$ was given a 20% probability to exist using a pseudorandom number generator. The posets were then transitively completed.

For each value of $n$ from $10,15,\dots,150$, we generated $n^{2}$ posets in this manner, and $n^{2}$ SIS estimates were performed on each poset. The relative variance of the estimates was calculated for each poset, and the results averaged for each value of $n$.

The results for the different importance functions were of differing orders of magnitude; therefore, they are compared in Figure 3 on a log-log scale. Linearity on a log-log scale suggests a power function relationship where the exponent is the slope of the line.

Numerical tests of Algorithm 2, the Recursive Connected Components algorithm, were performed in the same manner as for Algorithm 1. The results in Figure 4 show a modest reduction in the average relative variance for each of the three importance functions when recursion is used. The relative size of this reduction appears to diminish as the poset size increases.

These numerical results show that the variance of our method does not grow too quickly, and that our importance functions significantly reduce the variance as compared to a uniform importance function.

Acknowledgements

We thank Myra Deng and Amanda Crawford for work done on early versions of the algorithm.