A Simple Algorithm for Combinatorial n-Fold ILPs Using the Steinitz Lemma

The power of the linear algebra method has a long and storied history in the design and analysis of algorithms, [3, 32]. The question of reordering vectors so that all partial sums are bounded was posed by mathematicians Paul Lévy and George Riemann in the early 20th century. The question was first affirmatively answered by Ernst Steinitz [33, 15] showing that such a reordering of vectors, that sum to zero and are all in a unit ball with respect to any arbitrary norm, is indeed possible where all the partial sums can be bounded by a constant that depends only on $d$, the dimension of the vectors. This result, now eponymously known as the Steinitz Lemma, stated in Proposition 2, has become a cornerstone in our understanding of the geometry of numbers and has led to numerous algorithmic applications centered around problems that have an integer linear program (ILP) formulation [12]. In this paper, we study a simple Steinitz Lemma based non-augmenting algorithm for a class of combinatorial $n$-fold ILPs.

The usefulness of Steinitz Lemma in modern algorithmic research is perhaps best underscored by its applicability to analyze integer linear programs (ILP), as exhibited by Eisenbrand and Weismantel [12] who improved upon the longstanding best bound of Papadimitriou [29]. Papadimitriou’s seminal work gave a simpler proof of ILPs being in NP and the first pseudo-polynomial algorithm when the number of constraints is constant. This can be viewed as a pseudo-parameterized algorithm with respect to the number of constraints. Subsequently, Jansen and Rohwedder [20] further improved the running time by applying the Steinitz Lemma in an alternative way. In this paper, we are continuing this line of investigation by employing Steinitz Lemma on a special class of ILPs, namely combinatorial $n$-fold ILP, for the purpose of obtaining better bounds.

Notwithstanding these exciting possibilities, the bottom line is that solving an ILP is known to be NP-complete and that poses a computational challenge when we desire integral solutions only. While we know of certain classes of ILPs, such as totally unimodular, bimodular [2], binet [1] that can be solved in polynomial time, the task of mapping the ILP landscape into families that are tractable under modest restrictions is ongoing. In particular, in recent times researchers are motivated to design parameterized algorithms because they behave like polynomial-time algorithms when the parameter is constant-valued. An algorithm is said to be fixed-parameter tractable (FPT) with respect to parameters $k_1,\mathrm{\dots },k_{\mathrm{\ell }}$, if its running time is of the form $f(k_1,\mathrm{\dots },k_{\mathrm{\ell }})\cdot {|\text{input}|}^{𝒪(1)}$ for some computable function $f$. Steinitz Lemma has proven to be useful even in this direction of research, as exhibited in its usefulness in the study of block-structured ILPs.

Suppose that the given program (1.1) admits a solution $x$. From $x$, we obtain a sequence of vectors $v_1,\mathrm{\dots },v_q$ such that

1. 1.

   $q={\sum }_{i\in [n]}{b}^{(i)}$;

2. 2.

   each $v_j$ is a column of one of the top blocks $T^{(1)},\mathrm{\dots },T^{(n)}$;

3. 3.

   ${\sum }_{j\in [q]}{v}_j=b^{(0)};$

4. 4.

   the column $T^{(i)}[\cdot ,k]$ appears $x_k^{(i)}$ many times in the sequence.

Consider the $1\times q$ matrix $M$ with entries from $[n]$ where $M[j]=i$ if $v_j=T^{(i)}[\cdot ,k]$ for some $k\in [t]$, and each $i\in [n]$ is present exactly $b^{(i)}$ times in $M$. Informally, through $M$ we identify the block from which each vector is obtained. The first step of our approach is to “balance” this one-row matrix $M$.

Suppose that the symbol (or element) $e\in [n]$ appears $m_e$ times in $M$. Consider a random rearrangement of the matrix $M$. The expected number of occurrences of $e$ in the first $j$ columns, by the linearity of expectation, is $(j/q)\cdot m_e$. This is because the probability that a random rearrangement of symbols in $M$ contains $e$ in a particular index is $m_e/q$, and the linearity of expectation is applied over the first $j$ columns. Let ${\text{occ}}_M(e,j)$ denote the number of occurrences of symbol $e$ between columns $1$ to $j$ in $M$. The imbalance of a symbol $e$ at a column $j$ in $M$ is defined as

Thus, the imbalance of a symbol $e$ at a particular column is the difference between the number of occurrences of $e$ and the expected number of occurrences of $e$ until that column in a random rearrangement of the entries. We show that we can, in time $𝒪(nq)$, rearrange the symbols such that the matrix has bounded imbalance¹¹1We would like to note that a similar result can be obtained by applying the Steinitz Lemma on the $q\ast n$ matrix $M^{\prime }$ where for each $j\in [q]$, $M^{\prime }[\cdot ,j]$ is a vector containing $1$ in position $M[j]$ and $0$ in all other positions. However, such an application will result in a weaker imbalance upper-bound of $n$ and a slightly more complicated algorithm..

For each block $i\in [n]$, we do the following. Consider the $r\times b^{(i)}$ matrix formed by all the $b^{(i)}$-many vectors in the sequence corresponding to the block $i$, say $V^{(i)}$. By Proposition 3, there is a permutation of the vectors in $V^{(i)}$, say ${\pi }^{(i)}\in S_{b^{(i)}}$, such that each coordinate of the $j^{th}$ partial sum is at most $2r\mathrm{\Delta }$ away from the ${(j/b^{(i)})}^{th}$ fraction of the total sum of vectors in $V^{(i)}$. We construct a new $r\times q$ matrix $O$ from $M_{\sigma }$ by replacing the entries with vectors from those blocks. We place these vectors $V_{{\pi }^{(i)}}^{(i)}$ (in the same order) in $O$ based on where $i$ occurs in $M_{\sigma }$: the $w^{th}$ vector in $V_{{\pi }_{(i)}}^{(i)}$ is placed in the $w^{th}$ occurrence of $i$ in $M_{\sigma }$. We show that all partial sums in $O$ are bounded.

We construct a weighted directed acyclic graph $G$ based on the bounds we obtained in Lemma 5. The key property of $G$ is that any shortest path between two particular vertices in $G$ corresponds to an optimal solution of the program and vice versa. If there is no such path, then we assert that there is no solution. Otherwise, we compute an optimal solution corresponding to the path. Thus, we have our main result, Theorem 1.

In this section, we prove Lemma 4 in Section 3.1, followed by Lemma 5 in Section 3.2, and then present the overall algorithm in Section 3.3 which forms a proof of Theorem 1.

In this section, we generalize our main result by showing that programs which contain “certain” inequalities in the constraints can be solved in the same time that we take for (1.1) where we have $𝒜x=b$. This is towards the ease of expressing problems involving inequalities as $n$-fold integer programs which we will see later in Section 5.

We call the upper rows $(T^{(1)}{T}^{(2)}\mathrm{\cdots }{T}^{(n)})x=b^{(0)}$ as globally uniform constraints, and the lower rows $1^{⊺}{x}^{(i)}=b^{(i)}$ as locally uniform constraints. Consider the generalization of (1.1)

1. 1.

   where the globally uniform constraints may contain either of $\{\le ,=,\ge \}$;

2. 2.

   where the locally uniform constraints may contain either of $\{\le ,=\}$.

Let $A=(1,\mathrm{\dots },1)\in {\mathbb{Z}}^{1\times t}$, $T^{(1)},\mathrm{\dots },T^{(n)}\in {\mathbb{Z}}^{r\times t}$, $b\in {\mathbb{Z}}^{r+n}$, and $c\in {\mathbb{Z}}^{nt\times 1}$. Formally, we consider the following generalization.

We emphasize that the generalization does not allow locally uniform constraints to contain “$\ge$”. We would like to note that our approach at a high level is almost the same as that of [25], and the main difference lies in the fact that (unlike as in their model,) here the variables in the programs we consider are non-negative and have no explicit lower or upper bounds. The proof strategy is to construct an appropriate program of the form (1.1) by adding new (slack) variables. Let $\psi ≔2+2{\parallel c\parallel }_{\mathrm{\infty }}\cdot {\max }_{i\in [n]}\{b^{(i)}\}$. Observe that $\psi /2$ is strictly larger than the maximum value the objective function in (${𝐏}_{\mathrm{𝟐}}$) can take.

First, we deal with the locally uniform constraints in (${𝐏}_{\mathrm{𝟐}}$). We add $n$ variables: for each $i\in [n]$, we add $x_{t+1}^{(i)}$ and replace $T^{(i)}$ by $(T^{(i)}\mathbf{0})$ where $\mathrm{𝟎}$ is the $r\times 1$ matrix (i.e column) containing zeros. If the $i^{th}$ constraint has a “$\le$”, then we set $c_{t+1}^{(i)}=0$. Otherwise it contains an “$=$”, and we set $c_{t+1}^{(i)}=\psi$.

Next, we deal with the globally uniform constraints in (${𝐏}_{\mathrm{𝟐}}$). Without loss of generality, we assume that each constraint either contains a “$\le$” or an “$=$”. This is a safe assumption because if there is a constraint containing a “$\ge$”, then we can negate that constraint (which flips the inequality in that constraint to “$\le$”) and obtain a program that is still of the form we consider. We replace each $(T^{(i)}\mathbf{0})$ by $(T^{(i)}\mathbf{0}{I}_r)$, where $I_r$ is the $r\times r$ identity matrix. We also introduce a new ${(n+1)}^{st}$ block $(Z_r\mathbf{0}{I}_r)$ where $Z_r$ is the $r\times t$ matrix consisting of zeros, and $\mathrm{𝟎}$ is the $r\times 1$ matrix consisting of zeros. The coefficients of the new variables in the first $n$ blocks are set to $\psi$. That is, for each $i\in [n]$ and $j\in [t+2,t+r+1]$, we set $c_j^{(i)}=\psi$. For the ${(n+1)}^{st}$ block, we set the coefficients as follows. For each $j\in [t+r+1]$, we set $c_j^{(n+1)}=0$. Then, we set $b^{(n+1)}=2r\mathrm{\Delta }\cdot {\max }_{i\in [n]}\{b^{(i)}\}$, where $\mathrm{\Delta }$ is the largest entry in the constraint matrix $𝒜$ we started with. This concludes the construction of the program of the form (1.1).

Consider any feasible solution $y$ to the initial program (${𝐏}_{\mathrm{𝟐}}$). We define a solution $x$ to the constructed program of the form (1.1) of the same cost by setting the newly added variables appropriately. Firstly, we set $x_j^{(i)}=y_j^{(i)}$ for each $i\in [n]$ and $j\in [t]$. Recall that we assume w.l.o.g. that each globally uniform constraint contains either a “$\le$” or an “$=$”. For each $k\in [r]$, we deal with the $k^{th}$ globally uniform constraint in (${𝐏}_{\mathrm{𝟐}}$) as follows. We set $x_{t+1+k}^{(n+1)}=b^{(0)}[k]-𝒜[k,\cdot ]y$. Observe that $x_{t+1+k}^{(n+1)}$ is precisely the slack in the $k^{th}$ globally uniform constraint corresponding to solution $y$, and it is a number that is at most $2\mathrm{\Delta }\cdot b^{(k)}$. For each $i\in [n]$, we deal with the $i^{th}$ locally uniform constraint in (${𝐏}_{\mathrm{𝟐}}$) as follows. We set $x_{t+1}^{(i)}=b^{(i)}-{(1\mathrm{\dots }1)}^{⊺}{y}^{(i)}$. Finally, we set $x_{t+1}^{(n+1)}=b^{(n+1)}-{\sum }_{k\in [r]}{x}_{t+1+k}^{(n+1)}$. Similar to before, observe that $x_{t+1}^{(i)}$ is precisely the slack in the $i^{th}$ locally uniform constraint corresponding to solution $y$, and it is a number that is at most $b^{(i)}$. All other variables whose values have not yet been set in $x$, we set them to $0$.

We claim that $x$ is a solution to the constructed program (1.1) of the same cost. Firstly, it is a feasible solution by construction. Secondly, the cost is the same because any newly added variable that is set to a non-zero value corresponds to a “$\le$” constraint, and its coefficient is $0$ in (1.1). For the other direction, suppose that we are given a solution $x$ to (1.1). If the cost of $x$ is at least $\psi /2$ (which is strictly larger than the value of any feasible solution in (${𝐏}_{\mathrm{𝟐}}$)), then we can assert that (${𝐏}_{\mathrm{𝟐}}$) does not admit any solution. Otherwise, we obtain a solution to (${𝐏}_{\mathrm{𝟐}}$) of the same cost by discarding the newly added variables.

Now that we have established how to compute the solution of one program given that of the other, we move on to assess the time to construct (1.1) and the change in parameters. Each locally uniform constraint can be dealt with separately in time $𝒪(nt)$ by adding a column corresponding to it and updating $c$ appropriately. All globally uniform constraints can be dealt with together in time $𝒪(n\cdot rnt)$ by adding $r$ columns to each block, adding a new block, and then updating $b$ appropriately. In total, we can construct (1.1) in time $𝒪(n^{2}rt)$. Thus, given (${𝐏}_{\mathrm{𝟐}}$), we can construct (1.1) in time $𝒪(n^{2}rt)$.

Next, we assess the change in parameters: the number of blocks is $n+1$, the number of columns in a block is $t+r+1$, the number of globally uniform constraints is $r$, the parameter $q$ (which is the sum of the right hand sides of the locally uniform constraints) becomes $q+2r\mathrm{\Delta }{\max }_{i\in [n]}\{b^{(i)}\}$, and the largest entry $\mathrm{\Delta }$ remains the same since the newly added entries to the matrix are either $0$ or $1$. Applying Theorem 1, we can solve (1.1) in time