COMBINATORIAL OPTIMIZATION WITH D-WAVE OCEAN SDK

Solving Max-Cut Problems with D-Wave Quantum Annealing

Split the Network for Maximum Gain!

9 min readNov 24, 2024

Backstory

Imagine you’re a strategist for a Shōgun determined to weaken a network of clans across the region. While no single clan poses a threat, their web of alliances, varying in strength, could unite into a formidable opposition. Your mission is to disrupt these alliances by dividing the clans into two groups, eliminating the strongest connections to ensure they remain divided and unable to challenge the Shōgun’s rule!

The Max-Cut problem is a classic optimization problem in graph theory. Given a weighted graph, it seeks to partition vertices into two disjoint subsets to maximize the sum of the edge weights crossing between the two.

Example

We have vertices (nodes): V = {A, B, C, D} and edges with weights:

(A, B) = 1
(A, C) = 2
(B, D) = 3
(C, D) =4

We want to divide the vertices into two mutually exclusive groups to maximize the total weight of edges crossing between the two groups.

There are multiple possible groupings:

{A} and {B, C, D} → Cut edges (A, B) and (A, C) → Total weight: 1 + 2 = 3
{A, B} and {C, D} → Cut edges (A, C) and (B, D) → Total weight: 2 + 3 = 5
{A, B, C} and {D} → Cut edges (B, D) and (C, D) → Total weight: 3 + 4 = 7
{A, C} and {B, D} → Cut edges (A, B) and (C, D) → Total weight: 1 + 4 = 5
{A, D} and {B, C} → Cut edges (A, B), (A, C), (B, D), and (C, D) → Total weight: 1 + 2 + 3 + 4 = 10
so on …

In this example, we only have 4 vertices, so the total possible combinations are 2⁴ = 16 patterns. However, not all patterns are distinct partitions in Max-Cut because the subsets are unordered. The actual number of distinct cuts is 8. We’ll see an example later to clarify this point.

In any case, we don’t want to check all the permutations when the graph is larger and more complicated. It is an NP-hard problem because finding an exact solution requires exploring all possible partitions of the graph, and the number of possible partitions grows exponentially with the number of vertices.

Mathematical Definition

Variables

Let G = (V, E) be an undirected graph where:

V is the set of vertices.
E is the set of edges.
Each edge ( i, j ) ∈ E has weight w[i, j] ≥ 0 (0 means no edge).

Note: Since G is undirected, we treat ( i, j ) and ( j, i ) are the same thing. Also, we don’t consider ( i, i ) an edge.

Let x[i] ∈ {-1, 1} be a variable for vertex i ∈ V:

x[i] = 1 if vertex i is in subset S₁
x[i] = -1 if vertex i is in subset S₂

where S₁ ⋂ S₂ = ∅ (the subsets are disjoint). Therefore, each vertex can belong to either S₁ or S₂.

Objective Function

The goal is to partition the vertex set V into S₁ and S₂ such that the sum of weights of edges across two subsets is maximized:

Note: ( i, j ) ∈ E means we use each edge once.

Instead of the edge-based summation, we can use the vertex-based summation as below to achieve the same purpose:

Either of the above expressions effectively sums the weights of the edges crossing the cut because:

QUBO Formulation

We transform the variables and objective function to convert the max-cut problem into a QUBO (Quadratic Unconstrained Binary Optimization) formulation.

Convert binary variables from x[i] ∈ {-1, 1} to x[i] ∈ {0, 1}:

Press enter or click to view image in full size

Substituting this into x[i]x[j] in the objective function:

Changing the maximization to minimization objective:

Equivalently, in the vertex-sum:

This is a quadratic formula where:

-x[i] and -x[j] are linear terms (single binary variables)
2x[i]x[j] are quadratic terms (product of two binary variables)

So, we can express the minimization objective using an n ⨉ n symmetric matrix Q and vector x:

Diagonal Elements (i = j):

The diagonal element of xᵀQx is Q[i, i] x[i]², but x[i]² = x[i] since x[i] is either 0 or 1. Therefore, the diagonal element Q[i, i] represents the linear terms in the objective function.

For vertex i, we take the sum of all weights connected to it (and scale it by −1):

Here, the condition in the sum is j ≠ i (not j > i) because the contribution to i comes from -x[i] and -x[j].

Off-diagonal Elements (i ≠ j):

In the mathematical formulation, the coefficient for the quadratic term x[i]x[j] is 2w[i, j], as we consider each pair (i, j) only when j > i. But in the symmetric matrix representation, we distribute this coefficient equally between Q[i, j] and Q[j, i], such that:

Off-diagonal coefficients

Final QUBO Matrix:

Revisiting the Example

Let V be {A, B, C, D}:

Sorting out by linear and quadratic terms:

Substituting the weights:

Making it into a symmetric QUBO formulation:

Environment Setup

We’ll use D-Wave’s Ocean SDK to solve the combinatorial optimization.

Install Required Libraries

I’d recommend setting up a local environment for this project using venv or conda first. The below is for venv:

# create a project folder
mkdir ocean_test
cd ocean_test

# setup and activate a virtual environment 
python3 -m venv venv
source venv/bin/activate# upgrade to obtain the latest libraries
pip install --upgrade pip

Then, install the following Python libraries:

pip install dimod dwave-ocean-sdk matplotlib pandas

I’ll use Jupyter notebook (.ipynb) to write the following program.

Defining the Problem

We define the problem with edge weights:

# Define the edge weights of the graph
edge_weights = {
    ('A', 'B'): 1,
    ('A', 'C'): 2,
    ('B', 'D'): 3,
    ('C', 'D'): 4
}

Solution 1: Constrained Quadratic Model (CQM)

CQM is intended for constrained problems, and the Max-Cut is inherently unconstrained after QUBO transformation. So, it is not strictly necessary, but we’ll demonstrate it here anyway.

Setting Up Variables

We set up binary variables using a set of all vertices (nodes):

import dimod

# A set of all nodes in the graph
nodes = set([node for edge in edge_weights for node in edge])

# Create binary variables for each node
# (A dictionary of binary variables)
vars = {node: dimod.Binary(node) for node in nodes}

Setting Up Objective

This sets up the max-cut objective defined as a maximization problem:

# Define the Max-Cut objective
objective = 0
for (n1, n2), weight in edge_weights.items():
    # Add the edge contribution to the objective
    objective += weight * (vars[n1] + vars[n2] - 2 * vars[n1] * vars[n2])

(n1, n2) means an edge in the same way ( i, j ) is in the above discussion. So, we are handling each edge only once, as defined in edge_weights.

Initialize CQM

We initialize the CQM object with the minimization objective:

# Initialize a CQM
cqm = dimod.ConstrainedQuadraticModel()

# Add the objective to the CQM (maximize cut)
cqm.set_objective(-objective)  # Negate since ExactCQMSolver minimizes

Solving the CQM

Since it’s a small problem, we’ll use ExactCQMSolver:

# Solve the CQM
solver = dimod.ExactCQMSolver()
solutions = solver.sample_cqm(cqm)

ExactCQMSolver checks all possible grouping patterns. So, it is only for small problems like this one.

Examining the Solutions

Let’s examine all 16 patterns:

def display_all_max_cut_solutions(solutions):
    # Extract and display all solutions sorted by energy
    sorted_solutions = solutions.record
    sorted_solutions.sort(order='energy')  # Sort by energy

    # Display solutions
    print("All Solutions Sorted by Energy:")
    for sol in sorted_solutions:
        variables = dict(zip(cqm.variables, sol.sample))
        energy = sol.energy
        print(f"Solution: {variables}, Energy: {energy}")

display_all_max_cut_solutions(solutions)

So, we see two patterns yielding the minimum energy:

x_A: 0, x_B: 1, x_C: 1, x_D: 0
x_A: 1, x_B: 0, x_C: 0, x_D: 1

Both cases mean the same grouping of {A, D} and {B, C}. So, either solution works.

Solution 2: Binary Quadratic Model (BQM)

Setting Up BQM

Let’s solve the same problem using BQM:

# Solve the Max-Cut problem using the D-Wave dimod package
import dimod

# Create a Binary Quadratic Model
bqm = dimod.BinaryQuadraticModel('BINARY')

# Add terms to the BQM for Max-Cut
for (n1, n2), weight in edge_weights.items():
    # Linear terms (diagonal)
    bqm.add_linear(n1, -weight)
    bqm.add_linear(n2, -weight)

    # Quadratic interaction term
    bqm.add_interaction(n1, n2, 2 * weight)

So, we are adding the linear and quadratic contributions, as we discussed before. BQM automatically handles symmetry, so we don’t need to explicitly add both (n1, n2) and (n2, n1).

ExactSolver

We can use ExactSolver, and it gives the same solution:

# Solve the BQM
solver = dimod.ExactSolver()
solutions = solver.sample(bqm)

display_all_max_cut_solutions(solutions)

The result is the same as CQM.

SimulatedAnnealingSampler

Let’s use SimulatedAnnealingSampler, which starts with a completely random initial state (bit patterns) and iteratively updates variables (randomly flips one at a time) using the Metropolis-Hastings criterion over decreasing temperatures, mimicking a physical annealing process to minimize the energy.

It accepts a new state (pattern) if it has a lower energy than before.
Otherwise, the new state is accepted with probability, given in the formula above.

Note: SimulatedAnnealingSampler is a general-purpose classical solver with nothing to do with D-Wave quantum hardware.

We can execute the sampler as follows:

# Solve the BQM using simulated annealing
sampler = dimod.SimulatedAnnealingSampler()

# Sample the BQM
sampleset = sampler.sample(bqm, num_reads=5)

It does not check all permutations. I set the num_reads = 5, and it finds the best solution in the first read.

Each read updates temperature updates for a specific number of time, which is by default 1000 but you can specify it using num_sweeps as below:

# Sample the BQM
sampleset = sampler.sample(bqm, num_reads=1, num_sweeps=10)

display_all_max_cut_solutions(sampleset)

The temperature is specified in the beta_range, which defines the inverse temperature range [βstart, βend], where β = 1/T.

The sample method has more arguments that you can tweak. For details, see Ocean documentation.

Direct QPU for Quantum Annealing

It’s unnecessary, but you can solve it using D-Wave’s QPU (Quantum Processing Unit) for quantum annealing if you sign up with D-Wave Leap, which gives you a limited access time at no charge.

from dwave.system import DWaveSampler, EmbeddingComposite

# Use the D-Wave sampler to find the best solution
sampler = EmbeddingComposite(DWaveSampler())

# Run the sampler on QPU
sampleset = sampler.sample(bqm, num_reads=1)

Note: num_sweeps does not apply to quantum annealing. The QPU hardware determines the annealing schedule. You can specify annealing_time to control the total duration.

We can confirm the direct QPU usage as shown below: