A Spectral Bundle Method for Semidefinite Programming.
A central drawback of primal-dual interior point methods for semidefinite programs is their lack of ability to exploit problem structure in cost and coefficient matrices. This restricts applicability to problems of small dimension. Typically, semidefinite relaxations arising in combinatorial applications have sparse and well-structured cost and coefficient matrices of huge order. We present a method that allows us to compute acceptable approximations to the optimal solution of large problems within reasonable time.Semidefinite programming problems with constant trace on the primal feasible set are equivalent to eigenvalue optimization problems. These are convex nonsmooth programming problems and can be solved by bundle methods. We propose replacing the traditional polyhedral cutting plane model constructed from subgradient information by a semidefinite model that is tailored to eigenvalue problems. Convergence follows from the traditional approach but a proof is included for completeness. We present numerical examples demonstrating the efficiency of the approach on combinatorial examples.
Introduction
. The development of interior point methods for semidefinite
programming [19, 31, 1, 46] has increased interest in semidefinite modeling techniques
in several fields such as control theory, eigenvalue optimization, and combinatorial
optimization. In fact, interior point methods proved to be very useful and reliable
solution methods for semidefinite programs of moderate size. However, if the problem
is defined over large matrix variables or a huge number of constraints interior point
methods grow terribly slow and consume huge amounts of memory. The most efficient
methods of today [15, 23, 2, 32, 45, 29] are primal-dual methods that require, in each
iteration of the interior point method, the factorization of a dense matrix of order
equal to the number of constraints and one to three factorizations of the positive
semidefinite matrix variables within the line search. For a typical workstation this
restricts the number of constraints to 2000 and the size of the matrix variables to
500 if reasonable performance is required. For larger problems time and memory
requirements are prohibitive. It is important to realize that either the primal or the
dual matrix is generically dense even if cost and coefficient matrices are very sparse.
Very recently, a pure dual approach was proposed in [4] which offers some possibilities
to exploit sparsity. It is too early to judge the potential of this method.
In combinatorial optimization semidefinite relaxations where introduced in [27].
At that time they were mainly considered a theoretical tool for obtaining strong
bounds [11, 28, 40]. With the development of interior point methods hopes soared
high that these relaxations could be of practical value. Within short time several
approximation algorithms relying on semidefinite programming were published, most
of them based on the approach by Goemans and Williamson [8]. On the implementational
side [14, 16, 20] cutting plane approaches for semidefinite relaxations of
Konrad-Zuse-Zentrum fur Informationstechnik Berlin, Takustrae 7, 14195 Berlin, Germany.
helmberg@zib.de, http://www.zib.de/helmberg
y Universitat Klagenfurt, Institut f. Mathematik, Universitatsstr. 65-67, A 9020 Klagenfurt, Aus-
tria. franz.rendl@uni-klu.ac.at. Financial support through the Austrian FWF Project P12660-
MAT is greatfully acknowledged.
constrained quadratic 0-1 programming problems proved to yield solutions of high
quality. However, as mentioned above, they were very expensive to compute even for
problems of small size (a few hundred 0-1 variables). Problems arising in practical
applications (starting with a few thousand 0-1 variables) were out of reach. We believe
that the method proposed in this paper will open the door to problems of this
size.
Although combinatorial applications are our primary concern we stress that the
method is not restricted to this kind of problems. In fact it will be a useful alternative
to interior point methods whenever the number of constraints or the order of the
matrices is quite large.
We transform a standard dual semidefinite program into an eigenvalue optimization
problem by reformulating the semidefinite constraint as a non-negativity constraint
on the minimal eigenvalue of the slack matrix variable and lifting this constraint
into the cost function by means of a Lagrange multiplier. The correct value
of the Lagrange multiplier is known in advance if the primal feasible matrices have
constant trace. (This is the case for the combinatorial applications we have in mind.)
In this paper we develop a bundle method for solving the problem of minimizing
the maximal eigenvalue of an affine matrix function with an additional linear objective
term. These functions are well known to be convex and non-smooth. A very general
method for optimizing non-smooth convex functions is the bundle method, see e.g.
[21, 42, 17, 18]. In each step the function value and a subgradient of the function is
computed for some specific point. By means of the collected subgradients a cutting
plane model of the function is formed. The minimizer of the cutting plane model
augmented by a regularization term yields the new point. In the case of eigenvalue
optimization the subgradient is formed by means of an eigenvector to the maximal
eigenvalue. Extremal eigenvalues and associated eigenvectors of large symmetric matrices
can be computed efficiently by Lanczos methods (see e.g. [9]). Lanczos methods
need a subroutine that computes the product of the matrix with a vector. This allows
to exploit any kind of structure present in the matrix.
The polyhedral cutting plane model used in traditional bundle algorithms is up-dated
by new subgradient information such as to approximate well the subdifferential,
and thus the function itself, in the vicinity of the current point. For eigenvalue optimization
problems the subdifferential is generated by a semidefinite set, in particular
by the intersection of a simple affine constraint and a face of the semidefinite cone.
This suggests to use, instead of the traditional polyhedral cutting plane model, a
semidefinite cutting plane model that works with an approximation of this face of
the semidefinite cone. This specialization of the cutting plane model is the main
contribution of the paper.
The semidefinite bundle approach allows for an intriguing interpretation in terms
of the original semidefinite program. The cutting plane model requires that the dual
slack matrix of the semidefinite program is positive semidefinite only with respect to a
subspace of vectors, thus it may be interpreted as a relaxation of the dual semidefinite
program. In general the optimal solution of this relaxed semidefinite problem will
produce an indefinite dual slack matrix. One or more of the negative eigenvalues and
corresponding eigenvectors of the slack matrix are used to update the subspace in
order to improve the relaxation and the process is iterated.
This process trivially provides the optimal solution if the subspace grows to the full
space. However, we show that during the algorithm generically the dimension of the
subspace is bounded by (roughly) the square root of the number of constraints. If this
A SPECTRAL BUNDLE METHOD 3
is still considered too large the introduction of an aggregate subgradient guarantees
convergence for restricted bundle sizes. In the extreme the bundle may consist of one
new eigenvector to the maximal eigenvalue only.
In contrast, the 'classical' algorithms of Cullum, Donath, and Wolfe [6] and Polak
and Wardi [38] require in each iteration the computation of all eigenvectors to
eigenvalues within an "-distance of the maximal eigenvalue, thus close to the optimal
solution this number is at least as large as the multiplicity of the maximal eigenvalue
in the optimal solution. In the quadratically convergent algorithm of Overton [35]
each step is computed from a complete spectral decomposition of the matrix and a
guess of the exact multiplicity of the maximal eigenvalue in the optimal solution. In
recent work [33, 34] Oustry reinterprets the algorithm of Overton within the frame-work
of the U-Lagrangian introduced in [26] and embeds it in a first order method to
ensure global convergence. Again, for global convergence the approach relies on the
spectrum of all eigenvalues within "-distance of the maximal eigenvalue and makes
use of the entire spectral information to obtain local quadratic convergence.
Because of the restricted bundle size quadratic convergence is out of reach for
our algorithm, it is a first order method only. In principle convergence follows from
the traditional approach (see e.g. [21]) but we include a proof for completeness. We
also present a primal-dual interior point code for solving the quadratic semidefinite
programming problems associated with the semidefinite cutting plane models and
discuss efficiency aspects. The properties of the algorithm are illustrated on several
combinatorial examples.
In x2 some basic properties of semidefinite programs are stated. Then we transform
semidefinite programs into eigenvalue optimization problems. Section 3 introduces
the bundle method. The algorithm and the proof of convergence is given in
x4. The quadratic semidefinite subproblems arising in the bundle method can be
solved by interior point methods as explained in x5. Section 6 gives an outline of
the implementation and briefly discusses the computation of the maximal eigenvalue
and an associated eigenvector. Numerical examples on combinatorial problems are
presented in x7. We conclude the paper with a summary and possible extensions and
improvements in x8. For the convenience of the reader an appendix explaining the
notation and the symmetric Kronecker product is included at the end of the paper.
2. Semidefinite programs and eigenvalue optimization. We denote the
set of symmetric matrices of order n by Sn which we regard as a space isomorphic
to R ( n+1
As scalar product of A; B 2 Sn (or more general,
use the trace is the sum of the diagonal elements of a
square matrix. We will often use the same symbol for the canonical scalar product
of vectors a; b a, the appropriate space will be clear from the
context. The subset of positive semidefinite matrices S
n is a full-dimensional, non-
polyhedral convex cone in Sn and defines a partial order on the symmetric matrices
by A  B
n . Positive definite matrices are denoted by S ++
n or
A  0.
Consider the standard primal-dual pair of semidefinite programs,
(D)
Z  0:
linear operator and A T its adjoint operator, defined
by hAX; yi
ff for all X 2 Sn and y They are of the form
with A i 2 Sn , is the cost matrix, b 2 R m the right-hand-side
vector.
We assume some constraint qualification to hold, so that these problems satisfy
strong duality in the sense that for any optimal solution X   of (P) and any optimal
solution (y   ; Z   ) of (D) we have
The following assumption allows a simple reformulation of the dual (D) as an
eigenvalue optimization problem. We assume that
for some constant a ? 0. In this case we can add tr a as a redundant constraint
to the primal problem and obtain the following dual equivalent to (D)
Now a ? 0 implies X 6= 0 at the optimum, hence any optimal Z of this dual is
singular. Therefore all dual optimal solutions Z satisfy leading to
Thus we have shown that (D) is equivalent to min y amax
convenience we assume deal with the following problem.
The eigenvalue problem (E) is a convex, non-smooth optimization problem. It is well
studied in the literature. Here we only recall some basic facts. The function
is differentiable if and only if the maximal eigenvalue has multiplicity one. When
optimizing eigenvalue functions, the optimum is generically attained at matrices whose
maximal eigenvalue has multiplicity larger than one. In this case one has to consider
the subdifferential of max at X ,
(see e.g. [35]). In particular, for any v 2 R
n belonging to the eigenspace of the
maximal eigenvalue of X , contained in the subdifferential of max at X .
For the function of interest,
A SPECTRAL BUNDLE METHOD 5
the subdifferential of f at y can be derived by standard rules (see [17]),
Observe that the set of all subgradients is bounded.
Remark 2.1. Even though our assumption (2.2) might look artificial, it does
hold for SDP arising from quadratic 0-1 optimization. It also holds for many other
SDP derived as relaxations of combinatorial optimization problems, see for instance
[1, 12, 24].
3. The bundle method. In this section we develop a new method for minimizing
f . We use two classical ingredients, the proximal point idea, and the bundle
concept. The new contribution lies in the way that we derive the new iterate from
the 'bundle' of subgradient information collected from previous iterates. Since our
approach builds on several subtle ideas, we proceed in small steps and explain first,
how we derive a minorant of f from local information.
3.1. Minorizing f by "
f. Our first goal is to obtain a minorant "
f of f which
approximates f in the neighborhood of the current iterates reasonably well, and which
is easier to handle than f . Introducing the function
we can express f(y) as
This formulation shows that lower approximations of f can be obtained by constraining
W to a subset of all semidefinite matrices with tr
We propose the following choice for this subset. Let P be n \Theta r with P T
n with tr be two matrices. We restrict W to be contained in the
set
c
The
f , defined through P and W , now reads
Wg:
By definition, we have "
f(y)  f(y) 8y: If, for some "
W for some
eigenvector v to
f("y). This is e.g. the case if v is a
column of P or v is contained in the range space of P .
The intuitive idea behind our specific choice of c
W is as follows: the matrix P
contains subgradient information from the current point "
y, and perhaps from previous
iterates. We explain below in detail, how we propose to select and update the
matrix P . For computational efficiency, we would like to keep the number r of columns
of P small, independent of the multiplicity of the largest eigenvalue. Therefore we
collect indispensable subgradient information, that has to be removed from P , in an
aggregate subgradient. This aggregation is the final ingredient of our local model of f .
The matrix W plays the role of an aggregate subgradient. Again, we will discuss be-
low, how W is updated during the algorithm. The main point here is that instead of
optimizing over all semidefinite matrices W , we constrain ourselves to a small subset.
Remark 3.1. If we set use for the matrix P a set of eigenvectors
to the r largest eigenvalues at "
y, we would end up with a model closely related to the
approach from [6]. In this case it would be important to select r at least as large as the
multiplicity of the largest eigenvalue. In our present approach this is not necessary.
6 C. HELMBERG AND F. RENDL
3.2. Proximal point idea. The next goal is to minimize "
f instead of f . Since
f is built from local information from a few previous iterates, this model function is
unlikely to be reliable for points far from the current iterate. Therefore we use the
proximal point idea and add a penalty term for the displacement from the current
point. Thus we determine a new candidate y from the current iterate "
y by solving the
following convex problem, referred to as the augmented model. (Here u ? 0 is some
fixed real weight.)
min y
We note that this minimization problem corresponds to the Lagrangian relaxation of
Thus we replace the original function f by its minorant "
f
and minimize locally around "
y. The weight u controls (indirectly) the radius s of
the sphere around " y, over which we minimize. Substituting the definition of "
f , this
problem is the same as
min y
This problem can be simplified, because y is unconstrained. Note that
Therefore we obtain
min y
W2c W; b\GammaAW +u(y\Gamma"y)=0
W2c
The first equality follows from interchanging min and max (see Corollary 37.3.2 of
[41]) and using first order optimality for the inner minimization with respect to y,
The final problem is a semidefinite program with (concave) quadratic cost function.
We will discuss in x5 how problems of this kind can be solved efficiently. Its optimal
solution W k+1 gives the new trial point y by (3.3).
Remark 3.2. The choice of the weight u is somewhat of an art. There are several
clever update strategies published in the literature, see for instance [21, 42].
3.3. One iteration of the algorithm. The main ingredients of our approach
have now been explained, so we can give a formal description of a general iteration
k of the algorithm. To be consistent with the notation of the algorithm given in x4,
let us denote by x k what was called "
y in x3.2. The algorithm may have to compute
several trial points y k+1 , y keeping the same x progress
is not considered satisfactory (null step). For each y k+1 the function is evaluated and
a subgradient (eigenvector) is computed. This information is added to c
W k to form
an improved model c
W k+1 . Therefore, we assume that the current 'bundle'
A SPECTRAL BUNDLE METHOD 7
contains an eigenvector of its span (y k may or may not be equal
to x k ). Other than that, P need only The minorant of f in
iteration k is denoted by "
Here c
represents the current approximation to the set of all semidefinite matrices
of trace one, see (3.1). It will be convenient to introduce also the regularized version
of "
The new trial point y k+1 is obtained by minimizing f k (y) with respect to y. As
described above, this can be done as follows. First, solve by interior point methods
(see x5)
yielding a (not necessarily unique) maximizer W
use (3.3) to compute
To finish an iteration, we have to decide whether enough progress is made to perform
a serious step or not, i.e. whether we are going to set x
how to update P k and W k
If P k does not yet use the maximum number of columns allowed then the update
process is simple: orthogonalize the new eigenvector with respect to P k , add it as a
new column to form P k+1 and continue. In general, however, P k will already use the
maximum number of columns and so we have to make room for the new subgradient
information. Instead of simply eliminating some columns of P k we can do better by
exploiting the information available in ff   and V   .
Let Q\LambdaQ T be an eigenvalue decomposition of V   . Then the 'important' part of
the spectrum of W k+1 (the important subspace within the space spanned by P k ) is
spanned by the eigenvectors associated with the 'large' eigenvalues of V   . Thus we
split the eigenvectors of Q into two parts (with corresponding spectra   1
and   2 containing as columns the eigenvectors associated to 'large' eigenvalues
of V   and Q 2 containing the remaining columns,
Now the next P k+1 is computed to contain P k Q 1 and at least one eigenvector v k+1
to the maximal eigenvalue of
(The operator orth(.) indicates that we take an orthonormal basis of [P k
The next aggregate matrix is built in such a way that W k+1 2 c
contains only the important part of P k , given by P k Q 1 , we include the remaining
part of P k , given by P k Q 2 in W k+1
(ff   W k
Note that W k+1 is scaled to have trace equal one.
Proposition 3.3. Update rules (3.7) and (3.8) ensure that W k+1 2 c
W k+1 .
Proof. Let W k+1 be of the form (3.6). By (3.7) there is an orthonormal matrix
such that P k+1
W k+1 .
We summarize some easy facts, which will be used in the convergence analysis of
the algorithm.
since y k+1 is minimizer of f k . Because f k
Next let
Using the definition of y k+1 from (3.5) it follows easily that
(y
the augmented model of the next iteration will satisfy
(y) 8y:
Remark 3.4. While the choice for the update of P k is fairly natural, we could use
other update formulas, such as W . The main properties
guiding the update are that W k+1 2 c
ensuring (3.11) and that in y k+1 the model
is now supported by a subgradient of f pushing the model towards f in the vicinity of
the last minimizer.
4. Algorithm and convergence analysis. In the previous section we focused
on the question of doing one iteration of the bundle method. Now we provide a formal
description of the method and point out that except for the choice of the bundle, the
nature of the subproblem, and some minor changes in parameters the algorithm and
its proof are identical to the algorithm of Kiwiel as presented in [21]. To keep the
paper self-contained we present and analyze a simplified variant for fixed u. We refer
the reader to [21] for an algorithm with variable choice of u.
Algorithm 4.1.
Input: An initial point y to the maximal
eigenvalue of C \Gamma A T y 0 , an " ? 0 for termination, an improvement parameter mL 2
an upper bound R  1 on the number of columns of P .
1.
2. (Direction finding) Solve (3.4) to get y k+1 from (3.5). Decompose V   into
using
(3.8).
3. (Evaluation) Compute and an eigenvector v k+1 . Compute
P k+1 by (3.7).
4. (Termination) If f(x k
A SPECTRAL BUNDLE METHOD 9
5. (Serious step) If
then set x continue with Step 7. Otherwise continue with Step 6.
6. (Null step) Set x
7. Increase k by 1 and go to Step 2.
We prove convergence of the algorithm for If the algorithm stops after a
finite number of iterations then by
and thus by (3.5) 0 2 @f(x k ), so x k is optimal. Assume in the following that the
algorithm does not stop. First consider the case that only null steps occur after some
iteration K.
Lemma 4.2. If there is a K  0 such that (4.1) is violated for all k  K, then
Proof. For convenience we set Using the relations (3.10),
and (3.9), we obtain for all k  K
(y
Therefore the f k (y k+1 ) converge to some f    f(x) and
the computed gradient of f in y k+1 and observe
that the linearization
f of f in y k+1
ff
W k+1 . Thus
\Gamma\Omega
ff
The convergence of the f k (y k+1 ), the boundedness of the gradients and the fact that
imply that the last term goes to zero for k !1. So for all
there is an M 2 N such that for all k ? M
where '!' follows from (4.1) being violated for all k ? K. Thus the sequences f(y k+1 )
both converge to f(x). y k+1 is the minimizer of the regularized function
f k . On the one hand this implies that y k+1 ! x. On the other hand 0 must
be contained in the subgradient @f k (y k+1
Therefore there is a sequence h k 2 @ "
subgradients converging to zero. The
converge to f(x) and the y k+1 converge to x, hence zero
must be contained in @f(x).
We may concentrate on serious steps in the following. In order to simplify notation
we will speak of x k as the sequence generated by serious steps with all duplicates
eliminated. By f k (and the corresponding "
refer to the function whose
minimization gives rise to x k+1 .
The next lemma investigates the case that the f(x k ) remain above some value
f(~x) for some fixed ~
x.
Lemma 4.3. If
fixed ~
m and all k
then the x k converge to a minimizer of f .
Proof. First we prove the boundedness of the x k . To this end denote by g k+1 2
subgradient arising from the optimal solution of the minimization problem
observe that by (3.3)
Therefore the distance of x k+1 to ~
x can be bounded by
2\Omega ~
ff
2\Omega ~
ff
For any k ? K, a recursive application of the bound above yields
uX
By (4.1) the progress of the algorithm in each serious step is at least mL (f(x k
together with (4.2) we obtainX
Therefore the sequence of the x k remains bounded and has an accumulation point
x. By replacing ~ x by
x in (4.4) and choosing K sufficiently large, the remaining sum
can be made smaller than an arbitrary small ffi ? 0, thus proving the convergence of
the x k to  x. As the x k+1 converge to
x the g k+1 converge to zero by (4.3), and since
the sequence (f(x k
has to converge to zero as well, we conclude that
x is a minimizer of f .
The lemma also implies that f(x k there are no minimizers. We
summarize the discussion in the following theorem.
Theorem 4.4. [21] If the set of minimizers of f is not empty then the x k converge
to a minimizer of f . In any case f(x k
A SPECTRAL BUNDLE METHOD 11
Remark 4.5. We have just seen that the bundle algorithm works correctly even
if P contains only one column. In this case the use of the aggregate subgradient is
crucial.
To achieve correctness of the bundle algorithm without aggregate subgradients, it
suffices to store in P only the subspace spanning the eigenvectors corresponding to
non-zero eigenvalues of an optimal solution W k+1 of (3.2). Using the bound of [36] it
is not too difficult to show that in this case the maximal number of columns one has to
provide is the largest
plus the number of eigenvectors
to be added in each iteration (this is at least one). In our computational experiments
we found that this upper bound is hardly ever reached. In fact, typical values for the
maximal rank are around half this upper bound.
5. Solving the subproblem. In this section we concentrate on how the minimizer
of f k can be computed efficiently. We have already seen in x3 that this task
is equivalent to solving the quadratic semidefinite program (3.4). Problems of this
kind can be solved by interior point methods, see e.g. [7, 23]. Dropping the iteration
index k and the constants in (3.4) we obtain for
ff
ff  0; V  0:
Expanding into the cost function yields
ff
ff
ff  0; V  0:
Using the svec-operator (see the appendix for a definition and important properties
of svec and the symmetric Kronecker
product\Omega s ) to expand symmetric matrices
from S r into column vectors of length
we obtain the quadratic program (recall
that, for
I svec(V
where (after some technical linear algebra)
ff
ff
+\Omega C; W
At this point it is advisable to spend some thought on W . The algorithm is designed
for very large and sparse cost matrices C. W is of the same size as C. Initially it
might be possible to exploit the low rank structure of W for efficient representations,
but as the algorithm proceeds, the rank of W grows inevitably. Thus it is impossible
to store all the information of W . However, as we can see in (5.2) to (5.6), it suffices to
have available the vector AW 2 R
m and the
scalar\Omega C; W
ff to construct the quadratic
program. Furthermore, by the linearity of A(\Delta) and hC; \Deltai, these values are easily
updated whenever W is changed.
To solve (5.1) we employ a primal-dual interior point strategy. To formulate the
defining equations for the central path we introduce a Lagrange multiplier t for the
equality constraint, a dual slack matrix U  0 as complementary variable to V , a dual
slack scalar fi  0 as complementary variable to ff and a barrier parameter  ? 0.
The system reads
ts I \Gamma
I svec(V
The step direction (\Deltaff; \Deltafi; \DeltaU; \DeltaV; \Deltat) is determined via the linearized system
I svec(\DeltaV
In the current context we prefer the linearization
because it makes the system easy to solve for \DeltaV with relatively little computational
work per iteration. The final system for \DeltaV reads
ff
It is not too difficult to see that the system matrix is positive definite (because
suffices to show that Q
using
0). The main work per iteration is the factorization of this matrix
(with v 2 S r this is
it is not possible to do much better since Q 11
has to be inverted at some point. Because of the strong dominance of the factorization
it pays to employ a predictor corrector approach, but we will not delve into this here.
For strictly feasible primal starting point is
a strictly feasible dual starting point can be constructed by choosing t 0 sufficiently
negative such that
A SPECTRAL BUNDLE METHOD 13
Starting from this strictly feasible primal-dual pair we compute the first  by
compute the step direction (\Deltaff; \Deltafi; \DeltaU; \DeltaV; \Deltat) as indicated
above , perform a line search with line search parameter
strictly feasible, move to this
new point, compute a new  by
ae
oe
with
and iterate. We stop if (hU;
6. Implementation. In our implementation of the algorithm we largely follow
the rules outlined in [21]. In particular u is adapted during the algorithm. The first
guess for u is equal to the norm of the first subgradient determined by v 0 . The
scheme for adapting u is the same as in [21] except for a few changes in parameters.
For example the parameter mL for accepting a step as serious is set to
the parameter mR indicating that the model is so good (progress by the serious step
is larger than mR [f(x k
that u can be decreased is set to
The stopping criterion is formulated in relative precision,
in the implementation.
The choice of the upper bound R on the number of columns r of P and the
selection of the subspace merits some additional remarks. Observe that by Remark 4.5
it is highly unlikely the r violates the bound
even if the number of columns
of P is not restricted.
is also the order of the system matrix in (5.7) and is
usually considerably smaller than the size of the system matrix in traditional interior
point codes for semidefinite programming which is of order m. Furthermore the order
of the matrix variables is r as compared to n for traditional interior point codes.
Thus if the number of constraints m is roughly of the same size as n and a matrix of
order m is still considered factorizable then running the algorithm without bounding
the number of columns of P may turn out to be considerably faster than running an
interior point method. This can be observed in practice, see x7.
For huge n and m primal-dual interior point methods are not applicable any
more, because X , Z \Gamma1 , and the system matrix are dense. In this case the proposed
bundle approach allows to apply the powerful interior point approach at least on an
important subspace of the problem. The correct identification of the relevant subspace
in V   is facilitated by the availability of the complementary variable U   . U   helps to
discern between the small eigenvalues of V   (because of the interior point approach
we have V    0!). Eigenvectors v of V   that are of no importance for the optimal
solution of the subproblem will have a large value v T U   v, whereas eigenvectors, that
are ambiguous, will have both, a small eigenvalue v T V   v and a small value v T U   v.
In practice we restrict the number of columns of P to 25 and provide room for at
least five new vectors in each iteration (see below). Eigenvectors v that correspond
to small but important eigenvalues  of V
are added to W ; important eigenvectors are added to W only if
more room is needed for new vectors.
For large m the computation of (5.2) to (5.6) is quite involved. A central object
appearing in all constants is the projection of the constraint A i on the space spanned
14 C. HELMBERG AND F. RENDL
by P , P T A i P . Since the A i are of the same size as X which we assume to be huge, it
is important to exploit whatever structure is present in A i to compute this projection
efficiently. In combinatorial applications the A i are of the form vv T with v sparse and
the projection can be computed efficiently. In the projection step and in particular in
forming Q 11 the size of r is again of strong influence. If we neglect the computation of
the computation of Q 11 still requires 2m
flops. Indeed, if m is
large then for small r the construction of Q 11 takes longer than solving the associated
quadratic semidefinite program.
The large computational costs involved in the construction and solution of the
semidefinite subproblems may lead to the conviction that this model may not be worth
the trouble. However, the evaluation of the eigenvalue-function is in fact much more
expensive. There has been considerable work on computing eigenvalues of huge, sparse
matrices, see e.g. [9] and the references therein. For extremal eigenvalues of symmetric
matrices there seems to be a general consensus, that Lanczos type methods work best.
Iterative methods run into difficulties if the eigenvalues are not well separated. In our
context it is to be expected that in the course of the algorithm the largest eigenvalues
will get closer and closer till all of them are identical in the optimum. For reasonable
convergence block Lanczos algorithms with blocksize corresponding to the largest
multiplicity of the eigenvalues have to be employed. During the first ten iterations
the largest eigenvalue is usually well separated and the algorithm is fast. But soon
the eigenvalues start to cluster, larger and larger blocksizes have to be used, and
the eigenvalue problem gets more and more difficult to solve. In order to reduce the
number of evaluations it seems worth to employ powerful methods in the cutting plane
model. The increase in computation time required to solve the subproblem goes hand
in hand with the difficulty of the eigenvalue problem because of the correspondence
of the rank of P and the number of clustered eigenvalues.
Iterative methods for computing maximal eigenvectors generically offer approximate
eigenvectors to several other large eigenvalues, as well. The space spanned
by these approximate eigenvectors is likely to be a good approximation of the true
eigenspace. If the maximal number of columns for P is not yet attained it may be
worth to include several of these approximate eigenvectors as well.
In our algorithm we use a block Lanczos code of our own that is based on a
Fortran code of Hua (we guess that this is Hua Dai of [47]). It works with complete
orthogonalization and employs Chebyshev iterations for acceleration. The choice of
the blocksize is based on the approximate eigenvalues produced by previous evaluations
but is at most 30. Four block Lanczos steps are followed by twenty Chebyshev
iterations. This scheme is repeated till the maximal eigenvalue is found to the required
relative precision. The relative precision depends on the distance of the maximal to
the second largest eigenvalue but is bounded by 10 \Gamma6 . As starting vectors we use the
complete block of eigenvectors and Lanczos-vectors from the previous evaluation.
7. Combinatorial applications. The combinatorial problem we investigate is
quadratic programming in f\Gamma1; 1g variables,
In the case that C is the Laplace matrix of a (possible weighted) graph the problem
is known to be equivalent to the max-cut problem.
The standard semidefinite relaxation is based on the identity x T Cx
ff .
For all f\Gamma1; 1g n vectors, xx T is a positive semidefinite matrix with all diagonal elements
equal to one. We relax xx T to X  0 and and obtain the following
A SPECTRAL BUNDLE METHOD 15
primal-dual pair of semidefinite programs,
Z  0:
For non-negatively weighted graphs a celebrated result of Goemans and Williamson
[8] says, that there is always a cut within :878 of the optimal value of the relaxation.
One of the first attempts to approximate (DMC) using eigenvalue optimization
is contained in [39]. The authors use the Bundle code of Schramm and Zowe [42]
with a limited number of bundle iterations, and so do not solve (DMC) exactly. So
far the only practical algorithms for computing the optimal value were primal-dual
interior point algorithms. However these are not able to exploit the sparsity of the
cost function and have to cope with dense matrices X and Z \Gamma1 . An alternative
approach based on a combination of the power method with a generic optimization
scheme of Plotkin, Shmoys, and Tardos [37] was proposed in [22] but seems to be
purely theoretical.
In

Table

7.1 we compare the proposed bundle method to our semidefinite primal-dual
interior point code of [14] (called PDIP in the sequel) for graphs on
nodes that were generated by rudy, a machine independent graph generator written
by G. Rinaldi. Table 7.7 contains the command lines specifying the graphs. Graphs
G 1 to G 5 are unweighted random graphs with a density of 6% (approx. 19000 edges).
G 6 to G 10 are the same graphs with random edge weights from f\Gamma1; 1g. G 11 to G 13
are toroidal grids with random edge weights from f\Gamma1; 1g (1600 edges). G 14 to G 17
are unweighted 'almost' planar graphs having as edge set the union of two (almost
maximal) planar graphs (approx. 4500 edges). G to G 21 are the same almost planar
graphs with random edge weights from f\Gamma1; 1g. In all cases the cost matrix C is the
Laplace matrix of the graph divided by 4, i.e., let A denote the (weighted) adjacency
matrix of G, then
For a description of the code PDIP see [14], the termination criterion requires the
gap between primal and dual optimal solution to be closed to a relative accuracy
of
For the bundle algorithm, (DMC) is transformed into an eigenvalue optimization
problem as described in x2. In addition the diagonal of C is removed so that, in fact,
the algorithm works on the problem
min
with
This does not change problem (PMC) because the
diagonal elements of X are fixed to one. The offset 1e T (Ae \Gamma diag(A)) is added to
the output only and has no influence on the algorithm whatsoever, in particular it
has no influence on the stopping criterion. As starting vector y 0 we choose the zero
vector. All other parameters are as described in x6.
All computation times, for the interior point code PDIP as well as for the bundle
code, refer to the same machine, a Sun sparc Ultra 1 with a Model 140 UltraSPARC
CPU and 64 MB RAM. The time measured is the user time and it is given in the
leading zeros are dropped.
The first column of Table 7.1 identifies the graphs. The second and third refer
to PDIP and contain the optimal objective value produced (these can be regarded as
highly accurate solutions) and the computation time. The fourth and fifth column
give the same numbers for the bundle code.
On these examples the bundle code is superior to PDIP. Although the examples
do belong to the favorable class of instances having small m and relatively large n,
the difference in computation time is astonishing. Note that the termination criterion
used in the bundle code is quite accurate, except for G 11 which seems to be a difficult
problem for the bundle method. This deviation in accuracy is not caused by cancellations
in connection with the offset. The difficulty of an example does not seem to
depend on the number of nonzeros but rather on the shape of the objective function.
For toroidal grid graphs the maximum cut is likely to be not unique, thus the objective
function will be rather flat. This flatness has its effect on the distribution of the
eigenvalues in the optimal solution. Indeed, for G 11 more eigenvalues cluster around
the maximal eigenvalue than for the other problems. We illustrate this in Table 7.2,
which gives the largest eigenvalues of the solution at termination for problems G 1 ,
G 6 , G 11 , G 14 , and G

Table
Comparison of the interior point (PDIP) and the bundle (B) approach. sol gives the computed
solution value and time gives the computation time.
PDIP-sol PDIP-time B-sol B-time
G2 12089.43 1:19:14 12089.45 5:19
G6 2656.16 1:24:53 2656.18 3:57
G 11 629.16 1:28:41 629.21 45:26
G
G 17
G
G 19

Table

7.3 provides additional information on the performance of the bundle algorithm
on the examples of Table 7.1. The second column gives the accumulated time
spent in the eigenvalue computation, it accounts for roughly 90% of the computation
time. serious displays the number of serious steps, iter gives the total number of
iterations including both, serious and null steps. kgk is the norm of the subgradient
arising from the last optimal W k+1 before termination. For G 11 the norm is considerably
higher than for all other examples. Since the desired accuracy was not achieved
for G 11 by the standard stopping criterion it may be worth to consider an alternative
stopping criterion taking into account the norm of the subgradient as well. Column
max-r gives the maximal rank of P attained over all iterations. The rank of P would
A SPECTRAL BUNDLE METHOD 17

Table
The maximal eigenvalues after termination of examples G 1 , G6 , G11 , G14 , and G18 .
13 3.1190 3.2239 0.7651 1.0557 1.4047
14 3.1135 3.2181 0.7650 1.0515 1.4007
19 2.7214 2.7716 0.7647 1.0398 1.3725
22 2.6834 2.6756 0.7644 1.0341 1.3583
26 2.6274 1.8722 0.7636 1.0239 1.3480
28 2.5137 1.7974 0.7633 1.0211 1.3397
29 2.4840 1.4859 0.7630 1.0180 1.3345
have been bounded by 25, but this bound never came into effect for any of these
examples. Aggregation was not necessary. Observe that the theoretic bound allows
for r up to 39, yet the maximal rank is only half this number. The last column gives
the time when the objective value was first within 10 \Gamma3 of the optimum.
For combinatorial applications high accuracy of the optimal solution is of minor
importance. An algorithm should deliver a reasonable bound fast and its solution
should provide some hint on how a good feasible solution can be constructed. The
bundle algorithm offers both. With respect to computation time the bundle algorithm
displays the usual behavior of subgradient algorithms. Initially progress is very fast,
but as the bound approaches the optimum there is a strong tailing off effect. We
illustrate this by giving the objective values and computation times for the serious
steps of example G 6 (the diagonal offset is +77 in this example) in Table 7.4. After
one minute the bound is within 0:1% of the optimum. For the other examples see the
last column of Table 7.3.
With respect to a primal feasible solution observe that P k V k
successively
better and better approximation to the primal optimal solution X   . In case
too much information is stored in the aggregate vector AW k (remember that it is not
advisable to store W k itself), P k may be enriched with additional Lanczos-vectors
from the eigenvalue computation. The solution of this enlarged quadratic semidefinite
subproblem will be an acceptable approximation of X   . It is not necessary to

Table
Additional information about the bundle algorithm for the examples of Table 7.1. -time gives
the total amount of time spent for computing the eigenvalues and eigenvectors, serious gives the
number of serious steps, iter the total number of iterations including null steps. kgk refers to the
norm of the gradient resulting from the optimal solution of the last semidefinite subproblem. max-r
is the maximum number of columns used in P (the limit would have been 25). 0.1%-time gives the
time when the bound is within 10 \Gamma3 of the optimum in relative precision.
-time serious iter kgk max-r 0.1%-time
G1 3:12 22 33 0.1639
G4 2:38 19 27 0.08235 19 54
G 9
2:59
G13 17:24 43 78 0.218 15 6:17
G
G 19 11:24 41 71 0.1571 15 3:34
construct the whole matrix X   . In fact, the factorized form (P k
much more convenient to work with. For example the approximation algorithm of
Goemans and Williamson [8] requires precisely this factorization. A particular x ij
element of X   is easily computed by the inner product of row i and j of the n \Theta r
. In principle this opens the door for branch and cut approaches to
improve the initial relaxation. This will be the subject of further work.

Table

7.5 gives a similar set of examples for
A last set of examples is devoted to the Lov'asz #-function [27] which yields an
upper bound on the cardinality of a maximal independent (or stable) set of a graph.
For implementational convenience we use its formulation within the quadratic f\Gamma1; 1g
programming setting, see [24]. For a graph with k nodes and h edges we obtain a
semidefinite program with matrix variables of order
constraints. The examples we are going to consider have more than one thousand
nodes and more than six thousand edges. For these examples interior point methods
are not applicable any more because of memory requirements. It should be clear
from the examples of Table 7.5 that there is also little hope for the bundle method to
terminate within reasonable time. However, the most significant progress is achieved
in the beginning and for the bundle method memory consumption is not a problem.
We run these examples with a time limit of five hours. More precisely, the algorithm
is terminated after the first serious step that occurs after five hours of computation
time.
The graph instances are of the same type as above. The computational results are
displayed in Table 7.6. The new columns n and m give the order of the matrix variable
and the number of constraints, respectively. Observe that the toroidal grid graphs
A SPECTRAL BUNDLE METHOD 19

Table
Detailed account of the serious steps of example G 6 .
iter value time kgk max-r
8 2670.10 28 5.992 15
9 2666.07 34 4.173
14 2656.31 1:57 0.431
19 2656.19 3:33 0.07243

Table
Examples for
B-sol B-time -time serious iter kgk max-r
G22 14135.98 38:11 28:00 26 52 0.0781 23
G26 14132.93 34:45 26:37 31 48 0.3066 23
G 28 4100.81 29:41 21:08 23
G
G
G 36 8006.04 2:56:10 2:31:09 62 115 0.2634 24
G 37
G 38 8015.01 4:03:53 3:39:24 58 155 0.1937 22
G
G 48 and G 49 are perfect with independence number 1500; the independence number
of G 50 is 1440 but G 50 is not perfect. We do not know the independence number of
the other graphs. Except for G 48 and G 49 , which have '(G 48
perfectness, it is hard to judge the quality of the solutions. Tracing the development of
the bounds the last serious steps of examples G 43 to G 47 and G 51 to G 54 still produced
improvements of 0.5% to 1%. This and the rather large norm of the subgradient of

Table
Upper bound on the #-function after five hours of computation time.
serious iter kgk max-r
G 43
G 44 1001 10991 310.13 5:06:31 3:14:25
G48 3001 9001 1526.53 5:11:31 4:59:57 54 94 0.4062 15
G
G 50 3001 9001 1536.12 5:17:51 5:01:25 50 124 0.4728 15
G 53 1001 6915 463.86 5:08:36 4:36:34 41 104 2.593 25
G 51 and G 54 indicate that the values cannot be expected to be 'good' approximations
of the #-function. Also note, that the size of the subspace required for G 48 to G 50 is
still well below 25. In examples G 51 to G 54 the value of ff is almost negligible, but
for G 43 to G 47 the value of ff is roughly 1=3 at termination. Thus for these examples
the restriction to 25 columns became relevant.
The computational results of Table 7.6 demonstrate that the algorithm has its
limits. Nonetheless the bounds obtained are still useful and the primal approximation
corresponding to the subgradient is a reasonable starting point for primal heuristics.
8. Conclusions and extensions. We have proposed a proximal bundle method
for solving semidefinite programs with large sparse or strongly structured coefficient
matrices. The semidefinite constraint is lifted into the objective function by means of a
Lagrange multiplier a whose correct value is not known in general, except for problems
with fixed primal trace. In the latter case a is precisely the value of the trace. The
approach differs from previous bundle methods in that the subproblem is tailored
for semidefinite programming. In fact the whole approach can be interpreted as
semidefinite programming over subspaces where the subspace is successively corrected
and improved till the optimal subspace is identified. The set of subgradients modeled
by the semidefinite subproblem is a superset of the subgradients used in the traditional
polyhedral cutting plane model. Therefore convergence of the new method is a direct
consequence of previous proofs for traditional bundle methods. It is not yet clear
whether the specialized model admits stronger convergence results. The choice of u
is still very much an open problem of high practical importance.
For (constrained) quadratic f\Gamma1; 1g-programming the method offers a good bound
within reasonable time and allows to construct an approximate primal optimal solution
(of the relaxation) in compact representation. To improve the bound by a cutting
plane approach the algorithm must be able to deal with sign constraints on the
y-variables. In principle it is not difficult to model the sign constraints in the semidefinite
subproblem. However, as a consequence the influence of the sign constrained y
variables on the cost coefficients of the quadratic subproblem cannot be eliminated
any longer, rendering the method impractical even for a moderate number of cutting
planes. Alternatively one might consider active set methods but these entail the danger
of destroying convergence. Together with K.C. Kiwiel we are currently working
on alternative methods for incorporating sign constraints on y [13].
The backbone of the method is an efficient routine for computing the maximal
eigenvalue of huge structured symmetric matrices. Although our own implementation
A SPECTRAL BUNDLE METHOD 21

Table
Arguments for generating the graphs by the graph generator rudy.
G1 -rnd graph 800 6 8001
G2 -rnd graph 800 6 8002
G3 -rnd graph 800 6 8003
G4 -rnd graph 800 6 8004
G5 -rnd graph 800 6 8005
G6 -rnd graph 800 6 8001 -random
G7 -rnd graph 800 6 8002 -random
G8 -rnd graph 800 6 8003 -random
G9 -rnd graph 800 6 8004 -random
G12 -toroidal grid 2D 50
G13 -toroidal grid 2D 25
G14 -planar 800 99 8001 -planar 800
G17 -planar 800 99 8007 -planar 800
G22 -rnd graph 2000 1 20001
G23 -rnd graph 2000 1 20002
G24 -rnd graph 2000 1
G25 -rnd graph 2000 1 20004
G26 -rnd graph 2000 1 20005
G27 -rnd graph 2000 1 20001 -random
G28 -rnd graph 2000 1 20002 -random
G29 -rnd graph 2000 1 -times
G31 -rnd graph 2000 1 20005 -random
G32 -toroidal grid 2D 100 20 -random 0 1 -times
G35 -planar 2000 99 20001 -planar 2000
G36 -planar 2000 99 -planar 2000
G37 -planar 2000 99 20005 -planar 2000
G38 -planar 2000 99 20007 -planar 2000
-planar 2000 99 20001 -planar 2000
G40 -planar 2000 99 -planar 2000
-planar 2000 99 20005 -planar 2000 -times 2 -plus
G42 -planar 2000 99 20007 -planar 2000
G43 -rnd graph 1000 2 10001
G44 -rnd graph 1000 2 10002
G45 -rnd graph 1000 2 10003
G46 -rnd graph 1000 2 10004
G47 -rnd graph 1000 2 10005
G48 -toroidal grid 2D 50
-toroidal grid 2D
G50 -toroidal grid 2D 25 120
G52 -planar 1000 100 10003 -planar 1000 100 10004
G53 -planar 1000 100 10005 -planar 1000 100 10006
G54 -planar 1000 100 10007 -planar 1000 100 10008
(based on the code of Hua) seems to work sufficiently stable there is certainly much
room for improvement. A straight forward approach to achieve serious speed-ups is to
implement the algorithm on parallel machines, see for instance [43]. Rather recently
interest in the Lanczos method has risen again, see [25, 3, 5, 10, 30] and references
therein. Most of these papers are based on the concept of an implicit restart proposed
in [44] which is a polynomial acceleration approach that does not require additional
matrix vector multiplications. It will be interesting to test these new ideas within the
bundle framework.
We thank K.C. Kiwiel for fruitful discussions and C. Lemar'echal and an anony-
22 C. HELMBERG AND F. RENDL
mous referee for their constructive critisism that helped to improve the presentation.


Appendix

. Notation.
R
real column vector of dimension n
real matrices
real matrices
positive definite matrices
positive semidefinite matrices
A  0 A is positive definite
A  0 A is positive semidefinite
I , I n identity of appropriate size or of size n
e vector of all ones of appropriate dimension
maximal eigenvalue of A
tr A trace of A 2 M n;n , tr
product in Mm;n ,
dimensional vector representation of A 2 Sn
Kronecker product of A 2 M m;n
diag(A) the diagonal of A 2 Mn as a column vector
diagonal matrix with v on its main diagonal
Sn is isomorphic to R
via the map svec(A) defined by stacking the columns of
the lower triangle of A on top of each other and multiplying the offdiagonal elements
with
2,
a
The factor
for offdiagonal elements ensures that, for
The symmetric Kronecker
product\Omega s is defined for arbitrary square matrices
M n;n by its action on a vector svec(C) for a symmetric matrix C 2 Sn ,
Both concepts were first introduced in [2]. Here we use the notation introduced in
[45]. From the latter paper we also cite some properties of the symmetric Kronecker
product for the convenience of the reader.
1.
B\Omega s A
2.
(A\Omega s B)
3.
A\Omega s I is symmetric if and only if A is.
4.
5.
(A\Omega s
6. If A  0 and B  0 then
(A\Omega s B)  0
7.



--R

Interior point methods in semidefinite programming with applications to combinatorial optimization

Iterative methods for the computation of a few eigenvalues of a large symmetric matrix
Solving large-scale sparse semidefinite programs for combinatorial optimization
An implicitly restarted Lanczos method for large symmetric eigenvalue problems
The minimization of certain nondifferentiable sums of eigenvalues of symmetric matrices

Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming
Matrix Computations
A shifted block Lanczos algorithm for solving sparse symmetric generalized eigenproblems
Geometric Algorithms and Combinatorial Optimization
Fixing variables in semidefinite relaxations
Incorporating inequality constraints in the spectral bundle method

An interior-point method for semidefinite programming
Quadratic knapsack relaxations using cutting planes and semidefinite programming
Convex Analysis and Minimization Algorithms I

An interior-point method for minimizing the maximum eigenvalue of a linear combination of matrices
Solving graph bisection problems with semidefinite programming
Proximity control in bundle methods for convex nondifferentiable minimization
Efficient approximation algorithms for semidefinite programs arising from MAXCUT and COLORING

Connections between semidefinite relaxations of the max-cut and stable set problems
--TR

--CTR
Jiahai Wang, Letters: An improved discrete Hopfield neural network for Max-Cut problems, Neurocomputing, v.69 n.13-15, p.1665-1669, August, 2006
Samuel Burer , Renato D. C. Monteiro , Yin Zhang, Interior-Point Algorithms for Semidefinite Programming Based on a Nonlinear Formulation, Computational Optimization and Applications, v.22 n.1, p.49-79, April 2002
Abraham Duarte , ngel Snchez , Felipe Fernndez , Ral Cabido, A low-level hybridization between memetic algorithm and VNS for the max-cut problem, Proceedings of the 2005 conference on Genetic and evolutionary computation, June 25-29, 2005, Washington DC, USA
Gerald Gruber , Franz Rendl, The bundle method for hard combinatorial optimization problems, Combinatorial optimization - eureka, you shrink!, Springer-Verlag New York, Inc., New York, NY,
Tijl De Bie , Nello Cristianini, Fast SDP Relaxations of Graph Cut Clustering, Transduction, and Other Combinatorial Problems, The Journal of Machine Learning Research, 7, p.1409-1436, 12/1/2006
Kazuhide Nakata , Makoto Yamashita , Katsuki Fujisawa , Masakazu Kojima, A parallel primal-dual interior-point method for semidefinite programs using positive definite matrix completion, Parallel Computing, v.32 n.1, p.24-43, January 2006
Hernn Alperin , Ivo Nowak, Lagrangian Smoothing Heuristics for Max-Cut, Journal of Heuristics, v.11 n.5-6, p.447-463, December  2005
Alper Yildirim , Xiaofei Fan-Orzechowski, On Extracting Maximum Stable Sets in Perfect Graphs Using Lovsz's Theta Function, Computational Optimization and Applications, v.33 n.2-3, p.229-247, March     2006
Stephen Braun , John E. Mitchell, A Semidefinite Programming Heuristic for Quadratic Programming Problems with Complementarity Constraints, Computational Optimization and Applications, v.31 n.1, p.5-29, May       2005
Henry Wolkowicz , Miguel F. Anjos, Semidefinite programming for discrete optimization and matrix completion problems, Discrete Applied Mathematics, v.123 n.1-3, p.513-577, 15 November 2002
proximal bundle method;eigenvalue optimization;semidefinite programming;convex optimization;large-scale problems