Quantum Annealing in Python — Deep Dive

Ocean SDK Architecture

D-Wave’s Ocean toolkit is organized in layers:

dimod          → Problem definitions (BQM, CQM, DQM)
dwave-system   → Hardware interface and embedding
dwave-hybrid   → Hybrid classical-quantum workflows
dwave-samplers → Classical samplers for development
dwave-cloud    → Cloud API client

Installation and Setup

pip install dwave-ocean-sdk

# Configure API access
dwave config create
# Enter your API token from cloud.dwavesys.com

Problem Formulation

Binary Quadratic Model (BQM)

import dimod

# Direct QUBO construction
Q = {
    ('x1', 'x1'): -1,    # Linear bias for x1
    ('x2', 'x2'): -1,    # Linear bias for x2
    ('x1', 'x2'): 2,     # Quadratic coupling
}

bqm = dimod.BinaryQuadraticModel.from_qubo(Q)

# Or build incrementally
bqm = dimod.BinaryQuadraticModel('BINARY')
bqm.add_variable('x1', -1.0)
bqm.add_variable('x2', -1.0)
bqm.add_interaction('x1', 'x2', 2.0)

Constrained Quadratic Model (CQM)

More natural for real-world problems:

from dimod import ConstrainedQuadraticModel, Binary, Integer

cqm = ConstrainedQuadraticModel()

# Variables
x = [Binary(f'x_{i}') for i in range(10)]
y = Integer('y', lower_bound=0, upper_bound=100)

# Objective: minimize total cost
costs = [3, 7, 2, 5, 8, 1, 4, 6, 9, 3]
cqm.set_objective(sum(c * xi for c, xi in zip(costs, x)))

# Constraint: select exactly 4 items
cqm.add_constraint(sum(x) == 4, label='select_four')

# Constraint: budget limit
weights = [10, 20, 5, 15, 25, 8, 12, 18, 22, 7]
cqm.add_constraint(sum(w * xi for w, xi in zip(weights, x)) <= 50,
                    label='budget')

Discrete Quadratic Model (DQM)

For multi-valued variables:

from dimod import DiscreteQuadraticModel

dqm = DiscreteQuadraticModel()

# Variable with 5 possible values (e.g., color assignment)
dqm.add_variable(5, label='node_0')
dqm.add_variable(5, label='node_1')
dqm.add_variable(5, label='node_2')

# Penalize adjacent nodes having the same color
for color in range(5):
    dqm.set_quadratic_case('node_0', color, 'node_1', color, 10.0)
    dqm.set_quadratic_case('node_1', color, 'node_2', color, 10.0)

Solving Problems

Classical Samplers (Development)

# Exact solver (brute force, small problems only)
sampler = dimod.ExactSolver()
sampleset = sampler.sample(bqm)
print(sampleset.first.sample)  # Optimal solution
print(sampleset.first.energy)  # Optimal energy

# Simulated annealing (classical, any size)
sampler = dimod.SimulatedAnnealingSampler()
sampleset = sampler.sample(bqm, num_reads=1000)

Quantum Annealer (Hardware)

from dwave.system import DWaveSampler, EmbeddingComposite

# Direct QPU access with automatic embedding
sampler = EmbeddingComposite(DWaveSampler())

sampleset = sampler.sample(bqm,
    num_reads=1000,           # Number of annealing runs
    annealing_time=20,        # Microseconds per anneal
    chain_strength=None,      # Auto-calculated
    return_embedding=True     # Include embedding info
)

# Analyze results
print(f"Best energy: {sampleset.first.energy}")
print(f"Best solution: {sampleset.first.sample}")
print(f"Timing: {sampleset.info['timing']}")

Hybrid Solver (Production)

from dwave.system import LeapHybridCQMSampler, LeapHybridBQMSampler

# CQM hybrid solver — handles constraints natively
sampler = LeapHybridCQMSampler()
sampleset = sampler.sample_cqm(cqm, time_limit=10)

# Filter feasible solutions
feasible = sampleset.filter(lambda s: s.is_feasible)
print(f"Best feasible solution: {feasible.first.sample}")
print(f"Objective value: {feasible.first.energy}")

Embedding: Logical to Physical

D-Wave’s qubits are connected in a specific topology (Pegasus graph). Not every logical variable maps directly to one physical qubit. Embedding maps logical variables to chains of physical qubits.

from dwave.system import DWaveSampler
from minorminer import find_embedding

# Get the hardware graph
qpu = DWaveSampler()
target_graph = qpu.edgelist

# Find embedding of problem graph into hardware graph
source_graph = list(bqm.quadratic.keys())
embedding = find_embedding(source_graph, target_graph)

print(f"Variable 'x1' mapped to chain: {embedding['x1']}")
# A chain of physical qubits acts as one logical variable

Chain Breaks

Chains can break — physical qubits in a chain disagree on their value. This is a major source of errors:

from dwave.embedding import chain_breaks

# Analyze chain break frequency
for sample, energy, num_occ, chain_break in sampleset.record:
    if chain_break:
        print(f"Chain break detected in sample with energy {energy}")

# Strategies to reduce chain breaks:
# 1. Increase chain strength (but too high distorts the problem)
# 2. Use longer annealing times
# 3. Use postprocessing to fix broken chains
sampleset = sampler.sample(bqm,
    num_reads=1000,
    chain_strength=2.0,           # Manual chain strength
    postprocessing='optimization'  # Fix chain breaks
)

Real-World Problem: Vehicle Routing

from dimod import ConstrainedQuadraticModel, Binary
import numpy as np

def build_vrp(n_locations, n_vehicles, distances, demands, capacity):
    """
    Vehicle Routing Problem:
    Assign locations to vehicles and order visits to minimize total distance.
    """
    cqm = ConstrainedQuadraticModel()

    # x[i][v] = 1 if location i is served by vehicle v
    x = {}
    for i in range(n_locations):
        for v in range(n_vehicles):
            x[i, v] = Binary(f'x_{i}_{v}')

    # Objective: minimize total distance (simplified)
    objective = 0
    for i in range(n_locations):
        for j in range(n_locations):
            if i != j:
                for v in range(n_vehicles):
                    objective += distances[i][j] * x[i, v] * x[j, v]
    cqm.set_objective(objective)

    # Each location served by exactly one vehicle
    for i in range(n_locations):
        cqm.add_constraint(
            sum(x[i, v] for v in range(n_vehicles)) == 1,
            label=f'serve_{i}'
        )

    # Vehicle capacity constraints
    for v in range(n_vehicles):
        cqm.add_constraint(
            sum(demands[i] * x[i, v] for i in range(n_locations)) <= capacity,
            label=f'capacity_{v}'
        )

    return cqm

# Solve
n_loc, n_veh = 15, 3
distances = np.random.randint(1, 50, size=(n_loc, n_loc))
np.fill_diagonal(distances, 0)
demands = np.random.randint(1, 10, size=n_loc)

cqm = build_vrp(n_loc, n_veh, distances, demands, capacity=30)

sampler = LeapHybridCQMSampler()
sampleset = sampler.sample_cqm(cqm, time_limit=10)
best = sampleset.filter(lambda s: s.is_feasible).first

Real-World Problem: Portfolio Optimization

def portfolio_optimization(returns, covariance, budget, risk_aversion=0.5):
    """
    Select assets to maximize return while minimizing risk.
    Classic Markowitz portfolio theory as QUBO.
    """
    n_assets = len(returns)

    bqm = dimod.BinaryQuadraticModel('BINARY')

    for i in range(n_assets):
        # Linear term: expected return (negative because we minimize)
        bqm.add_variable(f'asset_{i}', -returns[i])

        # Quadratic terms: risk (covariance)
        for j in range(i + 1, n_assets):
            bqm.add_interaction(
                f'asset_{i}', f'asset_{j}',
                risk_aversion * covariance[i][j]
            )
        # Diagonal risk term
        bqm.add_variable(f'asset_{i}', risk_aversion * covariance[i][i])

    # Budget constraint as penalty
    penalty_strength = max(abs(r) for r in returns) * 10
    budget_terms = [(f'asset_{i}', 1.0) for i in range(n_assets)]
    bqm.update(dimod.generators.combinations(
        [f'asset_{i}' for i in range(n_assets)],
        budget,
        strength=penalty_strength
    ))

    return bqm

Annealing Schedule Customization

Fine-tune the quantum annealing process:

# Custom annealing schedule
# Format: list of (time_fraction, s_value) pairs
# s goes from 0 (full quantum) to 1 (full classical)

# Standard linear anneal
schedule = [[0.0, 0.0], [20.0, 1.0]]

# Pause in the middle (helps with hard problems)
schedule = [[0.0, 0.0], [10.0, 0.4], [60.0, 0.4], [80.0, 1.0]]

# Reverse anneal (start from a known solution)
initial_state = {'x1': 1, 'x2': 0, 'x3': 1}
sampleset = sampler.sample(bqm,
    initial_state=initial_state,
    anneal_schedule=[[0.0, 1.0], [10.0, 0.3], [20.0, 1.0]],
    reinitialize_state=True
)

Reverse annealing starts from a classical solution and briefly increases quantum fluctuations to explore nearby better solutions — useful for iterative refinement.

Performance Benchmarking

import time

def benchmark_problem(bqm, n_runs=5):
    """Compare classical vs quantum solving times."""

    # Simulated annealing (classical)
    sa_sampler = dimod.SimulatedAnnealingSampler()
    sa_times, sa_energies = [], []
    for _ in range(n_runs):
        start = time.time()
        result = sa_sampler.sample(bqm, num_reads=1000)
        sa_times.append(time.time() - start)
        sa_energies.append(result.first.energy)

    # QPU (quantum)
    qpu_sampler = EmbeddingComposite(DWaveSampler())
    qpu_result = qpu_sampler.sample(bqm, num_reads=1000)
    qpu_timing = qpu_result.info['timing']

    print(f"SA: avg energy={np.mean(sa_energies):.2f}, "
          f"avg time={np.mean(sa_times)*1000:.1f}ms")
    print(f"QPU: best energy={qpu_result.first.energy:.2f}, "
          f"qpu_time={qpu_timing['qpu_access_time']/1000:.1f}ms, "
          f"total_time={qpu_timing['total_real_time']/1000:.1f}ms")

Note: QPU access time (microseconds) vs. total time (includes network latency and preprocessing). The quantum computation itself is blazingly fast — the overhead is everything around it.

Production Patterns

Solution Validation

def validate_solution(cqm, sampleset):
    """Validate solutions against all constraints."""
    feasible = sampleset.filter(lambda s: s.is_feasible)

    if len(feasible) == 0:
        print("WARNING: No feasible solutions found")
        print("Constraint violations:")
        best = sampleset.first
        for label in cqm.constraint_labels:
            if not cqm.check_feasible(best.sample, label):
                print(f"  VIOLATED: {label}")
        return None

    return feasible.first

Iterative Refinement

def iterative_solve(bqm, iterations=5):
    """Use reverse annealing to refine solutions."""
    sampler = EmbeddingComposite(DWaveSampler())

    # Initial forward anneal
    result = sampler.sample(bqm, num_reads=100)
    best = result.first

    for i in range(iterations):
        # Reverse anneal from best known solution
        result = sampler.sample(bqm,
            initial_state=best.sample,
            anneal_schedule=[[0.0, 1.0], [5.0, 0.3], [15.0, 1.0]],
            num_reads=100,
            reinitialize_state=True
        )
        if result.first.energy < best.energy:
            best = result.first
            print(f"Iteration {i}: improved to energy {best.energy}")

    return best

One thing to remember: D-Wave’s Ocean SDK abstracts away the quantum physics — your job is formulating the right QUBO or CQM that accurately represents your problem. The quality of your formulation matters far more than any hardware tuning, and the hybrid solvers handle most of the complexity of bridging classical and quantum processing.