Computing the Edge Expansion of a Graph Using Semidefinite Programming

Abstract

Computing the edge expansion of a graph is a famously hard combinatorial problem for which there have been many approximation studies. We present two versions of an exact algorithm using semidefinite programming (SDP) to compute this constant for any graph. The SDP relaxation is used to first reduce the search space considerably. One version applies then an SDP-based branch-and-bound algorithm, along with heuristic search. The other version transforms the problem into an instance of a max-cut problem and solves this using a state-of-the-art solver. Numerical results demonstrate that we clearly outperform mixed-integer quadratic solvers as well as another SDP-based algorithm from the literature.

1 Introduction

Let G be a simple graph on n vertices. The (unweighted) edge expansion, also called the Cheeger constant or conductance or sparstest cut, of G is defined as

$$\begin{aligned} h(G)= \min _{\emptyset \ne S\subset V}\,\frac{|\partial S|}{\min \{|S|, |S^{\prime }| \}} = \min _{S\subset V}\,\left\{ \frac{|\partial S|}{|S|} :1 \le |S| \le \frac{n}{2}\right\} , \end{aligned}$$

(1)

where $\partial S = \{(i,j)\in E(G) : i \in S, j\in S^\prime \}$ is the cut-set associated with S, and $S^{\prime }=V{\setminus } S$. This constant is positive if and only if the graph is connected. A graph with $h(G)$ small is said to have a bottleneck. A threshold for good expansion properties is having $h(G)\ge 1$, which is desirable in many applications.

Edge expansions arise in the study of expander graphs, which have applications in network science, coding theory, cryptography, complexity theory, etc. [9, 13]. The famous Mihail-Vazirani conjecture [7, 20] in polyhedral combinatorics claims that the graph (1-skeleton) of any 0/1 polytope has edge expansion at least 1.

Computing the edge expansion is related to the uniform sparsest cut problem which is defined as

$$\begin{aligned} \phi (G)= \min _{\emptyset \ne S\subset V}\,\frac{|\partial S|}{|S|\cdot |S^{\prime }|} = \min _{S\subset V}\,\left\{ \frac{|\partial S|}{|S|\cdot |S^{\prime }|} :1 \le |S| \le \frac{n}{2}\right\} . \end{aligned}$$

(2)

Since $\frac{n}{2} \le |S^{\prime }| \le n$ it holds that $|S|\cdot |S^{\prime }| \le |S|\cdot n \le 2 |S|\cdot |S^{\prime }|$ and hence

$$\begin{aligned} h(G)\le n\cdot \phi (G)\le 2h(G)\end{aligned}$$

and a cut $(S,S^{\prime })$ that is $\alpha$-approx for $\phi (G)$ is a $2\alpha$-approx for $h(G)$.

Both $h(G)$ and $\phi (G)$ are NP-hard to compute [15], and the latter has received considerable attention from the approximation algorithms community.

The classical bounds on $h(G)$ are from the spectral relation due to Alon and Milman [1] who showed that $\frac{\lambda _2}{2} \le h(G)\le \sqrt{2\lambda _2\varDelta }$, where $\lambda _{2}$ is the second smallest eigenvalue of the Laplace matrix of G (spectral gap) and $\varDelta$ is the maximum degree of G. The best-known approximation for $\phi (G)$ is the famous $\mathcal {O}(\sqrt{\log n})$ factor by Arora et al. [2] which improved upon the earlier $\mathcal {O}({\log n})$-approximation [15]. Meira and Miyazawa [19] developed a branch-and-cut algorithm for $h(G)$ using LP relaxations and SDP-based heuristics. To the best of our knowledge, there is no other exact solver for the edge expansion.

Contribution and Outline. We develop an algorithm in two phases for computing the edge expansion of a graph. In the first phase, our algorithms splits 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 (B &B) algorithm, in the second version we transform the subproblem into an instance of a max-cut problem and solve this max-cut problem to optimality. We perform numerical experiments on different types of instances which demonstrate the effectiveness of our results. To the best of our knowledge, no other algorithms are capable of computing the edge expansion for graphs with a few hundred vertices.

In Sect. 2 we formulate the problem as a mixed-binary quadratic program and present an SDP relaxation. Section 3 investigates a related problem, namely the k-bisection problem. Our algorithm is introduced in Sect. 4, the performance of the algorithm is demonstrated in Sect. 5, followed by conclusions in Sect. 6.

Notation. The trace inner product for two real symmetric matrices X, Y is defined as $\langle X, Y \rangle = \text {tr}(XY)$ and the operator $\textrm{diag}(X)$ returns the main diagonal of matrix X as a vector. We denote by e the vector of all ones, and define $E=ee^\top$.

2 QP Formulation and a Semidefinite Relaxation

Consider a graph G with vertices $V = \{1,\dots ,n\}$ and its Laplacian matrix L, which is defined as $L = {\text {Diag}}(d) - A$, where A is the adjacency matrix of the graph and $d = (d_1,d_2,\dots ,d_n)$ is the vector of vertex degrees. Any binary vector $x \in \{0,1\}^n$ represents a cut in this graph and the value of this cut can be computed as $x^\top L x$. Hence, the expansion can be computed as

$$\begin{aligned} h(G) &= \min _{x \in \{0,1\}^n} \bigg \{\frac{x^\top L x}{e^\top x} :1 \le e^\top x \le \frac{n}{2}\bigg \}\\ &= \min _{\begin{array}{c} x \in \{0,1\}^n\\ y\in \mathbb {R} \end{array}} \bigg \{y \;:\dfrac{x^\top Lx}{e^\top x } \le y,\ 1 \le e^\top x \le \frac{n}{2} \bigg \}, \end{aligned}$$

which can be equivalently written as the mixed-binary quadratic problem

$$ h(G) = \min _{\begin{array}{c} x \in \{0,1\}^n\\ y \in \mathbb {R} \end{array}} \big \{y \;:x^\top L x - ye^\top x \le 0, 1 \le e^\top x \le \frac{n}{2}\big \}. $$

Standard solvers like Gurobi, CPLEX, Mosek, or CBC can handle this formulation but require a large computation time even for small instances. For example, Gurobi (version 11.0 with default parameter setting) terminated after 1.65 h/3 548 work units (resp. more than 24 h/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) [10].

A way to derive the spectral lower bound on $h(G)$ is via the SDP relaxation

$$\begin{aligned} \begin{array}{llllll} h(G) \ge &{}\min _{\widetilde{X}, k} ~ &{} \frac{1}{k} \langle L, \widetilde{X} \rangle &{}~ =~ \min _{X} ~ &{} \langle L, X \rangle &{} =~ \frac{\lambda _{2}(L)}{2},\\[0.2em] &{}\text {s.t.}~ &{} \text {tr}(\widetilde{X}) = k &{}\qquad \text {s.t.}~ &{} \text {tr}(X) = 1\\[0.1em] &{}&{} \langle E, \widetilde{X} \rangle = k^2 &{}&{} 1 \le \langle E, X \rangle \le \frac{n}{2} \\ &{}&{} 1 \le k \le \frac{n}{2} &{}&{} X \succcurlyeq 0\\ &{}&{} \widetilde{X} \succcurlyeq 0 \end{array} \end{aligned}$$

(3)

where $\tilde{X}$ models $xx^\top$ and we scale $X = \frac{1}{k}\widetilde{X}$ to eliminate the variable k. By considering the dual of the second SDP, it can be shown that the optimum is $\nicefrac {\lambda _{2}(L)}{2}$. To strengthen (3) we round down the upper bound to $\lfloor \frac{n}{2}\rfloor$ and add the following inequalities.

Lemma 1

The following are valid inequalities for X for all $1 \le i,\ j,\ \ell \le n$.

$$\begin{aligned} 0 \le {X}_{ij} & \le {X}_{ii} \end{aligned}$$

(4a)

$$\begin{aligned} {X}_{i\ell } + {X}_{j\ell } - X_{ij}& \le {X}_{\ell \ell } \end{aligned}$$

(4b)

$$\begin{aligned} {X}_{ii} + {X}_{jj} - {X}_{ij} \le \nicefrac 1 k &\le 1 \end{aligned}$$

(4c)

$$\begin{aligned} {X}_{ii} + {X}_{jj} + {X}_{\ell \ell } - {X}_{ij} - {X}_{i\ell } - {X}_{j\ell } \le \nicefrac {1}{k} & \le 1 \end{aligned}$$

(4d)

Proof

The inequalities result from scaling $\widetilde{X}$ in the facet-inducing inequalities of the boolean quadric polytope for $\widetilde{X}$ which are given as $0\le \widetilde{X}_{ij} \le \widetilde{X}_{ii}$, $\widetilde{X}_{i\ell } + \widetilde{X}_{j\ell } -\widetilde{X}_{ij} \le \widetilde{X}_{\ell \ell }$, $\widetilde{X}_{ii} + \widetilde{X}_{jj} - \widetilde{X}_{ij} \le 1$, $\widetilde{X}_{ii} + \widetilde{X}_{jj} + \widetilde{X}_{\ell \ell } - \widetilde{X}_{ij} - \widetilde{X}_{i\ell } - \widetilde{X}_{j\ell } \le 1$.    $\square$

Note, that in (4c) and (4d) we have to replace $\frac{1}{k}$ by its upper bound 1 in order to obtain a formulation without k.

Table 1 compares different lower bounds on the example on graphs of the grlex polytope, which is described in [10]. The first three columns indicate the dimension of the polytope, the number of vertices in the associated graph, and the edge expansion that is known to be one for these graphs [10]. In the fourth and fifth columns the spectral bound and the strenghened SDP bound (3) are displayed. Column 6 displays a 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 partition. In the last column, the minimum of the lower bounds $\ell _k$ for $1 \le k \le \lfloor \frac{n}{2} \rfloor$ is listed with $\ell _k$ being a bound related to the solution of (3) for k fixed. The definition of $\ell _k$ follows in Sect. 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 $\lfloor \frac{n}{2} \rfloor$ 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.

Table 1. Lower bounds for graphs from the grlex polytope in dimension d.

3 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 [8] and has the following formulation,

$$\begin{aligned} h_k = \frac{1}{k} ~ \min _{x\in \{0,1\}^n} \big \{x^\top L x :e^\top x = k\big \} \end{aligned}$$

(5)

which standard branch-cut solvers can solve in reasonable time only for small-sized graphs.

Since SDP-based bounds have been shown to be very strong for partitioning problems [cf 14, 18, 21, 22], we exploit these bounds by developing two kinds of solvers. We develop a tailored B &B 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$ (Sect. 3.1 and Sect. 3.2) as well as further ingredients of this exact solver (Sect. 3.3). An alternative to this B &B solver is presented in Sect. 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 [11].

3.1 SDP Lower Bounds for the Bisection Problem

A computationally cheap SDP relaxation of the bisection problem is

$$\begin{aligned} \ell _{\text {bisect}}(k) = \min _{X, x} \big \{ \langle L,X \rangle :\text {tr}(X) = k, \langle E,X \rangle = k^2, \text {diag}(X) = x, X \succeq xx^\top \big \}. \end{aligned}$$

(6)

Since the k-bisection of a graph has to be an integer, we get that

$$\ell _k = \frac{\lceil \ell _{\text {bisect}}(k) \rceil }{k}$$

is a lower bound on the scaled bisection $h_k$.

There are several ways to strengthen (6). In [22] 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.

$$\begin{aligned} {\begin{matrix} \min _X ~~&{}\frac{1}{2} \langle L, X^{11}+X^{22} \rangle \\ \text {s.t.}~~&{}\text {tr}(X^{11}) = k, ~ \langle { E}, X^{11} \rangle = k^{2}, \\ &{}\text {tr}(X^{22}) = n-k, ~ \langle { E}, X^{22} \rangle = (n-k)^{2}, \\ &{}\text {diag}(X^{12})={0},~ \text {diag}(X^{21})={0}, ~ \langle { E}, X^{12}+X^{21}\rangle = 2k(n-k), \\ &{} X= \begin{pmatrix} 1 &{} (x^1)^\top &{} (x^2)^\top \\ x^1 &{} X^{11} &{} X^{12} \\ x^2 &{} X^{21} &{} X^{22} \\ \end{pmatrix} \succeq {0}, ~~x^i=\text {diag}(X^{ii}), ~~ i=1,2,\\ &{}X \ge 0, \end{matrix}} \end{aligned}$$

(7)

where X is a matrix of size $(2n+1) \times (2n+1)$.

Meijer et al. [18] use this relaxation and strengthen it further by cutting planes from the boolean quadric polytope. Since this SDP cannot be solved by standard methods due to the large number of constraints, they present an alternating direction method of multipliers (ADMM) to (approximately) solve this relaxation even for graphs with up to 1000 vertices, and a post-processing algorithm is applied to guarantee a valid lower bound. Using this algorithm, we can compute strong lower bounds for each k with reasonable computational effort.

3.2 A Heuristic for the Bisection Problem

The graph bisection problem can be written as a quadratic assignment problem (QAP). To do so, we set the weight matrix W to be the Laplacian matrix L and the distance matrix D to the 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

$$\begin{aligned} \min _{\pi \in \varPi _n} \ \sum _{i = 1}^n \sum _{j=1}^n W_{i,j} D_{\pi (i),\pi (j)} = \min _{\pi \in \varPi _n} \ \sum _{i = 1}^k \sum _{j=1}^{k} L_{\pi ^{-1}(i),\pi ^{-1}(j)} = kh_k. \end{aligned}$$

To compute an upper bound $u_k$ on $h_k$, we use a simulated annealing heuristic for the QAP, as introduced in [6], and divide the solution by k.

3.3 A Branch-and-bound Algorithm for the Bisection Problem

We implement an open-source B &B solver to solve graph bisection problems of medium size to optimality using the ingredients described in the previous sections, namely the SDP bound as described in Sect. 3.1 and the upper bound as described in Sect. 3.2.

We base our branching decision on the solution of the relaxation of the subproblem. Namely, we branch on the node with corresponding value in $x^1$ being closest to 0.5. It turns out that we can write the subproblems again as problems of a similar form. In particular, if we set a variable $x_i$ to be 0, we can write the problem as the minimization problem $\min \{\bar{x}^\top \bar{L} \bar{x} : e^\top \bar{x} = k\}$, where $\bar{x}$ is obtained from x by deleting $x_i$ and $\bar{L}$ by deleting the i-th row and column of L. The subproblem where we set $x_i = 1$ can be written as $\min \{\bar{x}^\top \widetilde{L} \bar{x} + c : e^\top \bar{x} = k - 1\}$, with $\bar{x}$ again resulting from x by deleting $x_i$ and $\widetilde{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 Sect. 3.1 and Sect. 3.2 to compute bounds.

3.4 Transformation to a Max-Cut Problem

A different approach to solving the graph bisection problem is to transform it to a max-cut problem and use a max-cut solver, e.g. the open source parallel solver from [11]. To do so, we first need to transform the bisection problem into a quadratic unconstrained binary problem (QUBO).

Lemma 2

Let $\tilde{x} \in \{0,1\}^n$ such that $e^\top \tilde{x} = k$, and denote $\mu _k = \tilde{x}^\top L\tilde{x}$. Then

$$h_k = \frac{1}{k}\cdot \min _{x \in \{0,1\}^n} \big \{x^\top (L + \mu _k E) x -2\mu _k k e^\top x + \mu _k k^2\big \}.$$

Proof

First note that $x^\top Lx + \mu _k \Vert e^\top x - k\Vert ^2 = x^\top (L + \mu _k ee^\top )x - 2\mu _k k e^\top x + \mu _k k^2$. Let $x\in \{0,1\}^n$. For $e^\top x = k$ we have $x^\top Lx + \mu _k \Vert e^\top x - k\Vert ^2 = x^\top Lx$. And if $e^\top x \not = k$, then $x^\top Lx + \mu _k \Vert e^\top x - k\Vert ^2 > \mu _k$. Hence, for any infeasible $x \in \{0,1\}^n$, the objective is greater than the given upper bound $\mu _k$ and therefore the minimum can only be attained for $x \in \{0,1\}^n$ with $e^\top x = k$.    $\square$

It is well known that a QUBO problem can be reduced to a dense max-cut problem with one additional binary variable [cf. 5].

4 Split and Bound

We now assemble the tools developed in the previous section to compute the edge expansion of a graph by splitting the problem into $\lfloor \frac{n}{2} \rfloor$ 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.

4.1 Pre-elimination

The size k of the smaller set of the partition can theoretically be any value from 1 to $\lfloor \frac{n}{2}\rfloor$. 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 $\ell _k$ obtained by solving the SDP (6) in combination with the upper bound introduced in Sect. 3.2. We do not need to further consider values of k where the lower bound $\ell _k$ of the scaled bisection problem is already above an upper bound $u^*$ on the edge expansion. A pseudo-code of this pre-elimination step is given in Algorithm 1.

Algorithm 1

Pre-eliminate certain values of k

As it can be seen in Fig. 1, for a graph associated to a randomly generated 0/1 polytope and for a network graph, about 2/3 of the potential values of k can be excluded already by considering the cheap lower bound $\ell _k$.

Fig. 1.

Lower and upper bounds for each k.

We can further reduce the number of candidates for k by computing a tighter lower bound $\tilde{\ell }_k$ by solving the DNN relaxation (7) with additional cutting planes. Note that in our implementation we do not compute the tighter bound $\tilde{\ell }_k$ as part of the pre-elimination, since this bound is computed in the root node of the B &B tree in the algorithm from Sect. 3.3.

4.2 Computational Aspects

Stopping Exact Computations Earlier. For all values of k that are not excluded in the pre-elimination step, we have to compute the scaled bisection $h_k$. For some values of k, however, the optimum $h_k$ is greater than the threshold $u^*$. In that case, it makes no sense to compute $h_k$ but stop as soon as it is clear that we will not find the optimum with this choice of k.

Order of Considering Sub-Problems. If we find a smaller upper bound while computing $h_k$, this is also a better upper bound on the edge expansion. This affects all other open bisection problems since a better upper bound means that we can stop the B &B algorithms even earlier. Therefore, we do another 30 trials of simulated annealing for each k in $\mathcal {I}$. We expect the best further improvement on the upper bound for k with the smallest upper bound $u_k$ and therefore consider the subproblems in ascending order of their upper bound.

5 Numerical Results

All of our algorithms were written¹ in Julia². All computations were carried out on an AMD EPYC 7532 with 32 cores with 3.30 GHz and 1024 GB RAM. All max-cut problems were solved with the max-cut solver of BiqBin³. The SDPs to compute our cheap lower bounds $\ell _k$ are solved with MOSEK⁴ and MINLPs are solved with Gurobi⁵ using JuMP⁶.

5.1 Benchmark Instances

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 according to the uniform model in [16]. For any pair $(d,n_d)$ with $n_8 \in \{164, 189\}$, $n_9 \in \{153,178,203,228,253,278\}$, and $n_ {10} \in \{256, 281\}$, we generated 3 instances. Another class of polytopes we consider are the grlex and grevlex polytopes introduced and investigated in [10]. 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 [4] with at most 500 vertices. Additionally, we consider some more network graphs obtained from an online network repository⁷.

5.2 Discussion of the Experiments

Our numerical experiments indicate that the variant using BiqBin demonstrates superior performance compared to the B &B algorithm for bisection. For example, on the instance chesapeake from the DIMACS set, it took 2 s to compute the edge expansion compared to 9.8 s with the tailored B &B algorithm. Therefore, and due to space restrictions, we only report the results using BiqBin to solve the scaled bisection problems.

The detailed results are given in Table 2–3. In each of these tables, the first three columns give the name of the instance and the number of vertices and edges. In the columns 4–6 we report the optimal solution, i.e., the edge expansion of the graph, and the global lower and upper bound after the pre-elimination. The number of candidates for k after the pre-elimination is given in column 7. Column 8 lists the total number of B &B nodes in the max-cut algorithm for all values of k considered. The last two columns display the time spent in the pre-elimination and the total time (including pre-elimination) of the algorithm.

As reported in the tables, the pre-elimination phase only leaves a comparably small number of candidates for k to be further investigated. 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 B &B tree is of high quality: for many of the instances the gap is closed within the root node. This holds for all instances where the number of B &B nodes coincides with $|\mathcal {I}|$. As for comparing the run times of an instance for different values of k, no general statement can be derived. Typically, for k that is around the size where the optimum is attained and for k close to $\frac{n}{2}$ we experience the longest run time. 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, cf. columns titled h(G) and $u^*$. Overall, we solve almost all of the considered instances within a few minutes, for very few instances the B &B tree grows rather large and therefore computation times exceed several hours.

Table 2. Results of split & bound for graphs of random 0/1, grlex and grevlex polytopes.

Table 3. Results of split & bound for network instances.

6 Summary and Future Research

We developed a split & bound algorithm to compute the edge expansion of a graph. The splitting refers to separately considering different sizes k of the smaller partition. 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 B &B algorithm that solves the problem for k fixed. Through numerical results on various classes of graphs, we demonstrate that our algorithm outperforms other existing methods like the exact solver of [19] reporting an average run time of 2.7 h for instances with 60 vertices.

In some applications, one wants to check whether a certain value is a lower bound on the edge expansion, e.g., the Mihail-Vazirani conjecture. This verification is a straightforward modification of our algorithm and we are currently working on an implementation that enables this option. Another line of research is to replace the simulated annealing approach by a more sophisticated heuristic, e.g., in the spirit of the Goemans-Williamson rounding. This is necessary if one wants to obtain high-quality solutions for larger instances. We will also investigate convexification techniques by using recent results on fractional programming [12, 17] and on exploiting submodularity [3] of the cut function.