Practical integer-to-binary mapping for quantum annealers

Abstract

Recent advancements in quantum annealing hardware and numerous studies in this area suggest that quantum annealers have the potential to be effective in solving unconstrained binary quadratic programming problems. Naturally, one may desire to expand the application domain of these machines to problems with general discrete variables. In this paper, we explore the possibility of employing quantum annealers to solve unconstrained quadratic programming problems over a bounded integer domain. We present an approach for encoding integer variables into binary ones, thereby representing unconstrained integer quadratic programming problems as unconstrained binary quadratic programming problems. To respect some of the limitations of the currently developed quantum annealers, we propose an integer encoding, named bounded-coefficient encoding, in which we limit the size of the coefficients that appear in the encoding. Furthermore, we propose an algorithm for finding the upper bound on the coefficients of the encoding using the precision of the machine and the coefficients of the original integer problem. We experimentally show that this approach is far more resilient to the noise of the quantum annealers compared to traditional approaches for the encoding of integers in base two. In addition, we perform time-to-solution analysis of various integer encoding strategies with respect to the size of integer programming problems and observe favorable performance from the bounded-coefficient encoding relative to that of the unary and binary encodings.

A (Integer-to-binary encoding)

In this section, we aim to find the upper bounds on the coefficients of the integer encoding when we reduce an unconstrained integer quadratic programming (UIQP) problem to an unconstrained binary quadratic programming (UBQP) problem in the $\{ 0,1\}$ domain. Similar to what we discussed earlier, we assume that our UIQP problem has the form

$$\begin{aligned} \begin{array}{lll} \min &{} x^T Q x + q^t x,\\ \text {s.t.} &{}x_i \in \{0, 1, 2, \ldots , \kappa ^{x_i} \} &{} \quad \text {for } i=\{1,2, \ldots , n\}, \\ \end{array} \end{aligned}$$

with Q being symmetric.

We aim to represent the above problem as f(y) by substituting $x= Cy$, where C is the encoding matrix we had earlier, i.e.,

$$\begin{aligned} C= \left[ \begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} c^{x_1}_1&{}\,\, c^{x_1}_2 &{}\,\, \dots &{}\,\, c^{x_1}_{d^{x_1}} &{}\,\, 0 &{}\,\, 0 &{}\,\, \dots &{}\,\, 0 &{}\,\, 0 &{}\,\, \dots &{}\,\, \dots &{}\,\, \dots &{}\,\, 0 &{}\,\, 0&{}\,\, 0 &{}\,\, \dots &{}\,\, 0 \\ 0 &{}\,\, 0 &{}\,\, \dots &{}\,\, 0&{}\,\, c^{x_2}_1&{}\,\, c^{x_2}_2 &{}\,\, \dots &{}\,\, c^{x_2}_{d^{x_2}} &{}\,\, 0&{}\,\, \dots &{}\,\, \dots &{}\,\, \dots &{}\,\, 0 &{}\,\, 0&{}\,\, 0 &{}\,\, \dots &{}\,\, 0 \\ \vdots &{}\,\, \vdots &{}\,\, \ddots &{}\,\, \vdots &{}\,\, \vdots &{}\,\, \vdots &{}\,\, \ddots &{}\,\, \vdots &{}\,\, \vdots &{}\,\, \ddots &{}\,\, \ddots &{}\,\, \ddots &{}\,\, \vdots &{}\,\, \vdots &{}\,\, \vdots &{}\,\, \ddots &{}\,\, \vdots \\ 0 &{}\,\, 0 &{}\,\, \dots &{}\,\, 0&{}\,\, 0 &{}\,\, 0 &{}\,\, \dots &{}\,\, 0&{}\,\, 0&{}\,\, \dots &{}\,\, \dots &{}\,\, \dots &{}\,\, 0 &{}\,\, c^{x_n}_1&{}\,\, c^{x_n}_2 &{}\,\, \dots &{}\,\, c^{x_n}_{d^{x_n}} \end{array}\right] , \end{aligned}$$

(33)

and

$$\begin{aligned} y= & {} \left[ \begin{array}{ccccccccccccccc} y^{x_1}_1&\,\, y^{x_1}_2&\,\, \dots&\,\, y^{x_1}_{d^{x_1}}&\,\,y^{x_2}_1&\,\, y^{x_2}_2&\,\, \dots&\,\, y^{x_2}_{d^{x_2}}&\,\, \dots&\,\, \dots&\,\, \dots&\,\, y^{x_n}_1&\,\, y^{x_n}_2&\,\, \dots&\,\, y^{x_n}_{d^{x_n}} \end{array}\right] ^t \nonumber \\&\in \mathbb B^{\sum _{i=1}^n d^{x_i} }. \end{aligned}$$

(34)

After the substitution for x, we get the following equivalent binary formulation:

$$\begin{aligned} f(y) = y^t (C^t Q C) y + (C^tq)^t y. \end{aligned}$$

(35)

Unlike the spin variables for which the diagonal of $C^t Q C$ became constant, the diagonal of $C^t Q C$ is added to the linear term in this scenario since $y_i^2 = y_i$ for $y_i \in \{0,1\}$. Alternatively, we can represent f(y) as

$$\begin{aligned} f_B = y^t Q_B y, \end{aligned}$$

(36)

where

(37)

The diagonal entries of $Q_B$ are the linear terms, i.e., coefficients of the variables $y_j^{x_i}$. We use linear and quadratic terms for binary model f(y) instead of local fields and couplers, respectively. We also refer to inequalities (10) and (11) as the ratio inequalities for linear and quadratic terms, respectively. Notice that in these ratios, minimum and maximum coefficients are measured in magnitude, so in our discussion that follows, we consider merely the magnitude of the coefficients.

Considering the fact that $y_k^{x_i}y_l^{x_i}=y_l^{x_i}y_k^{x_i}$ and $y_k^{x_i}y_l^{x_j} = y_l^{x_j}y_k^{x_i}$, the coefficients of linear and quadratic terms are listed below:

$$\begin{aligned} y_j^{x_i}&:&\left( Q_{ii} \left( c_j^{x_i}\right) ^2 + q_i c_j^{x_i} \right) , \end{aligned}$$

(38)

$$\begin{aligned} y_k^{x_i}y_l^{x_i}&:&2\left( Q_{ii} c_k^{x_i} c_l^{x_i}\right) \quad \text {for } \ k,l \in \{1, \ldots , d^{x_i}\} \ \text { and } \ k<l, \end{aligned}$$

(39)

$$\begin{aligned} y_k^{x_i}y_l^{x_j}&:&2\left( Q_{ij} c_k^{x_i} c_l^{x_j}\right) \quad \text {for }\ k \in \{1, \ldots , d^{x_i}\} ,\ l \in \{1, \ldots , d^{x_j}\}, \ \text { and } \nonumber \\&\quad i,j \in \{1, \ldots , n\} : i<j. \ \ \ \end{aligned}$$

(40)

The difference between this case and what was presented in Sect. 3 is that the coefficients of $y_j^{x_i}$ are no longer linear in $c_j^{x_i}$ (compare with (17)); therefore, the smallest coefficient may no longer occur at $c_j^{x_i}=1$. Notice that $Q_{ii} \left( c_j^{x_i}\right) ^2 + q_i c_j^{x_i}$ intersects 0 at $c_j^{x_i}= 0$ and $c_j^{x_i} = \frac{-q_i}{Q_{ii}}$; therefore, if $\frac{-q_i}{Q_{ii}} < 1$, then $Q_{ii} \left( c_j^{x_i}\right) ^2 + q_i c_j^{x_i}$ takes its minimum at $c_j^{x_i} =1$ and is increasing afterward. When this is the case for all $i = 1, \ldots , n$, slight modification of Algorithm 2 is sufficient to find the $\mu ^{x_i}$’s. In the algorithm, the modification is at the initialization step of $\mu ^{x_i}$’s. Letting

$$\begin{aligned} m_l = \min _{i} \left\{ |Q_{ii} + q_{i}| \right\} \ \text { and } \ m_c = \min _{i,j} \left\{ |Q_{ii}| , |Q_{ij}| \right\} , \end{aligned}$$

(41)

the condition that needs to be satisfied for the linear coefficients is

$$\begin{aligned} \frac{m_l}{ |Q_{ii}| (\mu ^{x_i})^2 + \text {sgn}(q_iQ_{ii}) |q_i| \mu ^{x_i}} \ge \epsilon _l. \end{aligned}$$

(42)

Combined with the condition

$$\begin{aligned} \mu ^{x_i} \le \sqrt{\frac{m_c}{|Q_{ii} |\epsilon _c} } , \end{aligned}$$

(43)

we need to initialize $\mu ^{x_i}$ as

$$\begin{aligned} \mu ^{x_i}= \left\lfloor \min \left\{ \tilde{\mu }^{x_i}, \sqrt{\frac{m_c}{|Q_{ii}| \epsilon _c} } \right\} \right\rfloor , \end{aligned}$$

(44)

where

$$\begin{aligned} \tilde{\mu }^{x_i} = \frac{-\text {sgn}(q_iQ_{ii})|q_i| + \sqrt{|q_i|^2 +4 |Q_{ii}|\frac{m_l}{\epsilon _l}}}{2|Q_{ii}|}. \end{aligned}$$

(45)

It is worth mentioning the special cases where $\text {sgn}(Q_{ii}) = \text {sgn}(q_i)$ (resulting in $\frac{-q_i}{Q_{ii}} < 0$), and $Q_{ii}=0$ or $q_i= 0$ belong to the above category, where $Q_{ii} \left( c_j^{x_i}\right) ^2 + q_i c_j^{x_i}$ attains its minimum at 1.

In the general cases where there exists an i such that $\frac{-q_i}{Q_{ii}}\ge 1$, satisfying the ratio condition on the linear terms is more complicated than what we discussed above. Not only might the minimum coefficient no longer occur at 1, but also the maximum could occur at either $\mu ^{x_i}$ or $\frac{-q_i}{2Q_{ii}}$, where the derivative of the function $h(c_j^{x_i}) = Q_{ii} \left( c_j^{x_i}\right) ^2 + q_i c_j^{x_i}$ is zero. Our approach for these general cases is to first derive $\mu ^{x_i}$’s that satisfy the ratio condition for the quadratic terms and then adjust them accordingly to meet the ratio constraint on the linear terms.

Similar to the previous case, the minimum coefficient on the quadratic terms is

$$\begin{aligned} m_c = \min _{i,j} \left\{ |Q_{ii}| , |Q_{ij}| \right\} , \end{aligned}$$

(46)

and the ratio condition on quadratic terms enforces

$$\begin{aligned} \frac{m_c}{ |Q_{ii}| (\mu ^{x_i})^2 } \ge \epsilon _c, \ \text {and} \ \frac{m_c}{ |Q_{ij}| \mu ^{x_i}\mu ^{x_j} } \ge \epsilon _c, \end{aligned}$$

(47)

or, equivalently,

$$\begin{aligned} \mu ^{x_i}\le & {} \sqrt{\frac{m_c}{|C_{ii}| \epsilon _c} } , \end{aligned}$$

(48)

$$\begin{aligned} \mu ^{x_i} \mu ^{x_j}\le & {} \frac{m_c}{|C_{ij}| \epsilon _c} . \end{aligned}$$

(49)

We may now initialize $\mu ^{x^i}$ as $\left\lfloor \sqrt{\frac{m_c}{|Q_{ii}| \epsilon _c}} \right\rfloor$ and use the loop of Algorithm 2 to satisfy (49), i.e., Algorithm 3.

After $\mu ^{x_i}$’s are calculated to satisfy the conditions (48) and (49), we need to check the ratio condition for the linear terms. Although some of the integer values $\{ 1, 2, \ldots , \mu ^{x_i} \}$ may not appear in our encoding, knowing which integers will appear prior to finding $\mu ^{x_i}$’s is not trivial. The algorithm presented below guarantees that the ratio condition on the linear terms holds if any of these integer values appear in the encoding. Let us categorize the indices based on where the minimum and maximum linear coefficients occur; we introduce the following sets of indices for this purpose:

$$\begin{aligned} \mathcal I_0^m= & {} \left\{ i: \text {sgn}(Q_{ii} ) \ne \text {sgn}( q_i ) \text { and } \frac{-q_i}{Q_{ii}}=1 \right\} , \nonumber \\ \mathcal I_1^m= & {} \left\{ i: \text {sgn}(Q_{ii} ) = \text {sgn}( q_i ) \right\} \nonumber \\&\cup \left\{ i: \text {sgn}(Q_{ii} ) \ne \text {sgn}( q_i ) \text { and } 0\le \frac{- q_i}{Q_{ii}}< 1 \right\} \nonumber \\&\cup \left\{ i: \text {sgn}(Q_{ii} ) \ne \text {sgn}( q_i ) \text { but } \mu ^{x_i} < \left\lfloor \frac{- q_i}{Q_{ii}} \right\rfloor \right\} \nonumber \\&\cup \left\{ i: \text {sgn}(Q_{ii} ) \ne \text {sgn}( q_i ) \text { but } \frac{- q_i}{Q_{ii}} \text { is an integer greater than 1}\right\} , \nonumber \\ \mathcal I_2^m= & {} \left\{ i: \text {sgn}(Q_{ii} ) \ne \text {sgn}( q_i ) \text { and } \frac{- q_i}{Q_{ii}} > 1 \text { and } \frac{- q_i}{Q_{ii}} \text { is not an integer }\text { and } \mu ^{x_i} \ge \left\lfloor \frac{- q_i}{Q_{ii}} \right\rfloor \right\} . \nonumber \\ \end{aligned}$$

(50)

and

$$\begin{aligned} \mathcal I_1^M= & {} \left\{ i: \text {sgn}(Q_{ii} ) = \text {sgn}( q_i ) \right\} \nonumber \\&\cup \left\{ i: \text {sgn}(Q_{ii} ) \ne \text {sgn}( q_i ) \text { and } 0\le \frac{- q_i}{2Q_{ii}} < \frac{1}{2} \right\} \nonumber \\&\cup \left\{ i: \text {sgn}(Q_{ii} ) \ne \text {sgn}( q_i ) \text { but } \mu ^{x_i} \le \text {round}\left( \frac{- q_i}{2Q_{ii}}\right) \right\} , \nonumber \\ \mathcal I_2^M= & {} \left\{ i: \text {sgn}(Q_{ii} ) \ne \text {sgn}( q_i ) \text { and } \frac{- q_i}{2Q_{ii}} \ge \frac{1}{2} \text { and } \mu ^{x_i} > \text {round}\left( \frac{- q_i}{2Q_{ii}}\right) \right\} . \qquad \end{aligned}$$

(51)

Note that $\mu ^{x_i}$’s returned by Algorithm 3 are integers. The sets with m and M superscripts are formed to facilitate computing minimum and maximum linear coefficients, respectively. For indices $i \in \mathcal I^m_0$, the minimum linear term happens at 2; for indices in $\mathcal I^m_1$, it happens at 1, and for indices in $\mathcal I^m_2$, it happens at one of the two integers closest to $\frac{-q_i}{Q_{ii}}$. Similarly, for the indices in $\mathcal I^M_1$, the maximum linear coefficient happens at $\mu ^{x_i}$, whereas for indices in $\mathcal I^M_2$, it happens at either $\mu ^{x_i}$ or at the closest integer to $\frac{-q_i}{2Q_{ii}}$, i.e., $\text {round}\left( \frac{- q_i}{2Q_{ii}}\right)$.

After categorizing the indices in the above sets, we may form the arrays $v^m$ and $v^M$, which represent the minimum and maximum coefficients of the linear terms for each variable, respectively:

$$\begin{aligned} v^m_i= & {} {\left\{ \begin{array}{ll} \left( |4Q_{ii} + 2q_i |, 2, i \right) &{} \text {if } i \in \mathcal I^m_0 \\ \left( |Q_{ii} + q_i |, 1, i \right) &{} \text {if } i \in \mathcal I^m_1 \\ \left( \min _{x \in \mathcal S_i } |Q_{ii} x^2 + q_i x| , {{\,\mathrm{arg\,min}\,}}_{x \in \mathcal S_i } |Q_{ii} x^2 + q_i x| , i \right) : \mathcal S_i = \left\{ x \in \left\{ \left\lfloor \frac{- q_i}{Q_{ii}} \right\rfloor , \left\lceil \frac{- q_i}{Q_{ii}} \right\rceil \right\} : x\le \mu ^{x_i} \right\} &{} \text {if } i \in \mathcal I^m_2 \end{array}\right. } \nonumber \\ \end{aligned}$$

(52)

$$\begin{aligned} v^M_i= & {} {\left\{ \begin{array}{ll} \left( |Q_{ii} (\mu ^{x_i})^2 + q_i \mu ^{x_i} |, \mu ^{x_i}, i \right) &{} \text {if } i \in \mathcal I^M_1 \\ \left( \max _{x \in \mathcal S_i } |Q_{ii} x^2 + q_i x| , {{\,\mathrm{arg\,max}\,}}_{x \in \mathcal S_i} |Q_{ii} x^2 + q_i x| , i\right) : \mathcal S_i = \left\{ \text {round}\left( \frac{- q_i}{2Q_{ii}}\right) , \mu ^{x_i} \right\} &{} \text {if } i \in \mathcal I^M_2 \end{array}\right. } \nonumber \\ \end{aligned}$$

(53)

We then sort entries of $v^m$ (on the first entry) in increasing order; assume it results in vector $\bar{v}^m$, so $\bar{v}_1^m \le \bar{v}_2^m \le \cdots \le \bar{v}_n^m$ and sort $v^M$ in decreasing order; that is, $\bar{v}^M$ such that $\bar{v}_1^M \ge \bar{v}_2^M \ge \cdots \ge \bar{v}_n^M$. We check whether $\frac{\bar{v}_1^m}{\bar{v}_1^M} \ge \epsilon _l$. If this inequality holds, then we have obtained our set of $\mu ^{x_i}$’s; otherwise, we have the following options for improvement:

• $\bar{v}^m_1$ happens at an $i \in \mathcal I_2^m$, where we can either change $\mu ^{x_i}$ to increase the minimum coefficient, or change $\mu ^{x_i}$ to decrease the maximum coefficient;

• $\bar{v}^m_1$ happens at $i \in \mathcal I_1^m \cup \mathcal I_0^m$, in which case we can only change $\mu ^{x_i}$ such that the maximum coefficient decreases.

In the first scenario where we have the option to both increase the minimum or decrease the maximum coefficient, we use a greedy approach to make the decision. In other words, if the minimum coefficient is improved, it will be $v^m_{2,1}$ in the next iterate; similarly, the maximum coefficient will $\bar{v}^M_{2,1}$, if updated. In our approach, having the interval $[\bar{v}^m_{1,1} , \bar{v}^M_{1,1}]$ for the coefficients, we wish to update it to either $[\bar{v}^m_{2,1} , \bar{v}^M_{1,1}]$ or $[\bar{v}^m_{1,1} , \bar{v}^M_{2,1}]$. These two intervals will be $\left[ \frac{\bar{v}^m_{2,1}}{\bar{v}^M_{1,1}} , 1\right]$ or $\left[ \frac{\bar{v}^m_{1,1}}{ \bar{v}^M_{2,1}}, 1\right] ]$, respectively, after the rescaling. We choose the option that gives us better lower bound, i.e., if

$$\begin{aligned} \frac{\bar{v}^m_{2,1}}{\bar{v}^M_{1,1}}> \frac{v^m_{1,1}}{ \bar{v}^M_{2,1}} \ \equiv \ \bar{v}^m_{2,1} \bar{v}^M_{2,1} > \bar{v}^m_{1,1} \bar{v}^M_{1,1}, \end{aligned}$$

we attempt to improve the lower bound, thus decreasing $\mu ^{x_i}$ for $i = \bar{v}^m_{1,3}$. A formal presentation of what we have discussed is summarized in Algorithm 4.