Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Environment Setup

Quantum Annealing via D-Wave

David E. Bernal Neira
Davidson School of Chemical Engineering, Purdue University
Universities Space Research Association
NASA QuAIL

Pedro Maciel Xavier
Davidson School of Chemical Engineering, Purdue University

Open In Colab SECQUOIA

Environment Setup

If you open this notebook in Google Colab, make sure the runtime is set to Julia before running the next cell. The setup cell will clone SECQUOIA/QuIP into the Colab workspace when needed, activate notebooks_jl/Project.toml, and install the Julia packages used in this notebook.

import Pkg

IN_COLAB = haskey(ENV, "COLAB_RELEASE_TAG") || haskey(ENV, "COLAB_JUPYTER_IP")

function detect_quip_repo_dir()
    candidates = (pwd(), normpath(pwd(), ".."))

    for candidate in candidates
        if isfile(joinpath(candidate, "notebooks_jl", "Project.toml"))
            return candidate
        end
    end

    error("Could not locate the QuIP repository root from $(pwd()).")
end

QUIP_REPO_DIR = if IN_COLAB
    repo_dir = get(ENV, "QUIP_REPO_DIR", joinpath(pwd(), "QuIP"))

    if !isdir(repo_dir)
        run(`git clone --depth 1 https://github.com/SECQUOIA/QuIP.git $repo_dir`)
    end

    repo_dir
else
    detect_quip_repo_dir()
end

JULIA_NOTEBOOKS_DIR = joinpath(QUIP_REPO_DIR, "notebooks_jl")

Pkg.activate(JULIA_NOTEBOOKS_DIR)
Pkg.instantiate(; io = devnull)

cd(JULIA_NOTEBOOKS_DIR)

About this Notebook

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 JuMP, 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 Graphs.jl for network models/graphs.

Problem statement

We define a QUBO as the following optimization problem:

min⁡x∈{0,1}n∑(ij)∈E(G)Qijxixj+∑i∈V(G)Qiixi+β=min⁡x∈{0,1}nx′Qx+β\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}' \mathbf{Q} \mathbf{x} + \beta
(1)

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

Example

Suppose we want to solve the following problem via QUBO

min⁡x2x0+4x1+4x2+4x3+4x4+4x5+5x6+4x7+5x8+6x9+5x10s.t.[100111011110101011011100101011111]x=[111] x∈{0,1}11\begin{array}{rl} \displaystyle% \min_{\mathbf{x}} & 2x_0+4x_1+4x_2+4x_3+4x_4+4x_5+5x_6+4x_7+5x_8+6x_9+5x_{10} \\ \textrm{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{array}
(2)

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 = [
    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 = [1, 1, 1]
c = [2, 4, 4, 4, 4, 4, 5, 4, 5, 6, 5];

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

min⁡xc′xs.t.Ax=b x∈{0,1}11\begin{array}{rl} \displaystyle% \min_{\mathbf{x}} &\mathbf{c}' \mathbf{x} \\ \textrm{s.t.} & \mathbf{A}\mathbf{x} = \mathbf{b} \\ ~ & \mathbf{x} \in \{0,1 \}^{11} \end{array}
(3)

as follows:

min⁡xc′x+ρ(Ax−b)′(Ax−b)s.t.x∈{0,1}11\begin{array}{rl} \displaystyle% \min_{\mathbf{x}} & \mathbf{c}' \mathbf{x} + \rho (\mathbf{A}\mathbf{x}-\mathbf{b})' (\mathbf{A}\mathbf{x}-\mathbf{b}) \\ \textrm{s.t.} & \mathbf{x} \in \{0,1 \}^{11} \end{array}
(4)

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

ρ(Ax−b)′(Ax−b)=ρ(x′(A′A)x−2(A′b)x+b′b)\rho(\mathbf{A}\mathbf{x}-\mathbf{b})'(\mathbf{A}\mathbf{x}-\mathbf{b}) = \rho( \mathbf{x}'(\mathbf{A}'\mathbf{A}) \mathbf{x} - 2(\mathbf{A}'\mathbf{b}) \mathbf{x} + \mathbf{b}'\mathbf{b} )
(5)

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

using LinearAlgebra
ϵ = 1
ρ = sum(abs, c) + ϵ
Q = diagm(c) + ρ * (A'A - 2 * diagm(A'b))
β = ρ * b'b

display(Q)
println(β)
11×11 Matrix{Int64}: -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.

using Plots
plot(QUBOTools.SystemLayoutPlot(Q))
Loading...

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

using JuMP
using QUBO
using DWave
    CondaPkg Found dependencies: /home/pedroxavier/gits/DWave.jl/CondaPkg.toml


    CondaPkg Found dependencies: /home/pedroxavier/.julia/packages/PythonCall/1f5yE/CondaPkg.toml

    CondaPkg Found dependencies: /home/pedroxavier/.julia/packages/DWaveNeal/3cu1x/CondaPkg.toml

    CondaPkg Found dependencies: /home/pedroxavier/.julia/packages/PythonPlot/f591M/CondaPkg.toml

    CondaPkg Found dependencies: /home/pedroxavier/gits/DWave.jl/CondaPkg.toml

    CondaPkg Found dependencies: /home/pedroxavier/.julia/packages/PythonCall/1f5yE/CondaPkg.toml

    CondaPkg Dependencies already up to date

# Define empty model
qubo_model = Model()

# Define the variables
@variable(qubo_model, x[1:11], Bin)

# Define the objective function
@objective(qubo_model, Min, x' * Q * x + β)

# Print the model
print(qubo_model)
Min -46 x[1]² + 96 x[4]*x[1] + 96 x[5]*x[1] + 96 x[6]*x[1] + 96 x[8]*x[1] + 96 x[9]*x[1] + 96 x[10]*x[1] + 96 x[11]*x[1] - 44 x[2]² + 96 x[4]*x[2] + 96 x[6]*x[2] + 96 x[7]*x[2] + 96 x[9]*x[2] + 96 x[10]*x[2] + 96 x[11]*x[2] - 44 x[3]² + 96 x[5]*x[3] + 96 x[7]*x[3] + 96 x[8]*x[3] + 96 x[9]*x[3] + 96 x[10]*x[3] + 96 x[11]*x[3] - 92 x[4]² + 96 x[5]*x[4] + 192 x[6]*x[4] + 96 x[7]*x[4] + 96 x[8]*x[4] + 192 x[9]*x[4] + 192 x[10]*x[4] + 192 x[11]*x[4] - 92 x[5]² + 96 x[6]*x[5] + 96 x[7]*x[5] + 192 x[8]*x[5] + 192 x[9]*x[5] + 192 x[10]*x[5] + 192 x[11]*x[5] - 92 x[6]² + 96 x[7]*x[6] + 96 x[8]*x[6] + 192 x[9]*x[6] + 192 x[10]*x[6] + 192 x[11]*x[6] - 91 x[7]² + 96 x[8]*x[7] + 192 x[9]*x[7] + 192 x[10]*x[7] + 192 x[11]*x[7] - 92 x[8]² + 192 x[9]*x[8] + 192 x[10]*x[8] + 192 x[11]*x[8] - 139 x[9]² + 288 x[10]*x[9] + 288 x[11]*x[9] - 138 x[10]² + 288 x[11]*x[10] - 139 x[11]² + 144
Subject to
 
x[1] binary
 x[2] binary
 x[3] binary
 x[4] binary
 x[5] binary
 x[6] binary
 x[7] binary
 x[8] binary
 x[9] binary
 x[10] binary
 x[11] binary

# Use D-Wave's simulated annealer 'Neal'
set_optimizer(qubo_model, DWave.Neal.Optimizer)

set_optimizer_attribute(qubo_model, "num_reads", 1_000)

optimize!(qubo_model)

println("Minimum energy: $(objective_value(qubo_model))")
Minimum energy: 5.0

plot(QUBOTools.EnergyFrequencyPlot(qubo_model))
Loading...

Notice that this is the same example we have been solving earlier (via Integer Programming in the Quiz 2, via Ising model and QUBO in Notebook 2).

Now let’s solve this using Quantum Annealing!

First, we start by defining the DWAVE_API_TOKEN environment variable.

ENV["DWAVE_API_TOKEN"] = "<YOUR_KEY_HERE>";
# Use D-Wave's quantum annealer
set_optimizer(qubo_model, DWave.Optimizer)

set_optimizer_attribute(qubo_model, "num_reads", 1024)

optimize!(qubo_model)

println("Minimum energy: $(objective_value(qubo_model))")
Minimum energy: 5.0

plot(QUBOTools.EnergyFrequencyPlot(qubo_model))
Loading...
using Graphs
# Graph corresponding to D-Wave 2000Q
sampler = DWave.dwave_system.DWaveSampler(token = ENV["DWAVE_API_TOKEN"])
function get_topology(sampler)
    if string(sampler.solver.id) == "DW_2000Q_6"
        return DWave.dwave_networkx.chimera_graph(
            16,
            node_list=sampler.nodelist,
            edge_list=sampler.edgelist,
        )
    elseif string(sampler.solver.id) == "Advantage_system1.1"
        return DWave.dwave_networkx.pegasus_graph(
            16,
            node_list=sampler.nodelist,
            edge_list=sampler.edgelist,
        )
    elseif string(sampler.solver.id) == "Advantage_system4.1"
        return DWave.dwave_networkx.pegasus_graph(
            16,
            node_list=sampler.nodelist,
            edge_list=sampler.edgelist,
        )
    else
        error("Unknown solver id '$(sampler.solver.id)'")

        return nothing
    end
end

function draw_topology(sampler)
    χ = get_topology(sampler)

    g = Graphs.grpah
end

draw_topology(sampler)
sol = QUBOTools.sampleset(unsafe_backend(qubo_model))
149-element QUBOTools.SampleSet{Float64, Int64}: [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1] [-1, -1, -1, -1, -1, -1, -1, -1, 1, -1, -1] [-1, -1, -1, -1, -1, -1, -1, -1, -1, 1, -1] [1, -1, -1, -1, -1, -1, 1, -1, -1, -1, -1] [-1, -1, 1, -1, -1, 1, -1, -1, -1, -1, -1] [-1, -1, 1, 1, -1, -1, -1, -1, -1, -1, -1] [-1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1] [-1, 1, -1, -1, 1, -1, -1, -1, -1, -1, -1] [1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1] [-1, -1, -1, -1, -1, -1, -1, 1, -1, -1, -1] ⋮ [-1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1] [-1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1] [-1, -1, -1, 1, 1, -1, -1, -1, 1, -1, -1] [1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1] [-1, -1, -1, 1, -1, -1, 1, -1, 1, -1, -1] [-1, 1, 1, -1, -1, -1, -1, 1, -1, -1, 1] [-1, 1, 1, -1, 1, -1, -1, -1, -1, -1, 1] [-1, 1, 1, 1, -1, -1, -1, -1, 1, -1, -1] [1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1]
data = QUBOTools.metadata(sol)
info = get(data, "dwave_info", nothing)
Python: {'timing': {'qpu_sampling_time': 88043.52, 'qpu_anneal_time_per_sample': 20.0, 'qpu_readout_time_per_sample': 45.44, 'qpu_access_time': 103802.29, 'qpu_access_overhead_time': 4434.71, 'qpu_programming_time': 15758.77, 'qpu_delay_time_per_sample': 20.54, 'post_processing_overhead_time': 2044.0, 'total_post_processing_time': 2044.0}, 'problem_id': 'cb0dd99d-cc2d-487d-9161-75ab2c607cca'}

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

References

D-Wave
Quantum Annealing via D-Wave
Benchmarking
About this Notebook