Quantum Annealing via D-Wave¶

David E. Bernal Neira
Davidson School of Chemical Engineering, Purdue University

Pedro Maciel Xavier
Davidson School of Chemical Engineering, Purdue University

Benjamin J. L. Murray
Davidson School of Chemical Engineering, Purdue University
Undergraduate Research Assistant

Quantum Annealing via D-Wave (Python)¶

This notebook gives a first interaction with D-Wave’s quantum annealer. It uses the QUBO modeling problem introduced earlier, defines it with D-Wave’s dimod package, solves it classically with dwave-samplers, and then connects to the D-Wave system tools for quantum annealing. We also leverage NetworkX for graph-based visualizations.

Environment and execution notes

For local execution from the repository root, run uv sync --group qubo before launching Jupyter so the Python D-Wave stack is available. Before running the quantum annealing cells locally, configure and verify your D-Wave credentials with uv run --group qubo dwave setup and uv run --group qubo dwave ping.

In Google Colab, the next code cell installs the D-Wave Ocean packages into the active Python runtime. Later, the quantum annealing section runs dwave setup and dwave ping to configure and verify access to your D-Wave account before calling DWaveSampler().

Problem statement¶

We define a QUBO as the following optimization problem:

\min_{x \in \{0,1 \}^n} \sum_{(ij) \in E(G)} Q_{ij} x_i x_j + \sum_{i \in V(G)} Q_{ii} x_i + \beta = \min_{x \in \{0,1 \}^n} \mathbf{x}^\top \mathbf{Q} \mathbf{x} + \beta

(1)

where we optimize over binary variables $x \in \{0,1\}^n$ on a constrained graph $G(V,E)$ defined by a weighted adjacency matrix $\mathbf{Q}$ . We also include an arbitrary offset $\beta$ .

Example¶

Suppose we want to solve the following problem via QUBO $$ \min_{\mathbf{x}} 2𝑥_0+4𝑥_1+4𝑥_2+4𝑥_3+4𝑥_4+4𝑥_5+5𝑥_6+4𝑥_7+5𝑥_8+6𝑥_9+5𝑥_{10} \ s.t. \begin{bmatrix} 1 & 0 & 0 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1\ 0 & 1 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 1 & 1\ 0 & 0 & 1 & 0 & 1 & 0 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}\mathbf{x}=

\begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix}

(2)

\mathbf{x} \in {0,1 }^{11} $$

# If using this on Google collab, we need to install the packages
try:
  import google.colab
  IN_COLAB = True
except:
  IN_COLAB = False

if IN_COLAB:
    !pip install dwave-ocean-sdk

# Import the D-Wave packages dimod and the local simulated annealer
import dimod
from dwave.samplers import SimulatedAnnealingSampler
# Import Matplotlib to generate plots
import matplotlib.pyplot as plt
# Import numpy and scipy for certain numerical calculations below
import numpy as np
from scipy.special import gamma
import math
from collections import Counter
import pandas as pd
from itertools import chain
import time
import networkx as nx

First we would write this problem as a an unconstrained one by penalizing the linear constraints as quadratics in the objective. Let’s first define the problem parameters

A = np.array([[1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1],
            [0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1],
            [0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1]])
b = np.array([1, 1, 1])
c = np.array([2, 4, 4, 4, 4, 4, 5, 4, 5,6, 5])

In order to define the $\mathbf{Q}$ matrix, we write the problem

\min_{\mathbf{x}} \mathbf{c}^\top \mathbf{x}\\ \text{s.t. } \mathbf{A}\mathbf{x}=\mathbf{b} \\ \mathbf{x} \in \{0,1 \}^{11}

(3)

as follows:

\min_{\mathbf{x}} \mathbf{c}^\top \mathbf{x} + \rho(\mathbf{A}\mathbf{x}-\mathbf{b})^\top (\mathbf{A}\mathbf{x}-\mathbf{b}) \\ \mathbf{x} \in \{0,1 \}^{11}

(4)

Exploiting the fact that $x^2=x$ for $x \in \{0,1\}$ , we can make the linear terms appear in the diagonal of the $\mathbf{Q}$ matrix.

\rho(\mathbf{A}\mathbf{x}-\mathbf{b})^\top (\mathbf{A}\mathbf{x}-\mathbf{b}) = \rho( \mathbf{x}^\top (\mathbf{A}^\top \mathbf{A}) \mathbf{x} - 2\mathbf{b}^\top \mathbf{A} \mathbf{x} + \mathbf{b}^\top \mathbf{b} )

(5)

For this problem in particular, one can prove that a reasonable penalization factor is given by $\rho = \sum_{i=1}^n |c_i| + \epsilon$ with $\epsilon > 0$ .

epsilon = 1
rho = np.sum(np.abs(c)) + epsilon
Q = rho*np.matmul(A.T,A)
Q += np.diag(c)
Q -= rho*2*np.diag(np.matmul(b.T,A))
beta = rho*np.matmul(b.T,b)
print(Q)
print(beta)

We can visualize the graph that defines this instance using the Q matrix as the adjacency matrix of a graph.

G = nx.from_numpy_array(Q)
nx.draw(G, with_labels=True)

Let’s define a QUBO model and then solve it via simulated annealing.

model = dimod.BinaryQuadraticModel.from_qubo(Q, offset=beta)

def plot_enumerate(results, title=None):

    plt.figure()

    energies = [datum.energy for datum in results.data(
        ['energy'], sorted_by=None)]
    
    if results.vartype == 'Vartype.BINARY':
        samples = [''.join(c for c in str(datum.sample.values()).strip(
            ', ') if c.isdigit()) for datum in results.data(['sample'], sorted_by=None)]
        plt.xlabel('bitstring for solution')
    else:
        samples = np.arange(len(energies))
        plt.xlabel('solution')

    plt.bar(samples,energies)
    plt.xticks(rotation=90)
    plt.ylabel('Energy')
    plt.title(str(title))
    print("minimum energy:", min(energies))


def plot_energies(results, title=None):
    energies = results.data_vectors['energy']
    occurrences = results.data_vectors['num_occurrences']
    counts = Counter(energies)
    total = sum(occurrences)
    counts = {}
    for index, energy in enumerate(energies):
        if energy in counts.keys():
            counts[energy] += occurrences[index]
        else:
            counts[energy] = occurrences[index]
    for key in counts:
        counts[key] /= total
    df = pd.DataFrame.from_dict(counts, orient='index').sort_index()
    df.plot(kind='bar', legend=None)

    plt.xlabel('Energy')
    plt.ylabel('Probabilities')
    plt.title(str(title))
    plt.show()
    print("minimum energy:", min(energies))

Let’s now solve this problem using Simulated Annealing

simAnnSampler = SimulatedAnnealingSampler()
simAnnSamples = simAnnSampler.sample(model, num_reads=1000)

plot_enumerate(simAnnSamples, title='Simulated annealing in default parameters')
plot_energies(simAnnSamples, title='Simulated annealing in default parameters')

Now let’s solve this using Quantum Annealing!¶

Before running the next cells, configure access to your D-Wave account. For local execution, run uv run --group qubo dwave setup and uv run --group qubo dwave ping from the repository root. In Colab, use the setup and ping cells below before calling DWaveSampler().

# Let's setup the D-Wave connection
if IN_COLAB:
    !dwave setup

!dwave ping

import dwave_networkx as dnx
from dwave.system import (DWaveSampler, EmbeddingComposite,
                          FixedEmbeddingComposite)
from pprint import pprint

# Graph corresponding to D-Wave Model
qpu = DWaveSampler()
qpu_edges = qpu.edgelist
qpu_nodes = qpu.nodelist
# pprint(dir(qpu))
if qpu.solver.id == "DW_2000Q_6":
    print(qpu.solver.id)
    X = dnx.chimera_graph(16, node_list=qpu_nodes, edge_list=qpu_edges)
    dnx.draw_chimera(X, node_size=1)
    print('Number of qubits=', len(qpu_nodes))
    print('Number of couplers=', len(qpu_edges))
else:
    print(qpu.solver.id)
    X = dnx.pegasus_graph(16, node_list=qpu_nodes, edge_list=qpu_edges)
    dnx.draw_pegasus(X, node_size=1)
    print('Number of qubits=', len(qpu_nodes))
    print('Number of couplers=', len(qpu_edges))

DWavesampler = EmbeddingComposite(DWaveSampler())
DWaveSamples = DWavesampler.sample(bqm=model, num_reads=1000, 
                                   return_embedding=True, 
                                  #  chain_strength=chain_strength, 
                                  #  annealing_time=annealing_time
                                   )

print(DWaveSamples.info)

embedding = DWaveSamples.info['embedding_context']['embedding']
if qpu.solver.id == "DW_2000Q_6":
  dnx.draw_chimera_embedding(X, embedding, node_size=2)
else:
  dnx.draw_pegasus_embedding(X, embedding, node_size=2)

plot_enumerate(DWaveSamples, title='Quantum annealing in default parameters')
plot_energies(DWaveSamples, title='Quantum annealing in default parameters')

Now we can play with the other parameters such as Annealing time, chain strength, and annealing schedule to improve the performance of D-Wave’s Quantum Annealing.

Quantum Annealing via D-Wave (Python)

Quantum Annealing via D-Wave¶