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.

Quantum Annealing via D-Wave (Julia)

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

Open In Colab SECQUOIA

Environment setup

For local execution from the repository root, run uv sync --group qubo and make setup-julia NOTEBOOK=notebooks_jl/4-DWave.ipynb before launching Jupyter. This notebook reuses the repo-local Python environment for the D-Wave Ocean stack instead of relying on Julia’s CondaPkg resolver.

This notebook should open directly with the Julia runtime in Google Colab. If it does not, switch the runtime before running the next cell. The setup cell will clone SECQUOIA/QuIP into the Colab workspace when needed, activate notebooks_jl/envs/4-DWave/Project.toml, install the Python D-Wave Ocean packages into the Colab Python runtime, and print progress for clone, activate, instantiate, and package-load steps. Package precompilation still happens as needed during setup. The first Colab run can still take several minutes.

function load_notebook_bootstrap()
    candidates = (
        joinpath(pwd(), "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "..", "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "QuIP", "scripts", "notebook_bootstrap.jl"),
        joinpath("/content", "QuIP", "scripts", "notebook_bootstrap.jl"),
    )

    for candidate in candidates
        if isfile(candidate)
            include(candidate)
            return nothing
        end
    end

    in_colab = haskey(ENV, "COLAB_RELEASE_TAG") || haskey(ENV, "COLAB_JUPYTER_IP") || isdir(joinpath("/content", "sample_data"))
    if in_colab
        repo_dir = get(ENV, "QUIP_REPO_DIR", joinpath(pwd(), "QuIP"))
        if !isdir(repo_dir)
            println("[bootstrap] Cloning SECQUOIA/QuIP into $repo_dir")
            run(`git clone --depth 1 https://github.com/SECQUOIA/QuIP.git $repo_dir`)
        end
        include(joinpath(repo_dir, "scripts", "notebook_bootstrap.jl"))
        return nothing
    end

    error("Could not locate scripts/notebook_bootstrap.jl from $(pwd()).")
end

load_notebook_bootstrap()

BOOTSTRAP = QuIPNotebookBootstrap.bootstrap_notebook("4-DWave")

QUIP_REPO_DIR = BOOTSTRAP.repo_dir
JULIA_NOTEBOOKS_DIR = BOOTSTRAP.notebooks_dir
JULIA_PROJECT_DIR = BOOTSTRAP.project_dir
IN_COLAB = BOOTSTRAP.in_colab

Quantum Annealing via D-Wave (Julia)

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.

if !@isdefined(LinearAlgebra)
    using LinearAlgebra
end

ϵ = 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.

if !@isdefined(Plots)
    using Plots
end

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

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

if !@isdefined(JuMP)
    using JuMP
end

if !@isdefined(QUBO)
    using QUBO
end

if !@isdefined(DWave)
    using DWave
end

# 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...
if !@isdefined(Graphs)
    using Graphs
end

# 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.solution(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 (Python)
Benchmarking
Benchmarking (Python)