Solving Max-Cut with QAOA in Qiskit: A Practical Guide

What Is Max-Cut?

The Max-Cut problem asks: given a graph with weighted edges, how do you partition the vertices into two sets such that the total weight of edges crossing the partition is maximized? It is NP-hard in general, making it a natural target for quantum optimization algorithms.

Max-Cut is also a canonical benchmark for QAOA (Quantum Approximate Optimization Algorithm). It maps naturally to an Ising Hamiltonian, the cost function is easy to evaluate, and classical approximate algorithms (like Goemans-Williamson) provide a well-known comparison baseline.

This tutorial builds a complete solver using qiskit-optimization, from graph construction through result visualization. You will see why QAOA approximates rather than solves exactly, how circuit depth affects solution quality, and how to use the classical exact solver as a ground-truth reference.

Setting Up

You need the following packages:

pip install qiskit qiskit-optimization qiskit-algorithms networkx matplotlib scipy

Import everything up front:

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from qiskit_optimization.applications import Maxcut
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit_algorithms import QAOA, NumPyMinimumEigensolver
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import StatevectorSampler as Sampler

Constructing the Graph

Start with a small weighted graph. Seven nodes is enough to see interesting behavior while keeping circuit depth manageable.

# Requires: qiskit_optimization
def build_weighted_graph():
    """Create a weighted graph for Max-Cut experiments."""
    G = nx.Graph()
    
    # Add nodes
    n_nodes = 7
    G.add_nodes_from(range(n_nodes))
    
    # Add weighted edges
    edges = [
        (0, 1, 1.0),
        (0, 2, 2.0),
        (1, 3, 1.5),
        (2, 3, 1.0),
        (2, 4, 2.5),
        (3, 5, 1.0),
        (4, 5, 1.5),
        (4, 6, 2.0),
        (5, 6, 1.0),
        (1, 4, 0.5),
    ]
    G.add_weighted_edges_from(edges)
    return G


def visualize_graph(G, partition=None, title="Graph"):
    """Draw the graph, optionally coloring nodes by partition."""
    pos = nx.spring_layout(G, seed=42)
    fig, ax = plt.subplots(figsize=(8, 6))
    
    if partition is not None:
        colors = ['steelblue' if partition[i] == 0 else 'coral' 
                  for i in range(G.number_of_nodes())]
    else:
        colors = 'steelblue'
    
    edge_weights = nx.get_edge_attributes(G, 'weight')
    
    nx.draw_networkx(G, pos, node_color=colors, node_size=600,
                     font_color='white', font_weight='bold', ax=ax)
    nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_weights, ax=ax)
    
    if partition is not None:
        cut_value = sum(d['weight'] for u, v, d in G.edges(data=True)
                        if partition[u] != partition[v])
        ax.set_title(f"{title} (Cut = {cut_value:.1f})")
    else:
        ax.set_title(title)
    
    plt.tight_layout()
    plt.savefig(f"maxcut_{title.lower().replace(' ', '_')}.png", dpi=150)
    plt.show()


G = build_weighted_graph()
visualize_graph(G, title="Input Graph")
print(f"Graph: {G.number_of_nodes()} nodes, {G.number_of_edges()} edges")

Formulating as a QuadraticProgram

Qiskit Optimization provides the Maxcut application class that converts a NetworkX graph directly into a QuadraticProgram. The decision variable for each node is binary: 0 for one partition, 1 for the other.

# Requires: qiskit_optimization
maxcut_instance = Maxcut(G)
qp = maxcut_instance.to_quadratic_program()
print(qp.export_as_lp_string())

The LP output shows the objective: maximize the sum of w(i,j) * (x_i + x_j - 2x_ix_j) over all edges. For binary variables, this equals the edge weight when nodes i and j are in different partitions and zero when they are in the same partition. That is exactly the cut value.

Classical Exact Solution

Before running QAOA, get the optimal solution classically. This is your benchmark.

# Requires: qiskit_optimization
exact_solver = MinimumEigenOptimizer(NumPyMinimumEigensolver())
exact_result = exact_solver.solve(qp)

print("=== Exact Solution ===")
print(f"Optimal cut value: {-exact_result.fval:.4f}")  # negated because we minimized
print(f"Partition: {exact_result.x}")

# Decode partition
optimal_partition = {i: int(exact_result.x[i]) for i in range(G.number_of_nodes())}
visualize_graph(G, partition=optimal_partition, title="Optimal Cut")

Note: Qiskit Optimization minimizes by default, so the Max-Cut is formulated as a minimization of the negated cut value.

QAOA at p=1

QAOA uses two parameterized layers: a cost layer encoding the problem Hamiltonian and a mixing layer of X rotations. The number of repetitions p controls the approximation quality.

# Requires: qiskit_optimization
def run_qaoa(qp, p_layers, initial_point=None, shots=4096):
    """Run QAOA with p layers and return the OptimizationResult."""
    sampler = Sampler()
    optimizer = COBYLA(maxiter=300)
    
    if initial_point is None:
        rng = np.random.default_rng(42)
        initial_point = rng.uniform(0, np.pi, 2 * p_layers)
    
    qaoa = QAOA(
        sampler=sampler,
        optimizer=optimizer,
        reps=p_layers,
        initial_point=initial_point,
    )
    
    solver = MinimumEigenOptimizer(qaoa)
    result = solver.solve(qp)
    return result


# p=1
print("Running QAOA p=1...")
result_p1 = run_qaoa(qp, p_layers=1)
partition_p1 = {i: int(result_p1.x[i]) for i in range(G.number_of_nodes())}
cut_p1 = sum(d['weight'] for u, v, d in G.edges(data=True)
             if partition_p1[u] != partition_p1[v])
print(f"p=1 cut value: {cut_p1:.4f}")
visualize_graph(G, partition=partition_p1, title="QAOA p=1")

Comparing p=1, p=2, and p=3

The approximation ratio r = QAOA_cut / optimal_cut measures how close QAOA gets to optimal. Theoretically, p=1 QAOA guarantees r >= 0.6924 for unweighted Max-Cut (the Farhi-Goldstone-Gutmann bound). In practice, weighted problems and finite shots reduce this.

# Requires: qiskit_optimization
def evaluate_all_p_levels(qp, G, exact_result):
    """Compare QAOA across p levels."""
    optimal_cut = -exact_result.fval
    results = {}

    for p in [1, 2, 3]:
        print(f"\nRunning QAOA p={p}...")
        res = run_qaoa(qp, p_layers=p)
        partition = {i: int(res.x[i]) for i in range(G.number_of_nodes())}
        cut = sum(d['weight'] for u, v, d in G.edges(data=True)
                  if partition[u] != partition[v])
        ratio = cut / optimal_cut
        
        results[p] = {
            'cut': cut,
            'ratio': ratio,
            'partition': partition,
        }
        print(f"  Cut value: {cut:.4f}")
        print(f"  Approximation ratio: {ratio:.4f}")
        visualize_graph(G, partition=partition, title=f"QAOA p={p}")

    return results


results = evaluate_all_p_levels(qp, G, exact_result)
optimal_cut = -exact_result.fval

print("\n=== Summary ===")
print(f"Optimal cut value: {optimal_cut:.4f}")
for p, data in results.items():
    print(f"QAOA p={p}: cut={data['cut']:.4f}, ratio={data['ratio']:.4f}")

Visualizing QAOA Circuit Structure

Understanding what QAOA builds is important for reasoning about its depth and gate count:

# Requires: qiskit_optimization
from qiskit_algorithms import QAOA
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import StatevectorSampler as Sampler

# Build a p=2 QAOA circuit without running it
sampler = Sampler()
qaoa_p2 = QAOA(sampler=sampler, optimizer=COBYLA(), reps=2)

# Convert problem to Ising operator for circuit inspection
from qiskit_optimization.converters import QuadraticProgramToQubo
from qiskit_optimization.translators import to_ising

converter = QuadraticProgramToQubo()
qubo = converter.convert(qp)
operator, offset = to_ising(qubo)

# Get the ansatz circuit
qaoa_p2._check_operator_ansatz(operator)
circuit = qaoa_p2.ansatz
circuit = circuit.decompose()
print(f"Circuit depth: {circuit.depth()}")
print(f"Gate counts: {circuit.count_ops()}")
print(circuit.draw('text', fold=100))

Sampling the Full Solution Distribution

QAOA outputs a probability distribution over all 2^n bitstrings, not just one solution. Examining this distribution reveals the landscape:

# Requires: qiskit_optimization
from qiskit.primitives import StatevectorSampler as PrimitiveSampler
from qiskit_algorithms import QAOA
from qiskit_algorithms.optimizers import COBYLA
from qiskit_optimization.converters import QuadraticProgramToQubo
from qiskit_optimization.translators import to_ising

def get_qaoa_distribution(qp, G, p=2, optimal_params=None):
    """Get the full probability distribution from QAOA."""
    converter = QuadraticProgramToQubo()
    qubo = converter.convert(qp)
    operator, offset = to_ising(qubo)

    sampler = PrimitiveSampler()

    # Use pre-optimized parameters if provided
    if optimal_params is not None:
        gamma_beta = optimal_params
    else:
        rng = np.random.default_rng(0)
        gamma_beta = rng.uniform(0, np.pi, 2 * p)

    qaoa = QAOA(sampler=sampler, optimizer=COBYLA(), reps=p, initial_point=gamma_beta)

    # Access the sampler distribution directly
    qaoa._check_operator_ansatz(operator)
    bound_circuit = qaoa.ansatz.assign_parameters(gamma_beta)
    bound_circuit.measure_all()

    shots = 8192
    job = sampler.run([(bound_circuit,)], shots=shots)
    counts = job.result()[0].data.meas.get_counts()

    # Convert to cut values
    n = G.number_of_nodes()
    cut_distribution = {}
    for bitstring, count in counts.items():
        partition = [int(bitstring[-(i + 1)]) for i in range(n)]
        cut = sum(d['weight'] for u, v, d in G.edges(data=True)
                  if partition[u] != partition[v])
        prob = count / shots
        cut_distribution[round(cut, 2)] = cut_distribution.get(round(cut, 2), 0) + prob

    return cut_distribution


# Plot distribution of cut values
dist = get_qaoa_distribution(qp, G, p=2)
fig, ax = plt.subplots(figsize=(10, 5))
cuts = sorted(dist.keys())
probs = [dist[c] for c in cuts]
ax.bar(cuts, probs, width=0.1)
ax.axvline(x=optimal_cut, color='red', linestyle='--', label=f'Optimal ({optimal_cut:.1f})')
ax.set_xlabel("Cut Value")
ax.set_ylabel("Probability")
ax.set_title("QAOA p=2 Cut Value Distribution")
ax.legend()
plt.tight_layout()
plt.savefig("qaoa_distribution.png", dpi=150)
plt.show()

Interpreting the Results

Several practical points emerge from this experiment:

p=1 often achieves a good approximation ratio on small graphs, typically 0.70-0.85. The theoretical guarantee of 0.6924 is a worst-case lower bound.

Increasing p helps but with diminishing returns. Going from p=1 to p=2 often improves the ratio significantly. Going from p=2 to p=3 typically gives smaller gains while roughly tripling the two-qubit gate count.

Shot noise limits small-p gains. At p=1, the optimal landscape is relatively smooth and COBYLA converges reliably. At p=3, the parameter space is 6-dimensional and classical optimization becomes harder. More sophisticated optimizers like SPSA or gradient-based methods become necessary.

The problem size matters more than p for real hardware. On current NISQ hardware, circuit noise dominates for p > 3 on graphs with more than 20 nodes. Hardware-efficient approaches that map the problem graph to the device topology are essential for larger instances.

Max-Cut with QAOA remains an active research area. The computational advantage over classical approximation algorithms has not been demonstrated, but the framework is a template for applying QAOA to other combinatorial optimization problems: portfolio optimization, vehicle routing, and scheduling all admit similar formulations.