Strong SDP based bounds on the cutwidth of a graph

Given a linear ordering of the vertices of a graph, the cutwidth of a vertex $v$ with respect to this ordering is the number of edges from any vertex before $v$ (including $v$) to any vertex after $v$ in this ordering. The cutwidth of an ordering is the maximum cutwidth of any vertex with respect to this ordering. We are interested in finding the cutwidth of a graph, that is, the minimum cutwidth over all orderings, which is an NP-hard problem. In order to approximate the cutwidth of a given graph, we present a semidefinite relaxation. We identify several classes of valid inequalities and equalities that we use to strengthen the semidefinite relaxation. These classes are on the one hand the well-known 3-dicycle equations and the triangle inequalities and on the other hand we obtain inequalities from the squared linear ordering polytope and via lifting the linear ordering polytope. The solution of the semidefinite program serves to obtain a lower bound and also to construct a feasible solution and thereby having an upper bound on the cutwidth. In order to evaluate the quality of our bounds, we perform numerical experiments on graphs of different sizes and densities. It turns out that we produce high quality bounds for graphs of medium size independent of their density in reasonable time. Compared to that, obtaining bounds for dense instances of the same quality is out of reach for solvers using integer linear programming techniques.


Introduction
Several graph parameters are determined by finding an arrangement of the vertices of a graph on a straight line having a certain objective in mind. Depending on the objective, these parameters are, for instance, the treewidth, the pathwidth, the bandwidth or the cutwidth of a graph. Computing these graph parameters is necessary for several applications, e.g., when a certain layout has to be found (VLSI design [27]), in automatic graph drawing [9] or in many versions of network problems arising in energy or logistics. All these applications ask for algorithms to practically compute these parameters. However, this leads to NP-hard optimization problems and therefore algorithms for approximating these parameters in terms of lower and upper bounds are required. The parameter that is of our interest in this work is the cutwidth of a graph.
Definitions. The minimum cutwidth problem (MCP) can be defined as follows. We consider an undirected graph G = (V, E) with vertex set V and edge set E and assume without loss of generality that the set of vertices V is V = {1, . . . , n}. Furthermore, the set of all permutations of {1, . . . , n} is denoted by Π n . In any permutation π ∈ Π n of the vertices of G, the position of a vertex v ∈ V in π is given by π(v). Note that any linear ordering of the vertices V of G can be represented by a permutation π ∈ Π n .
The cutwidth CW π (v) of a vertex v with respect to the permutation π is the number of edges {u, w} ∈ E such that π(u) ≤ π(v) < π(w) holds, i.e., the cutwidth of v in π is the number of edges from any vertex before v (including v) to any vertex after v (not including v) in a linear ordering of the vertices according to the permutation π. Furthermore, the cutwidth CW π (G) of a graph G with respect to π is the maximum cutwidth of any vertex with respect to the permutation π, so CW π (G) = max v∈V CW π (v) holds. Finally, the cutwidth CW (G) of a graph G is the minimum cutwidth of G with respect to π over all possible permutations π, i.e., CW (G) = min π∈Πn CW π (G).
An obvious lower bound on the cutwidth is given by where ∆(G) is the maximum degree of any vertex in V . Indeed, let us denote by v max a vertex of G with degree ∆(G). Then, for every linear ordering π, every vertex in the neighborhood of v max is counted either in CW π (v max ) or in CW π (u max ), where u max is the vertex directly preceding v max in π. It follows that either CW π (v max ) or CW π (u max ) has to be greater or equal to ∆(G) 2 , and due to the integrality of CW (G) we obtain the lower bound ∆(G)+1

New bounds for the minimum cutwidth problem
The aim of this section is to introduce a new SDP relaxation for the MCP and to present ways to strengthen this relaxation.

Our new basic SDP relaxation
We follow the approach of Buchheim, Wiegele and Zheng [4] for the quadratic linear ordering problem in order to derive an SDP formulation of the MCP. This is a promising endeavor, as the feasible region of both the quadratic linear ordering problem and the MCP consist of permutations.
Towards this end, let A = (a ij ) 1≤i,j≤n be the adjacency matrix of the graph G, i.e., a ij = 1 if and only if the edge {i, j} is in the set of edges E of the graph G. To represents a permutation π, the authors of [4] introduce a binary variable χ π ij for any 1 ≤ i < j ≤ n that indicates whether π(i) < π(j) holds, so in total they have a vector χ π consisting of n 2 binary variables. They define the linear ordering polytope LO(n) as the convex hull of all vectors χ π , formally which is a subset of R ( n 2 ) . As a consequence, the elements of the set LO(n) ∩ {0, 1} ( n 2 ) are exactly all vectors representing permutations of {1, . . . , n}.
representing a permutation π ∈ Π n . By definition, the binary variable x ij is equal to 1 if and only if i is before j in π. Then, it can be checked that for any vertex v of G, holds. The first four terms of this expression count the number of edges from any vertex before v to any vertex after v in the permutation, both not including v. These four terms are necessary, as only one of the variables x uv (if u < v) and x vu (if u > v), and one of the variables x wv (if w < v) and x vw (if w > v) exist. The last two terms of (1) count the edges from v to any vertex after v. By expanding (1), grouping the constant, linear and quadratic terms together, and by combining three times two sums as one, we can rewrite (1) as and as a result, the MCP can be written as This is an integer program with a linear objective function and quadratic constraints, which is perfectly suited for deriving an SDP relaxation. To do so, we introduce the matrix variable X = (X ij,k ) 1≤i<j≤n 1≤k< ≤n . In particular, X ij,k represents the product x ij x k .
Then the following SDP is a relaxation of the MCP.
Note that x = (x ij ) 1≤i<j≤n is a vector of dimension n 2 andX is a square matrix with ( n 2 + 1) columns and rows. Thus, the SDP (4) has a matrix variable of dimension ( n 2 +1) and n+ n 2 +1 constraints. The n constraints (4b) make sure that α is at least the cutwidth of each vertex. The n 2 + 1 constraints (4c) together with the constraints (4d) and (4e) represent the relaxation of X − xx = 0 to X − xx 0 and taking the Schur complement. Due to the fact that x is binary, X − xx = 0 also implies that (4c) has to hold.
Furthermore, note that adding the non-convex constraint rank(X) = 1 to the SDP relaxation (4), the optimal objective function value becomes CW (G). This is also the case if we add integer conditions for X.

Strengthening the SDP relaxation
Next, we investigate several options to improve the SDP relaxation (4) of the MCP.

3-dicycle equations
In the basic SDP relaxation (4), no specific information about the fact that x should represent a permutation is used. One possible way to include such information is to model the transitivity by so-called 3-dicycle equations, as it is done by Buchheim, Wiegele and Zheng [4]. These 3-dicycle equations make sure that if i is before j and j is before k, then i must be before k. For a vector x = (x ij ) 1≤i<j≤n , they can be written as so we can include as additional constraints into our SDP relaxation (4) of the MCP. If we add all possible constraints of the form (6), we include n 3 equality constraints. Observe, that the authors of [4] show that if the variables are binary, the expression on the left hand-side of (5) is always non-negative. So in this case, one needs to consider only the inequality ≤ 0. However, this does not hold anymore for our SDP relaxation, as the variables are not necessarily binary.

Triangle inequalities
Another way to strengthen the SDP relaxation (4) is adding the triangle inequalities The inequalities (7a) and (7b) were introduced by Lovász and Schrijver [24] and the constraints (7c), (7d) and (7e) originate from Gruber and Rendl [16]. As a consequence, (7a), (7b) and (7c) each yield n It can be checked easily that all the inequalities of (7) are satisfied for each matrix X = xx for any x ∈ LO(n)∩{0, 1} ( n 2 ) , so for any vector x that represents a permutation. Thus, the inequalities (7) are valid inequalities for the SDP relaxation (4) for the MCP. Next, we give an alternative reasoning that the inequalities (7) are satisfied.

Inequalities obtained from the squared linear ordering polytope
Adams, Anjos, Rendl and Wiegele [1] introduced so-called exact subgraph constraints (ESC), which were later computationally exploited for the stable set problem, the Max-Cut problem and the coloring problem by Gaar and Rendl [12,13].
For the stable set problem the ESCs are defined in the following way. Let G = (V, E) be a graph with vertex set V = {1, 2, . . . , n}. The squared stable set polytope ST AB 2 (G) of G is defined as and the ESCs for the stable set problem (ESCSS) are where I ⊆ V , G I is the induced subgraph of G on the vertices I, X is the matrix variable of an SDP relaxation of the stable set problem, and X I is the submatrix of X that corresponds to G I . As detailed in [14], the ESCSS are equivalent to This implies that only the squared stable set polytope for a graph with k vertices and no edges needs to be considered.
Furthermore, the authors of [14] describe that the ESCSS can be represented by inequalities for X I , and therefore as inequalities for X. In [11], these inequalities that represent the ESCSSs are stated explicitly for subgraphs with 2 and 3 vertices. In fact, the triangle inequalities (7) are exactly the inequalities representing the ESCSS for subgraphs of order 2 ((7a), (7b) and (7c)) and for subgraphs of order 3 ((7d) and (7e)). When deriving the inequalities that represent the ESCSS for G 0 k , any binary vector of dimension k is feasible for the corresponding stable set problem in G 0 k . As a consequence, For the MCP, only specific (and not all) binary vectors of dimension n 2 induce a permutation of the n vertices. Thus, it makes sense to deduce the ESCs specifically for the linear ordering problem as they might be more restrictive. Towards this end, we introduce the quadratic linear ordering polytope which is a subset of R ( n 2 )×( n 2 ) . From our considerations above, it follows that The vertices of the polytope LO 2 (n) can easily be enumerated. With the help of the software PORTA [5] it is possible to determine the inequalities that represent LO 2 (n) for small values of n. In particular, any matrix X is in LO 2 (n) for n = 3 if and only if it is symmetric and all facet defining inequalities and equations are satisfied. Note, that (8a), (8c) and (8d) each yield n 3 potential inequalities and (8b) offers the option of including 4 n 3 constraints. It is easy to see that (8d) coincides with the 3-dicycle equations already considered as (6). Furthermore, (8a), (8b) and (8c) are a subset of the triangle inequalities (7a), (7b) and (7c), respectively. It turns out that the following holds. (8), then X also fulfills all inequalities of (7) whenever |{i, j, k, , u, v}| ≤ 3 holds.

Lemma 2.1. If a symmetric X satisfies
Proof. Assume a symmetric X satisfies (8).
Note that Lemma 2.1 is indeed no surprise: With the intuition of the ESCs for the stable set problem, the triangle inequalities (7) (which represent the membership to ST AB(G 0 3 )) must be satisfied for any 3 × 3 submatrix of X, as the only condition is that the corresponding submatrix must be formed by binary vectors. Instead, the constraints (8) (which represent the membership to LO 2 (3)) are only valid for 3 × 3 submatrices of X whose indices that correspond to the three pairs ij, ik, jk for any 1 ≤ i < j < k ≤ n. Thus, the inequalities (8) capture more structure, but are valid for fewer submatrices.
With the help of LO 2 (4) we were able to find the next valid inequalities, which are (to the best of the knowledge of the authors) not known as valid inequalities for the linear ordering problem so far.

Lemma 2.2. The inequalities
are valid inequalities for the SDP relaxation (4) for the MCP.
Proof. In order to prove this lemma, it is enough to show that for all i < j < k < and for all X = χ π χ π for π ∈ Π n the inequality (11) is satisfied. To do so, we distinguish two cases for fixed i < j < k < and fixed π.
As a result, in any case (11) is satisfied, which finishes the proof.

Liftings of inequalities from the linear ordering polytope
Buchheim, Wiegele and Zheng [4] derived their quadratic 3-dicycle equations (5) as an alternative to the linear 3-dicycle inequalities One could also use the standard reformulation linearization technique (RLT), of which the foundations were laid by Adams and Sherali [2]. Using this approach, we can multiply the inequalities (13) by x uv and (1−x uv ) for every pair (u, v) with u < v, and then replace products of the form x ij x k by X ij,k . In this way we obtain the set of valid inequalities which yields 4 n 3 n 2 potential additional constraints. However, some of them are already implied by the inequalities previously introduced. Indeed, (7a) and (6), (7b) and (6), • (14c) for (u, v) = (i, j) and (u, v) = (j, k) is implied by (7b) and (6), • (14d) for (u, v) = (i, j) and (u, v) = (j, k) is implied by (7b) and (7c), and • (14d) for (u, v) = (i, k) is implied by (7c) and (6).
Thus, at least one of u and v needs to be different to i, j and k such that an inequality of (14) has the potential to bring new information.

Algorithms
In this section we describe in detail our algorithm that utilizes our new SDP relaxation for the MCP and its strengthenings derived in the last section. Furthermore, we present a new upper bound obtained by a heuristic utilizing the optimizer of the SDP relaxation.

Algorithm for computing our new lower bound
Adding all the previously described inequality and equality constraints at once to the basic SDP relaxation (4) would render it too computationally expensive to solve. Therefore, a cutting-plane approach is used to obtain a tight lower bound for the MCP.

Separating constraints
The total number of possible inequality and equality constraints to add to the basic SDP relaxation (4) is O(n 6 ). While it is possible to exhaustively enumerate them and keep the ones with the largest violation, it would take a long time and does not guarantee the best tightening of the relaxation. A heuristic to find violated constraints is thus preferred.
The algorithm used to find a single violated constraint for a given solutionX of the SDP relaxation (4) is a simulated annealing heuristic, in which the current solution is represented by a tuple of indices ((i, j, k, ), (u, v), (q, w)), with i < j < k < , u < v and q < w, coupled with the constraint type. The possible constraint types are the 3-dicycle equations (6); the triangle inequalities (7a), (7b), (7c), (7d), (7e); the inequality from the squared linear ordering polytope of order four (10), and the inequalities obtained from lifting the linear ordering polytope (14a), (14c), (14b), (14d) and are given in the array constraintsToTest. Note that the tuple of indices is defined in such a way that the required indices of all possible constraints can be extracted from it. For example, the three pairs of indices of (7e) can be extracted as (i, j), (u, v) and (q, w) from the current solution.
Whenever a current solution is given, a neighbor solutions are obtained with ran-dom_neighbor_indices() by randomly replacing one index from the tuple of indices of the current solution before ordering it again. The violation of a solution (i.e., a fixed type of constraint for a fixed tuple of indices) for the solutionX of the SDP relaxation (4) can be computed via violation(), and is positive if the inequality or equality constraint is violated. The pseudocode of this simulated annealing heuristic can be found in Algorithm 1.

Algorithm 1: Simulated annealing to separate constraints
Parameters: T init , f t , maxIter, maxLenPlateau Input: solutionX of the SDP relaxation (4), array of types of constraints constraintsToTest

Outline of the overall algorithm
After detailing how to find violated inequality and equality constraints with simulated annealing in Algorithm 1, we are now able to describe our main algorithm to derive a lower bound for the MCP.
Our algorithm starts by solving the basic SDP relaxation (4), providing a first lower bound α init and the associated solutionX init . Some constraints violated by the current solutionX init , are then determined with Algorithm 1. In particular, at each iteration, Algorithm 1 is run 2numCuts times, and the numCuts most violated constraints are added. These violated constraints are then added to the SDP, which is solved again to obtain a new improved lower bound α and a new current solutionX. These steps are repeated for a fixed number of iterations maxIter, or until the improvement of the bound does not reach a certain threshold.
To reduce the computational effort, at each iteration the constraints that seem to be not necessary for obtaining the bound with the SDP are removed. To do so, the mean γ mean of all absolute values of all dual values γ associated with each added inequality is computed. The constraints that are associated with a dual value |γ| < 0.01γ mean are then removed. The pseudocode of our algorithm can be found in Algorithm 2.
The upadate of constraintsToTest is done in the following way to improve the efficiency. Empirically, the constraints added in the first two iterations are mostly 3dicycle equations (6). This is consistent with the fact that they are the ones adding the most to the structure of the problem. From the second or third iteration on (depending mostly on the size of the instances), the triangle inequalities (7) represent the most part of the added cuts. Thus, in the first two iterations of our algorithm only violated 3dicycle equations are added. Triangle inequalities are then added to our pool of potential constraints constraintsToTest for the third and fourth iterations, and all constraints are considered for the remaining iterations of the algorithm.

Our new heuristic for computing an upper bound
We can utilize the SDP relaxation (4) of the MCP not only to derive lower bounds as seen in Section 3.1, but also to obtain upper bounds on the cutwidth. In particular, our aim is to derive a feasible solution of the MCP, and thus an upper bound, from the SDP solutionX returned by Algorithm 2. For that purpose, the first column x ofX is considered, without its first element (equal to 1 by definition ofX). So x = (x ij ) 1≤i<j≤n is a vector of size n 2 , with entries between 0 and 1. For any i ∈ V = {1, 2, . . . , n}, we compute the relative position p i of vertex i as For any vector x ∈ LO(n) ∩ {0, 1} ( n 2 ) encoding a linear ordering π, p i is equal to the number of vertices before the vertex i in π, i.e., to the position of i in the linear ordering. In the general case of a vector x obtained from the SDP relaxation (4) (with or without Algorithm 2: Cutting-plane algorithm to obtain a lower bound for the MCP Parameters: maxIter, improvementMin, numCuts, minViolation Input: adjacency matrix A, array of types of constraints constraintsToTest 1 solve the SDP relaxation (4) to obtainX init and the lower bound α init 2X ←−X init 3 α ←− α init 4 iter ←− 0 5 improvement ←− +∞ 6 addedCuts = ∅ 7 while iter < maxIter and improvement > improvementMin do 8 iter ←− iter + 1 added cuts), and not necessarily feasible for the MCP, the p i are first sorted in ascending order: let (i 1 , . . . , i n ) be an ordering of V such that p i 1 ≤ p i 2 ≤ . . . ≤ p in . A linear ordering π can then be retrieved by assigning to each vertex i the position of p i in the sorted list, so π(i 1 ) = 1, π(i 2 ) = 2, . . . , π(i n ) = n. The pseudocode of this algorithm can be found in Algorithm 3. This provides a feasible linear ordering for the MCP, which we obtained directly from a feasible SDP solution. Thus, the cutwidth of this ordering is an upper bound for the MCP. This upper bound is then improved by running a simulated annealing heuristic, that uses the obtained feasible solution as starting solution.

Computational experiments
In this section we evaluate the quality of the bounds obtained by our algorithms through numerical results. We compare the bounds to those obtained by modeling the problem as MILP.

Computational setup
Our algorithm was implemented in Python, using the graph library networkx. The SDP relaxations were solved using the MOSEK Optimizer API [26] version 10. We use CPLEX [20] version 22.1 to solve the MILP in order to compare bounds and computation times of our approach to an MILP approach from the literature. We ran the tests on an Intel Xeon W-2125 CPU @ 4.00GHz with 4 Cores / 8 vP with 256 GB RAM. The code and the instances are available as ancillary files from the arXiv page of this paper.

Instances
The benchmark instances commonly used in the literature for linear ordering problems are extremely sparse. MILP based approaches are thus very efficient on them, as demonstrated by the results obtained by David Coudert in [7]. However, it is expected that denser graphs are much more challenging to solve by such approaches.
Among the graphs considered in [7], the densest instance has 22 vertices and 49 edges, i.e., a density of approximately 21%. Most of the instances have densities much below 20% or even below 10%. As we want to evaluate our bounds for graphs with various densities, and in particular for dense graphs, we generated instances according to well established models: Erdős-Rényi graphs and random geometric graphs.
A random geometric graph U (n, d) is generated by first placing n vertices uniformly at random in the unit cube. Two vertices are then joined by an edge if and only if the Euclidian distance between them is at most d, see for example [21]. Our set created this way is composed of 45 graphs with n ∈ {20, 30, 40, 50}, the distances d are chosen from the set {0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}. Figure 1 shows the evolution of the bound computed by our algorithm over seven iterations considering different sets of constraints on the example of an Erdős-Rényi graph on n = 20 vertices. Note that in the initial iteration no constraints are added and in the first iteration always only the 3-dicycle equations (6) are considered. The circles in Figure 1 show that the bound does not move anymore from iteration three on if only considering constraints (6). Adding the triangle inequalities (7) to the pool of possible constraints leads to a remarkable increase of the bound, cf. the triangles in the plot. And offering additionally the inequalities obtained from lifting the linear ordering polytope (14) and the inequalities from the squared linear ordering polytope of order four (10) further pushes the bound above with neglectable extra cost, cf. the stars in Figure 1. The plots are similar for all the instances we consider. Hence, these experiments confirm our choice of using all the constraints (6), (7), (10) and (14) as potential strengthenings within our algorithm.

Numerical results
In Tables 1-4 we give the numerical results for the above mentioned instances, comparing our SDP bounds with those obtained when using an MILP solver. The column label n refers to the number of vertices of the graph, p and d relates to the edges of the graph as specified in Section 4.2 above.
In Tables 1 and 3 columns "LB init" list the lower bounds obtained when solving the initial basic SDP relaxation without any cutting planes, while in column "LB final" the lower bound we finally obtained is displayed. "UB" is the upper bound obtained through our heuristic and "gap final" is computed as (UB − LB)/UB, where LB is the final lower bound. In the column with the "time" label, we report the total time for obtaining our lower bound, and how this total time is split into solving the SDPs and separating the violated inequality and equality constraints. Moreover, the time for computing the upper bound is given. The last column of these two tables indicate the number of cutting planes added when the algorithm stops. Tables 2 and 4 give the details when using CPLEX to obtain an (approximate) solution.
In the tables we see that adding the valid inequalities and equalities significantly improves the initial bound, in Table 1 the value of the bound increases by roughly 30%, for the random geometric graphs (Table 3) it is around 20 to 30%. Also the upper bounds we obtained are of good quality, overall the gap between our final lower bounds and the upper bounds is from 30% to 46% for the Erdős-Rényi graphs and between 27% and 50% for the random geometric graphs.
The time spent to compute the lower bounds ranges from a few seconds for the instances with 20 vertices up to 115 minutes for instances with 50 vertices. For the small instances, solving the SDPs can be done rather quickly and therefore, the time spent in the separation takes the bigger ratio. For larger instances, however, separating violated inequalities and equalities typically uses less than 20% of the overall time. Computing the upper bounds is done within less than 40 seconds.
We now turn our attention to comparing our results with using the MILP solver CPLEX, see Tables 2 and 4. As all our SDP based results are obtained within 2 hours, we set this as time limit for CPLEX. On small and sparse instances CPLEX performs very well. However, for the Erdős-Rényi graph with n = 40 and 70% density we are left with a gap of 45% after two hours, compared to the gap of 34% that we obtain with our algorithm in less than 30 minutes. The random geometric graphs with d equal to 0.3 or 0.4 can be solved for graphs with up to 40 vertices, for dense instances with at least 40 vertices a gap of more than 35% remains after two hours run time. In general, the density has a huge impact on the runtime for the MILP solver. For example, it increases from 3.45 seconds for p = 0.3 to 1206.53 seconds for p = 0.9 for Erdős-Rényi graphs on 20 vertices. Our SDP based bounds do not show significant differences for sparse and dense instances for the Erdős-Rényi graphs, and only a minor increase of the runtime for the random geometric graphs.
Overall, the MILP approach is clearly outperformed by our method for dense graphs but also for graphs having 40 or 50 vertices, independent of their density.

Conclusion and outlook
We presented an SDP based approach to compute lower and upper bounds on the cutwidth. In order to obtain tight bounds, we derive several classes of valid inequalities and equalities and a heuristic for separating those inequalities and equalities and add them iteratively in a cutting-plane fashion to the SDP relaxation. The solution obtained from the SDP relaxation also serves to compute an upper bound on the cutwidth. Our experiments show that we obtain high quality bounds in reasonable time. In particular, our method is by far not as sensitive to varying densities as MILP approaches.
As the number of vertices gets larger, solving the SDPs becomes the time consuming part and thus the bottleneck. This is due to the increasing size of the matrix, but even more due to the huge number of constraints. Therefore, in our future research we will investigate using alternating direction methods of multipliers (ADMM) instead of an interior-point method, as ADMM proved to be successful in solving SDPs with many constraints, see for example [8]. Also including our bounds in a branch-and-bound framework and thereby having an exact solution method is a promising future research direction.
In our experiments we observe that we were not able to push the gap significantly below 0.30. We want to further investigate this to find theoretical evidence of this behavior. And finally, we want to apply our approach to computing related graph parameters, i.e., those that can also be formulated in terms of orderings of the vertices, like the treewidth or the pathwidth of a graph.