BiqBin: A Parallel Branch-and-bound Solver for Binary Quadratic Problems with Linear Constraints

With the advent of data driven economy, a close cooperation between data science and mathematical optimization became a crucial driver for developing (i) new data science methods that are capable to reveal new knowledge hidden in the data and (ii) new algorithms for implementing these data science methods.

Studying how to group data instances according to their inner similarity, i.e., data clustering analysis [21, Ch.10.3], is a traditional and one of the most studied problems in data science. When we know that the data should be decomposed into a fixed number of groups, say K, and we want to find these K groups, we have the K-clustering problem, which is NP-hard [12].

Another interesting problem from data science, similar to the K-clustering problem, is the problem where we want to find a set of k vertices in a simple undirected weighted graph G such that the sum of the weights on the edges connecting these k vertices is maximum. This problem is called the densest k-subgraph problem [4] and can be formulated as an optimization problem with a quadratic objective function, one linear constraint, and binary ( \(0/1\) ) decision variables:Here, W is the weighted adjacency matrix of the underlying graph and \(\mathbf {e}\) is the all one vector.

The densest k-subgraph problem can be seen as a generalization of the Max-Clique problem [41] and also as a special case of the quadratic knapsack problem [34]. Even though it is easy to understand and it appears to be simple at first sight, it is not. In fact, it is one of the NP-hard problems for which even the approximability is not well understood. There is a huge gap between the best approximation algorithm and known inapproximability results [3].

The problems mentioned above are all special instances of binary quadratic problems with linear constraints, which are formally defined in Section 1.3. In practice, we typically solve such problems only approximately using heuristic algorithms. For example, in [6, 40, 42] one can find a vast number of such heuristics for the case of data clustering. However, to evaluate the performance of the heuristic algorithms, we still need the ground truth, i.e., the optimum solutions of the original problems. Therefore, solving (to optimality) constrained binary quadratic problems, like (DkS), on as large as possible problem instances is highly needed.

Many other optimization problems with clear real-life applications can also be represented in a similar way as a non-convex optimization problem in binary variables subject to linear constraints, e.g., the quadratic assignment problem, the stable set problem, and the graph coloring problem (see, e.g., [37] for definitions). Again, appropriate (meta)heuristic algorithms are used to approximately solve instances of big size, while optimum solutions for test instances of small or medium size are still needed to evaluate these heuristics.

We will use the following notation: \(\mathcal {S}_n\) denotes the space of symmetric \(n\times n\) matrices, \({\mathbb {R}}^n\) is the n-dimensional space of all n-tuples of real numbers and \({\mathbb {R}}^{m\times n}\) is the space of all \(m\times n\) real matrices. By \(\mathbf {e}_n\) we denote the vector of length n with all entries equal to one. Usually, its dimension is clear from the context, so we write only \(\mathbf {e}\) . Similarly, \(\mathbf {0}\) denotes the zero vector or the zero matrix. Given \(\mathbf {x}\in {\mathbb {R}}^n\) , \(\operatorname{Diag}(\mathbf {x})\) is the \(n \times n\) diagonal matrix with \(\mathbf {x}\) on its diagonal; and \(\operatorname{diag}(X)\) is the vector with the diagonal elements of matrix X. By \(\operatorname{rk}(X)\) we denote the rank of matrix X.

For an optimization problem \(\bf {P}\) , we refer to its optimum value by \({\mathrm{OPT}}_{\bf {P}}\) . For the SDP relaxations of (Max-Cut) we use slightly different notations for their optimum values. For example, the optimum value of the SDP relaxation (MCSDP) is denoted by \({\mathrm{OPT}}_{\bf {SDP}}\) .

The central problem that we consider in the paper is the linearly constrained binary quadratic problem (BQP), which can be formulated as follows:with given data \(F \in \mathcal {S}_n\) , \(c\in {\mathbb {R}}^{n}\) , \(A\in {\mathbb {R}}^{m\times n}\) , and \(\mathbf {b}\in {\mathbb {R}}^{m}\) . We assume to have only linear equality constraints, in case of inequalities we may add a slack variable and decompose it into a weighted sum of boolean variables. This problem encompasses also the densest k-subgraph problem and is thus an NP-hard problem.

The mathematical optimization community is interested in BQP as a problem per se. In this context, the main challenge is developing new algorithms for it by (i) exploring and exploiting new properties of the problem, (ii) using new results from other (non-optimization) areas (such as algebraic geometry), (iii) making new combinations of existing algorithms with best practical or theoretical performance, and (iv) exploiting best available high-performance hardware and software.

Lasserre [27] has proved that any instance of (BQP) can be transformed into an instance of the Max-Cut problem, which is an NP-hard optimization problem [12, 22] on graphs. This problem is among the most studied combinatorial optimization problems, it has connections to various fields of discrete mathematics and models a wide range of applications. The transformation is based on an exact penalty approach, that was further explored and advanced in a recent paper of two co-authors of this paper [16].

The Max-Cut problem can be defined as follows. Suppose a weighted undirected graph \(G=(V,E)\) is given, where V is the set of vertices, E is the set of edges, and each edge \(e\in E\) has the weight \(w_e\in {\mathbb {R}}\) . The Max-Cut problem asks to find a partition of V into two parts \((S, V\backslash S)\) such that the sum of the weights of the edges having one endpoint in S and the other one in \(V\backslash S\) is maximized.

If \(W=(w_{ij})\) is the weighted adjacency matrix of G, i.e., \(w_{ij} = w_{ji} = w_e\) for \(e=\lbrace i,j\rbrace \in E\) , and the Laplacian matrix of the graph associated with W isthen computing the maximum cut amounts to solving the following binary quadratic problem in variables from \(\lbrace - 1,1\rbrace\) :However, it is straightforward to transform it into a \(0/1\) problem by using the following simple linear substitution \(\mathbf {x}= 2\mathbf {z}- \mathbf {e}\) and by observing thatsince \(L\mathbf {e}=0\) . Hence, we obtain the equivalent formulation of the Max-Cut, which has a structure of (BQP) (with no linear constraints):

The reason why the community is interested in reformulating (BQP) into (Max-Cut) is related to the fact that for the latter problem, there exists a wide range of approximate and exact methods and solvers and we want to employ them in solving the former problem. We are particularly interested in the methods computing global optima of both problems.

The main contribution of this paper is the BiqBin solver for BQP, which outperforms the existing solvers on some special cases of BQP. The core of the BiqBin solver is the exact penalty algorithm EXPEDIS [16], meaning that we first transform every instance of (BQP) to a corresponding instance of (Max-Cut) and then solve the resulting instance of () BiqBin is coded in C and non-trivially improves existing solvers by

Additionally, we demonstrate practical efficiency of BiqBin by providing an extensive benchmarking with BiqCrunch [25], GUROBI [14], and SCIP [11] on the list of four special cases of BQP, including the Max-Cut problem, the unconstrained binary quadratic problem, the densest k-subgraph problem and randomly generated binary quadratic problems with linear constraints.

We can observe that BiqBin is performing very well on the instances with a small number of linear constraints, while on the instances with many linear constraints, the solvers which directly exploit the structure of the constraints are very competitive. Additionally, as we show on the benchmark instances, BiqBin is also highly scalable.

The BiqBin solver is a Max-Cut based solver for (BQP) instances that solves their reformulation to (Max-Cut) instances using a B&B algorithm.

The main ingredients of BiqBin are:

The overall performance of the algorithm is determined by the quality of the lower and upper bounds, computed in each B&B node, by the computational efficiency of computing these bounds, and by the strategy of exploring the B&B tree. The main reason why we cannot solve large instances on personal computers is that the B&B tree grows too big, i.e., the pruning is too slow and a single processor is not capable to explore all the generated nodes in the tree.

Exact penalty reformulation. The procedure which reformulates every instance of (BQP) into an instance of (Max-Cut) is described in Section 4.

Bounding procedure. The starting point of the algorithm is the strengthened SDP relaxation ((MCHYP)). Since the number of inequalities is too large to solve this SDP by standard solvers, we will use a bundle method to find an approximate solution to the partial Lagrangian dual. By dualizing only the inequality constraints, we obtain the nonsmooth convex partial dual functionwhere \(\gamma\) are the nonnegative dual variables associated with the constraints \(\mathcal {A}_{\bf {HYP}}(X)\le \mathbf {e}\) . Evaluating the dual function \(f(\gamma)\) and computing the subgradient amounts to solving an SDP of the form (MCSDP), which can be efficiently computed using an interior-point method tailored for this problem. It provides us with the matching pair \((X_{\gamma }, \gamma)\) such that \(f(\gamma) = \mathcal {L}(X_{\gamma },\gamma)\) . Moreover, the subgradient of f at \(\gamma\) is given by \(\partial f(\gamma) = \mathbf {e}- \mathcal {A}_{\bf {HYP}}(X_{\gamma }).\) For obtaining an approximate minimizer of the dual problemwe use the bundle method [23, 35]. Let the current iterate be \(\hat{\gamma }\) . Suppose we have evaluated f at \(k \ge 1\) points \(\gamma _1,\ldots ,\gamma _k\) with matching pairs \(X_1,\ldots ,X_k\) and subgradients \(\mathbf {e}- \mathcal {A}_{\bf {HYP}}(X_i)\) for \(i\in \lbrace 1,\ldots ,k\rbrace\) . The bundle method combines the following two ideas:

In summary, the bundle method finds a new trial point by minimizing, for some prescribed parameter \(t \gt 0\) , the functionover the nonnegative orthant. As the method terminates, we obtain the approximate minimizer of the dual function, as well as the convex weights \(\lambda \ge 0,\sum _i \lambda _i=1\) , determining the matrix \(\mathcal {X}\lambda\) .

In order to obtain a tight upper bound for an instance of () we use the cutting-plane approach, where multiple k-gonal inequalities are added and purged in the course of running the bundle algorithm. First, the optimum solution of the basic semidefinite relaxation (MCSDP) is computed using an interior-point method followed by separating a set of triangle inequalities. After doing a few bundle iterations with this set of constraints, we purge all inactive constraints. Next, invoking the separation routine described in Section 3, new violated k-gonal inequalities are added. The problem with the new set of constraints is solved and the process is iterated as long as there is a significant decrease of the upper bound.

Since the primal solution matrix X is not available, we use the same idea as in BiqMac to purge some inequality constraints. We look at the values of the corresponding dual multipliers in the vector \(\gamma\) . If the value of some dual multiplier is close to zero, this indicates that the corresponding constraint is not active and we remove it.

After each iteration of computing an upper bound and separating triangle inequalities, higher order k-gonal inequalities are added in order to further decrease the bound. We monitor the maximum violation of triangle inequalities \(r_{\text{tri}}=\max (\mathcal {A}_3(\mathcal {X} \lambda) - \mathbf {e})\) and as soon as the number is sufficiently small, the heuristic from Section 3 is used to add some strongly violated pentagonal inequalities to the relaxation. Similarly, as the maximum violation of pentagonal inequalities \(r_{\text{pent}}\) drops below some threshold, new heptagonal inequalities are separated and added to the relaxation. In our numerical tests, we have used the thresholds \(r_{\text{tri}} \lt 0.2\) and \(r_{\text{pent}}\lt 0.4\) .

In order to improve the performance of BiqBin, we can stop the bounding routine when we detect that we will not be able to prune the current node in the B&B tree. We again borrow an idea from BiqMac. After some cutting plane iterations, we make a linear (and hence optimistic) forecast to decide whether it is worth doing more iterations. If the gap cannot be closed, we terminate the bounding routine, branch the current node, and start evaluating new subproblems. This is especially important in the parallel solver, since its efficiency depends on how quickly idle workers receive subproblems.

Furthermore, at the beginning of the algorithm, the number of bundle iterations should be small. We start with three iterations. Then this number is increased after each separation of new cutting planes, until a limit is reached. In the BiqBin solver, this limit value is set to 15. This is motivated by the fact, that in the beginning it will take a while until we have identified the right cutting planes and we do not waste time by trying hard to decrease the bound if the current set of hypermetric inequalities does not allow much progress.

Branching strategies. BiqBin uses two branching strategies which are based on the bundle of matrices obtained from the bounding routine. Once the bundle method terminates, the last column \(\mathbf {x}\) of the matrix \(\mathcal {X}\lambda\) is extracted. Due to the diagonal and positive semidefiniteness constraints on the feasible matrices of ((MCHYP)), all the entries of \(\mathbf {x}\) lie in the interval \([-1,1]\) . Similarly to BiqCrunch, the decision on which variable \(x_i\) the branching should be performed is based on the following two strategies:

In the \(0/1\) formulation, these rules are usually referred to as most-fractional and least-fractional rules, respectively. The difficult first rule is set as the default strategy for the BiqBin solver.

When branching, two new subproblems, in which the chosen branching variable is fixed accordingly, are created and the corresponding nodes are added in the B&B tree, i.e., to the priority queue of unexplored problems. Priority is based on the upper bound obtained from the bundle method. When selecting the next subproblem, a node with the worst upper bound is evaluated first.

Rounding heuristic. For generating high quality feasible solutions of the Max-Cut problem, we apply the Goemans-Williamson rounding hyperplane technique [13]. Let X be an optimal solution of some SDP relaxation of Max-Cut. By computing the Cholesky factorization \(X = V^\top V\) with column vectors \(v_i\) of V and selecting some random vector r, the cut \((S,V\backslash S)\) is obtained by setting \(S = \lbrace i \mid v_i^\top r \ge 0 \rbrace .\) Since in our case we are working with the partial Lagrangian, the information about the primal X is not available. Instead, we use a convex combination of bundle matrices \(\mathcal {X}\lambda\) as the input. The cut vector x obtained from this heuristic is then further improved by flipping the vertices and using a convex combination of \(\mathcal {X}\lambda\) and the cut matrix \(xx^\top\) . To summarize, for generating good cuts, we use the following iterative scheme:

Interestingly, in most of our numerical experiments, this heuristic finds the optimum solution already in the root node.

Strategy for faster enumeration of the B&B tree. As described in Section 3, simulated annealing is used to heuristically separate k-gonal inequalities. Adding these k-gonal inequalities to the model in each B&B node is not necessary, especially if one can not prune the node. In that case all the work done and time are wasted, since for bigger graphs we usually need to reach a certain depth in the B&B tree in order to prune the nodes, even if the bounding routine produces tight bounds. It is beneficial to check (before including cutting planes when processing the node) whether there is hope to prune the node or whether it is better to branch and produce smaller subproblems. We propose the following strategy.

In the root node we compute: (i) the bound \({\mathrm{OPT}}_{\bf {SDP}}\) of the basic SDP relaxation (MCSDP), which is not strong but is quick to compute; and (ii) the bound \({\mathrm{OPT}}_{\bf {HYP}}\) by iteratively including violated triangle, pentagonal, and heptagonal inequalities. Letand let \({\mathrm{LB}}\) denote the current lower bound. Then, at all other nodes, we first compute only the basic SDP bound \({\mathrm{OPT}}_{\bf {SDP}}\) . If the conditionis satisfied, this means, we are already close to the lower bound, and so we add cutting planes to compute the tighter bound \({\mathrm{OPT}}_{\bf {HYP}}\) in order to increase the probability of pruning the node. With this idea, we can efficiently traverse the B&B tree and only invest time into the bounding routine when it is needed. Numerical results show that overall this strategy produces more B&B nodes than necessary, but the performance of the algorithm improves; in particular, it has a positive impact on the parallel version described in Section 5.2.

Numerical illustration of the bounding procedure. When using a B&B algorithm, one has two options regarding the quality of the upper bound. On one hand, a strong upper bound can be computed by iteratively adding and purging multiple cutting planes. In this way, more work is done within a node but overall this approach produces fewer B&B nodes. On the other hand, one can efficiently compute slightly weaker upper bounds, hence the whole tree grows larger but if the time spent for evaluating each node is small it can be traversed faster. In BiqBin we use the first approach.

We take the Beasley bqp250.8 problem from the BiqMac library and plot the convergence curve for the bounding routines in the solvers BiqBin and BiqCrunch in the root nodes of the B&B trees. Figure 1 depicts the decrease of the dual function values in the course of a bound computation. We note that by adding higher order k-gonal inequalities, our bounding routine attains a tighter bound compared to BiqCrunch. Consequently, BiqBin creates a smaller B&B tree, which consists of 81 nodes, while BiqCrunch terminates after traversing 325 nodes.

In this section, we describe how the algorithmic ingredients from sequential B&B algorithm are combined into a parallel solver which utilizes distributive memory parallelism. This is done in a similar fashion as in the parallel B&B solver for the stable set problem [18], introduced by some of the authors of this paper.

The load coordinator–worker paradigm with distributive work pools is applied, in which the rank 0 process becomes the master process carefully managing the status of each worker (idle or busy), while different workers concurrently explore branches of the B&B tree. Each worker has its own local queue of subproblems and the work is shared when one of them becomes idle. The master node knows the status of each worker and acts as a load coordinator receiving messages and, based on their content, replying in an appropriate manner.

At the beginning of the algorithm, the load coordinator reads and broadcasts the original graph to the workers and initializes the solution. It is important that every process has the knowledge of the original graph, since construction of subproblems via branching and encoding the MPI messages is done based on this information. All the data about the B&B nodes is encoded as an MPI structure, which is used in communication between different workers in order to efficiently exchange and construct the subproblems.

Next, the master process evaluates the root node and distributes the best lower bound. After the bounding step, two new subproblems are generated, which are sent to the first two idle processes. Afterwards, its job is restricted to monitoring the status of the workers, counting the number of B&B nodes, and distributing the best solution found so far. In order to solve multiple SDP relaxations in parallel, the workers need to send and receive subproblems and keep load balance. In this way, the work is shared and the whole dynamic B&B tree is enumerated faster.

After the initialization phase, the master process waits for three types of messages sent by the workers. Firstly, if the worker’s local queue is empty, the message is sent informing the master that the process is idle and can receive further work.

Secondly, the master process receives messages regarding the load balance and sharing of work. Throughout the algorithm, the master node manages a bool array containing the statuses of the workers. They specify whether the process is active or idle. The load coordinator uses this information to reply to the requests received from the workers during the branching step, in which the process wants to share one of the newly generated subproblems with some other free worker. The load coordinator sends a message specifying which workers are free and the value of best lower bound.

And thirdly, during the execution of the algorithm, the working processes compute multiple candidates for optimum solutions. The master node is keeping track of the currently best value and the corresponding solution. When a new solution is received, the value is compared, updated if necessary, and distributed back during the communication phase.

After a worker computes lower and upper bounds, it compares these values to see whether this branch of enumeration tree can be safely pruned or further branching is needed and construction of new subproblems takes place. In the latter case, a request message is sent to the load coordinator, asking for idle processes to share one of the newly generated subproblems or subproblems left in the queue from previous branching processes. If no idle worker is available, the generated subproblems are placed in the worker’s queue and the work continues locally. Otherwise, subproblems are encoded and sent to available idle workers. This is also where the exchange of the best lower bound happens.

When all the workers become idle, the master process sends a message to finish and the algorithm terminates. The algorithms for the load coordinator and the workers are summarized in Algorithms 1 and 2.

Lastly, we explain how the parallel version benefits from the strategy using the variable “ \({\mathrm{diff}}\) ” (see Equation (6) on page ). If the number of available workers is large, we need to reach a certain depth in the B&B tree in order for the processes to receive the work. Until this happens in the algorithm, the workers are idle. To fully exploit all HPC resources available, we need a strategy for the worker processes to start evaluating the nodes as soon as possible. This is where the idea from Section 5 helps. After the load coordinator evaluates the root node, the best lower bound found and the variable \({\mathrm{diff}}\) are distributed to all the workers. When the first two idle processes evaluate the generated subproblems, typically the value of the basic SDP relaxation is such that condition (7) is not satisfied. Hence, on the first few levels of the B&B tree, the workers compute only the basic SDP bound. Then branching of new subproblems takes place and idle workers quickly receive the generated subproblems. After the bounds are such that condition (7) is valid, hypermetric inequalities are added to compute tighter bounds \({\mathrm{OPT}}_{\bf {HYP}}\) approximately. This implies that more nodes are pruned, meaning that the size of the B&B tree decreases, and thus the algorithm terminates faster.

Max-Cut instances. The first family of (BQP) benchmark instances consists of Max-Cut instances. We selected benchmark instances from the BiqMac library that have been considered as hard by other solvers (see [39] for more details). Furthermore, we used rudy, a graph generator written by Giovanni Rinaldi [36], to construct new (hard) Max-Cut instances with 180 nodes, similar to the ones in the BiqMac library.

Unconstrained BQP instances. The second family of (BQP) instances consists of two sets of unconstrained BQP. The instances of the first set, due to Beasley [2], have density 0.1 with values uniformly chosen from \([- 100,100]\) . The second set of unconstrained BQP instances were generated by Billionet and Elloumi [5] according to [32]. They have size \(n \in \lbrace 100,120,150,200\rbrace\) and different densities. The diagonal coefficients are in the range \([ - 100,100]\) , while the off-diagonal ones are in from \([ - 50,50 ]\) . Similar instances with size 250 and density 0.1 were taken from the BiqMac library. For more details on these instances we refer to [39].

Densest k-subgraph instances. The third family of benchmark problems consists of instances of the (DkS), described in Section 1.1. It is a binary quadratic problem with one linear constraint, also called the cardinality constraint. These problems, also known as cardinality Boolean quadratic programming problems, include a wide range of applications in telecommunications or chemistry, see [29] for further details. We solve the benchmark instances that can be found at [26]. They have different sizes \(n \in \lbrace 120,140,160\rbrace\) , densities \(d \in \lbrace 0.25, 0.50, 0.75\rbrace\) , and values for the parameter \(k \in \lbrace n/4, n/2, 3n/4\rbrace\) . For testing the parallel version of BiqBin, we also created similar instances with sizes \(n \in \lbrace 180,200\rbrace\) . They can be found at the BiqBin web page.⁹

Randomly generated instances of (BQP). The fourth family of instances consists of randomly generated instances with a varying number of linear constraints. These instances have size \(n = 100\) and up to 15 constraints, their description can be found in [16] and they are available at [15].

In this section, we compare the sequential version of BiqBin with BiqCrunch [25], GUROBI [14], and SCIP [11] solvers on the four families of problems introduced in the previous section. We note that for the first three families of benchmark problems GUROBI and SCIP did not solve any instance within three hours, hence we report results for these problems only for BiqBin and BiqCrunch.

We present the results of the comparisons of these two solvers on the Max-Cut instances in Tables 1 and 2, the results for other unconstrained BQPs are in Tables 3 and 4, and the results for the densest k-subgraph problem are in Tables 5 and 6.

Before we list the tables, we summarize the results by the performance profile diagrams on Figure 3. These diagrams were created using the same raw data as we summarised in Tables 3 and 7. There were large time spans for each solver and each data set, so the diagrams are quite coarse-grained. Nevertheless, we can see that on the Max-Cut, unconstrained BQP and randomly generated instances of BQP our solver BiqBin is outperforming the other solvers. On the instances of the densest k-subgraph problem BiqCrunch is performing slightly better then BiqBin.

To highlight the difference between the solvers BiqBin and BiqCrunch, we also provide a relative performance profile in Figure 4. It shows, for each factor \(\alpha \in \lbrace 2^{-3},2^{-2},\ldots ,2^3\rbrace\) , the percentage of instances from the Max-Cut, unconstrained BQP and densest k-subgraph families of instances, respectively, for which the computation time required by BiqCrunch was shorter than \(\alpha\) times the computation time required by BiqBin. For the Max-Cut instances, this plot reveals that BiqCrunch does not require less than 1/2 of the time needed by BiqBin for any instance (i.e., it is never more than two times faster than BiqBin), but it is faster than BiqBin on about 13% of the Max-Cut instances, i.e., it is slower than BiqBin for about 87% of the Max-Cut instances. For the densest k-subgraph instances the situation is reversed: for approx. 84% of the instances, BiqCrunch is faster than BiqBin, while for the unconstrained BQP instances, both solvers perform similarly.

The first three columns in each of these tables contain the data about the instances: for each group of instances, we used several instances, e.g., for the Max-Cut problem and for the unconstrained BQP, we used 10 instances for every combination of size n and density d, and for the densest k-subgraph we used five instances. The second set of columns in each table contains the minimum, the maximum, and the average number of nodes in the B&B trees generated by BiqBin and BiqCrunch, respectively. Similarly, the third set of columns reports the minimum, the maximum, and the average computing times (in seconds) needed for the instances, and the last set of columns contains the minimum, the maximum, and the average initial gaps (in %).

For the Max-Cut problem, we can see in Tables 1 and 2 that the (sequential) BiqBin is always approximately two times faster than BiqCrunch and produces much fewer nodes in the B&B tree. One reason for this lies in the fact that the initial gap in the root node (the difference between the upper and lower bound, divided by the lower bound) is smaller for BiqBin.

For the unconstrained BQP (Tables 3 and 4), we can see that both, BiqBin and BiqCrunch, perform very similar. However, for BiqBin the B&B trees are much smaller, so there is potential for improving the computing times by increasing the computational efficiency in each B&B node.

For the densest k-subgraph problem (Tables 5 and 6), we can observe that BiqCrunch is in most cases slightly better (sometimes also much better) than BiqBin. The reason for this is that BiqCrunch has a specific version adapted to solving this problem, where the basic semidefinite relaxation is re-enforced with triangle and product constraints.

Turning to the randomly generated instances, GUROBI and SCIP performed well, while the BiqCrunch solver showed very poor performance. Thus, we report in Table 7 only the results for BiqBin, GUROBI and SCIP. The randomly generated instances have different intervals from which the entries of A, b, and F are selected, thus for every combination, we have 15 instances with an increasing number of constraints, from 1 to 15. Table 7 contains in the first column the name of the family of instances, the next three columns contain data about the numbers of instances (always 15), the sizes of problem n (always 100), and the densities of the data (the number of non-zero elements in F divided by \(n^2\) ). For each of the three solvers, we report the number of instances solved by the solver in the time limit of three hours. These instances are also considered in the performance profile 3d.

From Table 7 and from the performance profile 3d it can be observed that, overall, BiqBin is outperforming all the other solvers. It solves for each family at least two out of the 15 instances, and on average it solves 60% of the instances, while GUROBI and SCIP solve 48% and 16% of the instances, respectively.

We can see that BiqBin has worst performance on instances starting with 100_A_0_3_ in the middle part of the table. We do not have an unambiguous explanation for this, but the reason is very likely that the linear constraints enable generating good cutting planes, thus helping GUROBI and SCIP. Obviously, the linear constraints for the instances starting with 100_A_0_3_ yield more efficient cutting planes which enable these two solvers to more efficiently prune the B&B tree and therefore more often finish in a given time limit. Some of the 100_A_0_3 instances are also infeasible, while the instances having in the middle _b_0_ are always feasible (they have zeros on the right hand, i.e., \(\mathbf {b}= \mathbf {0}\) , hence \(\mathbf {x}= \mathbf {0}\) is always a feasible solution). GUROBI and SCIP are solving problems in the (BQP) form, while BiqBin first reformulates the instances into the(Max-Cut) format. This suggests to improve the performance of BiqBin by enhancing an early stopping condition for the case of infeasibility using Theorem 1.

In this section, we demonstrate how BiqBin is scaling across the high-performance computer that we used within the MPI framework. We measured wall-times needed to solve each instance using the sequential solver and the parallel solver with \(3, 6,\ldots ,48\) CPU cores (i.e., MPI processes). Note that one of the CPU cores was always reserved for the coordinator’s tasks, while the others were used as workers.

In the “sequential” columns of Table 8, we report the computational times needed by the sequential algorithm and the times needed for the computations in the root node of the B&B tree. The former are used to compute the speed-up factor while the latter are considered as the times for the non-parallelizable parts of BiqBin and are used to compute the upper bounds for the speed-up factors.

Pairs of columns denoted by “3 cores”, “6 cores”, etc. contain the sizes of the B&B trees and the wall-times of the parallel BiqBin using \(3, 6,\ldots\) CPU cores, respectively. For each instance inst from Table 8, we compute a vector of speed-up factors \({{\tt {speed-up} }_{{\tt inst}}}\) , defined as (see, e.g., [Equation (3.4)][31])where \({\tt time}_1\) is the time needed by the sequential solver while \({\tt time}_i\) denotes time needed by parallel solver using i CPU cores. These factors are aggregated in Table 9 for each family of instances and are depicted in Figure 5.

Times from the “root time” column of Table 8 are used to estimate the proportion of the BiqBin solver that is non-parallelizable. We denote this estimate by s and compute it as the sum of the 6th column divided by the sum of the 5th column (result is \(s=0.0136\) ). The upper-bounds for the speed-up factors, are computed according to Amdahl’s law by formula ([31, Equation (3.6)])where n corresponds to the number of CPU cores available. These upper-bounds are depicted by the green curve in Figure 5. We can see that speed-up factors are close to the theoretical bound for the problems g05 and deviate a lot for the problems pw09 and w05.

Table 8 also contains the numbers of nodes in the B&B trees generated by the parallel solver with different numbers of CPU cores. We can see that these numbers are different from the sequential solver and are varying with the number of CPU cores. The reason is mainly in the fact that parallel computations in different CPU cores are no more deterministic. Generating cutting planes in the process of solving ((MCHYP)) includes random numbers. We fix the seed for the random number generator at the beginning of the computation in each CPU core. However, when we vary the number of CPU cores, the amount and the order of the computational work in the cores change. Consequently, different cutting planes might be applied to the same B&B node and this results in slightly different bounds \({\mathrm{OPT}}_{\bf {HYP}}\) and finally in different sizes of B&B trees.

In this section, we report numerical results obtained by parallel BiqBin on the Max-Cut, unconstrained BQP and the densest k-subgraph instances that were to the best of our knowledge not solved so far. We used 200 CPU cores and the results are collected in Tables 10 and 12. We can see that the hardest families of instances are again g05 and w09. To solve them, we need to compute more than one million B&B nodes, which took at least six hours of wall-time.