Parallelized Conflict Graph Cut Generation

At the preprocessing stage, given a $\mbox{MIP}=(\mathcal{M},\mathcal{N},\mathcal{I},A,b,c,\circ,\ell,u)$, we perform a one-round simple presolve to detect set packing constraints and conflicting knapsack constraints for subsequent conflict graph cuts generation, which is described in Alg 1. This applies well-studied techniques from the CG literature (see e.g. Achterberg et al. [2]).

All techniques mentioned in line 3 can be found in [2, Sec 3.1 and 3.2]. For $A_{i\cdot}x\geq b_{i}$, we transfer it as $-A_{i\cdot}x\leq-b_{i}$. And for $A_{i\cdot}x=b_{i}$, we treat it as two inequalities $A_{i\cdot}x\leq b_{i}$ and $A_{i\cdot}x\geq b_{i}$.

In line 4 of Alg 1, if $A_{i\cdot}x\leq b_{i}$ is a set packing constraint, we call it an original set packing (OSP) constraint and collect all such constraints in the set $\mathcal{S}_{osp}$ (line 6).

In line 10, for $x_{j}$ with $j$ such that $A_{ij}\neq 0$, the lower or upper bounds are tightened as follows: if $A_{ij}>0$, $u_{j}:=\min\{u_{j},b_{i}/A_{ij}\}$; otherwise, $\ell_{j}:=\min\{\ell_{j},b_{i}/A_{ij}\}$.

In line 11 of Alg 1, if $\mbox{PBC}(i)$ is a set packing constraint, we call it a inferred set packing constraint (ISP) and collect all such constraints in the set $\mathcal{S}_{isp}$.

In line 14 we collect all (i.e. both original and PBC-inferred) conflicting knapsack constraints (CK) in the constraint set $\mathcal{S}_{ck}$ (line 15).

The complexity of Alg 1 is $O(NNZ)$, where $NNZ$ is the number of nonzero elements in $A$. We implement this procedure solely in serial since it executes very quickly in practice.

Following PBC generation, we proceed to detect maximal cliques from conflicting knapsack constraints with the Clique Detection method of Brito and Santos [10, Algorithm 1]. We modify the method, described as (serial) Alg 2: namely, instead of returning all detected maximal cliques in a single outputted set, we separately return the first detected maximal clique (see line 7 in Alg 2) in $\mathcal{S}_{org}$ (the original maximal clique) and all other cliques in $\mathcal{S}_{other}$ (other maximal cliques). This is used for more intensive clique extension and merging applied specifically to the original cliques, described in Section 2.4. Alg 2 for $\mathcal{S}_{ck}$ is, in turn, called in parallel via Alg 3.

Furthermore, in line 2 of Alg 3, we randomly shuffle and then partition $\mathcal{S}_{ck}$ into subsets of equal cardinality (modulo the last processor) from Alg 1 as a simple and fast heuristic for load balancing. Moreover, this can be efficiently parallelized [6, 33]. The same trick is used in Alg 5 and Alg 7. Note that the general problem of optimal load balancing—dividing up tasks with known computational load as evenly as possible across $k$ cores—is an NP-hard partitioning problem [11, 12]. The shuffling heuristic is justified both by average-case analysis as well as computational experiments indicating high parallel efficiency on hard instances.

The complexity for Algorithm 2 is $O(n\log n)$ (on average) and $O(n^{2})$ (worst case) (see [10, Page 4, Paragraph 2-3]), where $n$ is the number of elements in the constraint. The average-case complexity (average over random shuffles) for Algorithm 3 is $O(mn\log n/k)$, where $m$ is number of constraints, and worst case is in $O(mn^{2}/k)$ (see Proposition A.1 in Appendix A).

Following clique detection from PBC-derived conflicting knapsack constraints, we proceed to CG construction. Suppose there are $n_{\mathcal{B}}$ binary variables in $\mathcal{B}$, and so $n_{\mathcal{B}}$ complementary variables. We represent the CG with a sparse matrix $G\in\{0,1\}^{2n_{\mathcal{B}}\times 2n_{\mathcal{B}}}$, with the first set of rows $j\in\{1,...,n_{\mathcal{B}}\}$ representing the original variables $x_{j}$ and the next set $j\in\{n_{\mathcal{B}}+1,...,2n_{\mathcal{B}}\}$ representing the complements $\bar{x}_{j-n_{\mathcal{B}}}$. We build $G$ by Alg 5 in parallel.

We choose a sparse adjacency matrix here instead of a clique table in order to enable faster clique extension (see Section 2.4) by avoiding redundant computation. A clique table, however, is substantially more efficient with memory, so as a workaround to the sparse data structure we set limits on the number of nonzeros considered (see Section 2.6). Moreover, on cliques that we choose not to extend (namely $C_{other}$), we adopt the more memory-efficient data structure described in Section 2.3 of [10].

The initialization of the clique set $\mathcal{C}$ in line 1 includes: the trivial pairwise conflicts between variables and their complements (i.e. $G_{j,j+n_{\mathcal{B}}}=1$ for all $j=1,...,n_{\mathcal{B}}$); the clique of $\mathcal{S}$ from set packing constraints obtained by Alg 1; and cliques extracted in Alg 3.

As with Alg 3, in line 2 random shuffle and partition is applied as a load balancing heuristic.

In line 5, each $G(i)$ is updated in parallel using the clique subset $\mathcal{C}^{i}$ via Alg 4.

In lines 7-12, we reduce the individual $G(i)$ to the global $G$ via binary combination (see Fig 1). In line 9, the OR logical merge is applied to every element in CGs, e.g., $G(i)_{j_{1},j_{2}}:=\max\{G(i-2^{d-1})_{j_{1},j_{2}},G(i)_{j_{1},j_{2}}\}$.

For a clique set $\mathcal{C}$, suppose for each clique $Q\in\mathcal{C}$, whether a variable $x_{j}\in Q$ is a Bernoulli distribution with the probability $p$. Then the average runtime complexity of Alg 5 is $O(mn_{\mathcal{B}}^{2}p^{2}/k+\log k\cdot n_{\mathcal{B}}^{2})$ (see Proposition A.4 in Appendix A). Furthermore, the worst case runtime complexity of Alg 5 is $O(mn_{\mathcal{B}}^{2}/k+\log k\cdot n_{\mathcal{B}}^{2})$ (see Corollary A.5 in Appendix A).

After generating the conflict graph, we apply clique strengthening [1, 2, 10]—in particular, we modify the Clique Extension of Brito and Santos [10]: a greedy algorithm used to generate one strengthened clique based on the original clique and the CG (see [10, Algorithm 2]). Our modification (see Alg 6) involves (potentially) multiple clique strengthenings, and the increase in computational load is made practical by parallelization.

In line 1 of Alg 6, we generate a list $L$ including all variables $u$ that do not belong to the clique, but for which $u$ conflicts with all variables in the clique. In lines 19-20, we isolate a largest extended clique $\mbox{Clq}^{\prime}$ from the clique set $\mathcal{C}$. We distinguish this longest extended clique from all other cliques in order to perform cut management, described in Section 2.5.

Strengthening each clique is an independent procedure, and so we propose to apply strengthening in parallel as Alg 7.

Complexity for Alg 6 is $O(n_{\mathcal{B}}^{2})$ (see Lemma A.6 in Appendix A) for one clique; thus $O(mn_{\mathcal{B}}^{2})$ for $m$ cliques. Alg 7 has complexity $O(mn_{\mathcal{B}}^{2}/k)$ (see Proposition A.7 in Appendix A).

After obtaining extended cliques from Alg 7, we check for domination between extended cliques and remove dominated cliques. For a clique defined as a set pack (Equation (2)), we will only store the index set $\mathcal{S}$; thus, to check the domination between two cliques $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, we only need to check whether $\mathcal{S}_{1}\subseteq\mathcal{S}_{2}$ or $\mathcal{S}_{2}\subseteq\mathcal{S}_{1}$, which can be performed in $O(n_{\mathcal{B}})$ due to sparse data structures.

Given $m$ cliques, the domination checking process is in $O(m^{2}n_{\mathcal{B}})$. The parallel version of this, Alg 8, has complexity $O(m^{2}n/k)$.

Cut management applies the results from previous subsections (which are conducted in presolve) throughout the branch-and-cut tree.

In Alg 1, three constraint sets are generated: $\mathcal{S}_{osp},\mathcal{S}_{isp},\mathcal{S}_{ck}$. Subsequently, cliques are extracted from $\mathcal{S}_{ck}$ using Alg 3, yielding $\mathcal{C}_{org}$ and $\mathcal{C}_{other}$. Applying Alg 5 on the four sets $\mathcal{S}_{osp},\mathcal{S}_{isp},\mathcal{C}_{org},\mathcal{C}_{other}$, the CG is constructed. Cliques are then strengthened/extended from the three sets $\mathcal{S}_{osp},\mathcal{S}_{isp},\mathcal{C}_{org}$ with parallel Alg 7. As a result, this process yields $7$ clique sets: $\mathcal{C}^{long}_{osp},\mathcal{C}^{other}_{osp}$ (from $\mathcal{S}_{osp}$), $\mathcal{C}^{long}_{isp},\mathcal{C}^{other}_{isp}$ (from $\mathcal{S}_{isp}$), $\mathcal{C}^{long}_{org},\mathcal{C}^{other}_{org}$ (from $\mathcal{C}_{org}$), and $\mathcal{C}_{other}$. However, the number of cliques from these sets can be quite large, and so judicious management is critical for practical application. For instance, incorporating all corresponding inequalities at the root node is impractical both due to resource limitations on linear programming solves as well as potential numerical issues from excessive cuts.

Our cut management procedure applies only a few key inequalities from the longest cliques are added at the root node and apply the remaining cuts throughout the search tree:

First, we replace $\mathcal{S}_{org}$ in the MIP formulation (deleting it in line 5 of Alg 1) with the strengthened constraints from $\mathcal{C}^{long}_{osp}$.

Second, all cliques from $\mathcal{C}_{osp}^{other}$, $\mathcal{C}_{isp}^{other}$, $\mathcal{C}_{org}^{other}$, and $\mathcal{C}_{other}$ are incorporated via user cuts, which are placed in the cut pool and may be added to the model at any node in the branch-and-cut search tree to cut off relaxation solutions (see e.g. [16, Page 769, Lazy Attribute]).

Third, for cliques from $\mathcal{C}_{isp}^{long}$ and $\mathcal{C}_{org}^{long}$, we either add them to the root node or else place them as user cuts depended on following rules:

1. 1.

   If the size of $\mathcal{S}_{org}$ is excessively large (greater than 20000) or small (less than 100), we add cliques from $\mathcal{C}_{isp}^{long}$.

2. 2.

   If $2000<|\mathcal{S}_{org}|\leq 20000$ and $|\mathcal{C}_{org}^{long}|/|\mathcal{S}_{osp}|<0.02$, add cliques from $\mathcal{C}_{org}^{long}$ to the root node; otherwise, place them as user cuts.

3. 3.

   If $100<|\mathcal{S}_{org}|\leq 2000$ and $|\mathcal{C}_{org}^{long}|/|\mathcal{S}_{osp}|<1$, add cliques from $\mathcal{C}_{org}^{long}$ to the root node; otherwise, place them as user cuts.

4. 4.

   If the size of $\mathcal{S}_{org}$ is excessively large (greater than 12000) or small (less than 1500), add cliques from $\mathcal{C}_{isp}^{long}$ to the root node; otherwise, place them as user cuts.

This heuristic attempts to prioritize some inequalities and place them directly in the formulation at the root node, leaving the remaining inequalities for Gurobi to manage in the search tree via user cuts.

CG Management is time-consuming (see e.g. [2, Page 491 Paragraph 2]) and memory-consuming. For instance, in CG construction, the sparse adjacency matrix costs $O(n^{2})$ memory for the clique with length $n$. However, we are not obliged to consider every possible variable, nor identify every clique, etc.; as such, to avoid excessive memory or time costs, we set the following limits.

For CG construction in Alg 4, if the length of the clique $Q$ is greater than $3000$, we randomly select $3000$ elements from $Q$ to construct CG. For Clique Extension in Alg 7, we set each thread at most process $10^{6}$ nonzero terms and return generated extended cliques to the main thread. For Clique Merging in Alg 8, if the number of cliques is more than $10^{5}$, we give up to conduct Alg 8.

In the following experiments we analyze the parallel efficiency of our parallel CG management algorithms (excluding branch-and-cut solver time). We select the $163/173$ cases for which the CG procedure can be run with a single thread within a time of $600$ seconds. We furthermore exclude in the subsection $55$ cases which can be trivially handled in less than $0.1$ seconds with a single thread. On the remaining $108/173$ cases, parallel speed-ups (serial runtime divided by the parallel runtime) are presented in Fig. 2. The parallel efficiency, naturally, depends on the total work available. When CG takes more than 50 seconds we observe substantial speedups up to 64 threads, and on more modest instances there is a tail-off at 32 threads.

In this section we consider the impact of our CG procedure on total solver time, using Gurobi as our benchmark. As aforementioned, we focus on 173 cases from MIPLIB that have constraints applicable to our CG procedure. We exclude from comparison in this subsection an additional 74 cases for which our method cannot derive additional inequalities beyond standard CG techniques, namely those instances for which $\mathcal{S}_{isp}\cup\mathcal{S}_{ck}=\emptyset$. For these 74 excluded cases our procedure will simply reduce to a standard (albeit parallelized) CG procedure that is already incorporated into most solvers’ presolve routines (including Gurobi). This leaves us with 99 cases; among these 99, 3 more are excluded in the aggregated results as they exceed the excessive clique size caused issues in either JuMP (1 instance) or our procedure (2 instances). Note that a simple maximum clique-size limit would allow for skipping the CG procedure with minimal effect on runtime in these 3 instances. Thus $96/173$ cases are considered for comparison purposes to analyze the impact of our CG routine on Gurobi.

We present results in two parts: one with the $96$ cases excluding CG management/overhead time, and one with $96$ cases including CG management time and varying the number of threads used. Since our CG implementation will duplicate certain similar efforts in Gurobi’s CG itself, we expect a fully-integrated implementation to have less overhead than what is reported; hence, the two times presented provide optimistic (without overhead), and conservative (with overhead) bounds on potential impact, respectively.