Practical Integer-to-Binary Mapping for Quantum Annealers

TL;DR

This paper introduces a new integer-to-binary encoding method for quantum annealers that enhances noise resilience and expands their application to integer quadratic programming problems.

Contribution

We propose a bounded-coefficient encoding for integer variables in quantum annealing, improving noise robustness and applicability to broader optimization problems.

Findings

The bounded-coefficient encoding is more noise-resilient than traditional binary encoding.

The proposed algorithm effectively determines coefficient bounds based on machine precision.

Experimental results demonstrate improved performance on quantum annealers.

Abstract

Recent advancements in quantum annealing hardware and numerous studies in this area suggests 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. Finally, 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.

Tables

Table 1. Table 1: Resilience with Bounded-Coefficient Encoding.

	0.001	0.002	0.003	0.004	0.005	0.006	0.007	0.008	0.009	0.01
convex-1	1.0	1.0	0.6	0.4	0.8	0.5	0.2	0.1	0.1	0.0
convex-2	0.7	0.5	0.4	0.1	0.2	0.0	0.0	0.2	0.0	0.0
convex-3	1.0	0.7	0.6	0.5	0.2	0.1	0.0	0.1	0.0	0.0
convex-4	1.0	0.9	0.6	0.3	0.2	0.1	0.1	0.0	0.1	0.0
convex-5	1.0	1.0	0.5	0.6	0.3	0.4	0.2	0.2	0.4	0.1
$(U_2,U_{200})$	1.0	1.0	1.0	1.0	0.9	1.0	0.9	0.8	1.0	0.8
$(U_5,U_{200})$	1.0	0.9	0.9	0.9	0.8	0.4	0.7	0.8	0.5	0.6
$(U_{10},U_0)$	1.0	0.8	0.4	0.3	0.4	0.1	0.4	0.2	0.1	0.2
$(U_5,U_{10})$	1.0	1.0	1.0	0.9	0.7	0.9	0.4	0.4	0.6	0.7
$(U_5,U_{100})$	1.0	1.0	1.0	1.0	1.0	1.0	1.0	0.7	0.7	0.6
average	0.97	0.88	0.7	0.6	0.55	0.45	0.39	0.35	0.35	0.3

Table 2. Table 2: Resilience with Binary Encoding.

	0.001	0.002	0.003	0.004	0.005	0.006	0.007	0.008	0.009	0.01
convex-1	0.1	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
convex-2	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
convex-3	0.1	0.0	0.0	0.0	0.0	0.1	0.0	0.0	0.0	0.0
convex-4	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
convex-5	0.1	0.0	0.0	0.1	0.0	0.0	0.0	0.0	0.0	0.0
$(U_2,U_{200})$	0.6	0.4	0.2	0.1	0.3	0.1	0.3	0.1	0.4	0.2
$(U_5,U_{200})$	0.5	0.7	0.3	0.4	0.5	0.3	0.2	0.5	0.2	0.2
$(U_{10},U_0)$	0.4	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
$(U_5,U_{10})$	0.3	0.1	0.1	0.1	0.1	0.0	0.1	0.1	0.1	0.0
$(U_5,U_{100})$	0.5	0.3	0.1	0.4	0.3	0.4	0.3	0.1	0.4	0.1
average	0.26	0.15	0.07	0.11	0.12	0.09	0.09	0.08	0.11	0.05

Equations

$$\underset{s\in \{\pm 1{\}}^n}{\min }{s}^{t}JS+h^{t}s,$$

$$\begin{array}[]{lll}\min&x^{T}Qx+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}$$

$$x_i=\sum _{j=1}^{d^{x_i}}{c}_j^{x_i}{y}_j^{x_i}=(c^{x_i}{)}^t{y}^{x_i},$$

$$y_j=\frac{1}{2}\left(s_j+1\right),$$

$$x_i=\frac{1}{2}\left(\sum _{j=1}^{d^{x_i}}{c}_j^{x_i}+\sum _{j=1}^{d^{x_i}}{c}_j^{x_i}{s}_j^{x_i}\right),$$

$$x=c^{t}y,$$

$$c^x=\left[2^0,2^1,\dots ,2^{\lfloor \log ({\kappa }^x)\rfloor -1},{\kappa }^x-\sum _{i=1}^{\lfloor \log ({\kappa }^x)\rfloor }{2}^{i-1}\right]$$

$$c_i^x=\{\begin{array}{ll}{2}^{i-1}& \text{ for }i=1,\dots ,\rho \\ {\mu }^x& \text{ for }i=\rho +1,\dots ,\rho +\eta \\ \nu -\eta {\mu }^x& \text{ for }i=\rho +\eta +1\text{ if }\nu -\eta {\mu }^x\ne 0\end{array}$$

$$d^x=\{\begin{array}{ll}\rho +\eta +1& \text{ if }\nu -\eta \mu \ne 0\\ \rho +\eta & \text{ otherwise }\end{array}.$$

$$\begin{array}[]{rcll}\sum_{j=1}^{\rho}c_{j}y_{j}&\geq&(2^{\rho}-\mu^{x})y_{i}&\text{for }i=\rho+1,\\ y_{i}&\geq&y_{i+1}&\text{for }i=\rho+1,\ldots,d^{x}-1,\\ \sum_{j=1}^{\rho}c_{j}y_{j}&\geq&(2^{\rho}-c_{d^{x}})y_{d^{x}}.\end{array}$$

$$x=\frac{1}{2}\left({\kappa }^x+\sum _{j=1}^{d^x}{c}_j^x{s}_j^x\right).$$

$$\frac{{\min }_i|h_i|}{{\max }_i|h_i|}\ge {ϵ}_l.$$

$$\frac{{\min }_{i,j}|J_{ij}|}{{\max }_{i,j}|J_{ij}|}\ge {ϵ}_c.$$

$$\underset{i,j}{\min }\left|h_i-h_j\right|\ge {ϵ}_l,$$

$$\underset{\begin{array}{c}i,j\\ l,k\end{array}}{\min }\left|J_{ij}-J_{lk}\right|\ge {ϵ}_c,$$

$$f(x)=x^{t}Qx+q^{t}x,$$

$$C=\text{setcounter}MaxMatrixCols17\left[\begin{array}{ccccccccccccccccc}{c}_1^{x_1}& c_2^{x_1}& \dots & c_{d^{x_1}}^{x_1}& 0& 0& \dots & 0& 0& \dots & \dots & \dots & 0& 0& 0& \dots & 0\\ 0& 0& \dots & 0& c_1^{x_2}& c_2^{x_2}& \dots & c_{d^{x_2}}^{x_2}& 0& \dots & \dots & \dots & 0& 0& 0& \dots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & \ddots & \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 0& 0& 0& \dots & 0& 0& \dots & \dots & \dots & 0& c_1^{x_n}& c_2^{x_n}& \dots & c_{d^{x_n}}^{x_n}\end{array}\right],$$

$$s=\text{setcounter}MaxMatrixCols15{\left[\begin{array}{ccccccccccccccc}{s}_1^{x_1}& s_2^{x_1}& \dots & s_{d^{x_1}}^{x_1}& s_1^{x_2}& s_2^{x_2}& \dots & s_{d^{x_2}}^{x_2}& \dots & \dots & \dots & s_1^{x_n}& s_2^{x_n}& \dots & s_{d^{x_n}}^{x_n}\end{array}\right]}^t.$$

$$x=\frac{1}{2}\left(\kappa +Cs\right),$$

$$f(s)=\frac{1}{4}{s}^t\left(C^{t}QC-\Delta (\delta (C^{t}QC))\right)s$$

$$s_j^{x_i}$$

$$s_k^{x_i}{s}_l^{x_i}$$

$$s_k^{x_i}{s}_l^{x_j}$$

$$m_l=\underset{i}{\min }\left\{\left|{\left[Q\kappa +q\right]}_i\right|\right\},$$

$$m_c=\underset{i,j}{\min }\left\{\left|Q_{ii}\right|,\left|Q_{ij}\right|\right\}.$$

$$\frac{m_l}{\left|{\left[Q\kappa +q\right]}_i\right|{\mu }^{x_i}}\ge {ϵ}_l.$$

$$\frac{m_c}{|Q_{ii}|({\mu }^{x_i}{)}^2}\ge {ϵ}_c,\text{and}\frac{m_c}{|Q_{ij}|{\mu }^{x_i}{\mu }^{x_j}}\ge {ϵ}_c.$$

$${\mu }^{x_i}$$

$${\mu }^{x_i}{\mu }^{x_j}$$

$$\mu =\underset{i}{\min }\{{\mu }^{x_i}\}.$$

$$f(s)=\frac{1}{2}{q}^t\left(\kappa +Cs\right)=\frac{1}{2}{q}^t\kappa +\frac{1}{2}{\left(C^{t}q\right)}^{t}s.$$

Full text

1 Introduction

Adiabatic quantum computation (AQC) has been proposed as a successful technique for solving certain classes of optimization problems (see albash2016 ; farhi01 ; McGeoch2013 ; optimization_AQO ). In particular, quantum annealers (QA) such as the ones manufactured by D-Wave Systems Inc. DwaveNature approximate the adiabatic evolution in the presence of various sources of noise and finite temperature to solve Ising models of the form

[TABLE]

where diagonal entries of $J$ are zeroes, and the nonzero entries of $J$ create a subgraph of a sparse graph called a Chimera graph chimera . Many researchers have explored the applicability of quantum annealing to more-general optimization problems. Solving Ising models without the Chimera graph structure on a QA is discussed in AidanRoy2014 . Degree reduction techniques, such as the one described in Ishikawa2011 , enable us to solve unconstrained binary polynomial programming problems of higher order using a QA. Several studies Venturelli2015 ; Gili ; Dominic2015 ; Zick have employed quantum annealing for solving real-world applications, and karimi_GD discusses the possibility of using QAs for solving constrained problems.

Many optimization problems, on the other hand, involve integer-valued variables beyond binary values. Examples of such problems could be the number of vehicles or products travelling along each route of a supply chain (see supplyChain ), or the number of traded assets in portfolio optimization (see portfolioOpt ). This paper focuses on solving quadratic problems on bounded integer domains, alternatively referred to as unconstrained integer quadratic programming (UIQP) problems, i.e., problems of the following form:

[TABLE]

where $\kappa^{x_{i}}\in\mathbb{Z}_{+}$ is an upper bound on $x_{i}$; here, $\mathbb{Z}_{+}$ denotes the set of non-negative integers. Note that in problems where $x_{i}\in\{\alpha,\alpha+1,\alpha+2,\ldots,\alpha+\kappa^{x_{i}}\}$, we may shift $x_{i}$ and substitute it with $x_{i}-\alpha$; hence, without loss of generality, we may assume that the domain of $x_{i}$ starts at 0. It is worth mentioning that, although current QAs are limited to representing Ising models, an alternative approach to encoding bounded integer variables in terms of spin variables is to consider physical implementations of quantum processors as representing Potts-Ising models AQC_qudits ; katzgraber_qudits . Our approach is different and is based on reformulating a UIQP problem as an Ising model. To this end, we represent each integer variable as a linear combination of several binary variables, i.e.,

[TABLE]

where $c_{j}^{x_{i}}\in\mathbb{Z}_{+}$, $y_{j}^{x_{i}}\in\{0,1\}$ for $j\in\{1,\ldots,d^{x_{i}}\}$, and the superscript $x_{i}$ is used for clarity to denote that $c^{x_{i}}$ and $y^{x_{i}}$ correspond to the encoding of variable $x_{i}$. This representation is referred to as integer encoding. Some of the well-known integer encodings are binary and unary encodings in which $c_{j}^{x_{i}}=2^{j-1}$ and $c_{j}^{x_{i}}=1$, respectively. The width of an integer encoding, denoted by $d^{x_{i}}$ in (1), refers to the number of binary variables required for encoding integer variable $x_{i}$. The width of binary and unary encodings are $\left\lfloor\log_{2}(\kappa^{x_{i}})\right\rfloor+1$ and $\kappa^{x_{i}}$, respectively. Since a binary variable $y_{j}$ could be represented with a spin variable $s_{j}$ via the affine transformation

[TABLE]

integer variable $x_{i}$ could be encoded into several spin variables as

[TABLE]

which enables us to successfully represent a UIQP as an Ising model.

After derivation of an Ising model equivalent to (1), heuristic methods including quantum annealing could be employed to find the ground state of the problem. The performance of the heuristic methods highly depends on the energy landscape of the problem. In particular, landscapes with tall barriers are challenging for the hill-climbing heuristics, and having wide barriers could impact the performance of methods like quantum annealing that could potentially benefit from quantum tunnelling (see katzgraber2015 and references therein for more details). Having coefficients that are different in orders of magnitude is undesirable, as they could create energy landscapes with tall and wide barriers. Moreover, D-Wave quantum annealers have only a low precision of approximately $10^{-2}$ for couplings’ strengths and local fields’ biases (i.e., entries of $J$ and $h$ in (Ising)), scaling $J$ and $h$ to a very limited range, (e.g., $[-2,2]$). Several sources of noise such as thermal excitations and control errors contribute to the low precision of the machines; in albash2014 , Albash et al. propose a noise model for D-Wave devices that includes the control noise of the local field and couplings of the chip. Zhu et al. zoshk2016 show that increasing the classical energy gap beyond the intrinsic noise level of the machine can improve the success of the D-Wave Two QA. In zoshk2016 , resilience is defined as a metric for measuring the resistance of a problem to noise. More precisely, resilience is the probability that the ground state does not change under random field fluctuations found on the chip. In this paper, we argue that restricting the range of the coefficients of a problem could improve resilience; and our results, presented in Section 4 of this paper, support this argument

Studies in biasedGSsampling ; GSstat suggest that QAs could be biased in finding degenerate ground states; in other words, sampling with QAs does not return different degenerate solutions with uniform probability. Moreover, degenerate ground states are easier to reach; therefore, benchmarks created for QAs tend to avoid degeneracy and have a unique ground state katzgraber2015 . While there is a need to rigorously study whether highly degenerate low-energy excited states could impact reaching the ground state, avoiding degeneracy could be an alternative. Recall that in integer encoding, an integer variable, $x$, is substituted with $c^{t}y$, where we assume $c\in\mathbb{Z}_{+}^{d}$ and $y\in\mathbb{B}^{d}$. Any binary vector $y$ returns an integer value. If the total number of binary combinations, i.e., $2^{d}$, is larger than the summation of entries of $c$, i.e., $\sum_{i=1}^{d}c_{i}$, some of the integers occur at more than one binary combination. The code-word of an integer value $\chi$, for $\chi\in\{1,2,\ldots,\kappa\}$, refers to the cardinality of the set $\{y:c^{t}y=\chi\}$. Code-words more than one is what we refer to as redundancy in this paper. Under the assumption that the QA is inclined to observe integers with higher code-words, it is ideal in an integer encoding to have a unique code-word for all integers. Binary encoding, for example, keeps the code-word of each integer uniformly at one, but the code-word in unary encoding is highly variant and the redundancy has its peak at $\left\lceil\frac{\kappa}{2}\right\rceil$ with word count $\kappa\choose\lceil\frac{\kappa}{2}\rceil$.

This work presents an integer encoding with minimal width that creates noise-resilient Ising models. Since one or several qubits is assigned to each spin variable and the number of qubits on a chip is limited, it is desirable to keep the size of the Ising model, and thus the width of the encoding, as small as possible. Binary encoding has the minimum width and uniform code-word; however, the coefficients $c_{j}$’s in the binary encoding could get arbitrarily large and result in an error-prone (noisy, non-resilient) Ising model. On the other hand, unary encoding does not expand the range of the coefficients of the problem, but it has an exceedingly large degree and redundancy.

The notation used in this paper is as follows. Generally, we reserve the upper-case letters, lower-case letters, and Greek alphabet for matrices, vectors, and scalars, respectively, with only a few exceptions; these exceptions, however, should be clear from the context. An all-ones vector and the identity matrix are denoted by $e$ and $I$, respectively. Entries of a vector or matrix are differentiated with subscript indices. $\lceil\cdot\rceil$ and $\lfloor\cdot\rfloor$ denote the ceiling and floor of a number, respectively. The $\log(\cdot)$ function is in base two; and $\text{sgn}(\alpha)$ returns the sign of $\alpha$. The function $\delta(A):\mathbb{R}^{n\times n}\rightarrow\mathbb{R}^{n}$ returns a vector with the diagonal entries of $A$; and the function $\Delta(a):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n\times n}$ returns a diagonal matrix with entries of $a$ on the diagonal.

This paper is organized as follows. In the next section, we present bounded-coefficient encoding under the assumption that an upper bound on the coefficients of the integer encoding is available. In Section 3, we present how, given the precision, we may find an upper bound on the coefficients of the integer encoding. We test our algorithm by comparing binary and bounded-coefficient encodings in Section 4. Finally, we conclude our discussion in Section 5.

2 Bounded-Coefficient Encoding

Let $x$ be an integer variable with upper bound $\kappa^{x}$, and $\mu^{x}$ be an upper bound on the coefficients of the integer encoding. In other words, we represent $x$ as

[TABLE]

where $y\in\mathbb{B}^{d^{x}}$, $c\in\mathbb{Z}_{+}^{d^{x}}$, and $c_{i}\leq\mu^{x}$ for $i=1,\ldots,d^{x}$. The encoding we propose in this paper, which we refer to as bounded-coefficient encoding, is summarized in the following algorithm:

Notice that if $\kappa^{x}<2^{\lfloor\log\mu^{x}\rfloor+1}$, then the binary encoding respects the upper bound on the coefficients. These cases are captured in (5), and the width of encoding in these instances is $\lfloor\log(\kappa^{x})\rfloor+1$.

When $\kappa^{x}\geq 2^{\lfloor\log\mu^{x}\rfloor+1}$, the encoding is derived using (6) and the width of the bounded-coefficient encoding is

[TABLE]

Here are a few demonstrative examples of bounded-coefficient encoding:

• •

$\kappa=12$ and $\mu=8$, the bounded-coefficient encoding is $c=\left[1,\ 2,\ 4,\ 5\right]^{t};$

• •

$\kappa=20$ and $\mu=6$, the bounded-coefficient encoding is $c=\left[1,\ 2,\ 4,\ 6,\ 6,\ 1\right]^{t}.$

We assume $\kappa^{x}\gg 2^{\lfloor\log\mu^{x}\rfloor+1}$ when we refer to the bounded-coefficient encoding, and our general statements are focused on these cases where we have multiple coefficients of size $\mu^{x}$. Propositions 1 and 2, also, refer to these cases described by (6); these propositions hold true for binary encoding; hence, they hold for the cases that are derived with (5), such as our first example above.

Definition 1

An encoding $c$ is called $\kappa$-complete if it can encode only and all integers $\{0,1,\ldots,\kappa\}$.

Proposition 1

The bounded-coefficient encoding, generated by Algorithm 1, is $\kappa^{x}$-complete.

By the first part of our encoding, i.e., $1,\ldots,2^{\rho-1}$, we can encode all integers $\{0,1,\ldots,\mu^{x}\}$. By adding $\mu^{x}$ to those binary combinations, we can generate all integers from $\mu^{x}$ to $2\mu^{x}$; if we add two $\mu^{x}$ factors, we get integers $2\mu^{x}$ to $3\mu^{x}$, and so on. Finally, we are guaranteed that we get integers $\eta\mu^{x}$ to $\kappa^{x}$ because $\kappa^{x}-\eta\mu^{x}<\sum_{i=1}^{\rho}2^{i-1}-\mu^{x}$ by the fact that $\nu-\eta\mu^{x}<\mu^{x}$.

Definition 2

A sub-encoding $\tilde{c}$ of an encoding $c$ is a choice of entries of $c$, with a fixed ordering arranged in a column vector $\tilde{c}$.

Proposition 2

Bounded-coefficient encoding, generated by Algorithm 1, has minimum width among all $\kappa^{x}$-complete integer encodings of the form (1) that respect an upper bound on their coefficients.

Proof

Our arguments rely on the optimality of binary encoding, or, more rigorously, the following points:

1. (a)

binary encoding has minimum width; in other words, the binary encoding of values $\{0,1,\ldots,\mu^{x}\}$ has width $d_{B}=\lfloor\log\mu^{x}\rfloor+1$, and no other $\mu^{x}$-complete integer encoding could have a width less than $d_{B}$. 2. (b)

binary encoding of width $d_{B}$, i.e., $\left(2^{0},\ 2^{1},\ 2^{2},\ \ldots,\ 2^{d_{B}-1}\right)$, encodes integers of maximum value $2^{d_{B}}-1$; i.e., $\sum_{i=0}^{d_{B}-1}2^{i}=2^{d_{B}}-1$. Notice that by (a), this is the largest number that can have an encoding of width $d_{B}$.

It is worth mentioning that if $\kappa^{x}<2^{\lfloor\log\mu^{x}\rfloor+1}$, then the bounded-coefficient encoding has width $\lfloor\log\kappa^{x}\rfloor+1$, and by (a) has minimum width.

For more-general cases, our proof is by contradiction. Suppose $c$ is the bounded-coefficient encoding derived by Algorithm 1, and $f$ is a bounded-coefficient encoding of smaller width; i.e., $d_{c}>d_{f}$. Let $f_{\mu}$ be the minimal sub-encoding of $f$ that encodes $\{0,1,\ldots,\mu^{x}\}$, and let $d_{f_{\mu}}\leq d_{f}$ be the width of $f_{\mu}$. Similarly, $c_{\mu}$ in the bounded-coefficient encoding derived by Algorithm 1 refers to the sub-encoding of $c$ required to encode integers less than or equal to $\mu^{x}$. As implied by Algorithm 1 and (6), $d_{c_{\mu}}=\rho$. By (a), we conclude $d_{f_{\mu}}\geq\rho$. We can, now, consider two cases: $d_{f_{\mu}}=\rho$ and $d_{f_{\mu}}>\rho$.

Consider the case where $d_{f_{\mu}}=\rho$. By (a) and (b), $\sum_{\alpha\in c_{\mu}}\alpha\geq\sum_{\alpha\in f_{\mu}}\alpha$; therefore, $\kappa^{x}-\sum_{\alpha\in f_{\mu}}\alpha\geq\nu$, where $\nu$ is as defined in Algorithm 1. Since all of the coefficients need to be bounded by $\mu^{x}$, and by the fact that $\frac{\kappa^{x}-\sum_{\alpha\in f_{\mu}^{x}}\alpha}{\mu^{x}}\geq\frac{\nu}{\mu^{x}}$, we conclude that having $d_{f}<d_{c}$ is impossible.

Consider now the case that $d_{f_{\mu}}\geq\rho+1$. As all of our coefficients are bounded by $\mu^{x}$, $\sum_{\alpha\in f_{\mu}}\alpha<2\mu^{x}$; otherwise, this contradicts that $f_{\mu}$ is the minimal sub-encoding required to encode integers less than or equal to $\mu^{x}$. While $\sum_{\alpha\in f_{\mu}}\alpha<2\mu$, $\sum_{\alpha\in c_{\mu}}\alpha+\mu\geq 2\mu$. Similar to the previous case, we conclude that $\kappa^{x}-\sum_{\alpha\in f_{\mu}}\alpha>\kappa^{x}-\sum_{\alpha\in c_{\mu}}\alpha+\mu^{x}$; hence, $d_{f}<d_{c}$ would not be possible. ∎

As mentioned earlier, a potentially advantageous property of an integer encoding is uniform or low-variant redundancy. In general, the bounded-coefficient encoding has redundancy because we have more than one coefficient with a value of $\mu$. It is, however, worth mentioning that forcing the following constraints on the binary variables of (4) makes the word count of each integer value unique:

[TABLE]

Note that $(2^{\rho}-\mu^{x})\leq\mu^{x}$; however, $(2^{\rho}-c_{d^{x}})$ may not necessarily be less than or equal to $\mu^{x}$. If we wish to limit the coefficients of our constraints to $\mu^{x}$, we can substitute $(2^{\rho}-c_{d^{x}})y_{d^{x}}$ in the right-hand side of the last inequality with $\mu^{x}y_{d^{x}-1}y_{d^{x}}+(2^{\rho}-\mu^{x}-c_{d^{x}})y_{d^{x}}$ whenever $2^{\rho}-c_{d^{x}}>\mu^{x}$. Handling constraints, despite being possible, is non-trivial for QA; it may introduce exceedingly large coefficients (see discussion in karimi_GD ), and could be contradictory to the purpose of this work.

As we presented earlier in (2) and (3), we may directly encode integer variables into spin variables. Moreover, by the fact that in the bounded-coefficient encoding $\sum_{j=1}^{d^{x}}c_{j}=\kappa^{x}$, the encoding of each integer variable into spin variables would be

[TABLE]

In the next section, where we try to find the upper bound on the coefficients of the encoding, i.e., $\mu^{x}$, we use the integer-to-spin transformation.

3 Finding the Upper Bound on the Coefficients of the Encoding

We explained earlier, in Section 1, that the QA scales the coefficients of an Ising model to a limited range and has low precision. We also mentioned that resilience to noise could be helpful in overcoming the issue of precision. We propose that restricting the range of the coefficients of the problem could improve resilience. To this end, we wish to find upper bounds on the coefficients of the encoding, i.e., $\mu$, such that after the encoding of the integer problem (1) to an Ising model (Ising), the ratio of the smallest local field in magnitude to the largest one exceeds a threshold, i.e.,

[TABLE]

Similarly, we wish to bound the ratio of the smallest coupler to the largest one:

[TABLE]

Alternatively, one may desire to have the local fields biases and couplings’ strengths well-separated to avoid the effect of noise, or in other words, having

[TABLE]

and

[TABLE]

could be desirable. In this paper, however, we focus merely on inequalities (10) and (11).

Suppose we have a quadratic function

[TABLE]

where $Q$ is symmetric and the domain of $f(x)$ is $x=\left[x_{1},x_{2},\ldots,x_{n}\right]^{t}$ for $x_{i}\in\{0,\ 1,\ \ldots,\ \kappa^{x_{i}}\}$. Note that this is the function which we aim to minimize in (1). Also, note that since $x_{i}x_{j}=x_{j}x_{i}$ for two integers $x_{i}$ and $x_{j}$, we may substitute $Q$ with $\frac{1}{2}(Q+Q^{t})$ when $Q$ is not symmetric. Of course, everything presented here is under the assumption that the coefficients of the problem (before any encoding) respect inequalities (10) and (11), i.e., unary encoding satisfies them. Let us define the encoding matrix denoted by $C$ as

[TABLE]

and the vector of spin variables as

[TABLE]

Using (9), we conclude that the encoding of integer variables, $x$, into spin variables, $s$, refers to the following substitution:

[TABLE]

where $\kappa=\left[\kappa^{x_{1}},\ \kappa^{x_{2}},\ \ldots,\ \kappa^{x_{n}}\right]^{t}$ is the vector of all of the upper bounds on the integer variables. Using this substitution in (12), we get

[TABLE]

We are excluding the diagonal entries of $C^{t}QC$ in the quadratic term because they correspond to the square of spin variables, i.e., $s_{i}^{2}$, and $s_{i}^{2}$ is a constant equal to one. The local field bias corresponding to each spin variable and the coupling strength corresponding to each pair of spin variables, derived in (16), is summarized below:

[TABLE]

where $\left[Q\kappa+q\right]_{i}$ in the first line denote the $i$-th entry of vector $Q\kappa+q$, as expected, and (18) and (19) incorporate the fact that $s_{k}^{x_{i}}s_{l}^{x_{i}}=s_{l}^{x_{i}}s_{k}^{x_{i}}$ and $s_{k}^{x_{i}}s_{l}^{x_{j}}=s_{l}^{x_{j}}s_{k}^{x_{i}}$.

Since all of the local fields biases and couplings’ strengths share the $\frac{1}{2}$ factor and after rescaling our concern is their ratio, we will drop the $\frac{1}{2}$ factor for our calculations in the rest of this section. From the previous section and the derivation of the bounded-coefficient encoding, we conclude that the smallest coefficient used in the encoding is 1; hence, the minimum absolute value of the local fields biases is

[TABLE]

and the minimum absolute value of couplings’ strengths is

[TABLE]

The maximum absolute value of the local fields biases among all spin variables corresponding to an integer variable $x_{i}$ occurs at $\mu^{x_{i}}$—to be found—for each $i\in\{1,\ldots,n\}$. Therefore, by (10), we wish to have:

[TABLE]

Similarly, (11) enforces the following conditions for the couplers:

[TABLE]

Using (20) and (21), we wish to obtain $\mu^{x_{i}}$’s that satisfy the following set of inequalities:

[TABLE]

We may solve a feasibility problem to find a solution to the above set of inequalities, i.e., an optimization or auxiliary objective function along with these inequalities. A well-justified objective function could be $\max\ \min_{i}\{\mu^{x_{i}}\}$. Solving such problems is normally costly because (24) is non-convex and $\mu^{x_{i}}$’s are required to be integers. Alternatively, we present an algorithm below that heuristically finds $\mu^{x_{i}}$’s.

Note that any set of $\mu^{x_{i}}$’s satisfying the above inequalities will guarantee that (10) and (11) hold. Our proposed algorithm for finding $\mu^{x_{i}}$’s first initializes $\mu^{x_{i}}$’s using inequalities (22) and (23). If (24) is satisfied for all $i$ and $j$ at this step, it terminates. Otherwise, it greedily decrease $\mu^{x_{i}}$ or $\mu^{x_{j}}$ for an $i$ and $j$ pair for which the failure happened. In this process, we take $\frac{\kappa^{x_{i}}}{\mu^{x_{i}}}$ as an estimate on the width of the encoding, so when we want to decrease either $\mu^{x_{i}}$ or $\mu^{x_{j}}$, we choose the one that gives a lower combined width. See Algorithm 2 for a formal presentation of this algorithm.

It is worth mentioning that if we wish to keep $\mu^{x_{i}}$’s for all integer variables equal, we may use the value

[TABLE]

Moreover, if in (12) the quadratic terms were missing, i.e., $f(x)=\ q^{t}x$, we would get

[TABLE]

Ignoring the $\frac{1}{2}$ factor (since it is common for all spin variables), we get $m_{l}=\min_{i}\ \{\left|q_{i}\right|\}$, and our ratio condition reduces to

[TABLE]

therefore,

[TABLE]

Finally, we would like to point out that unlike QA, in several heuristic methods, variables tend to be binary, i.e., in the $\{0,1\}$ domain referred to as binary. In Appendix A, we present a modification of the method in this section that could be employed in this case.

4 Numerical Experiment

In this section, we test the bounded-coefficient encoding. We compare binary and bounded-coefficient encodings on ten randomly generated instances. To diversify our instances, however, we use several procedures to generate them. Let us define the set $U_{\alpha}=\{0,\ \pm 1,\ \pm 2,\ldots,\ \pm\alpha\}$. In all of our instances, we have five integer variables, i.e., $n=5$; the upper-bound on all of the integer variables is 50, i.e., $\kappa^{x_{i}}=50$ for $i=1,\ldots,5$; and the matrix $Q$ in the quadratic term has a sparsity around 50%. In half of our instances, the generated $Q$ is positive definite, and in the rest it is indefinite; these two categories are referred to as convex and non-convex instances, respectively.

For the convex instances, entries of matrix $Q$ are initially drawn from $U_{2}$; then, $\lambda I$ is added to $Q$ to make $Q$ positive definite. Our choice of $\lambda$ is $\left\lceil|\min\{\lambda_{\min},0\}|+r\right\rceil$, where $\lambda_{\min}$ is the minimum eigenvalue of $Q$ and $r$ is a random number between 0 and 1. A feasible sparse integer vector, i.e., $x^{\ast}$, where $x^{\ast}_{i}\in\{0,1,\ldots,50\}$ if $x^{\ast}_{i}\neq 0$, is then generated, and we set $q=-2Qx^{\ast}$. Note that by the fact that $Q$ is positive definite and by our choice of $q$, $x^{\ast}$ is the unique optimal solution to the generated instance of (1). These instances are shown with the prefix convex in Tables 1 and 2. In the non-convex instances, entries of $Q$ and $q$ are from separate $U_{\alpha}$’s. We denote these instances with pairs $(U_{Q},U_{q})$, where $U_{Q}$ and $U_{q}$ are distributions for the quadratic terms ($Q$) and linear terms ($q$), respectively. Our data sets are (U_{2},U_{200}),\ (U_{5},U_{200}),\ (U_{5},U_{10}),\ (U_{5},U_{100}),\ and $(U_{10},U_{0})$.

As mentioned earlier, we take resilience as the measure of success. More specifically, we hypothesized that the bounded-coefficient encoding is a technique for representing a UIQP problem as a UBQP problem that is more robust against noise, and resilience directly measures the robustness of an instance of UBQP against noise. In zoshk2016 , the resilience $R$ of an instance is defined as

[TABLE]

where $n_{\text{same}}$ is the number of times, among all $n_{\text{trial}}$, that the original ground state does not change with random noise perturbations. In other words, $n_{\text{trial}}$ different noise matrices with entries drawn from the normal distribution $\mathcal{N}(0,\epsilon)$ are generated; each is added to a scaled Ising model and the perturbed Ising model is solved exactly. $n_{\text{same}}$, then, refers to the number of times that the perturbed Ising model returns the same ground state as the unperturbed one. The number of trials, $n_{\text{trial}}$, in our experiment is set to 10. We take $\epsilon_{l}$ and $\epsilon_{c}$ to be 0.01, and we scale Ising models to $J\in[-1,1]$ before adding the noise. After finding the upper bound on the coefficients of the encoding, with Algorithm 2 and $\epsilon_{l}=\epsilon_{c}=0.01$, we obtain $\text{Ising}_{\text{bounded}}$ through the bounded-coefficient encoding. Also, by using the binary encoding, we derive $\text{Ising}_{\text{binary}}$. Then, we measure the resilience of each of these Ising models at $\epsilon\in\{0.001,0.002,\ldots,0.01\}$. The resilience of $\text{Ising}_{\text{bounded}}$ and $\text{Ising}_{\text{binary}}$ are summarized in Tables 1 and 2, respectively. The bounded-coefficient encoding has significantly outperformed the binary encoding; $\text{Ising}_{\text{bounded}}$ is five times more resilient to noise than $\text{Ising}_{\text{binary}}$, on average. Also, in most of our test cases, the convex instances and $(U_{10},U_{0})$, the resilience is almost zero with the binary encoding, except for a few exceptions at $\epsilon=0.001$, whereas for bounded-coefficient encoding, the resilience stays above zero and is considerably large for $\epsilon\leq 0.005$. In order to visualize the difference between the two encodings, we plot the average resilience (over all instances) with respect to $\epsilon$, the standard deviation of the added noise, in Figure 1. It is obvious from this plot that the resilience with binary encoding stays marginally above 0, but for the bounded-coefficient encoding it decreases at a slower rate and remains considerably above 0.

Note that in zoshk2016 , it is suggested to set $\epsilon=\frac{2}{J_{\max}}$, where $J_{\max}$ is the largest $J_{ij}$. Here, our approach slightly differs. We assume that we know the precision of the machine, and we derive the upper bound on the coefficients of our encoding to accommodate this restriction. The precision of the current D-Wave machine is in order of $10^{-2}$, which motivates our choice of $\epsilon=10^{-2}$.

5 Conclusion and Discussion

In this paper, we presented an encoding to represent an unconstrained integer quadratic programming problem as an Ising model. To deal with the low precision of quantum annealers, we suggested bounding the value of the coefficients in the encoding. This restricts the range of the local fields biases and couplings’ strengths in the derived Ising model, thereby creating Ising models that are more robust against noise after scaling. Resilience is used as a metric for the robustness of the Ising models. We compared bounded-coefficient encoding with the binary encoding. We infer from our results that bounding the coefficients of the encoding, as in the bounded-coefficient encoding, significantly improves the resilience of the model. The drawbacks are that the size of the derived Ising model is larger and it introduces redundancy.

In our proposed technique, we forced the ratios of the minimum absolute value of the local fields biases and the couplings’ strengths to their respective maximums exceed a certain threshold tolerance that is correlated or equal to the quantum annealer’s precision. However, all of our calculations are embedding-free, i.e., they ignore the Chimera graph structure. In theory, having an embedding with equal chain length for all variables and equal chain connectivity for any present quadratic terms will respect our calculation; however, finding such an embedding is nontrivial. Notice that, after the encoding, all spin variables corresponding to an integer variable are connected; additionally, all spin variables corresponding to $x_{i}$ and $x_{j}$ are connected if the term $x_{i}x_{j}$ appears in the integer formulation; so, the underlying graph of the Ising model is quite dense. Even for our small instances, the resulting Ising models either exceed the size that can be embedded on the current chip, or the chains’ lengths differ in orders of magnitude (e.g., minimum and maximum chain lengths in orders 1 and 10, respectively). Our algorithm can be tested on the future generation of the quantum annealers with improved connectivity.

According to our results in the previous section, using $\epsilon_{l}=\epsilon_{c}=0.01$ to find the upper bound on the coefficients of the encoding, the derived Ising model stays robust against noise, with an average resilience above 0.5, and minimum resilience above $0.1$ for $\epsilon=0.005$, i.e., $\frac{\epsilon_{c}}{2}$. This could suggest using $2\epsilon_{QA}$ as $\epsilon_{l}$ and $\epsilon_{c}$ for finding the upper bounds on the coefficients of the encoding, where $\epsilon_{QA}$ is the precision of the quantum annealer. This, however, requires more experimentation and could be an avenue for future study.

Appendix 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

[TABLE]

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.,

[TABLE]

and

[TABLE]

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

[TABLE]

Unlike the spin variables for which the diagonal of $C^{t}QC$ became constant, the diagonal of $C^{t}QC$ 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

[TABLE]

where

[TABLE]

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:

[TABLE]

The difference between this case and what was presented in Section 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 afterwards. 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

[TABLE]

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

[TABLE]

Combined with the condition

[TABLE]

we need to initialize $\mu^{x_{i}}$ as

[TABLE]

where

[TABLE]

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}}\geq 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

[TABLE]

and the ratio condition on quadratic terms enforces

[TABLE]

or, equivalently,

[TABLE]

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 (46), i.e., Algorithm 3.

After $\mu^{x_{i}}$’s are calculated to satisfy the conditions (45) and (46), 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:

[TABLE]

and

[TABLE]

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:

[TABLE]

[TABLE]

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}\leq\bar{v}_{2}^{m}\leq\cdots\leq\bar{v}_{n}^{m}$; and sort $v^{M}$ in decreasing order; i.e., $\bar{v}^{M}$ such that $\bar{v}_{1}^{M}\geq\bar{v}_{2}^{M}\geq\cdots\geq\bar{v}_{n}^{M}$. We then whether $\frac{\bar{v}_{1}^{m}}{\bar{v}_{1}^{M}}\geq\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 $[\frac{\bar{v}^{m}_{2,1}}{\bar{v}^{M}_{1,1}},1]$ or $[\frac{\bar{v}^{m}_{1,1}}{\bar{v}^{M}_{2,1}},1]$, respectively, after the rescaling. We choose the option that gives us better lower bound, i.e., if

[TABLE]

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.

Acknowledgements

We are thankful to Helmut G. Katzgraber and Gili Rosenberg for insightful discussions and helpful feedback; and Marko Bucyk for editing the manuscript.