N-fold integer programming in cubic time

N-fold integer programming is a fundamental problem with a variety of natural applications in operations research and statistics. Moreover, it is universal and provides a new, variable-dimension, parametrization of all of integer programming. The fastest algorithm for $n$-fold integer programming predating the present article runs in time $O(n^{g(A)}L)$ with $L$ the binary length of the numerical part of the input and $g(A)$ the so-called Graver complexity of the bimatrix $A$ defining the system. In this article we provide a drastic improvement and establish an algorithm which runs in time $O(n^3 L)$ having cubic dependency on $n$ regardless of the bimatrix $A$. Our algorithm can be extended to separable convex piecewise affine objectives as well, and also to systems defined by bimatrices with variable entries. Moreover, it can be used to define a hierarchy of approximations for any integer programming problem.


Introduction
N-fold integer programming is the following problem in variable dimension nt, where is an (r + ns) × nt matrix which is the n-fold product of a fixed (r, s) × t bimatrix A, that is, of a matrix A consisting of two blocks A 1 , A 2 , with A 1 its r × t submatrix consisting of the first r rows and A 2 its s × t submatrix consisting of the last s rows. This is a fundamental problem with a variety of natural applications in operations research and statistics which, along with extensions and variations, include multiindex and multicommodity transportation problems, privacy and disclosure control in statistical databases, and stochastic integer programming. We briefly discuss some of these applications in Section 5 (Corollaries 5.1 and 5.2). For more information see e.g. [4,5,7,9,10,11,12,14,15,20], [17,Chapters 4,5], and the references therein.
Moreover, n-fold integer programming is universal [6] and provides a new, variabledimension, parametrization of all of integer programming: every program is an n-fold program for some m over the bimatrix A := A(m) with first block A 1 the 3m × 3m identity matrix and second block A 2 the (3 + m) × 3m incidence matrix of the complete bipartite graph K 3,m . We make further discussion of this in Section 7.
The fastest algorithm for n-fold integer programming predating the present article is in [10] and runs in time O n g(A) L with L = w, b, l, u the binary length of the numerical part of the input, and g(A) the so-called Graver complexity of the bimatrix A. Unfortunately, the Graver complexity is typically very large [2,18]: for instance, the bimatrices A(m) mentioned above have Graver complexity g(A(m)) = Ω(2 m ), yielding polynomial but very large n Ω(2 m ) dependency of the running time on n.
In this article we provide a drastic improvement and establish an algorithm which runs in time O (n 3 L) having the cubic dependency on n which is alluded to in the title, regardless of the fixed bimatrix A. So the Graver complexity g(A) now drops down from the exponent of n to the constant multiplying n 3 . This is established in Section 3 (Theorem 3.9). Moreover, our construction can be used to define a natural hierarchy of approximations for (1) for the bimatrices A(m) with m variable, and therefore, by the universality theorem of [6], for any integer programming problem. These approximations are currently under study, implementation and testing, and will be discussed briefly in Section 7 and in more detail elsewhere.
Our algorithm extends, moreover, for certain nonlinear objective functions: using results of [16] on certain optimality criteria, we provide in Section 4 an optimality certification procedure for separable convex objectives whose time complexity is linear in n (Theorem 4.1) and an algorithm for solving problems with separable convex piecewise affine objectives whose time complexity is again cubic in n (Theorem 4.2). Furthermore, the algorithm also leads to the first polynomial time solution of n-fold integer programming problems over bimatrices with variable entries (Theorem 6.1).
So ⊑ extends the coordinate-wise partial order ≤ on the nonnegative orthant R n + to all of R n . By a classical lemma of Gordan, every subset Z ⊆ Z n has finitely-many ⊑-minimal elements, that is, x ∈ Z such that no other y ∈ Z satisfies y ⊑ x.
We have the following fundamental definition from [8].
Definition 2.1 The Graver basis of an integer m × n matrix A is defined to be the finite set G(A) ⊂ Z n of ⊑-minimal elements in {x ∈ Z n : Ax = 0, x = 0}.
For instance, the Graver basis of the matrix A := (1 2 1) consists of 8 vectors, Consider the general integer programming problem in standard form, A feasible step for feasible point x in (3) is any vector v such that x + v is also feasible, that is, Graver has shown that a feasible point x in (3) is optimal if and only if there is no element g ∈ G(A) in the Graver basis of A which is an augmenting step for x.
This suggests the following simple augmentation scheme: start from any feasible point in (3) and iteratively augment it to an optimal solution using Graver augmenting steps g ∈ G(A) as long as possible. While the number of iterations in this simple scheme as is may be exponential, it was recently shown in [10] that if in each iteration the best possible augmenting step of the form γg with γ positive integer and g ∈ G(A) is taken, then the number of iterations does become polynomial. In what follows, we call an augmenting step which is at least as good as the best possible augmenting step γg with g ∈ G(A), a Graver-best augmenting step.
It was shown in [5] that for fixed bimatrix A, the Graver basis G A (n) of the n-fold product of A can be computed in time polynomial in n. Thus, to find a Graver-best augmenting step of the form γg for an n-fold integer program (1), it is possible, as shown in [10], to check each element g ∈ G A (n) , and for each, find the best possible step size γ. However, as we explain below, the Graver basis G A (n) is very large. Therefore, in this article, we do it the other way around. For each of O(n) critical positive integer potential step sizes γ, we determine an augmenting step γh which is at least as good as the best possible augmenting step γg with g ∈ G A (n) . We then show that the best among these steps over all such γ is a Graver-best step.
In preparation for this, we need to review some material on Graver bases of nfold products. Let A be a fixed integer (r, s) × t bimatrix. For any n we write each vector x ∈ Z nt as a tuple x = (x 1 , . . . , x n ) of n bricks x i ∈ Z t . It has been shown in [1], [18], and [13], in increasing generality, that for every bimatrix A, the number of nonzero bricks appearing in any element in the Graver basis G A (n) for any n is bounded by a constant independent of n. So we can make the following definition.
Definition 2.2 The Graver complexity of an integer bimatrix A is defined to be the largest number g(A) of nonzero bricks g i in any element g ∈ G A (n) for any n.
This was used in [5] to show that the Graver basis G(A (n) ) has a polynomial number O(n g(A) ) of elements and is computable in time O(n g(A) ) polynomial in n. Thus, the computation of a Graver-best augmenting step in [10] was done by finding the best step size γ for each of these O(n g(A) ) elements of G(A (n) ), resulting in polynomial but very large O(n g(A) ) dependency of the running time on n. In Section 3 we show how to find a Graver-best step without constructing G(A (n) ) explicitly in quadratic time O(n 2 ) regardless of the bimatrix A and its Graver complexity.
We conclude this subsection with some remarks about complexity and finiteness. The binary length of an integer number z is the number of bits in its binary encoding, which is O(log |z|), and is denoted by z . The binary length z of an integer vector z is the sum of binary lengths of its entries. We denote by L the binary length of all numerical part of the input. In particular L = w, b, l, u for problem (1). All numbers manipulated by our algorithms remain polynomial in the binary length of the input and our algorithms are polynomial time in the Turing machine model. But we are mostly interested in the number of arithmetic operations performed (additions, multiplications, divisions, comparisons), so time in our complexity statements is the number of such operations as in the real arithmetic model of computation.
For simplicity of presentation we assume throughout that all entries of the bounds l, u in program (1) are finite and hence the set of feasible points in (1) is finite. This is no loss of generality since, as is well known, it is always possible to add suitable polynomial upper and lower bounds without excluding some optimal solution if any.

The algorithm
We now show how to decide if a given feasible point x = (x 1 , . . . , x n ) in (1) is optimal in linear time O(n), and if not, determine a Graver-best step γg for x in quadratic time O(n 2 ). This is then incorporated into an iterative algorithm for solving (1).
We begin with a lemma about elements of Graver bases of n-fold products.
Then for any n, any g ∈ G A (n) and any I ⊆ {1, . . . , n}, we have i∈I g i ∈ Z(A).
Proof. Consider any Graver basis element g ∈ G A (n) for some n. Then A (n) g = 0 and hence n i=1 A 1 g i = 0 and A 2 g i = 0 for all i. Therefore (see [17,Chapter 4]) each g i can be written as the sum g i = k i j=1 h i,j of some elements h i,j ∈ G(A 2 ) for all i, j. Let m := k 1 + · · · + k n and let h be the vector h := (h 1,1 , . . . , h 1,k 1 , . . . , h n,1 , . . . , h n,kn ) ∈ Z mt .
Then i,j A 1 h i,j = 0 and A 2 h i,j = 0 for all i, j and hence A (m) h = 0. We claim that moreover, h ∈ G A (m) . Suppose indirectly this is not the case. Then there is anh ∈ G A (m) withh ⊏ h. But then the vectorḡ ∈ Z nt defined byḡ i := k i j=1h i,j for all i satisfiesḡ ⊏ḡ contradicting g ∈ G A (n) . This proves the claim. Therefore, by Definition 2.2 of Graver complexity, the number of nonzero bricks h i,j of h is at most g(A). So for every I ⊆ {1, . . . , n}, we have that i∈I g i = i∈I k i j=1 h i,j is a sum of at most g(A) nonzero elements h i,j ∈ G(A 2 ) and hence i∈I g i ∈ Z(A).
Example 3.2 Let A := A(3) be the (9, 6) × 9 bimatrix mentioned in the introduction, which arises in the universality of n-fold integer programming discussed further in Section 7, having first block A 1 = I 9 the 9 × 9 identity matrix and second block the following 6 × 9 incidence matrix of the complete bipartite graph K 3,3 , Since A 2 is totally unimodular, its Graver basis G(A 2 ) consists of the 30 vectors in {0, ±1} 9 supported on circuits of K 3,3 with alternating ±1, see [17]. Also, it is known that the Graver complexity of this bimatrix is g(A) = 9, see [2,18]. Therefore, the set Z(A) in (4) which corresponds to A, consists of all sums of at most 9 such circuit vectors, and turns out to be comprised of 42931 vectors in Z 9 , such as We now define a dynamic program, that is, a weighted digraph, which will enable to find a Graver-best step γg for a feasible point x of (1) or detect that none exists.
is precisely the weight of the dipath h. The following lemma relates this dynamic program to Graver augmentations.
Proof. Let h be a minimum weight dipath from S 0 to S n and let g := g(h) be the vector associated with h. We claim that γg is the desired feasible step for x.
We begin with the complexity statement. Since A is fixed, so is g(A), and hence so is each S i . As the digraph is acyclic, the minimum weight dipath from S 0 to h i ∈ S i decomposes into a minimum weight dipath from S 0 to some h i−1 ∈ S i−1 plus the arc from h i−1 to h i . Thus, we have to check at most a constant number Repeating this for each of the sets S 1 , . . . , S n one after the other takes O(n) time.
We next show that γg is a good feasible step. Since dipath from S 0 to S n with weight wḡ = n i=1 w iḡi and associated vector g(h) =ḡ. Since h is a minimum weight dipath, wg ≤ wḡ and so w(x + γg) ≤ w(x + γḡ). The next lemma shows that we can quickly find a Graver-best augmentation. Lemma 3.6 A feasible step γg for x satisfying w(x + γg) ≤ w(x +γḡ) for any feasible stepγḡ withγ ∈ Z + andḡ ∈ G A (n) can be found in quadratic time O(n 2 ).
Proof. If Z(A) = {0} then G A (n) = ∅ by Lemma 3.1, so γg := 0 will do. Otherwise, construct a set Γ of O(n) positive integers in O(n) time as follows: for every i = 1, . . . , n and every z ∈ Z(A) \ {0} determine the largest positive integer γ such that l i ≤ x i + γz ≤ u i and include it in Γ. Now, for each γ ∈ Γ, construct and solve the corresponding dynamic program, resulting in total of O(n 2 ) time by Lemma 3.4. Let γg be that feasible step for x which attains minimum value wγg among the best steps obtained from all these dynamic programs, and letγḡ be that feasible step for x which attains minimum value wγḡ amongγḡ withḡ ∈ G A (n) if any. Assume that wγḡ < 0 as otherwise we are done since wγg ≤ wγ0 = 0. Then γ is the largest positive integer such that l ≤ x +γḡ ≤ u since otherwise the step (γ + 1)ḡ will be feasible and better. So for some i = 1, . . . , n, it must be thatγ is the largest positive integer such that l i ≤ x i +γḡ i ≤ u i . Sinceḡ ∈ G A (n) , it follows from Lemma 3.1 thatḡ i ∈ Z(A). Thereforeγ ∈ Γ. Now letγĝ be the best step attained from the dynamic program ofγ. Then wγg ≤ wγĝ by choice of γg and wγĝ ≤ wγḡ by Lemma 3.4. Therefore w(x + γg) ≤ w(x +γḡ) as claimed.
We next show, following [10], that repeatedly applying Graver-best augmenting steps, we can augment an initial feasible point for (1) to an optimal one efficiently.
Lemma 3.7 For any fixed bimatrix A there is an algorithm that, given n, w, b, l, u, and feasible point x for (1), finds an optimal solution x * for (1) in time O(n 3 L).
Proof. Iterate the following: find by the algorithm of Lemma 3.6 a Graver-best step γg for x; if it is augmenting then set x := x + γg and repeat, else x * := x is optimal.
To bound the number of iterations, following [10], note that while x is not optimal, and x * is some optimal solution, we have that x * − x = k i=1 γ i g i is a nonnegative integer combination of Graver basis elements g i ∈ G A (n) all lying in the same orthant, and hence each x + γ i g i is feasible in (1). Moreover, by the integer Carathéodory theorem of [3,19], we can assume that k ≤ 2(nt − 1). Letting γ i g i be a summand attaining minimum wγ i g i , and letting γg be a Graver-best augmenting step for x obtained from the algorithm of Lemma 3.6, we find that So the Graver-best step provides an improvement which is a constant fraction of the best possible improvement, and this can be shown to lead to a bound of O(nL) on the number of iterations to optimality, see [10] for more details. Since each iteration takes O(n 2 ) time by Lemma 3.6, the overall running time is O(n 3 L) as claimed.
We next show how to find an initial feasible point for (1) with the same complexity. We follow the approach of [5] using a suitable auxiliary n-fold program. Proof. Construct an auxiliary n-fold integer program min wz :Ā (n) z = b ,l ≤ z ≤ū , z ∈ Z n(t+2r+2s) (5) as follows. First, construct a new fixed (r, s) × (t + 2r + 2s) bimatrixĀ with We can now obtain the main result of this article.
Theorem 3.9 For every fixed integer (r, s) × t bimatrix A, there is an algorithm that, given n, vectors w, l, u ∈ Z nt and b ∈ Z r+ns having binary encoding length L := w, b, l, u, , solves in time O(n 3 L) the n-fold integer programming problem Proof. Use the algorithm of Lemma 3.8 to either detect infeasibility or obtain a feasible point and augment it by the algorithm of Lemma 3.7 to optimality.

Extensions to nonlinear objectives
Here we extend some of our results to programs with nonlinear objective functions, In [16] it was shown that Graver bases provide optimality certificates for problem (6) with separable convex functions as well: a feasible point x is optimal if and only if there is no feasible Graver step g for x which satisfies f (x + g) < f (x). This was used in [10] to provide polynomial time procedures for optimality certification and solution of problem (6) with separable convex functions f . However, this involved again checking each of the O(n g(A) ) elements of G A (n) .
Our results from Section 3 can be extended to provide linear time optimality certification for separable convex functions and a cubic time solution of (6) for separable convex piecewise affine functions. We discuss these respectively next.

Optimality certification for separable convex objectives
Here we assume that the objective function f is presented by a comparison oracle that, queried on two vectors x, y, asserts whether or not f (x) ≤ f (y). The time complexity now measures the number of arithmetic operations and oracle queries.
Theorem 4.1 For any fixed bimatrix A, there is an algorithm that, given n, b, l, u, separable convex f presented by comparison oracle, and feasible point x in program (6), either asserts that x is optimal or finds an augmenting step g for x which satisfies f (x + g) ≤ f (x +ḡ) for any feasible stepḡ ∈ G A (n) , in linear time O(n).
Proof. Given the feasible point x, set a dynamic program similar to that in Definition 3.3, with γ := 1, with the only modification that the weight of arc Then, for every dipath h and its associated vector g := g(h), we now have We now claim that the desired step is the vector g := g(h) associated with a minimum weight dipath h in this dynamic program. Indeed, an argument similar to that in the proof Lemma 3.4 now implies that for any feasible stepḡ ∈ G A (n) we have f (x + g) − f (x) ≤ f (x +ḡ) − f (x) and therefore f (x + g) ≤ f (x +ḡ).

Optimization of separable convex piecewise affine objectives
In [10] it was shown that problem (6) can be solved for any separable convex function in polynomial time, but with very large dependency of O(n g(A) ) of the running time on n, with the exponent g(A) depending on the bimatrix A. Here we restrict attention to separable convex objective functions which are piecewise affine, for which we are able to reduce the time dependency on n to O(n 3 ) independent of A.

So we assume again that
Moreover, we now also assume that for some fixed p, each f i j is p-piecewise affine, that is, the interval between the lower bound l i j and upper bound u i j is partitioned into at most p intervals with integer end-points, and the restriction of f i j to each interval k is an affine function w i j,k x i j + a i j,k with all w i j,k , a i j,k integers. We denote by f the binary length of f which is the sum of binary lengths of all interval end-points and w i j,k , a i j,k needed to describe it. The binary length of the input for the nonlinear problem (6) is now L := f, b, l, u .
Proof. We need to establish analogs of some of the lemmas of Section 3 for such objective functions. First, for the analog of Lemma 3.4, proceed as in the proof of Theorem 4.1 above: given a feasible point x and positive integer γ, set again a dynamic program similar to that in Definition 3.3, with the weight of arc A similar argument to that in the proof of Lemma 3.4 now shows that γg with g := g(h) the vector associated with the minimum weight dipath h is a feasible step for x which satisfies f (x+γg) ≤ f (x+γḡ) for any feasible step γḡ withḡ ∈ G A (n) .
For the analog of Lemma 3.6, we construct again a set Γ of critical step sizes γ as follows. First, as in the proof of Lemma 3.6, we collect the critical step sizes due to the lower and upper bound constraints by finding, for every i = 1, . . . , n and every z ∈ Z(A) \ {0}, the largest positive integer γ such that l i ≤ x i + γz ≤ u i , and include it in Γ. However, in contrast to Lemma 3.6, we are now dealing with the more general class of piecewise affine objective functions f i j . So we must add also the following values γ to Γ: for every i = 1, . . . , n, every z ∈ Z(A) \ {0}, and every j = 1, . . . , t, if x i j + γz j and x i j + (γ + 1)z j belong to different affine pieces of f i j , then γ is included in Γ. Since the number p of affine pieces is constant, the number of such values for each i, z and j is also constant. So the total number of elements of Γ remains linear and it can be constructed in linear time O(n) again. Now we continue as in the proof of Lemma 3.6: for each γ in Γ, using the analog of Lemma 3.4 established in the first paragraph above, we solve the corresponding dynamic program in O(n) time, resulting in total of O(n 2 ) time again. Let γg be that feasible step for x which attains minimum value f (x + γg) among the best steps obtained from all these dynamic programs, and letγḡ be that feasible step for x which attains minimum value f (x +γḡ) amongγḡ withḡ ∈ G A (n) if any. It now follows from the construction of Γ that ifγḡ is augmenting, namely, if f (x +γḡ) − f (x) < 0, thenγ ∈ Γ, as otherwise the step (γ + 1)ḡ will be feasible and better. Letγĝ be the best step attained from the dynamic program ofγ. Then where the first inequality follows from the choice of γg and the second inequality follows from the analog of Lemma 3.4. Therefore f (x + γg) ≤ f (x +γḡ).
For the analog of Lemma 3.7, we use the results of [10] incorporating the optimality criterion of [16], which assure that the number of iterations needed when using a Graver-best augmenting step at each iteration, is bounded by O(n f ) = O(nL), resulting again in overall time complexity O(n 3 L) for augmenting an initial feasible point to an optimal solution of (6). Since an initial feasible point if any can be found by Lemma 3.8 as before in the same complexity, the theorem now follows.

Some consequences
Here we briefly discuss two of the many consequences of n-fold integer programming which, now with our new algorithm, can be solved drastically faster than before.

Nonlinear multicommodity transportation
The multicommodity transportation problem seeks minimum cost routing of l commodities from m suppliers to n consumers subject to supply, consumption and capacity constraints. For l = 1 this is the classical transportation problem which is efficiently solvable by linear programming. But already for l = 2 it is NP-hard. Here we consider the problem with fixed (but arbitrary) number l of commodities, fixed (but arbitrary) number m of suppliers, and variable number n of consumers. This is natural in typical applications where few facilities serve many customers.
The data is as follows. Each supplier i has a supply vector s i ∈ Z l + with s i k its supply in commodity k. Each consumer j has a consumption vector c j ∈ Z l + with c j k its consumption in commodity k. The amount of commodity k to be routed from supplier i to consumer j is an integer decision variable x j i,k . The total amount l k=1 x j i,k of commodities routed on the channel from i to j should not exceed the channel capacity u i,j , and has cost f i,j l k=1 x j i,k for suitable univariate functions f i,j . We can handle standard linear costs as well as more realistic, convex piecewise affine cost functions f i,j , which account for channel congestion under heavy routing.
As a corollary of Theorem 4.2, for any fixed numbers l of commodities and m of suppliers, the problem can be solved in time cubic in the number n of consumers.
Corollary 5.1 For every fixed l commodities, m suppliers, and p, there exists an algorithm that, given n consumers, supplies and demands s i , c j ∈ Z l + , capacities u i,j ∈ Z + , and convex p-piecewise affine costs f i,j : Z → Z, solves in time O(n 3 L), with L := s i , c j , u i,j , f i,j , the integer multicommodity transportation problem Proof. Introduce new variables y j i and equations y j i = l k=1 x j i,k for all i, j. Then the objective function becomes i,j f i,j (y j i ) which is separable convex p-piecewise affine in the new variables, and the capacity constraints become y j i ≤ u i,j which are upper bounds on the new variables. Use u i,j as an upper bound on x j i,k and 0 as a trivial lower bound on y j i for all i, j, k. As shown in [11], arranging the original and new variables in a tuple z = (z 1 , . . . , z n ) of n bricks z j ∈ Z m(l+1) defined by

Privacy in statistical databases
A common practice in the disclosure of sensitive data contained in a multiway table is to release some of the table margins rather than the entries of the table. Once the margins are released, the security of any specific entry of the table is related to the set of possible values that can occur in that entry in all tables having the same margins as those of the source table in the database, see [7,20] and the references therein. In particular, if this set is small or consists of a unique value, that of the source table, then this entry can be exposed. Thus, it is desirable to compute the minimum and maximum integer values that can occur in an entry, which in particular are equal if and only if the entry value is unique, before margin disclosure is enabled.
Consider (d + 1)-way tables of format m 0 × · · · × m d , that is, arrays v = (v i 0 ,...,i d ) indexed by 1 ≤ i j ≤ m j for all j, with all entries v i 0 ,...,i d nonnegative integers. Our results hold for arbitrary hierarchical margins, but for simplicity we restrict attention to disclosure of d-margins, that is, the d + 1 many d-way tables (v i 0 ,...,i j−1 , * ,i j+1 ,...,i d ) obtained from v by collapsing one factor 0 ≤ j ≤ d at a time, with entries given by v i 0 ,...,i j−1 , * ,i j+1 ,...,i d := The problem is then to compute the minimum and maximum integer values that can occur in an entry subject to the margins of the source table in the database.
This problem is NP-hard already for 3-way tables of format n × m × 3, see [6]. However, as a corollary of Theorem 3.9, if only one side n of the table is variable, the problem can be solved in cubic time regardless of the other sides m i as follows. Proof. Let u be the maximum value of any entry in the given margins and use it as an upper bound on every variable. As shown in [5], this problem can be modeled as an n-fold integer programming problem, resulting in solution in large running time O n g(A) L , with exponent which depends on m 1 , . . . , m d . Theorem 3.9 now enables to solve it in cubic time independent of the table dimensions m 1 , . . . , m d .
We note that long tables, with one side much larger than the others, often arise in practical applications. For instance, in health statistical tables, the long factor may be the age of an individual, whereas other factors may be binary (yes-no) or ternary (subnormal, normal, and supnormal). Moreover, it is always possible to merge categories of factors, with the resulting coarser tables approximating the original ones, making the algorithm of Corollary 5.2 applicable.
We also note that, by repeatedly incrementing a lower bound and decrementing an upper bound on the entry x k 0 ,...,k d , and computing its new minimum and maximum values subject to these bounds, we can produce the entire set of values that can occur in that entry in time proportional to the number of such values.

Solvability over bimatrices with variable entries
The drop of the Graver complexity from the exponent of n to the constant multiple also leads to the first polynomial time solution of n-fold integer programming with variable bimatrices. Of course, by the universality of n-fold integer programming, the variability of the bimatrices must be limited. In what follows, we fix the dimensions r, s, t of the input bimatrix A, and let the entries vary. We show that, given as part of the input an upper bound a on the absolute value of every entry of A, we can solve the problem in time polynomial in a, that is, polynomial in the unary length of a. This holds for linear as well as separable convex piecewise affine objectives.
We have the following theorem, with L := f, a, b, l, u the length of the input. Theorem 6.1 For any fixed r, s, t and p, there is an algorithm that, given n, a, (r, s) × t bimatrix A with all entries bounded by a in absolute value, b, l, u, and separable convex p-piecewise affine f , in polynomial time O(a 3t(rs+st+r+s) n 3 L), solves Proof. Let G(A 2 ) be the Graver basis of the s × t second block A 2 of A, let p := |G(A 2 )| be its cardinality, and arrange its elements as the columns of a t × p matrix G 2 . Since r, s, t are fixed, it follows from bounds on Graver bases (see e.g. [17,Section 3.4]) that every g ∈ G(A 2 ) satisfies g ∞ = O(a s ) and hence p = O(a st ). Now, it is known (see [13,18] or [17,Section 4.1]) that the Graver complexity g(A) of A is equal to the maximum value v 1 of any element v in the Graver basis G(A 1 G 2 ) of the r ×p matrix A 1 G 2 . Since the entries of A 1 G 2 are bounded in absolute value by O(a s+1 ), the bounds on Graver bases (see again [17,Section 3.4]) imply that v 1 = O(p · (a s+1 ) r ) for every v ∈ G(A 1 G 2 ) and hence g(A) = O(a rs+st+r ). Now, consider again the following set defined in (4)  To find an initial feasible point, one can use the algorithm of Lemma 3.8, but then t would have to be replaced by t + 2r + 2s for the auxiliary bimatrixĀ, resulting in a somewhat larger exponent for a in the running time. However, it is possible to find an initial feasible point in an alternative, somewhat more involved way, keeping the original system with the bimatrix A, as follows. First find an integer solution to the system of equations only (without the lower and upper bounds) using the Hermite normal forms of the blocks A 1 and A 2 . Second, relax the bounds so as to make that point feasible. Third, minimize the following auxiliary objective function which is separable convex 3-piecewise affine, with If the optimal value is zero then the optimal auxiliary solution is feasible in the original problem, whereas if it is positive then the original problem is infeasible. Since this minimization can be done using the separable convex piecewise affine analog of Lemma 3.7 described in the proof of Theorem 4.2 in the same complexity O(|Z(A)| 3 n 3 L), the overall running time is O(a 3t(rs+st+r+s) n 3 L) as claimed.

Parametrization and approximation hierarchy
We conclude with a short discussion of the universality of n-fold integer programming and the resulting parametrization and simple approximation hierarchy for all of integer programming. As mentioned in the introduction, every integer program is an n-fold program for some m over the bimatrix A(m) having first block the identity matrix I 3m and second block the (3 + m) × 3m incidence matrix of K 3,m . It is convenient and illuminating to introduce also the following description.
Consider the following special form of the n-fold product operator. For an s × t matrix D, let D [n] := A (n) where A is the (t, s) × t bimatrix A with first block A 1 := I t the t × t identity matrix and second block A 2 := D. We consider such m-fold products of the 1 × 3 matrix (1 1 1). Note that (1 1 1) [m] is precisely the (3 + m) × 3m incidence matrix of the complete bipartite graph K 3,m . For instance, The following theorem was established in [6].
The Universality Theorem [6] Every (bounded) integer programming problem min{cy : y ∈ Z k + , V y = v} is polynomial time equivalent to some integer program This theorem provides a new, variable dimension, parametrization of integer programming: for each fixed value of the parameter m, the resulting programs above with n variable live in variable dimension 3mn and include natural models such as those described in Section 5, and can be solved in cubic time O(n 3 L) by Theorem 3.9; and when the parameter m varies, every integer program appears for some m.
Our new algorithm suggests a natural simple approximation hierarchy for integer programming, parameterized by degree d, as follows. Fix any d. Then given any m, let A := A(m), so t = 3m, A 1 = I 3m , and A 2 is the incidence matrix of K 3,m . Define the approximation at degree d of the set Z(A) in equation (4) Since A 2 is totally unimodular, G(A 2 ) consists of the O(m 3 ) vectors in {0, ±1} 3m supported on circuits of K 3,m with alternating ±1 and hence |Z d (m)| = O(m 3d ). Now, given a feasible point x in the universal program above, and positive integer γ, set a dynamic program similar to that in Definition 3.3, with the only modification that the sets S i are defined using the approximation Z d (m) of Z(A). Weaker forms of Lemmas 3.4 and 3.6 now assert that in time O(|Z d (m)| 3 n 2 ) = O(m 9d n 2 ), which is polynomial in both m and n, we can find a good feasible step γg. We use this iteratively to augment an initial feasible point to one which is as good as possible and output it. However, Lemma 3.1 no longer holds and not all brick sums of elements of the Graver basis G(A (n) ) lie in Z d (m). So the bounds on the number of iterations and total running time are no longer valid and the output point may be non optimal.
By increasing the degree d we can get better approximations at increasing running times, and when d = g(A) we get the true optimal solution. These approximations are currently under study, implementation and testing. They show promising behavior already at degree d = 3 and will be discussed in more detail elsewhere.
For m = 3, discussed in Example 3.2, for which the universal problem is equivalent to optimization over 3-way n × 3 × 3 tables, the approximation Z 3 (3) at degree d = 3 contains only 811 vectors out of the 42931 vectors in the true Z(A), such as −3 2 1 2 −3 1 1 1 −2 .