A range division and contraction approach for nonconvex quadratic program with quadratic constraints

This paper presents a novel range division and contraction approach for globally solving nonconvex quadratic program with quadratic constraints. By constructing new underestimating linear relaxation functions, we can transform the initial nonconvex quadratic program problem into a linear program relaxation problem. By employing a branch and bound scheme with a range contraction approach, we describe a novel global optimization algorithm for effectively solving nonconvex quadratic program with quadratic constraints. Finally, the global convergence of the proposed algorithm is proved and numerical experimental results demonstrate the effectiveness of the proposed approach.

programming problems can be converted into this mathematical modeling, and the solutions of many nonlinear optimization problems can be approximated or obtained by solving a sequence of nonconvex quadratic programs with quadratic constraints. Moreover, from research point of view, since the nonconvex quadratic programs with quadratic constraints possess multiple local minimum points which are not globally minimum point, they exist significant computational and theoretical challenges. Therefore, it is necessary to put forward a new global optimization approach for solving the nonconvex quadratic program with quadratic constraints.
In past several decades, various algorithms have been developed for solving the nonconvex quadratic program with quadratic constraints and its special form, which are given as follows. Based on Newton method, branching rule and cutting plane, Vandenbussche and Nemhauser presented a branch-and-cut method for nonconvex quadratic programs with box constraints (Vandenbussche and Nemhauser 2005). Used different d.c. decomposition method to construct the relaxation of quadratic objective function, and used the "optimal level solution" parametrical approach to solve relaxation problem, Cambini and Sodini (2005) proposed a decomposition method for solving nonconvex quadratic programs over a compact polyhedral feasible region. By decomposing a large-scale quadratic programming into a serial of small-scale ones and approximating the solution of the large-scale quadratic programming via the solutions of these smallscale ones, Li and Zhang (2006) presented a decomposition algorithm for solving largescale quadratic programming problems. Through decomposing quadratic objective function into a separable equivalent function, then by constructing linear under-estimator of the corresponding objective function,  presented a decomposition and linearization method for globally solving quadratic programs with linear constraints. Vavasis (1992) presented an approximation algorithm for indefinite quadratic programming, and concluded that such an approximation solution can be found in polynomial time. Based on D.C. decomposition, Cholesky factorization and convex relaxation, Yajima and Fujie (1998) proposed a decomposition-and-relaxation algorithm for the general quadratic problems with box constraints. Ye (1992) proposed an affinescaling algorithm to solve nonconvex indefinite or negative definite quadratic programs problems. By utilizing cutting plane technique, Gao and Deng (2008) presented a branch and bound method mixed with cutting plane technique for solving concave quadratic programming problems. Using lagrangian underestimates and interval Newton method, Voorhis (2002) proposed a global optimization algorithm for quadratic programs. Based on simplicial branch-and-bound scheme, Raber (1998) presented a simplicial branchand-bound method for solving nonconvex quadratic programs. Based on parametric linear relaxation and new linearizing technique, Jiao et al. (2015), Jiao and Chen (2013) proposed two branch and bound algorithms for globally solving nonconvex quadratic programs. Using duality bounds and branch-and-bound scheme, Thoai (2000) presented a duality bound method for the general quadratic programming problem with quadratic constraints. Based on linear relaxation method and branch-and-bound framework with rectangle reducing technique, Shen and Liu (2008), Gao et al. (2005a, b), Jiao et al. (2014) proposed four branch-and-reduce algorithms for solving nonconvex quadratic programs, respectively. Based on branch-and-bound scheme, An and Tao (1997), An and Tao (1998) presented two D.C. algorithms for solving nonconvex quadratic programs. Apart from the above reviewed quadratic programs algorithms, many algorithms (see Shen and Jiao 2006;Wang and Liang 2005;Wang et al. 2004;Shen 2005;Shen and Li 2013;Shen and Bai 2013;Shen et al. 2009) for solving geometric programming can be also used to solve the nonconvex quadratic program with quadratic constraints proposed in this paper. In addition, some recent artificial intelligent optimization algorithms (Zhang et al. 2013(Zhang et al. , 2014a(Zhang et al. , b, 2016 are developed, which are also used to obtain local optimal solution of the nonconvex quadratic program with quadratic constraints.
The above most deterministic algorithms for solving nonconvex quadratic programs are based on relaxation technique and branch-and-bound scheme, since the exhaustiveness of branching rule leads to a significant increase in the computational burden for solving nonconvex quadratic programs, until today, there is short of more effective algorithm for solving such problems, so it is necessary to establish a good algorithm for the NQPQC. Therefore, the main motivation for the paper is to construct a novel linearized technique and a range contraction approach, and based on these techniques develop an effective algorithm for solving the NQPQC.
In this paper, a novel range division and contraction approach will be proposed for globally solving the nonconvex quadratic program with quadratic constraints. The main features of the proposed algorithm in this paper are as follows. Firstly, a novel linearized technique is constructed for systematically converting the NQPQC into a sequence of linear program problems, and the solutions of the converted linear program problems are used to approximate the solution of the NQPQC, which can be computed by using simplex approach. Secondly, in order to accelerate the running speed of our approach, a range contraction method is constructed and employed in our algorithm. Next, use branch-and-bound framework with a range contraction approach, a new global optimization algorithm is established. Finally, the global convergence of the proposed approach is proved and numerical experimental results demonstrate the computational efficiency of our algorithm.
The remaining contents of the paper are stated as follows. Second section describes a new linearized method and the linear program relaxation problem of the initial NQPQC is derived. In third section, a new range division and contraction approach is established and its global convergence is proved. Fourth section describes numerical results for some examples, which be computed by using the proposed algorithm. Finally some conclusions are drawn.

New linearized approach
In this section, a new linearized approach is constructed for establishing the linear program relaxation problem of the NQPQC. Without loss of generality, let X k = {x ∈ R n | l k = (l k 1 , . . . , l k n ) T ≤ x ≤ u k = (u k 1 , . . . , u k n ) T } ⊆ X 0 , and let i be the minimum eigenvalue of Q i . Set In the following, for any i ∈ {0, 1, . . . , m}, for any x ∈ X k , define where ρ is an arbitrary positive real number.
Theorem 1 For any x ∈ X k , consider the functions F i (x) and F L i (x), we have: Proof (1) Consider the function x 2 j over [l k j , u k j ], according to the mean value theorem, we can get . . , n. By the expressions of the functions F i (x) and F L i (x), we can get that Therefore, we have Since therefore, we get The conclusion is drawn. From the above Theorem 1, we can establish the linear program relaxation problem (LPRP) of the NQPQC in X k as follows.
where Obviously, by the construction method of linear program relaxation problem (LPRP), we get that the feasible region of the LPRP contain all feasible points of the NQPQC in X k , and the optimal value of the LPRP is not more than that of the NQPQC in X k .

Range contraction approach
To accelerate running speed of our algorithm, a range contraction approach is formulated as the following Theorem 2. The proposed range contraction approach aims at contracting the investigated rectangle X without pruning any global optimum point of the initial NQPQC.
Without loss of generality, for any , we can rewrite the function F L i (x) as follows: Let UB be a known currently upper bound of the global optimum value of the NQPQC(X 0 ), and set Theorem 2 For any sub-rectangle X ⊆ X 0 , the following conclusions hold: . . , m}, then, for any p ∈ {1, 2, . . . , n}, if c ip > 0, then there does not contain global optimal solution of the NQPQC(X 0 ) over X = ( X j ) n×1 ; if c ip < 0, there does not contain global optimal solution of the NQPQC(X 0 ) over X = ( X j ) n×1 .
Proof The proof of the Theorem can be similarly given by the Theorem 3 in Jiao et al. (2014), therefore, here is omitted.
From the above Theorem 2, a new range contraction approach is presented as follows:

Range contraction approach
Without loss of generality, assume that X = (X j ) n×1 with X j = [l j , u j ] (j = 1, . . . , n) be any sub-rectangle of X 0 , and for each i = 0, 1, . . . , M, compute LB i . If there exist some i ∈ {1, . . . , M}, such that LB i > b i , then delete the whole rectangle X; Otherwise, for each i ∈ {1, . . . , M}, p ∈ {1, . . . , n}, compute If LB 0 > UB, then delete the whole rectangle X; Otherwise, for each p ∈ {1, . . . , n}, compute UB− LB 0 +min{c 0p l p ,c 0p u p } c 0p . If c 0p > 0 and UB− LB 0 +min{c 0p l p ,c 0p u p } c 0p < u p , then let u p = UB− LB 0 +min{c 0p l p ,c 0p u p } c 0p ; else if c 0p < 0 and UB− LB 0 +min{c 0p l p ,c 0p u p } c 0p > l p , then let l p = UB− LB 0 +min{c 0p l p ,c 0p u p } c 0p . From Theorem 2, by utilizing the above range contraction approach to compress the investigated rectangle region, or delete a part of the investigated rectangle region, in which there does not contain the global optimum solution of the NQPQC. Therefore, we can improve the computational speed or computational efficiency of the algorithm.

Range division and contraction algorithm
In this section, a new range division and contraction algorithm is presented for globally solving the initial NQPQC. In the algorithm, one of the most important operation is the choice of a suitable range division method. Here, we choose an w − division approach, which is sufficient to ensure the global convergence of the proposed algorithm, this is because that the selected range division method drives all interval to zero for all variables. This range division method is described as follows.
For any investigated rectangle X where ω ∈ (0, 1). At the same time, all other interval remain unchanged. Consequently, the investigated rectangle X ′ is partitioned into two new sub-rectangles X ′ 1 and X ′ 2 Without loss of generality, we denote LB(X k ) and x k = x(X k ) as the optimum value and optimal solution of the linear program relaxation problem over the sub-rectangle X k , respectively. Combining the former linear program relaxation problem, range division method and range contraction approach together, a new global optimization method is constructed for effectively solving the NQPQC, the main steps of the proposed algorithm are described as follows.

Range division and contraction algorithm
Step 0. (Initializing) Set the initial iteration number k = 0, the initial active node set 0 = {X 0 }, the given termination error ǫ > 0, the initial upper bound UB 0 = +∞ and the initial feasible point set = ∅, respectively.
Step 2. (Division) Use the proposed range division method to divide X k into two new sub-rectangles, and denote the set of new partitioned sub-rectangle as X k .
Step 3. (Contraction) For each investigated sub-rectangle X ∈ X k , use the presented range contraction approach to compress its range, and still denote the remaining rectangle part and the remaining partitioning set by X and X k , respectively.
Step 4. (Bounding) If X k is not empty, for each X ∈ X k , compute the LPRP(X), and denote its optimum value and optimum solution by LB(X) and x(X), respectively. If LB(X) > UB k , let X k := X k \ X; else if the midpoint x mid of X is a feasible solution of the NQPQC(X 0 ), then let := ∪ {x mid }, and if x(X) is a feasible solution of the NQPQC(X 0 ), then let � := � ∪ {x(X)}.
If is not empty, denote the new upper bound UB k := min x∈� F 0 (x), and denote the best feasible point x * := argmin x∈� F 0 (x).
Denote the new remaining partition set and the new lower bound by � k := (� k \X k ) ∪ X k and LB k := inf X∈� k LB(X), respectively.
The above algorithm either terminates after finite iteration or generates an infinite iteration sequence, from the exhaustiveness of the used range division approach, we can follow that all intervals of all variables must shrink to a singleton, i.e., �u k − l k � → 0. At the same time, the Theorem 1 guarantee that the linear program relaxation problem (LPRP) (X k ) will infinitely approaches the problem NQPQC(X k ) as �u k − l k � → 0.

Theorem 3
The above algorithm either finishes finitely at the global optimum point x * of the initial NQPQC, or produces an infinite iteration sequence {x k }, of which any accumulation point will be the global optimum point of the initial NQPQC.
Proof If the above algorithm finishes finitely at some iteration k, obviously, when the algorithm ends, we can get UB k = v * = LB k . Therefore, from step of the algorithm, we get that x k is a global optimum point of the initial NQCQP.
If the above algorithm does not finish at finite step, then it must produce an infinite sub-rectangle iteration sequence {X k }, from the exhaustiveness of division approach, we get that the sub-rectangle sequence {X k } must converge to a point. From the characteristics of our algorithm, we can get that {UB k } and {LB k } are nonincreasing and nondecreasing sequences, respectively. So that {UB k − LB k } is a monotonic non-increasing sequence. By the conclusion of Theorem 1, we get that the sequence {UB k − LB k } must be convergent to zero. From the structure of our algorithm, for each k, we have LB k ≤ v * ≤ UB k . Therefore, we can get that lim k→∞ UB k = lim k→∞ LB k = v * . From steps of our algorithm, we know that x k is always a feasible point of the initial NQCQP, we can get UB k = G 0 (x k ). By the continuity of constrained functions, we get that the limitation point x * of {x k } is also a feasible point of the initial NQCQP, and we can get the global optimum value v * = G 0 (x * ). Hence, the conclusion is proved.

Numerical experiments
In this section, in order to verify the performance of the algorithm proposed in this paper (our algorithm), several numerical examples in recent references are implemented on a Intel(R) Core(TM)2 Duo CPU microcomputer. Although these numerical examples have a relatively small number of variables, they are also quite challenging. The proposed algorithm program is coded in C++ language, all linear program relaxation problems are solved by using simplex approach, and the convergence errors are all set as ǫ = 10 −6 in our experiments. For all examples, the numerical results of optimal solutions, optimal values and number of iterations obtained by our algorithm and other approaches (Jiao and Chen 2013;Thoai 2000;Shen and Liu 2008;Gao et al. 2005b;Jiao et al. 2014;Shen and Jiao 2006;Wang and Liang 2005;Wang et al. 2004;Shen 2005;Shen and Li 2013) are illustrated in Table 1. The numerical experimental results show that our algorithm can globally solve the NQCQP problem. In Table 1, the notation "Iter. " represents "number of iterations". Example 1 (Jiao and Chen 2013;Thoai 2000).
For Example 1, from numerical results in Table 1, compared with the algorithms in Jiao and Chen (2013), Thoai (2000), use the same logic of the algorithm, our algorithm can obtain the same optimal solution (5.0, 1.0), but our algorithm spends less number of iterations.
For Example 2, by Table 1, compared with the algorithms in Shen and Liu (2008), Jiao et al. (2014), use the same logic of the algorithm, our algorithm can obtain the same optimal solution with less number of iterations; and compared with the algorithms in Wang and Liang (2005), our algorithm can obtain the better optimal solution with less number of iterations.
For Example 4, by the numerical results in Table 1, compared with the algorithms in Jiao and Chen (2013), Gao et al. (2005b), Jiao et al. (2014), our algorithm can obtain better or at least as good as optimal solution and optimal value with less number of iterations.
For Example 5, by the numerical results of Table 1, compared with the algorithms in Jiao and Chen (2013), Shen (2005), our algorithm can obtain the same optimal solution (1.5, 1.5), but our algorithm spend less number of iterations.
For Example 6, from the numerical results in Table 1, our algorithm obtains the optimal solution (1.177124344, 2.177124344) after 24 iterations, but the algorithm in Shen and Jiao (2006) obtains the optimal solution (1.177124327, 2.177124353) after 24 iteration, obviously our algorithm has higher computational efficiency than the algorithms in Shen and Jiao (2006); compared with the algorithm in Jiao and Chen (2013), our algorithm can obtain the same optimal solution (1.177124344, 2.177124344), but our algorithm spend less number of iterations.
Example 7 (Jiao and Chen 2013;Shen and Jiao 2006). For Example 7, from the numerical results in Table 1, compared with the algorithms in Jiao and Chen (2013), Shen and Jiao (2006), our algorithm can obtain the same optimal solution (2.0, 1.0), but our algorithm spend less number of iterations, obviously our algorithm has higher computational efficiency than the algorithms in Jiao and Chen (2013), Shen and Jiao (2006).
For Example 8, from the numerical results in Table 1, our algorithm obtains the optimal solution (1.0, 0.181818133, 0.983332175) and the optimal value −11.363635790 after 352 iterations, but the algorithm in Jiao and Chen (2013) obtains the optimal solution (1.0, 0.181818470, 0.983332113) and optimal value −11.363636364 after 420 iteration, and the algorithm in Shen and Li (2013) obtains the optimal solution (0.998712, 0.196213, 0.979216) and optimal value −10.35 after 1648 iteration, obviously our algorithm not only has higher efficiency but also get the better optimal solution and optimal value. In all, from numerical results for Examples 1-8, compared with Jiao and Chen (2013), Thoai (2000), Shen and Liu (2008), Gao et al. (2005b), Jiao et al. (2014), Shen and Jiao (2006), Wang and Liang (2005), Wang et al. (2004), Shen (2005), Shen and Li (2013), the presented algorithm in this paper can globally solve nonconvex quadratic program with quadratic constraints.

Conclusion
In this article, a new range division and contraction algorithm is proposed for globally solving nonconvex quadratic program with quadratic constraints (NQPQC). The linear program relaxation problem of the initial NQPQC is constructed by utilizing new linearizing method, which is derived by underestimating all quadratic objective function and constraint functions with linear relaxation functions. By applying the current upper bound and linear program relaxation problem of the NQPQC, a range contraction technique is introduced for improving the computational speed of our algorithm. By successive partition of the initial rectangle region, and by subsequently solving a sequence of linear program relaxation problems, the proposed algorithm converges to the global optimum point of the original NQPQC. Finally, Numerical computational results demonstrate the effectiveness and robustness of our algorithm.

Authors' contributions
This work was carried out in collaboration with all authors. All authors have a good contribution to design the new range division and contraction algorithm, and to perform the numerical analysis of this research work. All authors read and approved the final manuscript.