Approximating Integer Solution Counting via Space Quantification for Linear Constraints

1 Introduction

Solution Counting, the problem of computing the number of solutions for given constraints, is an important problem with various applications in artificial intelligence including probabilistic inference [Roth, 1996; Chavira and Darwiche, 2008], probabilistic logic programming [Aziz et al., 2015], planning [Domshlak and Hoffmann, 2007], knowledge representation [Lee and Wang, 2018] and privacy/confidentiality verification [Fredrikson and Jha, 2014].

For propositional formulas, the counting problem is known as model counting or #SAT, to which probabilistic inference is easily reducible. #SAT is the canonical #P-complete problem, thus exact counting is intractable for large instances. Recently, the hash-based approximate counters [Chakraborty et al., 2016b; Ge et al., 2018a] attract much attention for improving scalability and providing rigorous guarantees. For Satisfiability Modulo Theories (SMT) formulas, the counting problem is known as #SMT [Chistikov et al., 2015; Chakraborty et al., 2016a; Ge et al., 2018b], which is also an emerging topic in this realm.

As the fundamental constraints that are used in many areas, linear constraints (LCs) have been studied thoroughly.

∗Corresponding Authors

The integer solution counting over LCs is also an important problem and has many applications, such as counting-based search [Zanarini and Pesant, 2007; Pesant, 2016], simple temporal planning [Huang et al., 2018] and probabilistic program analysis [Geldenhuys et al., 2012; Luckow et al., 2014]. This problem is proved to be #P-hard [Valiant, 1979]. In practice, state-of-the-art tools for counting integer solutions on LCs, like LattE [Loera et al., 2004], often have difficulty when the number of dimensions is greater than 10 (i.e., when there are more than 10 numeric variables). As a result, they are usually not sufficient for practical applications.

To overcome such problem, we would like to investigate the relation between the integer solution count and the volume of the solution space. Consider the following set of LCs:

The volume of its solution space is 16291 which is close to the number of integer solutions, 16400. From more experiments, we observe that the execution time of space quantification (means volume computation and volume estimation) tools is usually much smaller than the counters, while their output values are close (for details, see Section 5). Naturally we wonder if there exists a relation between the space volume and the solution count, so that space quantification tools can be used to approximate integer solution counts.

In this paper, the first theoretical bound of approximating integer counts via space quantification for LCs is proposed and also proved. Based on these theoretical analysis, we present an approximate counting algorithm for LCs, which can give results with rigorous guarantees. This technique is orthogonal with respect to the current design of model counters and also space quantification tools. As a result, it is easy to apply. We evaluated our approach over random and structural cases generated from various applications. Experiments show that our approach achieves more than 20x speedup compared with state-of-the-art integer solution counters. Moreover, the advantage increases with the number of dimensions. The results also show that our approach could serve as a basis for many applications.

2 Background

The state-of-the-art tools for integer solution counting and space quantification are listed in Table 1. LattE [Loera et al., 2004] is an implementation of Barvinok's algorithm [Barvinok, 1993] dedicated to the counting of integer solutions for LCs. For solving space quantification problems, there are more candidates: qCORAL [Borges et al., 2014; 2015], PolyVest [Ge and Ma, 2015] and Vinci [B¨ueler et al., 2000]. qCORAL is a volume estimation tool, based on direct Monte-Carlo sampling. It can handle arbitrary constraints with real domain, such as non-linear polynomial constraints, trigonometric constraints, etc. PolyVest is also a volume estimation tool. It only handles LCs, but provides rigorous guarantees. It is an implementation of Multiphase Monte-Carlo algorithm [Dyer et al., 1991; Lov´asz and De´ak, 2012] with various improvements, whose time complexity is O∗ (n4 ) (O∗ indicates that we suppress factors of logn). Different from previous two tools, Vinci is a tool to compute the exact volume of solution space. It consists of several algorithms, but only Lasserre's method is called in this paper as others require inputs of convex hull.

Our experiments on LattE, Vinci and PolyVest show that the integer solution counts are usually very close to the volume for LCs (for details, see Section 5). Similar to a numeric integration procedure, it is trivial to estimate the volume of space with an algorithm based on integer solution counting. The estimation can be arbitrarily precise as we keep differentiating the space. On the contrary, there is no guarantee to approximate the number of integer solutions with space quantification. For example, a very 'thin' rectangle whose both sides are parallel to the coordinates and the short side lies in interval (0, 1). There is no integer point in it, but its volume can be arbitrarily large as the long side stretches. The theoretical result in [Kannan and Vempala, 1997] also shows such phenomenon. So it is necessary to find a bound to quantify the quality of an approximation.

3 Analysis

In this section, we propose and prove our main theoretical results, Theorem 1, 2 and Corollary 1, which provide bounds for approximating counts via volume, on LCs.

3.1 Preliminaries

We first provide the formal definitions of the approximation problem considered in this paper.

Definition 1. A LC can be written in the form a1x1 +a2x2 + · · · + anxn op a0, where xi are numeric variables, ai are constant coefficients, and op ∈ {<, ≤, >, ≥, =}.

Given a set of LCs F with real variables.

• Let #F denote the number of integer solutions of F.
• Let vol(F) denote the volume (Lebesgue measure) of solution space of F .

Problem 1. Given a set of LCs with real variables, we would like to obtain a bound c such that |vol(F) -#F | ≤ c.

When variables take values from integer domains, the definition of our approximation problem is slightly different. We consider the relaxed formula relaxed(F), where relaxed(F) is F whose variables take real values.

Problem 2. Given a set of LCs with integer variables, we would like to obtain a bound c such that |vol(relaxed(F))- #F | ≤ c.

A set of LCs also represents a convex polytope P. Similarly, let #P and vol(P) denote the number of integer points in P and the volume, respectively.

Problem 3. Given a convex polytope P, we would like to obtain a bound c that |vol(P) -#P | ≤ c.

Note that Problem 1, 2 and 3 are all equivalent. Without loss of generality, in this paper, we only use the definitions of problems on LCs with real variables and convex polytopes. In this section, we mainly discuss on convex polytopes. There are several more notations about polytopes.

• Let B(P) denote the boundary (facets) of P.
• Let mi(P) denote ⌈min{xi|x ∈ P}-1⌉. For example, min{xi|x ∈ P} = 1, then mi(P) = 0.
• Let Mi(P ) denote ⌊max{xi|x ∈ P} + 1⌋ similarly.
• Let Ph denote the intersection of P with a hyperplane h. For example, Px1=0 represents the cross-section of P with x1 = 0, i.e., Px1=0 = {Ax ≤ b, x1 = 0}.
• Let projh(P) denote the projection from P onto a hyperplane h.

Definition 2. An integer-cube is a unique unit-cube all of whose corners are integer points.

• Given an integer point α, we use c(α) to denote the integer-cube, where c(α) = {xi ∈ [αi, αi + 1], i = 1, . . . , n}.

• We use C(P) to denote the set of integer-cubes which intersects with P. For example, C(P) is the set of orange squares and grey squares in Figure 1. Similarly, we use C(B(P)) and C(projh(P)) to denote sets of integer-cubes that intersect with B(P) and projh(P) respectively.

3.2 Main Results

The difference between vol(P) and #P is intuitively caused by points that are close to the boundary of P . Moreover, these "controversial" points can be covered by integer-cubes that intersect with B(P), e.g., orange squares in Figure 1. So we draw the following theorem by this intuition.

Proof. Given an integer-cube c ∈ C(P) \ C(B(P)). Since c ∈ C(P), c has intersection with P. On the other hand, since c ̸∈ C(B(P)), c does not has intersection with the boundary B(P). So we know c ⊂ P, i.e., C(P) \ C(B(P)) is a set of integer-cubes inside P . Note that integercubes only share points on their boundaries and the volume of an integer-cube is 1. We have ∣∣C(P) \ C(B(P))∣∣ = ∑c∈C(P)\C(B(P)) vol(c) ≤ vol(P). We can map cubes in C(P) \ C(B(P)) to different integer points in P, e.g., pick a specific corner point of a cube. Therefore, the number of cubes ∣∣C(P) \ C(B(P))∣∣ is less or equal to #P. Note that ∣∣C(P ) \ C(B(P ))∣∣ = |C(P )| - ∣∣C(B(P ))∣∣. Then we obtain

In Figure 1, the size of the orange area is equal to ∣∣C(B(P))∣∣, which consists of integer-cubes in C(B(P)). Theorem 1 is the basis of theoretical guarantees, but lacks the bound of ∣∣C(B(P))∣∣. So it is our goal in the following analysis.

Definition 3. We call c(α){xn=αn}, which is a crosssection, the bottom of the integer-cube c(α), denoted as bottom(c(α)).

Definition 4. We call a polytope P cuts an integer-cube c if B(P) ∩ c ̸= ∅ and B(P) ∩ bottom(c) = ∅, i.e., B(P) intersects with c but does not intersects with bottom(c).

Lemma 1. Consider a set of integer-cubes Sn(α) = {c(αγ)| αγ = (α1, . . . , αn-1, γ), γ = ⌊mn⌋, . . . , ⌈Mn⌉}. Then P can cut at most two integer-cubes in Sn(α).

Proof. Assume there are three integer-cubes c1, c2, c3 ∈ Sn(α) such that B(P)∩ci ̸= ∅ and B(P)∩bottom(ci) = ∅. Let α1, α2 and α3 denote the points such that c(αi) = ci, i = 1, 2, 3. Note αis are same except the nth elements. Let vi represent the value of nth element in αi, i.e., αi =

1 Note that ∣∣C(B(P))∣∣ represents the cardinality of the set C(B(P )).

(α1, ..., αn-1, vi). Without loss of generality, let v1 < v2 < v3. Then c1 is above c2 and c2 is above c3.

Since B(P) ∩ c1 ̸= ∅ and B(P) ∩ c2 ̸= ∅, we can pick two points a ∈ B(P ) ∩ c1 and b ∈ B(P ) ∩ c2. The segment ab is inside P due to the convexity of P. In addition, ab ∩ bottom(c1) ̸= ∅ as c1 is above c2. So bottom(c1) ∩ P ̸= ∅. Furthermore, we have bottom(c1) ⊂ P, otherwise, there exists a point in bottom(c1) that also lies in B(P), which contradicts with B(P) ∩ bottom(c1) = ∅. Similarly, we can prove bottom(c2) ⊂ P , because B(P )∩c2 ̸= ∅, B(P )∩c3 ̸= ∅, and c2 is above c3.

From the assumption c2 ∩ B(P ) ̸= ∅, one can pick a point p ∈ c2∩B(P ), where p = (p1, . . . , pn). Consider two points p′ = (p1, . . . , pn-1, v1) and p′′ = (p1, . . . , pn-1, v2), obviously, p′ ∈ bottom(c1) and p′′ ∈ bottom(c2). Moreover, bottom(c1) ⊂ P and B(P) ∩ bottom(c1) = ∅, so we have p′ ∈ P \B(P ). Similarly, we have p′′ ∈ P \B(P ). From the convexity, the segment p′p′′ ⊂ P \B(P). We know that p is on segment p′p′′ , then p is not in the boundary of P, which is a contradiction.

Intuitively, we could obtain the bound of ∣∣C(B(P))∣∣ by mapping integer-cubes in C(B(P)) to integer-cubes which intersect with a bounding cuboid of P. For example, orange squares can be mapped to green squares in Figure 1. In detail, we prove the following theorem.

Theorem 2.

Proof. (Mathematical Induction)

Basis: Show that the Equation (1) holds for 1-dimensional problems. A 1-dimensional polytope is an interval. Its boundaries are the two end-points of this interval. Then ∣∣C(B(P ))∣∣ ≤ 2 and ∣∣C(proj{x=0}(P ))∣∣ = 1.

Inductive Step: Show that if Equation (1) holds for 1, 2, . . . , n - 1 dimensions, then it also holds for n dimensions.

Given a polytope P, we have the length of xi, i.e., [mi,Mi]. Then we subdivide P by constraints Xn(w) = {w ≤ xn ≤ w + 1}, where w = ⌊mn⌋, . . . , ⌈Mn⌉ - 1 are integers. It is easy to see that

where i = 1, . . . , n - 1.

Note that P{xn=w} is an (n - 1)-dimensional polytope. According to the induction hypothesis, there is

Since P{xn=w} ⊆ P ∩Xn(w), it is easy to see that

We know that projh(A ∩ B) ⊆ projh(A) ∩ projh(B) and proj{xi=0}(P) ∩ proj{xi=0}(Xn(w)) = proj{xi=0}(P) ∩ Xn(w), so we have

Then Equation (4) becomes

In order to have a better understanding of this proof, we also present a simple 3-dimensional example in Figure 2.

Then Fw is the number of integer-cubes that B(P) has intersections with, but no intersection with their bottoms, in C (B(P) ∩ Xn(w)). For example, in Figure 2, B(P{z=w}) are the 3 sides of △EFG and B(P) ∩ Z(w) are the 3 quadrangles EFIH, FGJI and EGJH. It is easy to see that an integer-square s ∈ C (B(P{z=w})) corresponds to a unique integer-cube c ∈ C (B(P) ∩ Z(w)) that s = (c ∩ (z = w)). Then s = bottom(c) and (s ∩ B(P)) ̸= ∅. It implies that bottom(c) intersects with B(P). On the other hand, we can also find that each integer-cube in C (B(P) ∩ Z(w)), whose bottom intersects with B(P ), corresponds to a unique integersquare in C (B(P{z=w})). So ∣∣C (B(P{z=w}))∣∣ is equal to the number of integer-cubes in ∣∣∣C (B(P) ∩ Z(w))∣∣∣ whose bottoms have intersection with B(P ).

Combining Equation (2), (3) and (5), we have

From Lemma 1, we can map each integer-cube, which B(P ) has intersections with, but no intersection with their bottoms, to a cube in C(proj{xn=0}(P)) by at most two times. So

Since computing the projection is not an easier problem, we would like to consider a simpler and practical version of Theorem 2, which is the following corollary.

Corollary 1.

We omit the proof of Corollary 1 which is plain, due to the limited space. This corollary provides a bound for ∣∣C(B(P ))∣∣ by volumes of surfaces of a bounding cuboid for P . To obtain the values of Mi(P) and mi(P), it requires at most 2n times of linear programming. In practice, the overhead is negligible compared with the cost of volume computation or estimation.

4 Algorithm

Based on Theorem 1 and Corollary 1, we present an approximation algorithm for integer solution counting with rigorous guarantees on LCs. The pseudo-code of our approach Vol2Lat is presented in Algorithm 1. The input is a convex polytope P (or a set of LCs) and the output is a triple that consists of the value, the lower bound and the upper bound of approximation.

We first apply Gaussian Elimination to reduce equalities, as they usually cause degeneration and lead to zero volume. Then we compute Mi and mi for Pk by linear programming and errk via the formula in Corollary 1. Since there are volume computation tools and volume estimation tools for convex polytopes, whose outputs are different, we present two procedures for them. For the volume v obtained by computation, the lower and upper bound are v - err and v + err respectively. Different from volume computation, the estimation routine, e.g., PolyVest, returns results in (ϵ, δ)-bound. It means that the estimation ˜ v of vol(P) lies in the interval [(1 + ϵ)-1 vol(P), (1 + ϵ)vol(P)] with at least probability (1 - δ). Such interval is also called the confidence interval. Let [Il, Iu] denote such (1 - δ) confidence interval. Then the bound is thus equal to [Il -errk, Iu +errk], which is also an (1 - δ) confidence interval.

Algorithm 1 Approximate Counting via Volume

Function Vol2Lat(P )

Apply Gaussian Elimination; Compute Mi(P) and mi(P), i = 1, . . . , n; err ← 2 n∑ i=1 ∏ i̸=j (⌈Mj (P )⌉ - ⌊mj (P )⌋); if enable volume computation then Call the volume computation routine to obtain the exact volume v of P ; return (v, v - err, v + err); else Call the volume estimation routine to obtain the estimation ˜ v with the confidence interval [Il, Iu]; return (˜ v, Il - err, Iu + err);

5 Evaluation

In order to evaluate our approach, we implemented a tool called Vol2Lat 2 . in C++ based on Algorithm 1 and a #SMT(LA) solver, which can handle LCs combined with propositional logic and also disequalities (̸=). We used a timeout of half an hour, and chose ϵ = 0.2 and δ = 0.1 for PolyVest (volume estimation). Experiments were conducted on a laptop workstation with 2.40GHz Intel Core i74700MQ CPU and 8GB memory. The suite of benchmarks consists of three families:

• Random Polytopes: For each LC, we first randomly pick a number l. Then we select l variables and generate coefficients (from -1000 to 1000) in uniform distribution. We filter out unsat cases and obtain 96 polytopes ranging from 2 to 15 dimensions at last.
• Multiple Knapsack: We generated 100 instances from the classic multiple knapsack problem [Ab´ıo and Stuckey, 2014]. The coefficients belong to [0, 10], the domain size is 64 and m = n, n ∈ {3, 4, 5}.
• Simple Temporal Planning: We randomly generated 150 instances based on the structure of Simple Temporal Networks (STNs) [Huang et al., 2018].
• Instances from program analysis: We analyzed 7 programs ('cubature', 'gjk', 'http-parser', 'muFFT', 'SimpleXML', 'tcas' and 'timeout') ranging from 0.4k to 7.7k lines of source code via a symbolic execution bugfinding tool called CAnalyze [Xu et al., 2014] 3 and generated 3803 instances in total.

Experimental Results

Figure 3 presents the average running time of three tools for each dimension. Note that the y-axis is in logarithmic coordinate. Table 2 presents more details, but only part of results, due to the limited space. #dims is the number of dimensions. L(F) and U(F) are the lower and upper bound computed by our approach. 'TIMEOUT' and 'MEMOUT' represent running out of time and memory (8GB) respectively.

2 Our tool and benchmarks are available at: https://github.com/ bearben/VOL2LAT

3 In order to analyze practical programs, bug-finding tools employ simplifications. For example, they unwind a loop by once (or several times), do not instantiate every element in an array, etc.

The results show that PolyVest and Vinci significantly outperform LattE. The running time of LattE exponentially increases as dimension increases. Vinci is very fast but suffers from tremendous memory consumption. However, the running time of PolyVest increases mildly. Due to the properties of the algorithms of LattE and Vinci, their running times are heavily affected by the structure of LCs. On some structural cases, LattE and Vinci can handle LCs with more dimensions and usually outperform PolyVest. However, since PolyVest employs a polynomial-time algorithm, it will eventually gain the upper hand. Besides the time-consumption, we can also observe that the values of vol(F) and ˜ vol(F) are very close to #F from Table 2. It is the key observation that motivates us to approximate integer counts via space volume. In addition, the bounds are usually close to #F. Recall that in Theorem 1, we compute bounds by subtracting or adding |C(B(P))| from#P. So, the lower bounds are worse and sometimes may be equal to zero.

To evaluate the quality of approximation bounds in applications, we conducted more experiments over constraints generated from multiple knapsack problems, STNs and program analysis. The average results for each family of benchmarks are listed in Table 3. Note that ¯ e, ¯ el and ¯ eu are the average values of e = |#F-vol(F)| #F , el = L(F) #F and eu = U(F) #F respectively. In general, the relative errors of approximations via volume are smaller than 10% and the bounds are often close to the exact counts #F, which are useful in many applications, such as counting-based search, knowledge representation, program analysis, etc. Since in these cases, we usually only need to find the polytope which has larger solution count instead of obtaining the exact count. For example, comparing the solution counts of branches to guide decision procedure during the search.

6 RelatedWorks

Kannan and Vempala [Kannan and Vempala, 1997] proposed an algorithm for sampling integer points in a convex polytope.

The proof of this algorithm shows that if a convex polytope P contains within it a ball of radius n√logm, then there is a constant c such that c · #P < vol(P) < 1c #P. Intuitively, the bound of c can be possibly approximated by inverting this sampling algorithm. However, there is no expression nor algorithm proposed to compute the constant c.

The Ehrhart Polynomial, to our best knowledge, is the result that is closest to our bounds. It says that #P can be expressed as a linear combination of the volume of the faces of P, i.e., #P = ∑F (vol(F) ∗ e(cone(P, F)), where F is a face of P, cone(P, F) is the supporting cone of P at F. Barvinok [Barvinok, 1994] proved that for integer polytopes, the coefficient e(cone(P, F )) can be computed using O(mk ) calls for the volume computation oracle, where m is the number of faces and F is k-dimensional. For rational polytopes, coefficients are considered to be more difficult to compute.

Dyer [Dyer et al., 1991] introduced a FPRAS algorithm for estimating the volume of convex bodies which is called Multiphase Monte-Carlo by its follow-up works. PolyVest is an implementation of such algorithm, which is restricted to polytopes for efficiency. But the goal of our work is different since we consider the integer solution counts instead of volume. Recall the discussion in Section 2, discretization will not work on integer counting as the length of grids cannot be arbitrarily small, although it is a natural way to handle continuous problems like volume.

7 Conclusion

Linear constraints appear frequently in many areas. Yet the associated counting problem has not received enough attention. We introduced an approximate counting algorithm via space quantification for LCs. It provides approximation with theoretical guarantees which are useful for users to make decisions of keeping or discarding approximations. Evaluation shows the scalability and high-efficiency of our approach. We also observe that the theoretical bounds are sometimes too loose compared with the practical errors. Therefore, improving bounds, especially the lower bound, is an interesting and also challenging direction of future research.

Acknowledgments

This work has been supported by Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS) under Grant QYZDJ-SSW-JSC036. Feifei Ma is also supported by the Youth Innovation Promotion Association, CAS. The authors would like to thank the anonymous reviewers for their comments and suggestions.