A clique algorithm for standard quadratic programming

Discrete Applied Mathematics 156 (2008) 2439–2448 www.elsevier.com/locate/dam A clique algorithm for standard quadratic programming Andrea Scozzari, Fabio Tardella ∗ Università di Roma “La Sapienza”, Dip. di Matematica per le Decisioni Economiche, Finanziarie ed Assicurative, Roma, Italy Received 7 September 2007; accepted 20 September 2007 Available online 20 February 2008 Dedicated to the memory of Peter L. Hammer Abstract A standard Quadratic Programming problem (StQP) consists in minimizing a (nonconvex) quadratic form over the standard simplex. For solving a StQP we present an exact and a heuristic algorithm, that are based on new theoretical results for quadratic and convex optimization problems. With these results a StQP is reduced to a constrained nonlinear minimum weight clique problem in an associated graph. Such a clique problem, which does not seem to have been studied before, is then solved with an exact and a heuristic algorithm. Some computational experience shows that our algorithms are able to solve StQP problems of at least one order of magnitude larger than those reported in the literature. c 2007 Elsevier B.V. All rights reserved. Keywords: Standard quadratic optimization; Indefinite quadratic programming; Maximum clique; Nonlinear programming 1. Introduction and notation A standard quadratic optimization problem (StQP) consists of finding (global) minimizers of a quadratic form over the standard simplex. It can be written in the form f ∗ = min{ f (x) = x T C x : x ∈ ∆}, (1) where C is a symmetric matrix; a T denotes transposition; and, with a slight abuse of notation, ∆ denotes the standard simplex in a Euclidean space Rn of any dimension: ∆ = {x ∈ Rn : eT x = 1, x ≥ 0}, where e = [1, . . . , 1]T . The standard quadratic optimization problem has been singled out by Bomze [3], who introduced the name and described several properties and applications of this class of problems. Several other authors have recently worked on this topic presenting exact or heuristic solution algorithm, applications, and connections with other problems (see [2–5,7,8,10,15,19] and the references therein). ∗ Corresponding author. Tel.: +39 064 976 6275; fax: +39 064 9766 765. E-mail addresses: andrea.scozzari@uniroma1.it (A. Scozzari), fabio.tardella@uniroma1.it (F. Tardella). 0166-218X/$ - see front matter c 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.dam.2007.09.020 2440 A. Scozzari, F. Tardella / Discrete Applied Mathematics 156 (2008) 2439–2448 Among other applications, the StQP has been used as a tool to solve maximum weight clique problems on the basis of a well-known reduction described by Motzkin and Straus [14] for the unweighted case, and by Gibbson et al. [12] for the weighted case. Of course these reductions also imply NP-completeness of StQP in the general case. In this work we establish a converse reduction from a StQP to a nonlinear constrained minimum weight clique problem. We then use this reduction to provide an exact and a heuristic algorithm for solving a StQP by solving the associated clique problem. To describe the reduction, first note that the vertices of the standard simplex are the unit vectors ei , and its graph is the complete graph on n vertices. We associate with the objective function f a convexity graph G f = (N , E), where N = {1, . . . , n}, and an edge (i, j) is in E if and only if f is strictly convex on the edge of ∆ joining the vertices ei and e j , i.e., if cii + c j j − 2ci j > 0. The intersection of a polyhedron P with one of its supporting hyperplanes is a face of P. The standard simplex ∆ has 2n faces, which are also standard simplices in lower dimensions. To every subset I of N we associate the face ∆ I = {x ∈ Rn : (e I )T x I = 1, x ≥ 0}, where, for a vector z ∈ Rn , z I denotes the subvector formed by its entries with indices in I . Clearly, the problem of minimizing f on ∆ I coincides with the problem of minimizing the function f I (x) = x T C I x on ∆, where C I is the submatrix of C formed by its entries with row and column indices in I . Let ϕ(I ) = min f I (x). x∈∆ Then ϕ is antitone on 2 N , i.e., ϕ(I ) ≥ ϕ(J ) whenever I ⊆ J . Hence, f ∗ = ϕ(N ) = min ϕ(I ). I ⊆N In the following section we show that the global minimum f ∗ is also equal to the minimum of ϕ(I ) on some subsets of the family of all (maximal) cliques of the convexity graph G f . Furthermore, we prove that the global minimizer belongs to the last face of an increasing sequence of faces of ∆ where f is strictly convex and has interior global minimizers. In Section 3 we devise an exact and a heuristic algorithm for a constrained nonlinear minimum weight clique problem that can be applied to the type of problems arising from the reduction of StQPs described in Section 2. The computational results described in Section 4 show that the size of the StQPs that can be solved with our algorithms is at least one order of magnitude larger than the ones solved in the previous works. We conclude with some final remarks and perspectives for future work. 2. Fundamental results We present here two results that are at the basis of the validity of the exact and heuristic clique algorithms that will be described in Section 3. The first result, which is a special case of more general results described by Tardella in [22, 23], is a Quadratic Programming version of the fundamental theorem of Linear Programming. We provide a proof here for the sake of completeness. Theorem 1. A quadratic function f that is bounded below on a pointed polyhedron P attains its minimum on P in the relative interior of a face where f is strictly convex. Proof. We first show that for any point x in the relative interior of a face F where f is not strictly convex, there exists a point y in a lower-dimensional face such that f (y) ≤ f (x). Indeed, if f is not strictly convex on F, there exists a line through x where f is concave. Since P is pointed, the intersection of this line with F is either a segment with endpoints y 1 , y 2 , or a half-line with origin y 3 . Clearly, y 1 , y 2 , y 3 belong to faces of P of dimension smaller than F. Furthermore, concavity implies that min{ f (y 1 ), f (y 2 )} ≤ f (x), while concavity and boundedness from below imply that f (y 3 ) ≤ f (x). Hence, in either case, we have a point y with f (y) ≤ f (x). Since every function is trivially strictly convex on zero-dimensional faces (vertices), by iterating the above argument we can conclude that for any point x in the relative interior of a face F where f is not strictly convex, there exists a point y in a lower-dimensional face where f is strictly convex such that f (y) ≤ f (x). A. Scozzari, F. Tardella / Discrete Applied Mathematics 156 (2008) 2439–2448 2441 Let Fc denote the set of all faces of P where f is strictly convex. Then from the previous argument we deduce that inf f (x) = min inf f (x). x∈P F∈Fc x∈F Since strictly convex quadratic functions are coercive (i.e., limkxk→+∞ f (x) = +∞), they attain their minimum on all closed sets. Hence we can conclude that the global minimum of f on P is attained in a face where f is strictly convex.  Note that when f is linear, the only faces of P where f is strictly convex are the vertices of P. So the Fundamental Theorem of LP is a special case of Theorem 1. Furthermore, the existence part of Theorem 1 is a well-known result of Frank and Wolfe [11]. We now generalize the construction described in Section 1 for the standard simplex. Let P be a polytope and N be the set of its vertices. We associate to a quadratic function f on P a graph G f = (N , E), called convexity graph of f , by setting (i, j) ∈ E if and only if either the line segment joining vertices i and j in P is not an edge of P, or if it is an edge of P and f is strictly convex on it. From Theorem 1 we then easily derive the following result: Corollary 2. A quadratic function f attains its minimum on a polytope P in the relative interior of a face whose vertices induce a clique in the associated convexity graph. Proof. The proof follows from Theorem 1 and from the definition of the convexity graph by observing that a strictly convex function on a face F is also strictly convex on all the edges of F.  Let us now focus on standard quadratic optimization problems, where P = ∆ and f (x) = x T C x. Let rint(S) denote the relative interior of a set S in Rn , C denote the set of all cliques of the convexity graph associated to f on ∆, and let   1 C ≻ = I ∈ C : C I ≻ 0 and x̄ = C I−1 e ∈ rint(∆ I ) 2 be the set of all cliques corresponding to those faces of ∆ where f is strictly convex and has an interior minimizer. Furthermore, for any subset A of 2 N , let A M denote the subset formed by all its (inclusionwise) maximal elements. We can then restrict the search for a global minimizer of f on ∆ to the (maximal) faces of ∆ where f is strictly convex and has an interior global minimizer. Corollary 3. min f (x) = min ϕ(I ) = min≻ ϕ(I ). x∈∆ I ∈C ≻ I ∈C M Proof. Equality minx∈∆ f (x) = min I ∈C ≻ ϕ(I ) is a rephrasing of Corollary 2. Furthermore, we have min I ∈C ≻ ϕ(I ) = min I ∈C M≻ ϕ(I ) by the antitonicity of ϕ(I ).  Note that ϕ(I ) can be easily evaluated on all cliques I ∈ C ≻ . Indeed, if x̄ I is the global minimizer of f in the relative interior of ∆ I , then, by first-order optimality conditions, we must have x̄ I = 21 C I−1 e and ϕ(I ) = f (x̄ I ) = 1 ⊤ −1 e C I e. 4 (2) We now present our second basic result in a setting that is more general than the one required for the algorithmic developments in the next section. We recall that a function f is strictly quasi-convex on a convex set X if f (x) < max{ f (x 1 ), f (x 2 )} for every pair of points x 1 , x 2 ∈ X and every point x in the line segment joining x 1 and x 2 . Let F be a face of P where f is strictly convex and let x F∗ denote the (unique) global minimizer of f on tng F, the tangential hull of F, defined as the smallest affine subspace containing F. Clearly, x F∗ need not belong to F in general. If x F∗ belongs to rint(F), we say that f has an interior global minimizer on F. A face of P of dimension dim(P) − 1 is called a facet of P. 2442 A. Scozzari, F. Tardella / Discrete Applied Mathematics 156 (2008) 2439–2448 Theorem 4. Let f be a strictly quasi-convex differentiable function that attains its minimum on a polytope P at an interior point x ∗ of P. Then there exists a facet F of P such that the minimum of f on F is attained at a point in its relative interior. Proof. Assume that P = {x ∈ Rn : ai · x ≤ bi , i = 1, . . . , m}, where ai ∈ Rn , and bi ∈ R. We can also assume that P has m facets described by Fi = {x ∈ P : ai · x = bi }, and that a j 6= tak for all j 6= k and for every t ∈ R+ . Let Hi = tng Fi = {x ∈ Rn : ai · x = bi }, for i = 1, . . . , m, and let x ′ be a solution of the problem min min f (x). 1≤i≤m x∈Hi (3) Then we must have x ′ ∈ P. Indeed, if x ′ 6∈ P the line segment joining x ∗ and x ′ intersects one of the faces of P at a point x ′′ where f (x ′′ ) < max{ f (x ∗ ), f (x ′ )}. From strict quasi-convexity of f we also have that f (x ∗ ) < f (x ′ ), so that f (x ′′ ) < f (x ′ ) contradicting the definition of x ′ . Since x ′ ∈ P, we must have x ′ ∈ F j for some index j. We now prove that x ′ ∈ rint(F j ). Indeed, if this is not the case we must have x ′ ∈ F j ∩ Fk for some k 6= j, so that f (x ′ ) = minx∈H j f (x) = minx∈Hk f (x) by definition (3). From first-order optimality conditions we then get ∇ f (x ′ ) = λ j a j = λk ak for some λ j , λk ∈ R+ contradicting the assumption a j 6= tak for all t ∈ R+ .  With repeated applications of the previous theorem we easily obtain the existence of a nested sequence of faces with interior global minimizers. Corollary 5. Let f be a strictly quasi-convex differentiable function that attains its minimum on a polytope P at an interior point x ∗ of P. Then there exists a nested sequence of faces F 1 ⊂ F 2 ⊂ · · · ⊂ F n such that F n = P, dim(F i ) = i, and x F∗ i ∈ rint(F i ) for i = 1, . . . , n. We can now specialize the previous result to the quadratic case. Theorem 6. For a quadratic function f on a polyhedron P there exists a nested sequence of faces F 1 ⊂ F 2 ⊂ · · · ⊂ F k of P where f is strictly convex, and such that x F∗ k is a global minimizer of f on P, dim(F i ) = i and x F∗ i ∈ rint(F i ) for i = 1, . . . , k. Proof. By Theorem 1 there exists a global minimizer x ∗ of f on P in the relative interior of a face F of P where f  is strictly convex. The proof then follows by applying Corollary 5 to polytope F. 3. An algorithm for nonlinear constrained minimum weight clique problems In view of Corollary 3, the problem of minimizing the quadratic form f (x) = x T C x on the standard simplex ∆ can be transformed into a nonlinear constrained minimum weight clique problem. Unfortunately, we could not find in the literature any reference to this type of problems. So we developed the new algorithm described below, which is an extension of the approaches of Östergard [20,21] to the standard unconstrained maximum (weight) clique problem. We address the problem of minimizing a (nonlinear) function ϕ(I ), defined on all subsets I of the set of vertices N of a graph G = (N , E), on a subset D of the family C of all cliques of G. We make the following assumptions: 1. ϕ(I ) is antitone; 2. ϕ(I ) can be computed efficiently on D, and for all sets I at least a lower bound L B(ϕ(I )) on ϕ(I ) can also be obtained efficiently; 3. the set D contains all single-vertex cliques {{1}, . . . , {n}}; 4. for every I in D there is a sequence I1 ⊂ · · · ⊂ Ik in D such that |I j | = j for j = 1, . . . k, and Ik = I . We may also use an upper bound U B(ϕ ∗ ) on ϕ ∗ = min I ∈D ϕ(I ), if available. Otherwise we just set U B(ϕ ∗ ) = +∞. For our purpose we use U B(ϕ ∗ ) = min1≤i≤n ϕ({i}) = min1≤i≤n f (ei ). In the case where ϕ(I ) = minx∈∆ f I (x) we take D = C ≻ on the basis of Corollary 3. Note that ϕ(I ) can be evaluated efficiently on D as described in (2). Furthermore, several efficiently computable lower bounds are available for ϕ(I ) = minx∈∆ f I (x) for every I ⊆ N (see [7]). In our implementation we used the bound ℓref Q described in Section 2 of [7]. This is the tightest known bound 2 among those computable in O(n ) time. A. Scozzari, F. Tardella / Discrete Applied Mathematics 156 (2008) 2439–2448 2443 We need to introduce some notation for the description of the algorithm. For i = 1, . . . , n, let Si = {1, . . . , i}, and let N (i) denote the set of all neighbors of vertex i. For increasing values of i, starting from i = 1, the algorithm Exact Clique below searches for the minimum weight clique in D with vertices in Si . More precisely, at each step i the global variable min contains the value of such minimum weight clique in D. The algorithm uses the recursive function Clique (I, U ) described below. For a clique I and a set of vertices U that are adjacent to all vertices of I , Clique (I, U ) updates the current value of min to the value of the minimum weight clique in D that is contained in I ∪ U when the latter value is less than min. E XACT C LIQUE() 1 for i = 1..n 2 do min ← U B(ϕ ∗ ) 3 Clique({i}, Si ∩ N (i)) 4 return C LIQUE(I, U ) 1 if |U | = 0 2 then 3 if ϕ(I ) < min 4 then min ← ϕ(I ) 5 return 6 while |U | 6= 0 7 do if L B(ϕ(I ∪ U )) ≥ min 8 then return 9 else Take k ∈ U 10 U ← U \ {k} 11 I ← I ∪ {k} 12 if I ∈ D 13 then 14 if ϕ(I ) < min 15 then 16 min ← ϕ(I ) 17 18 Clique(I, U ∩ N (k)) 19 else 20 I ← I \ {k} FUNCTION The function Clique (I, U ) is evaluated through implicit enumeration, with a pruning strategy based on the lower bound L B(ϕ(I ∪ U )). To establish its correctness, observe that when |U | = 0, the function Clique (I, U ) correctly updates the value of min when ϕ(I ) < min. On the other hand, when |U | > 0 and L B(ϕ(I ∪ U )) ≥ min, the function Clique (I, U ) terminates without updating min. Indeed, by antitonicity of ϕ, no clique contained in I ∪ U can have a value lower than L B(ϕ(I ∪ U )). Finally, when |U | > 0 and L B(ϕ(I ∪ U )) < min, the procedure tests for each vertex k ∈ U whether I ∪ {k} (in steps 12–16) or any clique containing I ∪ {k} (in the recursive call at step 18) belongs to D and improves the value of min. The efficiency of the algorithm is clearly highly dependent on the order of the vertices, and on the way of choosing vertex k in step 9 of the function Clique (I, U ). Hence, a vertex-coloring heuristic is performed before the start of the algorithm to find a better reordering of the vertices. This tries to place the best candidate vertices for the minimum weight clique in the first positions. Using a coloring that can be found in reasonable time, the vertices are ordered so that those belonging to the same color class are grouped. To get such a coloring, we use the following greedy heuristic due to Biggs [1]: The color classes are determined one at a time. When determining a new color class, the graph induced by the uncolored vertices is first constructed. Then, as long as there exists a vertex that can be added to the 2444 A. Scozzari, F. Tardella / Discrete Applied Mathematics 156 (2008) 2439–2448 color class, such a vertex with largest degree is added. Thus, the vertices are labeled 1, 2, . . . , n in the order they are added to a color class. In addition to the exact algorithm, we also developed a faster heuristic algorithm, called Clique Heuristic(), for the same problem. The procedure starts by evaluating the function ϕ({i}) at each vertex (single-vertex clique) of the graph G. It then proceeds to build the set Cq+1 of all the more promising feasible cliques of dimension q + 1, by trying to enlarge the cliques of Cq . More precisely, for a clique I ∈ Cq and a vertex k 6∈ I , the procedure checks if I ∪ {k} ∈ D, and if the function ϕ(I ∪ {k}) is smaller than the current minimum increased by a tolerance parameter εq that depends on the iteration and on the difference between an upper and a lower bound on the values of ϕ(I ) on D. C LIQUE H EURISTIC() 1 C1 ← {{1}, . . . , {n}} 2 min ← min1≤i≤n ϕ({i}) 3 for q = 1..(n − 1) 4 do Cq+1 = ∅ 5 while Cq 6= ∅ 6 do Take I ∈ Cq 7 Cq ← Cq \ {I } 8 while N (I ) 6= ∅ 9 do Take k ∈ N (I ) 10 N (I ) ← N (I ) \ {k} 11 if I ∪ {k} ∈ D and ϕ(I ∪ {k}) < min +εq 12 then 13 Add I ∪ {k} to Cq+1 14 if ϕ(I ∪ {k}) < min 15 then 16 min ← ϕ(I ∪ {k}) 17 18 return (min) Note that the correctness of both the exact and the heuristic algorithms are based on assumption 4, which guarantees that every clique I in D can be reached with an increasing sequence of cliques in D of all dimensions less than |I |. In the case where D = C ≻ and ϕ(I ) is defined as in (2), the validity of assumption 4 is established by Theorem 6. Let ω(G) denote the clique number of G, i.e., the size of a maximum clique of G. Since both the exact and the heuristic algorithms may analyze each clique of G at most once, the time complexity of both the algorithms are bounded by a polynomial in n and 2ω(G) multiplied by the time required for testing whether a clique I belongs to D, and to evaluate ϕ(I ). Hence the complexity may be exponential in n in the worst case. However, we may easily modify the heuristic algorithm in order to achieve polynomial time complexity. This is done, e.g., by restricting Cq to contain at most M elements for all q. In particular, we may keep only the M elements I in Cq with smallest values of ϕ(I ). In this case, the construction of Cq+1 requires at most Mn executions of step 11, which is usually the most expensive step of the algorithm (depending on how the set D and the function ϕ are defined). Hence, the whole algorithm Clique Heuristic() requires at most Mn 2 executions of step 11. In the case where D = C ≻ and ϕ(I ) is defined as in (2), condition I ∪ {k} ∈ D can be checked in O(n 3 ) time using Cholesky decomposition, and then ϕ(I ∪ {k}) can be evaluated in O(n 2 ) time by solving two triangular systems of linear equations. Hence, in this case the overall time complexity of Clique Heuristic() is O(Mn 5 ). 4. Experimental results In this section we provide some experimental results to compare the efficiency of the exact and heuristic algorithms. Since the algorithms are implemented and executed in the Matlab environment, we also compared our results with the ones obtained by the procedure quadprog available in the Optimization Toolbox of Matlab. This procedure is based on the algorithms described in [9,13], and guarantees only local optimality. Since no standard test problems seem to be available for StQP, we made numerical experiments on random ones. However, to allow comparison with previous studies, such problems have been generated according to the procedures 2445 A. Scozzari, F. Tardella / Discrete Applied Mathematics 156 (2008) 2439–2448 Table 1 Random test problems with known density of the convexity graph Problem n d Algorithm Sec. Cliques Heuristic Sec. Cliques %MRE Matlab Sec. %MRE 10 30 50 100 200 500 1000 1500 10 30 50 100 200 500 1000 1500 10 30 50 100 200 500 1000 1500 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.04 0.10 0.28 1.74 14.74 423.43 – – 0.02 0.47 2.70 50.56 874.71 – – – 0.08 3.92 61.35 1776.6 – – – – 7.6 107.3 379.9 2 683.3 24 194 576 001.2 – – 25.8 879.6 5 507.7 85 100 2262 540.5 – – – 98 10 899 130 250.7 6052 252.3 – – – – 0.03 0.12 0.24 0.70 2.56 55.03 113.91 477.22 0.06 0.28 0.66 2.47 8.81 96.13 1 472.2 2 783.7 0.09 0.60 1.77 6.32 24.6 237.3 5 028.6 26 017 5.3 54.1 105 223.1 503.8 1262.2 1865 2541 17 75.5 138.2 290.3 371 1001 2626.5 1368 22.1 89.7 165.4 300.1 434.6 746 2161 3490 0 0 0 0.24 0 0 – – 0 0 0 0 0.87 – – – 0 0 0 0 – – – – 0.05 0.13 0.28 1.36 18.1 743.4 – – 0.05 0.20 0.34 1.04 14.11 – – – 0.06 0.14 0.34 1.68 – – – – 40.1 42.9 38.3 29.8 36.1 122 – – 24.8 32.4 28.1 23.0 11.0 – – – 19.1 21.5 11.0 13.1 – – – – described in the papers of Nowak [17–19], and of Bomze and De Klerk [5]. We note that Nowak gives computational results concerning his algorithms for finding lower bounds or optimal solutions for StQPs only for dimensions up to n = 30, while Bomze and De Klerk [5] provide some computational results for a polynomial time approximation algorithm for StQPs of dimensions up to n = 40. In fact, most computational experiences reported on StQP concern only local solutions and/or are based on problems obtained as reformulations of maximum (weight) clique problems. We did not use these problems here because it is clearly not very efficient to use an algorithm devised for constrained nonlinear minimum weight clique problems to solve (reformulations of) unconstrained standard maximum (weight) clique problems. Our algorithms have been implemented in Matlab 7.0 under Windows, running on a 2GHz P4, and using standard Matlab functions for the Cholesky decomposition and for the solution of linear systems of equations. We report average results of 10 runs for each problem. In Tables 1–3, n is the problem size, d refers to the density of the convexity graph, Sec. gives the CPU time in −Opt seconds, ME = f min − Opt and %MRE = 100 ∗ fminOpt report on the Maximum Error and the (percentage of the) Maximum Relative Error, respectively, where Opt is the optimal solution found with the exact algorithm, and f min is the solution provided by either the heuristic algorithm or the quadprog procedure of Matlab. As a measure of performance of the exact and heuristic algorithms, we also report the average number of Cliques of order at least two evaluated. We observed that, for a given problem and for a given procedure, the variance of the number of cliques evaluated is close to zero. Hence we did not report it in the tables. Furthermore, in Table 2, s/n represents the fraction of negative eigenvalues of the Hessian matrix C of the quadratic form f with respect to the instance size n. The test problems in Table 1 are generated with the procedure rand qps described in [18,19], which generates StQP problems having specified densities of the associated convexity graphs with unknown global minimizer. In particular, we chose to generate problems with densities 0.25, 0.5, and 0.75 as in [18,19]. The test problems in Table 2 are generated with the procedure random qp described in [17,18], which produces problems with specified number of negative eigenvalues and known global minimizer belonging to a face of specified dimension. Following [17,18], in our experiments we considered 0.2, 0.5 and 0.8 as values of s/n. Furthermore, we require that the global minimizer lies in a face of dimension k = 0.1n, which coincides with half the maximum possible dimension of a face containing 2446 A. Scozzari, F. Tardella / Discrete Applied Mathematics 156 (2008) 2439–2448 Table 2 Random test problems with known optimal solution Problem n s/n d Algorithm Sec. Cliques Heuristic Sec. Cliques %MRE Matlab Sec. %MRE 10 30 50 100 200 500 1000 1500 10 30 50 100 200 500 1000 1500 10 30 50 100 200 500 1000 1500 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.83 0.58 0.68 0.69 0.80 0.81 0.75 0.79 0.70 0.48 0.62 0.60 0.81 0.85 0.81 0.76 0.64 0.46 0.55 0.74 0.73 0.81 0.80 0.78 0.09 0.16 1.64 2.20 644.45 – – – 0.05 0.18 0.90 1.10 610.4 – – – 0.03 0.13 1.19 1.25 601.75 – – – 139.4 61.1 3 349.2 916.6 1046 162 – – – 68 119.8 928.5 864.6 1045 634 – – – 31.9 46.1 2 269.4 930 1047 200 – – – 0.07 0.36 1.50 2.30 12.28 115.14 2 910.4 13 691 0.05 0.40 1.31 2.30 11.0 83.73 1 738.6 11 277 0.02 0.30 1.38 2.38 13.87 58.35 1 578.3 11 888 3.9 92.8 242.7 462.4 1178.3 2139.2 5453 7446.5 3.6 94.1 244.4 453.2 1035 1498.5 4273.5 7735 3 90.5 242.8 462.8 1293.3 1509.8 4425.5 6703 0 0 0 0 0 0.004 0 0.007 0 0 0 0 0 0 0.009 0 0 0 0 0 0 0 0 0 0.15 0.16 0.41 2.07 23.1 – – – 0.05 0.16 0.37 2.0 23.3 – – – 0.05 0.16 0.37 2.03 22.15 – – – 0 0 0 0 0 – – – 0 0 0 0 0 – – – 89.6 0 0 0 0 – – – Table 3 Test problems with matrix entries uniformly distributed in [0, 1] Problem n 10 20 30 40 50 100 200 400 500 1000 1500 2000 3000 4000 Algorithm Sec. 0.02 0.05 0.07 0.17 0.23 3.58 99.81 2310.70 3190.30 – – – – – Cliques 11.3 51.4 60.5 215.3 260.3 4 665.5 139 120 2112 986 2799 900 – – – – – Heuristic Sec. 0.05 0.13 0.27 0.34 0.50 1.70 5.80 17.20 26.10 105.68 311.71 522.76 1449.70 2517.90 Cliques 13.2 33 64.8 87.1 101.7 236 490.6 768.7 1 033.3 2 021.2 3 766 5 104 8 420 10 336 %MRE 0 0 0 0 0 0 0 0 0 – – – – – ME 0 0 0 0 0 0 0 0 0 – – – – – Matlab Sec. 0.04 0.10 0.14 0.22 0.35 1.80 21.81 353.24 812.60 – – – – – %MRE ME 3.6 × 103 0.19 0.22 0.17 0.18 0.23 0.16 0.17 0.14 0.18 – – – – – 8.2 × 103 2.1 × 103 1.4 × 104 5.9 × 103 5.1 × 103 5.8 × 104 3.1 × 104 3.2 × 105 – – – – – a global minimizer when s/n = 0.8. Only for smaller values of n (n = 10, 30, 50) we choose the maximum possible dimension k = 0.2n. Finally, the test problems in Table 3 are generated by taking the entries of the Hessian matrix C uniformly distributed in [0, 1] as described in [5]. For this class of problems, the density d of the convexity graph is always very close to 0.5. For the problems in Tables 1 and 3, we compared the results found by Matlab and by the heuristic algorithm with the optimum found by the exact algorithm, when available. In Table 1, the heuristic algorithm always found an optimal solution except for two cases where it found solutions with very small relative error, while the maximum relative error of the solutions provided by Matlab was between 11% and 122%. In Table 2 we were able to compare 2447 A. Scozzari, F. Tardella / Discrete Applied Mathematics 156 (2008) 2439–2448 Table 4 Accuracy of the heuristic and quadprog procedure Problem P1-0.25 P1-0.50 P1-0.75 P2-0.20 P2-0.50 P2-0.80 P3 ME Heuristic Matlab %MRE Heuristic Matlab Opt Heuristic Matlab 0 0.16 0.35 0 0 0 0 3.25 2.23 1.46 0 0 0 0.25 0 2.63 3.25 0 0 0 0 56.98 32.84 21.67 0 0 0 1.9 × 105 1000 998 990 1000 1000 1000 1000 25 24 47 1000 1000 1000 49 the solutions found by Matlab and the heuristic algorithm with the known global optimum for all dimensions. In this case both Matlab and the heuristic algorithm almost always found the global minimizer, although Matlab was unable to solve the larger problems within 2 hours of CPU time. We suspect however that this favorable behavior is mostly due to the existence of very few non-global local minimizers for this class of test problems. We observe that the performance of our algorithms on this class of problems does not seem to be significantly affected by the percentage of negative eigenvalues. An even better performance of the heuristic algorithm has been recorded in Table 3 for all the test problems where the global minimizer has been found by the exact algorithm. Here the heuristic always found the global minimizer, while this was never the case for Matlab. Clearly, we do not have a performance guarantee for the solutions provided by the heuristic algorithm for large problems in Tables 1 and 3 when the exact algorithm is unable to find a solution. However, based on the evidence with smaller dimensions, we believe that they should be quite accurate. To reinforce this evidence, in Table 4 we report the results concerning 1000 instances of dimension n = 100 of all types of problems considered in Tables 1–3, which are solved with our exact and heuristic algorithms, and with the quadprog procedure in Matlab. In Table 4, P1, P2, and P3 refer to the problems of Table 1, 2, 3, respectively. For P1 and P2 we considered problems with 0.25, 0.5, and 0.75 as graph densities, and with 0.2, 0.5, and 0.8 as fraction of negative eigenvalues, respectively. Here opt counts how many times the heuristic and the quadprog procedure found the optimal solution. Again we observed that the heuristic algorithm almost always found the global minimizer, and that in the few exceptions the gap was very small. We did not compare the results of the heuristic algorithm with known lower bounds of StQP such as those described in [7], because the ones that can be computed efficiently are not very tight, while the best lower bounds, which are based on semidefinite programming, are too hard to compute for n greater than 1000. A selection of StQP test problems of the type described in Tables 1 and 3 together with the global minimum and a global minimizer are available in the web page http://www.myjavaserver.com/˜asco/test StQP for comparison with future computational studies on this problem. 5. Final remarks and perspectives for future research Our approach to StQP consists in restricting the search for the global minimizer of f to a face of the standard simplex associated to a clique of the convexity graph of f , and then performing an implicit search on all such cliques. This approach is based on Theorem 1, that is an extension to QP of the Fundamental Theorem of LP, which is used to restrict the search for global minimizers of an LP to the vertices of the feasible region. Our approach could, in principle, be extended to general quadratic programs. However, its straightforward extension is hampered by the need to determine all vertices and faces of a general polytope, which is a difficult task that might have an exponentially long output. Nevertheless, when a polytope P has a limited number m of vertices, the problem of minimizing a quadratic function over P can be easily transformed into a StQP of dimension m as described in [7]. Thus we can solve such problems with the clique algorithms described in this paper. Examples of problems of this type include problems with few constraints, or the problem of minimizing a quadratic function over the ℓ1 ball ( ) X n B1 = x ∈ R : |xi | ≤ 1 , i 2448 A. Scozzari, F. Tardella / Discrete Applied Mathematics 156 (2008) 2439–2448 which is considered, e.g., in [6,16]. Indeed, the vertices of B1 are the 2n points ±ei for i = 1, . . . , n. Hence, the problem of minimizing a quadratic function on B1 can be transformed into a StQP of dimension 2n (see [6,7] for details on the transformation). Thus the sizes of problems of this kind tractable with our approach should be roughly half the ones reported in Tables 1 and 2. Finally we mention that, while the results we obtained compare favorably with the ones reported in the literature, there is still much space for optimizing the implementation of our algorithms. Some possible directions include using more sophisticated linear algebra techniques for testing positive definiteness of projections of quadratic forms on a sequence of nested linear subspaces; testing different ordering of the vertices and different lower bounds in the Exact Clique algorithm; experimenting with variable size restrictions on the sets Cq in the Clique Heuristic algorithm; and rewriting the code in more efficient programming environments. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] N. Biggs, Some heuristic in graph coloring, in: R. Nelson, R.J. Wilson (Eds.), Graph Colourings, Longman, New York, 1990, pp. 87–96. I.M. Bomze, Global escape strategies for maximizing quadratic forms over a simplex, Journal of Global Optimization 11 (1997) 325–338. I.M. Bomze, On standard quadratic optimization problems, Journal of Global Optimization 13 (1998) 369–387. I.M. Bomze, M. Dür, E. de Klerk, A. Quist, C. Roos, T. Terlaky, On copositive programming and standard quadratic optimization problems, Journal of Global Optimization 18 (2000) 301–320. I.M. Bomze, E. de Klerk, Solving standard quadratic optimization problems via linear, semidefinite and copositive programming, Journal of Global Optimization 24 (2002) 163–185. I.M. Bomze, F. Frommlet, M. Rubey, Improved SDP bounds for minimizing quadratic functions over the l 1 -ball, Optimization Letters 1 (2007) 49–59. I.M. Bomze, M. Locatelli, F. Tardella, New and old bounds for standard quadratic optimization: Dominance, equivalence and incomparability, Mathematical Programming (2007), doi:10.1007/s10107-007-0138-0. I.M. Bomze, L. Palagi, Quartic formulation of standard quadratic optimization problems, Journal of Global Optimization 32 (2005) 181–205. T.F. Coleman, Y. Li, A reflective Newton method for minimizing a quadratic function subject to bounds on some of the variables, SIAM Journal on Optimization 6 (1996) 1040–1058. E. de Klerk, D.V. Pasechnik, A Linear Programming Reformulation of the Standard Quadratic Optimization Problem, CentER Discussion Paper Series No. 2005-24, Tilburg University, The Netherlands, Available at SSRN: http://ssrn.com/abstract=683106, 2005. M. Frank, P. Wolfe, An algorithm for quadratic programming, Naval Research Logistics Quarterly 3 (1956) 95–110. L.E. Gibbons, D.W. Hearn, P.M. Pardalos, M.V. Ramana, Continuous characterizations of the maximum clique problem, Mathematics of Operations Research 22 (1997) 754–768. P.E. Gill, W. Murray, M.H. Wright, Practical Optimization, Academic Press, London, UK, 1981. T.S. Motzkin, E.G. Straus, Maxima for graphs and a new proof of a theorem of Turán, Canadian Journal of Mathematics 17 (1965) 533–540. Y.E. Nesterov, Global quadratic optimization on the sets with simplex structure, Discussion paper 9915, CORE, Katholic University of Louvain, Belgium, 1999. Y.E. Nesterov, H. Wolkowicz, Y. Ye, Nonconvex quadratic optimization, in: H. Wolkowicz, R. Saigal, L. Vandenberghe (Eds.), Handbook of Semidefinite Programming, Kluwer Academic Publishers, Dordrecht, 2000, pp. 361–416. I. Nowak, Some heuristics and test problems for nonconvex quadratic programming over a simplex, Preprint 98–17, Humboldt University, Berlin, 1998. I. Nowak, A global optimality criterion for nonconvex quadratic programming over a simplex, Preprint 98–18, Humboldt University, Berlin, 1998. I. Nowak, A new semidefinite programming bound for indefinite quadratic forms over a simplex, Journal of Global Optimization 14 (1999) 357–364. P.R.J. Östergard, A new algorithm for the maximum-weigth clique problem, Nordic Journal of Computing 8 (2001) 424–436. P.R.J. Östergard, A fast algorithm for the maximum clique problem, Discrete Applied Mathematics 120 (2002) 197–207. F. Tardella, On the equivalence between some discrete and continuous optimization problems, Annals of Operation Research 25 (1990) 291–300. F. Tardella, Connections between continuous and combinatorial optimization problems through an extension of the fundamental theorem of linear programming, Electronic Notes in Discrete Mathematics 17 (2004) 257–262.