An adiabatic quantum algorithm for the Frobenius problem

TL;DR

This paper introduces an adiabatic quantum algorithm tailored for the NP-hard Frobenius problem, translating it into a QUBO formulation to leverage D-Wave quantum hardware.

Contribution

It presents a novel quantum algorithm for the Frobenius problem that utilizes the Apéry set and is optimized for D-Wave quantum annealers.

Findings

Algorithm successfully encodes the Frobenius problem as a QUBO.

Designed specifically for the D-Wave 2X quantum annealer.

Provides a new approach to solving a classical NP-hard problem using quantum computing.

Abstract

The (Diophantine) Frobenius problem is a well-known NP-hard problem (also called the stamp problem or the chicken nugget problem) whose origins lie in the realm of combinatorial number theory. In this paper we present an adiabatic quantum algorithm which solves it, using the so-called Ap\'ery set of a numerical semigroup, via a translation into a QUBO problem. The algorithm has been specifically designed to run in a D-Wave 2X machine.

Tables

Table 1. Table 1: Combinatorial invariants for some examples of semigroups

$S=⟨A⟩$	$m\left(S\right)$	$e\left(S\right)$	$G\left(S\right)$	$g\left(S\right)$	$f\left(S\right)$
$⟨3,7⟩$	$3$	$2$	$\{1,2,4,5,8,11\}$	$6$	$11$
$⟨4,9⟩$	$4$	$2$	$\{1,2,3,5,6,7,10,11,14,15,19,23\}$	$12$	$23$
$⟨5,8,11⟩$	$5$	$3$	$\{1,2,3,4,6,7,9,12,14,17\}$	$10$	$17$
$⟨5,7,9⟩$	$5$	$3$	$\{1,2,3,4,6,8,11,13\}$	$8$	$13$

Equations

$$\sum _{j=1}^n{\lambda }_i{a}_i\text{ with }{\lambda }_i\in {\mathbb{Z}}_{\ge 0}.$$

$$S=⟨A⟩=⟨a_1,...,a_n⟩=\{{\lambda }_1{a}_1+...+{\lambda }_n{a}_n|{\lambda }_i\in {\mathbb{Z}}_{\ge 0}\}$$

$$\begin{array}{r}⟨3,7⟩=\{0,3,6,7,9,10,12,\to \}\\ ⟨4,9⟩=\{0,4,8,9,12,13,16,17,18,20,21,22,24,\to \}\\ ⟨5,8,11⟩=\{0,5,8,10,11,13,15,16,18,\to \}\\ ⟨5,7,9⟩=\{0,5,7,9,10,12,14,\to \}\end{array}$$

$$m(S)=a_1\text{ and }e(S)=n$$

$$m(S)=\min S^{\ast }$$

$$G(S)={\mathbb{Z}}_{\ge 0}\backslash S.$$

$$g(S)=|G(S)|,$$

$$f(S)=\max G(S),$$

$$Ap(S,s)=\{x\in S|x-s\notin S\}.$$

$$S=⟨5,8,11⟩=\{0,5,8,10,11,13,15,16,18,\to \}.$$

$$\begin{array}{rl}{\omega }_0& =\min \{x\in S|x\equiv 0\mathrm{mod}5\}=0\\ {\omega }_1& =\min \{x\in S|x\equiv 1\mathrm{mod}5\}=11\\ {\omega }_2& =\min \{x\in S|x\equiv 2\mathrm{mod}5\}=22\\ {\omega }_3& =\min \{x\in S|x\equiv 3\mathrm{mod}5\}=8\\ {\omega }_4& =\min \{x\in S|x\equiv 4\mathrm{mod}5\}=19\end{array}$$

$$e(S)\le m(S).$$

$$f(S)=\max Ap(S,s)-s$$

$$g(S)=\frac{1}{s}\left(\sum _{\omega \in Ap(S,s)}\omega \right)-\frac{s-1}{2}$$

$$\sum _{i=1}^n{\lambda }_i{a}_i=t.$$

$$H(s_1,\dots ,s_n)=\sum _{i=1}^n{h}_i{s}_i+\sum _{i=1}^{n-1}\sum _{j=i+1}^n{J}_{ij}{s}_i{s}_j,$$

$$H(s_1,\dots ,s_n)=\sum _{i=1}^n{h}_i{s}_i+\sum _{i=1}^{n-1}\sum _{j=i+1}^n{J}_{ij}{s}_i{s}_j$$

$$Ap(S,s)=\{x\in S|x-s\notin S\}.$$

$${\omega }_i=\min \{x\in S:x\equiv i\mathrm{mod}s\}.$$

$$\begin{array}[]{rl}\text{min}&\displaystyle\sum\limits_{j=1}^{n}a_{j}x_{j}\\ \textit{subject to:}&\displaystyle\sum\limits_{j=1}^{n}a_{j}x_{j}=i+sk,\\ \end{array}$$

$$65=0\times (11)+1\times (19)+2\times (23).$$

$$f(S)=\max \{Ap(S,s)\}-s.$$

$$f(S)\le (a_1-1)(a_n-1)-1.$$

$$t_j=1+\lfloor {\log }_2\left(\frac{(a_1-1)(a_n-1)+s-1}{a_j}\right)\rfloor$$

$$u=1+\lfloor {\log }_2\left(\frac{(a_1-1)(a_n-1)+s-i-1}{s}\right)\rfloor$$

$$\begin{array}[]{rll}\text{min}&\displaystyle\sum\limits_{j=1}^{n}a_{j}\displaystyle\sum\limits_{l=0}^{t_{j}}2^{l}x_{jl}\\ \textit{subject to:}&\displaystyle\sum\limits_{j=1}^{n}a_{j}\displaystyle\sum\limits_{l=0}^{t_{j}}2^{l}x_{jl}=i+s\sum\limits_{m=0}^{u}2^{m}k_{m}\end{array}$$

$$Q(x_1,\dots ,x_n)=c_0+\sum _{i=1}^n{c}_i{x}_i+\sum _{1\le i<j\le n}{q}_{ij}{x}_i{x}_j,$$

$$\sum _{j=1}^n{a}_j{x}_j=i+sk,$$

$$\begin{array}[]{rl}\text{min}&\displaystyle\sum\limits_{j=1}^{n}a_{j}x_{j}+\lambda^{(\nu)}\left(\displaystyle\sum\limits_{j=1}^{n}a_{j}x_{j}-i-sk\right)^{2}\\ \end{array}$$

$${\lambda }^{(\nu +1)}={\lambda }^{(\nu )}+\alpha \cdot \left|\sum _{j=1}^n{a}_j{x}_j-i-sk\right|$$

Full text

1 Numerical semigroups and the Frobenius problem

The study of numerical semigroups has its origins at the end of the 19th Century, when English mathematician James Joseph Sylvester (1814 – 1897) and German mathematician Ferdinand Georg Frobenius (1849 – 1917) were both interested in what is now known as the Frobenius problem, which we proceed to enunciate.

Definition 1.1**.**

Let $a_{1},a_{2},\ldots,a_{n}\in\mathbb{Z}_{\geq 0}$ with $\gcd(a_{1},a_{2},\ldots,a_{n})=1$, the Frobenius problem, or FP, is the problem of finding the largest positive integer that cannot be expressed as an integer conical combination of these numbers, i.e., as a sum

[TABLE]

This problem, so easy to state, can be extremely complicated to solve in the majority of cases, as will be seen later. It can be found in a wide variety of contexts, being the most famous the problem of finding the largest amount of money which cannot be obtained with a certain set of coins: if, for example, we have an unlimited amount of coins of 2 and 5 units, we can represent any quantity except 1 and 3.

In order to understand the relationship between this problem and numerical semigroups, we shall first define the latter.

Definition 1.2**.**

A semigroup is a pair $(S,+)$, where $S$ is a set and $+$ is a binary operation $+:S\times S\rightarrow S$ that is associative.

Definition 1.3**.**

A numerical semigroup $S$ is a subset of the non-negative integers $\mathbb{Z}_{\geq 0}$ which is closed under addition, contains the identity element [math], and has a finite complement in $\mathbb{Z}_{\geq 0}$.

From now on, we shall denote numerical semigroups as $S$, taking for granted that they are commutative and that their associated operation is the addition. As it can be easily noted, a numerical semigroup is trivially a semigroup. In other words, a numerical semigroup is a semigroup that, additionally, is a monoid (i.e., it also has an identity element) and has finite complement in $\mathbb{Z}_{\geq 0}$. In order to work with numerical semigroups, it will be necessary to characterize them somehow. For that, let us set forth the following lemma. The proof of this result (and of any other result in this section, unless otherwise stated) can be found in [24].

Lemma 1.4**.**

Let $A=\{a_{1},...,a_{n}\}$ be a nonempty subset of $\mathbb{Z}_{\geq 0}$. Then,

[TABLE]

is a numerical semigroup if and only if $\gcd(a_{1},...,a_{n})=1$.

The previous lemma tells us that, drawing from a set $A\subseteq\mathbb{Z}_{\geq 0}$, it is possible to generate a semigroup $S=\langle A\rangle$ as long as the elements of $A$ satisfy a certain condition. In this context, any set $A$ such that $S=\langle A\rangle$ for a certain numerical semigroup $S$ is called a system of generators of $S$. Even more, it can be proved that any numerical semigroup can be expressed that way, as will be shown next.

Theorem 1.5**.**

Every numerical semigroup $S$ admits a unique minimal system of generators, which can be calculated as $S^{*}\setminus(S^{*}+S^{*})$ with $S^{*}=S\setminus\{0\}$.

Corollary 1.6**.**

Let $S$ be a numerical semigroup generated by $A=\{a_{1},\ldots,a_{n}\}$ with $0\neq a_{1}<a_{2}<...<a_{n}$. Then, $A$ is a minimal system of generators of $S$ if and only if $a_{i+1}\notin\langle a_{1},a_{2},\ldots,a_{i}\rangle$, for all $i\in\{1,\ldots,n-1\}$.

It will be stated later in this section that the minimal system of generators of a numerical semigroup is in fact finite. Hereinafter, if we say that $S=\langle A\rangle$ with $A=\{a_{1},a_{2},\ldots,a_{n}\}$ is a numerical semigroup, then we shall assume without loss of generality that $a_{1}<a_{2}<\cdots<a_{n}$, $\gcd(a_{1},a_{2},\ldots,a_{n})=1$, and that $A$ is the minimal system of generators of $S$. Some examples of numerical semigroups, which will be used for the rest of the section, are:

[TABLE]

Where $\rightarrow$ means that all integers thenceforth are also included in the numerical semigroup. Having thus defined and characterized what numerical semigroups are, we proceed to describe some of their combinatorial invariants.

Definition 1.7**.**

Let $S=\langle a_{1},a_{2},\ldots,a_{n}\rangle$ be a numerical semigroup, then

[TABLE]

are called respectively the multiplicity of $S$ and the embedding dimension of $S$.

Lemma 1.8**.**

Let $S$ be a numerical semigroup, then

[TABLE]

Definition 1.9**.**

The set of gaps of a numerical semigroup $S$ is defined as

[TABLE]

Its cardinal,

[TABLE]

is called the genus of $S$; and its maximum,

[TABLE]

is called the Frobenius number of $S$.

In other words, as by definition any numerical semigroup $S$ has a finite complement in $\mathbb{Z}_{\geq 0}$, we can define the maximum of such complement as $f(S)$, known as the Frobenius number. In fact, the Frobenius problem described at the beginning of this chapter in Definition 1.1 can be enunciated as the problem of finding $f(S)$ for a certain numerical semigroup $S$. We shall expand the concepts related to the difficulties that surround the calculation of the Frobenius number later on. Table 1 shows the combinatorial invariants associated for the previously given examples of semigroups.

Greek-French mathematician Roger Apéry (1916 – 1994), better known for proving in 1979 the irrationality of $\zeta(3)$ [5], also laid the background in the context of the resolution of singularities of curves [4] for an important set associated to a numerical semigroup $S$ and one of its elements.

Definition 1.10**.**

The Apéry set of a numerical semigroup $S$ with respect to a certain $s\in S^{*}$ is defined as

[TABLE]

A possible characterization of the elements of the Apéry set one by one, which will come to special use in the last section, is given by the following lemma.

Lemma 1.11**.**

Let $S$ be a numerical semigroup and let $s\in S^{*}$. Then, $Ap(S,s)=\{\omega_{0},\omega_{1},...,\omega_{s-1}\}$ where $\omega_{i}$ is the least element of $S$ congruent with $i$ modulo $s$, for all $i\in\{0,...,s-1\}$. Consequently, $|Ap(S,s)|=s$.

By means of an example on how to calculate the Apéry set of a semigroup with respect to a certain element, let

[TABLE]

Then $Ap(S,5)=\{\omega_{0},\ldots,\omega_{4}\}$, where

[TABLE]

We proceed to hint how the Apéry set proves some elementary (but nevertheless important) results of numerical semigroups.

Proposition 1.12**.**

The minimal system of generators of a numerical semigroup $S$ is finite.

Proof.

Let $s\in S^{*}$. Then, it is easy to see that for every $t\in S$ there exists a unique pair $(u,v)\in\mathbb{Z}_{\geq 0}\times Ap(S,s)$ such that $t=us+v$. Thus, $S$ is generated by $A=Ap(S,s)\ \cup\ \{s\}$ and, as $A$ is finite, the unique minimal system of generators must be finite. ∎

Lemma 1.13**.**

Let $S$ be a numerical semigroup, then

[TABLE]

Proof.

Let $a=m(S)$ and let $A=Ap(S,a)\ \cup\ \{a\}$. Thus, as $S$ can be generated by $A\setminus\{0\}$ and $|A\setminus\{0\}|=a$, we can conclude that $e(S)\leq m(S)$. ∎

The Apéry set is noteworthy in the context of the Frobenius problem as there exists a relationship between its members (regardless of the element $s\in S^{*}$ we choose) and the Frobenius number, which we proceed to enunciate and prove.

Theorem 1.14**.**

(A. Brauer – J. E. Shockley, 1962) [10]* Let $S$ be a numerical semigroup and let $s\in S^{*}$. Then*

[TABLE]

[TABLE]

Now we exemplify the relationship between numerical semigroups and combinatorial optimization, as one of the most important problems in the latter branch of mathematics, known as the knapsack problem or rucksack problem, and more concretely one of its variants [21] (p. 374), can be seen as the problem of deciding if a given integer $t$ belongs to a certain numerical semigroup $S$.

Definition 1.15**.**

The numerical semigroup membership problem, or NSMP, is the problem of determining if, given a certain integer $t\in\mathbb{Z}_{\geq 0}$ and a numerical semigroup $S=\langle a_{1},...,a_{n}\rangle$, the integer $t$ is contained in $S$. That is to say, if there exist non-negative integers $\lambda_{1},\ldots,\lambda_{n}\in\mathbb{Z}_{\geq 0}$ such that

[TABLE]

The numerical semigroup membership problem is in NP-complete, as shown in [21]. This fact was used by Jorge Ramírez-Alfonsín in 1996 [22] to finally prove that the Frobenius problem is in NP-hard (under Turing reductions). For that, he define a polynomial algorithm $\Lambda_{\text{NSMP}}$ for solving the NSMP which uses as a subroutine an unknown algorithm $\Lambda_{\text{FP}}$ that solves the Frobenius problem. Thus, he proved that the NSMP can be Turing reduced to the FP. As the NSMP is in NP-complete, he concluded that the FP is in NP-hard.

2 The Ising spin problem in adiabatic quantum computing

Our tackle of the Frobenius problem uses adiabatic quantum computing, more specifically that related to the Ising spin problem. This problem has a story where NP-hard problems are no strangers.

A form of adiabatic quantum computation [19] is quantum annealing, where a known initial configuration of a quantum system evolves towards the ground state of a Hamiltonian that encodes the solution of an NP-hard optimization problem. The Canadian company D-Wave Systems announced in 2011 the first commercially available quantum annealer, composed of arrays of eight superconducting flux quantum bits with programmable spin–spin couplings, and published their results [18]. Subsequently in 2013, S. Boixo et al. [7] published their experimental results on the 108-qubit D-Wave One device. Their last chip, released in 2017 and called D-Wave 2000Q, has 2,048 qubits in a Chimera graph architecture [2], and can be seen as a computer that solves the Ising spin problem, a particular type of quantum annealing which we proceed to describe.

The Ising spin model, originally formulated by physicist Wilhelm Lenz and first solved by his student, Ernst Ising [17], consists of a model of ferromagnetism in statistical mechanics in which we have to find the ground state of a system of $n$ interacting spins. If we represent the spins of these particles as binary variables $s_{i}\in\{-1,1\}$ with $i\in\{1,\ldots,n\}$, then the Ising Spin problem can be expressed as an integer optimization problem whose objective is to find the spin configuration of the system that minimizes the function

[TABLE]

where $h_{i}\in\mathbb{R}$ is the energy bias acting on particle $i$ (i.e. the external forces applied to each of the individual particles) and $J_{ij}\in\mathbb{R}$ is the coupling energy between the spins $i$ and $j$ (i.e. the interaction forces between adjacent spins). This problem was proved to be in NP-hard by Francisco Barahona [6], and can be effectively solved with the hardware implemented by D-Wave, whose chip permits to program independently the values of $h_{i}$ and $J_{ij}$ [18, 7]. It will be shown in the next section how to embed a certain optimization problem inside the D-Wave machine.

3 An adiabatic quantum algorithm for the Frobenius problem

The next algorithm shows the possibilities of calculating the Apéry set and the Frobenius number of a numerical semigroup with an actual adiabatic quantum computer. As described in the previous section, current D-Wave quantum annealers solve a certain kind of mathematical optimization problems known as Ising spin problems, namely

[TABLE]

with $s_{i}\in\{-1,1\}$, where the objective is to find the minimum of $H(s_{1},\ldots,s_{n})$.

On the other hand, we recall that the Apéry set with respect to $s\in S\setminus\{0\}$, where $S=\langle a_{1},\ldots,a_{n}\rangle$ is a numerical semigroup, is defined as

[TABLE]

The question is, how could we be able to translate this problem to the Ising spin model solved by the D-Wave machine? We have already shown in Lemma 1.11 that $Ap(S,s)=\{\omega_{0},\ldots,\omega_{n-1}\}$, where

[TABLE]

How about transforming the definition of $\omega_{i}$ into the answer of a mathematical optimization problem?

Definition 3.1**.**

[21]* An integer linear program, or ILP, is defined in its canonical form as the optimization problem:*

[TABLE]

where $x\in\mathbb{Z}_{\geq 0}^{n}$, $A\in\mathbb{Z}^{n}\times\mathbb{Z}^{m}$, $b\in\mathbb{Z}^{m}$ and $c\in\mathbb{Z}^{n}$.

Thus, it is straightforward to see that the calculation of each of the $\omega_{i}$ can be redefined as the following ILP:

[TABLE]

with $x_{j}\in\mathbb{Z}_{\geq 0}$ for all $j=1,...,n$ and $k\in\mathbb{Z}_{\geq 0}$.

This mathematical optimization problem represents a way of calculating the Apéry set in its own right and, although integer linear programming is in NP-hard and its recognition version (deciding whether $Ax=b$ has a feasible solution or not regardless of its optimality) is in NP-complete [21], this computational complexity arises in the general case. In our context, it may prove to be an easier problem; however, details on this remain to be worked out. This approach for the calculation of the Apéry set first appeared in [16], although Greenberg’s work dealt with the direct calculation of the Frobenius number by means of $Ap(S,a_{1})$, as will be discussed later.

We have tested this algorithm using state-of-the-art optimization software for integer optimization. In our case, the optimization problem for $\omega_{i}$ was modeled using AMPL [13, 14], an algebraic modeling language for solving mathematical optimization problems. AMPL has a clear advantage, as our problem is entirely parameterized and AMPL allows us to describe the generic problem in a .mod file while defining the actual values for the parameters in a separate .dat file. This way, we have just to change the parameters inside the .dat file in order to solve a new instance of the Apéry set. The .mod file we propose is shown in Figure 1 (the names of the parameters and variables of the problem are maintained so that no further explanation is required).

Let us suppose that we are interested in the numerical semigroup $S=\langle 11,19,23\rangle$, and that we want to calculate the Apéry set of $s=30$ (i.e., $Ap(S,30)$). The calculation of, for example, $\omega_{5}\in Ap(S,30)$, will have the associated .dat file depicted in Figure 2. In order to compute $\omega_{5}$, we also need a .run file that will load the model and the data of the problem, and which will also call the solver for solving the problem. In our case, the solver we have chosen is Gurobi [3] which, among other things, solves integer linear problems. The .run file we have used is shown in Figure 3, and the corresponding output we obtain can be seen in Figure 4. This output tells us that $\omega_{5}=65$ and also that a representation of 65 with respect to the generators of the semigroup is

[TABLE]

We can also write an alternative .run file that will directly calculate and display the whole content of the Apéry set for a certain $s\in S$. First, we have to drop the param i := 5; line from the .run file, as it will be changed in each iteration of the main loop, and modify the .run so that it will look like the one in Figure 5. Thus, we obtain $Ap(S,30)=\{0,11,19,22,23,33,34,38,42,44,45,46,55,56,$ $57,61,65,66,67,69,77,78,80,84,88,89,92,100,103,111\}$ (see Figure 6 for the actual output).

This algorithm also provides a way for calculating the Frobenius number of a numerical semigroup. We recall that, for any numerical semigroup $S$ and any integer $s\in S\setminus\{0\}$, then

[TABLE]

Thus, as the difficulty of obtaining $f(S)$ this way increments with respect to the number $s$ we choose (we have to solve $s$ ILPs), the smartest way of proceeding is by solving it in the case where $s=\min(S\setminus\{0\})$ (i.e., $s=a_{1}$, as done by [16]). The AMPL file that represents this approach is shown in Figure 7. In our case, the output (Figure 8) tells us that $f(S)=81$. All these files can be found in the public GitHub repository [20]; however, a license for both AMPL and Gurobi is needed in order to run them and obtain the same results (or any result at all).

What we have shown is a classical algorithm for obtaining the Apéry set and the Frobenius number in a general manner that depends on a black box that solves an ILP with global optimality. From now on we will explain the steps we have followed in order to solve this ILP with an adiabatic quantum computer and, specifically, with a D-Wave 2X machine (and also the obstacles we have encountered in that path). In order to transform this ILP problem into the Ising model solved by the D-Wave hardware, one step further involves changing its integer variables into a new set of binary variables with at most a polynomial cost.

Definition 3.2**.**

[21]* A binary linear program, or 0-1LP, is defined in its canonical form as:*

[TABLE]

where $x\in\{0,1\}^{n}$, $A\in\mathbb{Z}^{n}\times\mathbb{Z}^{m}$, $b\in\mathbb{Z}^{m}$ and $c\in\mathbb{Z}^{n}$.

Both problems are polynomially equivalent, as shown in [21] (Theorem 13.6), where an upper bound for the number of binary variables representing each integer variable from the original ILP problem is given. However, we can tight this number of binary variables by using our knowledge of the problem and one of the lower bounds of the Frobenius number given in [23] (ideally, we could use $f(S)$). In 1935, Russian mathematician Issai Schur proved in a lecture in Berlin [9, 23] the following result:

Theorem 3.3**.**

(I. Schur, 1935) Let $S=\langle a_{1},...,a_{n}\rangle$ be a numerical semigroup. Then,

[TABLE]

Thus, if we define

[TABLE]

for every $j\in\{1,\ldots,n\}$, and

[TABLE]

we have given an upper bound for the number of bits needed to describe every variable $x_{j}$ in the worst case (we can also run every known bound for the Frobenius number, and stick to the minimum of them, but for now let us just use Schur’s result). It follows that the transformation of our problem from ILP form to 0-1LP delivers the minimization problem shown below:

[TABLE]

with $x_{jl}\in\{0,1\}$ and $k_{m}\in\{0,1\}$ for all $j$, $l$ and $m$ contained in the summation indices.

Let us come back to the definition of the Ising model. In it, we have binary variables which can have the values -1 or 1. On the other hand, in our 0-1LP for calculating the elements of the Apéry set the variables can be in 0 or 1. It is usually more advantageous to redefine the Ising spin problem as a quadratic unconstrained binary optimizacion problem, or QUBO, which consists of the minimization of the objective function

[TABLE]

with $x_{i}\in\{0,1\}$, thus changing the possible states of its variables from $1$ and $-1$ to $1$ and [math]. The transformation follows from $s_{i}=1-2x_{i}$, and can be easily checked. In fact, the documentation from D-Wave allows to program our problem directly into a QUBO formulation [1], besides the Ising model formulation. One detail remains to be figured out: QUBO problems have no constraints, but ours have one. How to proceed?

The most common way to transform a constrained optimization problem into an unconstrained one is by penalizing the constraints, putting them in the objective function. In this manner, we can force the fulfillment of the constraints of the problem by punishing the error committed in them. In our example, as we have a single equality constraint, namely

[TABLE]

we can force it to the objective function this way:

[TABLE]

where $\lambda^{(\nu)}$ is updated after every iteration in the following way, until no change in $\lambda$ is obtained:

[TABLE]

with $\alpha>0$.

In 0-1LP form, we obtain the following QUBO problem:

[TABLE]

with its corresponding penalty updating formula:

[TABLE]

with $\alpha>0$.

Now, we are going to test the performance of this new procedure for obtaining the Apéry set (and the Frobenius number) via a classical solver. For that, we have just to modify the original apery_set_member.mod file into a new one with an unconstrained ILP, with the constraint penalized in the objective function as previously exaplained. This new AMPL file can be found in Figure 9, along with its corresponding .run file in Figure 10. Additionaly, we are going to keep track of the number of iterations of the algorithm (i.e., $\nu$), and on the values of $\lambda^{(\nu)}$ needed for the algorithm to stop. This way, if we start with $\lambda^{(0)}=0$, the output for this algorithm is found in Figure 11. On the other hand, if we start with $\lambda^{(0)}=100$, the output is shown in Figure 12.

Foreseeably, the Apéry set we obtain for $S$ and $s=30$ is correct, but that should not be a surprise as both optimization problems are equivalent to the original one. What is of interest here is that a low starting point for $\lambda$ implies a relatively large number of iterations, which translates in more calls to our integer linear programming solver than in the original constrained version. However, if we start with a large enough value of $\lambda$, we can obtain the same results without additional calls to the solver and just one iteration for each of the elements of the Apéry set. On the other hand, if we exceed in the value for $\lambda$, we would also need just one iteration per $\omega_{i}$, but the solver would need more time in order to find the global optimum. The latter behavior can be observed when taking into account the number of simplex iterations and branch-and-cut nodes that Gurobi needs for each of the problems.

Nevertheless, this unconstrained reformulation of the problem is not meant to be solved by Gurobi or any other integer linear programming solver, but rather by an ideal adiabatic quantum computer, and more concretely by one of the few D-Wave machines that currently exist worldwide. Our unconstrained problem, in QUBO form, is shown below (please note that, if $x\in\{0,1\}$, then $x=x^{2}$):

[TABLE]

[TABLE]

where

[TABLE]

Right now, an ideal quantum annealer should be able to solve our QUBO problem for finding the members of the Apéry set, provided that such a computer has a sufficiently large enough amount of qubits, and also that the graph connecting the qubits is complete. However, to date there are no quantum annealers that fulfill those requirements (they may be available in the future, though). The latest quantum annealers commercially available are the D-Wave 2X, with 1152 and 3360 couplers (connections between adjacent qubits), and the D-Wave 2000Q, which has 2048 qubits and 6016 couplers. Their graphs are far from being complete, as every qubit in their grids are at most connected with six other qubits, as shown in Figure 15 (it may be less than that, as some qubits may be off after the last recalibration of the machine).

The importance of the completeness of the graph is problem dependant. Our QUBO instance, for example, has a complete connectivity graph between its variables. With an ideal quantum annealer we would have no problem, but with the D-Wave machine it is mandatory to transform our problem graph into an alternative graph that could be embedded into the Chimera graph. In other words, we have to solve an instance of the subgraph isomorphism problem, which happens to be in NP-complete. Even more, our problem may not be embeddable inside the D-Wave (for example, for the 1152 qubit Chimera grid of the D-Wave 2X, the largest complete graph that can be embedded into it is believed to be the $K_{33}$). Research on the subject of embedding a problem graph into D-Wave’s Chimera graph can be found in [11] and [12], where the concepts of embedding and parameter setting are explained.

There is a way to skip these limitations, by solving subinstances of our graph instead of the complete graph. For that, D-Wave released a graph partitioning open source library called qbsolv [8]. Its corresponding executable needs a certain kind of file format (called .qubo), to work. For example, if our QUBO instance is defined by

[TABLE]

where $x_{0},x_{1},x_{2}\in\{0,1\},$ its corresponding .qubo file is the one shown in Figure 13. The format is quite: the first line always starts with p qubo 0, followed by the number of variables, the number of nonzero diagonals, and the number of nonzero couplings. The output obtained by qbsolv for this file can be seen in Figure 14.

Please note that, in the previous paragraph, we have shown the solution that qbsolv obtains for a certain instance of the problem. This is because qbsolv has an auxiliary internal optimization solver based on the tabu search [15], which gives a solution for the subproblems and then unifies all the solutions into the complete one for all the variables. However, this classical method does not guarantee a global solution, and is of no help to us in the general case. We can decide if qbsolv tries to solve the problem with its tabu search, of if we prefer to connect to a D-Wave machine.

Finally, our QUBO subproblem is solved with the D-Wave 2X via one of the possible inputs allowed. The first one is shown in Figure 15, while the second one is again a .qubo file obtained internally with qbsolv. As part of the project Joint Research Unit Repsol-ITMATI (code file: IN853A 2014/03), we had the opportunity to try the D-Wave 2X machine based on the University of Southern California. However, due to the amount of subproblems needed to solve a proper instance of the Frobenius problem or the Apéry set, it was impractical to do so with the amount of time given and the current size of the graph of the D-Wave 2X machine.

4 A word on conclusions

Regarding the algorithms for the Apéry set and the Frobenius number, there are two aspects that need to be improved prior to completing a study of its feasibility and performance.

First, the current graph (i.e., Chimera graph) architecture of the available adiabatic quantum computers (i.e., the D-Wave machines) extremely obstruct the resolution of problems that have an almost full connectivity index between its variables, as in this case (the use of qbsolv is just a temporary workaround, or it should be as so).

And second, current adiabatic quantum computers do not guarantee global optimality, as opposed to the theoretical result deducted from the adiabatic theorem; in reality, the D-Wave just make a few runs of the process (instead of just one, as would be in the theoretical case) and, following a certain probability distribution, try to guarantee that the best of the solutions obtained by those runs is in fact the global solution to the problem. This solution, however, cannot be proven to be the global optimum (as global optimality is not known to be in the class NP), which makes useless our attempt to find those two combinatorial invariants via current adiabatic quantum computers.

In the future, with more reliable quantum annealers, however, this solution may prove to be effective and faster, but for now it is just a theoretical method.