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
Computer Science & Systems Engineering Program, Federal University of Rio de Janeiro

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

Quantum Annealing via D-Wave

This notebook will give the first interaction with D-Wave’s Quantum Annealer. It will use the QUBO modeling problem introduced earlier and will define it using D-Wave’s package dimod, and then solve them using neal’s implementation of simulated annealing classicaly and D-Wave system package to use Quantum Annealing. We will also leverage the use of Networkx for network models/graphs.

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 + c_Q = \min_{x \in \{0,1 \}^n} x^\top Q x + c_Q \]

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

Example

Suppose we want to solve the following problem via QUBO

\[\begin{split} \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} \\ \mathbf{x} \in \{0,1 \}^{11} \end{split}\]

# 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 Dwave packages dimod and neal
import dimod
import neal
# 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 \(Q\) matrix, we write the problem

\[\begin{split} \min_{\mathbf{x}} \mathbf{c}^\top \mathbf{x}\\ s.t. \mathbf{A}\mathbf{x}=\mathbf{b} \\ \mathbf{x} \in \{0,1 \}^{11} \end{split}\]

as follows:

\[\begin{split} \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} \end{split}\]

Exploting the fact that \(x^2=x\) for \(x \in \{0,1\}\), we can make the linear terms appear in the diagonal of the \(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} ) \]

For this problem in particular, one can prove the the penalization factor is given by \(\rho > \sum_{i=1}^n |c_i|\), therefore we choose this bound + 1.

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))
cQ = rho*np.matmul(b.T,b)
print(Q)
print(cQ)

[[ -46    0    0   48   48   48    0   48   48   48   48]
 [   0  -44    0   48    0   48   48    0   48   48   48]
 [   0    0  -44    0   48    0   48   48   48   48   48]
 [  48   48    0  -92   48   96   48   48   96   96   96]
 [  48    0   48   48  -92   48   48   96   96   96   96]
 [  48   48    0   96   48  -92   48   48   96   96   96]
 [   0   48   48   48   48   48  -91   48   96   96   96]
 [  48    0   48   48   96   48   48  -92   96   96   96]
 [  48   48   48   96   96   96   96   96 -139  144  144]
 [  48   48   48   96   96   96   96   96  144 -138  144]
 [  48   48   48   96   96   96   96   96  144  144 -139]]
144

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=cQ)

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 = neal.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')

minimum energy: 5.0

minimum energy: 5.0

Now let’s solve this using Quantum Annealing!

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

!dwave ping

Using endpoint: https://cloud.dwavesys.com/sapi/
Using region: na-west-1

Using solver: Advantage_system4.1

Submitted problem ID: f2abd2ae-64f2-4203-812d-a2d700e9b27c

Wall clock time:
 * Solver definition fetch: 550.347 ms
 * Problem submit and results fetch: 1043.683 ms
 * Total: 1594.029 ms

QPU timing:
 * post_processing_overhead_time = 1.0 us
 * qpu_access_overhead_time = 0.0 us
 * qpu_access_time = 15860.54 us
 * qpu_anneal_time_per_sample = 20.0 us
 * qpu_delay_time_per_sample = 20.58 us
 * qpu_programming_time = 15783.56 us
 * qpu_readout_time_per_sample = 36.4 us
 * qpu_sampling_time = 76.98 us
 * total_post_processing_time = 1.0 us

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))

Advantage_system4.1

Number of qubits= 5627
Number of couplers= 40279

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)

{'timing': {'qpu_sampling_time': 175620.0, 'qpu_anneal_time_per_sample': 20.0, 'qpu_readout_time_per_sample': 135.04, 'qpu_access_time': 191383.56, 'qpu_access_overhead_time': 543.44, 'qpu_programming_time': 15763.56, 'qpu_delay_time_per_sample': 20.58, 'post_processing_overhead_time': 1.0, 'total_post_processing_time': 1.0}, 'problem_id': '429a98d2-cd92-49d3-9cb9-c6b4890dcd8e', 'embedding_context': {'embedding': {3: (3721, 423), 0: (394,), 1: (468,), 4: (3707, 559), 2: (3647, 498), 5: (3662,), 6: (438, 3571), 7: (409, 3646), 8: (3736, 453), 9: (3692, 544), 10: (3677,)}, 'chain_break_method': 'majority_vote', 'embedding_parameters': {}, 'chain_strength': 159.83033942960208}}

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')

minimum energy: 5.0

minimum energy: 5.0

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.