A computational study of exact subgraph based SDP bounds for Max-Cut, stable set and coloring

The “exact subgraph” approach was recently introduced as a hierarchical scheme to get increasingly tight semidefinite programming relaxations of several NP-hard graph optimization problems. Solving these relaxations is a computational challenge because of the potentially large number of violated subgraph constraints. We introduce a computational framework for these relaxations designed to cope with these difficulties. We suggest a partial Lagrangian dual, and exploit the fact that its evaluation decomposes into several independent subproblems. This opens the way to use the bundle method from non-smooth optimization to minimize the dual function. Finally computational experiments on the Max-Cut, stable set and coloring problem show the excellent quality of the bounds obtained with this approach.


Introduction
The study of NP-hard problems has led to the introduction of various hierarchies of relaxations, which typically involve several levels.Moving from one level to the next the relaxations get increasingly tighter and ultimately the exact optimum may be reached, but the computational effort grows accordingly.
Among the most prominent hierarchies are the polyhedral ones from Boros, Crama and Hammer [5] as well as the ones from Sherali and Adams [32], Lovász and Schrijver [25] and Lasserre [22] which are based on semidefinite programming (SDP).Even though on the starting level they have a simple SDP relaxation, This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 764759 and the Austrian Science Fund (FWF): I 3199-N31 and P 28008-N35.We thank three anonymous referees for their constructive comments which substantially helped to improve the presentation of our material.An extended abstract containing the theoretical foundations of this work appeared as [13].This article now contains more details, additional information and many new computational results.
already the first nontrivial level in the hierarchy requires the solution of SDPs in matrices of order n 2 and on level k the matrix order is n O(k) .Hence they are considered mainly as theoretical tools and from a practical point of view these hierarchies are of limited use.
Not all hierarchies are of this type.In [5] a polyhedral hierarchy for the Max-Cut problem is introduced which maintains n 2 variables at all levels, with a growing number of constraints.More recently, Adams, Anjos, Rendl and Wiegele [1] introduced a hierarchy of SDP relaxations which act in the space of symmetric n × n matrices and at level k of the hierarchy all submatrices of order k have to be "exact" in a well-defined sense, i.e. they have to fulfill an exact subgraph constraint (ESC).
It is the main purpose of this paper to describe an efficient way to optimize over level k of this hierarchy for small values of k, e.g.k 7, and demonstrate the efficiency of our approach for the Max-Cut, stable set and coloring problem.These investigations were started in [12,13] and here we offer the full picture.
Maintaining n k possible ESCs in an SDP in matrices of order n is computationally infeasible even for k = 2 or k = 3, because each ESC creates roughly k 2 additional equality constraints and at most 2 k additional variables.
We suggest the following ideas to overcome this difficulty.First we proceed iteratively, and in each iteration we include only (a few hundred of) the most violated ESCs.More importantly, we propose to solve the dual of the resulting SDP.The structure of this SDP with ESCs admits a reformulation of the dual in the form of a non-smooth convex minimization problem with attractive features.First, any dual solution yields a valid bound for our relaxations, so it is not necessary to carry out the minimization to optimality.Secondly, the dual function evaluation decomposes into two independent problems.The first one is simply a sum of max-terms (one for each ESC), and the second one consists in solving a "basic" SDP, independent of the ESCs.The optimizer for this second problem also yields a subgradient of the objective function.With this information at hand we suggest to use the bundle method from non-smooth convex optimization.It provides an effective machinery to get close to a minimizer in few iterations.
As a result we are able to get near optimal solutions where all ESCs for small values of k (k 7) are satisfied up to a small error tolerance.Our computational results demonstrate the practical potential of this approach.
The paper is organized as follows.In Section 2 we briefly describe the Max-Cut, the stable set and the coloring problem along with their semidefinite relaxations, which are well-studied in the literature.Section 3 recalls the exact subgraph hierarchy, described in [1].We introduce a unified setting for all these problems and take a look at their structural properties.In Section 4 we reformulate the SDP and consider a partial Lagrangian dual.It results in many subproblems, separating the basic SDP part from the ESC part.The bundle method from non-smooth optimization is described in Section 5 as an attractive algorithmic framework to deal with the subproblems in the partial Lagrangian dual.In Section 6 we describe our algorithm in order to obtain exact subgraph based SDP bounds.We argue in Section 7 that standard SDP solvers are only of limited use when dealing with our ESC hierarchy and present extensive computational results.Finally we close with conclusions and future work in Section 8.
We finish this introductory section with some notation.We denote the vector of all-ones of size n with 1 n and 2 Combinatorial Problems and Semidefinite Relaxations

The Max-Cut Problem
In the Max-Cut problem a symmetric matrix L ∈ S n is given and c ∈ {−1, 1} n which maximizes c T Lc should be determined.
If the matrix L corresponds to the Laplacian matrix of a (edge-weighted undirected) graph G, this is equivalent to finding a partition of the vertices of G into two subsets such that the total weight of the edges joining these two subsets is maximized.Such an edge set is also called a cut in G.
Partitions of N into two subsets can be expressed as c ∈ {−1, 1} n where the two subsets of N correspond to the entries of c with the same sign.Given c ∈ {−1, 1} n we call C = cc T a cut matrix.The convex hull of all cut matrices (of order n) is denoted by or simply CUT if the dimension is clear from the context.Since c T Lc = L, cc T the Max-Cut problem can also be written as the following (intractable) linear program z mc = max{ L, X : X ∈ CUT}.
CUT is contained in the spectrahedron is a basic semidefinite relaxation for Max-Cut.This model is well-known, attributed to Schrijver and was introduced in a dual form by Delorme and Poljak [8].It can be solved in polynomial time to a fixed prescribed precision and solving this relaxation for n = 1000 takes only a few seconds.It is well-known that the Max-Cut problem is NP-hard.On the positive side, Goemans and Williamson [14] show that one can find a cut in a graph with nonnegative edge weights of value at least 0.878z mc in polynomial time.

The Stable Set Problem
In the stable set problem the input is an unweighted graph G.We call a subset of the vertices stable, if no two vertices are adjacent.Moreover we call a vector s ∈ {0, 1} n a stable set vector if it is the incidence vector of a stable set.The convex hull of all stable set vectors of G is denoted with STAB(G).In the stable set problem we want to determine the stability number α(G), which denotes the cardinality of a largest stable set in G, hence Furthermore we denote with the convex hull of all stable set matrices ss T .Then with the arguments of Gaar [12] it is easy to check that which is known as the theta body in the literature.Therefore is a relaxation of the stable set problem.The Lovász theta function ϑ(G) was introduced in a seminal paper by Lovász [24].We refer to Grötschel, Lovász and Schrijver [15] for a comprehensive analysis of ϑ(G).
Determining α(G) is again NP-hard.Contrary to Max-Cut, which has a polynomial time .878-approximation,for every ε > 0 there can be no polynomial time algorithm that approximates α(G) within a factor better than O(n 1−ε ) unless P = N P , see Håstad [17].

The Vertex Coloring Problem
The coloring problem for a given graph G consists in determining the chromatic number χ(G), which is the smallest t such that N can be partitioned into t stable sets.Let S = (s 1 , . . ., s k ) be a matrix where each column s i is a stable set vector and the corresponding stable sets partition N into k sets.Let us call such matrices S stable-set partition matrices (SSPM) and denote by |S| the number of columns of S or equivalently the number of stable set vectors of S. The n × n matrix X = SS T is called coloring matrix.The convex hull of the set of all coloring matrices of G is denoted by COL(G) = conv {X : X is a coloring matrix of G} .
We also need the extended coloring polytope The difficult set COL ε can be relaxed to the easier spectrahedron and we can consider the semidefinite program t * (G) = min t : Obviously t * (G) χ(G) holds because the SSPM S consisting of χ(G) stable sets yields a feasible coloring matrix X = SS T with objective function value χ(G).It is in fact a consequence of conic duality that t * (G) = ϑ(G) holds.
It is NP-hard to find χ(G), to find a 4-coloring of a 3-colorable graph [16] and to color a k-colorable graph with O(k log k 25 ) colors for sufficiently large k, [20].

Definition of the Hierarchy
In this section we discuss how to systematically tighten the relaxations (1), ( 2) and (3) with "exactness conditions" imposed on small subgraphs.We obtained the relaxations by relaxing the feasible regions CUT, STAB 2 and COL of the integer problem to simple spectrahedral sets.Now we will use small subgraphs to get closer to the feasible regions of the original problems again.
For I ⊆ N let k I = |I| be the cardinality of I. Furthermore let G I be the induced subgraph of G on the set of vertices I.If X is the n × n matrix from the relaxations (1), (2) or (3), then we denote with X I the principal k I × k I submatrix of X corresponding to the rows and columns in I.Note that X I is the submatrix of X corresponding to G I .
We first look at the exact subgraph relaxations for Max-Cut.Adams, Anjos, Rendl and Wiegele [1] introduced additional constraints for the Max-Cut relaxation (1) in the following way.The exact subgraph constraint (ESC) for I ⊆ N requires that the matrix X I corresponding to the subgraph G I lies in the convex hull of the cut matrices of G I , that is The ESC for I can equivalently be phrased as for some λ ∈ ∆ t I where C I i is the i-th cut matrix of the subgraph G I and t I is the total number of cut matrices.If X is a solution of (1) that fulfills the ESC for some I we say that X is exact on I and X is exact on G I .Now we want the ESCs to be fulfilled not only for one but for a certain selection of subgraphs.We denote with J the set of subsets I, on which we require X to be exact, and get the following SDP relaxation with ESCs for Max-Cut.max{ L, X : Before we give theoretical justification that ( 4) is worth to be investigated, we present the ESCs for the other problems.We start with the stable set problem on a graph G and its relaxation (2).In this case the ESC for I ⊆ N , and hence for the subgraph G I , requires that X I ∈ STAB 2 (G I ) holds and the SDP with ESCs for the stable set problem is max{trace(X) : X ∈ X S , X I ∈ STAB 2 (G I ) ∀I ∈ J}. ( Turning to the coloring problem, we analogously impose additional ESCs of the form X I ∈ COL(G I ) to obtain the SDP with ESCs min t : We now want to investigate the properties of (4), ( 5) and (6).Towards that end we define the k-th level of the exact subgraph hierarchy according to [1] by using J = {I ⊆ N : |I| = k} in the SDPs (4), ( 5) and (6), respectively.We denote the corresponding objective function values with z k mc , z k ss and z k c .So in other words the k-th level of the exact subgraph hierarchy is obtained by forcing all subgraphs on k vertices to be exact in the basic SDP relaxation.
Note that holds for all k ∈ {2, . . ., n}, see [1,12].Hence (4) and ( 5) are relaxations of Max-Cut and the stable set problem.Furthermore it can be verified that holds for all k ∈ {2, . . ., n}, so for the coloring problem we do not necessarily reach χ(G) at the n-th level.However, the following holds.Let z k cε be the optimal objective function value if we add the inequalities t where λ I ∈ ∆ t I is a variable for the convex combination for each subgraph G I to the SDP for z k c .Then z n cε = χ(G) holds.Hence z k c is a relaxation of z k cε , which is in turn a relaxation of the coloring problem.As a result it is clear that it makes sense to investigate (4), ( 5) and (6).
Note that in the case of the stable set and the coloring problem the polytopes STAB 2 (G I ) and COL(G I ) depend on the subgraph G I , whereas in Max-Cut the polytope CUT |I| only depends on the number of vertices of G I .
Finally let us mention that an important feature of this hierarchy is that the size of the matrix variable remains n or n + 1 on all levels of the hierarchy.On higher levels the ESCs are included into the SDPs in the most natural way through convex combinations.Hence on higher levels of the exact subgraph hierarchy new variables and linear constraints representing convex hull conditions are added to the SDP of the basic SDP relaxation.
Therefore it is possible to approximate z k mc , z k ss and z k c by forcing only some subgraphs of order k to be exact.This is our key ingredient to computationally obtain tight bounds on z mc , α(G) and χ(G) and also a major advantage over several other SDP based hierarchies [22,25,32] for NP-hard problems.

Structural Differences of the Three Problems
The focus of this paper lies in computational results, so we omit further extensive theoretical investigations, but we want to draw the attention to a major structural difference between the Max-Cut problem and the stable set and the coloring problem.Towards this end we consider one graph from the Erdős-Rényi model G(n, p) with n = 100 and p = 0.15.A graph from this model is a random graph of order n, in which each edge appears with probability p.
We compute the optimal solutions of the basic relaxations ( 1), ( 2) and ( 3) and denote them by X * .Then for each subgraph G I of order k ∈ {2, 3, 4, 5} we compute the projection distance δ I mc , δ I ss and δ I c of the submatrix X * I of the corresponding X * to CUT k , STAB(G I ) and COL(G I ), respectively.So for example where .denotes the Frobenius norm.We consider a subgraph G I as violated, if the projection distance is larger than the small tolerance 5 • 10 −5 .In Table 1 one sees that the number of violated subgraphs is much higher in the case of the Max-Cut problem than for the stable set and the coloring problem. Figure 1, 2 and 3 show the distribution of the projection distances of the violated subgraphs.They are normalized in such a way that 1 is the total number of violated subgraphs.Here it becomes obvious that for the Max-Cut problem most of the violated subgraphs have a large violation, whereas most of the violated subgraphs for the coloring problem have a small violation and an even smaller violation for the stable set problem.Therefore in the case of the Max-Cut problem there are very many violated subgraphs, and typically all of them have a large projection distance.On the other hand for the stable set and the coloring problem only very few subgraphs have a large projection distance, the majority of the subgraphs is either not violated at all or only violated a little bit.Hence finding significantly violated subgraphs is much more difficult for the stable set and the coloring problem, than it is for the Max-Cut problem.
A possible explanation for this consists of the following dimension argument.Let G be a graph on n vertices with m edges.The SDP relaxation for Max-Cut starts out with a matrix variable of size n and n equations, while the evaluation of ϑ(G) requires a matrix of size n + 1 and n + m + 1 equations and in the computation of t * (G) there is a matrix of size n + 1 and 2n + m equations.Hence the Max-Cut, stable set and coloring relaxation are contained in a n 2 , n 2 + n − m and n 2 − m + 1 dimensional space, and it makes sense that Max-Cut has the most and coloring has the least violated ESCs, just as we see it in Table 1.Furthermore in the stable set and the coloring relaxation the additional row and column together with the positive semidefiniteness constraint effect all entries of X, even if they are not directly addressed by any constraint.Therefore it is plausible that the violations for the Max-Cut problem are much larger than those for the stable set and the coloring problem.
For our computations that means that there is the hope that fewer ESCs are necessary to tighten the basic relaxation.This intuition is indeed confirmed in our computational experiments in Section 7.

Partial Lagrangian Dual
We are interested in solving relaxations (4), ( 5) and ( 6) with a potentially large number of ESCs, where using interior point solvers is too time consuming.In this section we will first establish a unified formulation of the relaxations (4), ( 5) and (6).Then we will build the partial Lagrangian dual of this formulation, where only the ESCs are dualized.
In order to unify the notation for the three problems observe that the ESCs X I ∈ CUT |I| , X I ∈ STAB 2 (G I ) and X I ∈ COL(G I ) can be represented as where C I i is the i-th cut, stable set or coloring matrix of the subgraph G I and t I is their total number.
A formal description of ESC in (7) requires some additional notation.First we introduce the projection P I : S n → S k I , mapping X to the submatrix X I .Second we define a map A I : S k I → R t I , such that its adjoint map A I : R t I → S k I is given by A I (λ) = t I i=1 λ i C I i and produces a linear combination of the cut, stable set or coloring matrices.Thus we can rewrite (7) as The left-hand side of this matrix equality is a symmetric matrix, of which some entries (depending on which problem we consider) are zero for sure, so which are b I + 1 equalities and t I inequalities.In consequence all three relaxations (4), ( 5) and ( 6) have the generic form where C, X , A I , M I and b I have to be defined in a problem specific way.Furthermore X = X in the case of Max-Cut and stable set and X = t 1 T 1 X for coloring, but for the sake of understandability we will just use X in the following.
The key idea to get a handle on problem ( 9) is to consider the partial Lagrangian dual where the ESCs (without the constrains λ I ∈ ∆ t I ) are dualized.We introduce a vector of multipliers y I of size b I for each I and collect them in y = (y I ) I∈J and also collect λ = (λ I ) I∈J .The Lagrangian function becomes and standard duality arguments (Rockafellar [31,Corollary 37 For a fixed set of multipliers y the inner maximization becomes max This maximization is interesting in at least two aspects.First, it is separable in the sense that the first term depends only on X and the second one only on the separate λ I .Moreover, if we denote the linear map A I M I : R b I → R t I with the matrix D I , maximizing the summands of the second term is easy, because the feasible region is a simplex.Hence the explicit solution of maximizing a summand of the second term is max In order to consider the first term in more detail, we define the following function.Let b = I∈J b I be the dimension of y.Then h : R b → R is defined as where X * is a maximizer over the set X for y fixed.Note that h(y) is convex but non-smooth, but (12) shows that is a subgradient of h with respect to y I .With (11) and ( 12) we reformulate the partial Lagrangian dual (10) to The dual formulation (14) of the original semidefinite relaxation (9) has the form of a convex minimization problem over the set of multipliers y.The evaluation of the function h at a given y requires solving a "simple" SDP, independent of the number of ESCs included in the relaxation.The function evaluation also provides a subgradient of h at y, given in (13).Hence we propose to use the bundle method from convex optimization to solve (14).The details are given in the subsequent section.

The Bundle Method
The bundle method is a well established tool in convex optimization to minimize a non-smooth convex function.We refer to the recent monograph Bonnans, Gilbert, Lemaréchal and Sagastizábal [4] for a nice introduction.In our setting we want to use the bundle method in order to solve an SDP.Helmberg and Rendl [18] were the first to use a bundle method to solve SDPs in 2000.Later Fischer, Gruber, Rendl and Sotirov [10] and Rendl and Sotirov [29] used the bundle method for SDPs in order to get good relaxations for the Max-Cut and the equipartition problem and the quadratic assignment problem, respectively.The bundle method setting described by Frangioni and Gorgone in [11], which is set up to handle max terms explicitly, is best suited for our purposes, so we apply it to our problem (14).
The bundle method is an iterative procedure.It maintains the current center y, representing the current estimate of the optimal solution, and the set B = {(y 1 , h 1 , g 1 , X 1 ), . . ., (y r , h r , g r , X r )}, which is called bundle, throughout the iterations.Here y 1 , . . ., y r are the points which we use to set up our subgradient model.Moreover h j = h(y j ), g j is a subgradient of h at y j and X j is a maximizer of h at y j as in (12).
At the start we select y 1 = y = 0 and evaluate h at y, which yields the bundle B = {(y 1 , h 1 , g 1 , X 1 )}.A general iteration consists of first determining the new trial point, then evaluating the function at this new point, and finally updating the bundle B. In the literature evaluating the function is referred to as calling the oracle.The subgradient information of the bundle B translates into the subgradient model h(y) h j + g j , y − y j for all j = 1, . . ., r.
It is common to introduce e j = h(y) − h j − g j , y − y j for j = 1, . . ., r and to define e = (e j ) j=1,...,r .With h = h(y) the subgradient model becomes h(y) max h − e j + g j , y − y .
The right-hand side above is convex, piecewise linear and minorizes h.In each iteration of the bundle method we minimize the right-hand side of (15) instead of h, but ensure that we do not move too far from y by adding a penalty term of the form 1 2 µ y − y 2 2 for a parameter µ ∈ R + to the objective function.We introduce auxiliary variables w ∈ R and v I ∈ R for all I ∈ J to model the maximum terms.With q = |J| and v = (v I ) I∈J ∈ R q we end up with min y,w,v This is a convex quadratic problem in 1 + q + b variables with r + I∈J t I linear inequality constraints which is often referred to as the bundle master problem.Its solution ( y, w, v) provides the new trial point y.In the following section we will discuss computational issues and present a practically efficient approach starting with its dual, see below.The second step in each bundle iteration is to evaluate the function h at y which means solving the basic SDP relaxation as introduced in Section 2 with a modified objective function.In the case of Max-Cut this function evaluation can be done very quickly (solve an SDP with n simple equations).For the stable set and the coloring problem the resulting SDP is computationally more demanding, as there are also equations for each edge in the graph.The bundle iteration is finished by deciding whether y becomes the new center (serious step, roughly speaking if the increase of the objective function is good enough) or not (null step).In either case the new point is included in the bundle, some other elements of the bundle are possibly removed, the bundle parameter µ is updated and a new iteration starts.

The Dual of the Bundle Master Problem
In the bundle method it is commonly proposed to solve the dual problem of ( 16), hence next we derive the dual of (16).Towards this end we collect the subgradients g i in the matrix G = (g 1 , . . ., g r ).After exchanging min and max by using strong duality the dual of ( 16) becomes max α 0 β 0 min y,w,v L(y, w, v, α, β).
Since ∇ w L = 0, ∇ v I L = 0, and ∇ y I L = 0 has to hold for all I ∈ J at the dual optimum, we get α ∈ ∆ r , β I ∈ ∆ t I and In consequence the dual of ( 16) simplifies to max This is a convex quadratic problem with r + I∈J t I variables and 1 + q simple equality constraints, asking that the respective block of variables adds up to one.Now instead of solving ( 16) within the bundle method directly, we solve its dual (18) to get the multipliers α and β and recover y using (17).

Our Bundle Method
So far we have sketched how to use our bundle method in order to obtain a solution y of ( 14), but actually we are interested in a solution X of ( 9).One can use the bundle B = {(y 1 , h 1 , g 1 , X 1 ), . . ., (y r , h r , g r , X r )}, which is updated in each iteration, in order to obtain a good approximate solution for X.In particular it follows from the convergence theory of the bundle method that under mild conditions converges to the optimal values of X and λ I of (9), see for example Robinson [30] for the general theory and Gaar [12] for the convergence in our particular setting.
We are now able to present our version of the bundle method.Note that there is no need of keeping y j in the bundle explicitly by computing and updating e in a proper way, so we drop y j from the bundle B. Algorithm 1 summarizes the main computational steps of our bundle method to get an approximate optimal solutions of ( 9) and ( 14).
Algorithm 1 Solving ( 9) and ( 14) for a given J 1: procedure Our Bundle(h, DI ) 2: y = y1 = 0 3: Evaluate h at y1 to get h1 = h(y1), a maximizer X1 and a subgradient g1 4: Set the bundle to B = {(h1, g1, X1)} 5: while "stopping condition is not satisfied" do 6: Solve the convex quadratic problem (18) to get α and β = (βI )I∈J 7: Determine y using (17) 8: Determine X and λ using (19) 9: Evaluate h at y to get h( y), a maximizer X and a subgradient g 10: Decide whether y becomes the new center y (serious step) or not (null step) 11: Update the bundle B and the bundle parameter µ 12: end while 13: return y, X, λ 14: end procedure The generic description of our bundle method in Algorithm 1 leaves some flexibility to the user.We will present implementation details in Section 6.3.

The Overall Algorithm
The goal of this paper is to get good bounds on the optimal Max-Cut value z mc , the stability number α(G) and the chromatic number χ(G) by including ESCs into the basic SDP relaxations (1), ( 2) and (3) in order to improve the bounds from the basic SDP relaxations.We will call bounds obtained in this way exact subgraph bounds (ESB).In other words ESBs are attained by solving (4), ( 5) and (6) or, in the generic form, by solving (9).
Up to now we have concentrated on the most subtle part of retrieving good ESBs, which consists in solving the SDP relaxation (9) with a given set J of ESCs.Our ultimate goal however is to reach ESBs where all ESCs of order k are (nearly) satisfied for small values of k like k 7.
We propose to reach this goal by proceeding iteratively.Starting with k = 3 in the Max-Cut case (as there are no violated ESCs of order 2) and k = 2 in the other cases we search for violated ESCs of order k and include only the most violated ESCs that we find into J.After solving the SDP (9), we follow an extreme strategy and remove any ESC that has become inactive.As we typically still find further badly violated ESCs this allows us a quick exploration of the entire space of ESCs.Once we do not find ESCs of order k with significant violation, we increase k and continue.We call each such iteration a cycle.
In each cycle so we keep some information, such as the current dual variables y i and the bundle B, appropriately modified to reflect possibly deleted and added new constraints.In particular we delete from all y i the positions corresponding to deleted ESCs, extend all y i with zeros for the newly added ESCs and deduce the update of all other variables.This choice allows us to reuse the bundle B. Our procedure to compute ESBs is sketched in Algorithm 2.

Algorithm 2 Computation of an exact subgraph bound
while "stopping condition is not satisfied" do 5: Get an approximate solution X of ( 9) and y of (14) 6: Update the ESB to the objective function value of y of (14) 7: Remove inactive ESCs from J 8: Include most violated ESCs of X with order k into J 9: if "not enough violated subgraphs found" then 10: end while 13: return ESB 14: end procedure The typical behavior over a set of cycles for one stable set instance can be seen in Figure 4.After only a few cycles with k = 2 we move to k = 3.Here it takes 16 cycles to reach a point with all ESCs nearly satisfied.The Figure clearly shows a continuing improvement of the ESB over the cycles.Note that the ESB computed in Algorithm 2 is indeed a valid bound, because any y is feasible for (14) and hence its dual objective function value is a valid bound on the primal optimal objective function value (9), which in turn is a valid bound on the optimal objective function value of the combinatorial optimization problem.Hence it is not necessary to solve ( 9) and ( 14) to optimality to obtain valid bounds.Of course we want to use our bundle method, Algorithm 1, in order to obtain the approximate solutions in line 5 of Algorithm 2.

Finding Violated Exact Subgraph Constraints
The key ingredients of Algorithm 2 are on the one hand Algorithm 1, which was detailed in Section 5, and on the other hand the update of the set J.
The crucial point in order to do so is to find violated ESCs.Let G I be a subgraph of oder k I of G and X * be the current solution of ( 9) and let U be an arbitrary k I ×k I matrix.Clearly CUT k I , STAB 2 (G I ) and COL(G I ) are bounded polytopes, hence the inner product of any element of these polytopes with U is contained in a certain interval.Thus finding I such that the inner product of U with the submatrix X * I of X * is minimum identifies a potentially violated subgraph.
This minimization may be recast as a quadratic assignment problem consisting of the data matrices X * and the matrix U embedded in an n × n matrix.We repeatedly use a local search heuristic for different fixed U in order to obtain potentially violated subgraphs.Then we compute the projection distances of X * I to the corresponding polytope for all these subgraphs G I and include those into J which have the largest projection distances and hence are violated most.
Possible choices for U make use of hyperplanes for the respective target polytope, but other choices are possible.In our computations we use a collection of different matrices for U , for example matrices that induce facets of the corresponding polytope (if their computation for a particular k I is possible easily, which is the case for k I 6), extreme copositive matrices with {0, 1, −1} entries and random matrices.For each cycle we use at most 50 different matrices U .

Details of the Bundle Implementation
We now briefly discuss some details our implementation of Algorithm 1 when used in line 5 of Algorithm 2. First of all one needs to decide on a stopping condition.Ideally we would stop, once a subgradient equal to zero is found.In our case, we either stop once the norm of the new subgradient is small enough (in the case of Max-Cut), or once the difference of the value of the function at the current center point and the value of the subgradient model of the function at the new trial point is smaller than some tolerance (as it is done in [4], the tolerance is 0.005 in our implementations) or once we reach a maximum number of iterations (30 in our implementations).The third condition is motivated by the fact that we typically will continue adding new violated ESCs, so there is no real need to get the exact minimum of (14).Note however that it is important to come close to the optimal solution, because otherwise the resulting X does not have a high enough precision in order to be useful for finding new violated subgraphs.
For updating the bundle we always add the new trial point to the bundle, but remove subgradients from the bundle that have become inactive.This extreme choice of updating the bundle led to the best performance in our computational experiments.In order to update the bundle parameter µ we use a modification of an update proposed by Kiwiel [21].We perform a serious step whenever the improvement of the objective function value of the new trial point is at leas a certain fraction of the expected improvement.This is a standard criterion, see for example [19].We solve the bundle master problem as a rotated second-order cone program (see [2] for more details) with MOSEK.

Bundle Approach versus Interior Point Methods
We start our computational investigation with a comparison of our bundle method with an interior point method in order to solve (9).In our overall Algorithm 2 presented in Section 6 this has to be done in each cycle, so we are highly interested in fast running times.
From a theoretical point of view it is clear which method will win this competition: Assume we include q = 1000 ESCs (so q = |J|) for subgraphs of order k I = 5 in (9) for the stable set problem.Then we have t I 2 5 = 32 stable set matrices that potentially span STAB 2 (G I ), and up to b I k I 2 + k I = 15 equality and one inequality constraint for each ESC.In total we have up to 32000 variables that have to fulfill up to 16000 constraints in (9)-additionally to the variables and constraints of the basic SDP relaxation (2).It is clear that the number of constraints will be a challenge for an interior point solver.In particular an interior point solver has to solve this SDP with a large number of constraints at once, whereas our bundle method in Algorithm 1 "only" has to solve the basic SDP relaxation and the bundle master problem over several iterations.Therefore, we expect the bundle method to be the clear winner in this competition and refrain from a large scale comparison.
Instead, we compare the two methods only on some instances to confirm our theoretical inspection.In Table 2 we list the results for one Max-Cut and one stable set instance, both are taken from the Erdős-Rényi model G(n, p).We vary the number of included ESCs for subgraphs of order 3, 4 and 5, so we solve (4) and ( 5) for different J.We choose J such that the total number of equality constraints induced by the convex hull formulation of the ESCs b ranges between 6000 and 15000.On the one hand we solve the instances with two interior point solvers, namely MOSEK and SDPT3 [33,34] and list the running times in seconds.On the other hand we use our bundle method.In our context we are mostly interested to improve the upper bounds quickly, so we do not run Algorithm 1 until we reach a minimizer, but stop after 30 iterations.We list the running time for the oracle, i.e. the sum of the solution times of the basic SDP relaxation, and the overall running times.Additionally we present how much % of the MOSEK running time the bundle method needs and how close the solution found by the bundle method is to solution of MOSEK in % (100% means the solutions coincide).
In Table 2 one sees that the running times decrease drastically if we use the bundle method compared to interior point solvers.For b ≈ 15000 it takes the bundle method only around 8% of the MOSEK running time to get as close as 95% to the optimal value, which is sufficient for our purposes.One sees that our bundle method scales much better for increasing |J|, so for an increasing number of ESCs.Furthermore MATLAB requires 12 Gigabyte of memory with interior point solvers for b = 15000, showing also memory limitations.
To summarize our small computational investigation confirms our intuition that the bundle method is much better suited for our purposes.
We want to point out that the number of bundle iterations can be increased in order to get closer to the optimum.For the larger instances in Table 2 this will still result in significantly shorter running times.
Note that the bundle method has another advantage: A warm start with the bundle B and the solution y of the previous iteration in line 5 of Algorithm 2 is possible.Since many ESCs remain the same in J the problem to solve in line 5 does not change too much and a warm start can be very beneficial.
As a last remark we want to draw the attention to the running times for the oracle in Table 2.For the stable set problem the oracle needs over half of the running time, whereas in the Max-Cut problem the oracle evaluation is much faster.This is due to the fact that the basic SDP relaxation is a simpler SDP for the Max-Cut problem.
In the following we present several computational results for obtained ESB by using the bundle method.Note that we refrain from comparing the running times of our bundle method with the running time of interior point methods, because interior point methods would reach their limit very soon.

The Stable Set and the Coloring Problem
In this section we will extend the computational results from [13] for the stable set and the coloring problem.The computational investigations show that (i) the ESB obtained by including ESCs of fixed order k I improve for increasing k I and (ii) after including several ESCs for subgraphs of order k I the maximum projection distance of the violated subgraphs found decreases drastically.
We extend these computational results by deriving one final ESB for several instances with Algorithm 2. We stop as soon as we have performed 50 cycles and only include subgraphs of order k 8.We add at most 100 ESCs in each cycle and warmstart the bundle with the information of the previous cycle.We already saw the typical behavior of the ESB over the cycles in Figure 4.
As a first structural easy class of graphs we consider two-dimensional torus graphs which are constructed as follows.For given d, the graph T d has d 2 vertices which we label by (i, j) for i, j ∈ {1, . . ., d}.The vertical edges join vertices with neighboring i indices (and j fixed), yielding edges {(i, j), (i + 1, j)} modulo d, and similarly the horizontal edges join vertices with i fixed {(i, j), (i, j + 1)} modulo d.So there is a total of n = d 2 vertices and m = 2n edges.It is not hard to verify that in case of odd d = 2t + 1, we get α(T d ) = t(2t + 1) and if d = 2t we have α(T d ) = 2t 2 .The even case is not interesting, as ϑ(T d ) = α(T d ).For d odd we summarize some computational results in Table 3.We observe that for these graphs our ESB is substantially better than ϑ(G) and we close the integer gap for all instances with n 121.When considering the running times observe that the majority of the running time (given in seconds in Table 3) is used for the oracle, because the SDP to evaluate ϑ(G) given in (2) with a slightly modified objective function is nontrivial.We tried several solver to solve this SDP, among them the interior point solver MOSEK [26], and solvers based on alternating direction method of multipliers as DADAL [7] and SDPNAL+ [35].Both these solvers show very good results on computing ϑ(G), but as soon as objective function slightly changes they do not perform well anymore.Hence it will be future research to develop an SDP solver dedicated to these kind of instances.Note that the running time in order to perform Algorithm 2 is not very high and in particular only increases mildly for larger instances.As a second class of problems we consider random near-r-regular graphs, which we generate as follows.We select a perfect matching on nr vertices and then we identify consecutive groups of r vertices into a single vertex.This yields a regular multigraph on n vertices.We remove loops and multiple edges resulting in a near-regular graph.In Tables 4 and 5 we provide results for random graphs.We compare near-regular graphs with random graphs from the Erdős-Rényi model where the density p is chosen so that the number of edges roughly matches those of the regular graphs.We compute our ESB and use a heuristic to compute large stable sets.In the results the gap between ϑ(G) and α(G) seems to be bigger for regular graphs, but we see in both cases that the ESB reduce the gap between ϑ(G) and the cardinality of the largest stable set found in a nontrivial way.Concerning running times we observe the same behavior as before.As a last experiment for the stable set problem in Table 6 we consider instances from the literature, taken mostly from the DIMACS challenge [9].On some instances there is hardly any improvement of the bound, while other instances are solved to optimality.It requires future research to get a better understanding for the fluctuation in quality on these instances, but for almost all instances the bound improves by at least one integer value.
The computation times for these instances range from 200 to 500 seconds for the smaller instances (n 125) to several hours for the biggest graphs.As in the instances before a faster oracle would improve the running times substantially.
Note that in our computations we aim for getting as good bounds as possible.If one wants to use the bounds in a branch-and-bound setting, a much more aggressive strategy with increasing k I faster and stopping as soon as we do not expect to reach the next integer value is favorable.Results for a selection of coloring instances from [27] are provided in Table 7.As in the case of the stable set problem we use Algorithm 2 to obtain ESBs.We include at most 100 ESCs in each cycle, only include ESCs for subgraphs of order k 8 and perform at most 25 cycles.The results are similar in quality to those for stable set from Table 6, so for the most instances we are able to obtain bounds, which are one integer value better than the original bounds from t * (G).The large running times are due to the difficult basic SDP relaxation (3).

The Max-Cut Problem
Finally we are ready to present computational results for the Max-Cut problem.It is well known that in the basic SDP relaxation of Max-Cut (1) all ESCs of order 3 can equivalently be represented by the metric polytope [23].Optimizing over it gives the exact solution to Max-Cut on graphs not contractible to K 5 , in particular on planar graphs.It is also well known that optimizing over the metric polytope may lead to rather weak relaxations for general graphs.In contrast, the simple SDP relaxation (1) provides an upper bound at most 14% above the optimal value of Max-Cut for graphs with nonnegative edge weights, see [14].
In our computational experiments with Max-Cut we noted that the number of ESCs necessary to insure that all ESCs for a given value k are satisfied can be quite large (see Section 3.2), even for small values of n, such as n = 100.We therefore simplify the ESC relaxation further.If a subgraph G I violates the ESC, then instead of asking that X I ∈ CU T k , we generate a single linear inequality separating X I from CU T k and include it instead of the ESC.This weakens the relaxation, but also reduces the computational effort, so that the total number of ESCs in the model may be quite large, and we can still compute the ESB.The computational effort is quite moderate, requiring no more than about 120 seconds for each of the instances.We first consider random graphs on n vertices from the Erdős-Rényi model G(n, p).Each edge is then assigned the weight 1 or −1 (each with probability 1/2).In Table 8 we report our computational results for n ∈ {100, 150} and p ∈ {0.1, 0.25, 0.5}.We compare the ESB with k = 3 (column labeled 3) to the ESB with k = 7 (column labeled 7).The column labeled 3 provides the deviation (in %) of the ESB with k = 3 from z mc .Thus if p is the value in the column labeled 3, then the ESB is equal to (1 + p/100)z mc .The column labeled 7 is to be understood in an analogous way for k = 7.In all cases we note a substantial gap reduction going from k = 3 to k = 7.The last column contains the number of ESCs at termination.It ranges from about 3000 for n = 100 to about 4500 for n = 150 and justifies our strategy to represent each ESC through a single cutting plane.
Next we consider graphs from the Beasley collection [3] with n = 250.Rendl, Rinaldi and Wiegele [28] used 10 of these instances in a branch-and-bound setting.The "hardest" instance 250-08 reported in [28] resulted in 4553 nodes in the branch-and-bound tree and took several days of computation time.All the other 9 instances from this collection resulted in branch-and-bound trees having between 17 and 223 nodes with computation times in the order of hours, see Table 6 from [28].We recomputed the root bound for all these instances and present our root gap in Table 9.We find it remarkable that our new bounding procedure is strong enough to prove optimality for all these instances right at the root node with the exception of problem 250-08.For this problem the gap at the root node was 2.19%.We recomputed the root bound in our setting and came up with a root gap of only 0.5%, thus reducing the gap by 75%.
As a final experiment we consider Max-Cut instances on Chimera graphs.This class of graphs has found increased interest in connection with quantum annealing, see [6] for further details.In Table 10 we provide computational results with such graphs on n = 512 vertices.We compute our ESB and also use a heuristic to find a good cut.It turns out that our bounding approach works nicely on these graphs, leading to provably optimal solutions in 2 out of 5 instances and the smallest possible positive gap (of 1) in the remaining cases.The computation times for each of these (big) instances range from 700 to 900 seconds, which we consider remarkable when dealing with more than 20000 ESCs.
Table 9. Max-Cut results for graphs from the OR library graph opt.cut # branch-and-bound nodes [28] root gap [28]  We conclude that for Max-Cut our ESB constitute a substantial improvement compared to the previously used strongest bounds based on SDP with triangle inequalities [28].These correspond to the column labeled 3 in Table 8.

Conclusions and Future Work
Summarizing, we offer the following conclusions from the computational results.Our computational approach based on the partial Lagrangian dual is very efficient in handling also a large number of ESCs.The dual function evaluation separates the SDP part from the ESCs and therefore opens the way for large-scale computations.The minimization of the dual function is carried out as a convex quadratic optimization problem without any SDP constraints, and therefore is also suitable for a large number of ESCs.
Our computational results for stable set and coloring confirm the theoretical hardness results for these problems.Including ESCs of rather small size (k 8) yields a noticeable improvement of the bounds.
The limiting factor for stable set instances is the solution time of the oracle.Hence it is desirable to have a fast solver for these kind of instances.
On the practical side we consider the cutting plane weakening of the ESCs for Max-Cut a promising new way to tighten bounds for this problem.
It will be a future project to explore these bounds in a branch-and-bound setting in order to solve Max-Cut, stable set and coloring instances to optimality.

Fig. 1 .Fig. 2 .Fig. 3 .
Fig. 1.Histogram of δ I mc for all violated subgraphs GI of order k ∈ {3, 4, 5} we do not have to include all k I × k I equality constraints into the SDP.Let b I be the number of equality constraints we have to include.Note that b I = k I 2 , b I = k I 2 + k I − m I and b I = k I 2 − m I for the Max-Cut, stable set and coloring problem respectively, if m I denotes the number of edges of G I .This is because in the case of the stable set problem we also have to include equations for the entries of the main diagonal contrary to Max-Cut and the coloring problem.Then we define a linear map M I : R b I → S k I such that the adjoint operator M I : S k I → R b I extracts the b I positions, for which we have to include the equality constraints, into a vector.So we can rephrase (8) equivalently as

2 +
It will be notationally convenient to partition the matrix G into blocks of rows corresponding to the subsets I ∈ J, namely G = (G I ) I∈J where each G I has r columns and b I rows.Furthermore we make the subgradient model and maximum term constraints more compact by reformulating them to w1 h1 − e + I∈J G I (y I − y I ) and v I 1 D I (y I ).We denote by α ∈ R r the dual variables to the subgradient model constraints and with β I ∈ R t I the dual variables of the constraints involving v I for the maximum terms.Furthermore we define β = (β I ) I∈J as the collection of all β I .Hence we obtain the Lagrangian function L(y, w, v, α, β) = w + I∈J α, h1 − e − w1 + I∈J α, G I (y I − y I ) + I∈J β I , D I (y I ) − v I 1 .

Fig. 4 .
Fig. 4. Progress of the ESB over 50 cycles for one instance of Table 4 the dimension is clear from the context we may omit the index and write 1 and ∆.Furthermore let N = {1, 2, . . ., n}.A graph G on n vertices has vertex set N and edge set E. The complement graph G of a graph G has the same vertex set N and contains an edge {i, j} ⊆ N if and only if {i, j} ∈ E. S n is the set of n-dimensional symmetric matrices.A spectrahedron is a set that is obtained as the intersection of the cone of positive semidefinite matrices with some linear affine subspace.

Table 1 .
Percentage of violated subgraphs of order k for one random graph

Table 2 .
Running times for one Max-Cut and one stable set instance with different sets of ESCs, where the graphs of order n = 100 are from the Erdős-Rényi model

Table 3 .
Stable set results for torus graphs

Table 4 .
Stable set results for near-regular graphs

Table 5 .
Stable set results for graphs from the Erdős-Rényi model G(n, p)

Table 8 .
Max-Cut results for graphs from the Erdős-Rényi model G(n, p)

Table 10 .
Max-Cut results for Chimera graphs with n = 512