New Characterizations and Efficient Local Search for General Integer Linear Programming

We develop a new local search algorithm Local-ILP, which is efficient for solving general ILP that is validated on a variety of different problems (i.e., the MIPLIB dataset).

Our algorithm is developed with two new operators tailored for ILP and a three-mode search framework. The two operators, namely the tight move and lift move, are designed based on a characterization of the solution space of ILP that involves the concept of boundary solutions. We show that searching for boundary solutions is complete for feasible and finite ILP instances, i.e., all optimal solutions belong to boundary solutions. We also design a local search framework that switches between three modes, namely Search, Improve, and Restore modes. Depending on the state of the best-found solution and the current solution, the framework selects the process of the appropriate mode to execute, where appropriate operators are leveraged to improve the quality of the current solution according to different situations.

For the Search and Restore modes, the tight move operator is adopted, which jointly considers variables and constraints’ information and tries to make some constraint tight. To distinguish important constraints and help guide the search to find feasible and high-quality solutions, we design a tailored weighting scheme and scoring function to select operations. For the Improve mode, an efficient operator lift move is used to improve the quality of the objective function while maintaining feasibility. To drive lift move operation, we propose a way to compute a variable’s new candidate values called local domain reduction. Additionally, we also design a specialized scoring function for lift move to pick a good operation. By putting these together, we develop a local search ILP solver called Local-ILP.

Experiments are conducted to evaluate Local-ILP on the MIPLIB benchmark, in which the ILP instances labeled hard and open are selected. We compare our algorithm with the state-of-the-art non-commercial ILP solver SCIP, as well as the state-of-the-art commercial solver Gurobi. Experimental results show that, in terms of finding a good feasible solution within a reasonable time, Local-ILP is competitive and complementary with Gurobi, and significantly better than SCIP. Experiments also demonstrate that Local-ILP significantly outperforms another local search algorithm Feasibility Jump, which won 1st place in MIP 2022 Workshop’s Computational Competition. Moreover, Local-ILP establishes 6 new records for MIPLIB open instances by finding the new best solutions.

Finally, we perform a theoretical analysis of our operators and algorithm based on the concept of boundary solutions, which are more explicit and closely related to the original form of ILP than the integral hull in polyhedral theory. We show that all feasible solutions visited by our algorithm are boundary solutions, efficiently avoiding unnecessary regions.

The remainder of the paper is organized as follows: Section 2 introduces the ILP problem and basic concepts for the local search algorithm, and then presents some basic characterizations of ILP we leveraged to design our local search algorithm. Section 3 proposes our local search framework for solving general ILP. In sections 4 and 5, we introduce two new operators that are key techniques for different modes of the framework. Section 6 presents detailed descriptions of how the Local-ILP algorithm implements the framework. The experimental results of our algorithm on public benchmark instances are reported in Section 7. The analysis of our algorithm is presented in Section 8, followed by the conclusion and future work in Section 9.

An instance of generalized ILP has the following form:

over integer variables $\bm{x}$, where $\bm{A}\in\mathbb{R}^{m\times n}$, $\bm{b}\in\mathbb{R}^{m}$, $\bm{c}\in\mathbb{R}^{n}$, $\bm{x^{l}}\in(\mathbb{Z}\cup-\infty)^{n}$ and $\bm{x^{u}}\in(\mathbb{Z}\cup+\infty)^{n}$ are given inputs. The goal of ILP is to minimize the value of the objective function, while satisfying all constraints. We denote the i-th constraint in the constraint system $\bm{A}\bm{x}\leq\bm{b}$ by $con_{i}$: $\bm{A_{i}}\bm{x}\leq\bm{b_{i}}$, where $\bm{A_{i}}\in\mathbb{R}^{n}$, $\bm{b_{i}}\in\mathbb{R}$. The coefficient of $x_{j}$ in $con_{i}$ is $A_{ij}$ and $con_{i}$ contains $x_{j}$ if $A_{ij}\neq 0$. The variables’ bounds are denoted by $\bm{x^{l}}\leq\bm{x}\leq\bm{x^{u}}$; they are indeed parts of $\bm{A}\bm{x}\leq\bm{b}$, but will be treated separately from the coefficient matrix $\bm{A}$ in practical algorithms. The infinite value of $\bm{x^{l}}$ or $\bm{x^{u}}$ indicates no lower or upper bound on the corresponding variable.

A solution $\bm{\alpha}$ for an ILP instance $Q$ is a vector of values assigned for each variable. Specifically, $\alpha_{j}$ denotes the value of $x_{j}$. A solution $\bm{\alpha}$ of $Q$ satisfies $con_{i}$ in $\bm{A}\bm{x}\leq\bm{b}$ if $\bm{A_{i}}\cdot\bm{\alpha}\leq\bm{b_{i}}$, and $\bm{\alpha}$ is feasible if and only if it satisfies all constraints in $Q$. The value of the objective function of a solution $\bm{\alpha}$ is denoted as $obj(\bm{\alpha})$.

The local search algorithm is well used in solving combinatorial optimization problems. It explores the search space comprised of all solutions, each of which is a candidate solution. Normally, local search starts with a solution and iteratively alters the solution by modifying the value of one variable, in order to find a feasible solution with a high-quality objective function value. The key to local search is to decide, under different circumstances, which variables to modify and to what value, to get a new solution.

An operator defines how the candidate solution is modified, i.e., given variables to be modified, how to fix them to new values. When an operator is instantiated by specifying the variable to operate on, we obtain an operation. Given an operation $op$, the scoring function $score(op)$ is used to measure how good $op$ is. An operation $op$ is said to be positive if $score(op)>0$, which indicates that performing $op$ can improve the quality of the current solution.

Currently, the most common theoretical tool to analyze ILP is polyhedral theory [22], of which the key component is the convex hull of all feasible solutions of an ILP (i.e., the integral hull), and it has established the equivalence between searching for the optimal solution and finding this integral hull. This equivalence provides additional theoretical explanations for the cutting plane method, which is widely used in commercial solvers, although the integral hull is hard to compute and could be done in exponential time.

However, the characterizations by integral hull are difficult to use for analyzing local search for ILP, as there is no explicit form of this integral hull and also no direct relations with constraints in the original form of ILP. To facilitate the analysis of local search algorithms for ILP, we introduce some new way to characterize solutions of ILP based on its original form, and show our local search operators and combinations are suited for ILP. We briefly present here the intuitions that we leveraged to design our local search algorithm and will give precise formulations and theoretical arguments in Section 8.

An important feature of ILP, is the linearity of its objective function and constraints. Although the ILP does not have true convexity, it still has some properties of this nature, such as all feasible solutions staying in the feasible domain of its LP relaxation, which is convex. From this aspect, some simple facts could be observed. Let $J=\{1,...,n\}$, $\bm{e}_{j}$ be the unit vector with $1$ in $j$-th coordinate and $0$ in all other places. There is:

In other words, solutions that lie between two solutions that are different in only one dimension have the objective value that is also in between (or could be equal to) the two solutions.

Thus we could make the following observations about ILP feasible solutions and optimal solutions, which will be used later to design new operators and our three-mode framework:

(1) all feasible solutions lie within a region (the integral hull)

(2) all optimal solutions lie in the “boundary” of the region

Furthermore, we call a set of solutions a complete space to an ILP if all optimal solutions are contained in this space. We will formalize the concept of “boundary” in Section 8 and show that boundary solutions are complete for an ILP. We design our operators based on this concept such that all feasible solutions our algorithm visits are boundary solutions.

For a search algorithm, we can distinguish three different stages: (1) enter the integral hull (= find a feasible solution)

(2) find good solutions within the integral hull (= improve the quality of a feasible solution)

(3) get back to the integral hull when jumping out of it (= from an infeasible solution to get a good feasible solution)

The characters of the three stages motivate us to propose a three-mode framework, enabling our algorithm to have different behaviors in different situations, improving efficiency.

In the following, we will first give descriptions of our algorithm, and later present the formal definitions and analysis of our algorithm in Section 8.

For solving ILP that may have infinite domains, a naive operator of local search is to modify the value of a variable $x_{j}$ by a fixed incremental value $inc$, i.e., $\alpha_{j}:=\alpha_{j}\pm inc$. However, the fixed $inc$ has drawbacks that it is hard to choose and may be instance-dependent: if $inc$ is too small, it may take many iterations before making any violated constraints satisfied; if $inc$ is too large, the algorithm may even become so troublesome that it can never satisfy some constraints, and therefore be hard to find a feasible solution. Our tight move operator modifies variable values according to the current solution and constraints and thus automatically adapts to different situations, which is defined below.

Our tight move operator has two properties while keeping the $x_{j}$’s bounds satisfied:

(1) If $con_{i}$ is violated, it makes $con_{i}$ as close to being satisfied as possible while minimally influencing other constraints, as it selects the minimal change of $x_{j}$ to make $con_{i}$ do so.

(2)If $con_{i}$ is satisfied, it can push $x_{j}$ to its extreme value, which keeps $con_{i}$ satisfied. This explores the influence of the maximal change of $x_{j}$ to the objective function and also helps to jump out of the local optimum.

The tight move operator, which considers both the information of variables and constraints, takes the idea to modify a variable’s value to make a constraint tight. This idea resembles some previous works which also consider some form of tightening constraints. The Simplex algorithm [23] adopts similar operator and is commonly used to solve linear systems. Specifically, the tight move is inspired by an operator named critical move, which is proposed in the context of solving Satisfiability Modulo Integer Arithmetic Theories [24]. The critical move operator assigns a variable $x_{j}$ to the threshold value making a violated literal (similar to a constraint in ILP) true, while our tight move operator keeps variables satisfying their global bounds, and allows modifications from satisfied constraints, thus expanding the range of operations to choose from.

We will show in Section 8 that although ILP does not have classical convexity, our tight move operator still has good theoretical properties to produce promising solutions and avoid visiting unnecessary search space.

During the local search process, scoring functions are used to compare different operations and pick one to execute in each step. For the tight move, we propose a weighted scoring function. Our scoring function has two ingredients: a weighting scheme to distinguish important constraints, and score computations to select operations.

To maintain the feasibility of a feasible solution, we must ensure that it satisfies every constraint. Therefore, we propose the local domain reduction to compute such a range to change a variable.

In order to compute $lfd(x_{j},\bm{\alpha})$, we consider the feasible domain of $x_{j}$ in $con_{i}$, which is denoted as $ldc(x_{j},con_{i},\bm{\alpha})$ meaning the $x_{j}$ can vary within this interval while keeping the satisfiability of $con_{i}$, assuming other variables in $\bm{\alpha}$ keep unchanged, where $con_{i}$ is a constraint containing $x_{j}$. Specifically, let $\Delta=\bm{b_{i}}-\bm{A_{i}}\cdot\bm{\alpha}$, $ldc(x_{j},con_{i},\bm{\alpha})$ is derived according to the sign of $A_{ij}$:

(1) If $A_{ij}\ <\ 0$, then $ldc(x_{j},con_{i},\bm{\alpha})=\left[\alpha_{j}+\left\lceil{\Delta}/{A_{ij}}\right\rceil,+\infty\right)$

(2) Otherwise, $ldc(x_{j},con_{i},\bm{\alpha})=\left(-\infty,\alpha_{j}+\left\lfloor\Delta/{A_{ij}}\right\rfloor\right]$.

Then, with $ldc(x_{j},con_{i},\bm{\alpha})$, we calculate the local feasible domain of $x_{j}$ as follows:

Clearly, moving $x_{j}$ within integers of $lfd(x_{j},\bm{\alpha})$ does not break the feasibility of $\bm{\alpha}$ as long as the other variables are kept constant. So once the current solution $\bm{\alpha}$ is feasible, we can choose a reasonable integer in $lfd(x_{j},\bm{\alpha})$ for updating $x_{j}$ to improve the quality of the objective function. We create the lift move operator for this purpose, which is based on the local domain reduction.

If there are multiple variables in the objective function, then multiple $lm$ operations could be constructed. To guide the search in Improve mode, we customize a scoring function to select a $lm$ operation. Since all $lm$ operations will maintain the feasibility of the solution, we propose lift score to measure the improvement of the objective function.

In the Improve mode, we pick the operation with the best $score_{lm}(op)$.

In Search mode (Algorithm 3), the goal of the algorithm is to find the first feasible solution.

If the algorithm fails to find any positive $tm$ operations in violated constraints, it first updates the weights of constraints according to the weighting scheme described in Section 5 (line 6). Then, it picks a random violated constraint to choose a $tm$ operation with the greatest $score_{tm}$ (line 7).

The Improve mode (Algorithm 4) seeks to improve the quality of the objective function’s value while maintaining the feasibility of a feasible solution. If no such operation is found, the algorithm will break feasibility, obtain an infeasible solution, and enter Restore mode.

If the algorithm fails to find any positive $lm$ operation, it randomly picks a variable $x_{j}$ in the objective function, and performs a simple unit incremental move in $x_{j}$ according to its coefficient $c_{j}$ (line 4-5), specifically, if $c_{j}<0$, then $\alpha_{j}=\alpha_{j}+1$; otherwise $\alpha_{j}=\alpha_{j}-1$. If a unit incremental move will break the global bound of a variable, we randomly select another. A unit incremental move that keeps all global bounds must exist; otherwise, all variables are at the corresponding bound value that improves the objective function, which means the current solution is the optimal solution and we could finish the search.

For an infeasible solution obtained from the Improve mode, the Restore mode (Algorithm 5) focuses on repairing the feasibility to obtain a new high-quality feasible solution.

If the algorithm fails to find any positive $tm$ operations in both violated and satisfied constraints, it updates weights and picks a random a $tm$ operation similarly to the Search mode (lines 6-7).

The difference between Restore and Search modes is that Restore mode will pick $tm$ operations in satisfied constraints because the number of violated constraints in Restore mode is usually much less than in Search mode. More importantly, since $score_{improve}(op)$ is meaningful in Restore mode while a constant in Search mode, $score_{tm}$ used in Restore mode considers the quality of the objective function while Search mode does not.

The classical techniques of local search are applied in our algorithm, including random sampling, forbidding strategy, and restart mechanism.

Random Sampling: We use a sampling method called Best from Multiple Selections, dubbed as BMS [27]. For all constraints that need to be considered in Search mode and Restore mode, the algorithm randomly samples a certain number of them, and for all operations derived from the sampled constraints that could be selected to apply, the algorithm randomly samples a certain number of these operations and selects the one with the highest $score_{tm}$. There are in total five parameters in the algorithm to control the number of samples, including: $c_{v}$ for violated constraints, and $o_{v}$ for $tm$ operations in sampled violated constraints; $c_{s}$ for satisfied constraints, and $o_{s}$ for $tm$ operations in sampled satisfied constraints; $o_{r}$ for $tm$ operations in a random violated constraint.

Forbidding Strategy: We employ a forbidding strategy, the tabu strategy [28, 24], to address the cycle phenomenon. After an operation is executed, the tabu strategy forbids the reverse operation in the following $tt$ iterations, where $tt$ is called tabu tenure, and we used the same setting as [24] to set $tt=3+rand(10)$. The tabu strategy is directly applied to Local-ILP. If a $tm$ operation that increases (decreases, resp.) the value of an integer variable $x_{j}$ is performed, then it is forbidden to decrease (increase, resp.) the value of $x_{j}$ by $tm$ operation in the following $tt$ iterations.

Restart Mechanism: The search is restarted when $\bm{\alpha^{*}}$ is not updated for enough steps iterations, which is simply set as 1500000. At each restart, we crossover $\bm{\alpha^{*}}$ and random integers within variables’ bounds to reset $\bm{\alpha}$: $x_{j}$ is assigned to $\alpha^{*}_{j}$ or a random integer within $[x^{l}_{j},x^{u}_{j}]$ with 50% probability of each, respectively. Additionally, all weights are restored to 1 when restarting.

In the MIPLIB dataset, each instance may contain multiple types of constraints⁶⁶6https://miplib.zib.de/statistics.html, such as knapsack constraints, set covering constraints, and so on. We categorize all instances by the type of their main constraint class. We mark the instance as hybrid if it contains multiple main constraint classes.

The results of the comparison with all competitors in terms of the ability to find a feasible solution, the quality of the best-found solution, and the primal integral are shown in Tables 2, 3 and 4, respectively.

To analyze the effectiveness of our proposed new techniques in Local-ILP, we tested 6 variations of our Local-ILP algorithm as follows:

(1) To analyze the effectiveness of the tight move operator. We modify Local-ILP by replacing the tight move operator with the operator that directly modifies an integer variable by a fixed increment inc, leading to two versions v fix 1 and v fix 5, where inc is set as 1 and 5, respectively.

(2) To analyze the effectiveness of the lift move operator and Improve mode, we modify Local-ILP by removing the Improve mode from the framework and using only Search and Restore modes, leading to the version v no improve. To analyze the effectiveness of the Restore mode, we modify Local-ILP by removing the Restore mode from the framework and returning to Search mode when Improve mode stuck in a local optimum, leading to the version v no restore.

(3) To compare different ways to escape from the local optimum in Improve mode. We modify Local-ILP by replacing the unit incremental move with operators that have larger step sizes, leading to two versions v per bound and v per random, where the step size is set as the distance to the bound of the variable and a random size between 1 and the bound, respectively.

We compare Local-ILP with these modified versions on the benchmark, with 10s, 60s, and 300s time limits. As shown in Table 5, Local-ILP outperforms all other variations, confirming the effectiveness of the strategies.

To examine the stability of Local-ILP which involves randomness, we run Local-ILP 10 times with 10 different seeds on the benchmark for 10s, 60s, and 300s time limits.

For all 10 times, we denote the average primal integral of each time by $avg_{P(T)}$, and the standard deviation of the primal integral by $std_{P(T)}$. The experimental results presented in Table 6 demonstrate that, for each time limit, the values of ${std_{P(T)}}/{avg_{P(T)}}$ for Local-ILP are less than 1.1%, indicating Local-ILP exhibits stable performance.

The optimal solutions for the instances labeled open have not yet been reported. The current best solutions known for each open instance are available on the official website of MIPLIB 2017.

As shown in Table 7, Local-ILP established the new best-known objective values for 6 instances. These 6 instances contain different constraint types, which simultaneously demonstrate the strong solving power of Local-ILP and its applicability to diverse types of problems. The new records have been submitted to MIPLIB 2017 and have been accepted; the links to the website are denoted in the footnotes ⁷⁷7https://miplib.zib.de/instance details sorrell7.html ⁸⁸8https://miplib.zib.de/instance details supportcase22.html ⁹⁹9https://miplib.zib.de/instance details cdc7-4-3-2.html ¹⁰¹⁰10https://miplib.zib.de/instance details ns1828997.html ¹¹¹¹11https://miplib.zib.de/instance details scpn2.html, except for one to be published in the next updates.

We assume an instance of ILP of the form:

For the convenience of analysis, here we do not distinguish variables’ bounds from other constraints, thus all variables’ bounds are included in $\bm{A}\bm{x}\leq\bm{b}$, this is slightly different in Formula (1), and we take Formula (5) as the description of an ILP instance within Section 8. We assume at least one of the coefficients in the objective function is non-zero: $\exists j\in\{1,...,n\}$, $c_{j}\neq 0$, otherwise this ILP instance has a constant objective function $0$ and turns to be a satisfiability problem. Since our algorithm does not prove feasibility, we exclude this case in analyzing our algorithm. We denote the index set $I=\{1,...,m\}$, $J=\{1,...,n\}$, $J_{\neq 0}=\{j|c_{j}\neq 0,j=1,...,n\}$, $J_{\neq 0}\neq\emptyset$. Moreover, we assume that there is no variable that is free (does not appear in any constraint, has both infinite upper/lower bounds, and does not appear in $\bm{c^{\top}}\bm{x}$).

We assume the ILP has a finite optimal solution $\bm{x}^{*}$ and a finite optimal objective value $opt$. We assume $\bm{x}^{*}$ exists (while it does not have to be unique), as our algorithm aims to find good solutions and does not try to prove an ILP is infeasible. We assume there is a finite optimal objective value; otherwise there is no meaningful optimal objective value.

Now we present the concept of boundary solutions that will be used to analyze our operators, and show that all optimal solutions are boundary solutions.

For an ILP instance in the form of Formula (5), let polyhedron $P=\{\bm{x}\in\mathbb{R}^{n}|\bm{A}\bm{x}\leq\bm{b}\}$, then the set of feasible solutions of the ILP could be described as $P\cap\mathbb{Z}^{n}$, and all feasible solutions belong to the integer hull $P_{I}=conv(P\cap\mathbb{Z}^{n})$ where $conv(S)\subseteq\mathbb{R}^{n}$ for $S\subseteq\mathbb{Z}^{n}$ denotes the convex hull of a set of points $S$. Let $U=\cup\{\bm{e}_{j},-\bm{e}_{j}\}$, $j\in J$. We call $\bm{x}+\bm{d},\bm{d}\in U$ the neighbors of $\bm{x}$. Given $P=\{\bm{x}\in\mathbb{R}^{n}|\bm{A}\bm{x}\leq\bm{b}\}$, we define the set of boundary points of $P$: