Edge expansion of a graph: SDP-based computational strategies

Computing the edge expansion of a graph is a famously hard combinatorial problem for which there have been many approximation studies. We present two variants of exact algorithms using semidefinite programming (SDP) to compute this constant for any graph. The first variant uses the SDP relaxation first to reduce the search space considerably. One implementation of this variant applies then an SDP-based branch-and-bound algorithm, along with heuristic search. The second implementation transforms the problem into an instance of a max-cut problem and solves this using an SDP-based state-of-the-art solver. Our second variant to compute the edge expansion uses Dinkelbach's algorithm for fractional programming. This is, we have to solve a parametrized optimization problem and again we use semidefinite programming to obtain solutions of the parametrized problems. Numerical results demonstrate that with our algorithms one can compute the edge expansion on graphs up to 400 vertices in a routine way, including instances where standard branch-and-cut solvers fail. To the best of our knowledge, these are the first SDP-based solvers for computing the edge expansion of a graph.


Introduction
Let G = (V, E) be a simple connected graph on n ≥ 3 vertices and with m edges.A cut (S, S ′ ), for any ∅ ̸ = S ⊂ V and S ′ = V \ S, in G is a partition of its vertices, and the cut-set ∂S is the edges between S and S ′ .The (unweighted) edge expansion, also called the Cheeger constant or isoperimetric number or sparsest cut, of G is a ratio that measures the relative number of edges across any vertex partition.It is defined as where ∂S = {(i, j) ∈ E : i ∈ S, j ∈ S ′ } is the cut-set associated with any vertex subset S ⊂ V , and S ′ = V \ S.This constant is positive if and only if the graph is connected, and the exact value tells us that the number of edges across any cut in G is at least h(G) times the number of vertices in the smaller partition.A weighted definition of edge expansion called the conductance of a graph, is where vol(S) = v∈S deg(v), and the second equality is due to vol(S) + vol(S ′ ) = 2m.Edge expansions arise in the study of expander graphs, for which there is a rich body of literature with applications in network science, coding theory, cryptography, complexity theory, cf.[13,21,39].A graph with h(G) ≥ c, for some constant c > 0, is called a c-expander.A graph with h(G) < 1 is said to have a bottleneck since there are not too many edges across it.A threshold for good expansion properties is having h(G) ≥ 1, which is desirable in many of the above applications.The famous Mihail-Vazirani conjecture [11,32] in polyhedral combinatorics claims that the graph (1-skeleton) of any 0/1-polytope has edge expansion at least 1.This has been proven to be true for several combinatorial polytopes [22,32] and bases-exchange graphs of matroids [1], and a weaker form was established recently for random 0/1-polytopes [28].
Lattice polytopes were constructed in [14] with the property that in every dimension their graphs lie on the threshold of being good expanders (i.e., h(G) = 1).
Computing the edge expansion is related to the uniform sparsest cut problem which asks for computing a cut in the graph with the smallest sparsity, where sparsity is defined as the ratio of the size of the cut to the product of the sizes of the two partitions, , which implies that the edge expansion problem is related to the sparsest cut problem up to a constant factor of 2. In particular, any cut (S, S ′ ) that is α-approx for ϕ(G) (resp.h(G)) is a 2α-approx for h(G) (resp.ϕ(G)), because There are polynomial reductions between h(G), h vol (G) and ϕ(G) and they are all NPhard to compute [27], in contrast to the minimum-cut of a graph which can be computed in polynomial time.Hence, almost all of the literature on edge expansion is devoted to finding good theoretical bounds.These are generally associated with the eigenvalues of the Laplacian matrix of the graph and form the basis for the field of spectral graph theory (see the monograph [9]).There have also been many approximation studies on this topic [3,27,34,37], and semidefinite optimization (SDP) has been a popular tool in this regard.The best-known approximation for ϕ(G) is the famous O( √ log n) factor by Arora et al. [3] which improved upon the earlier O(log n)-approximation [27].The analysis is based on an SDP relaxation with triangle inequalities and uses metric embeddings and concentration of measure results.Meira and Miyazawa [31] developed a branch-and-bound algorithm for computing ϕ(G) using SDP relaxations and SDPbased heuristics.Recall that ϕ(G) is related to h(G) in the approximate sense (up to a factor 2) but not in the exact sense.To the best of our knowledge, there is no exact solution algorithm for h(G).

Contribution and outline
We adopt mathematical programming approaches for numerical computation of h(G).All our approaches make use of tight bounds obtained via semidefinite programming.The first algorithm works in two phases.In the first phase, we split the problem into subproblems and by computing lower and upper bounds for these subproblems, we can exclude a significant part of the search space.In the second phase, we either solve the remaining subproblems to optimality or until a subproblem can be pruned due to the bounds.For the second phase, we develop two versions.The first version implements a tailored branch-and-bound algorithm, in the second version we transform the subproblem into an instance of a max-cut problem and compute the maximum cut using an SDP-based solver.
The second algorithm we implement uses the idea of Dinkelbach's algorithm to solve fractional optimization problems.The main concept of this algorithm is to iteratively solve linearly constrained binary quadratic programs.We solve these problems again by transforming them into instances of max-cut and using an SDP-based solver to compute the maximum cut.
We perform numerical experiments on different types of instances which demonstrate the effectiveness of our approaches.To the best of our knowledge, no other algorithms are capable of computing the edge expansion for graphs with a few hundred vertices.
The rest of the paper is structured as follows.In § 2 we formulate the problem as a mixed-binary quadratic program and present an SDP relaxation.§ 3 investigates a related problem, namely the k-bisection problem.We introduce a branch-and-bound algorithm and describe all the ingredients in this section.The first algorithm (relying on the k-bisection problem) for computing h(G) is introduced in § 4, and another algorithm (following Dinkelbach's idea) in § 5.The performance of all algorithms is demonstrated in § 6, followed by conclusions in § 7.
Notation The set of n × n real symmetric matrices is denoted by S n .The positive semidefiniteness condition for X ∈ S n is written as X ⪰ 0. The trace of X is written as tr(X) and defined as the sum of its diagonal elements.The trace inner product for X, Y ∈ S n is defined as ⟨X, Y ⟩ = tr(XY ) and the operator diag(X) returns the main diagonal of matrix X as a vector.The vector of all ones is e and the matrix of all ones is J = ee ⊤ .
For a n-vertex graph G = (V, E), the minimum and maximum vertex degrees are δ(G) and ∆(G), the adjacency matrix is a binary matrix A ∈ S n having A ij = 1 if and only if (i, j) ∈ E, and the degree matrix is a n × n positive diagonal matrix D having D ii equal to the degree of vertex i ∈ V .The Laplacian matrix is L = D − A, and thus has its nonzero entries as

Formulations and SDP relaxations
To write an algebraic optimization formulation for cut problems in graphs, we represent a cut (S, S ′ ) in G by its incidence vector χ S ∈ {0, 1} n which has χ S i = 1 if and only if i ∈ S. The cut function is the size of a cut-set, also called the value of the cut, and is equal to Any binary vector x ∈ {0, 1} n represents a cut in this graph.Denote the set of all cuts with S containing at least one vertex and at most half of the vertices by Using the common expression x ⊤ Lx for the cut function, the edge expansion problem is x,y y : The last formulation is a mixed-binary quadratically constrained problem (MIQCP).Standard branch-and-cut solvers may require a large computation time with these formulations even for instances of small to medium size, as we will report in § 6.
Although the focus of this paper is on computing h(G), let us also mention for the sake of completeness that analogous formulations can be derived for the graph conductance (weighted edge expansion) h vol (G) that was defined in § 1, by optimizing over the set where d = diag(D) is the vector formed by the vertex degrees.For example, the same steps as in (1) yields the MIQCP

Semidefinite relaxations
A well-known lower bound for the edge expansion is the spectral bound.It is based on the second smallest eigenvalue of the Laplacian matrix of the graph, namely h(G) ≥ λ 2 (L) /2.One way to derive this bound is by considering the following SDP relaxation where X models xx ⊤ and we scale X = 1 k X to eliminate the variable k.Proposition 2.1.The optimal solution of the second SDP in (2) is λ 2 (L) /2.

Proof. First observe that 1
n I is a strictly feasible point of the primal SDP and the optimum value is finite, hence strong duality holds.The dual of the second SDP in (2) is max The eigenvalue of W with respect to the eigenvector e is −v 1 + n(v 2 − v 3 ).The other eigenvalues of W are then λ i (L) − v 1 for 2 ≤ i ≤ n.Therefore, we can write the dual as which is a linear program with optimal solution v 1 = λ 2 (L), v 2 = λ 2 (L) /n and v 3 = 0 and optimal value λ 2 (L) /2.
To strengthen the SDP relaxation (2) we round down the upper bound to ⌊ n 2 ⌋ and add the following facet-inducing inequalities of the boolean quadric polytope [35] resulting in the following valid inequalities for for all 1 ≤ i, j, ℓ ≤ n.Note, that in (4c) and (4d) we have to replace 1 k in the rhs by its upper bound 1 in order to obtain a formulation without k.Therefore, we cannot expect these inequalities to strengthen the SDP relaxation significantly.

Illustrative examples for motivation
We motivate our algorithm by considering the example of the graph of the grlex polytope, which is described in [14].Table 1 compares different lower bounds on h(G) for these graphs.The first column indicates the dimension of the polytope and the second column lists the number of vertices in the associated graph.The third column gives the edge expansion that is known to be one for these graphs in all dimensions [14].The spectral bound is displayed in the fourth column.Column 5 lists the optimal value of the SDP relaxation (2) strengthened by inequalities (4) derived from the boolean quadric polytope.Column 6 displays a lower bound that is very easy to compute: the minimum cut of the graph divided by the largest possible size of the smaller set of the bipartition of the vertices, that is ⌊ n 2 ⌋.In the last column, the minimum of the lower bounds ℓ k for 1 ≤ k ≤ ⌊ n 2 ⌋ is listed with ℓ k being a bound related to the solution of (2) for k fixed.The definition of ℓ k follows in § 3.1.
The numbers in the table show that some of these bounds are very weak, in particular, if the number of vertices increases.Interestingly, if we divide the edge expansion problem into ⌊ n 2 ⌋ many subproblems with fixed denominator (as we did to obtain the numbers in column 6) the lower bound we obtain by taking the minimum over all SDP relaxations for the subproblems seems to be stronger than the other lower bounds presented in Table 1.We will, therefore, take this direction of computing the edge expansion, namely, we will compute upper and lower bounds on the problem with fixed k.Using these bounds will

Fixing the size k: Bisection problem
If the size k of the smaller set of the partition of an optimum cut is known, the edge expansion problem would result in a scaled bisection problem.That is, we ask for a partition of the vertices into two parts, one of size k and one of size n − k, such that the number of edges joining these two sets is minimized.This problem is NP-hard [12] and has the following formulations for any k ∈ {1, 2, . . ., ⌊ n 2 ⌋}, but standard branch-and-cut solvers can solve these in reasonable time only for smallsized graphs.Since SDP-based bounds have been shown to be very strong for partitioning problems, cf.[23,30,41,42], we exploit these bounds by developing two kinds of solvers.In the first variant, we develop a tailored branch-and-bound algorithm based on semidefinite programming to solve the bisection problem.In the subsequent sections, we describe how to obtain lower and upper bounds on h k ( § 3.1 and 3.2) as well as further ingredients of this exact solver ( § 3.3).The second variant is presented in § 3.4, where we transform the bisection problem into an instance of a max-cut problem which is then solved using the state-of-the-art solver BiqBin [17].For completeness, a description of BiqBin is given in Appendix A.

SDP lower bounds for the bisection problem
After squaring the linear equality constraint in problem (5) and employing standard lifting and relaxation techniques, we obtain the following SDP relaxation that is generally computationally cheap to solve, Since the bisection for a given simple unweighted graph has to be an integer, we get the following lower bound on the scaled bisection h k , There are several ways to strengthen the above relaxation of the bisection problem.In [42] a vector lifting SDP relaxation, tightened by non-negativity constraints, has been introduced.In our setting, this results in the following doubly non-negative programming (DNN) problem, min where X is a matrix of size (2n+1)×(2n+1).This relaxation can be further strengthened by cutting planes from the Boolean Quadric Polytope.In particular, we want to add the inequalities as Meijer et al. [30] demonstrated that these inequalities are the most promising ones to improve the bound.
The DNN relaxation (8) cannot be solved by standard methods due to the large number of constraints.The additional cutting-planes (9) make the SDP relaxation extremely difficult to solve already for medium-sized instances.Meijer et al. [30] apply facial reduction to the SDP relaxation which leads to a natural way of splitting the set of variables into two blocks in the following way.They present a (2n + 1) × n matrix V such that X = V RV ⊤ for R ∈ S n .Due to this facial reduction, we consider the sets R BP = {R ∈ S n : R ⪰ 0} and X BP being the set of matrices in S 2n+1 satisfying all polyhedral constraints.The facially reduced DNN relaxation of the bisection problem then reads min with appropriate matrix L BP .Using an alternating direction method of multipliers (ADMM) provides (approximate) solutions to this relaxation even for graphs with up to 1000 vertices.The steps to be performed in this ADMM algorithm result in projections onto the respective feasible sets.For projections onto polyhedra, Dykstra's projection algorithm is used.A careful selection of non-overlapping cuts, warm starts, and an intelligent separation routine are further ingredients of this algorithm in order to obtain an efficient solver for the SDP (8) enhanced with inequalities (9).A post-processing algorithm is also introduced to guarantee a valid lower bound.Using this algorithm, we can compute strong lower bounds for each k with reasonable computational effort.

A heuristic for the bisection problem
The graph bisection problem can also be written as a quadratic assignment problem (QAP) [25].To do so, we set the weight matrix W to be the Laplacian matrix L of the graph and the distance matrix D to be a matrix with a top left block of size k with all ones and the rest zero.The resulting QAP for this weight and distance matrix is In this formulation, the vertices mapped to values between 1 and k by the permutation π are chosen to be in the set of size k in the partition.To compute an upper bound u k on h k , we can use any heuristic for the QAP and divide the solution by k.Simulated annealing is a well-known heuristic to compute an upper bound for the QAP, we implement the algorithm introduced in [8].We use a slightly different parameter setting that we determined via numerical experiments.That is, the initial temperature is set to k 2 /( n 2 ) • tr(L), the number of transformation trials at constant temperature is initially set to n and increased by a factor of 1.15 after each cycle, and the cooling factor is set to 0.7.After every trial, we additionally perform a local search strategy to find the local minimum.
Other well-performing heuristics for the QAP are tabu search, genetic algorithms, or algorithms based on the solution of the SDP relaxation like the hyperplane rounding algorithm.Some of these heuristics have the potential to be superior to simulated annealing.However, as we will see later in the study of our experiments, the bounds we obtain with the simulated annealing heuristic are almost always optimal for the size of instances we are interested in.

A branch-and-bound algorithm for the bisection problem
The branch-and-bound algorithm is a helpful technique to solve optimization problems to optimality by dividing the feasible region.Existing exact algorithms for the graph bisection problem are presented in [19,24] considering instances with up to 148 vertices, most of them being very sparse, and in [2] for large sparse graphs.Moreover, most of the computational results are presented for equipartition problems and not for general sizes of the partition.Implementations of these algorithms are not available and from the results in the papers these algorithms might not be successful on our problems.
We, therefore, implement an open-source solver to solve graph bisection problems of medium size to optimality using the ingredients described in the previous sections, namely, we implement an ADMM to obtain approximate solutions of the SDP bound as described in § 3.1 and we implement a simulated annealing procedure to obtain the upper bound as described in § 3.2.
Branching In the case of binary optimization, a natural way of branching is to fix a variable to 0 or 1.We base our branching decision on the solution of the relaxation of the subproblem.Namely, we branch on the node with the corresponding value in x 1 being closest to 0.5, i.e., most fractional.It turns out that we can write the subproblems again as problems of a similar form as the bisection problem but with fewer variables.In particular, if we set a variable x i to be 0, we can write the problem as min{x ⊤ Lx : e ⊤ x = k}, where x is obtained from x by deleting x i and we get L by deleting the i-th row and column of L. The subproblem where we set x i = 1 can be written as min{x ⊤ Lx + c : e ⊤ x = k − 1}, with x again resulting from x by deleting entry x i and L is obtained from L by adding the i-th row and column to the diagonal before deleting them and c = L ii .Note that for both types of subproblems, although they are no bisection problems anymore, we can still use the methods discussed in the Sections 3.1 and 3.2 to compute lower and upper bounds.

Strategy on adding strengthening inequalities
In our branch-and-bound algorithm, similar to the strategy in [17], in the root node we first compute the (DNN) bound (8) and then iteratively add violated triangle inequalities (9).The improvement attained by adding triangle inequalities is stored in diff and handed down to the child nodes.We use it in the child nodes for the decision of whether or not we add violated triangle inequalities to (DNN) since adding triangle inequalities is time-consuming.In case the solution of (DNN) plus diff does not suffice to prune the node, we refrain from adding triangle inequalities.In case we do continue adding violated triangle inequalities and are not able to prune the node, we update diff accordingly for the subsequent child nodes.
Lower bound verification/early stopping In preparation for § 4, we add the option to stop the branch-and-bound algorithm as soon as there is no open subproblem with a lower bound below a given threshold, that is, if the lower bound on the bisection problem is greater or equal than that threshold.This functionality can be used to verify a given lower bound without having to solve the bisection problem.
Before moving on to § 4, we present a second variant on how to solve the bisection problem.

Transformation to a max-cut problem
A different approach to solving the graph bisection problem is to transform it to a maxcut problem and use a max-cut solver, e.g. the open source parallel solver BiqBin [17], see also Appendix A. To do so, we first need to transform the bisection problem into a quadratic unconstrained binary problem (QUBO).
Proof.First note that Let x ∈ {0, 1} n .Then we have Note that e ⊤ x − k is integer for x ∈ {0, 1} n .Hence, for any infeasible x ∈ {0, 1} n , the objective is greater than the given upper bound x⊤ Lx, and therefore the minimum can only be attained for x ∈ {0, 1} n with e ⊤ x = k.
Barahona et al. [6] showed that any QUBO problem can be reduced to a max-cut problem by adding one additional binary variable.In our context, this means the following.
Corollary 3.2.Let G = (V, E) and let G ′ be the complete graph with vertex set V ∪{v 0 }.Let the weights c uw on the edges of G ′ be as follows.
Then the minimum bisection of G where one side of the cut has size k is equal to Since max-cut solvers can benefit from edge weights of the input graph being integer, a possible choice for µ k is an upper bound on the bisection problem plus 1 /4.Note that we choose µ k to be as small as possible by doing so.
The formulation of the max-cut problem in ±1 variables additionally requires x v = 1 to hold.Because of the symmetry of the cut, we can omit this constraint.Due to our choice of the penalty parameter, it holds that on one side of the maximum cut, there are exactly k + 1 vertices, including vertex v.These k vertices on the same side as v are the vertices in the set of size k in the optimum of the bisection problem.
To summarize, we can solve the bisection problem by solving a dense max-cut problem.With this, we now have all the tools needed for our new split & bound approach.

Split & bound
We now assemble the tools developed in § 3 to compute the edge expansion of a graph by splitting the problem into ⌊ n 2 ⌋ many bisection problems.Since the bisection problem is NP-hard as well, we want to reduce the number of bisection problems we have to solve exactly as much as possible.To do so, we start with a pre-elimination of the bisection problems.This procedure aims to exclude subproblems unnecessary for the computation of the edge expansion of the graph.Computing the edge expansion by considering the remaining values of k is summarized in Algorithm 2 below.We now explain the pre-elimination step and further ingredients of our algorithm.

Pre-elimination
The size k of the smaller set of the partition can theoretically be any value from 1 to ⌊ n 2 ⌋.However, it can be expected that for some candidates, one can quickly check that the optimal solution cannot be attained for that k.As a first quick check, we use the cheap lower bound ℓ k obtained by solving the SDP (6) in combination with the upper bound introduced in the § 3.2.We do not need to further consider values of k where the lower bound ℓ k of the scaled bisection problem is already above an upper bound on the edge expansion.A pseudo-code of this pre-elimination step is given in Algorithm 1.
Compute an upper bound u k using the heuristic from § 3.2; 3 Compute the lower bound ℓ k from ( 7) by solving the cheap SDP (6); The hope is that many values of k can be excluded from computing the edge expansion.Clearly, this heavily depends on the instance itself, as in the worst case, it might happen that for many different values of k, the value of h k is close to the optimum.
We can further reduce the number of candidates for k by computing a tighter lower bound lk by solving the DNN relaxation (8) with additional cutting planes.In our implementation we do not compute lk as part of the pre-elimination, since this bound is computed in the root node of the branch-and-bound tree in the algorithm from § 3.3.
Impact of pre-elimination on sample instances Figures 1 and 2 display the bounds associated with four different graphs.For the graph of the grevlex polytope in dimension 7, considering the bounds u k and lk the only candidates for k where the optimal solution can be found are 12 and 14.For the grevlex polytope in dimension 8, the sizes 17 and 18 remain as the only candidates.Also for a graph associated to a randomly generated 0/1-polytope and to a network graph, about 2/3 of the potential values of k can be excluded already by considering the cheap lower bound ℓ k .

Stopping exact computations early and updating u *
All values of k that are not excluded in the pre-elimination step have to be further examined.For those values we run our branch-and-bound algorithm (either for k-bisection or max-cut, depending on the variant chosen by the user) to compute the scaled bisection h k .We can stop the branch-and-bound algorithm as soon as the lower bound of all open nodes in the branch-and-bound tree is larger or equal to u * .(Remember that u * is an upper bound on the edge expansion of the graph but not necessarily an upper bound on h k .)A simple way to implement this stopping criterion is to initialize the branch-and-bound algorithm for the k-bisection problem with ⌈u * k⌉ as an "artificial" upper bound.
In case h k < u * , we can update u * which might lead to eliminating further values from I.This fact is also part of the motivation for the order of choosing k for computing h k , as described in the next section.

Order of selecting values k from I
We consider the order of computing h k in ascending order based on their upper bounds u k .The motivation for this choice is as follows.
Remember that u * = min k u k is a global upper bound on the edge expansion.The most promising values for k to even further improve this bound are those with small u k .Therefore, before starting the branch-and-bound algorithm for the values k left after preelimination, we do another 30 trials of simulated annealing for each of these to hopefully further improve the upper bound.
Moreover, we run the exact computation of h k in this order since also during the branch-and-bound algorithm, the upper bound might drop further and this will improve the global upper bound u * most likely for those candidates with small u k .
An improvement of the upper bound u * means that there is a possibility to further eliminate values k from I.But even for values k that can not be eliminated, we obtain smaller artificial upper bounds and hence the computation of these bisection problems may be stopped earlier.
We summarize all the steps in Algorithm 2.

Algorithmic verification of lower bound
We close this section by addressing the important consideration that we are not interested in the exact value of the edge expansion in some applications, but want to check whether certain values are valid lower bounds on h(G).A lower bound c ≤ h(G), for some constant c > 0, means that the graph is a c-expander.The value of this lower bound means that the graph expands by at least that much.This also arises in the context of the  12 compute h k using the max-cut solver, initialize the lower bound for max-cut as offset − ⌈u * k⌉; Mihail-Vazirani conjecture on 0/1-polytopes where one wants to check whether h(G) ≥ 1 where G is the graph of a 0/1-polytope.
Our split & bound algorithm can also be used to verify a lower bound.Proof.Assume we initialize u * = υ.If we find a better upper bound (or some computed value h k is smaller than υ), this is a certificate that the given value υ is not a valid lower bound since we found a better solution.Otherwise, if the upper bound never gets updated, this means the provided bound is indeed a valid lower bound on the edge expansion of the graph.

Parametric optimization
Another approach to compute h(G) is following a discrete Newton-Dinkelbach algorithm.Dinkelbach [10] gave a general classical framework to solve (non)-linear fractional programs.The program one aims to solve is min x∈F f (x), where the objective f is a fraction of (non)-linear functions.In our case this is The main component of the algorithm is to form the following parametrized objective function and the corresponding parametrized optimization problem This problem then has the following useful properties.
Proposition 5.1.P (0) = ζ min (G) and P is a strictly decreasing concave piecewise linear function over R + whose unique root is equal to h(G).Consequently, Proof.We have where the penultimate equality is from symmetry of the cut function ζ(S) = |∂S| = ∂S ′ .Finiteness of F and linearity of g γ in γ tells us that P is the pointwise minimum of finitely many affine (in γ) functions, and so P is a concave piecewise linear function.The strictly decreasing property was shown in [10, Lemma 3] for general nonlinear fractional problems.
This implies that the edge expansion of a graph can be computed using a root-finding algorithm for the function P .One evaluation of P for a given γ still means solving a binary quadratic problem with two linear inequalities.Hence, reducing this number of evaluations is crucial to compute h(G) in reasonable time.There are several strategies to do so, such as binary search.Our approach is to evaluate P starting with γ 1 equal to some good upper bound on h(G) (in our experiments, we used our heuristic from § 3.2).We are already done if we have found the optimum with our heuristic, that is when P (γ 1 ) = 0. Otherwise, there is some x 1 ∈ F such that g γ (x 1 ) < 0 and therefore f (x 1 ) < γ.This means that f (x 1 ) is a better upper bound than γ 1 .Hence, we now set γ 2 = f (x 1 ) and repeat until we find the optimum as described in Algorithm 3. Since P (γ) < 0 if and only if γ > h(G), the stopping criterion is checking whether P (γ) < 0 at the current iterate.
The superlinear convergence rate of Dinkelbach's algorithm was established in [40].We derive a similar convergence result for Algorithm 3.
We solve Q(γ) by transforming it into a max-cut problem.The first step towards achieving this is to obtain an exact formulation as a QUBO using binary encoding of the slack variables and the penalty parameter suggested in [18,Thm. 15].To aid our derivation, let us denote the integer n s and vector v ns ∈ R ns+1 by Proof.For any x ∈ F it holds that e ⊤ x = 1 + s = ⌊ n 2 ⌋ − t for some slack variables s and t with 0 ≤ s, t ≤ ⌊ n 2 ⌋ − 1.In fact, any upper bound on s and t greater or equal than ⌊ n 2 ⌋ − 1 is fine, since from e ⊤ x = 1 + s and s ≥ 0 it follows that e ⊤ x ≥ 1 and from e ⊤ x = ⌊ n 2 ⌋ − t and t ≥ 0 it follows that e ⊤ x ≤ ⌊ n 2 ⌋.The smallest possible value for n s is ⌈log 2 (⌊ n 2 ⌋)⌉ − 1, since this gives an upper bound of 2 ns+1 − 1 on s and t.
The problem Q(γ) can then be equivalently formulated as the following QUBO, with σ > γ n n.
Proof.Let g(x, α, β) denote the objective function of (11).For any feasible vector x ∈ F there exist uniquely defined α x , β x such that g(x, α x , β x ) = γ d x ⊤ Lx − γ n e ⊤ x.Thus, the objective function of Q(γ) and ( 11) coincide for x ∈ F.Moreover, for x ′ it holds that g(x ′ , α x ′ , β x ′ ) = 0.For x ∈ {0, 1} n \ F there do not exist α, β ∈ {0, 1} ns+1 such that both equalities e ⊤ x − α ⊤ v ns = 1 and e ⊤ x + β ⊤ v ns = n 2 are satisfied, as one of the slack variables has to be negative in order to fulfill the constraints.Additionally, since L is positive semidefinite we can conclude that Therefore, Q(γ) is an equivalent formulation of (11) The unconstrained binary quadratic program (11) can again be transformed to a max-cut problem, as explained in [6] for example.Applied to our problem we obtain the following result.

Benchmark instances
Randomly generated 0/1-polytopes The first class of graphs are the graphs of random 0/1-polytopes.The polytopes are generated by randomly selecting n d vertices of the polytope in dimension d, i.e., n d different 0/1-vectors in dimension d.To obtain the graph, we then solve a linear programming feasibility problem to check whether there is an edge for a given pair of vertices.For any pair (d, n d ) with d ∈ {8, 9, 10} and n 8 ∈ {164, 189}, n 9 ∈ {153, 178, 203, 228, 253, 278}, and n 10 ∈ {256, 281}, we generated 3 random 0/1-polytopes.
Grlex and grevlex graphs Another class of graphs we consider are the graphs of grlex and grevlex polytopes introduced and characterised by Gupte and Poznanović [14].The grevlex-d and grlex-d instances of our benchmark set are the corresponding graphs of the polytopes in dimension d.

DIMACS and Network graphs
The last category of graphs originates from the graph partitioning and clustering application.The set of DIMACS instances are the graphs of the 10th DIMACS challenge on graph partitioning and graph clustering [5] with at most 500 vertices.Additionally, we consider some more network graphs obtained from the online network repository [36].

Discussion of the experiments
We compare different algorithms for computing the edge expansion of a graph, namely 1. Split & bound Algorithm 2 using the k-bisection solver, 2. Split & bound Algorithm 2 using the max-cut solver, 3. Fractional programming using Discrete Newton-Dinkelbach's method in Algorithm 3, 4. Gurobi for solving the MIQCP.k-bisection solver vs. max-cut solver In our preliminary numerical experiments, the max-cut variant demonstrated superior performance compared to the branch-and-bound algorithm for k-bisection.For example for the instance chesapeake from the DIMACS set, it took 2 seconds to compute the edge expansion using the max-cut variant compared to 9.8 seconds with the k-bisection branch-and-bound algorithm.Therefore we only report the results for the max-cut variant to solve the scaled bisection problems.Additionally, the fact that the max-cut solver BiqBin is parallelized further motivates this selection.
Algorithm 2 vs. Algorithm 3 vs.Gurobi The detailed results of our experiments are given in Tables 2 to 6.In each of the tables, the first column gives the name of the instance followed by the number of vertices and edges.Column 4 reports the optimal solution, i.e., the edge expansion of the graph.
In the split & bound section of the table, the first two columns give the global lower and upper bound after the pre-elimination Algorithm 1.The number of candidates for k after the pre-elimination is given in column 3.In column 4 we report the number of indices k ∈ I we were able to eliminate after solving the root node of the branch-andbound tree.Column 5 lists the total number of branch-and-bound nodes in the max-cut algorithm for all values of k considered.The last two columns display the time spent in the pre-elimination the total time (including pre-elimination) of the algorithm.
In the section for Dinkelbach's algorithm, the first column gives the first guess for the edge expansion, i.e., the first trial for γ.As described before, we take the upper bound from the heuristic for this initialization.Note, that this first guess may differ from u * , since in the pre-elimination step of split & bound we perform 30 additional rounds of simulated annealing for all indices k ∈ I. Column 2 indicates how many parametrized problems P γ i have been solved, and column 3 gives the total number of branch-andbound nodes for solving all parametrized problems.The fourth column of the results of Dinkelbach's algorithm displays the total time, including running the heuristic to obtain the first guess.
The final column of the tables holds information about computing the edge expansion using Gurobi.For the graphs from the randomly generated polytopes, Gurobi did not succeed to solve any of these instances within a time limit of 3 hours, we therefore report the gap after this time limit in Table 2.In Tables 3 to 6 we report the time for computing the edge expansion.
We highlight in the tables, which of our two algorithms performs better.

Algorithm 2 (Split & bound)
As can be seen in all tables, the pre-elimination phase of split & bound only leaves a comparably small number of candidates for k to be further investigated.Remember that the number of potential candidates is ⌊ n 2 ⌋, whereas |I| is the number of candidates that remain after the pre-elimination.For the randomly generated 0/1-polytopes on average only 12% of the candidates have to be further examined, for the other instance classes we are left with approximately 20% on average.This indicates that already the cheap SDP bound is of good quality.
We also observe that the SDP bound in the root-node of the branch-and-bound tree is of high quality: in 48 out of the 67 instances the gap is closed within the root node for all candidates to be considered.
The heuristic for computing upper bounds also performs extremely well: for almost all instances the upper bound found is the edge expansion of the graph, see columns titled h(G) and u * .In fact, only for 3 instances the heuristic fails to find the optimal solution.
Overall, we solve almost all of the considered instances within a few minutes, for very few instances the branch-and-bound tree grows rather large and therefore computation times exceed several hours.Algorithm 3 (Discrete Newton-Dinkelbach) Whenever the heuristic already returns the value of the optimal solution, we only have to solve one parametrized problem to certify the optimality of this value.For most of the instances tested, this certificate is already obtained in the root node of the branch-and-bound tree.However, there are many instances where BiqBin terminates because of numerical problems even for the first parametrized problem, see Table 2.This in particular arises when γ d and γ n (see Section § 5.1) are large.

Solving the MIQCP with Gurobi
To compute the edge expansion using Gurobi, we input the last formulation of (1) adding the redundant constraint y ≥ 0. Without this constraint, Gurobi terminated only after 1.65 hours/3 548 work units (resp.more than 24 hours/59 000 work units) on a graph with 29 vertices and 119 edges (resp.37 vertices and 176 edges) corresponding to the grevlex polytope in dimension 7 (resp.8).
Adding the redundant constraint, Gurobi is very efficient for graphs with an expansion less than one, see Tables 5 and 6, but as soon as the expansion (and also the number of vertices of the graph) gets larger, Gurobi cannot solve the instance within a few hours, see Tables 2 and 4.

Conclusion
To summarize the results, we give a performance profile in Figure 3. Gurobi solves the MIQCP very efficiently for several instances, but fails to yield results for others within a time limit of 3 hours.It is the clear winner for instances with very small edge expansion.Comparing split & bound with the algorithm following the Discrete Newton-Dinkelbach method, we observe the following behavior.For the grlex instances, Dinkelbach performs extremely well compared to the split & bound approach, see Table 3. Whereas for the grevlex instances in Table 4, we observe that, except for dimension 13, the split & bound algorithm by far outperforms Algorithm 3. In addition to the already mentioned, there are some other instances where the difference in the runtimes between the two algorithms is significant.For example, on the instances rand01-9-153-0 and malaria genes HVR1 split & bound clearly dominates Algorithm 3, whereas the latter is significantly better on the instances rand01-9-2781 and sp-office.
The conclusion is that in general for graphs with larger edge expansion, the split & bound algorithm is best, and for graphs with small edge expansion Algorithm 3 has a better performance than split & bound, but there are a few exceptions, and the difference in the total time of solving an instance can be quite large.
Overall, we conclude that with our algorithms we can compute the edge expansion of various graphs of size up to around 400 vertices and no other algorithm can achieve this.The time for solving an instance varies, it can be a few seconds for very structured instances and it can exceed one hour, in particular for the instances coming from 0/1polytopes with rather large expansion.For standard branch-and-cut solvers like Gurobi these instances are out of reach.

Summary and future research
We developed a split & bound algorithm as well as an algorithm applying Dinkelbach's idea for fractional programming to compute the edge expansion of a graph.The split-   (Dinkelbach) and Gurobi for graphs of random 0/1polytopes.The last column displays the gap reported by Gurobi after a time limit of 3 hours.ting refers to the fact, that we consider the different values of k (k being the size of the smaller partition) separately.We used semidefinite programming in both phases of our algorithm: on the one hand, SDP-based bounds are used to eliminate several values for k and we use SDP-based bounds in a branch-and-bound algorithm that solves the problem for k fixed.I.e., we also developed a branch-and-bound solver for the k-bisection problem using bounds from semidefinite programming.Also, the algorithm following the Dinkelbach framework uses semidefinite programming in order to solve the underlying parametrized problems.Through numerical results on various graph classes, we demonstrate that our split-and-bound algorithm is a robust method for computing the edge expansion while Dinkelbach's algorithm and Gurobi are very sensitive with respect to the edge expansion of the graph.
In some applications, one is not interested in the exact value of the edge expansion but wants to check whether a certain value is a lower bound on the edge expansion, e.g., to check the Mihail-Vazirani conjecture on 0/1-polytopes.We implemented an option in our algorithm that enables this feature of verifying a given lower bound.
As a heuristic, we use a simulated annealing approach that works very well for the problem sizes we are interested in.However, if one wants to obtain high-quality solutions for larger instances, a more sophisticated heuristic will be needed.Tabu-search, genetic algorithms, or a heuristic in the spirit of the Goemans-Williamson rounding could be potential candidates.In future research, we will also investigate convexification techniques by using recent results on fractional programming [20] and on exploiting submodularity of the cut function as has been done for mixed-binary conic optimization [4].
A. Computing the maximum cut of a graph Some of our algorithms for computing h(G) rely on finding the maximum cut in a graph.For computing the value of the max-cut, we will use the SDP-based solver BiqBin [17].Note that the software BiqBin can not only compute the max-cut in a graph and solve quadratic unconstrained binary problems (QUBOs) but it is also applicable to linearly constrained binary problems with a quadratic objective function.However, we only need the max-cut solver in this work, and briefly describe the main ingredients in this section.
BiqBin is a branch-and-bound algorithm that uses a tight SDP relaxation as upper bound and the celebrated Goemans-Williamson rounding procedure to generate a highquality lower bound on the value of the maximum cut in a graph.To be more precise, the SDP max ⟨L, X⟩ : diag(X) = e, A(X) = b, X ⪰ 0 (12) where A(X) ≤ b models a set of triangle-, pentagonal-and heptagonal-inequalities is approximately solved using a bundle method.To do so, only the inequality constraints are dualized yielding the nonsmooth convex partial dual function Evaluating the dual function f (γ) and computing the subgradient amounts to solving an SDP that can be efficiently computed using an interior-point method tailored for this problem.It provides us with the matching pair (X γ , γ) such that f (γ) = b ⊤ γ + ⟨ 1 /4L − A ⊤ (γ), X γ ⟩.Moreover, the subgradient of f at γ is given by ∂f (γ) = b − A(X γ ).For obtaining an approximate minimizer of problem min γ {f (γ) : γ ≥ 0}, the bundle method is used.We refer to [38] for more details.Note that interior-point methods are far from computing a solution of (12) already for small graphs due to the number of constraints being too large and therefore already forming the system matrix is an expensive task or even impossible due to memory requirements.
BiqBin dominates all max-cut solvers based on linear programming and is comparable to the SDP-based solver BiqCrunch [26].Moreover, BiqBin is available as a parallelized version.

Figure 1 :
Figure 1: Lower and upper bounds for each k.

Figure 2 :Algorithm 2 : 3 Run heuristic from § 3 .2 and update u k ; 4 if u k < u * then 5 u * ← u k ; 6 update I; 7 8 if variant = k-bisection then 9 compute
Figure 2: Lower and upper bounds for each k.

Proposition 4 . 1 .
Let υ be a given scalar and suppose we initialise Algorithm 2 with u * = υ.Then υ is a valid lower bound on h(G) if and only if the algorithm terminates without updating u * .

Figure 3 :
Figure 3: Performance comparison of the exact algorithms.Note the different scale on the x-axis: the plot on the left displays the time range from 0 to 500 seconds, the plot on the right from 500 seconds to 3 hours.

Table 1 :
Comparison of lower bounds for graphs from the grlex polytope in dimension d.