Counting Integral Points in Polytopes via Numerical Analysis of Contour Integration

In this paper, we address the problem of counting integer points in a rational polytope described by P(y) = { x ∈R^m Ax = y, x ≥ 0}, where A is an n × m integer matrix and y is an n-dimensional integer vector. We study the Z-transformation approach initiated by Brion-Vergne, Beck, and Lasserre-Zeron from the numerical analysis point of view, and obtain a new algorithm on this problem: If A is nonnegative, then the number of integer points in P(y) can be computed in O(poly (n,m, y_∞) (y_∞ + 1)^n) time and O(poly (n,m, y_∞)) space.This improves, in terms of space complexity, a naive DP algorithm with O((y_∞ + 1)^n)-size DP table. Our result is based on the standard error analysis to the numerical contour integration for the inverse Z-transform, and establish a new type of an inclusion-exclusion formula for integer points in P(y). We apply our result to hypergraph b-matching, and obtain a O(poly( n,m,b_∞) (b_∞ +1)^(1-1/k)n) time algorithm for counting b-matchings in a k-partite hypergraph with n vertices and m hyperedges. This result is viewed as a b-matching generalization of the classical result by Ryser for k=2 and its multipartite extension by Björklund-Husfeldt.

Authors

• 9 publications
• 1 publication
• 6 publications
• Counting vertices of integer polytopes defined by facets

We present a number of complexity results concerning the problem of coun...
05/04/2021 ∙ by Heng Guo, et al. ∙ 0

• On the number of integer points in translated and expanded polyhedra

We prove that the problem of minimizing the number of integer points inp...
05/09/2018 ∙ by Danny Nguyen, et al. ∙ 0

• How to Find the Convex Hull of All Integer Points in a Polyhedron?

We propose a cut-based algorithm for finding all vertices and all facets...
10/25/2020 ∙ by S. O. Semenov, et al. ∙ 0

• Fast Witness Counting

We study the witness-counting problem: given a set of vectors V in the d...
07/16/2018 ∙ by Peter Chini, et al. ∙ 0

• Efficient Distributed Computation of MIS and Generalized MIS in Linear Hypergraphs

Given a graph, a maximal independent set (MIS) is a maximal subset of pa...
05/09/2018 ∙ by Fabian Kuhn, et al. ∙ 0

• Partite Turán-densities for complete r-uniform hypergraphs on r+1 vertices

In this paper we investigate density conditions for finding a complete r...
03/11/2019 ∙ by Klas Markström, et al. ∙ 0

• Inapproximability of counting hypergraph colourings

Recent developments in approximate counting have made startling progress...
07/12/2021 ∙ by Andreas Galanis, et al. ∙ 0

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Counting integer points in polytopes is a fundamental problem. There are numerous applications in various areas of mathematical science, and fascinating mathematics behind; see e.g., [2, 4]. This problem is computationally intractable, i.e., it is P-hard [18]. Approximate counting as well as exact counting under fixed parameter settings has been rich sources for developments in the theory of algorithms and computational complexity. A seminal work by Barvinok [1] showed that there is a polynomial time algorithm to count integer points in rational polytope when the dimension of is fixed. His algorithm computes a certain “compact” expression (Brion-Lawrence formula) of the generation function for multivariate indeterminate . This (extremely difficult) algorithm is now implemented in computer package LattE [15], and provides a useful tool to the study of geometric combinatorics.

A “dual” generating-function approach was initiated by Brion-Vergne [8], Beck [3], and Lasserre-Zeron [13, 14]; see [12]. Suppose now that the input polytope is given by

 P(y)={x∈Rm:Ax=y,x≥0}, (1)

for integer matrix and -dimensional vector . Let , and consider the Z-transform . Brion-Vergne [8] showed that admits a very simple closed formula , and that the wanted is recovered by the inverse Z-transformation, which is a multi-dimensional contour integration of . This reduces the counting problem to the residue computation of . By this approach, Lasserre-Zeron [13, 14] developed an -time algorithm to count integral points in , where is a function of matrix .

In this paper, we study the contour integration of the inverse Z-transformation from the numerical analysis point of view, and obtain a new algorithm to count integer points for an important class of polytopes. Our main result is as follows.

Theorem 1.

Suppose that is nonnegative. For , the number of integer points in the polytope can be computed in time and space.

Notice that there is a simple DP algorithm with the same time complexity: For , consider the matrix consisting of the first columns of , and the number of integer points of polyhedron . By the nonnegativity of , integers for are obtained from for in time. The resulting DP algorithm, however, requires an -space for the DP table. Thus our result is regarded as an improvement in terms of space complexity.

Our technique for proving Theorem 1 is as follows. Instead of the residue computation of , we apply the numerical integration to the inverse Z-transform of , where we use the trapezoidal rule, a basic and popular method of numerical integration. By the standard error analysis of the trapezoidal rule [17]

, we obtain an error estimate with respect to the number

of sampled points of the numerical integration and contour radius . Interestingly this estimate gives rise to a new inclusion-exclusion type formula for , and brings the algorithm in Theorem 1, which is quite simple and is easier to be implemented. A notable feature of our inclusion-exclusion is to use the cancellation structure of trigonometric function in the complex plane . This extends the usual inclusion-exclusion based on the cancellation of and in . Our algorithm computes the number of integer points in the expression for , and recovers the “true” value by algebraic computation on the group ring of cyclic group , which avoids numerical computation of .

Our result is applicable to packing-type polytopes, which ubiquitously arise from graph theory and combinatorial mathematics. Consider the particular case where is 0-1 valued and is the all-one vector . Then is the number of perfect matchings in the hypergraph corresponding to . Exact exponential time counting algorithms for matchings have been intensively studied in recent years [5, 6, 9, 11]. The formula of hypergraph matching derived from our result is viewed as a variant of that given by Björklund-Husfeldt [6]. The time complexity matches that of their algorithm. By exploiting special properties, Björklund-Husfeldt [6] improve the time complexity to for -partite hypergraphs. This result is viewed as a multipartite extension of the classical result of Ryser [16] for bipartite matching. We also exploit a special structure of our formula in -partite hypergraphs, and prove the following -matching generalization of the results of Ryser and Björklund-Husfeldt.

Theorem 2.

Let be a -partite hypergraph, and let . The number of perfect -matchings in can be computed in time and space.

Counting perfect -matchings of a -partite graph (i.e., -dimensional -matchings) has many applications in a wide range of mathematical sciences that include combinatorics, representation theory, and statistics; see e.g., [10]

. For example, counting multiway contingency tables with prescribed margins, an important problem for statistical analysis on contingency tables, is nothing but

-dimensional -matching counting.

The rest of this paper is organized as follows. In Section 2, we set up basic notation, and introduce Z-transformation, its inverse, and approximate inverse Z-transformation obtained by numerical integration. In Section 3, we present our algorithm to prove the main theorem. In Section 4, we discuss hypergraph matching and prove Theorem 2 in a further generalized form.

2 Preliminaries

2.1 Notation

Let denote the set of positive real numbers. Let denote the set of nonnegative integers. For a matrix , let denote the -th column vector of the matrix. For an integer vector and a complex vector , define by

 zy\coloneqqzy11⋯zynn.

For a positive integer , define by

 ωN(h)\coloneqqexp2πiNh(h∈{0,1,2,…,N−1}). (2)

The following relation is well-known:

 N−1∑j=0ωN(kj)={Nif k=0modN,0otherwise,(k∈Z). (3)

For a function and , let denote , such as

 exp(z) =(exp(z1),exp(z2),…,exp(zn))⊤, (4) lnz =(lnz1,…,lnzn)⊤, (5) ωN(z) =(ωN(z1),…,ωN(zn))⊤. (6)

The symbol is meant as “undefined.” We use when the function value is defined via integration or infinite summation, possibly not converging.

2.2 Z-transformation

For a function , define the Z-transform of by

 ^f(z)\coloneqq∑y∈Znf(y)z−y(z∈Cn). (7)

The inverse of the Z-transformation is given as follows. For a function and , define by

 Ir[g](y) \coloneqq1(2πi)n∮|z1|=r…∮|zn|=rg(z)zy−1dz1⋯dzn (8) =∫[0,1)…∫[0,1)g(rexp(2πit))r1⊤yexp(2πit⊤y)dt1⋯dtn(y∈Zn), (9)

where we change variables by in (9). Under an appropriate condition on and , map is actually the inverse of the Z-transformation:

 Ir[^f]=f. (10)

We do not go into details under which conditions (10) holds. Instead, we consider an approximate inverse Z-transform by the numerical integration applied to (9). Here we use the trapezoidal rule, which is a basic and popular method of numerical integration; see e.g., [7]. For a positive integer (the number of points in the numerical integration), define by

 IN,r[g](y)\coloneqq1Nn∑j∈{0,1,…,N−1}ng(rωN(j))r1⊤yωN(j⊤y)(y∈Zn). (11)

Recall notation . Our counting algorithm is based on .

3 Counting integral points in a polytope

Let be an integral matrix. We assume that there is no nonzero nonnegative vector with . This assumption ensures that the polytope is bounded for every . In the case where is a nonnegative matrix, this assumption is equivalent to the property that each column of has at least one nonzero entry.

As mentioned in the introduction, define function by

 fA(y)\coloneqq|{x∈Zm:Ax=y,x≥0}|(y∈Zn).

Our starting point is the following formula for .

Theorem 3 ([3, 8, 13]).
1. For with , the Z-transform is given by

 ^fA(z)=∑h∈Zm≥0z−Ah=m∏k=111−z−Ak, (12)

where the series absolutely converges.

2. For with , it holds

 fA(y)=1(2πi)n∮|z1|=s1⋯∮|zn|=sn^fA(z)zy−1dz1⋯dzn. (13)

We establish an approximate version of the above theorem as follows.

Theorem 4.

Suppose that is nonnegative. For , , and , it holds

 IN,r[^fA](y) (14) =fA(y)+∞∑k=1r−Nk|{x∈Zm≥0:y−Ax≥0,1⊤(y−Ax)=Nk}|. (15)
Proof.

From the assumption that is nonnegative and each column of has nonzero entry, for with , it holds . By the previous theorem, the Z-transform is given by

 ^fA(z)=∑h∈Zm≥0z−Ah. (16)

Substituting this expression to (in (11)), we have

 IN,r[^fA](y) =1Nn∑j∈{0,1,…,N−1}n∑h∈Zm≥0r1⊤(y−Ah)ωN(j⊤(y−Ah)) (17) =1Nn∑h∈Zm≥0r1⊤(y−Ah)∑j∈{0,1,…,N−1}nωN(j⊤(y−Ah)) (18) =1Nn∑h∈Zn≥0r1⊤(y−Ah)m∏l=1N−1∑jl=0ωN(jl(yl−(Ah)l)), (19)

where the summations are interchangeable, thanks to the absolute convergence (by ). By the relation (3), is if , and zero otherwise. Notice that cannot occur since . Thus we have

 IN,r[^fA](y)=∑u∈Zn≥0r−N1⊤u|{x∈Zm≥0:y−Ax=Nu}|.

Gathering with , we obtain (15). ∎

This proof is inspired by the standard argument to derive the exponential convergence of the trapezoidal rule applied to periodic functions, where the the cancellation technique using (3) in the proof is known as aliasing in numerical analysis; see [17, Theorem 2.1].

Notice that is a variant of the Ehrhart (quasi)polynomial, and is bounded by a polynomial in with degree . Hence the series in (15) actually absolutely converges for .

Corollary 5.

Suppose that is nonnegative. Let and . Then is equal to the coefficient of in

 1Nn∑j∈{0,1,…,N−1}n∑h∈{0,1,2,…,N−1}mr−1⊤AhωN(j⊤(y−Ah)) (20) =1Nn∑j∈{0,1,…,N−1}nωN(j⊤y)m∏l=1N−1∑hl=0r−1⊤AlhlωN(−j⊤Alhl). (21)

Here we regard as an indeterminate.

Proof.

In the formula (14), if has , then for some , and does not contribute to the constant term . The claim follows from this fact and the observation . ∎

Our goal is to compute the coefficient of in (21). Instead of numerical computation of the trigonometric function , we develop an algebraic algorithm. Let denote the group ring of the cyclic group of order . Namely, consists of polynomials with variable , rational coefficients, and degree at most , in which the multiplication rule is given by . Consider the bivariate polynomial ring with variables . Then is a ring homomorphism from to . Letting and in (20), we obtain a polynomial in :

 1Nn∑j∈{0,1,…,N−1}n∑h∈{0,1,…,N−1}mt1⊤Ahsj⊤(y−Ah). (22)

Let denote the coefficient of in (22). Then it holds

 ¯¯¯¯¯¯fA(ωN(1))=fA(y). (23)

Our algorithm first computes , and then computes .

Lemma 6.

can be computed in time and space.

Proof.

Let . From (21), we first consider the computation of the coefficient of in

 m∏l=1N−1∑hl=0t1⊤Alhls−j⊤Alhl (24)

for fixed . It suffices to compute the above polynomial (24) modulo , which can be written as

 m∏l=1skl0+skl1t1+⋯+skldtdmod(td+1). (25)

for some . The integers are obtained in time by computing and for . Expand (25) to the form , and obtain . This computation is done by the multiplication of polynomials with degree modulo , where their coefficients are polynomials with degree in multiplication rule . Now is the sum of over divided by . Thus is obtained in arithmetic operations (over ).

Finally we estimate the bit-size required for the computation. It suffices to estimate the size of , which is the sum of at most terms of form . Then the coefficients of have bit-length . Thus the required bit-size is at most . ∎

Next we consider how to compute from .

Lemma 7.

is written as

 ¯¯¯¯¯¯fA(s)=fA(y)+∑iKi(1+si+s2i+⋯+sN−i).

where the sum is taken over divisors of with some coefficient .

Proof.

Regard as . Then the map defined by is a group homomorphism. Therefore the image of is the cyclic group for some divisor of . Also the number of inverse images of each is given by . Then, for with , it holds

 ∑j∈{0,1,…,N−1}nt1⊤Ahsj⊤(y−Ah) (26) =t1⊤y∑j∈{0,1,…,N−1}nsj⊤(y−Ah)=t1⊤yJh(1+snh+s2nh+⋯+sN−nh). (27)

Thus we have

 ¯¯¯¯¯¯fA(s) =∑h∈{0,1,…,N−1}mJh(1+snh+s2nh+⋯+sN−nh) (28) =∑i:divisor of NKi(1+si+s2i+⋯+sN−i), (29)

where is the sum of over such that the image of is . Notice . ∎

According to this lemma, we obtain a simple algorithm to compute from as follows.

• Let .

• If for all , then output ; stop

• Choose the minimum index with . Let for each index that is the multiple of the index , and go to step 1.

The correctness of the algorithm is clear from the above lemma: The chosen index in step 2 is a divisor of with . Hence the algorithm computes . After at most iterations, the algorithm terminates and outputs the correct answer .

Example 8.

Consider the following matrix and vector :

 A=(113111),y=(53). (30)

Then the polytope has three integer points:

 ⎛⎜⎝111⎞⎟⎠,⎛⎜⎝201⎞⎟⎠,⎛⎜⎝021⎞⎟⎠. (31)

Let us count the integer points according to Corollary 5. The coefficient of in

 1625∑j1=05∑j2=0s5j1+3j3⎛⎝5∑h1=0t2h1s(−j1−j2)h1⎞⎠2⎛⎝5∑h2=0t4h2s(−3j1−j2)h2⎞⎠ (32)

is

 1625∑j1=05∑j2=0s5j1+3j3(5s−4j1−4j2+3s−5j1−3j2+s−6j1−j2) (33) = 1365∑j1=05∑j2=0(5sj1−j2+3+s−j1+2j2) (34) = 136(120+12s+12s2+12s3+12s4+12s5)=3+13(1+s+s2+s3+s4+s5). (35)

Therefore we obtain .

Example 9.

Consider the following matrix and vector :

 A=(1221),y=(75) (36)

Then the polytope has only one integer point:

 (32). (37)

Consider the coefficient of in

 1827∑j1=07∑j2=0s7j1+5j2⎛⎝7∑h1=0t3h1s(−2j1−j2)h1⎞⎠⎛⎝7∑h2=0t3h2s(−j1−2j2)h2⎞⎠. (38)

By calculation, is

 1827∑j1=07∑j2=0s7j1+5j2(s−8j1−4j2+s−7j1−5j2+s−6j1−6j2+s−5j1−7j2+s−4j1−8j2) (39) = 1647∑j1=07∑j2=0(s−j1+j2+1+sj1−j2+s2j1−2j2+s3j1−3j2) (40) = 164(104+24s+40s2+24s3+40s4+24s5+40s6+24s7) (41) = (42)

Thus we obtain .

4 Hypergraph Matching

We next show Theorem 2 in a generalized form. Let be an nonnegative integer matrix. For each column index , consider subset consisting of row indices with . Let denote the hypergraph on vertex set and hyperedge set . By a stable set of we mean a vertex subset such that every hyperedge meets at most one vertex in .

Theorem 10.

Suppose that we are given a stable set of . For , we can compute in time and space.

Proof.

Let , , and . As before, it suffices to compute (22) modulo . We show that (22) admits the following factorization:

 ∑j∑ht1⊤Ahsj⊤(y−Ah)=∑j′∈{0,1,…,N−1}n−νF0(j′)F1(j′)⋯Fν(j′), (43)

where each is computable modulo in time and space.

By arranging indices of , we can assume that , and that is regarded as a block matrix , where consists of columns such that the corresponding hyperedge meets . Accordingly, vector is also partitioned as

 h=⎛⎜ ⎜ ⎜ ⎜⎝h0h1⋮hν⎞⎟ ⎟ ⎟ ⎟⎠,Ah=A0h0+A1h1+⋯+Aνhν.

We suppose that an -dimensional vector is embedded to by filling to the first components. Each is uniquely represented as for , where is the -th unit vector. Then we have

 j⊤y = j′⊤y+j1y1+⋯+j