Quantum Integer Programming (QuIP) 47-779: Lecture Notes

David E. Bernal, Sridhar Tayur, Davide Venturelli

Fall 2020

Abstract

This lecture series on Quantum Integer Programming (QuIP) – created by Professor Sridhar Tayur, David E.
Bernal and Dr. Davide Venturelli, a collaboration between CMU and USRA, with the support from Amazon
arXiv:2012.11382v2 [quant-ph] 11 Jan 2021

Braket during Fall 2020 – is intended for students and researchers interested in Integer Programming and
the potential of near term quantum and quantum inspired computing in solving optimization problems.
Originally created for Tepper School of Business course 47-779 (at CMU), these were also used for the
course ID5840 (at IIT-Madras, by Professors Anil Prabhakar and Prabha Mandayam) whose students (listed
at the beginning of each lecture) were scribes. Dr. Vikesh Siddhu, post-doc in CMU Quantum Computing
Group, assisted during the lectures, student projects and with proof-reading this scribe.
Through these lectures one will learn to formulate a problem and map it to a Quadratic Unconstrained
Binary Optimization (QUBO) problem, understand various mapping and techniques like the Ising model,
Graver Augmented Multiseed Algorithm (GAMA), Simulated or Quantum Annealing and QAOA, and ideas
on how to solve these Integer problems using these quantum and classical methods.
The course website (with lecture videos and colab notebooks): [Link]
Keywords: Ising model, Integer Programming, Computational Algebraic Geometry, Graver Basis, Quan-
tum Annealing, Simulated Annealing, Combinatorial Optimization, Graph coloring, discrete nonlinear opti-
mization.

0-1
 47-779. Quantum Integer Programming
Mini-1, Fall 2020

Room: Zoom Online Time: Tuesday and Thursday 5:20pm-7:10pm

Instructors:

Sridhar Tayur Email: stayur@[Link] Office: 4216 Tepper Quad

Davide Venturelli Email: DVenturelli@[Link] Office: Online

David E. Bernal Email: bernalde@[Link] Office: 3116 Doherty Hall

Course Webpage:
• [Link]

Objectives: This course is primarily designed for graduate students (and advanced undergraduates) in-
terested in integer programming (with non-linear objective functions) and the potential of near-term quan-
tum and quantum-inspired computing for solving combinatorial optimization problems. By the end of the
semester, someone enrolled in this course should be able to:
• Appreciate the current status of quantum computing and its potential use for integer programming
• Access and use quantum computing resources (such as D-Wave Quantum Annealers)
• Set up a given integer program to be solved with quantum computing
• Work in groups collaboratively on a state-of-the-art project involving applications of quantum com-
puting and integer programming
This course is not going to focus on the following topics:
• Quantum Gates and Circuits
• Computational complexity theory
• Quantum Information Theory
• Analysis of speedup using differential geometry, algebraic topology, etc.

Prerequisite classes and capabilities: Although this class has no explicit prerequisites we consider a
list of recommended topics and skills that the student should feel comfortable with. An undergraduate-level
understanding of probability, calculus, statistics, graph theory, algorithms, and linear algebra is assumed.
Knowledge of linear and integer programming will be useful for this course. Programming skills are strongly
recommended. Basic concepts in physics are recommended but lack of prior knowledge is not an issue as
pertinent ones will be covered in the lectures. No particular knowledge in quantum mechanics or algebraic
geometry is required.

Students with backgrounds in operations research, industrial engineering, chemical engineering, electrical
engineering, physics, computer science, or applied mathematics are strongly encouraged to consider taking
this course.

0-2
Tentative Course Outline:

Part 1 - Integer programming (classical methods) . . . . . . . . . . . . . . . . . 1 week

• Integer Programming basics [1].
• Cutting plane theory and relaxations [1].
• Introduction to Test Sets [2, 3, 4, 5].
– Gröbner basis [6, 7, 8].
– Graver basis [9].
Part 2 - Ising, QUBO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 week
• Ising model basics [10, 11, 12].
• Simulated Annealing [13, 14].
• Markov-chain Monte Carlo methods [15, 16, 17, 18].
• Benchmarking classical methods [19, 20].
• Formulating combinatorial problems as QUBOs [21].
Part 3 - GAMA: Graver Augmented Multiseed algorithm . . . . . . . 1 week
• GAMA [22].
– Applications: Portfolio Optimization [22], Cancer Genomics [23]
– Quantum Inspired: Quadratic (Semi-)Assignment Problem [24].
Part 4 - Quantum methods for solving Ising/QUBO . . . . . . . . . . . . . . 1 week
• AQC, Quantum Annealing and D-Wave [25, 26, 27, 28, 29, 30, 31].
• QAOA: Quantum Approximate Optimization Algorithm [32, 33, 34].
Part 5 - Hardware for solving Ising/QUBO . . . . . . . . . . . . . . . . . . . . . . . . 1 week
• Graphical Processing Units [35, 36, 37].
• Tensor Processing Units [38].
• Complementary metal-oxide-semiconductors (CMOS) [39].
• Digital Annealers [40].
• Oscillator Based Computing [41, 42].
• Coherent Ising Machines [43, 44, 45, 46, 47, 48].
Part 6 - Other topics and project presentations . . . . . . . . . . . . . . . . . . . 1 week
• Compiling
– Quantum Annealing [49, 50].
– Gate-based Noisy Intermediate Scale Quantum (NISQ) devices [51].
• Adiabatic Quantum Computing and Algebraic Geometry
– Minimizing Polynomial Functions [52].
– Prime Factorization [53].

Project description: This project will be completed in groups of 2-4 students and will reflect the un-
derstanding of the students of the material covered in the lecture. The components of this project are the
following
• Identify a problem that can be posed as an Integer Program. Discuss the importance of this problem.
• Solve instances of the identifies problem using classical tools. Identify which are the sources of com-
plexity while solving this problem.
• Model the problem as a Quadratic Unconstrained Binary Optimization (QUBO). Verify that the re-
formulation of the problem is valid, in the sense that it represents the original problem.
• Solve the resulting QUBO using non-conventional methods, e.g. Quantum Annealing, QAOA, simu-
lated annealing in GPUs/TPUs, etc. Compare at least two different methods.
• Write a report outlining the different approaches used and highlighting the knowledge obtained while
developing the project.
• Hold a final presentation in front of the class reporting the findings of the project.

0-3
Highlights: The specific skills that students will gain that will make them “quantum ready” for industry
or further academic research in this course are:

1. Classical

(a) Given a practical problem (from supply chain or physics or anything else), formulate it as a non-
linear integer program. We will provide a few practical problems, but we encourage you to suggest
one that you are already working on or are interested in.
(b) Solve such formulations via classical solvers.
2. Quantum

(a) Reformulate the problem to be “quantum ready” by making it in the form of a QUBO.
(b) Solve the QUBO “brute force” through D-Wave or IBM (via QAOA).
3. Hybrid Quantum-Classical
(a) Reformulate the problem again in the form suitable for GAMA.
(b) Solve GAMA compatible formulation via D-Wave and/or via QAOA.

USRA Collaboration

1. Access to D-Wave systems might be available via written proposals to the University Space Research As-
sociation (USRA). See [Link] for terms and conditions. The course
will discuss proposal preparation.

2. Students of this course are encouraged to apply to the Feynman Academy Internship program https://
[Link]/quantum/qacademy that sponsors research projects at NASA Ames Research Center.

Amazon Braket Collaboration

1. Access to D-Wave, Rigetti and ionQ have been made available by Amazon Braket.

2. Students will be provided individual sub-accounts for use with pre-funded amounts to access the quan-
tum machines.

Casual References: There is no single text book for the course. This is a short list of various interesting
and useful books that will be mentioned during the course. You need to consult them occasionally.
• Georges Irfah, The Universal History of Computing, John Wiley & Sons, 2001.
• A. Das and B.K. Chakrabarti (Eds.). Quantum Annealing and Related Optimization Methods, Springer-
Verlag, 2005.
• Eleanor G. Rieffel and Wolfgang H. Polak, Quantum Computing: A Gentle Introduction, MIT Press,
2011.
• Richard J. Lipton and Kenneth W. Regan, Quantum Algorithms via Linear Algebra. A Primer , MIT
Press, 2014.

0-4
 Quantum Integer Programming September-December 2020

Lecture 1: Integer Programming

Lecturer: Prof. Sridhar Tayur, David E. Bernal Scribe: Vighnesh Natarajan
Slides, Notebook Video, password: Pˆ2GSWny

1.1 Introduction

The standard approach to solving real world optimizations is as follows:

1. Mathematically model the problem

• Quantify the real world objective as a mathematical function
• Locate variables on which this function depends
• List real world constraints as constraints on the variables
2. Solve the model, i.e., find the best objective value and variables
3. Interpret Results

Using the first step above, most real world problems are brought to a mathematical model of the form:
minimize f (x)
x
subject to: g(x) ≤ 0
h(x) = 0
In this model, we call f (x) the objective function, g(x) ≤ 0 an inequality constraint, and h(x) = 0 an
equality constraint. Such models are called often called programs. These programs are categorized based
on the objective function and constraints, and the nature of the variables. Some popular categories include
Linear Programs, Convex Programs, and Integer Programs. We shall discuss such programs in more depth,
especially Integer Programs.

1.2 Linear Programming

Suppose there is a company that produces two different products, A and B, which can be sold at different
values, $5.5 and $2.1 per unit, respectively. The company only counts with a single machine with electricity
usage of at most 17kW/day, and producing each A and B consumes 8kW/day and 2kW/day, respectively.
Besides, the company can only produce at most 2 more units of A than B per day. This real world problem
can be mathematically modelled as follows.
Assuming the units produced of A are x1 and of B are x2 we have
max 5.5x1 + 2.1x2
x1 ,x2

s.t.x2 ≤ x1 + 2
8x1 + 2x2 ≤ 17
x1 , x2 ≥ 0

1-1
1-2 Lecture 1: Integer Programming

Figure 1.1: Feasible region and solution of sample LP

The model mentioned above is an example of a linear program (LP). In an LP, the variables are real,
continuous, and satisfy affine constraints. The objective function depends linearly on these variables. As a
result, the general form of an LP is

min cT x
x
subject to : Ax ≤ b

In an LP, the feasible region, i.e., the set variables which satisfy the problem’s constraints, is a convex
polyhedron. Such polyhedron are an intersection of half-spaces and hyperplanes.
The standard ways to solve an LP include
• Simplex methods
• Interior point methods

1.2.1 The simplex method

This is a vertex hopping method. The method provides a worst-case exponential time complexity. A
k-simplex is, in essence, a convex hull of k + 1 points. The optimal solution occurs at the vertices of the
simplex, where some of the constraints are active. Thus we find the optima by vertex hopping. The standard
LP can be expressed as a matrix, which is then used to solve the problem

     
 1 −cT 0   1 −cTB −cTN 0   1 0 −cTN zb 
  −→   −→  
     
0 A b 0 I D b 0 I D b

Here, we refer to the basic variables, −cTN is the relative cost, zb is the objective value, the subscript B denotes
the basic variables values, while D denotes the non-basic variables. The simplex method is implemented in
Lecture 1: Integer Programming 1-3

a variety of solvers including Gurobi1 , IBM CPLEX2 , XPRESS3 , MindOpt4 , CBC5 , GLPK6 , among others.

1.2.2 Interior point methods

These methods start with a point inside the feasible region and iteratively taking steps to never leave
the feasible region. This contrasts with optimization algorithms like projected gradient descent, projected
steepest descent, primal-dual optimization, etc.; where a step can leave the feasible set, and later be projected
back onto the feasible set. We start from a feasible solution and get closer and closer to the edges by having
a monotonic barrier when close. Thus we never reach the edge, just get closer and closer to activating the
appropriate constraint. This method is also known as a barrier method and is very useful when problems
are dual degenerate.
Examples of interior-point solvers include Gurobi7 , IBM CPLEX8 , MOSEK9 , and COPT LP10 , among
others.

1.3 Duality

Every problem has a dual problem that can be expressed using a Lagrangian. The original problem, that
has been described until now is known as the primal. For this primal problem:

minimize f (x)
x
subject to: g(x) ≤ 0
h(x) = 0,

define a Lagrangian,

L(x, λ, µ) = f (x) + λT g(x) + µT h(x)

λ ≥ 0.

We define

φ(λ, µ) = argmin L(x, λ, µ)

x

This function bounds from below all feasible values of the primal objective f . The dual problem is to
maximize φ(λ, µ). However, a question arises how the maximum dual value φ∗ is related to minimum value
f ∗ of the primal problem. There are two possibilities -
• Weak duality: φ∗ ≤ f ∗
• Strong duality: φ∗ = f ∗
1 [Link]
2 [Link]
3 [Link]
4 [Link]
5 [Link]
6 [Link]
7 [Link]
8 [Link]
9 [Link]
10 [Link]
1-4 Lecture 1: Integer Programming

This can be proven using the Minimax theorem. However in cases of an LP, the Minimax theorem can be
used to prove strong duality.
For convex optimization problems, where the objectives and inequality constraints in standard form have
to be convex, and the equality constraints have to be affine, Slater’s conditions are a sufficient condition for
strong duality to hold.
A short example has been illustrated below for an LP.

minimize cT x
x
subject to: Ax ≤ b
x ≥ 0, x  Rn

The Lagrange function becomes L(x, λ) = cT x + λT (b − Ax), x ≥ 0

φ(λ) = min cT x + λT (b − Ax)

x≥0

= λT b + min (cT − λT A)x

x≥0

φ(λ) = λT b, if cT − λT A ≥ 0
= −∞, otherwise

Thus, our dual problem becomes maximizing φ(λ), which is λT b. And we have strong duality holding in
LPs. So,

λTopt b = cT xopt

With this, we can move on to non-linear optimizations.

1.4 MIP and MINLP

Re-visiting our sample problem on the company producing products, What if we have a constraint that
lets us produce only integer number of products. This is something we need to modify and include in our
formulation.

max 5.5x1 + 2.1x2

x1 ,x2

s.t. x2 ≤ x1 + 2
8x1 + 2x2 ≤ 17
x1 , x2 ≥ 0
x1 , x2 are integers

The usual method to solve an Integer program is to first relax to a LP problem, and then round off to the
nearest integer. However, there are better techniques that provide more accurate solutions.
Before proceeding, we need to understand the challenges and the need for good Integer solvers. How hard
can it be to just enumerate all solutions? This depends on the complexity of the problem. Integer programs
generally have exponential complexity. This complexity issue is clearly understood from the image below.
A Mixed-Integer Program (MIP) is an optimization problem of the form
Lecture 1: Integer Programming 1-5

Figure 1.2: Feasible region and solution of sample ILP

Figure 1.3: Enumeration of feasible solutions

minimize f (x)
x
subject to : g(x) ≤ 0
some or all x0i s are integers

A MIP is a NP complete problem, and hence for large problem sizes, we do not have the means to solve it.
Various techniques to find approximate solutions, and if possible exact solutions are
• Branch and Bound - Solution of each search node using linear programming
• Cutting plane methods - based on polyhedral theory
In branch and bound, we essentially ask a set of “yes - no” questions. Answers to these questions help divide
the feasibly region, and also help prune it for optimal solutions. There are various ways we can prune a set,
by integrality, by bound and by infeasibility. The following image is an example of branch and bound on
our sample MIP. Cutting plane methods involve cutting off parts of the feasible region of the LP relaxed
problem. This can be used for convex optimization problems also, with a convex boundary approximated
1-6 Lecture 1: Integer Programming

Figure 1.4: Branching and bounding on sample MIP [1]

by a set of linear boundaries.

Recent improvements in MIP solvers are significant, with improvements in choosing cuts, pre-processing,
branching variable selection, and heuristics. We see that CPLEX 11 (2007) is 29000 times faster than
CPLEX 1.2 (1991). Gurobi 1 is comparable to CPLEX 11, and Gurobi 8 (2018) is 91 times faster. Overall
we see roughly a total speedup of 2.6M times from 1991 to 2018. In addition, we have had huge growth in
hardware speeds, which have also led to speedups. These solvers, which are now state of the art, are now
exploring using specialized hardware such as parallel processing GPU techniques to push performance.
The real world is nonlinear. Most systems are linearized about operating points, from where we can apply
linear programming and integer programming. However, there are cases where this linearization does not
yield satisfactory answers, and we need to involve nonlinear characteristics. Thus our objective function is
no more of the form cT x, and our objective is no more of the form Ax = b.
This brings us into the realm of Mixed Integer Non-Linear Programming or MINLP. MINLP are NP-hard
problems. We can see that all the problems we have considered are all subsets of MINLP, which is the general
and much more challenging set of problems. Similar to MIP, branching, and bounding is a method used to
solve MINLP problems. Also, MIP based methods are used, with the nonlinear objectives linearized in two
ways. The first way is where the objectives are under-estimated. The second way is where the objectives
are overestimated.
Continuing with our example of the two products, if we add an extra constraint, which is non linear as shown
below, we can see how the feasible region changes.

max 5.5x1 + 2.1x2

x1 ,x2

s.t. x2 ≤ x1 + 2
8x1 + 2x2 ≤ 17
(x2 − 1)2 ≤ 2 − x1
x1 , x2 ≥ 0
x1 , x2 are Integers

It is easy to observe that this particular non-linearity preserves convexity. Hence this is still a straightforward
problem to solve. Things start to get tricky when the problem becomes non-convex. We shall see non-
Lecture 1: Integer Programming 1-7

(a) Under-estimate the cost function (b) Over-estimate feasible region

Figure 1.5: Linearization of MINLP

Figure 1.6: Feasible region and optimal solution of the convex MINLP

convexity in the next example. Solvers can still solve such non-convex problems. However, we can understand
what can go wrong and where solvers can get stuck in these non-convex situations
We add an extra constraint to the formulation, as follows.
max 5.5x1 + 2.1x2
x1 ,x2

s.t. x2 ≤ x1 + 2
8x1 + 2x2 ≤ 17
(x2 − 1)2 ≤ 2 − x1
(x2 − 1)2 ≥ 0.5 + x1
x1 , x2 ≥ 0
x1 , x2 are Integers

A state of the art nonconvex MINLP solver is Baron11 . This solver is benchmarked using MINLPLib with
instances with approximately 2500 variables and constraints. MIP solvers like Gurobi, IBM-CPLEX, and
XPRESS are increasing their capabilities to solve nonlinear problems also.
11 [Link]
1-8 Lecture 1: Integer Programming

Figure 1.7: Feasible region and optimal solution of the non-convex MINLP

1.5 Complexity Analysis

Analysis of problems involves studying the time and memory required to solve them. The main focus here
is the time complexity involved in solving problems. Problems are usually categorized by the worst-case
scenario performance; however, some algorithms have become very popular due to their exceptional average-
case performance, such as quick-sort. A standard notation involved here is the “Big - O” notation.
The ”Big-O” notation is used where we have two functions f : S − → R+ and g : S − → R+ , there exist
positive numbers M and xo , such that f (x) ≤ M g(x) ∀x ≥ xo we say f (x) = O(g(x)).

We say that f (x) is polynomial time if it is polynomially bounded by g(x), in order of complexity. For
example, selection sort, is f (n) = O(n2 ), however for quick sort, f (n) = O(n log(n)). If there exists a poly-
nomial time solution for a problem, it belongs to the class P. LP problems are in the class P(Polynomial),
because interior point algorithms solve LPs in polynomial time.
Decision problems are those that have only two possible answers, ‘Yes’ or ‘No’. NP (Non-deterministic
Polynomial) class contains problems that cannot be solved in polynomial time complexity. However the
problem solutions can be verified in polynomial time. That is, a ‘Yes’ or ‘No’ answer can be verified in
polynomial time. A problem Q can be called NP-hard if all problems in NP can be reduced to Q in poly-
time. Integer Programming is NP-hard. We can transform an NP problem into an integer program. Hence,
if we find a way to solve IP in polynomial time, we can solve all NP problems in polynomial time and
prove the result P=NP! Unfortunately, integer programs with quadratic constraints have been proven to be
undecidable. Even after a long time without finding a solution, we cannot conclude that no solutions exists.
A problem is said to belong to the complexity class BPP (Bounded-error Probabilistic Polynomial time) if
there is an algorithm that solves the problem in such a way that the algorithm
• is allowed to make random decisions.
• is guaranteed to run in polynomial time.
• given an incorrect answer with probability at most 1/3
BQP (Bound error Quantum Polynomial time) is a complexity class that is the quantum analog of BPP. The
hope is that there are problems that belong to BQP and not in BPP, which will help in realizing quantum
advantage [54].
 Quantum Integer Programming September-December 2020

Lecture 2: Test-set methods - Gröbner Basis (Part 1)

Lecturer: Prof. Sridhar Tayur, David E. Bernal Scribe: Srijan Gupta, Vighnesh Natarajan
Slides, Notebook Video, password: 6.gKTx4+

In any interger program (IP), we have a set of constraints. This set is defined using equality and inequality
constraints. However, the inequality constraints can be converted to equality constraints by adding new
variables. Here, we discuss a way to find all feasible solutions (‘points’) to the set of constraints.
The following section starts with the necessary preliminary mathematical definitions: fun and deep in their
respect.

2.1 Algebraic Geometry & Polynomial Ring

Algebraic geometry uses tools from algebra to study of geometric objects. These geometric objects are
defined by polynomial equations.

2.1.1 Field, Ring, Affine-n space

Definition 1. (Field) A field is a set K along with two operations defined on it: ‘addition’ (a + b) and
‘multiplication’ (a · b) for a, b ∈ K. Both these operations satisfy the usual field axioms: associativity,
commutativity, distributivity, existence of identity and existence of inverse.
Examples of fields: the set of rational numbers Q, and real numbers R.
Definition 2. (Ring) A ring is a set R which is closed under addition and multiplication. This set forms
an abelian group under addition, and it satisfies associativity for multiplication and distributivity.
Definition 3. (Affine n-space) The affine n-space over the field K, denoted by AnK or simply An , is the
set of all n-tuples pf elements of K. An element P = (a1 , a2 , . . . , an ) ∈ An is called a point, and ai ∈ K are
called the coordinates of P .

2.1.2 Monomials, Polynomials, Polynomial Ring

Let K be a field and x = x1 , x2 , . . . , xn be indeterminates.

Definition 4. (Monomial) A monomial xα is given by
xα := xα1 α2 αn
1 x2 · · · xn (2.1)
where αi ’s are non-negative integers.
Definition 5. The tuple α = (α1 , α2 , . P
. . , αn ) is called the multidegree or exponent vector of the monomial.
The degree of a monomial, deg α = i αi
Definition 6. (Polynomial Ring) The polynomial ring K[x] is (a ring of polynomials) the set of all
polynomials p(x) (finite linear combinations of the monomials xα ) with coefficients in the field K.
 X 
K[x] := p(x) p(x) = cα xα , cα ∈ K (2.2)
α

2-1
2-2 Lecture 2: Test-set methods - Gröbner Basis (Part 1)

Definition 7. The degree of a polynomial is the maximum of the degrees of its constituting monomials.

Definition 8. A term is given by cxα where c ∈ K

2.1.3 Polynomial and Monomial Ideals

First, what is an ideal?

Definition 9. (Ideal) Let R be a ring. A subset I of R is an ideal if

(i) 0R ∈ I,

(ii) if a, b ∈ I then a + b ∈ I, and

(iii) if a ∈ I and k ∈ R then ak = ka ∈ I

√
Definition 10. (Radical of an ideal) The radical of an ideal of a ring R, denoted by rad(I) or I, is
defined as √ n o
I := a ∈ R | an ∈ I for some n ∈ Z+ (2.3)
√
Note. I ⊆ I

Let S be a set of polynomials f such that S ⊆ R[x], where R[x] is a polynomial ring. The set S corresponds
to the set of equality constraints.
The system S generates a polynomial ideal I of the polynomial ring R[x], such that S ⊆ I, and, for
f1 , f2 ∈ I, f1 + f2 ∈ I; for f ∈ I and r ∈ R, rf ∈ I.
Ideal generated by a set of functions. For the polynomials
Ps f1 , . . . , fs ∈ S, define hf1 , . . . , fs i as the
set that consists of all polynomials that are obtained by i=1 hi fi , with hi ∈ R[x].
 s
X 
hf1 , . . . , fs i := f f= hi fi , hi ∈ R[x] (2.4)
i=1

Note that hf1 , . . . , fs i is an ideal, called the ‘ideal generated by S’.

Theorem 1. (Hilbert Basis Theorem) Every polynomial ideal in R[x] is finitely generated.

Definition 11. (Monomial ideal) A monomial ideal is a polynomial ideal that can be generated by
monomials.

A polynomial belongs to a monomial ideal I if and only if all of its terms ∈ I.

2.2 Duality between affine variety and coordinate ring

First the definitions:

Definition 12. (Affine variety) The affine variety V(S) is the set of common zeros z = (z1 , . . . , zn ) of the
polynomials in S.
V(S) := {z | f (z) = 0, ∀f ∈ S} (2.5)
Lecture 2: Test-set methods - Gröbner Basis (Part 1) 2-3

It’s not hard to see why we would be interested in the affine variety V(S). This is the set of all feasible
solutions to the set of constraints S.
Definition 13. (Coordinate ring) Given an affine variety V in an affine-n space over the field K, the
coordinate ring of V is the quotient ring

K[V] := K[x]/I(V) (2.6)

(basically the set K[x] mod I(V)) where I(V) is the ideal formed by all polynomials f (x) with coefficients
in K which are zero at all points of V

There is a duality between affine algebraic varieties (i.e., the loci of common zeros of system of polynomial
equations) and coordinate rings.

affine algebraic varieties ∼ coordinate rings (2.7)

Example 1. Let’s look at the following example of an affine algebraic variety

V = the unit circle in R2 (2.8)

The very same data (the set of points lying on the unit circle) is captured algebraically with the coordinate
ring

R/ x2 + y 2 − 1 (2.9)
2 2
= (polynomials in x and y) mod (x + y − 1) (2.10)

which contains polynomials such as (x − cos θ)(y − sin θ), θ ∈ [0, 2π).
√
Theorem 2. V(S) = V(I) = V( I)

Proof. Let zS be an element of V(S), i.e., ∀f ∈ S, f (zS ) = 0. Thus, for r ∈ R, rf (zS ) = r × 0 = 0, and for
f1 , f2 ∈ S (hence ∈ I), f1 (zS ) + f2 (zS ) = 0. Thus, each element of V(S) is an element of V(I). Also, since
∀f ∈ S also ∈ I, and for each element zI of V(I), √ f (zI ) needs to be 0, which makes zI an element of S.
Hence, V(S) = V(I). The statement V(I) = V( I) can be easily seen to be true.

The ideal I generated by S reveals hidden polynomials. For instance, if one of the hidden polynomials is
the constant polynomial 1 (i.e., 1 ∈ I), then the system S is√inconsistent (because 1 6= 0). To be precise,
the set of all hidden polynomials is given by the radical ideal I.
√
Theorem 3. (Nullstellensatz) I(V(I)) = I, where I(V(I)) is the ideal that contains all polynomials that
vanish for the variety V(I).

This duality, and the reformulation of the set S so as to

√ easily find it’s common
√ roots
√ zS are what we are
primarily interested in. It’s true that V(S) = V(I) = V( I), but S ⊆ I ⊆ I and I is infinite. We seem
to have made our problem worse. Be patient, in √ Sec 2.4 we shall see a finite set of functions (spoiler alert:
the Gröbner basis) which are the generator for I, and have very nice properties.

2.3 Monomial Ordering

Definition 14. (Monomial ordering) A monomial ordering on C[x] is a relation  on Zn+ (i.e., the
monomial exponents), such that:
2-4 Lecture 2: Test-set methods - Gröbner Basis (Part 1)

(i) The relation  is a total ordering.

(ii) If α  β, and γ ∈ Zn+ , then α + γ  β + γ

(iii) The relation  is a well-ordering (every nonempty subset has a smallest element).

There are several term orderings which are of interest to us, for instance:
• Lexicographic (”dictionary”). α lex β if the left-most nonzero entry of α − β is positive. Notice
that a particular order of the variables is assumed, and by changing this, we obtain n! nonequivalent
lexicographic orderings.
• Graded lexicographic. Sort first by total degree, then lexicographic, i.e., α grlex β if |α|> |β|, or if
|α|= |β| and α lex β
• Graded reverse lexicographic. α grevlex β if |α|> |β|, or if |α|= |β| and the right-most nonzero entry
of α − β is negative. This ordering, although somewhat nonintuitive, has some desirable computational
properties.

Example 2. Consider the polynomial ring R[x, y]. In the lexicographic ordering (≺lex ) discussed above, we
have:
1 ≺ y ≺ y 2 ≺ · · · ≺ x ≺ xy ≺ xy 2 ≺ · · · ≺ x2 ≺ x2 y ≺ x2 y 2 ≺ · · · .

The graded lexicographic order gives (≺grlex , and the graded reverse lexicographic order givesn ≺grevlex ).
In a special case consisting of two variables, the ordering coincide:

1 ≺ y ≺ x ≺ y 2 ≺ xy ≺ x2 ≺ y 3 ≺ xy 2 ≺ x2 y ≺ x3 ≺ · · · .

Definition 15. (Leading term) Given a term order on R[x] and a polynomial f in this ring, say the highest
monomial in f with respect to the term order is cα xα , then
• Leading term, LT (f ) = cα xα
• Leading monomial, LM (f ) = xα
• Leading coefficient, LC(f ) = cα

2.4 Gröbner Basis

Just a little more patience, we are getting to it.

2.4.1 Initial Ideal

Definition 16. (initial ideal) Consider an ideal I ⊂ R[x], and a fixed monomial ordering. The initial ideal
of I, denoted in(I), is the monomial ideal generated by the leading monomials of all the elements in I, i.e.,

in(I) := hLM (f ) | f ∈ I\{0}i (2.11)

Given an ideal I = hf1 , . . . , fs i, we can construct two monomial ideals associated with it. We have the initial
ideal in(I), previously defined. We can also consider the monomial ideal generated by the initial monomials
of the generators, i.e., hLM (f1 ), . . . , LM (fs )i. Although we always have hLM (f1 ), . . . , LM (fs )i ⊆ in(I) and
in general these two monomial ideals are not equal.
Lecture 2: Test-set methods - Gröbner Basis (Part 1) 2-5

2.4.2 Gröbner Basis

Finally here. As we saw in the previous section, hLM (f1 ), . . . , LM (fs )i ⊆ in(I). It is possible to produce a
set of generators polynomials for which these two ideals are same. This is the Gröbner basis.
Definition 17. (Gröbner Basis) Consider the polynomial ring R[x], with a fixed monomial ordering, and
an ideal I. A finite set of polynomials {g1 , . . . , gs } ⊂ I is a Gröbner basis of I if the in(I) is generated by
the leading terms of the gi , i.e.,
in(I) = hLM (g1 ), . . . , LM (gs )i (2.12)
Theorem 4. Every ideal I has a Gröbner basis G. Furthermore, I = hg1 , . . . , gs i.

Gröbner bases as defined are not unique. This can be easily fixed by refining the concept to the so-called
reduced Gröbner bases, which are uniquely defined. For simplicity, we use the term ‘Gröbner basis’ to refer
to the reduced Gröbner basis.
√
Continuing from where we left off in section 2.2, the radical ideal I √
has a finite generating set, the Gröbner
basis B, which one might take to be a triangularization of the ideal I. It is, in fact, the generalization of
Gaussian elimination in linear systems.

Properties
√
• V(S) = V(I) = V( I) = V(B)
Hence, the set of all feasible solutions of S is the same as the set of solutions for B.
• If 1 ∈ B then V(S) = φ
Hence, it can be quickly figured out whether there exists a feasible solution to the constraints defining
S. Other algorithms to solve IP problems take very long if there is no feasible solution. Gröbner
basis computation takes a long time when there are a lot of feasible solutions, in which case the other
algorithms are faster. So if we run Gröbner basis computation in parallel to some other algorithm, we
reach an answer much quickly one way or the other.
• The size of V(S) = V(B) can be obtained using staircase diagram (without solving equations).

Figure 2.1: The monomial ideal M = x4 , x3 y 4 , xy 7 , y 10 , with its generators (white circles), standard
monomials (black dots), and irreducible components (shaded circles) [50]

Example 3. Consider the system given by,

S = {x2 + y 2 + z 2 − 4 = 0, x2 + 2y 2 − 5, xz − 1 = 0} (2.13)

We want to solve S. One way to do so is to compute a Gröbner basis for S. We can do this directly by
using the Sympy12 package in Python, which uses the Buchberger’s algorithm. It was specified to use the
12 [Link]/
2-6 Lecture 2: Test-set methods - Gröbner Basis (Part 1)

lexicographic ordering such that x > y > z. The computed Gröbner basis is,

B = {x + 2z 3 − 3z, y 2 − z 2 − 1, 2z 4 − 3z 2 + 1} (2.14)

We can see that the initial system has been triangulized: The last equation contains only the variable z,
while the second has an additional variable, and so on. The variable z is said to be eliminated with respect
to the rest of the variables. For the defined lexicographic ordering, the program tries to isolate z first, then
z and y, and finally x, y, and z (all variables). It is clear that different orderings yield different Gröbner
bases.

2.4.3 How to compute the Gröbner basis?

The following three definitions are useful in defining the Buchberger algorithm, which computes the Gröbner
basis.
Definition 18. (S-Polynomials) Given two multivariate polynomials f, g ∈ R[x], considering the least
common multiple of their leading monomials with respect to an ordering ,

L = lcm(LM (f ), LM (g)) (2.15)

we define, the S-Polynomial as,

L L
S(f, g) = f− g (2.16)
LT (f ) LT (g)
Note. S(f, g) is defined this way so that the leading terms of f, g cancel out.
Definition 19. (Normal or reduced form) A polynomial f ∈ R[x] is reduced or normal with respect to
G = {g1 , . . . , gs } ⊂ R[x] if no monomial of f is contained in the ideal hLM (g1 ), . . . , LM (gs )i, i.e., no term
in f is a multiple of any LM (gi ).

Another formal definition of Gröbner Basis:

Definition 20. For all the pairwise elements of the Gröbner basis G,

∀gi , gj ∈ G, mod (S(gi , gj ), G) = 0 (2.17)

Buchberger algorithm

Input: A set of polynomials F that generates an ideal I

Output: A Gröbner basis G for I

Algorithm:

Initialize G := F
And a set, M := {{fi , fj } | fi , fj ∈ G & fi 6= fj }

while M 6= φ
{p, q} := an element of M
Update M , M := M \{p, q}
Lecture 2: Test-set methods - Gröbner Basis (Part 1) 2-7

Define the S-polynomial, S = S(p, q)

Reduce S with respect to G, h := reduced(S, G)

if h 6= 0:
Add the new pairs {g, h} to M , M := M ∪ {g, h}∀g ∈ G
Add h to the set G, G = G ∪ h (corresponding to which we have already updated M )

The output of the Buchberger algorithm may not be the reduced Gröbner basis.
Definition 21. (Reduced Gröbner basis) A Gröbner basis G is reduced if ∀g ∈ G

(i) LT (g) is not a factor of any element of G\{g}

(ii) LC(g) = 1

Special cases of Buchberger algorithm:

• If the system of equations is univariate, it reduces to Euclid’s division algorithm
• If the system of equations is linear it reduces to Gauss’ algorithm

2.4.4 Applications of Gröbner basis

• Algebraic Geometry
• Coding Theory
• Cryptography
• Invariant Theory
• Integer Programming
• Graph Theory
• Statistics
• Symbolic Integration
• Symbolic Summation
• Differential Equations
• Systems Theory

2.5 Example: Graph coloring using Gröbner basis

Given a graph G(V, E), where V is the set of vertices and E is the set of edges, given k colors, K = {1, . . . , k},
• Can the graph be colored s.t. no two vertices connected by an edge have the same color?
• If yes, what are the possible colorings?
This can be posed as an IP (actually as a SAT since the objective is irrelevant):
min 1 (2.18)
x
s.t.
xij ∈ {0, 1}, ∀j ∈ K, ∀i ∈ V (If V is colored j, xij = 1, else 0) (2.19)
X
xij = 1, ∀i ∈ V (each vertex must have one (and only one) color) (2.20)
j∈K

xuj + xvj ≤ 1, ∀j ∈ K, ∀(u, v) ∈ E (Either none or only 1 of the 2 V ’s of an E have a color j) (2.21)
2-8 Lecture 2: Test-set methods - Gröbner Basis (Part 1)

The constraints can be remodelled in a way which is more suitable for the
Gröbner basis computation as
follows. We map the k colors to the k-th roots of unity, i.e., for ζ = exp ι 2π k , we map the colors {1, . . . , k}
to {1, ζ, . . . , ζ k−1 }. Thus, the statement that every vertex must have one of the colors can be written as the
set of polynomial equations:
Sk = {xki − 1 = 0 : i = 1, 2, . . . , n} (2.22)
We also require that two adjacent vertices xi and xj are assigned different colors. From the previous
discussion we know that xki = 1 and xkj = 1, so xki = xkj or, equivalently, xki − xkj = 0. By factorization we
obtain that:
xki − xkj = (xi − xj ) · f (xi , xj ) = 0
where f (xi , xj ) is a bivariate polynomial of degree k − 1 in both variables. Since we require that xi 6= xj then
xki − xkj can vanish only when f (xi , xj ) = 0. This allows us to write another set of polynomial equations:

SG = {f (xi , xj ) = 0 : (i, j) ∈ E}

Next we combine Sk and SG into one system of equations S. The graph G(V, E) is k-colorable if the Gröbner
basis B of S is non-trivial, i.e., B 6= {1}. If this is not the case, then the graph isn’t k–colorable. Otherwise
the Gröbner basis gives us information about all possible k–colorings of G.

2.6 Algebraic Geometry in Combinatorial Optimization

The BPT method

Proposed by Bertsimas et al. [6].
Given an optimization problem,

argminx f (x) (2.23)

s.t. g(x) = 0 (2.24)
n
x ∈ {0, 1} (2.25)

where f (x), g(x) ∈ R[x], Algebraic Geometry appears naturally.

• The feasible solutions are the variety of the ideal generated by the constraints, additional polynomials
whose roots are {0, 1}, enforcing the variables to be binary: xi (x1 − 1) = x2i − xi and a polynomial
z − f (x) with a new variable z. The ideal is given by,

I = z − f (x), g(x), x2i − xi ⊂ R [z, x0 , · · · , xn−1 ]

• We compute it’s Gröbner basis B with the lexicographic ordering xi > z, which results in B containing
a polynomial g(z) which is just a function of z.
• Solving g(z) = 0 gives the set of all feasible values zsol = {z1 , . . . , zm } of the objective function f (x).
• Choose the optimal value of zsol , zopt = min(zsol )
• Solve rest of the polynomials in B by substituting z = zopt , which gives the solution xsol = xopt
Hence, we have obtained the solution xopt to the optimization problem.
This ends Part 1 of Lecture 2.
 Quantum Integer Programming September-December 2020

Lecture 3: Test-set methods - Gröbner Basis (Part 2)

Lecturer: Prof. Sridhar Tayur, David E. Bernal Scribe: Srijan Gupta, Vighnesh Natarajan
Slides, Notebook Video, password: RV=9Uˆs3

3.1 Conti and Traverso method

The Conti and Traverso method [7] or the CT method is a purely algebraic method of solving integer
programs. This method is an application of the Buchberger algorithm to integer programs. We want to
solve an IP problem with the following constraints.
minimize cT x
x
subject to: Ax = b, x ∈ Zn

A is the coefficient matrix, A ∈ Zm×n

+ , b is the right hand side vector, b ∈ Zm
+ and c is the cost vector with
n
c ∈ Z+ . Here, we must note that (1) b must be in the column space of A for a feasible solution; (2) the
constraint Ax = b gives us a polytope. We create a term ordering with respect to an inner product with
respect to the cost vector c.

Toric Ideal A Toric ideal is a mapping of fields generated by the difference of monomials. For a mapping
K2 {x, y} −→ K3 {x2 , xy, y 2 }, the Toric ideal is < z1 z3 − z22 >. There is a mapping to a variable z, to which
we shall also map the integer problem constraints.

The first step is to create the Toric ideal of the constraints. The constraints of Ax = b are mapped to
a similar variable z, of a Toric ideal, which gives us, for each constraint
ziai1 x1 +ai2 x2 ....ain xn = zibi
Thus, for all constraints, we get,
m Y
n n
a x
Y Y
(zi ij j ) = z bi
i=1 j=1 i=1
Qn aij xj
With this, we define a mapping φ : Q[w1 , ...wn ] −→ Q[z1 , ...zm , z1−1 ...zm
−1
] such that φ(wj ) = j=1 (zi ).
+ −
We can write the columns of the matrix A as aj = aj − aj . This lets us define our ideal,
− +
I = < z aj wj − z aj , 1 − z1 − ... − zm >
Before proceeding further, let us look at a couple of examples for computing the ideal.

Example 1
4x1 + 5x2 + x3 = 37
2x1 + 3x2 + x4 = 20

3-1
3-2 Lecture 3: Test-set methods - Gröbner Basis (Part 2)

It is clear that φ(w1 ) = z14 z22 , φ(w2 ) = z15 z23 , φ(w3 ) = z1 , φ(w4 ) = z2 . A set of feasible solutions satisfy
φ(w1 , w2 , w3 , w4 ) = z137 z220 . Let fi = φ(wi ). Our ideal of consideration becomes

I = < f1 − w1 , f2 − w2 , f3 − w3 , f4 − w4 >

Example 2 We can also consider negative coefficients in the constraints.

2x1 − 1x2 + x3 = 4
−1x1 + 2x2 = 5

It is clear that φ(w1 ) = z12 /z2 , φ(w2 ) = z22 /z1 , φ(w3 ) = z1 . A set of feasible solutions satisfy φ(w1 , w2 , w3 ) =
z14 z25 . Our ideal in this case becomes

I = < w1 z2 − z12 , w2 z1 − z22 , w3 − z1 >

− +
Considering I = hz aj wj − z aj , 1 − z1 − ... − zm i = hxu − xv : Au = Av, u, v ∈ Zn+ i, we define a monomial
relative ordering with respect to the cost c, >c , defined as

α >c β −→ cT α > cT β or
cT α = cT β and α > β

With a Gröbner basis with respect to monomial ordering, we get B ∩ Q[w] is a Gröbner basis of the ideal
I ∩ Q[w] = IA . IA is called the Toric ideal of A. We can observe that this is independent of the right hand
side of the constraints. Thus, for a problem, if we can find the Gröbner basis of the constraint matrix A,
then no matter what the right hand side is, we can quickly solve the problem.

Theorem Let B>c be the reduced Gröbner basis with respect to >c of the Toric ideal IA . Then for any
right-hand side and integer constraint matrix, term order, and non-optimal feasible solution z0 , there is some
pair u, v such that z1 = z0 − u + v is a better feasible solution.

We can use this to generate the reduced Gröbner basis either by serial reduction or by a direct reduc-
tion using python. Using the Gröbner basis, we get our solution of the Integer program. A link to the
examples on python has been added in the references.

3.2 Improvements and Benchmarks in computing Gröbner basis

The main step in these algorithms’ performance is the computation of the Gröbner basis. Sympy takes
around 15 seconds for a 12 vertex graph coloring problem. Mathematica is more optimized, giving the
answer within a second. As the size of the problem grows, we can see that we will get stuck.
There have been various improvements optimizing the Buchberger algorithm. It has a doubly exponential
time complexity and is a crucial step in optimization. Besides, there are optimizations called project and
lift, which optimize the basis vectors’ serial reduction. Optimizations by Charles Faugere are based on
a transformation into linear algebraic problems, where more efficient algorithms can be applied. For the
efficient computation of the Gröbner basis of the Toric ideals, efficient algorithms try to avoid the zero
reduction of the polynomial during back substitution. 4ti2 is the most efficient algorithm to compute this.
As constraints get complicated, we can’t simply apply the BPT and the CT(Conti Traverso) methods on
the entire formulation. Techniques used for complicated constraints are first to relax the problem to a linear
Lecture 3: Test-set methods - Gröbner Basis (Part 2) 3-3

problem and perform a walk back from the solution. This applies to all branch and bound solvers. A key
critical factor is that we always need to stay feasible to find the solution. Finding the Gröbner basis is the
best way to guarantee feasibility. This way lets you iterate through a lattice of feasible points. This leads
us to the concepts of test sets.

3.3 Test-set methods

Given an integer linear program, there exists a finite set called the test set denoted by T = {t1 , t2 ....tn }
that only depends on A and c, the constraint matrix and the cost vector, that assures a feasible solution
x∗ is optimal if and only if cT (x∗ + ti ) ≥ cT x∗ whenever x∗ + ti is feasible. The essential idea is to create
a separation of constraints, calculate the Gröbner basis, and do a walk back with the complex constraints
acting as a complicated oracle.

A clever observation is that the Gröbner basis of a Toric ideal is a test set. Hence we can see that the CT
method is a test set method. However, things can get more complicated for general constraints, and we will
need to perform walk backs for better optimization.
The following is a comparison between the classic divide and conquer methods and test set-based methods.

Methods based on divide-and-conquer Methods based on test-sets

• Branch-and-Bound algorithms • Algorithms based on “augmentation”
• Harness advances in polyhedral theory • Harness tools from algebraic geometry
• Global optimality guarantees • Global convergence guarantees
• Efficient codes and exponential complexity • Few implementations available, Polynomial or-
acle complexity once we have test-set

The requirements to work with test sets on an Integer program are

• An initial feasible solution
• An oracle to compare the objective function
• A test set (set of directions)
Given an objective, the test set will give us directions to improve the objective. If we cannot improve the
objective, we have reached optimality. The Gröbner basis test-set only depends on the constraints and
objective and can be computed for equality constraints with integer variables.

Example
   
1 1 1 1 10
Ax = 

x =  
  
0 1 2 3 15

Figure 3.1 is a 2D section of a 4D lattice. From a feasible point, we are using the test set to iterate to
optimality. The test set direction vectors g1, g2 and g3 have been obtained by finding the Gröbner basis of
the allied problem’s objective and constraint matrix. The computation of this test set is the NP-hard part.
The rest of the problem is polynomial, as discussed above.
3-4 Lecture 3: Test-set methods - Gröbner Basis (Part 2)

Figure 3.1: 2D slice of a 4D lattice. Improvement by following directions of the test set to optimality

Chance constraints are a particularly difficult constraint to deal with. This is tackled by relaxing these
constraints to find a gröbner basis and using the chance constraints as a complex oracle during the walk
back. Regular MINLP solvers fail on chance constraints. However, chance constraints can be reformulated
and tackled using ILP solvers. However, these ILP solvers struggle and are only competitive in verifying
optimality.

Example Process reliability

We have a series-parallel process, shown in Fig. 3.2, and we want to minimize cost, while satisfying a reliability
threshold. With any unit failing, we want to ensure the entire system does not fail with a certain reliability.

Figure 3.2: Series parallel process [55]

Hence, back-up units are added, which take up the cost, but also improve reliability. The formulation of this
problem is
minimize Σi,j cij xij
x
subject to: 0 ≤ xij ≤ uij
Σj xij ≥ 1
xij ∈ Z
n kj
Y Y
(1 − (1 − rij )xij ) = R(x) ≥ Ro
i=1 j=1
Lecture 3: Test-set methods - Gröbner Basis (Part 2) 3-5

Here, cij is the cost of unit i, j, xij is the number of backups of a particular unit, uij is the maximum number
of each unit, ri j is the reliability of a unit, and Ro is the reliability threshold. This problem is relaxed to an
LP

minimize Σi,j cij xij

x
subject to: xij + tij = uij
Σj xij − dij = 1
xij ∈ Z, i ∈ I, j ∈ J
k
Σni Σj j cij xij = b

An interesting part about this problem is that the Gröbner basis can be calculated analytically.

B<c = {xik di − tik bcik , xiq tip − xip tiq bciq −cip }
 Quantum Integer Programming September-December 2020

Lecture 4: Graver Basis

Lecturer: Prof. Sridhar Tayur, David E. Bernal Scribe: Parth S. Shah
Slides, Notebook Video, password: 49&WKKAR

4.1 Introduction

This lecture introduces the concept of Graver basis, and its application to the field of Integer programming as
a test set. Rather than a rigorous mathematical definition, we focus on a somewhat intuitive understanding
of Grave basis.

4.2 History

The original paper on the Graver basis was written in 1975 by Jack E. Graver [56]. Later, applications of
Graver basis to more general objectives in integer programming were demonstrated by Onn in 2010 [57].
Recently there has been increasing interest in this method of solving integer programs, which inspired [58].

4.3 Hilbert Basis

Before getting into Graver basis, it is useful to first understand the Hilbert basis. This other basis will be
used later to describe the Graver Basis.
There are 2 methods of understanding the Hilbert Basis:
• First method: Given a cone C, F is defined as the points within the cone with integral co-ordinates,
i.e., F = C ∩ Zn . Now, H is considered an integral generating set of F if for every x ∈ F, there exists:
n
X
{h1 , .., hn } | x = λi hi , λi ∈ Z+ .
i=1

If such a set of hi is minimal with respect to inclusion, then this set is an integral basis. If this set is
finite, then it forms a Hilbert basis.
• Second method: Considers solving a homogeneous system of linear equations like Ax = b. The form
of constraint here are the same as those used in our previous lecture on integer programming. We want
to find the kernel of the matrix A:
n
ker(A) = {x ∈ Z+ | Ax = 0}.

The Hilbert basis of the matrix A is the subset of the above kernel, where the kernel is minimal with
respect to inclusion (non-zero is implicit from the definition of the kernel).
From a more geometric point of view, the Hilbert basis is pretty much like a lattice basis for F. It is the set
of (integral co-ordinates) vectors in F whose linear sums can give any vector in the space. A good animation
of the same is provided in the video lecture as well [59].

4-1
4-2 Lecture 4: Graver Basis

In the past (see for instance [60]), integral basis have been used to solve integer programs using pivoting
methods. Such methods differ from branch and cut, which is used today. However, there is potential to
bring integral basis methods into the present day. These integral basis method are very efficient in verifying
optimality. However, the test-set (introduced in lecture 2 and required to test optimality) may be very large
and hard to calculate classically.

4.4 Graver Basis

Now we can understand the basic idea of a Graver basis. Let Oj be the j th orthant of a Rn vector space.
Consider a constraint matrix A defined in the same space. Let:

Hj (A) = H(A) ∩ Oj

Given this, the Graver basis is defined as:

[
G(A) = Hj (A)
j

Again, the geometric intuition here is that we consider the Hilbert basis over each of the orthants (every
orthant is a cone). The union of those gives us the Graver Basis [59].

4.4.1 Some notation

• Recollect how we had an ordering for the Gröbner basis in lecture 2. Similarly, we need to define an
ordering for the Graver basis as well.
• u, v ∈ Rn are conformal if both the vectors belong to the same orthant, i.e. uj vj ≥ 0 ∀ j = {1, 2, .., n}
• For u, v ∈ Rn , u v v if u and v are conformal and v is farther from the origin than u
• L(A) is basically the kernel of the matrix A. This has the property that any element z belonging to
it has a v representation in terms of the Graver basis:
X
z= αi gi , αi ∈ Z+ , gi ∈ G(A), gi v z
i

4.4.2 Graver and Gröbner basis

• It turns out that the Gröbner basis and the Graver basis can be related in some way. Given an objective
c ∈ Zn and a toric ideal I for the constraint matrix A, we denote the Gröbner basis as B(A, c)
• We also define the Universal Gröbner basis as the union over the set of all objective functions as:
[
U(A) = B(A, c)
c∈Zn

• Now it turns out that this Universal Gröbner basis is contained in the Graver basis upto a negating
vector, i.e.,
U(A) ⊆ G(A)
• Sometimes these two are the same, specifically in cases where the matrix A is totally unimodular, or
it is a Lawrence lifting matrix. A totally unimodular matrix is a square matrix with every entry in the
matrix: 0, -1, or +1, and the determinant being -1 or 1. A unimodular matrix just does not have a
restriction on the entries except that they are integers.
Lecture 4: Graver Basis 4-3

4.4.3 Usage of the relation between the Graver and Gröbner basis

Using the relation shown in the previous subsection, there is a method to find the Graver basis for a matrix by
computing its Gröbner basis instead. What is done is that the matrix is ”Lawrence lifted”. Mathematically
what we do is consider the matrix Λ made from the original constraint matrix A as:
 
A 0
Λ(A) = 



In In

where In is the Identity matrix of the appropriate order. Its toric ideal (described in the previous lecture) is
+ − − +
given by IΛ(A) = xu y u − xu y u where u ∈ L(A). The y terms actually come from the identity matrices.
Interestingly, for such a matrix we have:

G(Λ(A)) = U(Λ(A)) = B(Λ(A), c)

for any ordering c defined on the vector space spanned by x, y. To compute the Graver basis, one simply
computes the Gröbner basis for this new matrix Λ where yi is now 1(remember this comes from the identity
matrices).
Code for the above algorithm can be seen in the first part of this notebook. The code uses the CT method
described in lecture 2 for computing the Gröbner basis. As we saw in the previous lecture, calculating the
entire Gröbner basis for a large matrix is very hard to do efficiently on a classical machine. Hence we look
at other algorithms.

4.4.4 Normal form algorithm

The normal form algorithm returns the ”normal form” of a vector and a set belonging to the same Lattice.
The normal form r of s ∈ L and a set G ⊂ L is related to them as:
X
s= αi gi + r, αi ∈ Z+ , gi ∈ G, s v r, gi 6v r ∀gi ∈ L
i

Basically we try to to find ”residual” or remainder of the vector with respect to the set, and then sift through
this remainder to find elements that won’t leave any. The algorithm can is described below, but note that
it very heavy computationally:
Input:s, G
Output:r
Initialise: s 7→ r
Loop while ∃g ∈ G | g v r:
Do: r-g 7→ r
End of loop
Return: r
This is not an algorithm to get the Graver basis, but rather helps in the understanding and implementation
of the following algorithm.
4-4 Lecture 4: Graver Basis

4.4.5 Pottier’s Algorithm

What would happen if the Normal Form algorithm’s input is the kernel of a matrix? We would get the
Graver basis! Why? Because it returns the v-minimal elements of the set. Hence, we arrive at the following
algorithm:

Input:Generating initial set F ⊆ L = ker(A)

Output:Graver basis G of L
Initialise:
G = F ∪ (−F)
S = ∪f,g∈G f + g
Loop while C 6= ∅:
Do: ∀s ∈ C : r ← normalf orm(s, G) and C ← C|{s}
if r 6= 0 then
G ← G ∪ {r} and C ← C ∪ {r + gi } ∀gi ∈ G
End of if
End of loop
Return: G
This is a sort of completion algorithm and is pretty hard to do classically. Intuitively what we are doing
is that we are taking each element in the kernel and finding the residuals again. An element is added to the
Graver basis, only when the element has no residual. The main issue here is that the kernel may contain
many elements that are not v-minimal, and hence we end up spending a lot of unnecessary time on such
elements.
There is an improvement in this algorithm. The computation is not done on the entire set, but only a subset
that conforms to certain restrictions. We don’t need to take every element’s sum with every other element as
shown above, but only if two vectors have the same pattern in the signs of their components do we consider
their sum. This is called the project and lift algorithm.
Note that the most efficient implementation of this algorithm is in a software named 4ti2.

4.4.6 Graver Basis as a Test Set

Test sets were discussed in the previous lecture as well, although there we looked at how the Gröbner basis
may be used as a test set for integer programming. Here we look at how the Graver basis works as a test
set and for what kinds of objective functions. A quick recap of test sets:
• Assume we are given an integer linear program:

min f (x) | Ax = b, l ≤ x ≤ u, x ∈ Zn
x

• Then there exists a test set T = t1 , .., tn that depends only on the constraint matrix A.
• This test set can be used to find more optimal points starting from a non-optimal feasible point xo to
the program. This can be done because:

∃ α ∈ Z+ , i ∈ {1, .., n} | f (xo + αti ) < f (xo )

all the while xo + ff ti still being feasible.

Lecture 4: Graver Basis 4-5

Now, for what kind of objective functions is the Graver basis a test set? Here is a list of some such objective
functions:
• Separable convex minimisation: i fi (cT
P
i x) where all the fi are convex functions.
• Convex integer minimisation: −f (Wx where W ∈ Zdxn and f is convex.
• Norm p minimisation: f (x =|| x − x̂ ||p where x̂ is some vector defined in the domain of x.
• Quadratic minimisation: f (x = xT Qx where Q lies in the dual of the quadratic Graver cone of the
constraint matrix A. This also includes certain non-convex Q 6≥ 0
• Polynomial minimisation:f (x) = P(x) with P being a polynomial of degree d that lies on the dual of
the dth degree Graver cone of the constraint matrix A.
Note that the Graver basis itself does not depend on the objective. The objective just determines whether
or not the Graver basis is a test set for the same. This point is beneficial, it allows us to compute a single
(possibly non-optimal) feasible point and the Graver basis, and just use these to find an optimal point. This
is at the heart of Graver Augmented Multiseed Algorithm(GAMA - lecture 6).
One other very interesting problem for which the Graver basis is a test set is that of the N-fold integer
programming. These are defined as:

min wx | Ax = b, l ≤ x ≤ u, x ∈ Zn
x

where the constraint matrix A is of the form:

 
E1 E1 ... E1 
 
 
E 0 ... 0
 2 
 
A = E(N )
 
=
0 E2 ... 0 
 
 . .. .. .. 
 .. . . . 
 
 
 
0 0 ... E2
| {z }
N

Such problems are widely encountered while dealing with multi-commodity flows, privacy in statistical
databases, and closest string determination. But what makes it interesting is that the Graver basis for
such a problem can be determined in polynomial complexity! This is unlike the other cases where any
current algorithm would still take a worst-case non-polynomial complexity, often exponential.
Using this fact, we could use a similar ”lifting” process as in Lawrence lifting and map any constraint matrix
to such an N -fold matrix. Then we can compute the Graver basis in polynomial time for such a matrix. Of
course, though polynomial in time, all this is still pretty computationally expensive as the lifting to the right
dimension itself can take a long time. Hence it isn’t such a general procedure after all.
4-6 Lecture 4: Graver Basis

4.5 Graver-based Augmentation

Based on the idea of using Graver basis as a test set, we arrive at the Graver-based Augmentation algorithm,
which is sort of a precursor to GAMA, which will be introduced in lecture 6.

4.5.1 Contrasting with the current approaches

Before jumping into the algorithm, it may be useful to discuss a little about the current classical methods
for solving integer programs. These are based on the divide-and-conquer approach, and the one used most
popularly is the Branch-and-Bound algorithm. These sorts of algorithms make use of polyhedral theory
and give global optimality guarantees. Although they have exponential complexity, there is code available
that has been made pretty efficient still.
In contrast, test-set methods rely on ”augmentation”, wherein we start from a non-optimal feasible point
and continually move towards a more optimum point. These make use of the theory of algebraic geometry.
They provide a global convergence guarantee and have a polynomial oracle complexity once supplied with a
test set.

4.5.2 The Algorithm

The algorithm is as follows:

Inputs:
• The constraint matrix A along with the constraint itself Ax = b, l ≤ x ≤ u.
• The test set (Graver basis if applicable), G, corresponding to the constraint.
• The objective function f for the problem.
• One feasible point zo to the problem.
Output: The optimal point zmin to the problem.
Loop while: ∃ t ∈ G, α ∈ Z+ | zo + αt is feasible and f (zo + αt) < f (zo )
Do:
α, t = arg minα,t f (zo + αt)
zo = zo + αt
done
Return:zo

This algorithm is demonstrated in the second half of the code here. A walk-through example is also
provided in the video lecture corresponding to this lecture.
Lecture 4: Graver Basis 4-7

4.5.3 Results against Gurobi

[Link] N -fold programming

Here we showcase some results of using the Graver basis for computing the optimal value for the N -fold
programs (closest string program in specific). The results are compared against those from Gurobi:

Figure 4.1: Fig 3.1: Time vs. Size of problem Figure 4.2: Zoomed in version of figure 4.1

Figure 4.3: Comparison of Test set augmented methods vs. Gurobi[61]

A few key points to note here:

• aug total - total augmentation time is the time to calculate the optimal point given the test set. Note
that it stays fairly constant, or rather grows slowly.
• aug int - The time taken to compute the test set. Notice that this is what takes most of the time in
the augmentation algorithm.
Even though we say polynomial time, the computation is still expensive, mostly because of memory issues.
If we had a way to compute the Graver basis more efficiently, we would do much better!

[Link] Quadratic Cardinality Boolean

Here we consider the problem of the type:

min{cT x + xt Qx | 1T
n x = b}
x

It turns out, for such a problem statement, the Graver basis can be calculated analytically! Following are
some results using the Graver basis, again compared against Gurobi.
Figure 4.4 shows that for convex problems, regardless of what initial feasible point we choose, the optimal
is always reached. In figure 4.5 we see that for non-convex problems getting to the optimal point is trickier:
we land at different points, in the end, depending on our first feasible point. Note that these results are
actually for GAMA, but in the first case, we would get the same results for Graver-best augmentation. In
the second case, we have to use GAMA to get a good optimum.
4-8 Lecture 4: Graver Basis

Figure 4.4: Convex Problem Figure 4.5: Non-convex Problem

Figure 4.6: GAMA convergence from multiple feasible points for convex and non-convex problems[58]

4.6 Food for Thought

Some stuff to think about before the next lecture:

• We see that computing the test-sets, in general, is classically intractable. Hence just directly applying
the same algorithm to every problem instance will not help much. What we could try to do is try
tailoring the algorithm to specific problem instances instead.
• What if we don’t compute the entire Graver basis test set but only a part of it? How good a solution
can we reach?
• Are there (more) special structured matrices that allow for systematic Graver Basis calculation?
Lastly, it worth noting that the Graver basis has applications where current commercial solvers fail, mostly in
the case of not-well-behaved objective functions. Hence the potential for it is immense if we find faster ways
of solving for the basis or find that, with many feasible solutions, only a partial Graver basis is sufficient.
 Quantum Integer Programming September-December 2020

Lecture 5: Ising Model

Lecturer: Prof. Sridhar Tayur, David E. Bernal Scribe: Sashank Kaushik Sridhar
Slides, Notebook Video, password: jKYs*s8@

5.1 Introduction & Background

The Ising model was first introduced by Lenz (1920) to explain the ferromagnetic-paramagnetic phase tran-
sition. It was solved by his student Ising (1925) for the 1-D case, followed by Onsager (1940) for the 2-D
case. The 3-D case was reformulated as a non-planar graph MAXCUT problem, proving the 3-D case to
be NP-Complete. By universality of dimensions and symmetries, Ising models have been shown to have
important implications for non-magnetic phase transitions as well.

5.1.1 Mathematical Definition

X X
H(σ) = − Jij σi σj − µ hi σi
(ij)∈E(G) i∈V (G)

where,
• H: Hamiltonian
• σ ∈ {−1, +1}V (G) : Spin per lattice site
• G = (V, E): Graph or Lattice defining the interactions
• µ: Magnetic Moment
• Jij : Interaction between sites
• hi : External field applied per site
Note:- Depending on the signs of Jij and hi , the interactions can be ferromagnetic or antiferromagnetic. The
partition function gives the probability of each configuration (Z(β)),

P (σ, β) = e−βH(σ) /Z(β)

Figure 5.1: Typical 2-D Ising Lattice.

5-1
5-2 Lecture 5: Ising Model

where,
• β = (kBP
T )−1
• Z(β) = σ e−βH(σ)

5.1.2 Solutions to Ising Models

Solutions exist in 1-D for both circular and free-boundary conditions without external fields. However, in
such solitions np phase transition is observed. With the external field applied, the phase transition is observed
when J = h. In the 2-D case, Onsager’s solution predicts a phase transition. In 3-D, the non-planar graph
is NP-Complete via the MAXCUT formulation, but mean-field approximations exist. Different approaches
can be adapted depending on the problem’s objectives, be it finding important statistical properties or the
exact spin configurations.

5.2 Reformulations

As stated, the Ising model can be recast as a MAXCUT problem, a QUBO, and thereby as an Integer Linear
Program.
For the MAXCUT reformulation, we start with no external field,

X
H(σ) = − Jij σi σj
(ij)∈E(G)
X X X
=− Jij − Jij + Jij
(ij)∈E(V + ) (ij)∈E(V − ) (ij)∈δ(V + )
X X
=− Jij + 2 Jij
(ij)∈E(G) (ij)∈δ(V + )
X
=C +2 Jij
(ij)∈δ(V + )

Taking Wij = −Jij as the weights of the edges, we recast the minimization of the Hamiltonian as a maxi-
mum cut of the graph.
Integer Linear Program:

X X
minx∈{0,1}n Qij xij + Qii xi + c
(ij)∈E(G) i∈V (G)

s.t. xij ≥ xi + xj − 1
xij ≤ xi
xij ≤ xj
xi , xij ∈ {0, 1} ∀ (ij) ∈ E(G), i ∈ V (G)

P P P
with Qij = 4Jij , Qii = 2hi − j∈V (G) (2Jij + 2Jji ), c = i<j Jij − i∈V (G) hi
Lecture 5: Ising Model 5-3

5.3 Solving Ising Models

Ising models can be solved computationally using the Markov-Chain Monte Carlo Metropolis-Hastings Al-
gorithm which we describe next.

5.3.1 Metropolis-Hastings Monte Carlo (MHMC) Algorithm

1. Start with a known spin configuration σ i , with corresponding energy H(σ i ) and temperature T =
(kB β)−1 .

2. Randomly change spin configuration: σ j = σ i + δ

3. Calculate new H(σ j ).
4. Compare to the previous energy and,
- If H(σ j ) < H(σ i ), keep new position h i
j i
- If H(σ j ) > H(σ i ), keep new position if exp − H(σ k)−H(σ
BT
)
≥ Rand[0, 1].

5. Repeat steps 2-4 K times.

5.3.2 Simulated Annealing

The MHMC algorithm can be improved by providing a heuristic to lower the temperatures in a scheduled
fashion. This algorithm is based on a metallurgical process of annealing, hence the name simulated annealing.
If the temperature is lowered slowly enough, the algorithm is guaranteed toh reach the ground istate of the
ln(2) ln(100)
energy profile. The default limits for the temperature are given by β ∈ max{∆E} , min{∆E} . At every
temperature, we run the MHMC algorithm as a nested loop. The outer loop incrementally reduces the
temperature from an effectively high starting point.

5.3.3 Advanced Simulated Annealing

For large graphs, the update δ is a limiting factor in the algorithm’s runtime. If individual spins are updated,
the algorithm can take an exponential time to converge to the global optimum. Hence, we update the spins
in clusters.
• This must be done carefully not to violate energy conservation and ensure that the evolution is ergodic.
Any state can be reached from another via the Markov-Chain.
• To ensure convergence, we generate multiple replicas of the system at different temperatures. After
a fixed number of MHMC updates, the temperatures of two replicas r1 , r2 are exchanged with a
probability:
P (r1 ↔ r2 ) = min{1, exp[(β1 − β2 )(H(σ 1 ) − H(σ 2 ))]}
• Two temperatures are always exchanged if a replica at higher temperature has a lower energy than a
replica with a lower temperature.
• Otherwise, the exchange of the two temperatures is either accepted or rejected using the random
number between 0 and 1
5-4 Lecture 5: Ising Model

Figure 5.2: Traversing the energy landscape through simulated annealing.

5.4 Comparing Heuristics

Comparing different heuristics is generally difficult due to several parameters to tune and differences in the
underlying hardware. Two metrics to consider are Time and Solution Quality. A useful metric is given as
Time To Solution (TTS):
log(1 − s)
TTS(m) = mτ (m)
log(1 − p(m))
• m: Number of times run, or sweeps in Simulated Annealing
• s: Success probability after m sweeps
• p(m): Probability of success to achieve (usually high)
• τ (m): Time it takes to perform a single sweep
Several heuristics exist for MAXCUT and QUBO. These heuristics can be compared with each other
using this metric. Simulated Annealing algorithms, as highlighted below, seem to perform well among other
heuristics.
Lecture 5: Ising Model 5-5

Figure 5.3: Heuristic methods for addressing QUBO problems, with the ones using Simulated Annealing
highlighted [19]
 Quantum Integer Programming September-December 2020

Lecture 6: QUBO - Quantum Unconstrained Binary Optimization

Lecturer: Prof. Sridhar Tayur, David E. Bernal Scribe: Suriyaa Valliapan, EP17B027
Slides, Notebook Video, password:

6.1 Abstract

This lecture introduces us to QUBO and its equivalence to the Ising model. We then move on to reformulating
a general Integer Programming problem as a QUBO problem, ending with a demonstration.

6.2 QUBO Model

QUBO expands to Quadratic Unconstrained Binary Optimization. Quadratic, because it contains expres-
sions of order 2; unconstrained, because there are no constraints; and binary, since all variables are either
zero or one. Equation (1) represents the QUBO model.

P P
minx∈{0,1}n (ij)∈E(G) xi Qij xj + i∈V (G) Qii xi + c
(6.26)
>
= minx∈{0,1}n x Qx + c

The offset c above is irrelevant for the optimization since it is independent of the variables. The problem
above represents an integer program that is non-linear (although it could be linearised). If Q is not positive
semi-definite, the problem may be non-convex.

6.2.1 Quadratic Coefficient Matrix Q

The Qij in equation (1) represents the Quadratic Coefficient Matrix. There is no unique representation for
this matrix. The matrix Q can be upper triangular or complete because of the fact that xi xj = xj xi . We
can also note that the elements on the diagonal can be linearised since x2i = xi if x ∈ {0, 1}. This matrix
can also be interpreted as the adjacency matrix of a graph (weighted edge between node xi and xj ).

6.3 QUBO Ising mapping

The Ising problem can be converted to a QUBO as follows: Transform the variables xi ∈ {0, 1} to σi ∈ {−1, 1}
such that

σi = 2xi − 1,
(6.27)
σi σj = 4xi xj − 2xi − xj + 1.

6-1
6-2 Lecture 6: QUBO - Quantum Unconstrained Binary Optimization

In these new variables, the objective function becomes

P P
minσ∈{−1,+1}n (ij)∈E(G) Jij σi σj + i∈V (G) hi σi + cI =
P P
minx∈{0,1}n (ij)∈E(G) xi Qij xj + i∈V (G) Qii xi + cQ
(6.28)

P P P
Qij = 4Jij , Qii = 2hi − j∈V (G) (2Jij + 2Jji ) , cI = cQ + i<j Jij − i∈V (G) hi

where E(G) is the set of pairs of nodes on the graph G that are connected (Edges), and V(G) is the set of
vertices of the graph G. Similarly, a QUBO problem can be converted to an Ising problem as follows:

xi = (σi + 1) /2
(6.29)
xi xj = (σi σj + σi + σj + 1) /4

Once we have this variable transformation, the objective function becomes,

P P
minx∈{0,1}n (ij)∈E(G) xi Qij xj + i∈V (G) Qii xi + cQ =
P P
minσ∈{−1,+1}n (ij)∈E(G) Jij σi σj + i∈V (G) hi σi + cI
(6.30)

P P P
Jij = Qij /4, hi = Qii /2 + j∈V (G) (Qij /4 + Qji /4) , cI = cQ + i<j Qij /4 − i∈V (G) Qii /2

Thus, we have shown that the QUBO and Ising are equivalent models.

6.4 QUBO-Integer Programming Model

Our goal is to formulate a generic integer programming problem as a QUBO. Here, for the sake of complete-
ness, we talk about QUBO as an Integer programming model.
min x> Qx + c (6.31)
x∈{0,1}n

This non-integer program is solvable using commercial INLP solvers such as Gurobi. To simplify things,
we convert this to an ILP by adding new variables xij = xi xj . The non-linearity in these can be captured
by using linear inequalities as shown below. Experimental results show that this is the most efficient ILP
formulation of QUBO.

P P
minx∈{0,1}n (ij)∈E(G) Qij xij + i∈V (G) Qii xi + c

(6.32)
s.t. xij ≥ xi + xj − 1, xij ≤ xi , xij ≤ xj ∀(ij) ∈ E(G)

xi ∈ {0, 1} ∀i ∈ V (G), xij ∈ {0, 1} ∀(ij) ∈ E(G)

Lecture 6: QUBO - Quantum Unconstrained Binary Optimization 6-3

6.5 General Integer Programming Problem as a QUBO

A generic integer programming problem can be converted to a QUBO model as follows. We convert all
the integer variables to binary variables, a process called Binarization. All non-linearities are converted to
quadratic terms, a process called quadratization, and finally one performs unconstraining, i.e, removing all
the constraints by introducing penalties in the objective function.

6.5.1 Binarization

Our goal is to transform all integer variables to binary variables. Let us denote a single positive integer
variable by y which can take values in the range 0, 1, .....y 0 . Now, y can be written in terms of binary
variables x as,
X d
y= kj xj = k> x, kj ∈ Z+ , xj ∈ {0, 1} (6.33)
j=1

with d being the width of the integer encoding. The most obvious choice for expressing an integer variable
in terms of binary variable is the binary encoding where kj = 2j−1 , d = blog2 (ȳ)c. The main issue with
this is that the coefficients become very huge which our quantum annealer may not be able to handle with
a good precision.
The other choice is to use an unary encoding where we set all the kj = 1 and the width d becomes y 0 .
The problem with this is that the width becomes very large, again causing issues.
The choice which fits best here is the bounded encoding where we find an encoding with an upper bound
for the coefficients (µ << y 0 ). We will not go into further details of this although we just state it here for
the sake of completion.

if ȳ < 2blog(µ)c + 1
h i (6.34)
0 1 [log(ȳ)]−1
P[log(ȳ)] i−1
k = 2 ,2 ,...,2 , ȳ − i=1 2

and,

if ȳ > 2blog(µ)c + 1
j k (6.35)
Pρ i−1 v
ρ = blog µc + 1, v = ȳ − i=1 2 , and η = µ


2i−1




 for i = 1, . . . , ρ


ki = µ for i = ρ + 1, . . . , ρ + η (6.36)





 v − ηµ
 for i = ρ + η + 1 if v − ηµ 6= 0

Noise is a significant issue in quantum annealers. Karimi et al., in their paper on ’Practical Integer-to-
Binary mapping for quantum annealers’, showed that bounded encoding gave the best performance, as can
be seen from figures (1) and (2). The probability of getting the right solution decreases exponentially with
the noise’s standard deviation for the binary case, as seen in figure (1). In contrast, figure (2) shows that the
6-4 Lecture 6: QUBO - Quantum Unconstrained Binary Optimization

unary encoding takes much time for reaching the solution. It grows with increase in the standard deviation
of the noise. Hence, we can see that bounded encoding is the better of the three.

Figure 6.1: The scaling of success probability as a function of standard deviation of noise

Figure 6.2: The scaling of time to solution as a function of standard deviation of noise

6.5.2 Quadratization

We state a theorem which is essential for our further discussion.

[Link] Theorem

Any pseudo-Boolean function defined as f : {0, 1}n → R can be written uniquely as a sum of multi-linear
functions.

X X X
f (x) = a0 + ai xi + aij xi xj + aijk xi xj xk + . . . (6.37)
i ij ijk

We define new variables xij . We do that for all terms of order greater than 2 in the above function
continuously until we are left with just quadratic terms. Hence, by using the above theorem and this
Lecture 6: QUBO - Quantum Unconstrained Binary Optimization 6-5

observation, we can conclude that any such objective function can be converted to a quadratic polynomial.
In an earlier section, we had defined the xij as xi xj and imposed some additional constraints to capture the
non-linearity. Instead, here we define a function H:

H(x) = 3xij + xi xj − 2xji xi − 2xij xj (6.38)

which satisfies H(x) = 0 if and only if xij = xi xj . Hence we can add H(x) = 0 to our constraints instead of
the three inequalities, thus simplifying our problem further. The variables xij are called ancillary variables.
We could also represent these new constraints through graphs.

Figure 6.3: Representation of the constraint H(x) = 0 through a graph. Here the edges denote the
coefficient of the product of the variables it connects in the constraint function. The value inside the red
circle represents the coefficient of just the variable corresponding to that node.

We had written H(x) for the case when we had terms of the form xi xj which is equivalent to logical AND.
We can come up with similar expressions for OR and NOT gates, and once we have the constraints for these,
we could come with expressions for everything else. This idea can also be easily extended to higher orders,
as shown in figure (4).

Figure 6.4: Representation of constraints for higher order problems using graphs

6.5.3 Unconstraining

Previously, we had removed the constraints in our linear programming problem using the concept of La-
grange multipliers. Using that concept, we moved our constraints into the objective function by weighting/
penalizing them with Lagrange multipliers. To demonstrate this,

minx≥0 cT x
(6.39)
s.t Ax = b

can be written as,

minx≥0 cT x + λT (b − Ax) (6.40)
6-6 Lecture 6: QUBO - Quantum Unconstrained Binary Optimization

As explained in the previous lectures, we can turn this optimization problem into a problem of finding these
Lagrange multipliers. For linear programs, we had this very nice principle of strong duality. This principle
states that the optimal objective of our initial (primal) problem is the same as the optimal objective of the
dual problem,

maxλ L = maxλ λT b
(6.41)
T T
s.t λ A ≤ c

In other words,
λ∗T b = cT x (6.42)

However for integer programs, strong duality doesn’t hold. However, we can still use the multipliers to
convert the problem into an unconstrained one, as demonstrated in the following example. In this example,
we move constraints into the objective (note the signs of the multipliers)

minx∈{0,1}n cT x

s.t g(x) ≤ 0 (λ) (6.43)

h(x) = 0 (ρ)

can be written as,

minx∈{0,1}n cT x + λT g(x) + ρT h(x)

(6.44)
where λTi ≤0

6.5.4 Putting everything together

Now, using all the tools explained above, we demonstrate how to reformulate an integer programming model
as a QUBO.

• Step 1: Binarize

m minx∈{0,1}n cT x

s.t g(x) ≤ 0(λ) (6.45)

h(x) = 0 (ρ)

• Step 2: Quadratize (getting new variables y)

miny∈{0,1}n dT y

s.t. y> Gy ≤ 0 (6.46)

y> Hy = 0
Lecture 6: QUBO - Quantum Unconstrained Binary Optimization 6-7

• Step 3: Transform inequalities into equalities

Introduce slack variables s to convert the inequalities to equalities so that we get, (Note that by scaling
the constraints, the slack variables become integers)

miny∈{0,1}n ,s∈Z+ d> y

s.t. y> Gy + s = 0 (6.47)

y> Hy = 0 (ρ)

• Step 4: Unconstraining
Considering that we can only deal with discrete variable in QUBO, we cannot change the value of
the multipliers. So, we set the multipliers to be constant penalty factors. In our case, Integer linear
equalities become penalties by squaring the constraints (We square them to avoid the case where the
minimisation of the term leads to negative infinities). So, we get,
min d> y + ρ(y> Gy + s)> (y> Gy + s) + λ(y> Hy)> (y> Hy) (6.48)
y∈Z+ ,s∈Z+

6.5.5 Example

Example optimization problem:

minx≥0 cT x
(6.49)
s.t Ax = b
Let us represent the constraint as HA = 0 and the objective function by HB
m n
!2
X X
HA = ρ Aij xi − bj (6.50)
j=1 i=1

n
X
HB = ci x i (6.51)
i=1
Hence, the total energy function of the QUBO can be written as, H = HA + HB . The next question
we encounter is how to determine the penalty factor ρ. These factors reflect the weightage given to the
constraints. We want it to be such that any infeasible solution that has the lowest value of HB should
still have an overall higher value of H than the best feasible solution (satisfies the constraints and has the
lowest value of HB ). We express this quantitatively as follows: we want even a reasonably small enough
infeasibility to be heavily penalized, this is required because our quantum annealers are noisy they have a
limited precision and a limited range of fidelity
max
∆HB
ρ≥ min
(6.52)
∆HA
where,
max
Pn
∆HB = i=1 max {ci , 0}

(6.53)
Pn
= minσi ∈{0,1},j max 1, 21 i=1 (−1)σi Aij
min
 

∆HA
6-8 Lecture 6: QUBO - Quantum Unconstrained Binary Optimization

6.6 Demonstration

In the Notebook a couple of simple problems formulated as QUBO are solved using simulated annealing.
• Problem 1

minx 2x0 + 4x1 + 4x2 + 4x3 + 4x4 + 4x5 + 5x6 + 4x7 + 5x8 + 6x9 + 5x10
   
1 0 0 1 1 1 0 1 1 1 1
   
   
s.t. 
0 1 0 1 0 1 1 0 1 1 x = 1
  (6.54)
   
   
0 0 1 0 1 0 1 1 1 1 1

x ∈ {0, 1}11

In order to define the Q matrix, we write the problem

minx c> x

[Link] = b (6.55)

x ∈ {0, 1}11

as follows:

minx c> x + ρ(Ax − b)> (Ax − b)

(6.56)
x ∈ {0, 1}11

Exploiting the fact that x2 = x for x ∈ {0, 1} we absorb linear terms in the diagonal of the Q matrix.

ρ(Ax − b)> (Ax − b) = ρ(x> (A> A)x − 2b> Ax + b> b) (6.57)

Pn
For this problem in particular, one can prove the the penalization factor is given by ρ > i=1 |ci |,
therefore we choose this bound + 1.
Remarks:
– The annealer gets stuck in a local minima majority of the time.
– The constraints may not be strictly satisfied as we have just added these terms as penalties to
the objective function. The whole expression being minimized does not necessarily mean that
the constraints have to be satisfied. For constraints to be satisfied exactly, we have to invoke
tools from algebraic geometry (G”oebner basis). Consider the quadratic polynomial where the
goal is to solve (obtain its zeros) Hij as a QUBO. This technique will allow us to eliminate slack
variables.
Hij := Qi Pj + Si,j + Zi,j − Si+1,j−1 − 2Zi,j+1 (6.58)

Now, compute the Gr”oebner basis of the system,

S = {Hij } ∪ x2 − x, x ∈ {Pj , Qi , Si,j , Si+1,j−1 , Zi,j , Zi,j+1 }


(6.59)
Lecture 6: QUBO - Quantum Unconstrained Binary Optimization 6-9

and look for a positive quadratic polynomial.

X
+
Hij = at t (6.60)
t∈B|deq(t)<2

+
Global minima of Hij are the zeros of Hij . (Solve for at )
• Problem 2
G(V, E) has 12 vertices and 23 edges. We ask if the graph is 3–colorable. (Graph Coloring problem)
Remarks:
– Because of precision issues in the translation to QUBO, we obtain very tiny coefficients 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.

6.7 Summary
• Introduction to QUBO model and a brief description of the model.
• A QUBO model is equivalent to the Ising problem.
• QUBO problem is an INLP, and it may or may not be convex.
• A general Integer Programming Problem can be reformulated as a QUBO by binarization of the
variables, Quadratization of the Objective and the constraints, converting the inequalities to equalities
using slack variables, and then unsconstraining by adding the constraints as penalties to the Objective
function.
• A couple of examples and the limitations, difficulties, and the subtleties involved.
 Quantum Integer Programming September-December 2020

Lecture 7: GAMA - Graver Augmented Multiseed Algorithm

Lecturer: Prof. Sridhar Tayur, David E. Bernal Scribe: Gautham Umasankar
Slides, Notebook Video, password: 3&UUXKgf

7.1 Introduction - Hybrid Quantum-Classical Algorithms

A naive way of solving an integer program given by:

min f (x) (7.61)
Ax = b (7.62)
l≤x≤u (7.63)
n
x∈Z (7.64)
T
is to first quadratize f (x) as x Qx and add the constraint using a Lagrange multiplier to get:
min xT Qx + λ(Ax − b)T (Ax − b). (7.65)
However, this has balancing issues and makes the objective even more complicated by increasing its size.
We want to take advantage of the annealer, but more intelligently.

7.1.1 GAMA: A Hybrid Quantum-Classical Optimization Algorithm

GAMA is an algorithm where the classically expensive Graver Basis computation is done by Ising solvers
like the D-Wave Quantum Annealer, Simulated Bifurcation Machine, etc. The Ising solver is also used to
find many initial feasible solutions, which are then augmented using the Graver Basis to find the optimal
solutions.
The defining feature of this algorithm is that it uses the objective function just like an oracle call. The
objective function doesn’t enter the QUBO formulation at all. This allows GAMA to be applied to complex
objective functions. A broad outline of the steps involved in GAMA: [62]

1. Calculate the Graver Basis (Quantum, with Classical post-processing)

2. Find many initial feasible solutions or ”seeds” (Quantum)
3. Augment each feasible solution using the Graver Basis computed in step 3 (Classical)

7.2 Test sets in Optimization, a recap

Given an integer program:

min f (x) (7.66)
Ax = b (7.67)
l≤x≤u (7.68)
x ∈ Zn (7.69)

7-1
7-2 Lecture 7: GAMA - Graver Augmented Multiseed Algorithm

There exists a finite set T = {t1 , ..., tN that depends only on A such that given a feasible, non-optimal point
xo , it satisfies for some α ∈ Z+

1. f (xo + αti ) < f (xo )

2. xo + αti is feasible

7.2.1 Some Definitions

• The lattice integer kernel of A is given by:
L(A) = {x|Ax = 0, x ∈ Zn , A ∈ Zm×n } r {0}
• Partial ordering: ∀x, y ∈ Rn x v y iff xi yi ≥ 0 &|xi |≤ |yi | ∀ i = 1, 2, ...n. If this condition is satisfied, x
is said to be conformal to y.
The Graver Basis is a very general and widely used test set. The Graver basis of A is defined by:

G(A) ⊂ L(A) ⊂ Zn (7.70)

such that G(A) is v - minimal ie, @x, y ∈ G(A) s.t. x v y.

7.2.2 Applicability of the Graver Basis

The Graver Basis can be used to augment feasible solutions and reach the optimum for the following objective
functions f (x):
• Separable convex minimization: i fi (cTi x) with fi convex
P
• Convex integer maximization: −f (W x) where W ∈ Zd×n and f () is convex
• Norm p minimization: f (x) = ||x − x̂||p
• Quadratic minimization: f (x) = xT Qx where Q lies on the dual of the quadratic Graver cone of A
(This includes certain non-convex Q  0
• Polynomial minimization: f (x) = P (x) where P is a polynomial of degree d that lies on cone Kd (A),
the dual of the dth degree Graver cone of A

7.2.3 Heuristics: the power of ”Multi-seed”

Mathematically, we need the complete Graver basis, a single feasible solution and an objective function of
the form mentioned above to reach optimum points. However, Ising solvers might yield only a partial Graver
basis and we also apply our algorithm on more complex objective functions, other than those described in
the previous section. How does the Algorithm then work? Answer: Having multiple starting points or
”multi”-seeding. Its advantages:
• We might obtain only a partial Graver basis. Since we have a large number of starting points, this
increases the likelihood that at least one of them reaches the optimum.
• The objective function might be complicated. However, it can be approximated by our wide classes of
functions given above. With a large number of seeds, we may cover various parts of the feasible space,
and the odds of not reaching the minimum are quite low.
Lecture 7: GAMA - Graver Augmented Multiseed Algorithm 7-3

7.3 Quantum Classical Method to find the Graver Basis

Find a binary (or unary) encoding for the variables to convert form an IP to a QUBO

The integer program: (7.71)

n m×n
Ax = 0 x ∈ Z A∈Z (7.72)
min xT QI x, QI = AT A, x ∈ Zn (7.73)
(7.74)
The Encoding: (7.75)
T
x = [x1 x2 · · · xn ] xi ∈ Z (7.76)
xi = eTi Xi (7.77)
xTi = [Xi,1 Xi,2 · · · Xi,ki ] ∈ {0, 1} ki
(7.78)
eTi = 20 21 . . . 2ki
 
(7.79)

Set a lower bound to allow for negative kernel values and reformulate:

    
 Lx1   eT1 0 T
... 0  X1 
    
    
 L
 x2
 
  0 eT2 ... 0 
 X2 
x = L + EX = 
 .
+   (7.80)
 ..
  .. .. .. ..   .. 

 
  . . . . 

 . 

    
    
Lxn 0 0 ··· eTn X0

Reformulate the original QUBO using this:

X T QB X, QB = E T QI E + diag 2LT Q1 E


min (7.81)
X ∈ {0, 1}nk , QI = AT A (7.82)

After reformulation, solve it with an Ising solver and invert previous transformations to obtain the kernel of
A. After this, v-minimal filtering done classically will yield a (partial) Graver Basis.

7.3.1 Finding the Graver Basis in 10 steps - A Toy Example

This example will walk us through the abstract procedure that was described above.

[Link] Matrix to Quadratic Unconstrained Integer Optimization (QUIO)

Consider a Constraint Matrix A given by:

 
A= 1 2 1 (7.83)

Ax = 0 (7.84)
7-4 Lecture 7: GAMA - Graver Augmented Multiseed Algorithm

We optimize the quadratic version of this:

min xT QI x (7.85)
 
 1 2 1 
 
QI = AT A = 
 
 2 4 2 
 (7.86)
 
 
1 2 1

[Link] Encode

We can use a two bit encoding as shown:

 
e= 0
2 2 1 (7.87)

We will also shift one step to the left L = −1 to cover negative values.

[Link] Obtain Encoding Matrix

With all this in place, the encoding matrix is

 
 1 2 0 0 0 0 
 
 
E=
 0 0 1 2 0 0 
 (7.88)
 
 
0 0 0 0 1 2

[Link] Obtain Encoded Equation

Now, we have the transformation from x to X given by:

x = L + EX (7.89)
 
 X1 
 
 
     
 X2 

 x1   −1   1 2 0 0 0 0  
 


 
 
 
 

 X3 

 x  =  −1  +  0 0 1 2 0 0    (7.90)
 2      


 
 
 
 

 X4 

x3 −1 0 0 0 0 1 2  



 X5 

 
 
X6
Lecture 7: GAMA - Graver Augmented Multiseed Algorithm 7-5

[Link] Formulate QUBO

The QUBO can be formulated as follows:

min (L + EX)T QI (L + EX) (7.91)

(7.92)
T T T


min X E QI E + 2 diag(L QI E) X (7.93)

Which boils down to:

 >   
 x1   −7 2 2 4 1 2  X1 
    
    

 x2 


 2 −12 4 8 2 4 
 X2 

    
    

 x3 


 2 4 −12 8 2 4 
 X3 

min 








 (7.94)

 x4 


 4 8 8 −6 4 8 
 X4 

    
    

 x5 


 1 2 2 4 −7 2 
 X5 

    
    
x6 2 4 4 8 2 −12 X6

[Link] Map to Ising Variables

In this simple transformation, S = 2X − 1, the Ising variable is S.

[Link] Reformulate as Ising problem

By substituting the above in X T QX, we get a minimization of the form:

min S T JS + hT S (7.95)

[Link] Solve Ising model and convert back to QUBO variables

Now, we call our preferred Ising solver. In this lecture, we shall treat this as a black box which gives us
a large number of optimal and near optimal samples for the given Ising Hamiltonian. We want as many
zeros as we can get, so that we can extract more Graver elements. Hence, we can exploit the annealer (the
sampler) since it can sample widely from the solutions.
The DWave quantum annealer gives some optimal and sub-optimal solutions. Of the latter, around 90%
of them have small overall sum-errors. We can postprocess these by adding and subtracting different sub-
optimal vectors to get zero to maximize the number of Graver elements we procure. D-Wave performs really
well for embeddable binary matrices and for narrow truncated band variables.
7-6 Lecture 7: GAMA - Graver Augmented Multiseed Algorithm

From this, we can convert back to QUBO variables from Ising variables (X = (S + 1)/2, to get:
 
 0 0 1 1 0 0 1 
 
 

 0 0 0 0 1 1 1 

 
 

 0 0 0 1 0 1 0 

[x] = 


 (7.96)

 0 1 0 0 0 0 0 

 
 

 0 0 1 1 0 0 1 

 
 
1 0 1 0 1 0 0

[Link] Decode and recover Kernel

Invert the previously made transformations and recover x as a vector of integers.

 
 −1 −1 0 0 1 1 2 
 
 
x = L + EX = 
 0 1 −1 0 −1 0 −1 
 (7.97)
 
 
1 −1 2 0 1 −1 0

[Link] Convert Kernel to Graver Basis

Now we do v-minimal classical filteration to get G(A)

 
 0 1 1 2 
 
 
G(A) = 
 −1 −1 0 −1 
 (7.98)
 
 
2 1 −1 0

Negative basis elements are also part of the Graver basis, due to the symmetry in the constraint. So we can
add −G(A) to our collection of Graver elements.

7.4 Finding feasible solutions and Augmenting

Along the same lines, we can find several feasible solutions by similarly reformulating the constraints as
shown below:
Ax = b l≤x≤u (7.99)
T
min X QB X, (7.100)
T T T


QB = E QI E + 2 diag L QI − b A E (7.101)
nk T
X ∈ {0, 1} , QI = A A (7.102)
Lecture 7: GAMA - Graver Augmented Multiseed Algorithm 7-7

Solving this QUBO using an Ising solver will give us many feasible solutions and starting points. Adaptive
centering and adaptive encoding width can be used for the feasibility bound, which makes the process faster.
Once these points are found, we can augment each point in parallel using our (Partial) Graver Basis. These
augmentations will lead us to the optimum point with high likelihood.

7.5 Demonstration of a Capital Budgeting problem

In this notebook, the problem shown below was solved using GAMA:
v
n u n
X u1 − ε X
min − µi xi + t σi2 x2i (7.103)
i=1
ε i=1

Ax = b x ∈ {0, 1}n (7.104)

Key insights from the demonstration:

• If we use a very small subset of the Graver Basis, we don’t reach the optimum
• On the other hand, using anywhere between 10-100% of the Graver Basis yielded very similar results
over the entire range, provided we start at multiple feasible solutions
• Augmenting with the entire Graver Basis is time consuming and augmenting with a small subset leads
to erroneous solutions. An appropriate trade-off has to be made. However, the entire Graver basis
need not be used either.
• There are various ways of augmenting solutions like Greedy, Bisection methods, etc. They can be
applied based on the situation at hand.
Please refer to the notebook for more details about the implementation.

7.6 Surpassing Best in Class Classical Methods

Three areas where the D-Wave quantum computer is lacking are the number of qubits, their connectivity
with each other and the coupling precision.
• If the precision is increased, we can run QUBOs over a larger range of parameters
• If the number of qubits increases, we can demonstrate GAMA on problems in which we know the
state-of-art methods will not perform on par because of timing analysis for smaller problems.
• GAMA can also surpass classical methods in problems with complex convex and non-convex objective
functions.

7.7 Classical Implementation of GAMA

Suppose the matrix A has a special structure. In that case, we can construct the Graver Basis from first
principles and randomly generate many feasible solutions. This works very fast as compared to solvers like
Gurobi. Problems belonging to such classes include QAP, QSAP, and CBQP. This work [58] describes this
in more detail.
7-8 Lecture 7: GAMA - Graver Augmented Multiseed Algorithm

7.8 Summary and Conclusion

• In conventional methods of solving Integer Programs, we quadratize the objective and add constraints
to the objective with multipliers and optimize this entire mix.
• GAMA is a novel method where the objective function is used just as an oracle, and augmentation
sets are found using annealers.
• It is a Quantum-Classical hybrid algorithm consisting of two parts, finding the Graver Basis and
augmenting feasible solutions.
• The Graver Basis is created by solving for the Kernel using an annealer and classical post-processing
to create a (partial) Graver Basis set.
• Many initial feasible solutions are found using an annealer and are later augmented classically
• Since we have many starting points, it is likely that the partial Graver Basis algorithm will find optimum
values, even with complicated objective functions.
• There are some special cases where the Graver Bases can be computed classically, resulting in efficient
performance compared to commercial solvers.
• GAMA can thus be used in areas like Portfolio management, Cancer Genome Pathways, etc. [63],
where we have a complex objective function subject to linear constraints.
 Quantum Integer Programming September-December 2020

Lecture 8: Quantum Annealing

Lecturer: Prof. Sridhar Tayur, David E. Bernal Scribe: Shashank Kumar Ranu
Slides, Notebook Video, password: 2O8wq@Y%

8.1 Qubit

A qubit (or quantum bit) is the quantum-mechanical analog of a classical bit. In classical computing,
information is encoded in bits. A bit can have one of two values zero or one. In quantum computing,
information is encoded in qubits. A qubit is a two-level quantum system. Consider a system with 2 basis
states, |0i and |1i. The state of a qubit can be any superposition of these two basis states

α |0i + β |1i , |α|2 +|β|2 .

Using the constraints on α and β, a single qubit state can also be written as

θ θ
cos |0i + eiφ sin |1i .
2 2

The numbers 0 ≤ θ ≤ π and 0 ≤ φ ≤ 2π define a point on a unit three-dimensional sphere. This is called
a Bloch sphere. The Bloch sphere is a geometric representation of qubit states as points on the surface of
a unit sphere (see Fig. ‘8.1). Operations on single qubits that are commonly used in quantum information
processing can be neatly described within the Bloch sphere picture.

Figure 8.1: Geometrical representation of Bloch sphere.

8-1
8-2 Lecture 8: Quantum Annealing

8.2 Superconducting qubits

The ability to control individual quantum degrees of freedom and their interactions allows us to perform
quantum computation. Different technologies such as trapped atoms (ions or neutral), photons, single
electrons trapped in silicon heterostructures, and magnetic/electric moments of molecules have been used to
realise qubits.
Superconducting qubits are a leading approaches for realizing quantum logic elements and quantum co-
herent interactions with sufficiently high controllability and low noise. Their design can be understood from
the energy levels of an oscillator. The quantum harmonic oscillator shown in Fig. 8.3 is a resonant circuit
comprising of a capacitor and an inductor. For sufficiently low temperature and dissipation, the resulting
harmonic potential supports quantized energy levels. However, due to the equidistant level spacing, the
quantum harmonic oscillator by itself cannot be operated as a qubit. To remedy this situation, the circuit
potential is made anharmonic by introducing a nonlinear inductor – the Josephson junction. The imparted
anharmonicity leads to a non-equidistant spacing of the energy levels, enabling one to uniquely address each
transition, see Fig. 8.3. Typically, the two lowest levels are used to define a qubit, with |0i corresponding to
the ground state and |1i corresponding to the excited state. Large anharmonicity is generally favorable to
suppress unwanted excitations to higher levels. Different types of superconducting qubits exist in literature.

Figure 8.2: Energy spectrum of a quantum harmonic oscillator.

D-wave annealers use flux superconducting qubits, while Google and IBM use transmon qubits for designing
their quantum processors. The transmon qubit is a capacitively shunted variant of the Cooper pair box that
is largely insensitive to charge, resulting in improved reproducibility and coherence times. A flux qubit is
based on magnetic flux whose states correspond to clockwise and counter-clockwise currents flowing around
a loop interrupted by Josephson junctions. .
Lecture 8: Quantum Annealing 8-3

Figure 8.3: Tansmon and Flux qubits architecture.

8.3 The Quantum Adiabatic Algorithm (Quantum Annealing)

Time evolution of a quantum system is fairly straightforward for a time-independent Hamiltonian. As per
Schrodinger equation, the quantum state evolves as

|ψ(t)i = e−iHt |ψ(0)i , (8.105)

where H is the system Hamiltonian and |ψ(0)i is the initial state of the system. Hence, an eigenstate |Ei of
the Hamiltonian acquires a time-dependent overall phase.
Time evolution of a system having a time-dependent Hamiltonian is complex. However, if a Hamiltonian
changes slowly, the dynamics remain simple and can be understood using the quantum adiabatic theorem.
This theorem states that a quantum system which begins in a non-degenerate ground state of a time-
dependent Hamiltonian will remain in the instantaneous ground state provided the Hamiltonian changes
sufficiently slowly.

8.4 The Quantum Adiabatic Algorithm for Ising Machines

The quantum adiabatic theorem can be used to solve optimization problems. The first step is to map the
objective function of the problem to the Hamiltonian of the quantum Ising system, given as
X X
HP = Zi ⊗ Zj + hi Zi , (8.106)
Jij i

where Jij captures the nearest neighbour interaction strength, and Z is the Pauli Z matrix.
To apply the quantum adiabatic theorem: start from the minimum energy state of a problem with a known
solution. One such Hamiltonian is the transverse field Hamiltonian
X
HD = τ Xi . (8.107)
i

When acting upon n qubits, this Hamiltonian has an equal superposition of all the 2n possible binary
strings as its minimum eigenstate. In D-wave annealers, analog signals control the implementation of these
transverse and Ising models. We start with a transverse field acting on all the qubits and quench it slowly
to 0. In parallel, we increase the flux signal controlling the Ising model Hamiltonian. At some point in
time, both the signals co-exist, resulting in a complex dynamics. The analog signals are changed at a rate
which depends upon the optimization problem being implemented. The quantum adiabatic theorem does
not hold for the D-wave quantum annealer due to the system’s inherent noise. For example, mapping an
objective function to the Ising Hamiltonian is not perfect due to the intrinsic control errors. Benchmarking
and software mitigation of these errors is important to gain a proper understanding of the annealer output.
However, D-wave annealer is still an adiabatic inspired quantum computer.
8-4 Lecture 8: Quantum Annealing

8.5 Minor Embedding

Limited connectivity between the qubits is a restriction to employing the D-Wave quantum annealer for
real-world applications. Before solving an optimization problem, it is necessary to map a problem graph
onto a subgraph of the hardware graph. This process is called minor embedding. The problem graph is
defined as a graph in which the vertices and edges represent the logical variables and interactions between
them. The hardware graph is defined as a graph for which the vertices and edges represent the qubits and
interactions between them, respectively. For example, suppose we have a QUBO with all non-zero quadratic
coefficients. Implementing such a QUBO needs a fully connected annealer. Minor embedding procedure

Figure 8.4: An illustration of minor embedding.

consists of finding a set of connected subgraphs (logical qubits). These correspond to different colors in
Fig.??, of the original graph such that each logical qubit can be associated with a node in the original
graph. This association needs to be such that there is at least one edge between the qubits belonging to the
associated logical qubits for every two connected nodes. However, such an embedding increases the qubit
cost and changes the energy landscape. The problem of finding an optimal graph minor (i.e., minimizing the
number of required nodes) is itself NP-hard and is typically tackled with heuristic approaches. For many
graphs with a regular structure, an efficient embedding can be found systematically.

Figure 8.5: Change in energy landscape due to minor embedding.

8.6 Unembedding

Modeling the connectivity structure of a given problem instance thus necessitates the computation of a minor
embedding of the variables in the problem specification onto the logical qubits. One logical qubit consists of
several physical qubits” chained” together. After annealing, it is, however, not guaranteed that all chained
Lecture 8: Quantum Annealing 8-5

qubits get the same value (−1 or +1 for an Ising model, and 0 or 1 for a QUBO), and several approaches
exist to assign a final value to each logical qubit, a process called unembedding. The most straightforward
one is called majority vote. Suppose some logical qubit x is mapped onto a chain with m physical qubits.
After annealing, the values of m physical qubits are read, and x is assigned the most common value among
the m chained qubits (see Fig. 8.6).

Figure 8.6: Unembedding using majority voting.

8.7 Annealing Schedule Parameters

Two energy scales are relevant during quantum annealing, A(s) and B(s). A(s) represents the transverse, or
tunneling, energy. B(s) is the energy applied to the problem Hamiltonian. s is the ratio between the current
time during the anneal and the total annealing time and has a value between 0 and 1. Energy scales A(s)
and B(s) change during the quantum annealing process. In particular, a single, global, time-dependent bias
controls the trajectory of A and B. At any intermediate value of s, the ratio A(s)/B(s) is fixed. We can
choose the trajectory of one of A or B with time. The standard annealing schedule, s, produces a quadratic
growth in B(s). Typical values of A(s) and B(s) are shown in the Fig.8.7. There is a Reverse annealing
feature: this allows one to start (in a classical state with the strength of the problem Hamiltonian and
the driver Hamiltonians) at the end of the forward anneal, and to then increase the strength of the driver
Hamiltonian and reduce that of the problem Hamiltonian, following the schedule in reverse. One can pause,
and then repeat this forward-backward approach gain.
The pause feature, which allows one to pause the anneal, keeping the strengths of the driver Hamiltonian
and the problem Hamiltonian constant for extended periods of time before completing the default annealing
schedule (see Fig. 8.8. A pause can increase the performance by orders of magnitude when the pause occurs

Figure 8.7: Forward annealing.

8-6 Lecture 8: Quantum Annealing

Figure 8.8: Annealing with a pause.

within a well defined, relatively narrow region of the anneal but has little effect if placed outside that region.
 Quantum Integer Programming September-December 2020

Lecture 9: Quantum Annealing and QAOA

Lecturer: Prof. Sridhar Tayur, David E. Bernal, Dr. Davide Venturelli Scribe: Shrigyan Brahmachari
Slides, Notebook Video, password: ghE58r=E

• minor embedding
• implementing quantum annealing and comparisons
• benchmarking study-QAOA

9.1 Minor embedding

Quantum annealing works by constructing a physical situation where the energy of the system maps to the
objective function. Once the system to reaches a state of minimum energy, one obtains the minimum value
of the objective function.
A physical system’s mathematical abstraction is a graph ( the adjacency matrix can be used to get the
QUBO, as explained in the previous lectures).
However, the chimera graph ( subgraph of the chip) isn’t necessarily complete ( that is, not all pairs of
vertices are connected). In this case, we need to ”duplicate” certain variables into several qubits/nodes (
representing decision variables) This duplication step is non-trivial: we can either use a heuristic method or
solve a highly constrained problem. Here is an example to illustrate this:

zi ∈ 1, −1, (9.108)
min(z1 z2 + z2 z3 + z3 z1 ) (9.109)

needs to be mapped into a subgraph which ( at least locally) has a square lattice-like structure. Here is an
equivalent mapping using constraints.:

min(x1 x2 + x2 x3 + x3 x4 ) (9.110)
such that(x1 = x4 ) (9.111)

However constraints can’t be directly implemented. We need a QUBO form, so we use penalty functions:

min(x1 x2 + x2 x3 + x3 x4 + p(x1 − x4 )2 ), or (9.112)

2
min(x1 x2 + x2 x3 + x3 x4 − p(x1 x4 ) ), (9.113)

where the value of the chain strength p > 0 should be large enough so the constraint is satisfied. However
making p very large will dilute the effect of other coefficients in the objective function.

9-1
9-2 Lecture 9: Quantum Annealing and QAOA

9.2 Implementing quantum annealing and comparisons

Using classical Simulated Annealing with the default parameters increases the probability of finding the
optimal solution from (.511 ) to 0.25 and a feasible solution to almost 1. The ”infeasible” solutions are
heavily penalized, and random sampling is not an option.
Using Quantum Annealing with default parameters (annealing time, chain strength) results in the proba-
bility of finding the optimal solution of 5/10000 and feasible of 15/10000.
In Quantum Annealing, we analyzed two different factors, the chain strength and the annealing time. Our
primary concern is maximizing the probability of success (feasible and optimal).
From the plots below, it is observed that there is a massive fraction of chain-breaking for small values of
chain strength, i.e., the constraints aren’t satisfied. However, when we increase the chain strength (p), we
see that the number of optimal solution initially increases and then falls. Here, p can be decided by taking
both plots into account.
The following plots are for a shorter annealing time of 2 s, one-tenth of the previous set. While there
is a drastic decrease in the number of optimal solutions (the chain breaking seems unaffected), one should
consider whether the increased running time is a better strategy than to run the annealing several times,
somehow making sure different sets of solutions are found.
Lecture 9: Quantum Annealing and QAOA 9-3
9-4 Lecture 9: Quantum Annealing and QAOA
Lecture 9: Quantum Annealing and QAOA 9-5

9.3 How to do a benchmarking strategy

There are 3 kinds of quantum speedup that are considered:

• Provable Quantum Speedup (e.g. Grover)
• Strong Quantum Speedup (e.g. Shor)
• Quantum Speedup (potential, limited)
Note: it is possible to sometimes observe an apparent quantum speedup only because we may be dealing with
small sets of test cases and the algorithms perform suboptimally. This does not give us the true asymptotic
behavior.
Once we decide how to run, the benchmarking question is: What is the quality of a solution that I can
expect for a random new instance with a given confidence?
• Set up the quantum algorithm on the QPU with some initial parameters
• Run it several times and process the performance collecting the statistics of the distribution
• If performance is not acceptable, use the distribution to choose new parameters (might involve pro-
cessing) → Repeat 1-3 times until satisfaction
• Process final result and measure resources used (time, energy) → Repeat 1-4 for many benchmarking
instances and collect the distribution of performance.
• Compare (time, energy) against best classical method on available hardware
9-6 Lecture 9: Quantum Annealing and QAOA

9.4 Quantum Approximate Optimisation Algorithm

• A quantum computer is UNIVERSAL if its instruction set allows the implementation of any algorithm
allowed by quantum mechanics.
• The time-evolution of the Ising model in a transverse field (Quantum Annealing as implemented in
D-Wave) is NOT universal. However, the general AQC procedure is universal (need more complex
(Hp ) and (HD )). For quantum advantage in optimization heuristics, universality is not necessarily
required (the final state we are searching for is classical).
Why might you want a Universal Quantum Computer?
• Simulation of Quantum Systems
• Flexibility in implementation of multiple quantum algorithms (e.g., Grover/Shor)
• Exploit all the power of quantum mechanics
• Making sure that what you do is not classically simulatable efficiently
QAOA operates using the concept of quantum circuits. Essentially, it can be shown that any unitary
operation can be physically done to a qubit (or a set of entangled qubits) because unitary operators can be
shown to be an exponential of the time integral of a Hermitian, i.e., Hamiltonian operator. Using operators
that work on several subsystems (qubits), we can start with unentangled states into some desirable entangled
states. To understand the kind of circuits and therefore operators one needs to get the desirable quantum
entangled state, one uses linear algebra to decompose the transformation. Such decompositions are studied
in detail in introductory course on Quantum Computing. It can be shown that any unitary transformation
can be done by merely repeatedly using a recognized set of gates (X, Hadamard, CNOT) in an appropriate
sequence (see Neilsen and Chuang).
QAOA aims to implement adiabatic transitions coherent operations more flexibly than AQC (digitally).
For infinite circuit this is at least as powerful as AQC. For finite circuit its power is unknown in general.
The following is an illustration of a simple example:

We start with some unentangled qubits; after the circuit transformation, the entangled state is such that
when you measure the probability of a bit string, it depends both on γ and β in a non-linear way.
The primary difference to AQC is that the state that maps to the minimum energy is obtained as part of a
superposition with other states. So there is a certain probability of realizing it by observation. Calculating
this probability becomes very difficult.
 R-1

References
[1] Michele Conforti, Gérard Cornuéjols, Giacomo Zambelli, et al. Integer programming, volume 271.
Springer, 2014.

[2] Robert Weismantel. Test sets of integer programs. Mathematical Methods of Operations Research, 47
(1):1–37, 1998.

[3] Bernd Sturmfels. Grobner bases and convex polytopes, volume 8. American Mathematical Soc., 1996.

[4] Sridhar R Tayur, Rekha R Thomas, and NR Natraj. An algebraic geometry algorithm for scheduling
in presence of setups and correlated demands. Mathematical Programming, 69(1-3):369–401, 1995.

[5] Jesús A De Loera, Raymond Hemmecke, and Matthias Köppe. Algebraic and geometric ideas in the
theory of discrete optimization. SIAM, 2012.

[6] Dimitris Bertsimas, Georgia Perakis, and Sridhar Tayur. A new algebraic geometry algorithm for integer
programming. Management Science, 46(7):999–1008, 2000.

[7] Pasqualina Conti and Carlo Traverso. Buchberger algorithm and integer programming. In International
Symposium on Applied Algebra, Algebraic Algorithms, and Error-Correcting Codes, pages 130–139.
Springer, 1991.

[8] Serkan Hoşten and Bernd Sturmfels. Grin: An implementation of gröbner bases for integer programming.
In International Conference on Integer Programming and Combinatorial Optimization, pages 267–276.
Springer, 1995.

[9] Raymond Hemmecke, Shmuel Onn, and Robert Weismantel. A polynomial oracle-time algorithm for
convex integer minimization. Mathematical Programming, 126(1):97–117, 2011.

[10] Stephen G Brush. History of the Lenz-Ising model. Reviews of Modern Physics, 39(4):883, 1967.

[11] David Sherrington and Scott Kirkpatrick. Solvable model of a spin-glass. Phys. Rev. Lett., 35:1792–
1796, Dec 1975. doi: 10.1103/PhysRevLett.35.1792. URL [Link]
PhysRevLett.35.1792.

[12] P. Ray, B. K. Chakrabarti, and Arunava Chakrabarti. Sherrington-kirkpatrick model in a transverse

field: Absence of replica symmetry breaking due to quantum fluctuations. Phys. Rev. B, 39:11828–11832,
Jun 1989. doi: 10.1103/PhysRevB.39.11828. URL [Link]
39.11828.

[13] Scott Kirkpatrick, C Daniel Gelatt, and Mario P Vecchi. Optimization by simulated annealing. Science,
220(4598):671–680, 1983.

[14] Christos Koulamas, SR Antony, and R Jaen. A survey of simulated annealing applications to operations
research problems. Omega, 22(1):41–56, 1994.

[15] Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward
Teller. Equation of state calculations by fast computing machines. The Journal of Chemical Physics,
21(6):1087–1092, 1953.

[16] Alfred B Bortz, Malvin H Kalos, and Joel L Lebowitz. A new algorithm for monte carlo simulation of
ising spin systems. Journal of Computational Physics, 17(1):10–18, 1975.

[17] Matthias Troyer and Uwe-Jens Wiese. Computational complexity and fundamental limitations to
fermionic quantum monte carlo simulations. Physical Review Letters, 94(17):170201, 2005.
R-2

[18] A Peter Young, Sergey Knysh, and Vadim N Smelyanskiy. Size dependence of the minimum excitation
gap in the quantum adiabatic algorithm. Physical Review Letters, 101(17):170503, 2008.

[19] Iain Dunning, Swati Gupta, and John Silberholz. What works best when? a systematic evaluation of
heuristics for max-cut and qubo. INFORMS Journal on Computing, 30(3):608–624, 2018.

[20] Carleton Coffrin, Harsha Nagarajan, and Russell Bent. Evaluating ising processing units with inte-
ger programming. In International Conference on Integration of Constraint Programming, Artificial
Intelligence, and Operations Research, pages 163–181. Springer, 2019.

[21] Andrew Lucas. Ising formulations of many np problems. Frontiers in Physics, 2:5, 2014.

[22] Hedayat Alghassi, Raouf Dridi, and Sridhar Tayur. Graver bases via quantum annealing with application
to non-linear integer programs. arXiv preprint arXiv:1902.04215, 2019.

[23] Hedayat Alghassi, Raouf Dridi, A Gordon Robertson, and Sridhar Tayur. Quantum and quantum-
inspired methods for de novo discovery of altered cancer pathways. bioRxiv, page 845719, 2019.

[24] Hedayat Alghassi, Raouf Dridi, and Sridhar Tayur. Gama: A novel algorithm for non-convex integer
programs. arXiv preprint arXiv:1907.10930, 2019.

[25] Catherine C McGeoch. Theory versus practice in annealing-based quantum computing. Theoretical
Computer Science, 816:169–183, 2020.

[26] Tameem Albash and Daniel A Lidar. Adiabatic quantum computation. Reviews of Modern Physics, 90
(1):015002, 2018.

[27] Arnab Das and Bikas K. Chakrabarti. Colloquium: Quantum annealing and analog quantum com-
putation. Rev. Mod. Phys., 80:1061–1081, Sep 2008. doi: 10.1103/RevModPhys.80.1061. URL
[Link]

[28] Giuseppe E Santoro and Erio Tosatti. Optimization using quantum mechanics: quantum annealing
through adiabatic evolution. Journal of Physics A: Mathematical and General, 39(36):R393, 2006.

[29] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, Joshua Lapan, Andrew Lundgren, and Daniel Preda.
A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem.
Science, 292(5516):472–475, 2001.

[30] Tadashi Kadowaki and Hidetoshi Nishimori. Quantum annealing in the transverse ising model. Physical
Review E, 58(5):5355, 1998.

[31] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J.

Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky,
N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin,
J. Wang, B. Wilson, and G. Rose. Quantum annealing with manufactured spins. Nature, 473(7346):
194–198, May 2011. ISSN 0028-0836. doi: 10.1038/nature10012. URL [Link]
nature/journal/v473/n7346/full/[Link].

[32] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann. A quantum approximate optimization algorithm.
arXiv preprint arXiv:1411.4028, 2014.

[33] Stuart Hadfield, Zhihui Wang, Eleanor G Rieffel, Bryan O’Gorman, Davide Venturelli, and Rupak
Biswas. Quantum approximate optimization with hard and soft constraints. In Proceedings of the
Second International Workshop on Post Moores Era Supercomputing, pages 15–21, 2017.
 R-3

[34] Stuart Hadfield, Zhihui Wang, Bryan O’Gorman, Eleanor G Rieffel, Davide Venturelli, and Rupak
Biswas. From the quantum approximate optimization algorithm to a quantum alternating operator
ansatz. Algorithms, 12(2):34, 2019.

[35] A. Yavorsky, L. A. Markovich, E. A. Polyakov, and A. N. Rubtsov. Highly parallel algorithm for the
ising ground state searching problem, 2019.

[36] Chase Cook, Hengyang Zhao, Takashi Sato, Masayuki Hiromoto, and Sheldon X-D Tan. Gpu-based ising
computing for solving max-cut combinatorial optimization problems. Integration, 69:335–344, 2019.

[37] Joshua Romero, Mauro Bisson, Massimiliano Fatica, and Massimo Bernaschi. A performance study of
the 2d ising model on gpus. arXiv preprint arXiv:1906.06297, 2019.

[38] Kun Yang, Yi-Fan Chen, Georgios Roumpos, Chris Colby, and John Anderson. High performance monte
carlo simulation of ising model on tpu clusters. In Proceedings of the International Conference for High
Performance Computing, Networking, Storage and Analysis, pages 1–15, 2019.

[39] Masanao Yamaoka, Chihiro Yoshimura, Masato Hayashi, Takuya Okuyama, Hidetaka Aoki, and Hi-
royuki Mizuno. A 20k-spin ising chip to solve combinatorial optimization problems with cmos annealing.
IEEE Journal of Solid-State Circuits, 51(1):303–309, 2015.

[40] Maliheh Aramon, Gili Rosenberg, Elisabetta Valiante, Toshiyuki Miyazawa, Hirotaka Tamura, and
Helmut Katzgrabeer. Physics-inspired optimization for quadratic unconstrained problems using a digital
annealer. Frontiers in Physics, 7:48, 2019.

[41] Jeffrey Chou, Suraj Bramhavar, Siddhartha Ghosh, and William Herzog. Analog coupled oscillator
based weighted ising machine. Scientific Reports, 9(1):1–10, 2019.

[42] Tianshi Wang and Jaijeet Roychowdhury. Oim: Oscillator-based ising machines for solving combinato-
rial optimisation problems. In International Conference on Unconventional Computation and Natural
Computation, pages 232–256. Springer, 2019.

[43] Charles Roques-Carmes, Yichen Shen, Cristian Zanoci, Mihika Prabhu, Fadi Atieh, Li Jing, Tena
Dubček, Chenkai Mao, Miles R Johnson, Vladimir Čeperić, et al. Heuristic recurrent algorithms for
photonic ising machines. Nature Communications, 11(1):1–8, 2020.

[44] Takahiro Inagaki, Yoshitaka Haribara, Koji Igarashi, Tomohiro Sonobe, Shuhei Tamate, Toshimori
Honjo, Alireza Marandi, Peter L McMahon, Takeshi Umeki, Koji Enbutsu, et al. A coherent ising
machine for 2000-node optimization problems. Science, 354(6312):603–606, 2016.

[45] Andrew D King, William Bernoudy, James King, Andrew J Berkley, and Trevor Lanting. Emulating
the coherent ising machine with a mean-field algorithm. arXiv preprint arXiv:1806.08422, 2018.

[46] Ryan Hamerly, Takahiro Inagaki, Peter L McMahon, Davide Venturelli, Alireza Marandi, Tatsuhiro
Onodera, Edwin Ng, Carsten Langrock, Kensuke Inaba, Toshimori Honjo, et al. Experimental inves-
tigation of performance differences between coherent ising machines and a quantum annealer. Science
advances, 5(5):eaau0823, 2019.

[47] Egor S Tiunov, Alexander E Ulanov, and AI Lvovsky. Annealing by simulating the coherent ising
machine. Optics express, 27(7):10288–10295, 2019.

[48] Peter L McMahon, Alireza Marandi, Yoshitaka Haribara, Ryan Hamerly, Carsten Langrock, Shuhei
Tamate, Takahiro Inagaki, Hiroki Takesue, Shoko Utsunomiya, Kazuyuki Aihara, et al. A fully pro-
grammable 100-spin coherent ising machine with all-to-all connections. Science, 354(6312):614–617,
2016.
R-4

[49] David E Bernal, Kyle EC Booth, Raouf Dridi, Hedayat Alghassi, Sridhar Tayur, and Davide Ven-
turelli. Integer programming techniques for minor-embedding in quantum annealers. arXiv preprint
arXiv:1912.08314, 2019.
[50] Raouf Dridi, Hedayat Alghassi, and Sridhar Tayur. A novel algebraic geometry compiling framework
for adiabatic quantum computations. arXiv preprint arXiv:1810.01440, 2018.
[51] Raouf Dridi, Hedayat Alghassi, and Sridhar Tayur. Knuth-bendix completion algorithm and shuffle
algebras for compiling nisq circuits. arXiv preprint arXiv:1905.00129, 2019.
[52] Raouf Dridi, Hedayat Alghassi, and Sridhar Tayur. Minimizing polynomial functions on quantum
computers. arXiv preprint arXiv:1903.08270, 2019.
[53] Raouf Dridi and Hedayat Alghassi. Prime factorization using quantum annealing and computational
algebraic geometry. Scientific reports, 7:43048, 2017.
[54] Michael A Nielsen and Isaac Chuang. Quantum computation and quantum information, 2002.
[55] Jesús Gago-Vargas, I Hartillo, Justo Puerto, and JM Ucha. An improved test set approach to nonlinear
integer problems with applications to engineering design. Computational Optimization and Applications,
62(2):565–588, 2015.

[56] Jack E. Graver. On the foundations of linear and integer linear programming i, 1975. URL https:
//[Link]/article/10.1007/BF01681344.
[57] Shmuel Onn. Nonlinear discrete optimization, 2010. URL [Link]
Nachdiplom/.

[58] Hedayat Alghassi, Raouf Dridi, and Sridhar Tayur. Gama: A novel algorithm for non-convex integer
programs, 2019. URL [Link]
[59] Tepper School of Business at Carnegie Mellon University. Quantum integer programming, 2020. URL
[Link]

[60] Matthias Köppe & Robert Weismantel Utz-Uwe Haus. The integral basis method for integer program-
ming, 2001. URL [Link]
[61] Martin Koutecký Kateřina Altmanová, Dušan Knop. Evaluating and tuning n-fold integer programming,
2019. URL [Link]
[62] Hedayat Alghassi, Raouf Dridi, and Sridhar Tayur. Graver bases via quantum annealing with application
to non-linear integer programs, 2019.
[63] Hedayat Alghassi, Raouf Dridi, A. Gordon Robertson, and Sridhar Tayur. Quantum and quantum-
inspired methods for de novo discovery of altered cancer pathways. bioRxiv, 2019. doi: 10.1101/845719.
URL [Link]