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.

Quadratic Unconstrained Binary Optimization

QUBO & Ising Models

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

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

Open In Colab SECQUOIA

Quadratic Unconstrained Binary Optimization

This notebook will explain the basics of the QUBO modeling. In order to implement the different QUBOs we will use D-Wave’s packages dimod, and then solve them using neal’s implementation of simulated annealing. We will also leverage the use of D-Wave’s package dwavebinarycsp to translate constraint satisfaction problems to QUBOs. Finally, for Groebner basis computations we will use Sympy for symbolic computation in Python and Networkx 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+cQ=min⁡x∈{0,1}nx⊤Qx+cQ\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
(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 an adjacency matrix QQ. We also include an arbitrary offset cQc_Q.

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}=

[111]\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

# Let's install dimod, neal, and pyomo
if IN_COLAB:
    !pip install -q pyomo
    !pip install dimod
    !pip install dwave-neal
# Import the Pyomo library, which can be installed via pip, conda or from Github https://github.com/Pyomo/pyomo
import pyomo.environ as pyo
# 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 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\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.

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 using DWaves code for complete enumeration and simulated annealing (eventually with Quantum annealiing too!).

model = dimod.BinaryQuadraticModel.from_qubo(Q, offset=Beta)
print(model)

Since the problem is relatively small (11 variables, 211=20482^{11}=2048 combinations), we can afford to enumerate all the solutions.

exactSampler = dimod.reference.samplers.ExactSolver()
exactSamples = exactSampler.sample(model)
# Some useful functions to get plots
def plot_enumerate(results, title=None):

    plt.figure()

    energies = [datum.energy for datum in results.data(
        ['energy'], sorted_by='energy')]
    
    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_samples(results, title=None):
    plt.figure()
    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')

    counts = Counter(samples)
    total = len(samples)
    for key in counts:
        counts[key] /= total
    df = pd.DataFrame.from_dict(counts, orient='index').sort_index()
    df.plot(kind='bar', legend=None)

    plt.xticks(rotation=80)
    plt.ylabel('Probabilities')
    plt.title(str(title))
    plt.show()
    print("minimum energy:", min(energies))


def plot_energies(results, title=None, skip=1):
    # skip parameter given to avoid putting all xlabels
    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()
    ax = df.plot(kind='bar', legend=None)

    plt.xlabel('Energy')
    plt.ylabel('Probabilities')
    # Plot only a subset of xlabels (every skip steps)
    ax.set_xticklabels([t if not i%skip else "" for i,t in enumerate(ax.get_xticklabels())])
    plt.title(str(title))
    plt.show()
    print("minimum energy:", min(energies))
plot_enumerate(exactSamples, title='Enumerate all solutions')
plot_energies(exactSamples, title='Enumerate all solutions', skip=10)

Let’s now solve this QUBO via traditional Integer Programming.

# We do not need to worry about the tranformation to QUBO since dimod takes care of it
Q, c = model.to_qubo()

# Define the model
model_pyo = pyo.ConcreteModel(name='QUBO example as an IP, 47-779/785 QuIPML')

I = range(len(model))
J = range(len(model))
#Define the original variables
model_pyo.x = pyo.Var(I, domain=pyo.Binary)
# Define the edges variables
model_pyo.y = pyo.Var(I, J, domain=pyo.Binary)

obj_expr = c

# add model constraints
model_pyo.c1 = pyo.ConstraintList()
model_pyo.c2 = pyo.ConstraintList()
model_pyo.c3 = pyo.ConstraintList()
for (i,j) in Q.keys():
    if i != j:
        model_pyo.c1.add(model_pyo.y[i,j] >= model_pyo.x[i] + model_pyo.x[j] - 1)
        model_pyo.c2.add(model_pyo.y[i,j] <= model_pyo.x[i])
        model_pyo.c3.add(model_pyo.y[i,j] <= model_pyo.x[j])
        obj_expr += Q[i,j]*model_pyo.y[i,j]
    else:
        obj_expr += Q[i,j]*model_pyo.x[i]

# Define the objective function
model_pyo.objective = pyo.Objective(expr = obj_expr, sense=pyo.minimize)
# Print the model
model_pyo.display()

Let’s install the MIP solver GLPK

# Let's install the LP/MIP solver GLPK
if IN_COLAB:
    !apt-get install -y -qq glpk-utils
# Define the solver GLPK
if IN_COLAB:
    opt_glpk = pyo.SolverFactory('glpk', executable='/usr/bin/glpsol')
else:
    opt_glpk = pyo.SolverFactory('glpk')
# Here we could use another solver, e.g. gurobi or cplex
# opt_gurobi = pyo.SolverFactory('gurobi')
# We obtain the solution with GLPK
result_obj = opt_glpk.solve(model_pyo, tee=False)
model_pyo.display()

We observe that the optimal solution of this problem is x8=1,0x_{8} = 1, 0 otherwise, leading to an objective of 5. Notice that this problem has a degenerate optimal solution given that x10=1,0x_{10} = 1, 0 otherwise also leads to the same solution.

Ising model

This notebook will explain the basics of the Ising model. In order to implement the different Ising Models we will use D-Wave’s packages dimod and neal, for defining the Ising model and solving it with simulated annealing, respectively. When posing the problems as Integer programs, we will model using Pyomo, an open-source Python package, which provides a flexible access to different solvers and a general modeling framework for linear and nonlinear integer programs. The examples solved here will make use of open-source solver GLPK for mixed-integer linear programming.

Problem statement

We pose the Ising problem as the following optimization problem:

min⁡σ∈{−1,+1}nH(σ)=min⁡σ∈{−1,+1}n∑(ij)∈E(G)Jijσiσj+∑i∈V(G)hiσi+cI\min_{\sigma \in \{ -1,+1 \}^n} H(\sigma) =\min_{\sigma \in \{ -1,+1 \}^n} \sum_{(ij) \in E(G)} J_{ij}\sigma_i\sigma_j + \sum_{i \in V(G)}h_i\sigma_i + c_I
(6)

where we optimize over spins σ∈{−1,+1}n\sigma \in \{ -1,+1 \}^n, on a constrained graph G(V,E)G(V,E), where the quadratic coefficients are JijJ_{ij} and the linear coefficients are hih_i. We also include an arbitrary offset of the Ising model cIc_I.

Example

Suppose we have an Ising model defined from

h=[145.0122.0122.0266.0266.0266.0242.5266.0386.5387.0386.5],J=[00024242424242424240002402424242424240000240242424242400002448242448484800000242448484848000000242448484800000002448484800000000484848000000000727200000000007200000000000] and β=1319.5h = \begin{bmatrix} 145.0 \\ 122.0 \\ 122.0 \\ 266.0 \\ 266.0 \\ 266.0 \\ 242.5 \\ 266.0 \\ 386.5 \\ 387.0 \\ 386.5 \end{bmatrix}, J = \begin{bmatrix} 0 & 0 & 0 & 24 & 24 & 24 & 24 & 24 & 24 & 24 & 24\\ 0 & 0 & 0 & 24 & 0 & 24 & 24 & 24 & 24 & 24 & 24\\ 0 & 0 & 0 & 0 & 24 & 0 & 24 & 24 & 24 & 24 & 24\\ 0 & 0 & 0 & 0 & 24 & 48 & 24 & 24 & 48 & 48 & 48\\ 0 & 0 & 0 & 0 & 0 & 24 & 24 & 48 & 48 & 48 & 48\\ 0 & 0 & 0 & 0 & 0 & 0 & 24 & 24 & 48 & 48 & 48\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 24 & 48 & 48 & 48\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 48 & 48 & 48\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 72 & 72\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 72\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ \end{bmatrix} \text{ and } \beta = 1319.5
(7)

Let’s solve this problem

# These could also be simple lists and numpy matrices
h = {0: 145.0, 1: 122.0, 2: 122.0, 3: 266.0, 4: 266.0, 5: 266.0, 6: 242.5, 7: 266.0, 8: 386.5, 9: 387.0, 10: 386.5}
J = {(0, 3): 24.0, (0, 4): 24.0, (0, 5): 24.0, (0, 7): 24.0, (0, 8): 24.0, (0, 9): 24.0, (0, 10): 24.0, (1, 3): 24.0, (1, 5): 24.0, (1, 6): 24.0, (1, 8): 24.0, (1, 9): 24.0, (1, 10): 24.0, (2, 4): 24.0, (2, 6): 24.0, (2, 7): 24.0, (2, 8): 24.0, (2, 9): 24.0, (2, 10): 24.0, (3, 4): 24.0, (3, 5): 48.0, (3, 6): 24.0, (3, 7): 24.0, (3, 8): 48.0, (3, 9): 48.0, (3, 10): 48.0, (4, 5): 24.0, (4, 6): 24.0, (4, 7): 48.0, (4, 8): 48.0, (4, 9): 48.0, (4, 10): 48.0, (5, 6): 24.0, (5, 7): 24.0, (5, 8): 48.0, (5, 9): 48.0, (5, 10): 48.0, (6, 7): 24.0, (6, 8): 48.0, (6, 9): 48.0, (6, 10): 48.0, (7, 8): 48.0, (7, 9): 48.0, (7, 10): 48.0, (8, 9): 72.0, (8, 10): 72.0, (9, 10): 72.0}
cI = 1319.5

model_ising = dimod.BinaryQuadraticModel.from_ising(h, J, offset=cI)

Since the problem is relatively small (11 variables, 211=20482^{11}=2048 combinations), we can afford to enumerate all the solutions.

exactSamples = exactSampler.sample(model_ising)
plot_enumerate(exactSamples, title='Enumerate all solutions')
plot_energies(exactSamples, title='Enumerate all solutions')

Let’s now solve this Ising Model via traditional Integer Programming.

# We do not need to worry about the tranformation from Ising to QUBO since dimod takes care of it
Q, c = model_ising.to_qubo()

# Define the model
model_ising_pyo = pyo.ConcreteModel(name='Ising example as an IP, 47-779/785 QuIPML')

I = range(len(h))
J = range(len(h))
#Define the original variables
model_ising_pyo.x = pyo.Var(I, domain=pyo.Binary)
# Define the edges variables
model_ising_pyo.y = pyo.Var(I, J, domain=pyo.Binary)

obj_expr = c

# add model constraints
model_ising_pyo.c1 = pyo.ConstraintList()
model_ising_pyo.c2 = pyo.ConstraintList()
model_ising_pyo.c3 = pyo.ConstraintList()
for (i,j) in Q.keys():
    if i != j:
        model_ising_pyo.c1.add(model_ising_pyo.y[i,j] >= model_ising_pyo.x[i] + model_ising_pyo.x[j] - 1)
        model_ising_pyo.c2.add(model_ising_pyo.y[i,j] <= model_ising_pyo.x[i])
        model_ising_pyo.c3.add(model_ising_pyo.y[i,j] <= model_ising_pyo.x[j])
        obj_expr += Q[i,j]*model_ising_pyo.y[i,j]
    else:
        obj_expr += Q[i,j]*model_ising_pyo.x[i]

# Define the objective function
model_ising_pyo.objective = pyo.Objective(expr = obj_expr, sense=pyo.minimize)
# Print the model
model_ising_pyo.display()
# We obtain the solution with GLPK
result_obj = opt_glpk.solve(model_ising_pyo, tee=False)
model_ising_pyo.display()

We observe that the optimal solution of this problem is x10=1,0x_{10} = 1, 0 otherwise, leading to an objective of 5. Notice that this problem has a degenerate optimal solution given that x8=1,0x_8 = 1, 0 otherwise also leads to the same solution.

Let’s go back to the slides

We can also solve this problem using Simulated Annealing

simAnnSampler = neal.SimulatedAnnealingSampler()
simAnnSamples = simAnnSampler.sample(model_ising, num_reads=1000)
plot_enumerate(simAnnSamples, title='Simulated annealing in default parameters')
plot_energies(simAnnSamples, title='Simulated annealing in default parameters')
simAnnSamples.info

Let’s go back to the slides

Let’s solve the graph coloring problem in the slides using QUBO.

Vertex kk-coloring of graphs

Given a graph G(V,E)G(V, E), where VV is the set of vertices and EE is the set of edges of GG, and a positive integer kk, we ask if it is possible to assign a color to every vertex from VV, such that adjacent vertices have different colors assigned.

G(V,E)G(V, E) has 12 vertices and 23 edges. We ask if the graph is 3–colorable. Let’s first encode VV and EE using Julia’s built–in data structures:

Note: This tutorial is heavily inspired in D-Wave’s Map coloring of Canada found here.

# Let's install with dimod and neal
if IN_COLAB:
    !pip install dwavebinarycsp
    !pip install dwavebinarycsp[maxgap]
    !pip install dwavebinarycsp[mip]

import dwavebinarycsp
V = range(1, 12+1)
E = [(1,2),(2,3),(1,4),(1,6),(1,12),(2,5),(2,7),(3,8),(3,10),(4,11),(4,9),(5,6),(6,7),(7,8),(8,9),(9,10),(10,11),(11,12),(5,12),(5,9),(6,10),(7,11),(8,12)]
layout = {i: [np.cos((2*i+1)*np.pi/8),np.sin((2*i+1)*np.pi/8)] for i in np.arange(5,13)}
layout[1] = [-1.5,1.5]
layout[2] = [1.5,1.5]
layout[3] = [1.5,-1.5]
layout[4] = [-1.5,-1.5]
G = nx.Graph()
G.add_edges_from(E)
nx.draw(G, with_labels=True, pos=layout)
# Function for the constraint that two nodes with a shared edge not both select
# one color
def not_both_1(v, u):
    return not (v and u)

# Valid configurations for the constraint that each node select a single color, in this case we want to use 3 colors
one_color_configurations = {(0, 0, 1), (0, 1, 0), (1, 0, 0)}
colors = len(one_color_configurations)

# Create a binary constraint satisfaction problem
csp = dwavebinarycsp.ConstraintSatisfactionProblem(dwavebinarycsp.BINARY)

# Add constraint that each node select a single color
for node in V:
    variables = ['x'+str(node)+','+str(i) for i in range(colors)]
    csp.add_constraint(one_color_configurations, variables)

# Add constraint that each pair of nodes with a shared edge not both select one color
for edge in E:
    v, u = edge
    for i in range(colors):
        variables = ['x'+str(v)+','+str(i), 'x'+str(u)+','+str(i)]
        csp.add_constraint(not_both_1, variables)

Defining the Binary Quandratic model (QUBO) using the CSP library we have:

bqm = dwavebinarycsp.stitch(csp)
simAnnSamples = simAnnSampler.sample(bqm, num_reads=1000)
plot_enumerate(simAnnSamples, title='Simulated annealing in default parameters')
plot_energies(simAnnSamples, title='Simulated annealing in default parameters')

Because of precision issues in the translation to BQM, we may obtain very tiny coefficeints that should be zero. In any case, since this is a constraint satisfaction problem, any of the solutions with energy ~0 is a valid coloring.

# Check that a good solution was found
sample = simAnnSamples.first.sample     # doctest: +SKIP
if not csp.check(sample):           # doctest: +SKIP
        print("Failed to color map. Try sampling again.")
else:
        print(sample)
# Function that plots a returned sample
def plot_map(sample):
    # Translate from binary to integer color representation
    color_map = {}
    for node in V:
          for i in range(colors):
            if sample['x'+str(node)+','+str(i)]:
                color_map[node] = i
    # Plot the sample with color-coded nodes
    node_colors = [color_map.get(node) for node in G.nodes()]
    nx.draw(G, with_labels=True, pos=layout, node_color=node_colors)
    plt.show()
plot_map(sample)
Mathematical Programming
Environment Setup
QUBO
Environment Setup