Prime factorization using quantum annealing and
computational algebraic geometry
arXiv:1604.05796v2 [quant-ph] 16 Jun 2016

Raouf Dridi∗ and Hedayat Alghassi†
1QB Information Technologies (1QBit)
458–550 Burrard Street
Vancouver, British Columbia, Canada V6C 2B5
June 17, 2016

Abstract
We investigate prime factorization from two perspectives: quantum annealing and
computational algebraic geometry, specifically Gröbner bases. We present a novel
autonomous algorithm which combines the two approaches and leads to the factorization of all bi-primes up to just over 200 000, the largest number factored to date
using a quantum processor. We also explain how Gröbner bases can be used to
reduce the degree of Hamiltonians.

1

Introduction

Prime factorization is at the heart of secure data transmission because it is widely believed
to be NP-complete. In the prime factorization problem, for a large bi-prime M , the task
is to find the two prime factors p and q such that M = pq. In secure data transmission,
the message to be transmitted is encrypted using a public key which is, essentially, a large
bi-prime that can only be decrypted using its prime factors, which are kept in a private
key. Prime factorization also connects to many branches of mathematics; two branches
relevant to us are computational algebraic geometry and quantum annealing.
To leverage the problem of finding primes p and q into the realm of computational algebraic
geometry, it suffices to transform it into an
S. This is done
P algebraic system of equations
P
using the binary representation p = 1 + i=1..sp 2i Pi and q = 1 + i=1..sq 2i Qi , which is
plugged into M = pq and expanded into a system of polynomial equations. The reader is
invited to read the sections Methods 4.1 and 4.2 for the details of this construction. The
system S is given by this initial system of equations in addition to the auxiliary equations
∗
†

raouf.dridi@1qbit.com
hedayat.alghassi@1qbit.com

1

expressing the binary nature of the variables Pi and Qi , carry-on, and connective variables.
The two primes p and q are then given by the unique zero of S. In theory, we can solve
the system S using Gröbner bases; however, in practice, this alone does not work, since
Gröbner basis computation (Buchberger’s algorithm) is exponential.
The connection to quantum annealing can also be easily described. Indeed, finding p
and q can be formulated into an unconstrained binary optimization problem (P), where
the cost function f is the sum of the squares of polynomials in S. The unique zero of S
now sits on the unique global minimum of (P) (which has minimum energy equal to zero).
There are, however, a few non-trivial requirements we need to deal with before solving
the cost function using quantum annealing. These requirements concern the nature of
cost functions that quantum annealers can handle. In particular, we would like the cost
function of (P) to be a positive quadratic polynomial. We also require that the coefficients
of the cost function be rather uniform and match the hardware-imposed dynamic range.
In the present paper, we suggest looking into the problem through both lenses, and
demonstrate that indeed this approach gives better results. In our scheme, we will be using quantum annealing to solve (P), but at the same time we will be using Gröbner bases
to help us reduce the cost function f into a positive quadratic polynomial f + with desired
values for the coefficients. We will be also using Gröbner bases at the important step of
pre-processing f + before finally passing it to the quantum annealer. This pre-processing
significantly reduces the size of the problem. The result of this combined approach is an
algorithm with which we have been able to factorize all bi-primes up to 2 × 105 using the
D-Wave 2X processor. The algorithm is autonomous in the sense that no a priori knowledge, or manual or ad hoc pre-processing, is involved. We refer the interested reader to
Supplementary materials 5 for a brief description of the D-Wave 2X processor, along with
some statistics for several of the highest numbers that we embedded and solved. More
detail about the processor architecture can be found in [JAG+ 11]. Another important
reference is the work of S. Boixo et al. in [BAS+ 13], which presents experimental evidence that the scalable D-Wave processor implements quantum annealing with surprising
robustness against noise and imperfections. Evidence that, during a critical portion of
quantum annealing, the qubits become entangled and entanglement persists even as the
system reaches equilibrium is presented in [LPS+ 14].
Relevant to us also is the work in [PS01], which uses algebraic geometry to solve optimization problems (though not specifically factorization; see Methods 4.4 for an adaptation
to factorization). Therein, Gröbner bases are used to compute standard monomials and
transform the given optimization problem into an eigenvalue computation. Gröbner basis computation is the main step in this approach, which makes it inefficient, given the
high complexity of Gröbner basis computation. In contrast to that work, we ultimately
solve the optimization problem using a quantum annealing processor and pre-process and
adjust the problem with algebraic tools, that is, we reduce the size of the cost function
and adjust the range of its parameters. However, we share that work’s point of view of
using real algebraic geometry, and our work is the first to introduce algebraic geometry,
and Gröbner bases in particular, to solve quantum annealing-related problems. We think
that this is a fertile direction for both practical and theoretical endeavours.
2

Mapping the factorization problem into a degree-4 unconstrained binary optimization
problem is first discussed in [Bur02]. There, the author proposes solving the problem using
a continuous optimization method he calls curvature inversion descent. Another related
work is the quantum annealing factorization algorithm proposed in [SS10]. We will discuss
it in the next section and improve upon it in two ways. The first involves the addition
of the pre-processing stage using Gröbner bases of the cost function. This dramatically
reduces the number of variables therein. The second way concerns the reduction of the
initial cost function, for which we propose a general Gröbner basis scheme that precisely
answers the various requirements of the cost function. In Results 2.2, we present our
algorithm (the column algorithm) which outperforms this improved algorithm (i.e., the
cell algorithm). Using a reduction proposed in [SS10] and ad-hoc simplifications and
tricks, the paper [XZL+ 12] reports the factorization of bi-prime 143 on a liquid-crystal
NMR quantum processor, which until now was the largest reported bi-prime number
factored in any quantum computation realization.
This review is far from complete without mentioning Shor’s algorithm [Sho97] and Kitaev’s
phase estimation [Kit95], which, respectively, solve the factorization problem and the
abelian hidden subgroup problem in polynomial time, both for the gate model paradigm.
The largest number factored using a physical realization of Shor’s algorithm is 15 [MNM+ 16];
see [SSV13] also for a discussion about oversimplification in the previous realizations. Finally, in [Rau13], it has been proved that contextuality is needed for any speed-up in a
measurement-based quantum computation factorization algorithm.

2

Results

The binary multiplication of the two primes p and q can be expanded in two ways: cellbased and column-based procedures (see Methods 4.1 and 4.2 ). Each procedure leads to
a different unconstrained binary optimization problem. The cell-based procedure creates
the unconstrained binary quadratic programming problem
P
(
minZ2 ij Hij 2 ,
(2.1)
(P1 )
with Hij := Qi Pj + Si,j + Zi,j − Si+1,j−1 − 2 Zi,j+1 ,
and the column-based procedure results in the problem
P

minZ2 1≤i≤(sp +sq +1) Hi 2 ,


sq +1+i−m
sq
(P2 )
i
P
P
P i j−i

 with Hi :=
Qj Pi−j +
Zj,i − mi −
2 Zi,i+j .
j=0

j=1

(2.2)

j=1

The two problems (P1 ) and (P2 ) are equivalent. Their cost functions are not in quadratic
form, and thus must be reduced before being solved using a quantum annealer. The
reduction procedure is not a trivial task. In this paper we define, for both scenarios:
1) a reduced quadratic positive cost function and 2) a pre-processing procedure. Thus,
we present two different quantum annealing-based prime factorization algorithms. The
3

first algorithm’s decomposition method (i.e., the cell procedure, Methods 4.1 ) has been
addressed in [SS10], without pre-processing and without the use of Gröbner bases in the
reduction step. Here, we discuss it from the Gröbner bases framework and add the important step of pre-processing. The second algorithm, however, is novel in transformation of
its quartic terms to quadratic, outperforming the first algorithm due to its having fewer
variables.
We write R[x1 , . . . , xn ] for the ring of polynomials in x1 , . . . , xn with real coefficients and
V(f ) for the affine variety defined by the polynomial f ∈ R[x1 , . . . , xn ], that is, the set of
zeros of the equation f = 0. Since we are interested only in the binary zeros (i.e., xi ∈ Z2 ),
we need to add the binarization polynomials xi (xi − 1), where i = 1, ..., n, to f and obtain
the system S = {f, xi (xi − 1), i = 1, ..., n}. The system S generates an ideal I by taking
all linear combinations over R[x1 , . . . , xn ] of all polynomials in S; we have V(S) = V(I).
The ideal I reveals the hidden polynomials which are the consequence of the generating
polynomials in S. To
√ hidden polynomials is given by rthe so√ be precise, the set of all
called radical ideal √I, which is defined by I = {g ∈ R[x1 , . . . , xn ]| ∃r ∈ N : g ∈ I}.
In practice, the ideal I is infinite, so we represent such an √
ideal using a Gröbner basis B
which one might take to be a triangularization of the ideal I. In fact, the computation
of Gröbner√
bases generalizes Gaussian elimination
in linear systems. We also have V(S) =
√
V(I) = V( I) = V(B) and I(V(I)) = I. A brief review of Gröbner bases is given in
Methods 4.3.

2.1

The cell algorithm

Suppose we would like to define the variety V(I) by the set of global minima of an
unconstrained optimization problem minZn2 (f + ), where f + is a quadratic polynomial.
For instance, we would like f + to behave like f 2 . Ideally, we want f + to remain in
R[x1 , . . . , xn ] (i.e., not in a larger ring), which implies that no slack variables will be
added. We also want f + to satisfy the following requirements:
√
(i) f + vanishes on V (I) or, equivalently, f + ∈ I.
(ii) f + > 0 outside V (I), that is, f + > 0 over Zn2 − V (I).
(iii) Coefficients of the polynomial f + are adjusted with respect to the dynamic range
allowed by the quantum processor.
Let B be a Gröbner basis for I. We can then go ahead and define
X
f+ =
at t,
t∈B| deg(t)≤2

where the real
√ coefficients ai are subject to the requirements above; note that we already
+
have f ∈ I and thus the first requirement (i) is satisfied.
Let us apply this procedure to the optimization problem (P1 ) above. There, f = Hij
and the ring of polynomials is R[Pj , Qi , Si,j , Si+1,j−1 , Zi,j , Zi,j+1 ]. We obtain the following
Gröbner basis:
4















































t1

:= Qi Pj + Si,j + Zi,j − Si+1,j−1 − 2 Zi,j+1 ,

t2

:= (−Zi,j+1 + Zi,j ) Si+1,j−1 + (Zi,j+1 − 1) Zi,j ,

t3

:= (−Zi,j+1 + Zi,j ) Si,j + Zi,j+1 − Zi,j+1 Zi,j ,

t4

:= (Si+1,j−1 + Zi,j+1 − 1) Si,j − Si+1,j−1 Zi,j+1 ,

t5

:= (−Si+1,j−1 − 2 Zi,j+1 + Zi,j + Si,j ) Qi − Si,j − Zi,j + Si+1,j−1 + 2 Zi,j+1 ,

t6

:= (−Si+1,j−1 − 2 Zi,j+1 + Zi,j + Si,j ) Pj − Si,j − Zi,j + Si+1,j−1 + 2 Zi,j+1 ,

t7

:= (−Zi,j+1 + Zi,j+1 Zi,j ) Qi + Zi,j+1 − Zi,j+1 Zi,j ,

t8

:= −Si+1,j−1 Zi,j+1 + Si+1,j−1 Qi Zi,j+1 ,

t9

:= (−Zi,j+1 + Zi,j+1 Zi,j ) Pj + Zi,j+1 − Zi,j+1 Zi,j ,

(2.3)

t10 := −Si+1,j−1 Zi,j+1 + Si+1,j−1 Pj Zi,j+1 .

We have used the lexicographic order plex(Pj , Qi , Si,j , Si+1,j−1 , Zi,j , Zi,j+1 ); see Methods
4.3 for definitions. Note that t1 = Hij . We define
X
X
Hij+ =
at t, that is, Hij+ =
ai ti ,
1≤i≤6

t∈B| deg(t)≤2

where the real coefficients ai are to be found. We need to constrain the coefficients ai
with the other requirements. The second requirement (ii), which translates into a set
of inequalities on the unknown coefficients ai , can be obtained through a brute force
evaluation of Hij and Hij+ over the 26 points of Z62 . The outcome of this evaluation is a set
of inequalities expressing the second requirement (ii) (see Supplementary materials 5 ).
The last requirement (iii) can be expressed in different ways. We can, for instance, require that the absolute values of the coefficients of Hij+ , with respect to the variables
Pj , Qi , Si,j , Si+1,j−1 , Zi,j , and Zi,j+1 , be within [1 − , 1 + ]. This, together with the set of
inequalities from the second requirement, define a continuous optimization problem and
can be easily solved. Another option is to minimize the distance between the coefficients to
one coefficient. The different choices of the objective function and the solution of the corresponding continuous optimization problem are presented in Supplementary materials 5 .
Having determined the quadratic polynomial Hij+ ∈ R satisfyies the important requirements (i, ii, and iii) above, we can now phrase our problem
(P1 ) as the equivalent quadratic
P
unconstrained binary optimization problem minZ2 ij Hij+ . Notice that this reduction is
performed only once for all cases; it need not to be redone for different bi-primes M .
Before passing the problem to the quantum annealer, we use Gröbner bases P
to reduce the

size of the problem. In fact, what we pass to the quantum annealer is H = NFB Hij+ ,
where NF is the normal form and B is now the Gröbner basis cutoff, which we discuss in
the next section. The largest bi-prime number that we embedded and solved successfully
using the cell algorithm is ∼35 000. The following table presents some bi-prime numbers M that we factored using the cell algorithm, the number of variables using both the
customized reduction CustR and the window-based GB reduction, the overall reduction
5

percentage R%, and the embedding and solving status inside the D-Wave 2X processor
Embed.
Cell algorithm
M
31861
34889
46961
150419

2.2

p×q
CustR GB R% Embed
√
211 × 151
111
95 14
√
251 × 139
111
95 14
311 × 151
125
109 13
×
431 × 349
143
125 12
×

The column algorithm (factoring up to 200 000)

The total number of variables in the cost function of the previous method is 2sp sq , before
any pre-processing. Here we present the column-based algorithm where the number of
variables (before pre-processing) is bounded by 1 + sp sq + log2 (sp )(sp + sq ). Recall that
here we are phrasing the factorization problem M = pq as
X
(P2 ) : minP1 ,...,Psp ,Q1 ,...,Qsq ,Z12 ,Z23 ,Z24 ,...∈Z2
Hi 2 ,
i

where Hi , for 1 ≤ i ≤ sp , is
Hi =

sq
X
j=0

Li
i
X
X
Qj Pi−j +
Zj,i − mi −
2j−i Zi,i+j
j=1

(Q0 = P0 = m0 = 1, Li = sq + 1 + i − mi ).

j=1

The cost function is of degree 4 and, in order to use quantum annealing, it must be
replaced with a positive quadratic polynomial with the same global minimum. The idea
is to replace the quadratic terms Qj Pi−j inside the different Hi with new binary variables
Wi−j,j , and add the penalty (Qj Pi−j − Wi−j,j )+ to the cost function (now written in terms
of the variables Wi−j,j ). To find (Qj Pi−j − Wi−j,j )+ , we run Gröbner bases computation
on the system

Qj Pi−j − W,




 Qj 2 − Qj ,
(2.4)
2

P
−
P
,

i−j
i−j



Wi−j,j 2 − Wi−j,j .
Following the same steps as in the previous section, we get
(Qj Pi−j − Wi−j,j )+ = a(Qj Wi−j,j − Wi−j,j ) + b(Pi−j Wi−j,j − Wi−j,j ) + c(Pi−j Qj − Wi−j,j ),
with a, b, c ∈ R such that −a − b − c > 0, −b − c > 0, −a − c > 0, c > 0 (e.g., c = 1, a =
b = −2). The new cost function is now
X
X
H=
Hi (W )2 +
(Qj Pi−j − Wi−j,j )+ .
i

ij

6

We can obtain a better Hamiltonian by pre-processing the problem before applying the
W transformation. Indeed, let us first fix a positive integer cutoff ≤ (sp + sq + 1) and let
B ⊂ R [P1 , . . . , Psp , Q1 , . . . , Qsq , Z12 , Z23 , Z24 . . .] be a Gröbner basis of the set of polynomials
{Hi }i=1...cutoff ∪ {Pi (Pi − 1), Qi (Qi − 1), Zij (Zij − 1)}i,j .
In practice, the cutoff is determined by the size of the maximum subsystem of polynomials Hi on which one can run a Gröbner basis computation; it is defined by the hardware.
We also define a cutoff on the other tail of {Hi }, that is, we consider {Hi }i=2ndcutoff...(sp +sq +1) .
Notice that here we are working on the original Hi rather than the new Hi (W ). This is because we would like to perform the replacement Qj Pi−j → Wi−j,j after the pre-processing
(some of the quadratic terms might be simplified by this pre-processing). Precisely, what
we pass to the quantum annealer is the quadratic positive polynomial
X
2 X
H=
NFWi−j,j −LT(NFBc (Qj Pi−j )) (NFBc (Hi )) +
(Wi−j,j − LT (NFBc (Qj Pi−j )))+ .
ij

(2.5)
Here LT stands for the leading term with respect to the graded reverse lexicographic
order. The second summation is over all i and j such that LT (NFB (Qj Pi−j )) is still
quadratic. The outer normal form in the first summation refers to the replacement
LT (NFB (Qj Pi−j )) → Wi−j,j , which is again performed only if LT (NFB (Qj Pi−j )) is still
quadratic.
The columns of the following table present: a sample of bi-prime numbers and their prime
factors, the number of variables using each of a naı̈ve polynomial-to-quadratic transformation P2Q, our novel polynomial-to-quadratic transformation CustR, and our windowbased reduction GB after applying pre-processing. The overall reduction percentage R%
and the embedding and solving status in the D-Wave 2X processor Embed are also shown.
The adjacency matrix of the corresponding positive quadratic polynomial graph H and
its embedded pattern inside the Chimera graph of the D-Wave 2X processor for one of
the bi-primes are also depicted (see Figure 1). Details pertaining to use of the hardware
can be found in Supplementary materials 5.3.
Column Algorithm
M
150419
151117
174541
200099
223357

3

p×q
P 2Q CustR GB R% Embed
√
431 × 349 116
86
73 37
√
433 × 349 117
88
72 38
√
347 × 503 117
86
72 38
√
499 × 401 115
89
75 35
557 × 401 125
96
80 36
×

Discussion

In this work, factorization is connected to quantum annealing through binarization of the
long multiplication. The algorithm is autonomous in the sense that no a priori knowledge,
or manual or ad hoc pre-processing, is involved. We have attained the largest bi-prime
7

Figure 1: The column algorithm: the adjacency matrix pattern (left) and embedding into
the the D-Wave 2X quantum processor (right) of the quadratic binary polynomial for
M = 200 099.
factored to date using a quantum processor, though more-subtle connections might exist.
A future direction that this research can take is to connect factorization (as an instance of
the abelian hidden subgroup problem), through Galois correspondence, to covering spaces
and thus to covering graphs and potentially to quantum annealing. We believe that morerewarding progress can be made through the investigation of such a connection.

4
4.1

Methods
Column factoring procedure

Here we discuss the two single-bit multiplication methods of the two primes p and q.
The first method generates a Hamiltonian for each of the columns of the long multiplication expansion, while the second method generates a Hamiltonian for each of the
multiplying cells in the long multiplication expansion. The column factoring procedure of
p = 2sp Psp + 2sp−1 Psp−1 + ... + 2P1 + 1 and q = 2sq Qsq + 2sq−1 Qsq−1 ... + 2Q1 + 1 is depicted
in the following table:

8

···

Psp

msp+sq+1

...
···
···
···
···

Qsq−1Psp
Qsq Psp Qsq Psp−1
msp+sq msp+sq−1

Pi

···
Pi
···
Q1Pi−1
···
Q2Pi−2
..
···
· · · Qsq−1Pi−sq+1
···
Qsq Pi−sq
···
mi

Psp
Q1Psp−1
Q2Psp−2
..
Qsq−1Psp−sq+1
Qsq Psp−sq
msp

···
···
···
···
···
···
···
···

Psq
Qsq
Psq
Q1Psq−1
Q2Psq−2
..
Qsq−1P1
Qsq
msq

· · · P2 P1 P0 = 1
· · · Q2 Q1 Q0 = 1
· · · P2 P1
1
· · · Q1 P 1 Q1
· · · Q2
···
...
···

m2

m1 m0 = 1

The equation for an arbitrary column (i) can be written as the sum of the column’s multiplication terms (above) plus all previously generated carry-on terms from lower significant
columns (j < i). This sum is equal to the column’s bi-prime term mi plus the carry-ons
generated from higher significant columns. The polynomial equation for the i-th column
is
sq
P

Qj Pi−j +

j=0

i
P

Zj,i = mi +

j=1

Li
P

2j−i Zi,i+j (Q0 = P0 = m0 = 1) .

j=1

The above equation is used as the main column procedure’s equation Hi . The Hamiltonian
generation and reduction is discussed in detail in Results 2.2.

4.2

Cell bi-prime factoring procedure

In the cell multiplication procedure, the ultimate goal is to break each of the column equations discussed above into multiple smaller equations so that each equation contains only
one quadratic term. This not only simplifies the generation of quadratic Hamiltonians,
but also generates Hamiltonians with more-uniform quadratic coefficients in comparison
to the column procedure. The following table depicts the structure of the cell procedure:
.

...

Psp

· · · Pi+j · · ·

0
Q1 P n
Z1,n

Pn
Q1Pn−1
Z1,n−1

..

Z1,n+1
Q2Pn
Z2,n

S25
Q2Pn−1
Z2,n−1

S24
Q2P4
Z2,n−2

..

..

..

..
Zi,j+1

msq+sp+1

Zsq−1,sp+1
Qsq Psp
Zsq,sp
msq+sp

Zsq−2,sp+1
Qsq−1Psp
Zsq−1,sp

Ssq−1,sp−1
Qsq−1Psp−1
Zsq−1,sp−1

Ssq−1,sp−2
Qsq−1Psp−2
Zsq−1,sp−2

Ssq−1,sp−3
Qsq−1Psp−3
Zsq−1,sp−3

Ssq,sp−1
Qsq Psp−1
Zsq,sp−1
msq+sp−1

Ssq,sp−2
Qsq Psp−2
Zsq,sp−2
msq+sp−2

Ssq,sp−3
Qsq Psp−3
Zsq,sp−3
msq+sp−3

Ssq,sp−4
Qsq Psp−4
Zsq,sp−4
msq+sp−4

Si,j
Qi P j
Zi,j
Si+1,j−1
..

..
· · · mi+j · · ·

9

Psq
Qsq
Psq
Q1Psq
Z1,3

···
···

···

S2,sq−1
Q2Psq−1
Z2,sq−1

···

..

..

Ssq−1,1
Qsp−1P1
Zsq−1,1

Ssq−1,0
Qsq−1
0
msq−1

Ssq0
Qsq
0
msq

P3
Q3
P3
Q1P2
Z1,2

P2
Q2
P2
Q1 P 1
Z1,1

S2,1
Q2P1
Z2,1

S20
Q2
0
m2

...

P1
Q1
P1
Q11
0
m1

P0 = 1
Q0 = 1

m0 = 1

Each cell contains one of the total (sp − 1)(sq − 1) quadratic terms in the form of Qi Pj .
To chain a cell to its upper cell, one extra sum variable Si,j is added. Also, each carry-on
variable Zi,j in a cell is the carry-on of the cell directly to its right, so each cell contains
four variables. The sum of three terms Qi Pj , Si,j , and Zi,j is at most 3; thus, it generates an additional sum variable Si+1,j−1 and one carry-on variable Zi,j+1 . Therefore, the
equation for an arbitrary cell indexed (i, j), shown in the centre of the above table, is
Si,j + Qi Pj + Zi,j = Si+1,j−1 + 2Zi,j+1 .
As we can see, only six binary variables are involved in each cell equation and the equation
contains one quadratic term, so it can be transformed into a positive Hamiltonian without
adding slack variables. The Hamiltonian generation and reduction procedure is discussed
in detail in Results 2.1.

4.3

Gröbner bases

Good references for the following definitions are [Stu96, CLO98].
Normal forms. A normal form is the remainder of Euclidean divisions in the ring of
polynomials R[x1 , . . . , xn ]. Precisely, the normal form of a polynomial f ∈ R[x1 , . . . , xn ],
with respect to the set of polynomials B ⊂ R[x1 , . . . , xn ] (usually a Gröbner basis), is the
polynomial NF(f ) ∈ R[x1 , . . . , xn ], which is the image of f modulo B. It is the remainder
of the Euclidean of f by all g ∈ B.
Term orders. A term order on R[x1 , . . . , xn ] is a total order ≺ on the set of all monomials
xa = xa11 . . . xann , which has the following properties: (1) if xa ≺ xb , then xa+c ≺ xb+c
for all positive integers a, b, and c; (2) 1 ≺ xa for all strictly positive integers a. An
example of this is the pure lexicographic order plex(x1 , . . . , xn ). Monomials are compared
first by their degree in x1 , with ties broken by degree in x2 , etc. This order is usually
used in eliminating variables. Another example, is the graded reverse lexicographic order
tdeg(x1 , . . . , xn ). Monomials are compared first by their total degree, with ties broken by
reverse lexicographic order, that is, by the smallest degree in xn , xn−1 , etc.
Gröbner bases. Given a term order ≺ on R[x1 , . . . , xn ], then by the leading term (initial
term) LT of f we mean the largest monomial in f with respect to ≺. A (reduced) Gröbner
basis to the ideal I with respect to the ordering ≺ is a subset B of I such that: (1) the
initial terms of elements of B generate the ideal LT(I) of all initial terms of I; (2) for each
g ∈ B, the coefficient of the initial term of g is 1; (3) the set LT(g) minimally generates
LT(I); and (4) no trailing term of any g ∈ B lies in LT(I). Currently, Gröbner bases are
computed using sophisticated versions of the original Buchberger algorithm, for example,
the F4 algorithm by J. C. Faugère.

4.4

Factorization as an eigenvalue problem

In this section, for completeness, we describe how the factorization problem can be solved
using eigenvalues and eigenvectors. This is an adaptation of the method presented in

10

[PS01] to factorization, which is itself an adaption to real polynomial optimization of the
method of solving polynomial equations using eigenvalues in [CLO98].
Let H be in R[x1 , . . . , xn ] as in (2.5), where we have used the notation xi instead of the
P s, Qs, Zs, and W s. Define
X
Hα := H +
αi xi (xi − 1),
i

which is in the larger ring R[x1 , . . . , xn , α1 , . . . , αn ]. We also define the set of polynomials
C = {∂Hα /∂x1 , . . . , ∂Hα /∂xn , ∂Hα /∂α1 , . . . , ∂Hα /∂αn } .
The variety V(C) is the set of all binary critical points of H. Its coordinates ring is
the residue algebra A := R[x1 , . . . , xn , α1 , . . . , αn ]/C. We need to compute a basis for A.
This is done by first computing a Gröbner basis for C and then extracting the standard
monomials (i.e., the monomials in R[x1 , . . . , xn , α1 , . . . , αn ] that are not divisible by the
leading term of any element in the Gröbner basis). In the simple example below, we
do not need to compute any Gröbner basis since C is a Gröbner basis with respect to
plex(α, x). We define the linear map
mHα : A → A
g 7→ Hα g
Since the number of critical points is finite, the algebra A is always finite-dimensional by
the Finiteness Theorem ([CLO98]). Now, the key points are:
• The value of Hα (i.e., values of H), on the set of critical points V(C), are given by
the eigenvalues of the matrix mHα .
• Eigenvalues of mxi and mαi give the coordinates of the points of V(C).
• If v is an eigenvector for mHα , then it is also an eigenvector for mxi and mαi for
1 ≤ i ≤ n.
We illustrate this in an example. Consider M = pq = 5 × 3 and let
H = 2 + 7 x4 + 2 x3 + 2 x4 x3 − 2 x3 x 2 − x1 − 4 x4 x1 − 2 x3 x1 + x2 x1
be the corresponding Hamiltonian as in (2.5), where x1 = p2 , x2 = q1 , x3 = w2,1 , and
x4 = z2,3 . A basis for the residue algebra A is given by the set of the 16 monomials
{1, x4 , x3 , x4 x3 , x2 , x4 x2 , x3 x2 , x3 x2 x4 , x1 , x4 x1 , x3 x1 , x1 x3 x4 , x2 x1 , x4 x1 x2 , x1 x3 x2 , x1 x3 x2 x4 }.

11

The matrix mHα is

2 7

 0 9


 0 0


 0 0


 0 0

 0 0


 0 0


 0 0

mHα := 
 0 0


 0 0


 0 0


 0 0

 0 0


 0 0


 0 0

0 0

2

2

0 0 −2

0

4

0 0

4

9

0 0 −2

0

−1 −4 −2

0
−2

0

−5

0

0

0

1

0

0

−2 0

1

0

−3 −4 0

0

1

−7 0

0

0

0

0

0 13 0 0

0

−2

0

0

0

0

0

2 7

0

2

0

0

0

0

0 −4 −2

0

0

0 9

0

2

0

0

0

0

0 −4

0

0

0 0

2

9

0

0

0

0

0

0

−2

0

0

0 0

0

11

0

0

0

0

0

0

0

0

0

0 0

0

0

1

3

0

2

1

0

−2

0

0

0 0

0

0

0

4

0

2

0

1

0

0

0

0 0

0

0

0

0

1

5

0

0

−1

0

0

0 0

0

0

0

0

0

6

0

0

0

0

0

0 0

0

0

0

0

0

0

2

3

−2

0

0

0 0

0

0

0

0

0

0

0

5

0

0

0

0 0

0

0

0

0

0

0

0

0

0

0

0

0 0

0

0

0

0

0

0

0

0

0

0

0




0 


0 


1 


0 

−2 


−4 


−6 


0 


−2 


0 


−1 

2 


0 


5 

5

We expect the matrix’s smallest eigenvalue to be zero and, indeed, we get the following
eigenvalues for mHα :
{0, 1, 2, 4, 5, 6, 9, 11, 13}.
This is also the set of values which Hα takes on V(C). The eigenvector v which corresponds
to the eigenvalue 0 is the column vector
v := (1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0)T .
This eigenvector is used to find the coordinates of x̂ ∈ V(C) that cancel (minimize) Hα .
The coordinates of the global minimum x̂ = (x̂1 , . . . , x̂n ) are defined by mxi v = x̂i v, and
this gives x1 = x2 = x3 = 1, x4 = 0, and α1 = 2α2 = α3 = 2, α4 = 5.

12

5

Supplementary materials

5.1

Continuous optimization problems for the requirements (ii–
iii)

In Results 2.1, we describe how a positive quadratic polynomial Hij+ can be extracted
using Gröbner bases. Here we provide the details of the calculation.
The second requirement (ii) is equivalent to each of the following linear polynomials being
greater than zero:
a1 , −a1 + a3 , −a1 − a4 , −a1 + a5 , −a1 + a6 , 2 a1 + a3 , 2 a1 − a4 , −a2 − a1 , −a2 + 2 a1 ,
−2 a1 + a3 + 2 a5 , −2 a1 + a3 + 2 a6 , −2 a1 − a4 + 2 a5 , −2 a1 − a4 + 2 a6 , −a1 + a5 + a6 ,
a1 + a3 − a5 , a1 + a3 − a6 , a1 − a4 − a5 , a1 − a4 − a6 , −a2 − 2 a1 + 2 a5 , −a2 − 2 a1 + 2
a6 , −a2 + a1 − a5 , −a2 + a1 − a6 , −2 a1 + a3 + 2 a5 + 2 a6 , −2 a1 − a4 + 2 a5 + 2 a6 ,
a1 + a3 − a5 − a6 , a1 − a4 − a5 − a6 , −a2 − 2 a1 + a3 − a4 , −a2 − 2 a1 + 2 a5 + 2 a6 ,
−a2 + a1 − a5 − a6 , −a2 + 3 a1 + a3 − a4 , −a2 − 3 a1 + a3 − a4 + 3 a5 ,
−a2 − 3 a1 + a3 − a4 + 3 a6 , −a2 + 2 a1 + a3 − a4 − 2 a5 , −a2 + 2 a1 + a3 − a4 − 2 a6 ,
−a2 − 3 a1 + a3 − a4 + 3 a5 + 3 a6 , −a2 + 2 a1 + a3 − a4 − 2 a5 − 2 a6
For the third requirement (iii), the first choice for the objective function f : R5 → R is
2
2
f (a1 , . . . , a5 ) = (−a1 + a5 + a6 )2 − 1 + (−2 a1 + a3 + 2 a5 + 2 a6 )2 − 1
2
2
+ (a1 − a2 − a5 − a6 )2 − 1 + (a1 − a4 − a5 − a6 )2 − 1
2
2
2
2
+ 2 a2 2 − 1 + a1 2 − 1 + 2 a3 2 − 1 + 2 a4 2 − 1
2
2
2
2
+ 2 a5 2 − 1 + 2 a6 2 − 1 + 4 a5 2 − 1 + 4 a6 2 − 1
The solution is a1 = 0.214, a2 = −1.082, a3 = 0.514, a4 = −1.082, a5 = 0.314, and
a6 = 0.314.
The second choice for f is
2
+ (−2 a1 + a3 + 2 a5 + 2 a6 )2 − a2
2
2
+ (a1 − a2 − a5 − a6 )2 − a2 + (a1 − a4 − a5 − a6 )2 − a2
2
2
2
2
+ 2 a2 2 − a2 + a1 2 − a2 + 2 a3 2 − a2 + 2 a4 2 − a2
2
2
2
2
+ 2 a5 2 − a2 + 2 a6 2 − a2 + 4 a5 2 − a2 + 4 a6 2 − a2
(5.1)

f (a1 , . . . , a5 ) =

(−a1 + a5 + a6 )2 − a2

2

The solution is a1 = 1.0, a2 = −4.0, a3 = 4.0, a4 = −4.0, a5 = 2.0, and a6 = 2.0 (identical
to the solution given in [SS10]).

5.2

Basic description of the quantum annealing processor

Here we introduce the quantum annealing concept that ultimately solves a general Ising
(quadratic unconstrained binary optimization, or ”QUBO”) problem, then talk about the
13

important topic of embedding a QUBO problem into the specific quantum annealer (the
D-Wave 2X processor).
Quantum annealing (QA), along with the D-Wave processors, have been the focus of
much research. We refer the interested reader to [JAG+ 11, CCD15, BAS+ 13, BRI+ 14,
LPS+ 14]. QA is a paradigm designed to find the ground state of systems of interacting
spins represented by a time-evolving Hamiltonian:
S(s) = E(s)HP −
HP = −

X

1X
∆(s)σix ,
2 i

hi σix +

i

X

Jij σiz σjz .

i<j

The parameters hi and Jij encode the particular QUBO problem P into its Ising formulation. QA is performed by first setting ∆  E, which results in a ground state into which
the spins can be easily initialized. Then ∆ is slowly reduced and E is increased until
E  ∆. At this point the system is dominated by HP , which encodes the optimization
problem. Thus, the ground state represents the solution to the optimization problem.
An embedding is the mapping of the nodes of an input graph to the nodes of the destination graph. The graph representing the problem’s QUBO matrix needs to be embedded
into the actual physical qubits on the processor in order for it to solve the QUBO problem. The specific existing connectivity pattern of qubits in the D-Wave chip is called the
Chimera graph. Embedding an input graph (a QUBO problem graph) into the hardware
graph (the Chimera graph) is in general NP-hard ([Cho08]).
Figure 1–right shows an embedding of the (column algorithm) QUBO corresponding to
the bi-prime M = 200 099 into the Chimera graph of the D-Wave 2X chip consisting of
a 12 by 12 lattice of 4 by 4 bipartite blocks. The Chimera graph is structured so that
the vertical and horizontal couplers in its lattice are connected only to either side of each
bipartite block. Each node in this graph represents one qubit and each edge represents
a coupling between two qubits. Adjacent nodes in the Chimera graph can be grouped
together to form new effective (i.e., logical) nodes, creating nodes of a higher degree. Such
a grouping is performed on the processor by setting the coupler between two qubits to a
large negative value, forcing two Ising spins to align such that the two qubits end up with
the same values. These effective qubits are expected to behave identically and remain in
the same binary state at the time of measurement. The act of grouping adjacent qubits
(hence forming new effective qubits) is called chain creation or chain identification.
An embedding strategy consists of two tasks: mapping and identification. Mapping is the
assignment of the nodes of the input graph to the single or effective nodes of the destination
graph. Solving such problems optimally is in general NP-hard, but one can devise various
approximations and enhancement strategies to overcome these difficulties, for example,
using statistical search methods like simulated annealing, structure-based methods, or a
combination of both. For a better understanding of current embedding approaches, we
refer the reader to [Cho08, BCI+ 14, JWA14, TAA15]. In Figure 1–right, the blue lines
14

indicate the identified couplers, the yellow lines indicate the problem couplers (i.e., the
edges of the problem graph), and the grey lines indicate empty couplers.

5.3

Embedding and solving details

We have used one of the D-Wave 2X processors, DW2X SYS4, as our quantum annealing
solver. This processor operates at a temperature range of 26(±5) millikelvin (mK) and
has 1100 qubits with a 95.5-qubit yield. To utilize the processor, we used D-Wave’s
SAPI software development kit (version 2.2.1). To embed the problem graph into the
hardware graph we used the sapiFindEmbedding and sapiEmbedProblem modules, and
to solve the problems we used the sapiSolveIsing and sapiUnembedAnswer modules. For
all problems we opted for the maximum number of reads available (10 000) in order to
increase the fraction of ground state samples. The following table shows some statistics
of the embedding and solving stages for several of the highest numbers that we were able
to successfully embed and solve.
Embedding & Solving Statistics

M
31861
34889
150419
151117
174541
200099

n emT ry idC prC #qubits jRatio rT ime
95
33
848 721
815
10
3.52
95
27
803 740
833
10
3.52
73
1
941 830
902
64
3.52
72
7
1001 846
918
64
3.52
72
3
1004 897
966
64
3.52
75
5
884 824
897
64
3.52

In the above table, M stands for the bi-prime, n is the number of variables in the QUBO
problem, emTry is the number of block trials of the sapiFindEmbedding routine, idC
is the total number of identified couplers, prC is the total number of problem couplers,
max({|J |})
#qubits is the total number of (physical) qubits, jRatio is the ratio min({|Jijij|}) , and rTime
is the chip run time in seconds.

6

Acknowledgements

We would like to thank Pooya Ronagh for constructive discussions and helpful comments
on the paper. We also thank Marko Bucyk for proofreading the manuscript.

References
[BAS+ 13]

Sergio Boixo, Tameem Albash, Federico M. Spedalieri, Nicholas Chancellor,
and Daniel A. Lidar, Experimental signature of programmable quantum annealing, Nat Commun 4 (2013).

[BCI+ 14]

Zhengbing Bian, Fabian Chudak, Robert Israel, Brad Lackey, William G
Macready, and Aidan Roy, Discrete optimization using quantum annealing
on sparse ising models, Frontiers in Physics 2 (2014), no. 56.

15

[BRI+ 14]

Sergio Boixo, Troels F. Ronnow, Sergei V. Isakov, Zhihui Wang, David
Wecker, Daniel A. Lidar, John M. Martinis, and Matthias Troyer, Evidence for
quantum annealing with more than one hundred qubits, Nat Phys 10 (2014),
no. 3, 218–224.

[Bur02]

C.J.C. Burges, Factoring as optimization, Tech. Report MSR-TR-2002-83,
Microsoft Research, January 2002.

[CCD15]

Cristian S. Calude, Elena Calude, and Michael J. Dinneen, Guest column:
Adiabatic quantum computing challenges, SIGACT News 46 (2015), no. 1,
40–61.

[Cho08]

Vicky Choi, Minor-embedding in adiabatic quantum computation: I. the parameter setting problem, Quantum Information Processing 7 (2008), no. 5,
193–209.

[CLO98]

David A. Cox, John B. Little, and Donal O’Shea, Using algebraic geometry,
Graduate texts in mathematics, Springer, New York, 1998.

[JAG+ 11]

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 (2011), no. 7346, 194–198.

[JWA14]

Cai Jun, G. Macready William, and Roy Aidan, A practical heuristic for
finding graph minors, Preprint arXiv:1406.2741 (2014).

[Kit95]

A.Yu. Kitaev, Quantum measurements and the abelian stabilizer problem,
quant-ph/9511026 (1995).

[LPS+ 14]

T. Lanting, A. J. Przybysz, A. Yu. Smirnov, F. M. Spedalieri, M. H. Amin,
A. J. Berkley, R. Harris, F. Altomare, S. Boixo, P. Bunyk, N. Dickson, C. Enderud, J. P. Hilton, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, R. Neufeld, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva,
S. Uchaikin, A. B. Wilson, and G. Rose, Entanglement in a quantum annealing processor, Phys. Rev. X 4 (2014), 021041.

[MNM+ 16] Thomas Monz, Daniel Nigg, Esteban A. Martinez, Matthias F. Brandl,
Philipp Schindler, Richard Rines, Shannon X. Wang, Isaac L. Chuang, and
Rainer Blatt, Realization of a scalable shor algorithm, Science 351 (2016),
no. 6277, 1068–1070.
[PS01]

Pablo A. Parrilo and Bernd Sturmfels, Minimizing polynomial functions, DIMACS Series in Discrete Mathematics and Theoretical Computer Science
(2001).

[Rau13]

Robert Raussendorf, Contextuality in measurement-based quantum computation, Phys. Rev. A 88 (2013), 022322.
16

[Sho97]

Peter W. Shor, Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer, SIAM J. Comput. 26 (1997), no. 5,
1484–1509. MR 1471990

[SS10]

Gernot Schaller and Ralf Schutzhold, The role of symmetries in adiabatic quantum algorithms, Quantum Information & Computation 10 (2010),
no. 1&2, 109–140.

[SSV13]

John A. Smolin, Graeme Smith, and Alexander Vargo, Oversimplifying quantum factoring, Nature 499 (2013), no. 7457, 163–165.

[Stu96]

Bernd Sturmfels, Gröbner bases and convex polytopes, University Lecture
Series, vol. 8, American Mathematical Society, Providence, RI, 1996. MR
1363949

[TAA15]

Boothby Tomas, D. King Andrew, and Roy Aidan, Fast clique minor generation in chimera qubit connectivity graphs, Preprint arXiv:1507.04774 (2015).

[XZL+ 12]

Nanyang Xu, Jing Zhu, Dawei Lu, Xianyi Zhou, Xinhua Peng, and Jiangfeng
Du, Quantum factorization of 143 on a dipolar-coupling nuclear magnetic
resonance system, Phys. Rev. Lett. 108 (2012), 130501.

17

