Journal of Global Optimization 28: 1–15, 2004. © 2004 Kluwer Academic Publishers. Printed in the Netherlands. 1 Ellipsoidal Approach to Box-Constrained Quadratic Problems PASQUALE L. DE ANGELIS1 , IMMANUEL M. BOMZE2 and GERARDO TORALDO3 1 Institute of Statistic and Mathematics, University of Naples ‘PARTHENOPE’, Naples, Italy Department of Statistics and Decision Support Systems, University of Vienna, Viennna, Austria 3 CPS-CNR & Department of Agricultural Engineering and Agronomy, University of Naples ‘Federico II’, Naples, Italy 2 (Received 19 September 2001; accepted in revised form 19 November 2001) Abstract. We present a new heuristic for the global solution of box constrained quadratic problems, based on the classical results which hold for the minimization of quadratic problems with ellipsoidal constraints. The approach is tested on several problems randomly generated and on graph instances from the DIMACS challenge, medium size instances of the Maximum Clique Problem. The numerical results seem to suggest some effectiveness of the proposed approach. Key words: Nonconvex quadratic programming, Ellipsoidal constraints, Global optimization 1. Introduction We consider global optimization problems of the form 1 (1.1) min
f x = x⊤ Ax +b ⊤ x  x ∈ −1 1 n 2 where A = aij ∈ IRn×n is a symmetric n×n matrix and x b ∈ IRn with x⊤ denoting the transposed vector. Problems of such form appear in a number of applications. Moreover (1.1) arises as subproblem in the solution of more complex optimization problems. For the convex case, which is known to be polynomially solvable [12], many different approaches have been proposed. Also for the concave case, which is equivalent to the 0–1 quadratic problem, several algorithms exist based on the combinatorial nature of that problem. Here public domain packages like Q01SUBS [18] are available. However, there are very few algorithms for the general indefinite case, where the matrix A has eigenvalues of mixed sign (see, for example, a recent paper [20] or [7] for a survey on bound constrained quadratic programming). In general only problems with few variables can be handled with exact algorithms, and therefore, for medium-size and large problems, heuristics have to be devised, which iteratively produce a sequence of local solutions. Obviously, the hard task calling for heuristic approaches is to escape from an inefficient local solution towards an improving one. The paper is organized as follows. In Section 2 we first review some results taken from the trust region literature, about global minimization of a quadratic 2 PASQUALE L. DE ANGELIS, IMMANUEL M. BOMZE AND GERARDO TORALDO function over a ball. The key result for the purposes addressed in this paper is apparently new: the Lagrange multiplier for the global solution (which is known to be unique) depends continuously, and in an decreasing manner, on the radius of the ball. In fact, this dependence is even smooth, with the — possible, but from a generic viewpoint highly improbable — exception of one kink. We also show that a global solution can be chosen such that it depends continuously (again, generically, even smoothly) on the radius. Note that in degenerate cases, the global solution need not be unique. In Section 3 we present a heuristic algorithm for the problem (1.1), in which a local solver is used, using a ‘promising’ starting point whose computation is motivated by the results in Section 2. Such an algorithm can naturally be embedded in a branch and bound strategy. Finally, Section 4 presents some computational results on a wide variety of test problems. In particular the algorithm is tested on small and medium sized graph instances from the DIMACS Challenge. 2. Local and global minimizers for quadratic functions over a ball In this section we study the auxiliary problem 1 (2.2) min
f x = x⊤ Ax +b ⊤ x  x ∈ IRn with x⊤ x
 2 It is well-known in the trust-region literature, see, e.g., [14, 13, 19], that there is at most one local solution to (2.2) which is not global, and generically, the global solution to (2.2) is unique. Moreover, we exhibit a way how to obtain a global solution which depends continuously on the ball radius, or, equivalently, on  (as √ ⊤ usual, we denote by x = x x). It will turn out that this dependence is in fact smooth, with the possible exception of a single kink at some critical crit in degenerate cases. Similarly, the corresponding Lagrange multiplier depends (again, with the possible exception of one kink at crit ) smoothly and in a monotonically decreasing way on the radius. While most of the used technology is folklore by now, the monotonicity and smoothness results seem to be new, and in order to establish them we anyhow need the prerequisites which were used previously. Thus we decided to present also short proofs of known results, to make the presentation self-contained, and in particular to shed some light on the genericity conditions which exclude degeneration of the problem. We start with some second-order optimality conditions for Karush/Kuhn/Tucker (KKT) points x̄ of (2.2). Recall that such a point satisfies (i) A+In x̄ +b = 0 (ii) x̄⊤ x̄ − = 0 with (iii) x̄⊤ x̄
 and   0 (2.3) The last conditions (iii) are the primal (for x̄), and dual (for the Lagrange multiplier ) feasibility conditions, respectively, while (i) and (ii) represent first-order ELLIPSOIDAL APPROACH TO BOX-CONSTRAINED QUADRATIC PROBLEMS 3 optimality and complementarity. Note that as Slater’s constraint qualifications are satisfied for the problem (2.2), all local solutions must satisfy (2.3). First observe that if there is no local solution at the boundary of the ball, then from compactness of the feasible region there must be a local solution x̄ in the interior (namely at least the global one), and therefore the unrestricted secondorder optimality conditions familiar from elementary calculus imply that A must be positive semidefinite. Clearly, Ax̄ = −b but there there is no alternative solution to Ax = −b (otherwise we would obtain a global solution also at the boundary). Hence A cannot be singular, and therefore A must be positive-definite. Consequently, the optimal solution reads x⋆  = −A−1 b for all   b ⊤ A−2 b. For smaller , we again must have (local) solutions on the boundary of the ball. So without loss of generality we may and do investigate in the sequel only solutions that satisfy x̄⊤ x̄ = , regardless of the definiteness of A. Returning to the general case, the second-order condition necessary for local optimality of a KKT point x̄ can be expressed as a lower bound for  in terms of the spectrum of A. Denote by 1 = 2 = ··· = k < k+1
···
n the eigenvalues of A, where the smallest eigenvalue of A has multiplicity k (possibly k = 1). Observe that the case where all eigenvalues coincide, i.e., where A = 1 In , is trivial: unless b = 0 (where x⋆ = o is the global solution for any  > 0 if 1  0 while x⋆
is some x with x⊤ x =  if 1 < 0), we get x⋆ = −b/1 + where  = −1 + b ⊤ b/ if this quantity is positive, or else x⋆ = −b/1 in the interior of the ball. So again without loss of generality let us assume in the sequel that A has a nonzero spectral spread. The next three results were established by Martínez [14]. For the readers’ convenience, we provide a self-contained (and slightly simplified) presentation. In the sequel, we sometimes use coordinates relative to the orthogonal eigensystem u1 un of A, where ui satisfy u⊤i ui = 1 and Aui = i ui . Then A is replaced with U ⊤ AU = diag1  n , and b with c = U ⊤ b, while the constraint x⊤ x = vU ⊤ Uv = v⊤ v
 remains unchanged. THEOREM 1 (Lemma 2.4 in [14]) If x̄ is a local solution to (2.2), then x̄ satisfies (2.3), i.e., x̄ is a KKT point of (2.2), and the Lagrange multiplier  must necessarily satisfy   −2  (2.4) Proof. We adapt a similar proof for global optimality conditions in [11]. Put ! > 0 (to be chosen later) and let d be some direction with −! < d ⊤ x̄ < 0. Then for t = −2d ⊤ x̄/d ⊤ d, we get x̄ +td⊤ x̄ +td = x̄⊤ x̄ −4 d ⊤ x̄ ⊤ d ⊤ x̄2 ⊤ x̄ +4 d d= d d⊤ d d ⊤ d2 4 PASQUALE L. DE ANGELIS, IMMANUEL M. BOMZE AND GERARDO TORALDO so that x̄ +td is feasible to (2.2). Now if $ > 0 is so small that f x  f x̄ for every feasible x with x − x̄ < $, choose ! = $ d . Then, by construction of t, 3 we get for x = x̄ +td that x − x̄ < $ and thus f x  f x̄. But 2 f x̄ +td−f x̄ = td ⊤ Ax̄ +b+ t2 d ⊤ Ad 2 = td ⊤ −x̄+ t2 d ⊤ Ad   2 = t2 d ⊤ A+In d + 2 −2td ⊤ x̄ −t 2 d ⊤ d   2 = t2 d ⊤ A+In d + 2 x̄2 −x̄ +td2 2 = t2 d ⊤ A+In d (2.5) Thus we arrive at the condition d ⊤ A+In d  0 for all directions d with −! < d ⊤ x̄ < 0. Now let v be a direction perpendicular to x̄. Then v can be arbitrarily well approximated by such a direction d, and continuity of the quadratic form finally yields v⊤ A+In v  0 for all v ⊥ x̄ (2.6) Evidently, this cannot happen if A+In has more than one negative eigenvalue (including multiplicities). Hence the result.  Theorem 1 can be also derived from much more general second-order optimality conditions for nonlinear optimization problems, see, e.g., [6] (and observe that Abadie’s constraint qualification for M& =
x ∈ IRn  x =  if & =  > 0 are satisfied in the present context). However, for the readers’ convenience we here provided a concise direct proof. The global case is also treated in Theorem 2.1 of [11] (earlier references are Lemmas 2.2 and 2.3 in [14] or, e.g., [17]): THEOREM 2 Suppose that x̄ satisfies (2.3), i.e., x̄ is a KKT point of (2.2). Then x̄ is a global solution to (2.2) if and only if the Lagrange multiplier  satisfies   −1  (2.7) Proof. The assertion is equivalent to the requirement that A+In be positive semi-definite, as formulated in [11].  The next step deals with a non-generic case which complicates analysis, and shows that then no non-global solution may exist: LEMMA 3 (Lemma 3.2 in [14]). If b ∈
u1  uk solution to (2.2) can exist. ⊥ , then no non-global local ELLIPSOIDAL APPROACH TO BOX-CONSTRAINED QUADRATIC PROBLEMS 5 Proof. Let  be the Lagrange multiplier of a non-global local solution x̄ to (2.2). Then, by Theorems 1 and 2, −2
 < −1 . Thus necessarily 1 < 2 so that k = 1. Now by hypothesis, 0 = −b ⊤ u1 = u⊤1 A+In x̄ = 1 +u⊤1 x̄ and  = −1 imply x̄ ⊥ u1 . On the other hand, u⊤1 A+In u1 = 1 + < 0, contradicting (2.6). Hence the result.  As a consequence of the previous proof, existence of a local, non-global solution implies simplicity of the smallest eigenvalue 1 of A (i.e., k = 1). Next we attack the question how many local solutions to (2.2) may coexist: to this end the so-called secular function [19] 't = n  ci2 2 i=1 i +t plays an important role. This function is well-defined on D' = IR+ \
−1  −n , and ' is strictly convex on D' . Observe that if the multiplier  ∈ D' , then  determines the KKT point uniquely since A+In is then non-singular and x̄ = −A+In −1 b (2.8) as can be seen from the conditions (2.3)(i). But the remaining (ii) and (iii) also imply that for any  ∈ D' \
0 , we get x̄ = U v̄ with v̄i = − ci i +  where  satisfies ' = ni=1 v̄i 2 = . This establishes the significance of the secular function '. We next show that under conditions which are generically (w.r.t. A and b) satisfied, there are at most three local solutions to (2.2). THEOREM 4 Suppose that b 
u1  uk ⊥ ∪
uk+1  um ⊥ where
uk+1  um is the eigensystem to the eigenvalue k+1 with multiplicity m− k+1  1 (if the smallest two eigenvalues of A are simple, this condition reduces to b ⊤ u1 b ⊤ u2  = 0. Then there are at most three local solutions to (2.2) at the boundary of the ball. Every such local solution x̄ satisfies x̄ = U v̄ with v̄i = −ci /i +, where  > −k+1 satisfies  = −1 and ' = . Proof. In coordinates w.r.t. the eigensystem of A, the assumption b 
u1  uk ⊥ reads ci = b ⊤ ui = 0 for some i ∈
1  k . Without loss of generality we assume c1 = 0. But as diag 1 +  n +v̄ = U ⊤ A+In UU ⊤ x̄ = U ⊤ −b = −c 6 PASQUALE L. DE ANGELIS, IMMANUEL M. BOMZE AND GERARDO TORALDO we get 1 +v1 = −c1 = 0, so that  = −1 . Similarly we obtain from b 
uk+1  um ⊥ (again, after rearranging indices in case of m−k+1 > 1 if necessary) from ck+1 = 0 also  = −k+1 . Hence we obtain via Theorem 1 that  ∈ D' ∩ −k+1 + . On this region, the secular function ' has only one pole located at −1 . For t → 1 , we get 't → +, and also 't → + as t ց −k+1 while 't ց 0 as t ր +. From strict convexity of ' on this region, we infer that there are at most three distinct solutions 1 2 3 to the equation ' =  (at least one for every , two for a particular value of , and three for all larger ). To every such i there corresponds exactly one KKT point, and so there are at most three local solutions.  Note that similar but refined arguments are used by Lucidi et al. in [13] to bound the number of KKT points of (2.2). Also, the last result can be sharpened considerably for global solutions: THEOREM 5 Suppose that b 
u1  uk ⊥ (if the smallest eigenvalue of A is simple, this reduces to b ⊤ u1 = 0). Then there is exactly one global solution to (2.2) at the boundary of the ball. This solution x̄ satisfies x̄ = U v̄ with v̄i = −ci /i +, where  > −1 satisfies ' = . Proof. The arguments are virtually the same as in the previous theorem. Indeed, using Theorem 2, we see that now the relevant part of D' shrinks to the region IR+ ∩ −1 + , upon which the secular function is even smooth.  Now we prove the sharpening of Theorem 4 due to Martínez: THEOREM 6 (Lemma 3.2 in [14]) There is at most one local, non-global solution x̄ to (2.2), and this solution has a Lagrange multiplier  which satisfies '′   0 where '′ is the derivative of the secular function. Proof. We only have to establish the asserted inequality. The remainder follows from the arguments in Theorem 4 (use Lemma 3). Starting point of our reasoning is the relation (2.6). We first construct a basis of x̄⊥ . Consider the vectors wj = cj c1 u u1 − j + 1 + j j = 2  n Since c1 = 0 as shown in the proof of Theorem 4, these vectors
w2  wn are linearly independent and they all satisfy wj ⊥ x̄, because of x̄ = U v̄ = − ni=1 ci ui /i +, so that   n  cj ci c1 ⊤ ⊤ ⊤ u u = 0 wj x̄ = u1 ui − 1 + j i i=1 i + j + ELLIPSOIDAL APPROACH TO BOX-CONSTRAINED QUADRATIC PROBLEMS 7 Thus, forming the n×n−1 matrix B = w2  wn , we see that condition (2.6) can be rephrased as a positive-semidefiniteness condition for the symmetric n− 1×n−1 matrix B ⊤ A+In B. Now observe that A+In wj = cj 1 + c1 j + u1 − u j + 1 + j whence we deduce that the entry in the ith column and jth row of B ⊤ A+In B equals wi⊤ A+In wj = ci cj 1 + c 2  +2 + 1 i  i +j + 1 +2 ij where ij is the Kronecker symbol. Written in compact form, we thus obtain B ⊤ A+In B = D −ss ⊤ where D is a diagonal, positive-definite matrix, and s has coordinates
sj = 1 + cj j + for j = 2  n. Hence B ⊤ A+In B is positive-semidefinite if and only if In−1 −D−1/2 sD−1/2 s⊤ is so, with  1 + D−1/2 = diag
c1  j + j=2  n designating the inverse of the square-root factorization of D. But a symmetric rankone update of In−1 in this form, In−1 −vv⊤ is positive-semidefinite if and only if v⊤ v
1, as a direct argument via v⊥ and the Cauchy-Schwarz-Bunyakovsky inequality shows. Finally, we calculate 1  D−1/2 s⊤ D−1/2 s = s ⊤ D−1 s = − n cj2 1 +3  3 c12 j=2 j + which obviously is equivalent to '′   0.  Finally, we obtain a simple genericity condition guaranteeing that there are at most two local solutions: COROLLARY 7 If there is a local, non-global solution x̄ to (2.2), then there is no other local-nonglobal solution. If, furthermore, b ⊤ u1 = 0, then there are in total two local solutions, namely x̄ and the (unique) global solution x∗ . 8 PASQUALE L. DE ANGELIS, IMMANUEL M. BOMZE AND GERARDO TORALDO Proof. Existence of x̄ implies, as observed in the proof of Lemma 3, k = 1. Hence the result follows by Theorems 6 and 5.  At the end of this section, we investigate continuous dependence of the global solution from the parameter  (the square of the ball radius). For completeness, we also deal with the non-generic case where b ∈
u1  uk ⊥ , and provide an explicit construction both for the multiplier ⋆  and the optimal solution x⋆ . THEOREM 8 For varying , we can choose a global solution x⋆  to (2.2), which depends continuously upon . This dependence is smooth if b 
u1  uk ⊥ , whereas it has exactly one point of non-differentiability in case of b ∈
u1  uk ⊥ . Also the corresponding multiplier ⋆  exhibits the same smoothness and continuity behavior. Proof. The first assertion can be derived from Theorem 5, invoking the Implicit Function Theorem, applied to the equation F t  = 't− on the region −1 + × 0 + where F is well-defined and smooth. However, if b ∈
u1  uk ⊥ , we have to argue more carefully since there may be more than one global solution, and since the case  = −1 (the largest potential pole of ') is not ruled out a priori. However, in this case we get 't = ci2 2 i>k i +t  which is well defined for all t > −k+1 (in particular, for t = −1 ). Now put crit = '−1 . If  < crit , then there is a unique ⋆  =  > −1 such that ' = , since ' decreases strictly on −1 + . In this case, put vi⋆  = 0 if i
k whereas vi⋆  = − ci if i > k i + Then x⋆  = Uv⋆  is the only KKT point with multiplier ⋆   −1 , and thus the only global solution. Again invoking the Implicit Function Theorem, we obtain smooth dependence of the above quantities upon . If, however,   crit , then define ⋆  = −1 , a constant in , and put
ci if i > k vi⋆  = −'−1  /k if i
k whereas vi⋆  = − i −1 Then the continuity assertion is a consequence of ⋆  ց −1 as  ր crit , which in turn follows from strict monotonicity and continuity of ' on −1 + . Finally, we obtain via '′ −1  < 0 and the Implicit Function Theorem that the left-hand side derivative of ⋆  at crit satisfies .F −1 crit  −1 d⋆ . crit  = − .F <0 =− ′ d ' −1  −1 crit  .t ELLIPSOIDAL APPROACH TO BOX-CONSTRAINED QUADRATIC PROBLEMS so that there must be a kink at crit for ⋆ . 9  Here we should stress that it is not suggested to employ the above explicit construction in numerical application. It serves only for establishing the theoretical result. Only in the highly unlikely case of b ∈
u1  uk ⊥ one might be forced to return to the above (or similar) constructions. Observe that for testing the latter conditions, we need only information on the eigenspace for the smallest eigenvalue rather than complete spectral decomposition of A. Finally, we obtain an important monotonicity result for the multiplier: COROLLARY 9 The Lagrange multiplier ⋆  is (can be chosen as) decreasing in . If b ∈
u1  uk ⊥ , then ⋆  strictly decreases on the interval 0 crit (and is constant for larger ), while if b 
u1  uk ⊥ , then ⋆  decreases strictly on the whole IR+ . 3. The algorithm Once we have established continuous dependence of the global solution x⋆  upon the radius, we can proceed as follows: 3.1. algorithm ebh 1. Start with 0 = n (then the ball is containing the box B), determine 0 = ⋆ 0  and x0 = x⋆ 0  by solving the corresponding trust region problem. 2. If x0 ∈ B, stop and return x0 that is the exact solution of the BCQP; else determine the multiplier  = ⋆ 1 for the ball completely contained in B and put k = 1. 3. Calculate k = 0 + /2 xk = −A+2k In −1 b (then for k = xk⊤ xk we have xk = x⋆ k  and k = ⋆ k , according to (2.8) and Theorem 5). 4. If xk ∈ B and near the boundary .B of B, stop and return xk as the promising point (to be used as starting point for any local solver); else, if xk ∈ B then put  = k else put 0 = k . 5. Replace k with k+1 and repeat from 3. From Corollary 9 and Theorem 8 we can be sure that the above procedure is finite. The next step is a branch-and-bound approach with the help of the obtained promising point xk . We cut the box B in two halves by bisecting along that coordinate direction i which is most promising in the sense that the minimal deviation of the i-th coordinate of xk from ±1 is largest. Repeating this recursively, we generate 10 PASQUALE L. DE ANGELIS, IMMANUEL M. BOMZE AND GERARDO TORALDO 2m subboxes by m such cuts. Of course, subboxes for which the minimum of f x on the including ball is larger than the best current solution can be discarded from further analysis. This cannot be used as an exact algorithm, because of the prohibitive computational cost that it would require, so a heuristic stopping criterion must be used. In our computational experiments we used as a stopping criterion a prefixed number maxit of cuts. 4. Computational results A prototype Matlab version of algorithm EBH was tested on a set of test problems. For the minimization of f x over a ball, we used the code trustpen [13]. As a local solver for the box-constrained quadratic problem solved we used the Matlab routine quadprog [15]. As already pointed out, for a prefixed maximum number maxit of cuts, at most 2maxit problems of the form (2.2) had to be solved. 4.1. TEST PROBLEMS WITH KNOWN GLOBAL SOLUTION We first try test problems with known global, generated as in [10], moving from the sufficient global optimality condition originally proposed by Neumaier [16]. However, we found such problems very easy, in the sense that the Matlab local solver was almost always able to find the global solution. Therefore for constructing test problems, we used the sufficient global optimality conditions recently described by Beck and Teboulle [2] for quadratic optimization problems with binary constraints, that can be easily rephrased for concave box-constrained quadratic problems. THEOREM 10 Consider the problem (1.1), with A negative definite. Then x is a global solution for (1.1) if x is feasible and &max Ae  XAXe +Xb where e = 1  1 ⊤ ∈ IRn , eigenvalue of the matrix A. (4.9) X = diagxi , and ≪max A is the largest Based on Theorem 10, concave test problems were generated with A = HDH , with D a diagonal matrix whose elements are randomly generated from the interval −2 −1 ; H a Householder matrix of a random vector; the prefixed solution x randomly generated with components in
−1 +1 , and the vector b is computed so that the condition (4.9) holds. Test problems were generated with n = 50 100 200; for each dimension 10 different problems were generated. The value of maxit was set equal to 15 in our experiments. From the computational results (Table I), we see that for such problems, the ellipsoidal approach show to be very effective, since only for one instance of such problems the global solution was not found; to be fair, we note that problems ELLIPSOIDAL APPROACH TO BOX-CONSTRAINED QUADRATIC PROBLEMS 11 Table I. Failure rates for EBH n Failure 50 100 200 0 1 0 with known solution are usually not very difficult, because of the sufficient global optimality conditions they are based on. So are the problems we generated, for which the linear term is very strong. Nevertheless, to check that the problems are not trivial, for each of them we run the linear solver 10n times, with randomly generated starting points; we found that for none of them, the stationary point is unique. For example, for the problems of size 100, we found for each of them 25 stationary points on the average. Also, for the local solvers, the problems shown to be harder than the ones we generated as in [10]. 4.2. TEST PROBLEMS WITH UNKNOWN GLOBAL SOLUTION Next, we generated a first set of random test problems with unknown solution. The matrix A has the form A = HDH , where H is a Householder matrix of a random vector and D = diagd, with d a random vector such that di  ∈ 1 2 . For each matrix we generate five different problems, by taking b = d/10p , with p = 0 1 2 3 4 5. For p = 5 we have problems in which the quadratic part of the objective function dominates, and this is the opposite which happens for p = 0. Several problems were generated varying the size n and the percentage Neg of negative eigenvalues of the Hessian matrix. To check the quality of the solution that was found with the ellipsoidal approach, for each problem we run n times the two local solvers with randomly generated starting points, and then we compared the best solution (random solution), to the one we found using the ellipsoidal approach. As expected, the number of stationary points for the test problems with unknown global solution is much higher than for the test problems with known global solution. For the problems with 100 variables, in finding the random solution, 25 different stationary points were found on the average. From Table II, it is observed that the ellipsoid algorithm was able to find the best (known) solution in 24 out of the 30 test problems; actually, for 13 problems the solution that the algorithm was able to detect, was better than the one which was found by the random approach. We also point out that for the six problems in which the algorithm failed, just doubling the maximum number of iterations maxit, the EBH algorithm was eventually able to get the best solution in four cases. The number of failures (i.e., the problems in which the random solution was better than the one found by the algorithm) is very low; although there is of course no guarantee that the solution the algorithm was able to find is the global one, even the suboptimal solutions found were of considerable high quality (measured by 12 PASQUALE L. DE ANGELIS, IMMANUEL M. BOMZE AND GERARDO TORALDO Table II. Performance of EBH on random test problems (See text) n Neg Failures 50 50 50 100 100 100 25 50 75 25 50 75 1 1 1 0 1 2 the fraction of EBH-result over the maximum over all random solutions found). It seems that the effectiveness of the algorithm does not deteriorate when the number of negative eigenvalues increases (since the number of failures increases just slightly). On the other hand, for Neg = 75 and Neg = 50, the ellipsoidal approach beats the random approach in 11 problems out of 20; in the problems with Neg = 25, the two approaches almost always get the same solution. 4.3. MAXIMUM CLIQUE PROBLEMS The maximum clique problem is a classical problem in graph theory which can be formulated as follows. Let G = GV E, with vertex set V =
1 2  n and edge set E ⊆ V ×V , be an undirected graph. The maximum clique problem is the problem of finding a complete subgraph of G of maximal cardinality. This problem has many equivalent formulations as an integer problem, or a continuous nonconvex optimization problem, for a recent survey see [4]. Here we follow a recent continuous reformulation approach from [1]: starting from the following quadratic zero-one reformulation min f x = x⊤ Ax −e⊤ x  s.t. x ∈
0 1 n (4.10) where A is the adjacency matrix of the complement graph of G and e is the vector of all ones, consider the continuous relaxation problem min f x = x⊤ Ax −e⊤ x  s.t. x∈ 0 1 n  (4.11) Since the matrix A has zeroes on the diagonal, this problem has at least one 0–1 solution, then the global solution of (4.10) is equal to (one of the) global solution(s) of (4.11), and therefore a maximum clique can be computed by solving the problem (4.11). The EBH algorithm was then tested on problem (4.11) generated by a set of some benchmark graph instances from the DIMACS challenge [9]. We selected only graphs with less than 400 vertices. The results are shown in Table 3. On purpose, we did not include the results of the similar reformulation approach in [1], because not all instances in Table III were considered there. In short, the authors 13 ELLIPSOIDAL APPROACH TO BOX-CONSTRAINED QUADRATIC PROBLEMS Table III. Maximum clique problems from the DIMACS challenge (see text) Problem c-fat200-1 c-fat200-2 c-fat200-5 hamming6-2 hamming6-4 hamming8-2 hamming8-4 johnson8-2-4 johnson8-4-4 johnson16-2-4 p_hat300-1 p_hat300-2 p_hat300-3 MANN_a9 MANN_a27 keller4 brock200_1 brock200_2 brock200_3 brock200_4 san200_0.7_1 san200_0.7_2 san200_0.9_1 san200_0.9_2 san200_0.9_3 sanr200_0.7 sanr200_0.9 n Dens. (%) Clique size EBH (cuts) CBH SR AR 200 200 200 64 64 256 256 28 70 120 300 300 300 45 378 171 200 200 200 200 200 200 200 200 200 200 200 07.7 16.3 42.6 90.5 34.9 96.9 63.9 55.6 76.8 76.5 24.4 48.9 74.5 92.7 99.0 64.9 74.5 49.6 60.5 65.8 70.0 70.0 90.0 90.0 90.0 69.7 89.8 12 24 58 32 4 128 16 4 14 8 8 25 36 16 126 11 21 12 15 17 30 18 70 60 44 18 > 42 − 12 (1) 24 (1) 58 (1) 32 (1) 4 (1) 128 (1) 16 (1) 4 (1) 14(1) 8 (1) 8 (2) 24 (3) 32 (1) 16 (1) 125 (4) 9 (8) 20 (1) 10 (1) 11 (1) 15 (5) 15 (3) 12 (1) 46 (2) 40 (10) 34 (9) 16 (7) 35 (3) 12 24 58 32 4 128 16 4 14 8 8 25 36 16 121 10 20 12 14 16 15 12 46 36 30 18 41 12 24 58 32 4 128 16 4 14 8 6 22 32 12 117 7 17 8 9 12 15 12 45 36 32 14 37 12 24 58 32 4 128 16 4 14 8 8 25 35 16 117 8 19 10 13 14 15 12 45 39 31 16 41 report in [1] the same figures as in Table III for the first 10 lines. Their algorithm has equal performance as EBH on the Mannino graphs MANN_ax (x ∈
9 27 ) and the Sanchis graphs san200_0.9_x (x ∈
2 3 ), beats EBH by one node for san200_0.9_1, and hits the maximum clique in the Keller graph keller4. All other instances were unfortunately not covered in [1]. The columns labelled by CBH, SR, and AR, refer to different but related continuous-based heuristics for the maximum clique problem. Essentially, all these are based upon the Motzkin– Strauß program which reformulates the Maximum Clique Problem as a Standard Quadratic Optimization Problem, i.e. to find the maximum of a quadratic form over the standard simplex (rather than over a box). Here, the matrix generating the quadratic form is the adjacency matrix of the graph. The first method considered is the Continuous Based Heuristic (CBH) [8], which employs a parameterized version of the original Motzkin–Strauß program. The problem is divided into a series of subproblems with the simplex constraints relaxed into spherical ones (this is a related feature to our approach). A combinatorial post-processing phase is needed to round the solutions produced by the procedure that solves the subproblems back to the standard simplex. 14 PASQUALE L. DE ANGELIS, IMMANUEL M. BOMZE AND GERARDO TORALDO The second procedure is the Simple Replicator Dynamics approach (SR) [5] which basically finds local solution to a regularized version of the Motzkin-Strauß program (adjacency matrix plus 1/2 times the identity matrix) with the help of replicator dynamics, a discrete time (non-Euclidean gradient) dynamical system operating on the standard simplex (thus feasibility is guaranteed throughout the algorithm). The third algorithm is Annealed Replication (AR) [3]. It uses a parameterized quadratic form emerging from perturbing similar as in SR the adjacency matrix (here plus 8 times the identity matrix), and iteratively searches for local solutions via replicator dynamics. In order to escape inefficient local solutions, thereafter the parameter 8 is varied as in simulated annealing, but motivated by more principled arguments. 5. Conclusions Moving from some new theoretical results concerning the global minimization of a quadratic function over a ball, and more specifically the continuous dependence of the Lagrange multiplier for the global solution on the radius of the ball, we proposed an heuristic algorithm for the global solution of box constrained quadratic programming problems. The algorithm, embedded in a branching strategy, showed to be successful in many of the test problems we used in our computational experiences. In particular, our approach gave reasonable good results on a wide set of some benchmark graph instances from the DIMACS challenge, although no special care was taken to fit the algorithm to the maximum clique problem. We therefore believe that the proposed algorithm can be fruitfully used, at least, as an excellent procedure to give a very good starting point to any local solver. Acknowledgement The authors want to thank Laura Palagi for making the code trustpen available to them. G. T. was supported by the MURST National Project ‘Algorithms for Complex Systems Optimization’. References 1. J. Abello, S. Butenko, P.M. Pardalos and M.G.C. Resende (2001), Finding independent sets in a graph using continuous multivariable polynomial formulations. J. Global Optimization. 21, 111–137. 2. A. Beck and M. Teboulle (2000), Global optimality conditions for quadratic optimization problems with binary constraints. SIAM J. Optim. 11, 179–188. 3. I.M. Bomze, M. Budinich, M. Pelillo and C. Rossi (2001), Annealed replication: A new heuristic for the maximum clique problem, to appear in: Discrete Appl. Math. 4. I.M. Bomze, M. Budinich, P.M. Pardalos and M. Pelillo (1999), The maximum clique problem. In: D.Z. Du and P.M. Pardalos (Eds), Handbook of Combinatorial Optimization, suppl. Vol. A, 1–74. Kluwer Academic Publishers, Dordrecht. ELLIPSOIDAL APPROACH TO BOX-CONSTRAINED QUADRATIC PROBLEMS 15 5. I.M. Bomze, M. Pelillo and R. Giacomini (1997), Evolutionary approach to the maximum clique problem: Empirical evidence on a larger scale. In: I. M. Bomze, T. Csendes, R. Horst and P. M. Pardalos (Eds.), Developments in Global Optimization, pp. 95–108. Kluwer Academic Publishers, Dordrecht. 6. G. Danninger and I.M. Bomze (1994), Generalized convexity for second-order optimality conditions. In: S. Komlósi, T. Rapcsák, S. Schaible (Eds.), Generalized Convexity, pp. 137–144. Springer, Berlin. 7. P.L. De Angelis, P.M. Pardalos and G. Toraldo (1997), Quadratic Programming with box constraints. In: I.M.Bomze, T. Csendes, R. Horst and P.M. Pardalos (Eds.), Developments in Global Optimization, pp. 73–93. Kluwer Academic Publishers, Dordrecht/Boston/London. 8. L. E. Gibbons, D. W. Hearn and P. M. Pardalos (1996), A continuous based heuristic for the maximum clique problem. In: D. S. Johnson and M. Trick (Eds.), Cliques, Coloring, and Satisfiability—Second DIMACS Implementation Challenge, pp. 103–124. American Mathematical Society, Providence, RI. 9. J. Hasselberg, P.M. Pardalos and G. Vairaktarakis (1993), Test case generators and computational results for the maximum clique problem, Journal of Global Optimization 3, 463–482. 10. C.G. Han, P.M. Pardalos and Y. Ye (1992), On the solution of indefinite quadratic problems using an interior point algorithm, Informatica 3, 474–496. 11. J.B. Hiriart-Urruty (1995), Conditions for global optimality. In: R. Horst and P. M. Pardalos (Eds.), Handbook of Global Optimization, pp. 1–26. Kluwer Academic Publishers, Dordrecht. 12. M.K. Kozlov, S.P. Tarasov and L.G. Khacian (1979), Polynomial solvability of convex quadratic programming, Soviet. Math. Dokl. 20, 1108–1111. 13. S. Lucidi, L. Palagi and M. Roma (1998), On some properties of quadratic programs with a convex quadratic constraint, SIAM J. Optim. 8, 105–122. 14. J.M. Martínez (1994), Local minimizers of quadratic functions on Euclidean balls and spheres, SIAM J. Optim. 4, 159–176. 15. Matlab - Optimization Tolbox 2.1 (2001), The MathWorks, Inc. 16. A. Neumaier (1992), An optimality criterion for global quadratic optimization, J. Global Optimization 2, 201–208. 17. J. Moré and D.C. Sorensen (1983), Computing a trust region step, SIAM J. Sci. Stat. Comput. 4, 553–572. 18. P.M. Pardalos and G.P. Rodgers (1990), Computational aspects of branch and bound algorithms for quadratic zero-one programming, Computing 45, 131-144. 19. R.J. Stern and H. Wolkowicz (1994), Trust region problems and nonsymmetric eigenvalue perturbations, SIAM J. Matrix Anal. Appl. 15, 755–778. 20. Y. Yajima, M.V. Ramana and P.M. Pardalos (2001), Cuts and Semidefinite Relaxations for Nonconvex Quadratic Problems, In: N. Hadjisavvas, J.E. Martinez-Legaz and J.P. Penot (Eds.), ‘Generalized Convexity and Generalized Monotonicity,’ Lecture Notes in Economics and Mathematical Systems, Vol. 502, pp. 48–70. Springer, Berlin.