Globally solving quadratic programs with convex objective and complementarity constraints via completely positive programming

Quadratic programs with complementarity constraints (QPCC) are NP-hard due to the nonconvexity of complementarity relation between the pairs of nonnegative variables. Most of the existing solvers are capable of solving QPCC by finding stationary solutions, which are not able to be verified as global or local optimal ones. In this paper, we aim to globally solve QPCC by a branch-and-bound algorithm, in which the doubly nonnegative (DNN) relaxation in each node is efficiently solved via an augmented Lagrangian method. The method is practically efficient due to the fact that the augmented Lagrangian function can be decomposed into two easy-to-solve subproblems. Computational results demonstrate the effectiveness of the proposed algorithm, with a particular highlight in only a few nodes for some instances.

the relaxations via an augmented Lagrangian method, which is efficient in small and medium scale. Secondly, although the SDP relaxations are polynomial-time tractable, they are not required to be solved exactly in every node, except for the leaf node due to the finiteness of the algorithm, for the purpose of overall efficiency. Instead, a good lower bound is obtained in an efficient way if the relaxation in one node can not be fully solved by our method. Computational results show that the proposal algorithm can solve medium-sized instances within reasonable amount of time.
The rest of this paper is arranged as follows. Section 2 presents the completely positive reformulation for QPCC and derives the doubly nonnegative program (DNP) as an SDP relaxation. Section 3 proposes the augmented Lagrangian method for efficiently solving the DNP relaxations. Section 4 states the branchand-bound algorithm for globally solving QPCC and the implementation details for improving the efficiency. Results of numerical experiment are presented in Section 5, and the conclusions are given in Section 6.
The following notations are adopted in this paper. S n denotes the cone of n × n real symmetric matrices, S n + denotes the cone of n × n real positive semidefinite matrices, and N n denotes the cone of n × n real symmetric matrices with nonnegative components. C n denotes the cone of completely positive matrices, that is denotes the diagonal vector of matrix A. For two matrices A, B ∈ S n , the inner product between A and B is defined as A • B = n i=1 n j=1 A ij B ij where A ij and B ij are the (i, j)-th components of matrices A and B, respectively. For two vectors a, b ∈ n , the Hadamard product between a and b, denoted by a • b, is the component-wise product of a and b.

2.
Copositive reformulation and doubly nonnegative relaxation. In [4], Burer has shown that problem (QPCC) is equivalent to the following completely positive problem in terms of optimal value: Since the completely positive cone C n+1 is intractable (ref. [7]), problem (CPP) is NP-hard in general. Due to the fact C n+1 ⊆ S n+1 ∩ N n+1 , problem (CPP) naturally has the following SDP relaxation: which is a doubly nonnegative program (DNP) because the eigenvalues and components of matrix Y are both nonnegative. By the theory of interior point method (IPM), DNP can be solved in polynomial time to any fixed precision. However, due to the existence of doubly nonnegative constraint, the Newton system in every iteration of IPM is of size O(n 4 ). Hence, IPMs are limited to fairly small or medium problem sizes in practice because the CPU time will increase dramatically as the problem size n goes larger. In recent years, alternating direction methods of multipliers (ADMMs) are developed to solve large-scale SDPs (ref. [15], [22] and [23]). In general, ADMMs decompose the augmented Lagrangian function into several subproblems, each of which is relatively easy to handle. In particular, several ADMMs are specially designed for DNP problem. For example, Lu et al. [17] proposed a primal ADMM for the DNP relaxations of a class of mixed integer quadratic programming problem, whose domain is separable for each x i . This special structure of feasible domain is utilized in the design of algorithm so that the solutions of decomposed subproblems have closed form. Burer [5] proposed an augmented Lagrangian method for DNP relaxation of completely positive programs. We highlight here that Burer used coordinate descent method for solving the decomposed subproblem from the augmented Lagrangian function. It has been shown that this procedure is actually equivalent to the primal ADMM in [17]. Now, we derive a property of problem (DNP), which allowes us to establish the finite branch-and-bound algorithm in Section 4.
is an optimal solution of (DNP). Ifx ∈ H, thenx is an optimal solution of problem (QPCC).
Proof. Ifx ∈ H andx is part of an optimal solution of (DNP), thenx is a feasible solution of (QPCC). It follows thatŶ = 1x T xxx T is a feasible solution of (DNP) andQ •Ȳ ≤Q •Ŷ . On the other hand, following from the factX −xx T ∈ S n + , Note thatQ •Ȳ is a lower bound on the optimal value of (QPCC), it follows thatx is an optimal solution of (QPCC).
Based on Property 1, we can design a finite branch-and-bound algorithm in which x ∈ H is enforced gradually. The finiteness is ensured by fixingx i = 0 orx j = 0 wheneverx ixj = 0 in the branching process.
3. An augmented Lagrangian method. In this section, we investigate the special structure of doubly nonnegative relaxation (DNP) in details and design the decomposition method for augmented Lagrangian function.
Lemma 3.1 allows us to restate problem (DNP) as the following formulation: In the next subsection, we will design a decomposition method based on this new formulation.

3.2.
Augmented Lagrangian and decomposition. Define the following closed and convex cone: Let K * be the dual cone of K. It is easy to verify that K = {Y = N P N T |P ∈ S m−r + }. Hence, for any given matrix R ∈ S n+1 , the projection of R into the cone K is (ref. [21]) and Now, we are ready to derive the augmented Lagrangian function. Introducing an auxiliary variable Z and consider the following auxiliary problem: where matrix U = 1 u T u uu T with u being the upper bounds on x implied by the assumption that feasible region F is bounded. The idea of the decomposition is to relax the constraint Z = Y and then apply the standard augmented Lagrangian method. We introduce the multiplier S and the penalty parameter σ > 0 to yield the augmented Lagrangian function where · F is the Frobenius norm. The associated subproblem is then (AL-Sub) After we get an optimal solution (Y, Z) of problem (AL-Sub), we update the multiplier S in the standard way [23] In fact, the constraint Z ∈ K implies that S must be in K * , otherwise the optimal value of problem (AL-Sub) is unbounded. In the next subsection, we will discuss how to solve problem (AL-Sub) in an approximate but efficient way, and the approach for calculating the lower bound of problem (DNP).

3.3.
Solving the subproblem (AL-Sub) and calculating lower bound. As in [24], we adopt the block coordinate descent over Y and Z for solving problem (AL-Sub). When Z is held constant, problem (AL-Sub) over Y becomes where a constant term Z 2 F is ignored in the objective function. Obviously, problem (Y-Block) can be immediately decomposed into n(n+1)/2−|E|−1 one-dimensional strictly convex quadratic problems of the form min y {αy + σy 2 /2|0 ≤ y ≤ µ}, whose solution is analytic. We remark that the upper bound µ is finite and positive σ ensures the boundedness of each one-dimensional subproblem.
On the other hand, when Y is held constant, problem (AL-Sub) over Z becomes Completing the square and dividing by σ, the objective function of problem So the optimal solution of this problem is just the projection of Y − σ −1 S onto K, which can be calculated according to (2). The mathematical operations for completing this projection include matrix production and a projection onto S m−r + . In total, it requires O(n 2 (m − r) + (m − r) 3 ) = O(n 2 (m − r)). We point out that this is the most dominant computational burden in the proposed algorithm and we implemented it in C programming language and connected to MATLAB using MEX library.
In order to calculate a valid lower bound v for problem (DNP), for a given multiplier S ∈ K * , we consider the standard Lagrangian relaxation of problem (DNP-New), i.e., setting σ to be 0 in problem (AL-Sub). In this case, problem (AL-Sub) is separable over Y and Z. Moreover, problem (Z-Block) has optimal solution Z = 0 because S ∈ K * . Hence, problem (Y-Block) with σ = 0 and Z = 0 can be decomposed into O(n 2 ) one-dimensional linear optimization problems of the form min y {αy|0 ≤ y ≤ µ} = min{0, αµ}. In all, the total cost for calculating a lower bound v requires O(n 2 ) mathematical operations.
We state the augmented Lagrangian method for problem (DNP) as Algorithm 1. Some details are not specified here, such as how to initialize the parameter σ and the termination criterion. The reason is that these strategies are closely related to the branch-and-bound method that will be adopted. Therefore, we only state a basic framework of the algorithm here and leave the discussion on the details in Section 4. We end this section with two remarks. Solve problem (Y-Block) to get Y .

4:
Solve problem (Z-Block to get Z.

A branch-and-bound algorithm.
In this section, we present implementation details about the branch-and-bound algorithm for globally solving problem (QPCC). We also prove the finiteness and correctness of the proposed algorithm under the mild assumption.
In each node, some variables of x i are fixed to be zero, where (i, j) ∈ E or (j, i) ∈ E for some j ∈ {1, ..., n}. The relaxation problem in this node is then x i = 0, ∀i ∈ I X ij = X ji = 0, ∀i ∈ I, j = 1, ..., n, where I is the set of indices of the variables fixed to be zero. It is obvious that Property 1 and Algorithm 1 still work for the above problem. The upper bound U = 1 u T u uu T is obtained in the preprocessing phase, in which a sequence of linear programs (LP) in the form of u i = arg max{x i |Ax = b, x ≥ 0} are solved for each variable x i , i = 1, ..., n. We used warm start to speed up the solution of these LPs. Specifically, the row and column basis information of previous LP is used in the next one.
The practical performance of Algorithm 1 depends immensely on how σ and S is initialized. We discuss the rule for updating and initializing these parameters. The default initialization is σ 0 = max i,j=1,...,n {|Q ij |, |q i |} and S 0 = 0. We calculate the lower bound v every 10 iterations, and then update σ according to the rule proposed by Burer [5]: wherev is the best lower bound found so far. If the ratio 1 + (v −v)/(1 + |v) is nonpositive, then we do not update σ. We terminate Algorithm 1 if (i) The relative change in the last 5 updates of v is less than 10 −5 ; (ii) the maximal iteration 10,000 is achieved, or (iii) the relative error Y − Z F /( Y F + Z F ) is less than 10 −5 . At the end of Algorithm 1, the optimal σ and S are saved as the initialization for the child nodes. When I is a vertex cover of E (that is x i x j = 0 for all (i, j) ∈ E), the node becomes a leaf, which has to be fathomed for the seek of correctness of a branchand-bound algorithm. The problem in the leaf node has the following form whereQ,q,Ā andb are the data obtained by deleting corresponding rows and columns indexed by I. Obviously, problem (Leaf) is a convex quadratic program, whose feasible region is a polyhedron with nonempty relative interior. In this case, we switch to use the interior-point method (IPM) to solve the leaf node exactly. The optimal value of problem (Leaf) provides a valid global upper bound (GUB) for (QPCC). A node is fathomed if (i) the lower bound obtained by solving the relaxation problem is within a prespecified tolerance of GUB; (ii) the solution of relaxation problem satisfies the condition of Property 1; or (iii) the relaxation problem is infeasible.
If a node cannot be fathomed, then branching step involves enforcing the complementarity constraints x i x j = 0 for (i, j) ∈ E by fixing some variable to be zero. We choose a maximum violation rule to select the complementarity constraint to branch on. Specifically, for the current optimal solutionx, we choose (ĩ,j) = arg max{x ixj |(i, j) ∈ E,x ixj > 0}. Then, two children nodes are created: one enforces xĩ = 0 while the other one xj = 0.
More details about the branch-and-bound algorithm are stated as follows: • When selecting a node to solve, we choose the one with the smallest lower bound. • In each node, the MATLAB function fmincon is used to find a local optimal solution. The global upper bound GUB is then updated if the optimal value is less than GUB. • We use the relative optimality condition for fathoming. That is, for a given tolerance > 0, a node with lower bound v is fathomed if (GUB − v)/ max{1, |GUB|} < , and Algorithm 1 is stopped before the termination criteria are met.
Remark 3. The fmincon function in MATLAB can find a good upper bound most of the time, thus helping to obtain a good GUB and eventually reduce the number of nodes. Our preliminary computational results show that the number of nodes increases dramatically if fmincon is not used. Other efficient solution methods that can find a good feasible solution of problem (QPCC) could be alternatives of fmincon for updating the GUB.
We are now ready to formally state the proposed branch-and-bound method for problem (QPCC) as Algorithm 2. The main steps of the algorithm are presented here and the details of each step can be found in previous paragraphs. Proof. The finiteness of Algorithm 2 is ensure by the branching rule because there are finite possible vertex covers for the set E. In order to have a correct branchand-bound algorithm, the relaxation problem should have no gap at any leaf node.

Algorithm 2 A Branch-and-Bound Algorithm for Globally Solving Problem (QPCC)
Input: Data (Q, q, A, b, E). Output: A global solution of (QPCC).
1: Preprocessing Step: Calculate the upper bound for variable x.

2: Initialization
Step: The list L to be explored is initialized to contain the root node. Upper bound is set to +∞, and lower bound is set to −∞. 3: while L is not empty do 4:

Node Selection
Step: Select and remove the best-first node from L.

5:
Bounding Step: Solve the relaxation problem of the node by Algorithm 1. Update upper and lower bound if possible. If the node is not fathomed, go to the next step.

6:
Branching Step: Branch the node on the most violated complementarity constraint, generate two children nodes and add them to the list L. 7: end while This is ensured in Algorithm 2 because the problem in the leaf node is a quadratic convex program, which can be solved exactly when feasible.
5. Numerical experiment. In this section, we test the proposed algorithm on various types of instances. The algorithm has been implemented in MATLAB R2015a on a PC with Windows 7, 3.60GHz Intel Dual Core CPU processors and 4GB RAM.

Instances on MacMPEC.
We first solved some QPCC instances from MacMPEC [16]. Note that the instances are not in the standard form we presented in this paper, thus we transformed them by adding some auxiliary variables and constraints. Table 1 displays the computational results for six instances. The first column "Id" shows the problem name on the web site, the second column shows the size of the transformed problem, the column "nodes" shows the nodes explored by the proposed algorithm, and the forth column "time" displays the total CPU time in seconds. We confirmed, by our approach, that the global optimal values given on the MacMPEC were correctly identified. Especially, we can solve the largest instance "flp4-4" on MacMPEC in about one minute. Random instances. We also tested the proposed algorithm on randomly generated instances. QPECgen [13] is an instance generator for QPCC problems. It allows the user to control the properties of the generated QPCC instances and gives stationary points, which are often the global solutions of the instances. Table  2 contains the computational results for a collection of QPCC instances of small to medium size. For each given problem size, ten instances are generated, and the average number of explored nodes and average CPU time required by the proposed algorithm are displayed in the second and third column, respectively. The results clearly show that the proposed algorithm can efficiently solve the randomly generated instances in a reasonable time. It is worthwhile to point out that the average number of explored nodes is still small, even if n = 100 and |E| = 40.
Here, the matrix Q is usually not positive semidefinite. This problem is notable for embracing a wide range of applications, including set partition, resource allocation, and job assignment problems. Although the form of BQP is simple, it has been proved to be NP-hard by Pardalos and Jha [20]. One of the solving methods for BQP is quadratic convexification reformulation (QCR) [3]. Specifically, by using the fact that x i = x 2 i for i = 1, ..., n, BQP can be reformulated as where Λ is a diagonal matrix and e is a vector of all ones. QCR tries to find a Λ such that Q+Λ is positive semidefinite and the continuous relaxation of (RBQP) provides the best lower bound for BQP. Lu and Guo [18] proposed another QCR method by choosing Λ such that the trace of Λ is minimized while keeping the positive semidefiniteness of Q + Λ. Our preliminary experiment shows that the proposed algorithm works better with Lu and Guo's QCR method. Thus, we adopted their QCR method in this numerical experiment. Note that problem (RBQP) can be transformed to the form of (QPCC) by introducing an auxiliary variable y = e − x and adding complementarity constraints x i y i = 0 for i = 1, ..., n. The data sets used for the experiment are from the OR-Library [2]. All instances have 10% density and the elements of Q are integral values in the [-100, 100] interval. The problem name, optimal value, the number of explored nodes and CPU time are given in Table 3. The proposed algorithm can solve 8 out of 10 instances within half an hour. For the rest 2 instances ("bqp50-1" and "bqp50-10"), the algorithm terminated due to the time limit of 1,800 seconds, and the best upper bound, instead of optimal value, is reported in Table 3. It is notable that, for instance "bqp50-3", the number of explored nodes is only 69 while it took 114 iterations for CPLEX 12.6 to solve the same instance. 6. Conclusion. In this paper, we propose a branch-and-bound algorithm based on the doubly nonnegative relaxation for quadratic program with complementarity constraints. The doubly nonnegative relaxation in each node is efficiently solved by an augmented Lagrangian method due to the fact that the Lagrangian function can be decomposed into two subproblems that are easy to solve. The computational results show that the proposed algorithm can solve problems of small and medium size within a reasonable amount of time, with a particular highlight in a few nodes for some instances.