Space-Efficient Algorithm for Integer Programming with Few Constraints
Lars Rohwedder Karol Węgrzycki

A New First-Order Meta-Learning Algorithm
 with Convergence Guarantees
    
  El Mahdi Chayti

  Machine Learning and Optimization Laboratory (MLO), EPFL


  Martin Jaggi

  Machine Learning and Optimization Laboratory (MLO), EPFL

    September 9, 2024
================================================================================================================================================================

§ ABSTRACT 
  Integer linear programs min{c x : A x =
  b, x ∈ℤ^n_≥ 0}, where A ∈ℤ^m × n, b
  ∈ℤ^m, and c ∈ℤ^n, can be
  solved in pseudopolynomial time for any fixed number of constraints m = O(1).
  More precisely, in time (mΔ)^O(m)(I),
  where Δ is the maximum absolute
  value of an entry in A and I the input size.

  Known algorithms rely heavily on dynamic programming, which
  leads to a space complexity of similar order of magnitude as the running time.
  In this paper, we present a polynomial space algorithm that solves integer linear programs in
	(mΔ)^O(m (log m + loglogΔ))(I) time, that is,
  in almost the same time as previous dynamic programming algorithms.




§ INTRODUCTION


We consider integer linear programs (ILP) in equality form

    min{ c x : A x = b, x ∈^n }

where A ∈^m × n, b ∈^m, and c ∈^n. 
Integer linear programs are a widely used model to describe problems in combinatorial optimization and
the algorithms available for them are used as general purpose optimization suites.
Unlike its continuous counterpart, where x∈^n, integer linear programming is NP-complete even in the case of m=1. Hence, theoretical research is mainly focused on identifying structural assumptions, under which ILPs can still be
solved efficiently. Here, two dominant directions can be distinguished.

Based on the geometry of numbers, Lenstra <cit.> showed that ILPs of the form min{c x : A x ≤ b, x∈^n} can be solved in polynomial time for constant n, in fact, that they are fixed-parameter tractable (FPT) in parameter n. 
Note that this form is equivalent to (<ref>) by a standard transformation, but this would led to an increase in
parameter n.
Subsequent improvements by Kannan <cit.>
and most recently Reis and Rothvoss <cit.> lead to the currently best running time of (log n)^O(n)(I), where I is the input size.

The second direction is mainly based on dynamic programming approaches and goes
back to Papadimitriou <cit.>. He considered (<ref>) in the case when the absolute values of the entries of A and b are bounded by Δ and designed a pseudopolynomial time algorithm when m is fixed. This generalizes the classical pseudopolynomial time algorithm for Knapsack problems <cit.>. Later, his algorithm was improved by Eisenbrand and Weismantel <cit.>. The current best algorithm is due to Jansen and Rohwedder <cit.> with a running time of O(√(m)Δ)^2m(I).
In the latter two results the bound of Δ is only required on A, but not b.
Similar approaches were also applied successfully to constraint matrices A, where both dimensions n and m can
be large but satisfy certain structures of their non-zero entries, see e.g. <cit.>.



    Space efficiency of integer linear programs.

NP-completeness tells us that any algorithm that solves general ILPs needs to have an exponential running
time at least in some parameter unless P=NP. But not only running time, also space usage can be a limiting factor.
It is well known that for practical ILP solvers based on Branch-and-Bound it is often a bottleneck to store
the list of active nodes, see e.g. <cit.>.
This motivates the theoretical question of whether the aforementioned
results can also be achieved in polynomial space or, more generally, what tradeoffs between running
time and space one can achieve.

Frank and Tardos <cit.> showed that Lenstra's algorithm and the improvement by Kannan
with running time n^O(n)(I) can be implemented in polynomial space.
For the improvement by Reis and Rothvoss this is not known.

Regarding the algorithms based on dynamic programming, the space usage is also pseudopolynomial, which
is inherent in the method. Alternative, space-efficient methods are known for Subset Sum, a very
restricted variant of ILP, where the task is to find a solution to
{x ∈{0,1}^n : a x = b } with a ∈^n and b ∈. Note that here m=1 and there is
no objective function.
For the Subset Sum problem, Lokshtanov and Nederlof <cit.> introduced an algebraization technique and designed an Õ(n^3 b logb) time and Õ(n^2) space algorithm (see also <cit.>). This algorithm was later improved by Bringmann <cit.>, and the current best low-space algorithm for Subset Sum is due to Jin et al. <cit.>, which runs in Õ(nb) time and O(lognloglogn + logb) space.
Despite the intensive research on polynomial space algorithms for Subset Sum, little is known for
more general problems.
The seminal result of Lokshtanov and Nederlof <cit.> is the only one that considers the more general Knapsack problem, that is, max{c x : a x ≤ b, x∈{0,1}^n} with c, a∈^n, and b∈. Nevertheless, their method suffers from pseudopolynomial factor in both b and c_∞ in the running time, see <cit.>, whereas the classical dynamic programming algorithm only requires pseudopolynomial time in one of the two parameters. A (n, b) time and polynomial space algorithm is widely open for the Knapsack problem and Unbounded Knapsack (the variant with x∈^n) <cit.>.

Savitch's theorem can be used to avoid the high space usage of dynamic programming, but often comes at a significant
cost for the running time.
In the case of (<ref>), this method leads to an algorithm in time (mΔ)^O(m^2 log(mΔ))(I)
and polynomial space, see <Ref> for details.



    Our contribution.

In this paper, we show that integer programming (<ref>) can be solved in polynomial space
and a running time that nearly matches the best known running time of the dynamic programming approach by <cit.>.

	Let A ∈^m × n, b ∈^m, and c ∈^n and
	suppose that the absolute value of each entry of A is bounded by
	Δ > 1. Then, in time 
	
    (mΔ)^O(m(logm+ loglogΔ))·(I)

	and polynomial space, one can compute a solution
	to (<ref>) or conclude that such a solution does not exist.
	Here, I is the encoding size of A, b, and c.

Compared to the running time in <cit.>, this is an additional factor of
O(log m + loglogΔ) in the exponent of mΔ, whereas all known techniques lead to an additional factor of at
least O(mlogΔ), see <Ref>.



§ OVERVIEW OF PREVIOUS AND NEW TECHNIQUES
 
Let us briefly introduce some notation and the model of computation.
Throughout this paper, we assume all logarithms are base 2. We write Õ(f) = O(f ·(f)). Further, we make use of the shorthand notation [n] = {1, …, n}. For a vector x ∈^n, let x_p denote the ℓ_p norm of vector x. In the special case of p = 0, recall that x_0 = |(x)|, where (x) = { i ∈ [n] : x_i ≠ 0} are the non-zero coefficients.

For two vectors x,y ∈^n, we say that x is lexicographically smaller than y (x  y) if there exists an index i ∈ [n] such that
(i) x_j = y_j for every j < i, and (ii) x_i < y_j.

As in related works, see <cit.>, we work with a random-access model with word length Θ(logn + logW), where W is the largest absolute value of an input integer. We measure space complexity in terms of the total number of words in the working memory.
We say that an algorithm runs in polynomial space if the space complexity is bounded by a polynomial in the
input size, that is, in (I) := (n,m,log(Δ), log(b_∞), log(c_∞)).
Similarly, the running time of an algorithm is measured in the number of arithmetic operations
on words.

Before we explain the ideas behind <ref> we will briefly review
previous techniques used to achieve polynomial space dependency.

The standard tool to reduce space complexity usage are the ideas in Savitch's theorem, which transforms
algorithms with non-deterministic space complexity to slightly higher deterministic space complexity.
Here, the key technique is to decide (s, t)-reachability in a directed graph with N vertices in
N^O(logN) time and O(log N) space. This is straight-forward by recursively
guessing the middle vertex of the path from s to t. The space complexity comes from the O(log N)
vertices stored in the stack.
This can easily be adapted to finding the shortest
path in a weighted graph in the same time and space.
One can already apply this to ILP. Namely,
the algorithm of Eisenbrand and Weismantel <cit.>
constructs a directed acyclic graph with N = (m Δ)^O(m)·b_∞ vertices and
then computes the shortest path between two of the vertices. 
With standard polynomial time preprocessing described in <cit.>, one can reduce b_∞ to (mΔ)^O(m).
This shortest path query is the main source of the high space usage of <cit.> and
using Savitch's theorem, the ILP can therefore be solved in N^O(logN) =
(mΔ)^O(m^2 log(mΔ)) time and O(logN) = O(m log(mΔ)) space.

A more intricate technique is the algebraization technique <cit.>, which relies on the following coefficient test lemma.

	Consider a polynomial f(x) = ∑_i=0^d α_i x^i ∈𝔽[x] of degree at most d
	and assume there is an algorithm running in time T and space S for
	evaluating f(a) for a given a ∈𝔽. Then, there
	is an algorithm running in Õ(d T) time and Õ(S) space
	that for a given t ∈{0,…,d} decides if α_t = 0. 

The idea to solve the Subset Sum problem, is to consider the polynomial p(x) =
∏_i=1^n (1+x^a_i) and use <ref> to test whether the coefficient
α_t is nonzero, where a_1,…,a_n ∈ are the input integers,
t ∈ is the target, and α_t is a coefficient in front of monomial
x^t of polynomial p. Note that the polynomial p(x) can be evaluated in
(n, a_∞) time, which, combined with <ref> yields the result
of <cit.>. 

This technique can be applied to check feasibility of (<ref>)
by aggregating the m constraints of A into a single one a x = b' where
a_∞ = (mΔ)^O(m^2), see e.g. <cit.>,
and reducing the non-negative integer variables to binary ones,
which can be done by applying proximity techniques <cit.>.
We omit the details here, since
the method has two inherent disadvantages for our application:
it has a higher running time with an additional factor of O(m) in the exponent compared to <cit.> and
to use <ref> to solve
the optimization variant in (<ref>) one would need to additionally encode the
coefficient of the objective function c in a, which would yield a pseudopolynomial running time
in c_∞ that is not necessary in <cit.>.

Another straight-forward approach is based on the notable property of (<ref>) that it always has sparse optimal solutions.

	If (<ref>) is bounded and feasible, then there exists an optimal solution x^∗∈^n such that
	
		x^∗_0 ≤Γ,
	
	where Γ := 2(m+1)(log(m+1) + logΔ + 2). Moreover x^∗ is the lexicographically minimal solution.[Eisenbrand and Shmonin <cit.> did not state that x^∗ is the lexicographically minimal solution, but this can be derived easily: 
	one can replace the objective by c'_i = c_i + ϵ^i for each i∈ [n], where ϵ > 0 is sufficiently small. This simulates a lexicographical tie-breaking rule.]

This property is sometimes used in the following way: one can guess the support of an optimal solution
from the O(n^Γ) potential variable sets and then apply Lenstra's algorithm for integer linear programming
in small dimension. This idea was used for example in <cit.>.
Recall that Lenstra's algorithm and the improvement by Kannan can be implemented in polynomial space.

	There is an algorithm that solves
	
    max{c x : Ax ≤ b, x∈^n}

	in time n^O(n)(I) and polynomial space.

Note that (<ref>) can be translated to the form in <Ref> without increasing the number
of variables, hence the same result holds there.
Together with the bound on Γ, this yields
a running time of n^Γ·Γ^Γ·(I) ≤ (mΔ)^O(m^2 logΔ)·(I) and polynomial space, using that without loss of generality all columns of A are different and, in particular,
n ≤ (2Δ + 1)^m.

For our result we introduce two new ideas.
In <ref> we introduce the following branching algorithm that can be seen as an application of
Savitch's theorem to the support of the solution. Consider
the lexicographically minimal solution x^∗ and let x^∗ = x^∗_ℓ +
x^∗_r be such that x^∗_ℓ and x^∗_r have half the support each.
We will then guess a vector b_ℓ∈^m such
that Ax^∗_ℓ = b_ℓ and for each such guess, recurs into subproblems
with the right-hand sides b_ℓ and b - b_ℓ and then retrieve both x_ℓ^∗
and x_r^∗. Note, that in each recursive call, the support halves,
so the total depth of the recursion is O(logΓ) =
O(log m + loglogΔ). 
This idea already yields a (mΔ)^O(m^2 loglogΔ)(I) time and
polynomial space algorithm. However, the dependence on the number
of constraints m is much higher than what is promised in <ref>.
The bottleneck of the earlier mentioned idea of guessing the support and applying Lenstra's algorithm
is that one enumerates O(n^Γ) support sets.
In <ref> we show that the number of relevant support sets is much
lower and that they can be enumerated using the branching algorithm as a subroutine
on small instances where it runs more efficiently. This is then combined with running
Lenstra's algorithm on each support set.



§ BRANCHING ALGORITHM

As a subroutine in <ref>, we will first introduce an algorithm based on the idea
of branching on the support, which has the following properties.

	Let A ∈^m × n, b ∈^m, c ∈^n, and Δ≥A_∞.
	Assume that (<ref>) is bounded and feasible and let x^* be
	the lexicographically minimal optimal solution.

	For any σ∈ℕ we can in time (σΔ)^O(m logΓ)·(n)
	and polynomial space, compute x^* under the promise that x^*_1 ≤σ. If this promise is not satisfied, there is no guarantee of correctness.

Note that the running time matches <ref> except for the dependence on σ.
There are known bounds on h, see e.g., <cit.>, and
after the preprocessing in <cit.> one can guarantee a bound of σ = (mΔ)^O(m).
This would therefore give a polynomial space algorithm with running time (mΔ)^O(m^2logΓ)(I),
but this is larger than the running time claimed in <ref>. We will
apply it more efficiently in <Ref>.


Now, we present the algorithm behind <ref>. For the pseudocode, see <ref>.
This recursive algorithm has an additional parameter s, which is a bound on the support. 
On the initial call we use the bound Γ from <ref>.
The base case consists of the
scenario when s ≤ 1. If s = 0, we return the all-zero vector if and only
if b is an all-zero vector; otherwise, we return . If s = 1 let e_i, i
∈ [n], be the indicator vector, i.e., a vector with 1
at position i and 0 otherwise. We consider the set {z · e_i : z · A e_i = b, z
∈, i ∈ [n]} and return the last element
that minimizes the dot product with c (or  if no such element exists).
Observe that for each i, such an integer z can be determined in O(m)
time; hence, in total, the base case can be decided in O(mn) time.

When s > 1, we intuitively split the support of the solution into two (almost) equal-size parts
and guess the target b_ℓ for the first
part in the original solution. More precisely, we guess a
vector b_ℓ∈{-σΔ, …, σΔ}^m, recursively compute
x_ℓ and x_r, which are lexicographically minimal optimal solutions to A x_ℓ = b_ℓ and A x_r = b -
b_ℓ with half the support of the previous call (rounded up for one recursion and down for the other),
and if both of them exist,
return x_r + x_ℓ, breaking ties according to lexicographical order.
If all of the guesses for b_ℓ fail in at least one recursion, we return .
This concludes the description of the algorithm.


For the running time, observe
that the integer s decreases by a constant factor in each recursive call and starts with
	value Γ.
Thus, the depth of the recursion is O(logΓ). The number of possible b_ℓ
vectors at each recursion depth is (2 σΔ + 1)^m, so the total
size of the branching tree is (σΔ)^O(m logΓ).
The base case takes O(mn) time; hence, the total running time is as claimed.
	Similarly, the space complexity is (n, logΓ) as the algorithm needs to store the
stack of the recursion of height O(logΓ); the base case requires
	(n) space.

Next, we argue about the correctness. We prove by induction over s that
branch(A,b,c,σ,s) returns the lexicographically minimal optimal solution x^* ∈^n to
(<ref>), provided that x^*_1 ≤σ and x^*_0 ≤ s.

In the base case, s ≤ 1, the
algorithm straightforwardly checks over all possible non-zero coefficients of
x^* and returns a lexicographically minimal one. Note that if b is the all-zero vector then
assuming the integer program is bounded there cannot be a non-zero solution with negative objective
value. Hence, the all-zero solution is optimal and lexicographically minimal.

Now assume that s > 1. Let x^∗∈^n be the lexicographically
minimal solution to (<ref>) and assume that x^*_1 ≤σ and x^*_0 ≤ s. Let x_ℓ^∗, x_r^∗∈^n be
arbitrary vectors such that x_ℓ^∗ + x_r^∗ = x^∗, x_ℓ_0
≤x^∗_0/2, and x_r^∗_0 ≤x^∗_0/2.  For example, x_ℓ^∗ could be the first half
of the non-zero coefficients of x^∗, and x_r^∗ could be the remaining
non-zero coefficients.  Finally, let b^∗_ℓ = Ax^∗_ℓ. Observe
that x^*_ℓ_1 ≤x^*_1 ≤σ, x^*_r_1 ≤x^*_1 ≤σ, and b^∗_ℓ∈{-σΔ,…,σΔ}^m.
Furthermore, x^*_ℓ and x^*_r must be lexicographically minimal optimal solutions
for right-hand sides b^∗_ℓ and b - b^∗_ℓ since otherwise this would
contradict x^∗ being the lexicographically minimal optimal solution for b.

Hence, if in <ref> we guess b_ℓ = b_ℓ^∗, the induction hypothesis guarantees
that branch(A,b_ℓ,c,σ,s/2) and
ILP(A,b-b_ℓ,c,σ,s/2) return x^*_ℓ and x^*_r.
In particular, x^* = x^*_ℓ + x^*_r must be returned by the algorithm, since no
lexicographically smaller optimal solution exists that could dominate x^*.




§ CANDIDATE SUPPORT SETS AND MAIN THEOREM



A naive enumeration of all support sets of cardinality at most Γ would require
O(n^Γ). In this section we will design a more efficient enumeration and prove <Ref>.


	
	For a matrix A ∈^m × n, we say that a set 𝒞⊆{0,1}^n is a candidate support set if:
	
		
  *  x_1 ≤Γ for every x ∈𝒞, and
		
  *  for every b ∈^m and c ∈^n such that (<ref>) is bounded and feasible there exists an optimal solution x^∗ to (<ref>) with (x^∗) ∈𝒞.
	



	For every A ∈^m× n with A_∞≤Δ, there exists a candidate support set of size (m Δ)^O(m). Moreover, this set can be enumerated in time (m Δ)^O(m logΓ)·(I) and polynomial space.


We note that, with <ref>, the proof of <ref> is a straightforward application of known results.


	We enumerate a candidate support set 𝒞 of A using <ref>. Then for every c ∈𝒞, we consider the instance truncated to variables in (c) and solve it using the exact algorithm for integer programming by Lenstra and Kannan, i.e., <Ref>.

	The correctness follows from Property <ref> of 𝒞. For the running time, note that by <ref> we have |𝒞| ≤ (m Δ)^O(m) and it requires overhead of (m Δ)^O(m logΓ)(I) to enumerate 𝒞. Recall that the algorithm of Lenstra and Kannan <cit.> solves (<ref>) in Γ^O(Γ)(I) time and polynomial space for at most Γ variables. In total, the running time of the algorithm is
	
    (m Δ)^O(m logΓ)(I) + |𝒞| ·Γ^O(Γ)·(I)
    		   ≤ (mΔ)^O(m logΓ)·Γ^O(Γ)·(I) 
       ≤ (mΔ)^O(m (log m + loglogΔ))·(I).


It remains to prove <ref>.
This relies on the following property that Deng, Mao, and Zhong <cit.> noticed for Unbounded Knapsack
and which easily generalizes to ILP.

	Let x^* ∈^n be the lexicographically minimum solution to (<ref>). Then, there exists x' ∈{0,1}^n and b' ∈^m such that (i) x' is the lexicographically minimal optimal solution to min{c x' : A x' = b'} and (ii) (x') = (x^*).


	Let x' be the indicator vector of (x^*) and let b' = Ax'. Clearly, Property (ii) holds. Hence, it remains to prove (i). For the sake of contradiction, assume that y is a solution to the integer program with right-hand side b' and c y < c x'
Consider z := y + x^* - x'. Then z is a solution to {Ax = b', x∈^n}. But c z < c x^*, which contradicts the optimality of x^*.
Similarly, if there is a solution y with c y = c x' and y  x'. Then
z = y + x^* - x' is an optimal solution to min{c x, Ax = b', x∈^n} that is lexicographically smaller than x^*.


First, we iterate over every b' ∈^m such that b'_∞≤Γ·Δ. Note that the number of such b' is (2ΓΔ + 1)^m. Next, for each such b', we invoke <ref> with σ = Γ, which returns the lexicographically minimal optimal solution x^* to

    min{c x : Ax = b', x∈^n }

	if x^*_1≤σ.
	Finally, for each of these guesses, we output x if x ∈{0,1}^n. 
	This concludes the description of the algorithm. The algorithm runs in time
	
    (2ΓΔ+1)^m · (ΓΔ)^O(m logΓ)·(I) ≤ (m Δ)^O(m logΓ)·(I)

	and uses polynomial space.
	Property <ref> of <ref> holds because every x^* considered
	is lexicographically minimal and therefore by <ref> has support at most Γ. For Property <ref> consider some b∈^m, c∈^n for which (<ref>) is bounded and feasible. Let x^* be the lexicographically minimal optimal solution. By <ref> 
	we have x^*_0 ≤Γ. Furthermore, by <Ref> there is some x'∈{0,1}^n with (x') = (x^*), which is the lexicographically minimal optimal solution for b' = Ax'. Because b'_∞≤ΓΔ and x'_1 ≤σ the algorithm above will consider b' and correctly output x'.





§ CONCLUSION

In this paper, we show that an integer linear program can be solved
in time O(mΔ)^O(m(log m + loglogΔ)(I)
and polynomial space.
It remains open whether the additional factor of O(log m + loglogΔ) in the exponent
compared to the algorithm without space limitation can be avoided.
We note that our algorithm has two bottlenecks that lead to this running time:
 (i) the depth of O(log m + loglogΔ) for our branching and
(ii) because we use the algorithm of Lenstra and
Kannan that runs in time N^O(N)(I) on an instance with N = O(m log(mΔ)).
Thus, even if the improvement to (logN)^O(N)(I) by Reis and
Rothvoss <cit.> could be implemented in polynomial space, this would not immediately
lead to an improvement.
Already for the special case of Unbounded Knapsack it is unclear how to avoid
the exponent of loglogΔ and achieve a pseudopolynomial running time.

Finally, we want to mention the case of binary integer linear programs, i.e., where
we substitude the domain constraint for x∈{0,1}^n in (<ref>).
Here, the best known algorithm runs in time
(mΔ)^O(m^2)(I) <cit.> and naive application of
Savitch's theorem yields polynomial
space at significant higher running time of (mΔ)^O(m^3 log(mΔ))(I).
It is unclear how to come close to
an exponent of O(m^2) in polynomial space and the ideas in this paper do not seem to be
applicable, since the solutions need not be sparse.

plain


