Interior Point Methods in Semidefinite Programming with Applications to Combinatorial Optimization ^{∗}

Farid Alizadeh^{†}

Abstract We study the semidefinite programming problem (SDP), that is the problem of minimizing a linear function of a symmetric matrix X subject to linear constraints on X and the additional condition that X be positive semidefinite. We develop a duality theory for SDP similar to that of linear programming. We extend Ye’s projective interior point method for linear programming to solve the SDP problem. We also show that various eigenvalue optimization problems may be formulated as SDP problems. Finally, we show that many combinatorial optimization problems may be approximated (and even sometimes exactly expressed) by an SDP problem. In particular we show how our interior point approach may be applied to a new class of semidefinite matrix cuts for 0–1 integer programming developed by Lovász and Schrijver. We also apply our SDP interior point algorithm to the maximum stable set problem, maximum k–partite induced subgraph problem, and graph partitioning problems.

1 Introduction.

Consider the following optimization problem which we call the standard semidefinite programming problem (SDP):

min{C • X : A_{i} • X = b_{i} for i = 1, ⋯, m and X ⪰ 0}  (1.1)

where C, A_{i}’s and X are n × n matrices, and X is symmetric; the “•” operation is the inner product of matrices: A • B := ∑_{i,j} A_{ij}B_{ij} = trace A^{T} B; and the “inequality” constraint ⪰ indicates the Löwner partial order, that is, for real symmetric matrices A, B, A ⪰ B (respectively A ≻ B), whenever, A − B is positive semidefinite (respectively positive definite.)

It is well-known that linear programming (LP for short) is an effective and powerful method for solving hard combinatorial optimization problems, for instance as a tool to obtain bounds in a branch and bound scheme or as a means to generate sharper constraints in a cutting plane approach. Our goal in this article is to show that semidefinite programming, which is even more powerful tool for solving such integer programming problems, can be solved effectively using interior point methods. Generalization from linear programming to semidefinite programming is quite natural. First, the duality theory extends easily from LP to SDP and results in generalization of Farkas’s lemma, strong duality theory and complimentary slackness theorem. These extensions are already known and we review them in section 1.

Second, the positive semidefinite cone is studied extensively and its properties are well–understood. This helps us extend, in a mechanical way, all known interior point linear programming techniques to algorithms which solve semidefinite programs; proofs of convergence and polynomial time complexity also extend in a similar fashion. In section 2. we develop a particular interior point method which is essentially a word by word extension of Ye’s projective method for linear programming. Once this particular extension is understood, it should become clear how to extend other LP interior point methods similarly. Proofs of convergence and polynomial time complexity extend in the same fashion. In fact we show that the number of iterations of our algorithm is proportional to O(\sqrt{n}) if X is an n × n symmetric matrix. (This extension is the main contribution of this paper.) Thus at least as far as interior point methods are concerned, semidefinite programming is not significantly more difficult than linear programming.

∗Research supported in part by NSF grant no. CDA-9211106, by the Air Force Office of Scientific Research grant AFOSR-87-0127, the National Science Foundation grant DCR-8420935 and the Minnesota Supercomputer Institute.

†International Computer Science Institute, Berkeley CA 94704. E-mail: alizadeh@icsi.berkeley.edu

Third, we study various “eigenvalue optimization” problems and show that they are indeed especial cases of SDP. This is the most common form that SDP arises in applications.

Finally, we study applications of SDP to combinatorial optimization problems. The connection is usually spectral properties of graphs. In particular, we study applications to maximum stable sets, graph partitioning, maximum cut, and the largest multipartite subgraph problem. Furthermore, we examine the new approach of Lovàsz and Schrijver [LS91] to the general 0-1 programming and demonstrate that our interior point approach may be used in most cases to solve Lovàsz –Schrijver semidefinite relaxations efficiently and perhaps practically. Lovász and Schrijver present a convincing argument that semidefinite programming relaxations of 0-1 optimization problems provide much deeper cuts than the analogous linear programming relaxations (see section 5.1 below.) Thus, SDP interior point methods may be employed to obtain considerably sharper bounds on the solutions of 0-1 programs.

In addition, some times it is possible to show that interior point methods may be used to derive sublinear time parallel algorithms. The first such result was obtained by Goldberg et. al in [GPST91] where they used interior point linear programming techniques to obtain the maximum matching in a bipartite graph. In [Ali91b] we used interior point methods for eigenvalue optimization to show that maximum cliques in perfect graphs may be computed in sublinear parallel time. This result is studied in section 5.2.

The SDP problem is a convex programming problem and, given exact input and à priori bound on the norm of the solution, may be solved in polynomial time to given accuracy by the ellipsoid method see [GLS88]. Y. Nesterov and A. Nemirovsky were the first ones to show that one may use interior point methods to optimize over the positive semidefinite cone [NN92]. In this important work the authors develop a general approach for using interior point methods for solving convex programming problems which is based on the concept of p–selfconcordant barrier functions. Nesterov and Nemirovsky show that for any convex set K that is endowed with a p–selfconcordant barrier function, there is an interior point algorithm which optimizes a linear function on K. Furthermore, every O(\sqrt{p}) iteration of this algorithm results in an interior point with half the distance to the optimal solution. As a special case, Nesterov and Nemirovsky show that linear programs with p inequality constraints, quadratic programs with p convex quadratic constraints and semidefinite programs over p × p matrices all admit p-selfconcordant barriers. Therefore, the authors extend the revolutionary result of Karmarkar to a rather general class of convex programs.

Our point of view regarding semidefinite programming is somewhat in the opposite direction of that of Nesterov and Nemirovsky. They obtain their complexity theorems by specializing their general results to SDP. In fact, they show that − ln det X is a n-self–concordant barrier for the cone of n × n positive semidefinite matrices and from here obtain an O(\sqrt{n}) iteration algorithm for SDP. We arrive at this conclusion by taking interior point methods for LP, a special SDP, and generalizing it to SDP. (Note that LP is a special case SDP: simply impose the linear constraints that X is a diagonal matrix and we get LP.) As a result we show that all known LP interior point methods may be extended, in a mechanical way, to SDP.

We should also mention two other works dealing with special instances of SDP. The first one is an early work of Cullum, Donath and Wolfe [CDW75], where a nonsmooth optimization technique is used to minimize sum of the first k eigenvalues of a linear combination of matrices. Their motivation was the graph partitioning problem described at the end of this paper. The second is the work of M. Overton [Ove92] where he develops efficient, though not polynomial time, algorithms based on the sequential linear programming and sequential quadratic programming techniques in nonlinear programming to solve the eigenvalue optimization problem. Overton’s algorithms, when viewed as solving SDP problems, are active set methods. In this sense they may be regarded in the same class as the simplex method for linear programming in that they both traverse the boundary of the feasible region to get to the solution. An empirical comparison of Overton’s approach to interior point methods has not been done, and seems to be a worthwhile project.

Notation.

We use lower case boldface letters to name column vectors, and upper case letters to name matrices. 1 and 0 denote vector of all ones and the zero vector, respectively. I and 0 denote the identity and zero matrices, respectively. ℜ^{\frac{n×n}{2}} is the set of symmetric n × n matrices. The i^{th} largest eigenvalue of a symmetric matrix X is λ_{i}(X); its i^{th} largest eigenvalue absolute-value-wise is λ^{i}(X). The Löwner partial order ⪰ and the dot product “•” were defined above; the symbol “≥” is used for component-wise comparison between two matrices or two vectors. Diag(x) denotes the diagonal matrix made up of vector x; diag(X) is the vector made up of diagonal elements of X. For matrices, ∥X∥ and ∥X∥_{2} are the Frobonius and the spectral norms of X, respectively (recall that in case of symmetric matrices ∥X∥_{2} equals the spectral radius ρ(X) = |λ^{1}(X)|). For vectors, ∥x∥ and ∥x∥_{∞} are the Euclidean and the maximum norms of x; also ∥x∥_{p} := (∑ |x_{i}^{p}|)^{1/p} is the p–norm of x. If A is a p × q matrix then, vec A is a pq column vector made up of columns of A stacked on each other. If v is a pq-vector then Mat_{pq} v is a p × q matrix whose i^{th} column is made up of the entries at (i − 1)p + 1 through ip in v; if p and q are clear from the context we will drop them from the subscript. A ⊗ B is the Kronecker product of matrices (see [Gra81].) If I and J are subsets of integers from 1 to p and from 1 to q, respectively, then A_{I,J} is the submatrix of A whose rows are taken from those rows of A indexed by I, and whose columns are indexed by J. A_{I.} and A_{.J} indicate rows indexed by I and columns indexed by J, respectively. In this paper the semidefinite programming problem refers to any optimization problem with any mixture of (symmetric) matrix and scalar valued variables which has a linear objective function and any combination of linear equality or (either component–wise ‘≥’ or Löwner ‘⪰’) inequality constraints.

G = (V, E) is a simple undirected graph without loops or multiple edges. A stable set S in G is a subset of vertices which are mutually nonadjacent. A clique K in G is a subset of vertices that are all mutually adjacent. A k-partite graph is one whose vertices can be partitioned into k subsets V_{j}, for j = 1, ⋯, k, where each V_{j} is an stable set. A clique covering of G is a collection K_{j}, j = 1, ⋯, k of vertices, where each K_{j} is a clique, and ∪_{j}K_{j} = V.

2 Duality theory.

A duality theory quite similar to that of linear programming may be constructed for the semidefinite programming problem. In this section we state the theory for the standard form SDP problem, the result for the more general form follows exactly as in linear programming. This theory has been developed in a more general context in many works before. It is easy to see that any pointed convex cone K ⊆ ℜ^{n} (that is K ∩ (−K) = {0}) induces a partial order ≥_{K} on ℜ^{n}: x ≥_{K} y iff x − y ∈ K. For instance, the nonnegative orthant and the positive semidefinite matrices induce the component–wise “≥” and the Löwner “⪰” partial orders, respectively. The duality theory in the usual linear programming can be extended to a generalized linear programming problem where “≥_{K}” (induced by an arbitrary pointed convex cone K) replaces “≥”. Duffin in [Duf56] was the first one to study such generalized duality theories. Later Ben-Israel, Charnes and Kortanek [BICK69], Borwein and Wolkowicz [BW81b, BW81a], and Wolkowicz [Wol81] among others developed more general formulations of the duality theory. Here we state this theory as applied to the SDP problem for the sake of completeness. For more general treatment consult the references mentioned above. Papers of Overton and Womersley [OW91a] and Fletcher [Fle85] treat duality theory for the more special SDP problems.

To derive a duality theory for SDP it is necessary to assume that C and A_{i} are symmetric. There is no loss of generality in this assumption. If C is not symmetric then, since C • X = C^{T} • X we may replace C by \frac{C+C^{T}}{2} . The same argument holds for A_{i}’s. These assumptions of symmetry will allow us to formulate the pair of primal and dual problems:

\underline{Primal} min  C • X s.t.  A_{i} • X = b_{i} for i = 1, ⋯, m X ⪰ 0  \underline{Dual} max b^{T} y s.t.  C − ∑_{i=1}^{m} y_{i}A_{i} ⪰ 0  (2.2)

Notice the similarity of primal and dual SDP pair to the corresponding linear programming pair. First we state a weak duality lemma.

Lemma 2.1. Let X be any feasible matrix for Primal and y any feasible vector for Dual. Then C •X ≥ b^{T} y.

Proof. (This result is well-known; we provide a proof to underscore the similarity to linear programming.) We have:

C • X − ∑_{i=1}^{m} b_{i}y_{i} = C • X − ∑_{i=1}^{m}(A_{i} • X)y_{i}

= (C − ∑_{i=1}^{m} y_{i}A_{i}) • X
≥ 0.

The last inequality is true because the inner product of two positive semidefinite matrices is nonnegative due to self–polarity of the positive semidefinite cone, see [HJ85].

We now state generalizations of Farkas’s lemma. Here is the most common form of the lemma as given in Schrijver’s text [Sch86], and as extended to the positive semidefinite cone:

Lemma 2.2. extend Farkas lemma: Let b ∈ ℜ^{m} and A ∈ ℜ^{m×n^{2}} be a matrix such that its rows A_{i.}^{T} = vecA_{i} where A_{i} is symmetric for i = 1, ⋯, m. Furthermore, let there be an m–vector y such that Mat(A^{T} y) ≻ 0. Then there exists a symmetric matrix X ⪰ 0, with A vecX = b if and only if y^{T} b ≥ 0 for all y for which Mat(A^{T} y) ⪰ 0.

(Recall that Mat(A^{T} y) is the same as ∑ y_{i}A_{i}.) Proof. For the only if part we have,

b^{T} y = (A vecX)^{T} y = Mat(A^{T} y) • X ≥ 0.

(the last inequality is due to self-polarity of the positive semidefinite cone.) To prove the if part, Suppose that the system A vecX = b, and X ⪰ 0 is infeasible. Define

K := {A vecX : X ⪰ 0}

Then b ∉ K. We claim that the set K is closed and prove it shortly. Since K is a closed cone there must exist a linear space that separates b and K strongly, i.e. there exists some vector y such that b^{T} y < 0 and (A vecX)^{T} y ≥ 0 for all X ⪰ 0, by [Roc70], Theorem 11.7, pp.100. But this means that X • Mat(A^{T} y) ≥ 0 for all X ⪰ 0, which is equivalent to Mat(A^{T} y) ⪰ 0, and therefore the if part of the theorem is proved.

Now we prove that the cone K is indeed closed. Let P be the positive semidefinite cone of n × n matrices and L := {Mat(A^{T} y) : y ∈ ℜ^{m}}. The condition that there is some y such that Mat(A^{T} y) ≻ 0 means that

L ∩ Int P ≠ ϕ

Thus any translate of the linear space L also intersects P and its interior. This is equivalent to saying that every symmetric n × n matrix can be written as the sum of two matrices, one of which is positive semidefinite and the other belongs to L. Therefore, ℜ^{\frac{n×n}{2}} = P + L. Taking polar we have

{0} = P ∩ L^{⊥}

Here L^{⊥} is the set {X : A vecX = 0}. Therefore, we have that X = 0 is the only solution of the system A vecX = 0, and X ⪰ 0. Therefore, [Roc70], Theorem 9.1, pp. 73 implies that K is closed.

A strong duality theorem similar to linear programming holds for SDP. In the following we say the primal problem in (2.2) is feasible if the set \{X ∈ ℜ^{\frac{n×n}{2}} : AvecX = b, and X ⪰ 0\} is nonempty, otherwise it is infeasible. Feasibility is defined similarly for the dual in (2.2). Recall that infimum over the empty set is by definition +∞ and supremum over the empty set is −∞. Furthermore, the primal (respectively dual) problem in (2.2) is unbounded if the infimum (respectively supremum) over the feasible set is −∞ (respectively +∞).

Theorem 2.1. Let

z_{1} := inf {C • X : AvecX = b, and X ⪰ 0}
z_{2} := sup {b^{T} y : C − Mat(A^{T} y) ⪰ 0}

Assume that there is an m–vector y such that A^{T} y ≻ 0. If either Primal or Dual is infeasible the other one is either unbounded or infeasible. If both Primal and Dual are feasible then z_{2} = z_{1}.

Proof. Suppose, z_{2} < z_{1}. Then, the system

\begin{matrix}C • X = z_{2} \\ AvecX = b \\ X ⪰ 0 \end{matrix}

is infeasible. Therefore, by the extended Farkas lemma 2.2, there exists a scalar y_{0} and m-vector y such that

y_{0}C + Mat(A^{T} y) ⪰ 0 and z_{2}y_{0} + b^{T} y < 0.  (2.3)

Now:

1. If y_{0} = 0, (2.3) is equivalent to

Mat(A^{T} y) ⪰ 0 and b^{T} y < 0,

which by extended Farkas lemma implies that AvecX = b and X ⪰ 0 is infeasible and thus z_{1} = ∞.

2. If y_{0} > 0, then dividing both sides of (2.3) by y_{0} we get

C − Mat(A^{T} (−y/y_{0})) ⪰ 0 and z_{2} − b^{T} (−y/y_{0}) < 0

which means that z_{2} is not the optimal solution of the dual problem.

3. If y_{0} < 0, then dividing both sides of (2.3) by −y_{0} we get

−C + Mat(A^{T} (−y/y_{0})) ⪰ 0 and − z_{2} + b^{T} (−y/y_{0}) < 0

In fact, since we have strict inequality we must have

−C + Mat(A^{T} (−y/y_{0})) ⪰ 0 and − z_{2} + b^{T} (−y/y_{0}) < −ϵ

for some ϵ > 0. But also, by optimality of z_{2} there must exist a y^{∗} such that

C − Mat(A^{T} y^{∗}) ⪰ 0 and z_{2} − b^{T} y^{∗} < ϵ

Adding the last two relations we get

Mat(A^{T} (−y/y_{0} − y^{∗}) ⪰ 0 and b^{T} (−y/y_{0} − y^{∗}) < 0

which again by extended Farkas lemma implies the primal problem is infeasible and z_{1} = ∞.

So the assumption z_{2} < z_{1} results in contradiction. Since by weak duality lemma we have that z_{2} ≤ z_{1} we conclude that z_{2} = z_{1}.

The condition that there should be a vector y such that A^{T} y ≻ 0, is known as the Slater’s condition. Without this stipulation the strong duality theorem 2.1 may not always hold. Wolkowicz and Borwein [BW81b, BW81a] and Wolkowicz [Wol81] develop a more general duality theory which does not require Slater’s condition. Since in this article we focus on interior point methods, and therefore have to assume that both the primal and the dual problems in (2.2) have nonempty interiors, imposing Slater condition does not restrict our presentation. Furthermore, in section 3.4 we show how to transform any pair of primal and dual problems to an equivalent pair which does have nonempty interior both for the primal and the dual (this is achieved by adding two scalar variables and two scalar constraints.)

It is also possible to derive a “complementary slackness” theorem. First we prove the following lemma:

Lemma 2.3. Let A and B be symmetric n×n matrices. If A ⪰ 0, B ⪰ 0, and A•B = 0, then AB = 0.

Proof. Let B = UΩU^{T} be the eigenvalue decomposition of B, with Ω = Diag(ω_{i}) and ω_{i} ≥ 0 for i = 1, ⋯, n. Set C := U^{T} AU, thus C ⪰ 0, and in particular, its diagonal elements C_{ii} ≥ 0. We only need to show that CΩ = 0. From A • B = 0 we have C • Ω = 0 and therefore, ∑_{i=1}^{n} C_{ii}ω_{i} = 0. Since all the summands are nonnegative, it follows that they are all zero. Thus we have:

1. If ω_{i} > 0 then C_{ii} = 0, and by C ⪰ 0, the entire row and column i is zero.

2. If C_{ii} > 0, then ω_{i} = 0.

Now suppose (CΩ)_{ij} ≠ 0 for some i, j. Then C_{ij}ω_{j} ≠ 0, which by 1 above we must have that the entire column j is zero, and so C_{ij} = 0, a contradiction.

Now the “complementary slackness” theorem follows easily :

Theorem 2.2. Let X^{∗} be a feasible matrix for the primal, and y^{∗} a feasible vector for the dual in (2.2). Define S^{∗} := C − Mat(A^{T} y^{∗}). Then X^{∗} and y^{∗} are primal and dual optimal, respectively, if and only if

X^{∗}S^{∗} = 0.  (2.4)

Proof. In the proof of the weak duality lemma 2.1 we established that the duality gap for any pair of primal and dual feasible solutions X and y is given by:

C • X − b^{T} y = X • S.  (2.5)

Now the “if” part of the theorem follows by the fact the duality gap will be zero. The “only if” part follows by applying lemma 2.3 to (2.5).

Notice that in contrast with linear programming, component-wise multiplication in the complementary slackness theorem is replaced by the ordinary matrix multiplication. The complementary slackness for SDP can be restated in the following way which makes it quite similar to the LP case:

Corollary 2.1. Let X^{∗} be a feasible matrix for the primal problem in (2.2) with eigenvalues λ_{1}, ⋯, λ_{n}; and S^{∗} := C − Mat(A^{T} y^{∗}) be feasible for the dual problem with eigenvalues ω_{1}, ⋯, ω_{n}. Then X^{∗} and S^{∗} are primal and dual optimal, respectively, if and only if they commute and there is a permutation π of eigenvalues of S^{∗} such that

λ_{i}ω_{π_{i}} = 0 for i = 1, ⋯, n.

Proof. X^{∗} and S^{∗} are optimal if and only if X^{∗}S^{∗} = 0. Thus, X^{∗} and S^{∗} commute with each other and therefore, they share a system of eigenvectors. Let columns of U be a joint system of orthonormal eigenvectors of X^{∗} and S^{∗}, i.e

X^{∗} = UDiag(λ_{1}, ⋯, λ_{n})U^{T} and S^{∗} = UDiag(ω_{π_{1}}, ⋯, ω_{π_{n}})U^{T}

for some permutation π. The corollary follows immediately by multiplying the right hand sides of these two identities.

Similar to linear programming, the complementary slackness theorem 2.2 may be used as a basis for primal-dual algorithms. Indeed in this paper, our interior point algorithm is a primal-dual method which maintains a primal feasible X_{k} and dual feasible S_{k} and each iteration moves X_{k}S_{k} closer to the zero matrix. The Frobonius norm ∥X_{k}S_{k}∥ is an indication of how close our current solution is to the optimum. We should mention here that the complementary slackness theorem, for a more restricted version of SDP has been reported in Grötschel, Lovász and Schrijver [GLS84], and in a paper by Shapiro [Sha85].

As in linear programming, semidefinite programs may arise in a variety of forms; the standard form (2.2) is just one type. Sometimes we may have positive semidefinite constraints imposed on linear combinations of matrices (as in the Dual, for example). Sometimes we may have ordinary inequalities “≥” on scalar variables in addition to Löwner relations. We may have several Löwner relations on several matrices. Finally, we may have some or all of these. Of course, as in linear programming, it is possible to convert all such problems to the standard form, usually by introducing new scalar and matrix variables and new constraints. However, it is more convenient to apply duality directly, as in linear programs in general form. One may come up with many variants of extended Farkas lemmas, all of which are extensions of similar lemmas for the component-wise inequalities, as given for example in Schrijver’s text [Sch86]. All of these extensions assume some variant of Slater’s condition.

Lemma 2.4. Let A ∈ ℜ^{n^{2}×m} be a matrix whose columns are linearly independent and are of the form vecA_{i} for symmetric A_{i}, and B ∈ ℜ^{\frac{n×n}{2}} . Assume that there exists some symmetric matrix Y ≻ 0 such that (vecY)^{T} A = 0. Then Mat(Ax) ⪯ B has a solution in x if and only if B • Y ≥ 0 for all Y ⪰ 0 for which (vecY)^{T} A = 0.

Proof. To show the “if” part notice that

B • Y = B • Y − x^{T} (A^{T} vecY) = (B − MatAx) • Y ≥ 0

To show the “only if” part suppose the system Mat(Ax) ⪯ B does not have any solution. Then B ∉ K where

K := {Mat(Ax) + S : S ⪰ 0 and x ∈ ℜ^{m}}

We show that K is a closed cone; but notice that closedness of K implies the “only if” part of the lemma because by the separating hyperplane theorem for cones, ([Roc70], Theorem 11.7, pp.100), there must exist a matrix Y such that Y • B < 0 and Y • (Mat(Ax) + S) ≥ 0 for all S ⪰ 0 and all x ∈ ℜ^{m}. In particular, for x = 0 We have that Y • S ≥ 0 for all S ⪰ 0 which implies Y ⪰ 0; and for S = 0 we have that (vecY)^{T} Ax ≥ 0 for all x. Replacing x by −x we find that (vecY)^{T} Ax = 0 for all x which implies (vecY)^{T} A = 0.

To show K is closed notice that if P is the positive semidefinte cone, and L := \{Y ∈ ℜ^{\frac{n×n}{2}} : (vecY)^{T} A = 0} then we have that

L ∩ IntP ≠ ϕ

Thus any translate of L also intersects IntP, which implies ℜ^{\frac{n×n}{2}} = L + P. Taking polar we get

{0} = L^{⊥} ∩ P

which means S = 0 is the only matrix satisfying MatAx − S = 0, S ⪰ 0, which in turn implies that x = 0 since rows of A are linearly independent. Therefore, by [Roc70] theorem 9.1 pp.73, implies that K is closed.

Several other variants of extended Farkas lemma can be formulated and proved in a similar vain. We mention a few more:

Lemma 2.5. Let A ∈ ℜ^{m×n} and B ∈ ℜ^{m×m}. Suppose there exist some matrix Y such that A^{T} Y A ≻ 0. Then the system AXA^{T} = B and X ⪰ 0 has a solution iff for all symmetric matrices Y, A^{T} Y A ⪰ 0 implies that B • Y ≥ 0.

Lemma 2.6. Let A ∈ ℜ^{m×n} and B ∈ ℜ^{m×m}. Suppose there exist some matrix Y such that A^{T} Y A = 0 and Y ≻ 0. Then the system AXA^{T} ⪯ B has a solution iff for all symmetric matrices Y ⪰ 0 and A^{T} Y A = 0 implies that B • Y ≥ 0.

MIN  MAX
matrix or scalar, ≥ 0  ⟷ C  matrix or scalar, ≤
V  matrix or scalar, ≤ 0  ⟷ O  matrix or scalar, ≥
A  matrix, ⪰ 0  ⟷ N  matrix ⪯
R  matrix, ⪯ 0  ⟷ S  matrix, ⪰
matrix or scalar, unrestricted ⟷ T  matrix or scalar, =
C  matrix or scalar, ≥  ⟷  matrix or scalar, ≥ 0
O  matrix or scalar, ≤  ⟷ V  matrix or scalar, ≤ 0
N  matrix, ⪰  ⟷ A  matrix, ⪰ 0
S  matrix, ⪯  ⟷ R  matrix, ⪯ 0
T  matrix or scalar, =  ⟷  matrix or scalar, unrestricted

Figure 1: Duality rules for semidefinite programming.

Lemma 2.7. Let A ∈ ℜ^{m×n} and B ∈ ℜ^{m×m}. Suppose there exist some matrix Y such that A^{T} Y A ≻ 0. Then the system AXA^{T} ⪰ B and X ⪰ 0 has a solution iff for all symmetric matrices Y ⪰ 0 and A^{T} Y A ⪰ 0 implies that B • Y ≥ 0.

It is easy to show that the rules of obtaining the dual are a straightforward extension of these rules for the linear programming problem. The main extension is that constraints that involve semidefinite relations on matrix valued expressions give rise to matrix valued dual variables with semidefinite constraints on them. These rules are summarized in the table in figure 1; this table is an extension of a similar table in the text of Bazaraa, Jarvis and Sherali [BJS90].

3 An interior point Algorithm.

In this section we develop a potential reduction method for solving the primal problem so that, within O(\sqrt{n}| log ϵ|) iterations, we get an approximate solution with at least ϵ relative accuracy, if ϵ is sufficiently small. Our development closely follows Ye’s projective technique for linear programming [Ye90]. Ye’s complexity analysis is also extended to semidefinite programs.

3.1 Potential functions and projective transformations. First, recall that the interior of the cone of positive semidefinite matrices is the set of positive-definite matrices; therefore, all interior points are nonsingular. The boundary of the cone consists of singular semidefinite matrices and so, some of the eigenvalues of the boundary matrices are zero. In particular, optimal solutions of the primal problem in 2.2 are singular.

We assume that the primal and the dual problems have non-empty interiors, with given initial primal and dual points and with finite optimal solutions. Later, in section 3.4, we show how to transform any primal-dual pair to an equivalent one where an initial interior primal–dual solution is available. Let q > 0, \underline{z} be given constant known to be a lower bound on the optimal value z^{∗} of the primal problem in (2.2). Let X be an interior primal feasible matrix, y an interior dual feasible vector, and S := C − ∑_{i=1}^{m} y_{i}A_{i}; thus, X ≻ 0 and S ≻ 0. Define the primal potential function:

ϕ(X, \underline{z}) = q ln(C • X − \underline{z}) − ln det X,  (3.6)

and the primal-dual potential function:

ψ(X, S) = q ln(X • S) − ln det(XS).  (3.7)

For motivation, one may think of semidefinite constraints X ⪰ 0, expressed as λ_{i}(X) ≥ 0 for i = 1, ⋯, n. When standard logarithmic barrier is applied to these constraints we get: ∑_{i=1}^{n} ln λ_{i}(X) = ln det X.

The strategy of the algorithm is to generate a sequence of interior primal feasible matrices X_{k}, and a sequence of interior dual vector–matrix pairs (y_{k}, S_{k}) such that the sequence of ψ(X_{k}, S_{k}) decreases at least like an arithmetic progression. With an appropriate choice of q, this would imply that the “duality gap” C • X_{k} − b^{T} y_{k} decreases at least like a geometric progression with k; in particular it becomes a constant fraction of the original gap after O(\sqrt{n}) iterations.

Before describing the algorithm we state the following lemma which is a direct generalization of a similar lemma that appears in the analysis of most interior point linear programming methods. (Recall that ρ(X) is the spectral radius of matrix X, which equals its largest eigenvalue when X is positive semidefinite.)

Lemma 3.1. Let X be a symmetric n × n matrix. If 0 ≺ X ≺ I, then

ln det X ≥ trace X − n − \frac{trace (X − I)^{2}}{2[1 − ρ(X − I)]}

Proof. In most interior-point linear programming algorithms it is shown that if ∥x − 1∥_{∞} < 1 then

∑_{j=1}^{n} ln x_{j} ≥ (1^{T} x − n) − \frac{∥x − 1∥^{2}}{2(1 − ∥x − 1∥_{∞})}

which is easily proved by expanding ln x, (see for example, Karmarkar [Kar84] or Ye [Ye91].) Now to prove the lemma simply replace λ_{j}(X) for x_{j}.

We use a projective transformation to bring the current iterate to the center, except that the center here is the identity matrix (in contrast with linear programming in which the center is 1). An important point is that the transformation should map the set of symmetric matrices to itself. This is needed so that the transformed problem remains a meaningful SDP problem. To find a symmetry preserving projective transformation that maps X_{0} to the identity matrix I, let L_{0} be any matrix such that L_{0}L_{0}^{T} = X_{0}. There are infinitely many choices for L_{0}. For instance, it could be a Cholesky factor of X_{0}, or it could be its square root, X_{0}^{1/2}. We shall see shortly that it does not matter how we select L_{0} as it will not affect the algorithm’s behavior and performance. Fix integer r. Define T : ℜ^{\frac{n×n}{2}} → ℜ^{\frac{n×n}{2}} × ℜ^{r}, such that T (X) = (\overline{X}, \overline{x}). Then:

\overline{X} := \frac{(n + r)L_{0}^{−1}XL_{0}^{−T}}{r + X_{0}^{−1} • X} and \overline{x} := (\frac{n + r}{r + X_{0}^{−1} • X} ) 1  (3.8)

Also, the inverse transformations is given by:

X = T^{−1}(\overline{X}, \overline{x}) := \frac{L_{0}\overline{X}L_{0}^{T}}{∑ \overline{x}_{j}/r}  (3.9)

Under T , Primal is transformed into

min \overline{C} • \overline{X} + \overline{c}(\underline{z})^{T} \overline{x}
s.t. \overline{A}vec \overline{X} + \overline{A}\overline{x} = 0
trace \overline{X} + 1^{T} \overline{x} = n + r
\overline{X} ⪰ 0
\overline{x} ≥ 0

(3.10)

where

\overline{C} := L_{0}^{T} CL_{0}  (3.11)
\overline{c}(\underline{z}) := − (\underline{z}/r) 1  (3.12)
\overline{A}_{i} := L_{0}^{T} A_{i}L_{0}  (3.13)
\overline{A} := A(L_{0} ⊗ L_{0})  (3.14)
\overline{A} := (−1/r) b1^{T}  (3.15)

Note that \overline{A} is an m × r rank one matrix. The transformed problem may be viewed as a mixed linear and semidefinite program. We may define the following primal potential function for the transformed problem:

\overline{ϕ}(\overline{X}, \overline{x}, \underline{z}) := q ln [\overline{C} • \overline{X} + \overline{c}(\underline{z})^{T} \overline{x}] − ln det \overline{X} − ∑_{j=1}^{r} ln \overline{x}_{j}  (3.16)

The following invariant property holds for the potential functions under projective transformations:

Lemma 3.2. If \overline{x}_{1} = ⋯ = \overline{x}_{r}, and q = n + r then

ϕ(X, \underline{z}) − ϕ(X_{0}, \underline{z}) = \overline{ϕ}(\overline{X}, \overline{x}, \underline{z}) − \overline{ϕ}(I, 1, \underline{z}).  (3.17)

Also the following result is easily proved by expanding \overline{ϕ} and applying lemma 3.1; later we use it to prove the reduction in the potential function.

Corollary 3.1. For q = n + r we have

\overline{ϕ}(\overline{X}, \overline{x}, \underline{z}) − \overline{ϕ}(I, 1, \underline{z}) ≤ (n + r) ln (\frac{\overline{C} • \overline{X} + \overline{c}(\underline{z})^{T} \overline{x}}{trace \overline{C} + \overline{c}(\underline{z})^{T} 1}) + \frac{∥\overline{X} − I∥^{2} + ∥\overline{x} − 1∥^{2}}{2(1 − ∥\overline{X} − I∥ + ∥\overline{x} − 1∥)}  (3.18)

3.2 A potential reduction algorithm. Similar to linear programming, in (3.10) we replace the inequality constraints \overline{X} ⪰ 0 and \overline{x} ≥ 0 by an inscribed “ball” constraint, except that for the SDP problem the ball is centered at (I, 1). Therefore, (3.10) is replaced by the “ball optimization” problem:

min \overline{C} • \overline{X} + \overline{c}(\underline{z})^{T} \overline{x}
s.t. \overline{A}vec \overline{X} + \overline{A}\overline{x} = 0
trace \overline{X} + 1^{T} \overline{x} = n + r
∥X − I∥^{2} + ∥x − 1∥^{2} ≤ β^{2} < 1

(3.19)

where β is a fixed constant between 0 and 1 to be determined shortly. Once we solve this problem and map the result back to the original space, we get a point that serves as a candidate for the next iterate. The solution of (3.19) is given by

( \begin{matrix}vec & \overline{X} & 1 \\ \overline{x} & 1 \end{matrix} ) := ( \begin{matrix}vecI \\ 1 \end{matrix} ) − β \frac{P (\underline{z})}{∥P (\underline{z})∥},  (3.20)

and thus the candidate for the new primal iterate is given by:

X_{1} := T^{−1}(\overline{X}_{1}, \overline{x}_{1})  (3.21)

where

P (\underline{z}) := P_{A′} ( \begin{matrix}vec\overline{C}(\underline{z}) \\ \overline{c}(\underline{z}) \end{matrix} ) , A′ := ( \begin{matrix}\overline{A} & \overline{A} \\ (vec I)^{T} & 1^{T} \end{matrix} ) ,  (3.22)

and P_{A′}(u) is the projection of the (n^{2} +r)-vector u to the null space of A′. After expansion the projection P_{A′} in (3.22) becomes:

P (\underline{z}) = (I − \frac{[vecI|1][vecI|1]^{T}}{n + r} ) (I − [\overline{A}|\overline{A}]^{T} ([\overline{A}|\overline{A}][\overline{A}|\overline{A}]^{T})^{−1}[\overline{A}|\overline{A}]) ( \begin{matrix}vec\overline{C}(\underline{z}) \\ \overline{c}(\underline{z}) \end{matrix} )  (3.23)

Define:

y(\underline{z}) := ([\overline{A}|\overline{A}][\overline{A}|\overline{A}]^{T} )^{−1} [\overline{A}|\overline{A}] ( \begin{matrix}vec\overline{C}(\underline{z}) \\ \overline{c}(\underline{z}) \end{matrix} )
= (A(X_{0} ⊗ X_{0})A^{T} + (1/r) bb^{T} )^{−1} [Avec(X_{0}CX_{0}) + (\underline{z}/r) b]

(3.24)

and

S(\underline{z}) := C − Mat(A^{T} y(\underline{z}))  (3.25)

S(\underline{z}) and y(\underline{z}) serve as candidates for the new dual iterates. In terms of these quantities P (\underline{z}) may be written as:

P (\underline{z}) = ( \begin{matrix}vec(L_{0}^{T} S(\underline{z})L_{0}) \\ \frac{b^{T} y(\underline{z})−\underline{z}}{r} 1 \end{matrix} ) − \frac{C • X_{0} − \underline{z}}{n + r} ( \begin{matrix}vecI \\ 1 \end{matrix} )  (3.26)

Observe that X_{1}, S(\underline{z}) and y(\underline{z}) are all independent of L_{0}; in fact in actual computation we do not need to have L_{0} explicitly.

Now we show that either the primal candidate X_{1}, or the dual candidates S(\underline{z}) and y(\underline{z}) reduce the value of the primal-dual potential function ψ by a constant amount. First observe that P_{A′} is a projector, that is P_{A′}^{2} = P′_{A}. Therefore, from (3.20) we get:

\overline{C}(\underline{z}) • (\overline{X}_{1} − I) + \overline{c}(\underline{z})^{T} (\overline{x} − 1) = −β∥P (\underline{z})∥

Hence, noting that ln(1 + x) ≤ x, for nonnegative x, corollary 3.1 implies:

Corollary 3.2. Let q = n + r and \overline{X}_{1} and \overline{x} be as in (3.20). Then

\overline{ϕ}(\overline{X}, \overline{x}, \underline{z}) − \overline{ϕ}(I, 1, \underline{z}) ≤ −(n + r)β \frac{∥P (\underline{z})∥}{c(\underline{z})^{T} 1 + trace \overline{C}} + \frac{β^{2}}{2(1 − β)}

Let Δ_{0} be the size of the duality gap in the current iterate, that is

Δ_{0} := C • X_{0} − \underline{z}

and let

Δ_{1} := S(\underline{z}) • X_{0} = C • X_{0} − b^{T} y(\underline{z})

Thus Δ_{1} should be interpreted as the value of duality gap if we choose y(\underline{z}) as our new dual iterate. Before deriving the amount of reduction in the potential function we prove the following lemma:

Lemma 3.3. If there is some real number α with 0 < α < 1, such that

∥P (\underline{z})∥ ≤ α\frac{Δ_{0}}{n + r}

then S(\underline{z}) ≻ 0, and b^{T} y(\underline{z}) > \underline{z}. Furthermore,

∥L_{0}^{T} S(\underline{z})L_{0} − \frac{Δ_{1}}{n} I∥ ≤ \frac{Δ_{1}}{n} α\sqrt{\frac{n + n^{2}/r}{n + n^{2}/r − α^{2}}} ,  (3.27)

and

|\frac{n + r}{n} \frac{Δ_{1}}{Δ_{0}} − 1| ≤ \frac{α}{\sqrt{n + n^{2}/r}}  (3.28)

Proof. Suppose S(\underline{z}) ⊁ 0. Then L_{0}^{T} S(\underline{z})L_{0} is not positive definite and so some of its eigenvalues are less than or equal to 0. Thus, from (3.26) we have

∥P (\underline{z})∥ ≥ ρ (\frac{Δ_{0}}{n + r} I − L_{0}^{T} S(\underline{z})L_{0}) ≥ \frac{Δ_{0}}{n + r} ,

a contradiction. Also, If b^{T} y(\underline{z}) ≤ \underline{z} then from (3.26) we have

∥P (\underline{z})∥ ≥ \frac{Δ_{0}}{n + r} − \frac{b^{T} y(\underline{z}) − \underline{z}}{r} ≥ \frac{Δ_{0}}{n + r} ,

which is again a contradiction. Now from (3.26) we have

P (\underline{z}) = ( \begin{matrix}[vec(L_{0}^{T} S(\underline{z})L_{0}) − \frac{Δ_{1}}{n} I] − [\frac{Δ_{0}}{n+r} − \frac{Δ_{1}}{n} ] vecI \\ [\frac{Δ_{0}−Δ_{1}}{r} − \frac{Δ_{0}}{n+r} ] 1 \end{matrix} )

Since I • [(L_{0}^{T} S(\underline{z})L_{0}) − (Δ_{1}/n)I] = 0, we have

∥P (\underline{z})∥^{2} = ∥L_{0}^{T} S(\underline{z})L_{0} − \frac{Δ_{1}}{n} I∥^{2} + n (\frac{Δ_{0}}{n + r} − \frac{Δ_{1}}{n} )^{2} + r (\frac{Δ_{0} − Δ_{1}}{r} − \frac{Δ_{0}}{n + r} )^{2}

= ∥L_{0}^{T} S(\underline{z})L_{0} − \frac{Δ_{1}}{n} I∥^{2} + (n + \frac{n^{2}}{r} ) (\frac{Δ_{1}}{n} − \frac{Δ_{0}}{n + r} )^{2} .  (3.29)

If (3.27) is false we must have from (3.29) we have

∥P (\underline{z})∥^{2} = (\frac{Δ_{1}}{n} )^{2} α^{2} \frac{n + n^{2}/r}{n + n^{2}/r − α^{2}} + (n + \frac{n^{2}}{r} ) (\frac{Δ_{1}}{n} − \frac{Δ_{0}}{n + r} )^{2}

≥ α^{2} (\frac{Δ_{0}}{n + r} )^{2}  (3.30)

(The last inequality is proved by taking the right hand side of the equality in (3.30) as a quadratic function in Δ_{1}/n and minimizing it.) But (3.30) contradicts the assumption of the lemma, so (3.27) must be true. Finally (3.28) follows from (3.29).

Now we may prove the potential reduction theorem.

Theorem 3.1. Let X_{0} be any interior feasible point for the primal problem (1.1) and y_{0} interior feasible for the dual. Let also, r := ⌈\sqrt{n} ⌉ and q := n + r, S_{0} := C − ∑_{i=1}^{m} y_{i}A_{i}, \underline{z}_{0} := b^{T} y_{0}, X_{1} := T^{−1}(\overline{X}_{1}, \overline{x}_{1}), as in (3.9), y_{1} := y(\underline{z}_{0}), and S_{1} := S(\underline{z}_{0}). Then there exist an absolute constant δ such that either

ψ(X_{1}, S_{0}) ≤ ψ(X_{0}, S_{0}) − δ,

or

ψ(X_{0}, S_{1}) ≤ ψ(X_{0}, S_{0}) − δ,

Furthermore, if we set α = 0.55 and β := 0.3, then δ > 0.1.

Proof. If

∥P (\underline{z})∥ ≥ α\frac{Δ_{0}}{n + r}

then

ψ(X_{1}, S_{0}) − ψ(X_{0}, S_{0}) = ϕ(X_{1}, \underline{z}_{0}) − ϕ(X_{0}, \underline{z}_{0})
= \overline{ϕ}(\overline{X}_{1}, \overline{x}_{1}, \underline{z}_{0}) − \overline{ϕ}(I, 1, \underline{z}_{0})
≤ −βα + \frac{β^{2}}{2(1 − β)}

(the last inequality is true by corollary 3.2). Otherwise, the conditions of lemma 3.3 are satisfied. Also

applying lemma 3.1 to (n/Δ_{1})L_{0}^{T} S_{1}L_{0}, and setting γ := α\sqrt{\frac{n+n^{2}/r}{n+n^{2}/r−α^{2}}} we have:

n ln X_{0} • S_{1} − ln det X_{0}S_{1} = n ln (\frac{nX_{0} • S_{1}}{Δ_{1}} ) − ln det \frac{nX_{0}S_{1}}{Δ_{1}}
= n ln n − ln det \frac{nX_{0}S_{1}}{Δ_{1}}
≤ n ln n + \frac{∥nL_{0}^{T} S_{1}L_{0}/Δ_{1} − I∥^{2}}{2(1 − ∥nL_{0}^{T} S_{1}L_{0}/Δ_{1} − I∥)}
≤ n ln X_{0} • S_{0} − ln det X_{0}S_{0} + \frac{γ^{2}}{2(1 − γ)},  (3.31)

Algorithm SDP: Input:

An n × n matrix X_{0}, interior feasible for the primal problem in (2.2); an m-vector y_{0} interior feasible for the dual problem; a constant ϵ.

Output:

A primal feasible solution X and dual feasible solution y such that C • X − b^{T} y < ϵ.

Method:

1) Set \underline{z}_{0} = b^{T} y_{0}.

2) Set S_{k} := C − Mat(A^{T} y_{0}).
3) Set k = 0.
4) While C • X_{k} − b^{T} y_{k} ≥ ϵ do

begin

Compute S(\underline{z}_{k}) from (3.25) and P (\underline{z}_{k}) from (3.22).
If ∥P (\underline{z}_{k})∥ ≥ α(C • X_{k} − \underline{z}_{k})/(n + r) then

a)Find β^{∗} := argmin_{0≤β≤1}ψ(X_{k} − βL_{k}P (\underline{z}_{k})L_{k}^{T} , S_{k}),

using a line search procedure.

b) Set (\overline{X}_{k+1}, \overline{x}_{k+1}) = (I, 1) − β^{∗}P (\underline{z}), and set X_{k+1} := T^{−1}(\overline{X}_{k+1}, \overline{x}_{k+1}).

c) Set S_{k+1} := S_{k}, and \underline{z}_{k+1} := \underline{z}_{k}.

Else

d) Find \underline{z}^{∗} := argmin_{\underline{z}≤\underline{z}_{k}}ψ(X_{k}, S(\underline{z})) by a line search.

e) Set S_{k+1} = S(\underline{z}^{∗}).

f) Set X_{k+1} = X_{k}, and \underline{z}_{k+1} = b^{T} y(\underline{z}^{∗}).

Set k = k + 1. end.

Figure 2: A projective potential reduction algorithm.

where the last relation results from applying the arithmetic-geometric mean inequality to the eigenvalues of X_{0}S_{0} (which are all real.) By (3.28) of lemma 3.3 we have

Δ_{1} < (1 − \frac{r}{n + r} − \frac{n}{n + r} \frac{α}{\sqrt{n + n^{2}/r}} ) Δ_{0}.

Thus,

r ln \frac{X_{0} • S(\underline{z})}{X_{0} • S_{0}} = r ln \frac{Δ_{1}}{Δ_{0}} ≤ \frac{r^{2}}{n + r} (−1 + \frac{α}{\sqrt{r + r^{2}/n}}) .  (3.32)

Adding (3.31) and (3.32) we get

ϕ(X_{0}, S_{1}) − ϕ(X_{0}, S_{0}) ≤ \frac{r^{2}}{n + r} (−1 + \frac{α}{\sqrt{r + r^{2}/n}}) + \frac{γ^{2}}{2(1 − γ)}.  (3.33)

It is easily verified that choice of α = 0.55, β = 0.3 and δ = 0.1 is consistent with all the conditions of the theorem.

Based on this result we present the projective version of the algorithm displayed in figure 2.

3.3 Potential reduction and polynomial-time solvability. Now we show that starting from any pair of interior primal and dual feasible points, and a tolerance ϵ, we get a pair whose duality gap is less than ϵ by running the algorithm in figure 2 a number of times that depends polynomiality on ln ϵ, n and the error in the initial pair.

Theorem 3.2. Let X_{0}, y_{0} and S_{0} := C −Mat(A^{T} y_{0}) be given initial interior points for the primal and dual semidefinite programming problems in (2.2). Let also that r = ⌈\sqrt{n} ⌉ and q = n + r in the primal-dual potential function ψ, and assume that ψ(X_{0}, S_{0}) ≤ O(\sqrt{n}E) for some constant E. If an algorithm generates a sequence of interior primal and dual points X_{j}, y_{j} (and thus S_{j}) such that ψ(X_{j}, S_{j}) ≥ ψ(X_{j+1}, S_{j+1})+δ for some fixed number δ then, after k = O(\sqrt{n}| log ϵ|) iterations, for primal and dual solutions X_{k}, y_{k} and S_{k} we have

C • X_{k} − b^{T} y_{k} < 2^{E}ϵ.

Proof. Each iteration reduces the potential function by at least δ. Thus, if ψ(X_{0}, S_{0}) < O(\sqrt{n}E) then after O(\sqrt{n}| log ϵ|) iterations we have:

ψ(X_{k}, S_{k}) < (\sqrt{n}E − \sqrt{n}| log ϵ|)
= (\sqrt{n}[log 2^{E} − | log ϵ|]
≤ \sqrt{n}| log(2^{E}ϵ)|.

Therefore,

\sqrt{n} ln X_{k} • S_{k} < −n ln X_{k} • S_{k} + ln det X_{k}S_{k} + \sqrt{n}| log(2^{E}ϵ)|
< −n ln n + \sqrt{n}| ln(2^{E}ϵ)|

The last inequality comes from applying the arithmetic-geometric inequality to the eigenvalues of X_{k}S_{k}, which are real, as both matrices are positive definite. Thus, ln X_{k} • S_{k} < | log(2^{E}ϵ)|, and since X_{k} • S_{k} = C • X_{k} − b^{T} y_{k}, the theorem follows.

What we have just shown essentially says that if we start our potential reduction algorithm at a pair of primal and dual points where the initial error is small enough so that the value of the potential function is less than O(\sqrt{n}E), then after O(\sqrt{n}(E + | log ϵ|)) iterations we will have a solution with duality gap less than ϵ. Therefore for all ϵ < 2^{−E} the term | log ϵ| dominates E and so the number of iterations is bounded by O(\sqrt{n}| log ϵ|). Also observe that this proof solely depends on the reduction of the potential function ψ. We have to guarantee a reduction of at least δ in each iteration; but larger reductions may speed up the algorithm without affecting its worst case complexity. Therefore, in steps 4a and 4d of the algorithm in figure 2 we allow a line search to find a step length β^{∗} which maximizes the reduction in the potential function.

3.4 Feasibility, boundedness and polynomial–time computability. In this section we study some sufficient conditions for polynomial time computability of the general SDP problem. The situation is somewhat different from the case of linear programming. Let us assume that all entries in the primal and dual problems (2.2) are integers. First, in contrast with linear programming, the optimal solution of (2.2) is not necessarily a rational number. Therefore we need to specify an error tolerance, ϵ, and ask for a pair of primal and dual solutions X and S such that the duality gap X • S ≤ ϵ^{1}. If ϵ is also a rational number, define L, the size of the SDP problem, as the number of bits in the binary representation of ϵ and entries of C, A, and b. One might expect that the interior point method developed in the previous sections leads to an algorithm which runs in time polynomial in m, n and L. However, this may not be true as the solution itself may be exponentially large. To see this consider the optimization problem:

min{x_{n} : x_{1} = 2, and x_{i} ≥ x_{i−1}^{2} for i = 2, ⋯, n}  (3.34)

1Since X, S and y are solution of the algebraic system of equations: XS = 0, AvecX = b and A^{T} y + S = C, there are algebraic solutions among all optimal solutions of an SDP problem with integral input.

Clearly, x_{n} = 2^{2^{n}} is the solution of this problem which requires exponential number of bits. Now (3.34) can be written as the following semidefinite program:

min x_{n}
s.t. x_{1} = 2
( \begin{matrix}x_{i} & x_{i−1} \\ x_{i−1} & 1 \end{matrix} ) ⪰ 0 for i = 2, ⋯, n

This SDP can be easily turned into a standard form SDP whose input size (taking ϵ = 1, say) is polynomial in n and whose output requires more than exponential number of bits. So no algorithm can solve it in polynomial time^{2}

In many cases, including all of the combinatorial optimization problems described below, one may be able to put an à priori bound on the norms of the optimal solutions. For instance in special cases we may be able to prove that ∥X∥ = O(2^{nm}), and ∥y∥ = O(2^{nm}). In such cases we can show that the interior point algorithm developed earlier can produce primal and dual solutions in polynomial time whose duality gap is smaller than ϵ. Let L′ be the number of bits in the binary expansion of some rational number known to be a bound on ∥X^{∗}∥ and ∥y^{∗}∥. Then, similar to linear programming, one can always transform the pair of primal and dual problems (2.2) to another pair for which initial interior feasible points are readily available. We extend the construction suggested by Kojima, Mizuno and Yoshise in [KMY89] which in turn is based on Megiddo’s in [Meg89].

Consider the following pair of primal and dual problems:

min C • X + Mx_{1}
s.t. Avec(X) + [b − Avec(X_{0})]x_{1} = b
[Mat(A^{T} y_{0}) + S_{0} − C] • X + x_{2} = N
X ⪰ 0
x_{1}, x_{2} ≥ 0

(3.35)

and,

max b^{T} y − Ny_{1}
s.t. Mat(A^{T} y) + S + [C − Mat(A^{T} y_{0}) − S_{0}]y_{1} = C
[b − Avec(X_{0})]^{T} y + y_{2} = M
S ⪰ 0
y_{1}, y_{2} ≥ 0

(3.36)

where X_{0} and S_{0} are arbitrary positive definite n × n matrices, y_{0} an arbitrary m-vector, and M and N are large enough positive numbers to ensure that y_{2} ≥ 0 and x_{2} ≥ 0. Clearly X := X_{0}, x_{1} := 1 and x_{2} := N −(Mat(A^{T} y_{0})+S_{0} −C)•X_{0} are interior feasible for the primal (3.35) (with large enough N); and S := S_{0}, y := y_{0}, y_{1} := 1, and y_{2} := M − (b − Avec(X_{0}))^{T} y_{0} are interior feasible for the dual problem in (3.36) (for large enough M). By choosing X_{0} = S_{0} = I, x_{1} = y_{1} = 1, y_{0} = 0, it suffices to choose M and N such that

N > max (n − ∑_{i} C_{ii}, ∑_{i} X_{ii}^{∗} − C • X^{∗}) ,

M > max (0, b^{T} y^{∗} − ∑_{i} y_{i}^{∗}trace (A_{i})) .

For instance we may set N = M = 2^{L+L′}. It is easy to see that if the optimal value of x_{1} is not zero, then the original primal is unbounded (the proof is exactly like the one given in Kojima et al. in [KMY89]). Similarly if the optimal value of y_{1} is not zero, then the original dual is unbounded. Otherwise, the optimal X^{∗} and y^{∗} are also optimal for the original primal and dual problems, respectively. Furthermore, It is

2I am indebted to Joshi Ramana for bringing to my attention an error in [Ali91a, Ali92] where I had claimed that the norm of the solution to any SDP problem is bounded by 2^{L}. Joshi essentially provided this counter example.

LP  SDP
unknown vector: x  unknown symmetric matrix: X
inequality constraints: ≥  Löwner constraints: ⪰
dual variable: y  dual variable: y
dual slack vector: s  dual slack symmetric matrix: S
1  I
linear scaling:  linear scaling:
x → (x_{i}/(x_{0})_{i})_{i=1}^{n} = [Diag(x_{0})]^{−1}x X → L_{0}^{−1}XL_{0}^{−T} = Mat[(L_{0}^{−1} ⊗ L_{0}^{−1})vec(X)]
projective scaling:  projective scaling:
x → \frac{c_{1}[Diag(x_{0})]^{−1}x}{c_{2}+1^{T} [Diag(x_{0})]^{−1}x}  X → \frac{c_{1}L_{0}^{−1}XL_{0}^{−T}}{c_{2}+trace L_{0}^{−1}XL_{0}^{−T}}
barrier function:  barrier function:
∑ ln x_{i}  ln det X
norms:  norms:
∥x∥  ∥X∥
∥x∥_{∞}  ∥X∥_{2}
∥x∥_{p}  (∑ |λ_{i}(X)|^{p})^{1/p}

Figure 3: Correspondence between linear programming and semidefinite programming

easily verified that the value of the primal-dual potential function ψ at the initial point is bounded by O(\sqrt{n}(L + L′)). So, for the general SDP problem, any algorithm that reduces the primal-dual potential function ψ by a constant amount may find, in O(\sqrt{n} max(L, L′, | log ϵ|)) iterations, a pair of primal and dual feasible solutions whose duality gap is less than ϵ; if ϵ < 2^{−L−L′}, then the number of iterations is bounded by O(\sqrt{n}| ln ϵ|).

3.5 A correspondence between proofs in linear and semidefinite programming. The remarkable similarity between the algorithm presented here and Ye’s LP algorithm in [Ye90] suggests that other LP interior point methods may also be extended to SDP problems. All proofs of convergence and polynomialtime complexity may be extended as well. The correspondence is summarized in figure 2. Given any interior point algorithm for linear programming we may construct, in a mechanical way, an algorithm for the SDP problem by replacing any references to the entries under the LP column, with the corresponding entry under the SDP column. Proofs of convergence or polynomial time complexity may also be extended mechanically in the same manner. We have already verified this claim on the approaches of Gonzaga [Gon89], Ye [Ye91] (see [Ali92]), and Monteiro and Adler [MA88].

4 Eigenvalues as semidefinite programs.

In most cases semidefinite programs arise in the form of minimizing or maximizing a linear combination of eigenvalues of a symmetric matrix subject to constraints on the matrix. In this section we study problems of this form, and show that under appropriate assumptions they are indeed special case semidefinite programs. We give primal and dual characterization of each problem and examine the complementary slackness theorem as specialized to that problem.

4.1 Minimizing sum of the first few eigenvalues. First we consider minimizing sum of the first k eigenvalues of a symmetric matrix subject to linear constraint on the matrix, namely

min{λ_{1}(X) + ⋯ + λ_{k}(X) : AvecX = b}.  (4.37)

or

min ∑_{i=1}^{k} λ_{i}(A(x)) where A(x) = A_{0} + ∑_{i=1}^{m} x_{i}A_{i}  (4.38)

To show that these problems are indeed semidefinite programs, we use the following elegant characterization by Overton and Womersley [OW91b, OW91a].

Theorem 4.1. For the sum of the first k eigenvalues of a symmetric matrix A the following semidefinite programming characterization holds:

λ_{1}(A) + ⋯ + λ_{k}(A) = max A • U
s.t. trace U = k
0 ⪯ U ⪯ I

(4.39)

Proof. See Overton and Womersley [OW91b, OW91a].

It is worth mentioning that this result is based on a beautiful convex hull characterization which was known at least as early as 1971, (see [FW71],) but unfortunately has remained somewhat obscure. Here is the statement of this result:

Lemma 4.1. Let

S_{1} := {Y Y^{T} : Y ∈ ℜ^{n×k}, Y^{T} Y = I}

and

S_{2} := {W : W = W^{T} , trace W = k, 0 ⪯ W ⪯ I}.

Then

conv S_{1} = S_{2},

and S_{1} is exactly the set of extreme points of S_{2}.

For an historical account of this result, its connection to the well-known, but computationally less useful theorem of Ky Fan, and interesting connections to the theorem of Birkhoff and Von–Neumann concerning the convex hull of doubly stochastic matrices refer to Overton and Womersley [OW91a].

Now to express (4.37) as a semidefinite program we first derive another characterization of sum of the first k eigenvalues of A, by taking the dual of (4.39).

Theorem 4.2. For the sum of the first k eigenvalues of a symmetric matrix A the following semidefinite programming characterization holds:

λ_{1}(A) + ⋯ + λ_{k}(A) = min kz + trace V
s.t. zI + V ⪰ A
V ⪰ 0

(4.40)

Now, it is easy to incorporate the equality constraints into (4.40) by replacing A with X. So (4.37) is equivalent to

min kz + trace V
s.t. AvecX = b
zI + V − X ⪰ 0
V ⪰ 0

(4.41)

and taking dual again we have the following dual characterization:

max b^{T} y
s.t. U = Mat(A^{T} y)
trace U = k
0 ⪯ U ⪯ I

(4.42)

The complementary slackness result for primal feasible z^{∗}, X^{∗}, and V^{∗}, and dual feasible U^{∗} states that these are optimal if and only if

(z^{∗}I + V^{∗} − X^{∗})U^{∗} = (I − U^{∗})V^{∗} = 0

Similarly (4.38) may be expressed by the following primal and dual pair:

min kz + trace V s.t. zI − V − ∑ x_{i}A_{i} ⪰ A_{0} V ⪰ 0  max A_{0} • Y s.t. trace Y = k A_{i} • Y = 0 for i = 1, ⋯, m 0 ⪯ Y ⪯ I  (4.43)

When k = 1, these characterization become simpler, because in that case the constraint U ⪯ I is redundant. Therefore, the problem

min{λ_{1}(X) : AvecX = b}

may be expressed as the solution of the primal and dual SDP pair:

min z s.t. zI − X ⪰ 0 AvecX = b  max b^{T} y s.t. trace Mat(A^{T} y) = 1 Mat(A^{T} y) ⪰ 0  (4.44)

and the complementary slackness theorem indicates that for X^{∗} and y^{∗} to be primal and dual optimum solution for (4.44), they must satisfy–in addition to being primal and dual feasible:

Mat(A^{T} y^{∗})(λ_{1}(X^{∗})I − X^{∗}) = 0.

4.2 Minimizing weighted sums of eigenvalues. In this section we consider the weighted sum of eigenvalues of a matrix. Let m_{1} ≥ m_{2} ≥ ⋯ ≥ m_{k} ≥ m_{k+1} = 0 be a set of fixed numbers. We are interested in the following problem:

min{m_{1}λ_{1}(X) + ⋯ + m_{k}λ_{k}(X) : AvecX = b}  (4.45)

Note that without the condition m_{1} ≥ m_{2} ≥ ⋯ ≥ m_{k} > 0 (4.45) is not necessarily a convex program. To formulate this problem as a semidefinite program, we use a technique originally employed by Donath and Hoffman in [DH73]. They rewrote the sum as follows:

m_{1}λ_{1}(A) + m_{2}λ_{2}(A) + ⋯ + m_{k}λ_{k}(A) = (m_{1} − m_{2})λ_{1}(A) +
(m_{2} − m_{3})[λ_{1}(A) + λ_{2}(A)] + ⋯ +
(m_{k−1} − m_{k})[λ_{1}(A) + ⋯ + λ_{k−1}(A)] +
m_{k}[λ_{1}(A) + ⋯ + λ_{k}(A)]  (4.46)

and observed that the right hand side of (4.46) is a nonnegative combination of convex functions, and therefore, itself is convex. This formulation also allows us to write (4.45) as a semidefinite programming problem. For each of the partial sums of eigenvalues in (4.46) we may use the relations in the last subsection and obtain the primal:

min ∑_{i=1}^{k} iz_{i} + ∑_{i=1}^{k} trace V_{i}
s.t. z_{i}I + V_{i} − (m_{i} − m_{i+1})X ⪰ 0 for i = 1, ⋯, k
AvecX = b
V_{i} ⪰ 0 for i = 1, ⋯, k

(4.47)

and the dual

max b^{T} y
s.t. A^{T} y − ∑_{i=1}^{k}(m_{i} − m_{i+1})U_{i} = 0
trace U_{i} = i for i = 1, ⋯, k
0 ⪯ U_{i} ⪯ I for i = 1, ⋯, k

(4.48)

formulations of (4.45).

The complementary slackness condition for feasible X^{∗}, z_{i}^{∗}, V_{i}^{∗}, y^{∗}, and U_{i}^{∗} for i = 1, ⋯, k to be optimal may be stated as:

(z_{i}^{∗}I + V_{i} − (m_{i} − m_{i+1})X^{∗})U_{i}^{∗} = (I − U_{i}^{∗})V_{i}^{∗} = 0 for i = 1, ⋯, k

Notice that the primal and dual characterizations (4.47) and (4.48) contain 2k semidefinite constraints each involving n × n matrices, and therefore, the interior point methods discussed earlier require O(\sqrt{kn}) iterations for each new significant digit of accuracy.

4.3 Minimizing sums of absolute-value-wise largest eigenvalues. The results of the two preceding sections may be extended to the sum of the k absolute value-wise largest eigenvalues as well. Overton and Womersley derived a dual characterization similar to (4.39); applying duality to their result we obtain the following theorem: (Recall that λ^{i}(X) is the i^{th} largest eigenvalue of X absolute value wise)

Theorem 4.3. For a symmetric matrix A the sum |λ^{1}(A)| + ⋯ + |λ^{k}(A)| is equal to optimal solution of the pair of primal and dual semidefinite programs:

max A • Y − A • W s.t. trace (Y + W) = k 0 ⪯ Y ⪯ I 0 ⪯ W ⪯ I  min kz + trace V + trace U s.t. zI + V − A ⪰ 0 zI + U + A ⪰ 0 U ⪰ 0 V ⪰ 0  (4.49)

Now to solve the optimization problem

min{|λ^{1}(X)| + ⋯ + |λ^{k}(X)| : AvecX = b}  (4.50)

we may simply add the equality constraints to the min formulation in (4.49) and then take its dual and we get the following pair of primal and dual semidefinite programs:

min kz + trace V + trace U s.t. AvecX = b zI + V − X ⪰ 0 zI + U + X ⪰ 0 U ⪰ 0 V ⪰ 0  max b^{T} y s.t. A^{T} y = Y − W trace (Y + W) = k 0 ⪯ Y ⪯ I 0 ⪯ W ⪯ I  (4.51)

The complementary slackness theorem indicates that for primal feasible z^{∗}, V^{∗}, and U^{∗}, and dual feasible Y^{∗}, and W^{∗}, complementary slackness theorem indicates that these are optimal if and only if

(z^{∗}I + V^{∗} − X^{∗})Y^{∗} = (z^{∗}I + U^{∗} + X^{∗})W^{∗} = (I − Y^{∗})U^{∗} = (I − W^{∗})V^{∗} = 0.

Again these results may be generalized to the weighted sums of absolute-value-wise largest eigenvalues. In other words, the problem

min{m_{1}|λ^{1}(X)| + ⋯ + m_{k}|λ^{k}(X)| : AvecX = b}  (4.52)

may be expressed by a primal and dual pair of semidefinite programs. First, let us ignore the equality constraints AvecX = b, and assume that X is a fixed matrix A. Then, we have

Theorem 4.4. The sum m_{1}|λ^{1}(A)|+⋯+m_{k}|λ^{k}(A)|, where A is a symmetric matrix equals the optimal solution of the primal program:

min ∑_{i=1}^{k} iz_{i} + ∑_{i=1}^{k} trace (U_{i} + V_{i})
s.t. z_{i}I + U_{i} − (m_{i} − m_{i+1})A ⪰ 0 for i = 1, ⋯, k
z_{i}I + V_{i} + (m_{i} − m_{i+1})A ⪰ 0 for i = 1, ⋯, k
U_{i} ⪰ 0 for i = 1, ⋯, k
V_{i} ⪰ 0 for i = 1, ⋯, k

(4.53) and the dual program:

max ∑_{i=1}^{k}(m_{i} − m_{i+1})(A • Y_{i} − A • W_{i})
s.t. trace (Y_{i} + W_{i}) = i for i = 1, ⋯, k
0 ⪯ Y_{i} ⪯ I for i = 1, ⋯, k
0 ⪯ W_{i} ⪯ I for i = 1, ⋯, k

(4.54)

Now we may replace A by X and impose the equality constraints on the min characterization in (4.53). After taking the dual we will have the following pair of primal and dual formulation of (4.52):

min ∑_{i=1}^{k} iz_{i} + ∑_{i=1}^{k} trace (U_{i} + V_{i})
s.t. AvecX = b
z_{i}I + U_{i} − (m_{i} − m_{i+1})X ⪰ 0 for i = 1, ⋯, k
z_{i}I + V_{i} + (m_{i} − m_{i+1})X ⪰ 0 for i = 1, ⋯, k
U_{i} ⪰ 0 for i = 1, ⋯, k
V_{i} ⪰ 0 for i = 1, ⋯, k

(4.55)

and,

max b^{T} y
s.t. A^{T} y = ∑_{i=1}^{k}(m_{i} − m_{i+1})(Y_{i} − W_{i})
trace (Y_{i} + W_{i}) = i for i = 1, ⋯, k
0 ⪯ Y_{i} ⪯ I for i = 1, ⋯, k
0 ⪯ W_{i} ⪯ I for i = 1, ⋯, k

(4.56)

Finally, the complementary slackness theorem for problem (4.52) states that primal (4.55) feasible z_{i}^{∗}, V_{i}^{∗}, and U_{i}^{∗}, and dual (4.56) feasible Y_{i}^{∗}, and W_{i}^{∗}, for i = 1, ⋯, k are optimal if and only if

(z_{i}^{∗}I + V_{i}^{∗} − (m_{i} − m_{i+1})X_{i}^{∗})Y_{i}^{∗} = (z_{i}^{∗}I + U_{i}^{∗} + (m_{i} − m_{i+1})X_{i}^{∗})W_{i}^{∗} = (I − Y_{i}^{∗})U_{i}^{∗} = (I − W_{i}^{∗})V_{i}^{∗} = 0

for i = 1, ⋯, k.

The characterization (4.39), and the max part of (4.49) were given in Overton and Womersley [OW91b]. Also, Fletcher in [Fle85] derives a closely related result to (4.39) but the result was incorrect (Fletcher had 0 ⪯ S rather than 0 ⪯ S ⪯ I.) The min characterizations as well as the primal and dual formulation of the variants with equality constraints, we believe are new.

Similar formulations can be derived for maximizing (weighted) sums of the last few smallest eigenvalues of symmetric matrices, and minimizing the first few largest singular values of arbitrary matrices. We omit these formulations here. However, maximizing the last few smallest eigenvalues of a symmetric matrix absolute-value-wise and the last few smallest singular values of an arbitrary matrix cannot be formulated as SDP because these problems are not convex programs.

5 Applications in combinatorial optimization.

The semidefinite programming problem studied in the previous sections has applications in combinatorial optimization, especially in graph theory. The connection usually is the spectral properties of graphs. In the following sections we first examine a general approach of Lovász and Schrijver which applies semidefinite programming to general zero-one integer programming problem. Then we study other applications such as the maximum stable set, the maximum induced k-partite subgraph, and graph partitioning (in particular, graph bisection) problems.

5.1 Nonlinear relaxations of 0-1 programming. Consider the homogenized integer programming problem and its linear programming relaxation:

\underline{IP} max c^{T} x s.t. a_{i}^{T} x ≥ 0 for i = 1, ⋯, m x_{i} ∈ \{0, 1\} for i = 0, ⋯, n  \underline{LP} max c^{T} x s.t. a_{i}^{T} x ≥ 0  for i = 1, ⋯, m 0 ≤ x_{i} ≤ x_{0} for i = 0, ⋯, n  (5.57)

Let P be the convex cone which is the feasible region of the LP relaxation, and ℑ(P) its integer hull (that is, ℑ(P) is the feasible region of IP.) The LP relaxation serves as a first approximation of the solution of IP. In general, this first approximation may be nonintegral and far from the actual solutions. Most effective methods of integer programming consist of adding new “cutting planes” to the LP relaxation. It seems however, that little work has been done in generating “nonlinear” but convex cuts in P. Generally such cuts may produce far better approximations than planar cuts. An ingenious approach for creating a class of nonlinear cuts has been recently proposed by Lovász and Schrijver in [LS91]. The idea is to “lift” the space from n-vectors to n × n symmetric matrices. First, we decompose the set of constraints into two sets (with possible overlap); then multiply each inequality in the first set by each inequality in the second set to obtain quadratic constraints, then replace each occurrence of x_{i}x_{j} by a new variable x_{ij} to get linear constraints again, finally impose on the matrix X = (x_{ij}) positive semidefinite constraints. If P_{1} and P_{2} are the cones defined by the first and second sets of constraints, then P = P_{1} ∩ P_{2}, and the space of matrices just defined is denoted by M_{+}(P_{1}, P_{2}). More formally, let J_{1} and J_{2} be two subsets that cover the index set of the inequality constraints in LP. Define A_{1} := A_{J_{1}}, and A_{2} := A_{J_{2}}, and P_{i} the feasible region {x : A_{i}x ≥ 0} for i = 1, 2. We require that constraints 0 ≤ x_{i} ≤ x_{0} be in both subsets. Then

M_{+}(P_{1}, P_{2}) := \{X ∈ ℜ^{\frac{n×n}{2}} : X ⪰ 0, Xe_{0} = diag(X), and (A_{1} ⊗ A_{2})vec(X) ≥ 0\}  (5.58)

Also, let N_{+}(P_{1}, P_{2}) be the n-space made up of diagonals of matrices in M_{+}(P_{1}, P_{2}), that is

N_{+}(P_{1}, P_{2}) := {diag(X) : X ∈ M_{+}(P_{1}, P_{2})}.

Then Lovász and Schrijver show that ℑ(P) ⊆ N_{+}(P_{1}, P_{2}) ⊆ P.

It is clear that optimizing a linear function on M_{+}(P_{1}, P_{2}) and N_{+}(P_{1}, P_{2}) are mixed linear and semidefinite programming problems, and interior point techniques may be applied (as long as P is given by an explicit system of inequalities.) The process just described may be quite powerful in certain combinatorial optimization problems. For instance in a general branch and bound algorithm, one may use interior point algorithms to solve the optimization problem

max{c^{T} x : x ∈ N_{+}(P_{1}, P_{2})}.

The solution then may be used as a bound for any branch and bound algorithm. This bound is stronger than the one generated by the LP relaxation, although it may be considerably more expensive to obtain.

Lovász and Schrijver show that applying N_{+} operator to the LP relaxation of the stable set polytope of a graph G = (V, E) gives bounds that are already stronger than a combination of several well-known classes of linear cuts. Recall that a stable set in a graph G = (V, E) is a subset of vertices S where each pair of vertices i and j in S are nonadjacent. Let w be a weight vector on the vertices of G, such that w_{i} is the weight of vertex i. The weighted maximum stable set problem in graphs can now be formulated as the following 0-1 program:

max w^{T} x
s.t. x_{i} + x_{j} ≤ 0 for all {i, j} ∈ E
x_{i} ∈ {0, 1} for all i ∈ V

(5.59)

Now we homogenize (5.59) by adding a new variable x_{0}, then apply the N_{+} operator with

P_{1} := P = {x : x_{i} + x_{j} ≤ x_{0} for all i, j ∈ E, and 0 ≤ x_{i} ≤ x_{0} for all i ∈ V },

and

P_{2} := {x : x_{0} − x_{i} ≥ 0, and x_{i} ≥ 0},

Finally intersect the result with hyperplane x_{0} = 1. Let the resulting set be N_{+}(STAB G). Optimization over this set is a semidefinite program and can be done in polynomial time using interior point methods (Lovász and Schrijver use the ellipsoid method to establish polynomiality). Furthermore, Lovász and Schrijver show that

STAB G ⊆ N_{+}(STAB G) ⊆ E-STAB G

where STAB G is the convex hull of all 0-1 vectors that characterize some stable set of G, and E-STAB G is the polytope associated with the LP relaxation of (5.59) (that is the polytope obtain by replacing constraints x_{i} ∈ {0, 1} by 0 ≤ x_{i} ≤ 1.) The set N_{+}(STAB G) is convex, but generally nonpolyhedral. However, Lovász and Schrijver show that the set of points in STAB G and in N_{+}(STAB G) already satisfy the following classes of well-known valid inequalities for STAB G:

1. Clique constraints. Let K be a clique in G, that is a subset of vertices every pair of which is adjacent. Let S be a stable set in G. Then clearly |1_{S} ∩ 1_{K}| ≤ 1, where 1_{S} and 1_{K} are characteristic vectors of S and K, respectively. This observation implies that for all cliques in G the inequality

1_{K}^{T} x ≤ 0  (5.60)

is valid for STAB G. Let Q-STAB G be the polytope defined by the inequalities in (5.60) and x_{i} ≥ 0. Lovász and Schrijver show that

N_{+}(STAB G) ⊆ Q-STAB G

2. Odd hole constraints. For every cycle (hole) C with 2k + 1 edges and every stable set S we know that |C ∩ S| ≤ k. Thus, for all cycles C in G the constraint

1_{C}^{T}x ≤ k  (5.61)

is valid for STAB G. Let C-STAB G be the polytope defined by all (5.61) induced by all cycles of G and inequalities x_{i} ≥ 0. Then,

N_{+}(STAB G) ⊆ C-STAB G

3. Odd anti-hole constraints. Let C̄ be a graph whose edge complement set is an odd cycle. Then the maximum stable set in G has two vertices and therefore, |C̄ ∩ S| ≤ 2 for all stable sets S. Therefore, for all antiholes C̄ in G every inequality

1T ; C̄x ≤ 2  (5.62)

is valid for STAB G. Let C̄-STAB G be the polytope defined by all inequalities (5.62) and x_{i} ≥ 0. Then,

N_{+}(STAB G) ⊆ C̄-STAB G

4. Odd wheel Constraints. Let W be a graph with 2k vertices such that vertices 1, 2, ⋯, 2k − 1 induce a cycle and vertex 2k is adjacent to all other vertices. Then W is called an odd wheel. It can be shown (see [GLS88]) that for all wheels W in G, the inequality

∑_{i=1}^{2k−1} x_{i} + (k − 1)x_{2k} ≤ k − 1.  (5.63)

is valid for STAB G. It can be shown that if W -STAB G is the set of all inequalities (5.63) and x_{i} ≥ 0, then

N_{+}(STAB G) ⊆ W -STAB G

(See [LS91].) It follows that

STAB G ⊆ N_{+}(STAB G) ⊆ Q-STAB G ∩ C-STAB G ∩ C̄-STAB G ∩ W -STAB G

and N_{+}(STAB G) already provides sharper relaxation of STAB G than any of the polytopes defined above. Yet optimization over N_{+}(STAB G) is an SDP problem and the interior point methods developed in this paper may yield practical ways of achieving strong bounds on for the maximum stable set problem.

5.2 Maximum cliques in perfect graphs. A particularly nice application of semidefinite programming is to the solution of the maximum clique problem in perfect graphs. A graph G(V, E) is called perfect if for all induced subgraphs G′ of G, the size of the maximum clique, ω(G′), equals the size of minimum proper coloring, χ(G′). (A proper coloring of vertices of a graph is an assignment of colors to each vertex such that no two adjacent vertices have the same color.) It is clear that ω(G) ≤ χ(G) for all graphs, as one needs at least ω(G) colors just to cover the vertices of the maximum clique. Several interesting properties of perfect graphs should be noted. First, the perfect graph theorem of Lovász indicates that a graph is perfect if and only if its complement is perfect, [Lov72]. This statement is equivalent to saying that for all induced subgraphs G′ of G, α(G′) = ρ(G′), where α(G′) is the size of the largest stable set in G′, and ρ(G′) is the size of the smallest number of cliques that cover all vertices of G′. Thus, in effect studying cliques in perfect graphs is equivalent to studying stables sets and any algorithm for one is valid for the other one (by simply applying it to the complementary graph.) As a consequence of the perfect graph theorem one can show that equality of maximum cliques and minimum coloring extends to the weighted graphs. More precisely, let w ∈ N^{n} be a weight vector defined on the vertices of G. A proper w-coloring of G is an assignment of colors to the vertices of G such that each vertex has at least w_{i} colors and for two adjacent vertices, their color sets are disjoint. χ(G, w) is the minimum number of colors over all proper w-colorings of G. A maximum weighted clique in G is the clique whose sum of weights of vertices is maximum; this sum is denoted by ω(G, w). A graph is perfect if and only if for all weight vectors w ∈ N^{n}, ω(G, w) = χ(G, w). Restating this for the complements of graphs, we have that a graph is perfect if and only if α(G, w) = ρ(G, w), where, α(G, w) is the weight of the maximum weighted stable set in G, and ρ(G, w) is the minimum number of cliques required to cover vertices of G such that each vertex i is in at least w_{i} cliques.

Lovász in [Lov79] discovered an invariant of graphs, θ(G, w), which has two desirable properties: first it is polynomial time computable, and second it is simultaneously an upper bound for ω(G, w) and a lower bound for χ(G, w). This invariant can be defined by a pair of primal and dual semidefinite programs. Let

M := \{X ∈ ℜ^{\frac{n×n}{2}} : X_{ij} = 0 for all i, j ∈ E or i = j\}

and

M^{⊥} := \{Y ∈ ℜ^{\frac{n×n}{2}} : Y_{ij} = 0 for all i, j ∉ E\}

Then the weighted Lovász number of G is defined by the following pair of primal and dual semidefinite programs:

θ(G, w) := min{λ_{1}(X + W) : X ∈ M}
= max{W • Y : Y ∈ M^{⊥} and trace Y = 1}

(5.64)

where W := \sqrt{w}\sqrt{w}^{T} and \sqrt{w} is an n-vector whose i^{th} component is \sqrt{w}_{i}. This min–max equality is proved directly in [GLS88], and also follows easily from the duality theory stated earlier, see (4.44).

Lemma 5.1. For every vertex weighted graph G = (V, E),

ω(G, w) ≤ θ(G, w) ≤ χ(G, w)

and

α(G, w) ≤ ϑ(G, w) := θ(Ḡ, w) ≤ ρ(G, w)

See [GLS88] chapter 9 for a thorough treatment of Lovász number of graphs including several other characterization and many interesting properties. Now our interior point algorithm can compute θ(G, w) in polynomial time; however in case of perfect graphs we have

ω(G, w) = θ(G, w) = χ(G, w)

and

α(G, w) = ϑ(G, w) = ρ(G, w)

Therefore, in case of perfect graphs, one can compute ω(G, w) and α(G, w) in polynomial time using interior point methods. By the standard self-reducibility process one can compute the maximum clique and the maximum stable set in any perfect graph, (in [GLS88] the ellipsoid method was used for this purpose.)

We now show that interior point methods give us a slightly stronger result than the ellipsoid method. More precisely, we show that computing maximum cliques (and maximum stable sets) in perfect graphs can be accomplished in O^{∗}(\sqrt{n}) randomized parallel time using P-RAM model of computation if ∥w∥_{∞} = O(n^{c}) for some constant c; (O^{∗}(\sqrt{n}) means O(\sqrt{n} log^{k} n) for some constant k.) This is straightforward. First recall that we showed a standard SDP problem can be solved in O(\sqrt{n} max(L, L′, | ln ϵ|)) iterations, if L is the number of bits in the input SDP, L′ is an à priori bound on the norm of the solution, and ϵ is the accuracy required on the size of the duality gap. In case of perfect graphs we only need to set ϵ = 1/2; in fact, if z_{k} and Y_{k} are our current primal and dual solution where there is only one integer between z_{k} and W • Y_{k} then we can stop and declare θ(G, w) = ⌈z_{k}⌉ = ⌊W • Y_{k}⌋. Furthermore, L = O(log n) since all coefficients in the primal-dual characterization of θ(G, w) in (5.64) are either zero or one or w_{i}w_{j}. Finally, L′ = O(log n) because at most the whole graph may be a clique and so its weight is ∑ w_{i}. Thus computing θ(G, w) requires O^{∗}(\sqrt{n}) iterations. Each iteration essentially involves solving a system of linear equations which is already known to be in complexity class NC, that is requires O^{∗}(1) time with polynomial number of processors. Therefore, computing θ(G, w) for polynomially bounded w requires O^{∗}(\sqrt{n}) operations on a P-RAM model of computation.

It remains to show that computing the maximum clique itself can be accomplished in O^{∗}(\sqrt{n}). We cannot use the self reducibility process here since it may require O(n) time even on a P-RAM machine. However, observe that if the maximum clique is unique then we can compute it in O^{∗}(\sqrt{n}) parallel time. One could remove one vertex i of the graph and compute θ(G\i, w) for the remaining graph. The vertex i is in the unique maximum clique if and only if ω(G\i, w) < ω(G, w). Therefore, testing this simultaneously for all vertices we get the set of vertices in the maximum clique. When we do not have uniqueness, we may use the randomized perturbation scheme of Mulmuley, Vazirani and Vazirani, [MVV87]. First recall their isolating lemma:

Lemma 5.2. Let S = {x_{1}, ⋯, x_{n}} and F a family of subsets of S, that is F = {S_{1}, ⋯, S_{N}}. Further, let elements of S be assigned integer weights chosen uniformly and independently from [1, 2n]. Then,

Pr[There is a unique maximum weight set in F] ≥ \frac{1}{2}.

See [MVV87] for proof.

To get an (unweighted) maximum clique in a perfect graph we follow a procedure similar to the one adopted by Mulmuley, Vazirani and Vazirani for constructing the minimum weighted perfect matching in bipartite graphs. The idea is to assign weights to vertices randomly so that with high probability the maximum weighted clique is unique, but at the same time, this clique is among the maximum cliques before assigning weights.

First give a weight of 2n^{2} to each vertex i so that the weight of maximum weighted cliques is at least 2n^{2} more than the next largest clique weight. Then perturb weight of each vertex i by adding integer u_{i} uniformly and independently chosen from integers in [1, 2n]. So now each vertex has weight w_{i} = 2n^{2} + u_{i}. Notice that if a clique was not maximum before, then it is impossible for it to become maximum after assigning weights. Therefore, the maximum weighted clique is among one of the unweighted maximum cliques. The isolating lemma implies that this clique is unique with a probability at least 1/2 and we may use the scheme mentioned at the beginning of this section to find it in parallel.

We should mention that this scheme, in fact, results in a Las Vegas type randomized algorithm. No randomization is involved in computing the size oracle, ω(G, w); only constructing a maximum clique involves probabilistic choices. If the weights generated do not result in a unique maximum weighted clique, the scheme mentioned at the beginning of this section may return a set which is not a maximum clique. This can be checked in parallel and the algorithm will return a message of failure; any set returned by the algorithm is a genuine maximum clique with no possibility of error. We summarize these results in the following theorem:

Theorem 5.1. Let G = (V, E) be a perfect graph with weight vector w on its vertices. Let also that ∥w∥_{∞} = O(n^{c}) for some constant c. Then one can compute the maximum weighted clique and the maximum weighted stable set of G in O^{∗}(\sqrt{n}) Las Vegas randomized parallel time using a P-RAM model of computation.

5.3 The maximum induced k-partite subgraph problem. In [NM90] G. Narasimhan and R. Manber generalized the concept of the Lovász number of graphs as follows: Let G = (V, E) be a graph. Let α_{k}(G) be the size of the largest induced k-partite graph in G. Recall that ρ(G) is the minimum number of cliques that can cover all vertices of G. Then Narasimhan and Manber show that

α_{k}(G) ≤ ϑ_{k}(G) := min_{X∈M} ∑_{i=1}^{k} λ_{i}(X + J) ≤ kρ(G)  (5.65)

where J is the matrix of all 1’s. For k = 1 ϑ_{k} reduces to the Lovász number ϑ. It is clear now that computing θ_{k}(G) is an SDP problem and may be solved by interior point methods. Taking dual of (5.65) we get

ϑ_{k}(G) = max J • Y
s.t. trace Y = k
Y ∈ M^{T}
0 ⪯ Y ⪯ I

(5.66)

It is not difficult to extend the bound of Narasimhan and Manber to the weighted case. Let w be a weight vector over the vertices of G and α_{k}(G, w the maximum weight k-colorable induced subgraph of G.

Theorem 5.2. Let G = (V, E) be a graph with a weight vector w ∈ N^{n} defined on its vertices. Set W = (\sqrt{w})(\sqrt{w})^{T}. Then α_{k}(G, w) ≤ ϑ(G, w, where ϑ(G, w is defined as

ϑ(G, w) = min{∑_{1}^{k} λ_{i}(X + W) : X ∈ M}
= max{W • Y : Y ∈ M^{⊥} and trace Y = k}

(5.67)

Proof. The proof is identical to [GLS84] Theorem (4.20) and is omitted.

Let ℧ be the class of graphs for which α_{k}(G′) = θ_{k}(G′) for all induced subgraphs G′. Then the sublinear parallel time algorithm of theorem 5.1 may be extended to solve the largest induced k-partite subgraph problem for the graphs in class ℧. It remains an interesting open problem to fully characterize ℧.

5.4 The graph partitioning problem. Suppose we are given a set of integers m_{1} ≥ m_{2} ≥ ⋯ ≥ m_{k}, with ∑_{j} m_{j} = n. Let also that G = (V, E) be a complete edge-weighted graph with n vertices and each edge {i, j} with weight w_{ij}. We want to partition vertices of G into k subsets such that the j^{th} subset has cardinality m_{j}, and that sum of the weights of those edges whose endpoints are in different subsets is minimized. Let us denote this minimum number by π_{k}(G). Computing π_{k}(G) is of course NP-hard. Hoffman and Donath in [DH72] and [DH73] derive the following lower bound on the size of the minimum partition (see also Barnes and Hoffman [BH84]). Let A be a matrix with A_{ij} = w_{ij}. Then Donath and Hoffman prove the following relation [DH73]:

π_{k}(G) ≥ −\frac{1}{2} {min \atop 1^{T} x=a} ∑_{j=1}^{k} m_{j}λ_{j}(A + Diag x)  (5.68)

where a := − ∑ A_{ij}. Again it is clear that computing this bound is an SDP problem. Using results from section 4 and after considerable simplifications we may formulate it by the following pair of primal and dual SDP programs:

max ∑_{i=1}^{k}(m_{i} − m_{i+1})A • U_{i} + a s.t. trace U_{i} = i  for i = 1, ⋯, k ∑_{i=1}^{k}(m_{i} − m_{i+1})(U_{i})_{jj} = 1 for j = 1, ⋯, n 0 ⪯ U_{i} ⪯ I  for i = 1, ⋯, k ,  (5.69)

and

min ∑_{i=1}^{k} iz_{i} + 1^{T} x + ∑_{i=1}^{k} trace V_{i} s.t. z_{i}I + V_{i} − (m_{i} − m_{i+1})Diagx ⪰ (m_{i} − m_{i+1})A for i = 1, ⋯, k V_{i} ⪰ 0 for i = 1, ⋯, k .  (5.70)

Barnes and Hoffman in [BH84] describe how to use the eigenvectors associated with the k largest eigenvalues of the optimal matrix A + Diag x^{∗}, to generate a partition of the nodes of the graph. See also Barnes [Bar82a] and [Bar82b].

An important special case of the graph partitioning problem is the graph bisection problem. In this case, m_{1} = m_{2} = n/2 (assuming n is even), and the pair of primal and dual problems may further be simplified as follows:

min 2z + trace V + a
s.t. zI + V − Diag x ⪰ A
V ⪰ 0

(5.71)

And the dual is:

max A • Y + a
s.t. Y_{ii} = \frac{2}{n}  for i = 1, ⋯, n
0 ⪯ Y ⪯ I

(5.72)

R. Boppana in [Bop87] considers the graph bisection problem and using the primal characterization shows that one may get the optimal bisection with high probability. The graph bisection problem has important applications in the VLSI routing problem. Combining the SDP formulation of Hoffman and Donath, favorable average case analysis of Boppana, and the interior point technique developed in this paper may result in an effective and practical method for solving this problem.

Acknowledgment

Discussions with Stephen Boyd, Arkadii Nemirovsky, Michael Overton, Joshi Ramana, Rob Womersley and Yinyu Ye, have been most useful in preparing this article and removing several errors in the preliminary versions.

References

[Ali91a] F. Alizadeh. Combinatorial Optimization with Interior Point Methods and Semi-Definite Matrices. PhD thesis, University of Minnesota, Minneapolis, Minnesota, 1991.

[Ali91b] F. Alizadeh. A sublinear-time randomized parallel algorithm for the maximum clique problem in perfect graphs. In Proc. 2nd ACM-SIAM Symposium on Discrete Algorithms, 1991.

[Ali92] F. Alizadeh. Optimization Over Positive Semi-Definite Cone; Interior-Point Methods and Combinatorial Applications. In P. Pardalos, editor, Advances in Optimization and Parallel Computing. North–Holland, 1992.

[Bar82a] E. R. Barnes. An algorithm for partitioning the nodes of a graph. SIAM J. Alg. and Disc. Meth., 3, 1982.

[Bar82b] E. R. Barnes. Partitioning the nodes of a graph. In L. Lesniak, D. R. Lick, and C. E. Wall, editors, Graph Theory with Applications to Algorithms and Computer Science. John Wiley, [BH84] E. Barnes and A. J. Hoffman. Partitioning, spectra, and linear programming. In W. E. Pulley-

blank, editor, Progress in Combinatorial Optimization. Academic Press, 1984.

[BICK69] A. Ben-Israel, A. Charnes, and K. O. Kortanek. Duality and Asymptotic Solvability over Cones. Bulletin of American Mathematical Society, 75(2):318–324, 1969.

[BJS90] M. Bazaraa, J. Jarvis, and H. Sherali. Linear Programming and Network Flows. John Wiley, 1990.

[Bop87] R. B. Boppana. Eigenvalues and graph bisection: an average case analysis. In Proc. 28th IEEE Annual Symposium on Foundations of Computer Science, 1987.

[BW81a] J. Borwein and H. Wolkowicz. Characterization of optimality for the abstract convex program with finite dimensional range. J. Austral. Math. Soc. series A, 30:390–411, 1981.

[BW81b] J. Borwein and H. Wolkowicz. Facial reduction for a cone-convex programming problem. J. Austral. Math. Soc. series A, 30:369–380, 1981.

[CDW75] J. Cullum, W. E. Donath, and P. Wolfe. The minimization of certain nondifferentiable sums of eigenvalue problems. Mathematical Programming Study, 3:35–55, 1975.

[DH72] W. E. Donath and A. J. Hoffman. Algorithms for partitioning graphs and computer logic based on eigenvectors of connection matrices. IBM Technical Disclosures Bulletin, 15, 1972.

[DH73] W. E. Donath and A. J. Hoffman. Lower bounds for the partitioning of graphs. IBM Journal of Research and Development, 17, 1973.

[Duf56] R. J. Duffin. Infinite Programs. In H. W. Kuhn and A. W. Tucker, editor, Linear inequalities and Related Systems, pages 157–170. Princeton University Press, Princeton, NJ, 1956.

[Fle85] R. Fletcher. Semi-definite matrix constraints in optimization. SIAM J. Control Optim., 23:493– 513, 1985.

[FW71] P. A. Fillmore and J. P. Williams. Some convexity theorems for matrices. Glasgow Mathematical Journal, 12:110–117, 1971.

[GLS84] M. Grötschel, L. Lovász, and A. Schrijver. Polynomial algorithms for perfect graphs. In C. Berge and V. Chvatal, editors, Perfect Graphs. North Holland, 1984. Annals of Discrete Mathematics, 21.

[GLS88] M. Grötschel, L. Lovász, and A. Schrijver. Geometric Algorithms and Combinatorial Optimization. Springer Verlag, 1988.

[Gon89] C. C. Gonzaga. An Algorithm for Solving Linear Programming in O(n^{3}L) Operations. In N. Megiddo, editor, Progress in Mathematical Programming, pages 1–28. Springer Verlag, Berlin, 1989.

[GPST91] A. V. Goldberg, S. A. Plotkin, D. Shmoys, and É. Tardos. Interior-Point Methods in Parallel Computation. SIAM J. Comput., 21(1):149–150, 1991.

[Gra81] A. Graham. Kronecker Products and Matrix Calculus: with Applications. Ellis Horwood, London, 1981.

[HJ85] R. Horn and C. Johnson. Matrix Analysis. Cambridge University Press, 1985.

[Kar84] N. Karmarkar. A New Polynomial-Time Algorithm for Linear Programming. Combinatorica, 4:373–395, 1984.

[KMY89] M. Kojima, S. Mizuni, and A. Yoshise. A Primal-Dual Interior-Point Algorithm for Linear Programming. In N. Megiddo, editor, Progress in Mathematical Programming. Springer Verlag, Berlin, 1989.

[Lov72] L. Lovász. Normal hypergraphs and the weak perfect graph conjecture. Discrete Mathematics, 2, 1972.

[Lov79] L. Lovász. On the shannon capacity of a graph. IEEE Trans. Info. Theory, 25(1), January 1979.

[LS91] L. Lovász and A. Schrijver. Cones of Matrices and Setfunctions, and 0-1 Optimization. SIAM J. Optimization, 1(2), 1991.

[MA88] R. Monteiro and I. Adler. Polynomial-time primal-dual affine scaling algorithm for linear and convex quadratic programming and its power series extension. Technical Report ESRC 88-8, Industrial Engineering and Operations Research Department, University of California–Berkeley, March 1988.

[Meg89] N. Megiddo. Pathways to the optimal set in linear programming. In N. Megiddo, editor, Progress in Mathematical Programming. Springer Verlag, Berlin, 1989.

[MVV87] K. Mulmuley, U. V. Vazirani, and V. V. Vazirani. Matching is as Easy as Matrix Inversion. Combinatorica, pages 105–131, 1987.

[NM90] G. Narasimhan and R. Manber. A Generalization of Lovász’s θ Function. In W. Cook, and P. D. Seymour, editor, Polyhedral Combinatorics, volume 1 of DIMACS series in Discrete Mathematics and Theoretical Computer Science, pages 19–27. American Mathematical Society and Association for Computing Machinery, 1990.

[NN92] Y. Nesterov and A. Nemirovskii. Interior Point Polynomial Methods in Convex Programming: Theory and Applications. SIAM, 1992. To be published.

[Ove92] M. Overton. Large-scale optimization of eigenvalues. SIAM J. Optimization, 2(1):88–120, 1992.

[OW91a] M. Overton and R. Womersley. On the sum of the largest eigenvalues of a symmetric matrix. Report 550, NYU Computer Science Department, March 1991. to appear in SIAM J. Matr. Anal.

[OW91b] M. Overton and R. Womersley. Optimality conditions and duality theory for minimizing sums of the largest eigenvalues of symmetric matrices. Report 556, NYU Computer Science Department, June 1991.

[Roc70] T. Rockafellar. Convex Analysis. Princeton University Press, 1970.

[Sch86] A. Schrijver. Theory of Linear and Integer Programming. J. Wiley & Sons, 1986.

[Sha85] A. Shapiro. Extermal problems on the set of nonnegative definite matrices. Linear Algebra and its Applications, 67:7–18, 1985.

[Wol81] H. Wolkowicz. Some applications of optimization in matrix theory. Linear Algebra and its Applications, 40:101–118, 1981.

[Ye90] Y. Ye. A class of projective transformations for linear programming. SIAM J. Comput., 19(3), 1990.

[Ye91] Y. Ye. An O(n^{3}L) Potential Reduction Algorithm for Linear Programming. Math. Prog., 50(2),