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)
minf 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)
minf 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.