Efficient Solutions for Large-Scale Max-Cut Problems: A Hybrid Local Search Heuristic Approach and Comparative Analysis with Quantum Annealing

1. Introduction

The Quadratic Unconstrained Binary Optimization (QUBO) problem holds significant relevance across various applications, particularly as a unifying approach for many combinatorial optimization challenges [1]. The emergence of quantum computing, particularly quantum annealing solvers, has brought QUBO into the spotlight [2]. However, the complexities of modelling and analyzing very large-scale problems on conventional computing platforms often hinder effective solutions through mathematical programming methods.

While quantum computing platforms offer the potential to overcome hardware limitations, transforming real-world problems into the suitable qubit architecture remains challenging. This study addresses the gap between quantum and traditional computing platforms, focusing on the very large-scale Max-Cut problem (MCP)—a critical research area with applications spanning statistical physics [3], communication infrastructure [3], machine learning [4], and computer chip design [5].

The computational landscape of the Max-Cut problem has sparked diverse solution strategies, ranging from approximation algorithms leveraging semidefinite programming to metaheuristic methods and exact approaches. While approximation algorithms like the seminal work of Goemans and Williamson [6] and subsequent contributions by Karish et al. offer theoretical performance guarantees [7], empirical evaluations often reveal their inferiority to other methodologies.

Exact methods, exemplified by Krishnan and Mitchell’s cut and price approach [8] and Rendl et al.’s branch and bound method [9], have emerged as viable alternatives, particularly for smaller problem instances. However, their applicability remains constrained to scenarios with relatively few nodes, typically in the hundreds.

For more substantial Max-Cut instances, characterized by thousands of nodes, metaheuristic methods become imperative. Recent advancements in this domain include the rank-2 relaxation method, dubbed CirCut, by Burer et al. [10], Festa et al.’s hybrid randomized method [11], and Marti et al.’s scatter search method [12]. Computational evaluations consistently demonstrate the superiority of these techniques, with the scatter search method notably achieving outstanding performance across test problems featuring up to 3000 nodes. Additionally, Alidaee et al. proposed a diversification strategy based on sequential local improvement, demonstrating substantial enhancements in solutions for a significant portion of benchmark Max-Cut problems within reasonable CPU time constraints [13].

Further breakthroughs have been made to tackle even larger problem sizes. Kochenberger et al. introduced the DDT approach, delivering state-of-the-art solutions for test problems boasting up to 10,000 nodes [14].

Max-Cut problem is a popular testbed for Quantum Computing platforms and algorithms due to its natural formulations. Recent developments in quantum computing, especially quantum annealing solvers, put Max-Cut in the spotlight [15,16,17,18,19,20,21,22,23,24]. Alam et al. used the graph Max-Cut problem as a prototype while using the quantum approximate optimization algorithm (QAOA) where a quantum circuit and a classical optimizer operate in a closed loop solving hard combinatorial optimization problems [15]. Basso et al. also applied QAOA to Max-Cut on large-girth regular graphs [16]. Galda et al. showed that the transferability of the parameters among QAOA instances and convergence of the optimal parameters around specific values can be predicted based on local graph properties [17]. Bianchi et al. proposed a pooling operator for graph neural networks to generate coarser graphs while maintaining the overall topology of the graph and used a spectral algorithm approximating the Max-Cut solution [18]. Shaydulin et al. utilized a multi-start optimization technique in a QAOA framework to improve the performance of quantum machines on various graph clustering problems [19] and introduced QAOA Kit for the combination and standardization of previously known parameters for Max-Cut problems [20]. Shaydulin and Wild proposed a new method to accelerate the assessment of QAOA energy by using the problem symmetry and implemented the proposed method on the Max-Cut problem using a graph auto-morphism solver and tensor network simulator [21]. Umasankar et al. used a Time-Multiplexed Opto-Electronic Oscillator-based Coherent Ising Machine on Max-Cut Problems [22]. Wurtz and Love defined the Spanning Tree QAOA for Max-Cut and utilizes an “ansatz” whose structure results from an approximate classical solution and reaches the same performance level as the classical algorithm to outperform QAOA at low depth [23].

In this research, we introduce a local search strategy based on the r-flip strategy for QUBO [24]. The dynamic adjustment of the r value in the r-flip strategy strikes a balance between computational cost and solution quality. This approach effectively reduces computational time while reaching superior solutions [25]. Our proposed heuristic achieves results within a strict CPU time limit of 600 seconds, addressing Max-Cut instances ranging from 10,000 to 40,000 variables. Comparative evaluations against a quantum annealing solver based on multi-start Tabu Search [26] provide insights into the efficiency and effectiveness of our proposed approach.

By bridging the gap between quantum and conventional computing, our study contributes to the exploration of quantum computing’s potential in solving complex optimization problems, showcasing a novel local search strategy tailored for QUBO with promising results in the context of the challenging Max-Cut problem.

2. QUBO and MCP Solution

Before we present the r-flip strategy to solve QUBO and MCP in this study, the notations used are given as follows:

n The number of variables

x An initial feasible solution

x* The best solution found by the algorithm

K The largest value of k for r-flip, k$\le$ r

π(i) The i-th element of x in the order π(1)⋯π(n)

S {i:xi is tentatively chosen to receive a new value to produce a new solution xi’} restricting consideration to |S| = r

D The set of candidates for an improving move

$E\left(x_i\right)$ Derivative of f(x) with respect to $x_i$

$E\left(x\right)=\left(E\right(x_1),\dots ,E(x_n\left)\right)$ The vector of derivatives

x(.) A vector representing the solution of x

E(.) A vector representing the value of derivative $E\left(x_i\right)$

The general form of QUBO can be given as [27],

$$Max f\left(x\right)=\sum _{i=1}^n{q}_i{x}_i+\frac{1}{2}\sum _{i=1}^n\sum _{j\ne i}^n{q}_{i,j}{x}_i{x}_j, s.t. x_i\in \left\{\mathrm{0,1}\right\}, i=1,\cdots ,n$$

(1)

In (1), $\frac{1}{2}{q}_{i,j}$ is the i,j-th entry of a given n by n symmetric matrix Q and the linear term $q_i$ can be the entry in the diagonal of the matrix. Since ${x_i}^2=x_i$, and Q may be written as an upper triangular matrix by doubling each entry of the upper triangle part of the matrix and letting $q_{i,i}=q_i$, then we can write (1) as (2).

$$Max f\left(x\right)={\sum }_{i=1}^n{\sum }_{j\ge i}^n{q}_{i,j}{x}_i{x}_j=x^{T}Qx, s.t. x\_i\in \left\{\mathrm{0,1}\right\},i=1,\cdots ,n$$

(2)

Numerous combinatorial optimization problems lend themselves to reformulation as QUBO by employing a quadratic infeasibility penalty. However, determining the optimal penalty values poses a considerable challenge. In light of this, we focus our investigation on the very large-scale MCP as a prime candidate application, given its inherent and natural formulation within the QUBO framework.

2.1. Max-Cut problems

Max-Cut problem with many real-world implementations is an NP-hard problem which means that exact solution algorithms are not sufficient for large-scale problems, since obtaining a solution is very time-consuming [28]. Following the advances in quantum computers, max-cut problems have received much interest recently [29]. Two formulations, linear and nonlinear, exist for the MCP [27]. Initially, we present the natural QUBO model for MCP.

$$Max {\sum }_{i=1}^n{\sum }_{j=1,i<j}^n{q}_{ij}\left(x_i+x_j-2\ast x_i{x}_j\right)$$

(3)

$$\mathrm{s}.\mathrm{t}. x_i \in \left\{\mathrm{0,1}\right\}\forall i=1,\dots ,n$$

The QUBO model (3) can be transformed to a linear model using classical linearization techniques [30].

$$Max {\sum }_{i=1}^n{\sum }_{j=1,i<j}^n{q}_{ij}\left(x_i+x_j-2\ast y_{ij}\right)$$

(4)

$$\mathrm{s}.\mathrm{t}. y_{ij}\le x_i, \forall i,j=1,\dots ,n and i<j$$

(5)

$$y_{ij}\le x_j, \forall i,j=1,\dots ,n and i<j$$

(6)

$$y_{ij}\ge x_i+x_j-1, \forall i,j=1,\dots ,n and i<j$$

(7)

$$x_i \in \left\{\mathrm{0,1}\right\}, \forall i=1,\dots ,n$$

$$y_{ij}\in [0,1], \forall j=1,\dots ,n and i<j$$

The linear model for the MCP undergoes a substantial increase in size on very large-scale instances, both in terms of the number of variables and constraints. In this research, we specifically present solutions for the QUBO model, solely for the purpose of conducting a comparative analysis with a quantum annealing solver.

2.2. Solution Representation

Here Q is a symmetric matrix. Given a solution $x=(x_1,\dots ,x_n)$ define

$${\Delta }_i\left(x\right)=q(i,i)+{\sum }_{k\ne i}2x_k{q}_{ik}, i\in N$$

(8)

With respect to 1-flip neighborhood, x is locally optimal if and only if it satisfies

$${\Delta }_i\left(x\right)\le 0, \forall x_i=0, and {\Delta }_i\left(x\right)>0, \forall x_i=1$$

(9)

Then the objective function will satisfy

$$∆f={\sum }_{i\in S}(x{\text{'}}_i-x_i)E\left(x_i\right)+{\sum }_{i,j\in S,i<j}(x{\text{'}}_i-x_i)(x{\text{'}}_j-x_j)q_{i,j}$$

(10)

and

$$\forall j\in S, E\left(x_j\right)\leftarrow E\left(x_j\right)+\sum _{i\in S\backslash \left\{j\right\}}\left(x_i^{\text{'}}-x_i\right)q_{i,j}$$

(11)

The rationale behind the r-flip move is rooted in the notion that a larger value of moves (here, a larger value of r) offers the opportunity to explore a more extensive and diverse solution space. Moreover, when improvement within a specific neighborhood (e.g., a 1-flip move) is exhausted, and a solution x attains local optimality with respect to the 1-flip move, employing another neighborhood move (such as a 2-flip move) becomes viable for further exploration of the solution space. This strategic approach forms the basis for various meta-heuristics. Upon obtaining a locally optimal solution x concerning a 1-flip move, the elements of x are ordered based on the ascending absolute value of derivatives, as follows:

$$\left|E\right(x_{\pi \left(1\right)}\left)\right|\le \cdots \le \left|E\right(x_{\pi \left(n\right)}\left)\right|$$

(12)

Here, $\pi \left(i\right)$ means the i-th element of x in the order $\left(\pi \right(1),\cdots ,\pi (n\left)\right)$. Now, one at a time in the given order, check the summation in (13) for k=1,2,…,n. For m<n, define (m, n) to be the number of combinations of m elements out of n, and let $\phi =\underset{i,j\in N}{\mathrm{Max}}\left\{\left|q_{i,j}\right|\right\}$, and for the r-flip move let $M=\phi \ast (2,r)$. Let K be the largest value of k where the inequality (13) holds, and define the set D(n) (14).

$${\sum }_{i=1}^k\left|E\right(x_{\pi \left(i\right)}\left)\right|<M, \mathrm{for} k=\mathrm{1,2},3,\cdots ,n$$

(13)

$$D\left(n\right)=\{x_{\pi \left(1\right)},\cdots ,x_{\pi \left(K\right)}\}$$

(14)

To implement an r-flip local search, the main decision is to choose a set S with |S|=r at each step and change $x_i$to $1-x_i$for all $i\in S$if it improves the objective function. This search process is inherently computationally demanding. The strategies mentioned above, including Lin and Kernighan's approach and very large-scale neighborhood search, are inherently tailored to specific problem domains. For instance, these strategies have predominantly found application in graph theoretic problems. However, in this paper, we introduce a novel implementation of the r-flip strategy, coupled with a sequential diversification technique. This approach enables the exploration of diverse sets S, each comprising r elements, while strategically avoiding non-improving sets. The incorporation of sequential diversification fundamentally reduces the computational time required for QUBO.

2.3. Implementation Details

The solution implementation is detailed in Section 2.2 through Algorithm 1 and Algorithm 2, with a comprehensive comparison of very large-scale MCP instances presented in the subsequent section. For additional insights and theoretical underpinnings, we draw upon the work and proofs by Alidaee and Wang [31], particularly for the novel r-flip results in general QUBO problems. It’s worth noting that in Algorithm 2, the sets S and M dynamically evolve as the search progresses.

Algorithm 1. 1-flip Local Search.

Note that, the algorithm is similar to the greedy algorithm for a KP. Here, the value that each time is added to the knapsack is equal to. We used a random sequence to greedily choose items, however, another alternative that many researchers have used is, each time choosingse the item with the maximum contribution to F(x). Furthermore, due to the dynamicsm nature of the problem, when an is changed to it is possible the algorithm will come back to the same variable and change its value to what it was before (this might happen several times during the search possible). Also, note that, reaching a locally optimal solution by the greedy Algorithm 1, is NP-hard, (…).

In a 1-flip search, we only look at one variable at a time to possibly change its value. However, in r-flip search the process is similar to 1-flip search except that up to r variables may be considered for a change of values. The idea behind the r-flip search is that, the larger value of moves (here larger value of r) may provide an opportunity to explore a larger and thus more diverse solution space. In this process when a neighborhood search, for example, k-flip ( ) search, is exhausted, i.e., x is locally optimal with respect to k-flip move, another neighborhood move, for example (k+1)-flip search, may be implemented to explore the solution space. Such a strategy is the basis for several meta-heuristics.

Algorithm 2. r-flip Local Search.

An alternative way to implement this algorithm is that k=r throughout the process. The Algorithm 2 is simple but inefficient; it requires some mechanism to choose the set S, and a process mechanism to implement elements of the k-flip move. Below, we give several results that allow an efficient implementation of the r-flip search process.

Lemma 1.

Given a locally optimal solution x with respect to a 1-flip search, we have

$$(x{\text{'}}_i-x_i)E_i(x)\le 0,\mathrm{for} i=1,\dots ,n$$

(15)

Proof. 

Directly follows from the local optimality condition (6), which states

$$\left(E_i(x)\ge 0,iff{x}_i=1\right) \mathrm{and} \left(E_i(x)\le 0,iff{x}_i=0\right),\mathrm{for} i=1,\dots ,n$$

(16)

Lemma 2.

Let $S\subseteq N$, and |S|=r, given any solution x, and M as defined earlier, we have

$${\displaystyle {\sum }_{i,j\in S}(x{\text{'}}_i-x_i)(x{\text{'}}_j-x_j)q_{ij}}\le M$$

(17)

Proof. 

For each pair of elements, i,j in S, the left-hand-side (LHS) in (17) is either $q_{ij}$ or $-q_{ij}$.Using definition of M, the summation in LHS is at most equal to M.

Theorem 3.

Given a locally optimal solution x of f(x) with respect to 1-flip search, a subset $S\subseteq N$, with |S|=r, is an improving r-flip move if and only if we have

$${\displaystyle {\sum }_{i\in S}|E_i(x)|}<{\displaystyle {\sum }_{i,j\in S}(x{\text{'}}_i-x_i)(x{\text{'}}_j-x_j)q_{ij}}$$

(18)

Proof. 

Using (10), a subset S with r elements is an improving r-flop move if and only if we have (19)

$${\displaystyle {\sum }_{i,j\in S}(x{\text{'}}_i-x_i)(x{\text{'}}_j-x_j)}{q}_{ij}>-{\displaystyle {\sum }_{i\in S}(x{\text{'}}_i-x_i)E_i(x)}={\displaystyle {\sum }_{i\in S}\left|E_i(x)\right|}$$

(19)

Since x is a locally optimal solution with respect to a 1-flip search, it follows from Lemma 1 that

inequality (19) is equivalent to (20); which completes the proof.

$${\displaystyle {\sum }_{i,j\in S}(x{\text{'}}_i-x_i)(x{\text{'}}_j-x_j)}{q}_{ij}>-{\displaystyle {\sum }_{i\in S}(x{\text{'}}_i-x_i)E_i(x)}={\displaystyle {\sum }_{i\in S}\left|E_i(x)\right|}$$

(20)

Proposition 2.

Let x be any locally optimal solution of f(x) with respect to a 1-flip search. If a subset $S\subseteq N$, with |S|=r, is an improving r-flip move, then we must have${\displaystyle {\sum }_{i\in S}\left|E_i(x)\right|}<M$.

Proof. 

The result follows from Theorem 3 and Lemma 2.

A consequence of Theorem 3 is as follows. Given a locally optimal solution x with respect to a 1-flip search, if there is no subset of S with |S|=r that satisfies (18), then x is also locally optimal solution with respect to an r-flip search. In fact, if there is no subset S of any size that (18) is satisfied, then x is also a locally optimal solution with respect to an r-flip search for all $r\le n$. Similar statements are true also regarding Proposition 2. The result of Proposition 2 is significant in the implementation of an r-flip search. It illustrates that, after having a 1-flip search implemented, if an r-flip search is next served as a locally optimal solution, only those elements with the sum of absolute value of derivatives less than M are eligible for consideration. Furthermore, when deciding about the elements of an r-flip search, we can easily check to see if any element $x_i$ by itself or with a combination of other elements is eligible to be a member of an improving r-flip move.

We adopted the following notation for computational results in next section:

BFS Best found solution among 10 runs within the CPU time limit.

APD Average percentage difference: the measurement for the performance of individual algorithm across runs on the same instance.

RPD Relative percentage deviation: the measurement for the relative performance across algorithms/methods on the same instance

3. Computational Experiments

In this investigation, we conducted extensive computational experiments to assess the efficacy of the proposed r-flip strategy across varying problem sizes, densities, and r values. The implementation of Algorithm 2 is executed using the GNU C++ compiler and undergoes 10 runs on the AMD Ryzen 5 3500U processor equipped with 16GB RAM. Simultaneously, the Quantum Annealing Solver (QAS), employing a multi-start tabu search, is implemented via a Python client on the same processor.

For the hybrid 1-flip/r-flip strategy in Algorithm 2 applied to very-large-scale MCP instances, we set a CPU time limit of 600 seconds. To comprehensively evaluate the r-flip strategy’s performance on very large-scale instances, we generate a diverse set of planar graph instances tailored for the MCP. This dataset comprises 120 instances, ranging from Q matrices sized 10,000 by 10,000 to 40,000 by 40,000. Among these instances, three distinct weight types with discrete values are introduced: Type ‘a’ with values of 1 and -1, Type ‘b’ with random discrete values between -10 and 10, and Type ‘c’ with random discrete values between -1,000 and 1,000. The instance data, ranging in size from 10,000 to 40,000, is conveniently accessible in the supplementary materials provided.

4. Numerical Results

Table 1, Table 2 and Table 3 present the outcomes obtained from both the hybrid 1-flip/r-flip Algorithm 2 and QAS across the 120 instances. It’s crucial to note that the results from the QAS are derived using the sampler approach (sampler=10). In reporting the QAS results, we focus on Relative Percentage Deviation (RPD) as the quality of solutions hinges on the parameter of sufficient reads (num_reads) and the collection of samples for potential solutions.

Concurrently, for Algorithm 2, we provide the Average Percentage Deviation (APD) value to illustrate the consistent performance over 10 runs. This dual reporting approach allows for a comprehensive evaluation of the comparative performance between the hybrid 1-flip/r-flip Algorithm 2 and the Quantum Annealing Solver across the diverse set of instances.

The r-flip strategy holds the potential for integration into various other local search heuristics as an enhancement procedure. A judicious implementation of the hybrid 1-flip/r-flip strategy not only diminishes computing time but also enhances solution quality. The outcomes reported here represent the best solutions obtained from 10 independent runs for each instance using the hybrid 1-flip/r-flip strategy. Given the constraints on computing resources in this study, the solution deviation across these 10 runs is calculated within the specified short CPU time limits. Figure 1, Figure 2 and Figure 3 illustrate the performance of Algorithm 2 vs. QAS for 10, 20, and 40 thousand variables.

5. Conclusions

In this paper, we explore the application of r-flip search within the context of a very large-scale Max-Cut problem, treating it as a case of QUBO. Specifically, our investigation centers on assessing the local optimality of the r-flip search when a 1-flip search has already reached its local optimum. Our computational findings demonstrate a significant reduction in the pool of candidates for potential r-flip implementation, showcasing the efficiency of this novel strategy in solving very large-scale Max-Cut instances within a strict 600-second timeframe.

In our approach, we employed a 1-flip strategy to achieve local optimality. Yet, it’s well recognized that numerous problems, particularly those in constrained optimization, can reap substantial benefits from r-flip strategies as local procedures. Comparisons with alternative methodologies documented in the literature underscore the remarkable performance of our approach, surpassing leading competitors by a significant margin. Notably, our method consistently discovered previously undiscovered optimal solutions for the majority of the problem instances tackled. What distinguishes our achievement is that our solution framework isn’t tailored specifically for the Max-Cut problem but rather accommodates a broader class of problems. Our findings strongly indicate the potential for significant advancements in r-flip solution techniques, hinting at their capacity to tackle even larger-scale Max-Cut challenges with further refinement.

These results present notable advantages, particularly in scenarios where variable neighborhood strategies are being employed for very large-scale problems or sparse matrices. Furthermore, we benchmark our outcomes against those obtained from a quantum annealing solver based on multi-start tabu search. Notably, the quantum annealing solver has also proven effective in addressing the challenges posed by very large-scale MCP. This comparative analysis provides insights into the strengths and applicability of both the r-flip strategy and quantum annealing solver in tackling complex optimization problems.