Binary quadratic programming (BQP) problem is a kind of typical combinatorial optimization problem, and has a variety of applications in computer aided design, traffic management, cellular mobile communication frequency allocation, operations research and engineering. And many combinatorial optimization problems with constraints can be transformed into BQP problems form by certain transformation, so these problems can be on behalf of kinds of important problems in combinatorial optimization. Hammer  pointed out that any integer programming problems where the objective function is quadratic or linear and the constraint condition is linear can be described as BQP problems. Glover et al.  have successfully used this transformation method to transform the quadratic knapsack problem into BQP problem for solving. Considering the feasibility of this transformation, BQP problem has more extensive practical application. Due to the wide application background of this problem and the difficulties caused by NP (Non-Deterministic Polynomial) hard properties, it is of great academic value to study effective algorithms for solving such problems.
Binary quadratic programming can be uniformly written in the following general
where Q is a real symmetric matrix, and the objective function can be followed by another linear term, but it can be omitted because there is no substantial influence on the method introduced later. Of course, in the following discussion, we will mainly focus on the so-called 0 - 1 quadratic programming:
However, the continuous method proposed later in this paper is also applicable to the general binary quadratic programming (1), which can be easily accomplished by simple variable transformation. It is well known that 0 - 1 programming problem will increase exponentially with the increase of the scale of the problem, and within limited computer resources and time the traditional solution method will be difficult to realize. Because of the wide application prospect and difficult characteristic of 0 - 1 programming, how to solve this kind of combinatorial optimization problem effectively has been the focus of many scholars. The classical methods for solving 0 - 1 programming problems are mainly divided into two types: accurate algorithms and random search algorithms. The accurate algorithms include implicit enumeration method, branch and bound algorithm, cutting plane algorithm, dynamic programming method, etc. Random search algorithms include genetic algorithm, simulated algorithm, artificial neural network algorithm, etc.      . But continuous method has become a new tendency of the 0 - 1 programming problem in recent years, whose advantage is to avoid inherent characteristics of the combinatorial optimization problem, and no longer limited to the size of the problem. This method is to transform the combinatorial optimization Aproblem into equivalent continuous optimization problem for solving, thus effectively avoid the so-called combination explosion problem, which can solve some large combinatorial optimization problems efficiently. While this approach is still in the exploratory stage with few research work in this area, but we think it is a noteworthy new research direction. Interested readers can see the literature  -  .
In this article, we will first turn the binary constraints of the 0 - 1 programming problem into equivalent complementarity problem, that is to say, a mathematical programming problem with equilibrium constraints (MPEC); then NCP function is used to transform the complementary condition into equivalent equation. Finally, the non-smooth equation is smoothed by using the aggregate function and we finally get a nonlinear smooth optimization problem.
2. Continuous Formulation of BQP
The simplest and most intuitive way to solve 0 - 1 planning problem is to adopt the round integration technology, that is, to consider 0 - 1 variable as a continuous variable, and to replace the variable constraint in the original problem with the following interval constraint:
Then a relaxation optimization problem is obtained, and the components of the relaxation solution are rounded to the nearest discrete value. Although this method is easy to implement, but the theory of this technology is not strict. It cannot guarantee that the optimal solution obtained is the global optimal solution, or even a feasible solution.
The other method is the penalty function method, which is easy to see that is always equivalent to a complementarity condition:
Therefore, the punishment function can be constructed as follows:
where is a large enough penalty parameter. Adding this penalty function to the object function can obtain an equivalent continuous optimization problem:
Due to the strong concavity of , the object function in (6) is concave. The equivalence of (2) and (6) is based on the fact that the concavity function obtains the minimum value at a vertex, as well as has included and the sum is included. Although the vertex of the feasible domain is not necessarily the vertex of the unit hypercube, if is sufficiently large, the global minimum can only be obtained while .
However, the specific value of cannot be determined. If the value is too small, the penalty term will not play its role. If the value is too large, due to the strong concavity of the function, this penalty function method cannot effectively obtain the optimal solution of the original problem, so it is not an effective method. It is pointed out in literature  , but only a new viewpoint is provided, and no more effective method is given to solve the original problem. Therefore, it is not ideal to solve the original problem directly with (6). Great-Ming Ng later constructed a logarithmic barrier term:
This barrier term can be used to replace the boxed constraint conditions as well as to avoid falling into the local minimum during the iteration, so as to find the global optimal solution for the original problem. At this point, the following unconstrained optimization problems can be obtained:
where is a penalty parameter and is a barrier parameter. In a sense, (8) is an improvement on (6), but in essence they both replace constraints with penalty terms or logarithmic barrier terms.
For the BQP problem, another continuous method is proposed in this paper. Consider the following two sets of constraints:
Obviously, these two sets of constraint conditions and binary constraint conditions are equivalent. In fact, these two constraints skillfully reform complementary constraints:
The advantage of constructing constraint conditions in this way is that: the two constraint conditions of (9) work at the same time, and to each variable there is a unique complementary constraint condition
And more importantly, complementary constraints (10) and the general constraint conditions are different: for general complementary constraints, each component is a function about , and for (10), the value of every only has relationship with the values of , nothing to do with the other variables. Namely, each of the complementary constraint of (10) is independent of each other.
Through the above, the problem (2) can be equivalently transformed into the following mathematical programming problems with complementary constraints:
For the above complementary constraints, we can use NCP function  to replace them further. To better describe NCP functions, we first introduce the concept of them:
Definition (NCP function): If a function is set , and , then the function of is called an NCP function.
Therefore, we can replace the complementary constraint in (11) with a series of equation equations. Finally, we get the continuous optimization model as follows:
Thus, the global optimal solution of (12) is the exact solution of problem (2).
3. Smoothing Method for BQP
According to the definition of NCP function given above, functions satisfying conditions can take various forms. Several commonly used NCP functions can be given:
In the following section, we mainly choose the minimum value function to study.
It is easy to know , and for the non-differentiable maximum function:
we often use aggregate function to carry out smooth approximation to it   . The so-called aggregate function is a function in the following form:
where the smooth parameter is large enough, and are continuous differentiable (i.e., smooth) real value function. The aggregate function has very good properties  , which can be found from the following properties:
Lemma 3.1 The function which is defined by (15) satisfies the following conditions:
Proof 1) can be equivalent to deformation , because , so , and there is at least one indicator that makes , so exist . So we have
2) When , the above conclusion can be proved by formula (16).
Based on the above introduction, we can handle the minimum function as follows  :
where the smooth parameter is large enough, as shown by lemma 3.1, when , and are equivalent.
For ease of expression, let’s call:
At this point, the problem (12) can be further transformed into an equivalent continuous and smooth nonlinear constraint optimization problem:
By this continuous method, the global optimal solution of the original problem (2) can be solved by solving the corresponding nonlinear constrained optimization problem (19). Many mature optimization algorithms have been developed for solving nonlinear equality constraint optimization problems. This paper mainly uses the multiplier penalty function method to solve the above problems (19). Multiplier method is an optimization algorithm independently proposed by Powell and Hestenes in 1969 to solve equality constraint optimization problems. Later, it was extended to solve inequality constraint optimization problems by Rockfellar in 1973. The basic idea is to start from the Lagrangian function of the original problem and add an appropriate penalty function, so as to transform the original problem into a series of unconstrained optimization subproblems.
For convenience’s sake, let’s first:
The Lagrange function of problem (19) is  :
where is the Lagrangian multiplier vector, and is set as KT pairs of problem (19), then it can be known from the optimality condition:
In addition, it is not difficult to find any x in the feasible domain is satisfied:
The above equation shows that if the multiplier vector is known, the problem (19) can be equivalently transformed into:
Then the external penalty function method is considered to solve the problem (24). The augmented objective function is
In this way, we can fix and find a minimum of , then change the value of appropriately to find a new value , until we get the and that we want. Specifically, when solving the minimum of in the k-th iteration of the unconstrained subproblem, the necessary conditions for taking the extreme values are known:
And the KT-point of the original problem satisfies:
In order for and , after comparing the above two expressions, the updating formula of the multiplier sequence is:
It can be seen from Equation (28) that the sufficient and necessary condition for convergence is And then we proof that is also a necessary and sufficient condition for judging the KT-point .
Theorem 3.1 Let be the minimum point of the unconstrained optimization problem:
then being a KT-point to the (19) has a necessary and sufficient condition of .
Proof The necessity is obvious. The following proof is sufficient, since is a minimum of (29) and , therefore, for any feasible point:
that is, is also a minimum of (19). On the other hand, it is noted that is also the stable point of (29), therefore
The above formula indicates that is the Lagrangian multiplier vector with respect to , that is to say is also the KT-point of (19).
Based on the above discussion, before giving the detailed steps to solve the multiplier method of equation constraint problem (19), the following lemmas are given by studying some properties of .
Lemma 3.2 For each , the following formula holds:
Proof According to the properties of aggregate function:
thus . That is to say, . Therefore, .
Lemma 3.3 For any , there is
Proof Easy to obtain by lemma 3.2 that: , so
From the above equation we know:
Lemma 3.4 For each , is strictly convex on the interval ; When , is convex on the interval .
Proof According to the definition of , it can be seen that:
let . We could get,
So is strictly convex on the interval ; for , when all the satisfy , so when , is convex on the interval .
Lemma 3.5 For any and that are large enough, as long as the region where x satisfies the following equation: , the augmented Lagrange function defined by Equation (25) is convex in this region.
Proof Let then the hesse matrix of is
So we just have to satisfy , the matrix is positive, therefore for sufficiently large , is convex in the region that satisfies the above conditions. When the interval is is convex on the interval ; when the interval is .
According to lemma 3.5, we can use the augmented Lagrange penalty function to solve the problem (19). For a strictly monotonic increasing sequence
and , the solution of unconstrained optimization is
, the corresponding Lagrange multiplier is , and is modified according to (28) in the iteration.
Based on the previous analysis, the basic algorithm for solving the problem (19) is as follows  :
Step 1 Given parameters , , , and , . Select a starting point and set , k = 0;
Step 2 Using BFGS algorithm solve the unconstrained minimization problem
with the starting point , and denote by its optimal solution;
Step 3 If , set go to Step 6; else go to Step 4;
Step 4 If go to Step 5; else update the parameter , set , and go to Step 1;
Step 5 Set , , modify according to Equation (28), set , and go to Step 1;
Step 6 Finish.
The optimal solution sequence required by the above algorithm makes converge to 0 at a speed of at least 0.1. If this requirement is not met in an iteration, the penalty factor and smooth parameter should be increased automatically.
Then we introduce the BFGS algorithm mentioned in step 2. It is the most popular and effective quasi-newton methods at present. It was proposed by Broyden, Fletcher, Goldfarb, and Shanno independently in 1970, so it is called BFGS algorithm. We know that the basic idea of quasi-newton method is to replace the Hesse matrix with one of its approximate matrixs in the steps of basic Newton method, and the following relational expression is required:
where displacement , gradient difference , (42) is usually referred to as quasi-newton equation or quasi-newton condition.
Algorithm 2 (BFGS)
Step 1 Given parameters , select a starting point and positive matrices ( usually set as or the unit matrix ),and set , ;
Step 2 Calculate , if go to Step 6, output as approximate minimum point; else go to Step 3;
Step 3 Solve the linear equations and get the solution ;
Step 4 Suppose is the minimum non-negative solution that satisfies the following inequality:
set , ;
Step 5 Set ,by the correction formula:
identify the , set and go to Step 1;
Step 6 Finish.
5. Number Experiments
About BQP questions, the current relatively popular in the question bank is J. E. Beasleys OR-library (http://people.brunel.ac.uk/~%20mastjjb/jeb/orlib/bqpinfo.html), in this paper, we has carried on the numerical experiment on parts of question bank BQP questions, and the calculation results of our continuous algorithm are compared with those of the best solution of known, as shown in Table 1.
The first column Mu in the Table 1 is the question number; The second column m denotes the number of variables; The third is the numerical result of this paper; The fourth is the best known solution; The fifth column represents our numerical result as a percentage of the best solution.
Table 1. Comparison of numerical results.
It’s not hard to see from the Table 1 that in the process of solving these problems, the optimal solution obtained by our algorithm is basically close to the best known solution. Through the comparison in the table above, it can be seen that the proposed continuous algorithm is basically close to the known best solution for solving the general unconstrained 0 - 1 quadratic programming problem.
In this paper, we have reformulated an unconstrained BQP problem as a MPEC problem by the equivalent complementarity conditions of a binary vector. To seek a global minimizer of the resulting continuous optimization problem, we construct a global smoothing function and develop a global continuation algorithm via a sequence of unconstrained minimization. And the numerical results indicate that the continuous approach proposed is extremely promising, especially for those large problems, in terms of the quality of the optimal values generated and the computational work involved. In addition, this new approach can be extended to general nonlinear binary optimization problems, and we can leave it as a future research topic.
Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.
 Wolkowicz, H. and Anjos, M.F. (2002) Semidefinite Programming for Discrete Optimization and Matrix Completion Problem. Discrete Applied Mathematics, 123, 513-577.
 Glover, F. (2002) One-Pass Heuristics for Large-Scale Unconstrained Binary Quadratic Problems. European Journal of Operational Research, 137, 272-287.
 Pardalos, P.M. (1996) Continuous Approaches to Discrete Optimization Problems. Nonlinear Optimization and Applications. Plenum Publishing, New York, 313-328.