A Vectorized Positive Semidefinite Penalty Method for Unconstrained Binary Quadratic Programming

Due to the potential of the unconstrained quadratic programming problem, many authors have studied the optimality conditions and sought high-quality solution methods for the UBQP problem in recent years. Beck and Teboulle doi:10.1137/S1052623498336930 derived the sufficient and necessary global optimality conditions. Based on the work of Beck and Teboulle, Xia Xia2009NewOC obtained tighter sufficient optimality conditions and explored the relationship between the optimal solution of the problem and that of its continuous relaxation. These studies have provided us with a deeper understanding of UBQP problems, leading to the development of efficient solution strategies.

One of the main categories of methods for solving UBQP is exact solvers. Typically, this is achieved by a tree-search approach that incorporates a generalized branch-and-bound technique. In the 1980s, both Gulati et al. GULATI1984121 and Barahona et al. 1989Experiments described a branch-and-bound algorithm for UBQP. More recently, Krislock et al. Krislock2017BiqCrunchAS present BiqBin, an exact solver for binary quadratic problems. The approach is based on an exact penalty method that efficiently transforms the original problem into an instance of Max-Cut, and then solves the Max-Cut problem using a branch-and-bound algorithm. The numerical results show that BiqBin is a highly competitive solver.

Heuristic algorithms are also often used to solve UBQP problems exactly. The most common are tabu search inproceedings , quantum annealing PMID:35140264 . For example, Wang et al. 10.1007/978-3-642-29828-8_26 introduce a backbone-guided tabu search algorithm that involves a fundamental tabu search process alternating with a variable fixing/freeing phase, which is determined by the identification of strongly determined variables. Due to the NP-hardness, the optimal solution usually cannot be found in polynomial time. Therefore, it is necessary to find a fast method that is not strictly exact. Based on this, many scholars have proposed a number of inexact methods.

Among the inexact methods, there are mainly methods based on relaxation and penalty function methods. There are some penalty function methods based on vector variables. Liu et al. Liu2018ACA and Nayak et al. RePEc:spr:jcomop:v:39:y:2020:i:3:d:10.1007_s10878-019-00517-8 have adopted a continuous approach, converting binary constraints into continuous ones by using Fischer-Burmeister nonlinear complementarity problem functions, then it can be further carried on the smoothing processing by aggregate function. However, this method is relatively sensitive to the initial values, and the parameter updates are complex and may not serve as exact penalty. Furthermore, the exact penalty method is widely used in many methods, transforming the binary constraints into a set of inequality constraints. Yuan et al. yuan2017exact introduced an alternative approach by replacing the binary constraints with a Mathematical Programming with Equilibrium Constraints (MPEC). This transformation converts the UBQP problem into a biconvex problem. While this algorithm demonstrates high computational efficiency, the solution accuracy still needs to be improved for certain data sets. Luo et al. luo2019enhancing introduced a special penalty method called conditionally quasi-convex relaxation based on semidefinite relaxation and they showed that the constructed penalty function is exact. However, the subproblem of this algorithm still needs to solve an SDP subproblem, which is time-consuming.

Matrix-based relaxation methods are relatively common approaches, with typical examples including the reconstruction-linearization technique Sherali2007RLTAU and positive semidefinite relaxation (SDR) fujie1997semidefinite . SDR is a classical method for solving QCQP, which transforms quadratic terms of vector variables into matrix variables, drops the rank-one constraint, relaxes the problem to a semidefinite programming (SDP) problem, and utilizes SDP solvers for the solution. In many cases, SDR can obtain the exact solution to QCQP or prove to achieve a certain approximation ratio 2023Exact ; 5447068 . It also has wide applications in practice.

When the matrix solution of SDR is not rank one, the traditional strategies are to use projection or random sampling 5447068 . Recently, other methods have also been studied to further obtain high-quality rank one solutions from SDR solutions 10.1093/imanum/draa031 ; article . Gu et al. provided a broader perspective to SDR research by introducing a positive semidefinite-based penalty (PSDP) formulation 10.1093/imanum/draa031 . When the penalty factor is zero, it corresponds exactly to the semidefinite relaxation of QCQP problems. Moreover, if the penalty factor is sufficiently large, the penalty function problem is equivalent to QCQP. Therefore, by gradually increasing the penalty factor, the matrix solution of SDR can reach rank one. In addition, the authors used a proximal point algorithm (PPA) to solve the penalty function subproblems, which demonstrated remarkable performance when applied to UBQP, with a significant proportion yielding exact solutions for UBQP.

While the solution quality provided by PSDP is high, its significant computational cost as an SDR-based approach cannot be overlooked compared to vector methods. Therefore, a crucial research question is how to maintain the solution quality of PSDP while significantly reducing its computational cost, ideally to a form comparable to the computational efficiency of vector methods.

In the PSDP framework, we select a specified initial value for the penalty factor and find that the penalty function problem can depend only on the diagonal elements of the matrix. Following the penalty factor update rule of PSDP, the generated sequence of PSDP penalty problems still ensures that it is only related to the matrix diagonal elements, thus causing the matrix problem in PSDP to degrade to a vector problem, achieving the vectorization of PSDP. Furthermore, in the vectorized PSDP algorithm we propose, the algorithm for solving subproblem and penalty update rule maintain the same form as PSDP, while incorporating some improvements in the details. Therefore, the vectorized PSDP not only maintains the quality of solutions in PSDP, but also keeps the computational cost at the level of vector methods.

To further accelerate the convergence and improve the solution accuracy, we apply an algorithm accelerated by the projection alternating Barzilai-Borwein (BB) method 2005Projected . This method uses the projection step to effectively handle the 0-1 constraints, while the alternating BB method provides a robust and efficient framework for updating the search direction and step size. Additionally, we introduce a criterion to update the penalty parameter during the iterations. This ensures that the active set is appropriately updated throughout the iteration process, thereby ensuring the convergence of the algorithm. We tested our algorithm on one randomly generated dataset and five benchmark datasets, comparing its performance with four existing methods. Experimental results show that our algorithm performs well in terms of both time efficiency and accuracy.

In section 2, we will present our vectorized exact penalty function and its transformation process. In section 3, we will present our subproblem, two methods for solving the subproblem, and an approach for updating the penalty factor. Section 4 will provide the theoretical foundation of our algorithm, including the exact penalty function theory and the convergence of the subproblem algorithm. Finally, in section 5, we will demonstrate the comparison of our algorithm with several other methods across different datasets through numerical experiments.

For a general QCQP as in (4), PSDP 10.1093/imanum/draa031 adds a matrix variable $Z\succeq 0$ and a penalty term $P\cdot Z$ with $P\succeq 0$ to the objective, proposing the following penalty problem (5).

When $P$ is 0, (5) is equivalent to a standard SDR.

In the UBQP problem, with $Q_{0}=Q$, $Q_{i}=e_{i}e_{i}^{T}$ ($e_{i}$ is a vector with 1 in the $i$-th position and 0 elsewhere), $g_{0}=0.5b_{i}$, $g_{i}=-0.5e_{i}$ and $c_{i}=0$ as in (4), the penalty problem (5) corresponds to the following problem.

where $z_{i}$ is the $i$-th diagonal element of $Z$.

PSDP employed a PPA method to solve problem (6), with each PPA iteration involving the solution of an SDP problem, necessitating the repeated solving of SDPs. PSDP begins with $P=0$, which is equivalent to SDR, and then updates the penalty factor $P$ by $P=P+\alpha Z$.

Solving SDP is the major computational obstacle in PSDP, so we consider what kind of problem structure can make problem (6) easy to solve.

We note that when $Q+P$ is a diagonal matrix, the solution of $Z$ depends only on its diagonal elements. This is because the first constraint involves only the diagonal elements of $Z$, and when $Q+P$ is a diagonal matrix, the objective function also involves only the diagonal elements of $Z$. For the remaining positive semidefinite constraints on $Z$, as long as its diagonal elements are non-negative and the non-diagonal elements are set to zero, it naturally satisfies the constraints. Therefore, we only need to make $Q+P$ a diagonal matrix so that there exists a diagonal matrix solution for $Z$. In addition, since $Z$ is a diagonal matrix, $Q+P$ actually only updates the diagonal elements when updating $P=P+\alpha Z$. Therefore, the new penalty problem still has a diagonal matrix solution for $Z$.

The construction of $P$ is straightforward that we simply need to set the off-diagonal elements of $P$ to be the negation of the off-diagonal elements of $Q$. Supposing that we choose an appropriate $P$ such that $Q+P$ also becomes a diagonal matrix $P^{\star}$, problem (6) is thus transformed into the following form.

where $p_{i}$ is the $i$-th diagonal element of $P^{\star}$, the inequality $x_{i}^{2}-x_{i}\leq 0$ has been replaced with $0\leq x_{i}\leq 1$, and $Z$ has been eliminated.

Billionnet and Elloumi Billionnet2007UsingAM have discussed the form of problem (7), exploring how to construct an appropriate $p$ to preserve the convexity of problem (7) while also preserving equivalence with problem (2) as much as possible. In contrast, we continue to follow the PSDP approach of updating $p$, allowing the quadratic problem (7) to be non-convex, and solving it using iterative algorithms.

In PSDP, when applying the PPA to solve (6), the $k$-th iteration is carried out in the following form.

where $G_{1}=\{(x,Z)|x_{i}^{2}-x_{i}+z_{i}=0,i=1,2,\dots,n,Z\succeq 0\}$ and $H\left(x^{k}+d,Z,P\right)=x^{T}Qx+b^{T}x+(Q+P)\cdot Z$. Its proximal term is $d^{T}Pd$. For the vectorized problem (7), we observe the implementation of this PPA as follows:

where $F_{1}=\{x|0\leq x_{i}\leq 1\}$, $h\left(x,p\right)=x^{T}(Q-P^{\star})x+(b+p)^{T}x$, and $p$ is the vector composed of diagonal elements of $P^{\star}$.

By utilizing the nature of $Q+P=P^{\star}$, it can be observed that the objective function of (9) is a linear function of $d$. It is known that minimizing a linear function over box constraints leads to all components of the optimal solution being taken on the boundaries. This causes $x$ to reach the critical point of (7) in one iteration. However, it is typically a low-quality local solution, since problem (7) is non-convex. Therefore, we have to modify the proximal term to better solve (7) using PPA.

Building on (9), we introduce an additional quadratic term $d^{T}H_{p}d$ to the existing proximal term, where $H_{p}\succ 0$ to ensure strict convexity of the new problem and avoid linearity. Furthermore, we need to ensure that the new proximal term $(H_{p}+P)$ is positive semidefinite since we do not strictly require $P\succeq 0$ here. Lastly, it is necessary for $H_{p}$ to be a diagonal matrix to facilitate the computation of the subproblems. Therefore, we define $H_{p}$ in the following form.

where $\epsilon>0$ is an arbitrarily small number. Thus, $H_{p}\succ 0$ and

where the final step utilizes the fact that a diagonally dominant matrix is positive semidefinite.

Setting the proximal term to be $(H_{p}+P=H_{p}-Q+P^{\star})$, we obtain the $k$-th iteration as follows:

In fact,

Thus, (10) is separable in variables, equivalent to solving n one-dimensional problems. Without considering box constraints, the optimal solution is as follows:

where $\tilde{d}^{k}_{j}$ is $j$-th element of $\tilde{d}^{k}$, $\nabla h\left(x^{k},p\right)^{j}$ is the $j$-th element of $\nabla h\left(x^{k},p\right)$ and $h_{p}^{j}$ is the $j$-th diagonal element of $H_{p}$. Then $d^{k}$ can be obtain by projecting $\tilde{d}^{k}$ onto box set $(F_{1}-x^{k})$, that is,

Therefore, the next iterative point $x^{k+1}=x^{k}+d^{k}$ is given by

where $\mathcal{P}$ is a projection operator defined by

We present the algorithmic framework of this PPA in Algorithm 1, which is used to solve the penalty problem (7) with a given penalty factor $p$.

While the previously proposed algorithm is straightforward and computationally efficient, it exhibits relatively long computation time for local convergence. Consequently, we opt to use the projection alternating BB method introduced by Dai and Fletcher 2005Projected to speed up the computational process.

The projection alternating BB method strategically alternates between utilizing long and short BB step sizes, as outlined below:

Here,

where $s_{k-1}=x_{k}-x_{k-1}$ and $y_{k-1}=g^{k}-g^{k-1}$. In the context of problem (7), we have $g^{k}:=\nabla h\left(x^{k},p\right)$. Then the update of $x^{k+1}$ is expressed as

We present the framework of this projection BB method in Algorithm 2.

For the optimization problem defined by (7), the termination criterion is established based on the Karush-Kuhn-Tucker (KKT) conditions.

The Lagrangian function is given by

where $u$ and $v$ correspond to the dual variables. Then the first-order optimality conditions can be expressed as follows:

By introducing the index sets $I=\{i|0<x_{i}<1\}$, $K=\{i|x_{i}=1\}$, $J=\{i|x_{i}=0\}$, and defining $\bar{Q}=Q-P^{\star}$ and $\bar{b}=b+p$, the conditions can be rewritten as:

Consequently, we present the termination criterion as

where

Due to the adoption of the penalty factor update from PSDP, specifically $P^{\star}=P^{\star}+\alpha Z$, where both $P^{\star}$ and $Z$ are diagonal matrices in the vectorized PSDP algorithm, only the diagonal elements of $P^{\star}$ need to be updated, that is,

Here, $z^{s}=(x_{1}^{s}-x_{1}^{s^{2}},x_{2}^{s}-x_{2}^{s^{2}},\dots,x_{n}-x_{n}^{s^{2}})^{T}$, $x^{s}$ denotes the best value of $x$ after the $s-1$ th update of $p$, and $p^{s}$ denotes the penalty parameter after s updates.

The choice of the parameter $\alpha$ significantly impacts the overall efficiency of the algorithm. When $\alpha$ is too small, it may lead to smaller iteration steps and an increase in the number of iterations, thus increasing the runtime and slowing down convergence. Conversely, if $\alpha$ is too large, it can easily cause $x$ to rapidly approach the boundaries during the update process, potentially halting at a low-quality local solution. Therefore, selecting an appropriate value for $\alpha$ is crucial.

We opt for an appropriate $\alpha$ from the perspective of the active set. Ideally, after updating $p$, we desire the new optimal solution $x$ of (7) to have a different active set compared to before the update, while ensuring that the change is not too drastic. In the following lemma, we provide an $\alpha$ that guarantees a change in the active set.

Within this section, we introduce two initialization techniques: reformation and random perturbation.

Reformation. Without loss of generality, we initialize the vectorized PSDP by setting $p=0$. For the UBQP problem (2), leveraging the binary nature of $x\in\{0,1\}^{n}$, we can reformulate the objective in an equivalent form by introducing an arbitrary vector $\gamma$, yielding

Although the new $Q$ (i.e. $Q+Diag(\gamma)$) and $b$ (i.e. $b-\gamma$) do not change the solution of problem (2), the solution of problem (7) subject to the constraint $x\in[0,1]^{n}$ is affected by $\gamma$. Therefore, $\gamma$ will affect the initial solution of $x$ in our algorithm.

As noted in the previous subsection, it is imperative to prevent the iterative point sequence from approaching the boundary too rapidly, as this could impede the discovery of high-quality solutions. Consequently, the selection of an appropriate $\gamma$ is crucial to ensure that the initial solution for $x$ resides within the interior of the feasible region. Notably, this goal can be achieved by setting $\gamma$ sufficiently large. This is due to the inactive nature of the boundary defined by $0\leq x_{i}\leq 1$, where the first-order condition is satisfied as

Analysis of (21) reveals that as $\gamma$ approaches infinity, $x_{i}$ converges to $1/2$. To maintain the constraint $0<x_{i}<1$, an appropriately large $\gamma$ can be chosen. For example, setting $\gamma_{i}>\sum_{j=1}^{n}2|Q_{ij}|+|b_{i}|$ for $i=1,\ldots,n$ in the selection of $\gamma$ ensures $0<x_{i}<1$, as can be deduced from (21).

Random Perturbation. For integer problems (where the elements of $Q$ and $b$ are integers), certain components of the algorithm may occasionally stall around 1/2. Therefore, following the approach of Tang and Toh tang2023feasible , we introduce a symmetric random perturbation from a normal distribution.

where $\mu_{ij}\sim\mathbb{N}(0,\sigma)$ and $\mu_{ij}=\mu_{ji}$.

Due to the continuity of problem (7), small perturbations do not alter the optimal solution. By incorporating random perturbations, the iterations no longer stagnate around 1/2, thereby enhancing the robustness and stability of the algorithm.

Having outlined all the details of the vectorized PSDP algorithm, we now present its framework in Algorithm 3.

In this subsection, we prove that the penalty problem (7) is exact, which means that its optimal solution is also an optimal solution of problem (3), when all the components of $p$ are larger than a certain constant.

To simplify the description of the problem, we rewrite the problem (3) as

where $f(x)=x^{T}Qx+b^{T}x$ and $F=\{x|x\in\{0,1\}^{n}\}$. Also we rewrite the problem (7) as

where $h(x,p):=f(x)+\sum_{i=1}^{n}p_{i}(x_{i}-x_{i}^{2})$ and $G=\{x|0\leq x\leq 1\}$.

For the convenience of our proof, let us first propose two lemmas.

By leveraging Lemma 2 and 3, we can determine the exact penalty of problem (24), as demonstrated in Theorem 4.1.

In Algorithm 1, we apply a proximal point algorithm to solve the subproblem (7), then we present the convergence of the Algorithm 1 in Lemma 4 and Theorem 4.2.

In Algorithm 2, an alternative BB method is adopted, which is non-monotonic, allowing for instances where $h\left(x^{k+1},p\right)$ may increase during certain iterations. Dai and Fletcher 2005Projected provided counterexamples illustrating that the Projected BB and Projected Alternating BB methods, without line searches, could cyclically oscillate between multiple points, and such behavior is atypical. Furthermore, Li and Huang li2023note theoretically proved that the alternative BB method without a monotone line search converges at least R-linearly with a rate of $1-1/\kappa$ under certain properties for convex unconstrained quadratic problems. However, we have not yet found a relevant theoretical result for such a non-convex problem. In the subsequent section, numerical experiments will be conducted to demonstrate that the PABB method does not exhibit cycling phenomena.

Due to the absence of any line search, it is challenging to provide the theoretical convergence analysis for Algorithm 2. Here, we validate its practical convergence through numerical experiments. We present one instance from each of the six datasets, and Fig.1 illustrates the number of iterations used by Algorithm 2 in solving the subproblems. It can be observed from the figure that the number of iterations for the subproblems is mostly within 1000. For the instance with the QPLIB ID 3565, there is a small probability that the number of iterations for the subproblems exceeds 2000. However, even in such cases, all iteration counts remain below 6000, not reaching the maximum of 10000 iterations set for the experiments. This indicates to some extent that the PABB method without line search does not exhibit cycling, confirming the convergence of Algorithm 2 in practice.

Next, we validate the overall convergence of the two vectorized PSDP algorithms (PPA and PABB) in numerical experiments. Fig.2 illustrates the convergence process of the two algorithms in reducing the number of components violating the 0-1 constraints over six test cases. While not strictly monotonically decreasing, the overall trend is toward reduction and convergence to zero. This indicates that our algorithms ultimately converge to feasible points satisfying the 0-1 constraints.

The UBQP problem to be addressed is formulated as follows:

where Q is a symmetric matrix of order $n$. For this problem, we conducted experiments on four datasets, namely: random instances, QPLIB, BE, and GKA.

(Random Instances) We randomly generate the symmetric matrix $Q$ with dimensions of 50, 100, 150, 200, and 250. For each dimension, we create five test cases. Due to the stochastic nature of some algorithms, we run each algorithm 10 times for each test case and report the average GAP and CPU time. In cases where the runtime of the ALM algorithm exceeds 100 seconds or yields a non-optimal trivial solution of 0, we classify such instances as failures.

The detailed results are summarized in Table 1. We observe that the five algorithms do not exhibit significant differences in terms of runtime, with the SDR algorithm running relatively slower. In terms of solution quality, PPA achieves the smallest GAP with the lower bound, closely followed by PABB. However, ALM fails on average 1 to 4 times in each test case.

Additionally, we consider some lower-dimensional matrices with dimensions of 20, 40, 60, and 80, and apply the branch-and-bound method to obtain the optimal solution as the lower bound (lb). As shown in Table 2, PABB is relatively faster in terms of time, with a smaller GAP compared to the optimal value, and even the maximum GAP is only 1.02%. PPA is more accurate, achieving an average GAP of just 0.02% when solving examples with 80 dimensions. In contrast, ALM has a probability of failure during the solving process, and even when successful, the GAP is relatively high. The MPEC solution is also accurate, but with longer run times. On the other hand, the SDR algorithm has slower solving speeds and relatively lower solution quality.

This observation suggests that the two vectorized PSDP methods have a clear advantage in solution quality over other approaches. PABB generally outperforms PPA in speed, effectively speeding up the process.

(QPLIB) The QPLIB is a library of quadratic programming instances, which contains 319 discrete and 134 continuous instances of different characteristics. We conduct testing on all binary instances, totaling 23 cases, with dimensions ranging from 120 to 1225 in the QPLIB dataset.

As shown in Table 3, MPEC emerges as the fastest approach, with an average of 0.49 seconds per instance. PABB and ALM follow closely with an average of about 1 second per instance. PPA and SDR take longer. In terms of solution quality, PPA and PABB maintain their superiority.