# Quantum annealing with inequality constraints: the set cover problem

Hristo N. Djidjev  Institute of Information and Communication Technologies,
Bulgarian Academy of Sciences, Sofia, Bulgaria  Los Alamos National
Laboratory, Los Alamos, NM 87545, USA

###### Abstract

This paper presents two novel approaches for solving the set cover problem
(SCP) with multiple inequality constraints on quantum annealers. The first
method uses the augmented Lagrangian approach to represent the constraints,
while the second method employs a higher-order binary optimization (HUBO)
formulation. Our experimental analysis demonstrate that both approaches
outperform the standard approach with slack variables for solving problems
with inequality constraints on D-Wave quantum annealers. The results show that
the augmented Lagrangian method can be successfully used to implement a large
number of inequality constraints, making it applicable to a wide range of
constrained problems beyond the SCP. The HUBO formulation performs slightly
better than the augmented Lagrangian method in solving the SCP, but it is less
scalable in terms of embeddability in the quantum chip. These findings could
impact the use of quantum annealers for solving constrained optimization
problems.

Keywords: Set cover problem, Quantum annealing; Augmented Lagrangian method;
Quadratic penalty method, D-Wave; Ising problem; QUBO; HUBO

##  1 Introduction

###  1.1 Quantum annealing

Quantum annealers, such as those produced by D-Wave Systems Inc., leverage
quantum mechanical phenomena, including entanglement and tunneling, to tackle
NP-hard optimization problems that are inherently challenging for classical
computers. While D-Wave’s machines are currently the most powerful quantum
devices available, boasting over 5000 qubits, the state of Noisy Intermediate-
Scale Quantum (NISQ) technology is still limited in its abilities and cannot
yet outperform classical computers in solving practical problems due to high
levels of noise and decoherence. However, each new generation of D-Wave
quantum annealers is improving hardware to reduce noise and increase coherence
time, enabling the machines to accurately solve progressively more complex
problems.

To solve a problem on a quantum annealer, it has to be formulated as a problem
of the type

$\mbox{minimize}\;\mathit{Is}(\bm{x})=\sum_{i<j}J_{ij}x_{i}x_{j}+\sum_{i}h_{i}x_{i},$
(1)

where $\bm{x}=\\{x_{1},\dots,x_{n}\\}$, $J_{ij}$ and $h_{i}$ are real numbers,
and variables $x_{i}$ are either in $\\{-1,1\\}$, in which case the problem is
called an Ising problem, or in $\\{0,1\\}$, when the formulation is called a
quadratic unconstrained binary optimization (QUBO) problem. The two
representations are equivalent and can be easily converted into each other by
using a linear variable transformation.

Problem 1, which is a quadratic function of the variables $x_{i}$, is known to
be NP-hard [1] and many important optimization problems can be easily
formulated as Ising or QUBO problems [21]. Such formulations are usually
constructed in two steps. In the first step, the problem of interest is stated
as a 0-1 quadratic programming problem, i.e., a problem to minimize a
quadratic form of $n$ binary variables subject to linear equality or
inequality constraints. In the second step, that constrained problem is
converted into an unconstrained one, which is necessary since problem 1
formulation doesn’t allow constraints. Next we discuss methods to convert
constrained problems into unconstrained ones.

###  1.2 Handling constraints

The penalty method is the most commonly employed technique for dealing with
constrained problems, in which the constraints are included in the objective
function as penalty terms. For instance, a constraint
${\bm{a}}^{\intercal}\bm{x}=b$ can be added to the objective as a penalty term
$P(\bm{x})=\mu(\bm{a}^{\intercal}\bm{x}-b)^{2}$, where $\mu>0$ is a penalty
constant. If $\bm{x}$ satisfies the constraint, then $P(\bm{x})=0$ and the
penalty term doesn’t change the value of the objective. But if $P(\bm{x})\neq
0$, then $P(\bm{x})>0$ and if the constant $\mu$ is chosen large enough, then
the penalty term will prevent $\bm{x}$ to be a minimum of the combined
objective.

If we have an inequality constraint, we can first convert it into an equality
one and then proceed as described above. For instance, a constraint
${\bm{a}}^{\intercal}\bm{x}\leq b$ can be represented as
${\bm{a}}^{\intercal}\bm{x}+d=b$, where $d\geq 0$ is a new slack variable. But
since problem 1 accept only binary variables and $d$ can potentially be as
large as $b$, integer variable $d$ has to be encoded using $\lfloor\log
b\rfloor+1$ binary variables. This may lead to a large number of new
variables, especially if the problem has multiple inequality constraints. But
with the size of the problem increasing, the limited availability of qubits
and the reduced precision of current generation quantum devices can pose
significant challenges [22].

Another issue with the penalty method is that a large penalty constant $\mu$
can significantly impact the accuracy of quantum annealing. This is because
all $J_{ij}$ and $h_{i}$ coefficients from 1 are normalized before being
submitted to the annealer to satisfy hardware-imposed restrictions of the
D-Wave device. Consequently, a large value for $\mu$ leads to some
coefficients of the resulting problem being very small in absolute value. This
poses a challenge due to the analog nature of the quantum device and the
finite precision of the digital-to-analog converter, resulting in further
degradation of the accuracy of the quantum annealing.

###  1.3 The set cover problem

The problem we are considering in this paper, whose formulation involves
dealing with multiple inequality constraints, is the set cover problem (SCP).
SCP is a classical optimization problem with multiple applications including
in scheduling, resource allocation, logistics, and bioinformatics [30]. The
weighted version of the problem is, given a set $U$ of $n$ elements
$\\{1,\dots,n\\}$ and a set $M$ of $m\geq 2$ sets $S_{i}\subset U$ with
positive weights $\mathit{wt}_{i}$ such that $\bigcup_{S\in M}S=U$, to find
$M^{*}\subseteq M$ for which

$\bigcup_{S_{i}\in M^{*}}S_{i}=U$

that minimizes

$\sum_{S_{i}\in M^{*}}\mathit{wt}_{i}.$

In its unweighted version, all weights are one, and both versions are NP-hard
[18]. In this paper we consider the weighted version of the SCP.

###  1.4 Our objectives

In this work, we propose two new approaches for more accurate constraint
handling in solving larger SCP problems on current quantum annealers. The
first approach, described in the next section, uses the augmented Lagrangian
method to represent inequality constraints, as an alternative to the penalty
method. The second approach formulates the SCP as a HUBO problem, which is
similar to the QUBO formulation presented in 1 but permits higher degree
monomials. These methods are intended to address the limitations imposed by
the restricted number of qubits and the reduced precision of quantum annealers
when handling a large number of variables.

This work makes three main contributions. Firstly, we demonstrate that the
augmented Lagrangian method can effectively handle a large number of
inequality constraints when solving constrained problems with quantum
annealing. While our application focuses on the SCP, this approach has broad
applicability to other problem types. Secondly, we show that a HUBO
representation can be leveraged to solve the SCP on a quantum annealer. This
method can be extended to other problems with inequality constraints, although
the problem must have certain structural characteristics. Lastly, our proposed
methods produce QUBO or HUBO problems with only $m$ variables, and our
experimental results indicate that the current D-Wave Advantage machines can
solve problems with up to 400 sets (variables).

The paper is structured as follows. In Section 2, we give a brief literature
review of related results. In Section 3, we provide a brief introduction to
the augmented Lagrangian method and the HUBO formulation, and detail our
proposed algorithms for handling constraints in quantum annealing. In Section
4, we present the experimental results, including a comparison of the
accuracies achieved by our methods and those achieved by the penalty method.
Finally, we conclude the paper with a summary of our findings and suggestions
for future research directions.

##  2 Previous work

Given the set cover problem’s significance, many classical heuristic
algorithms have been proposed for its solution. For the unweighted version of
the problem, Johnson [16] and Lovász [20] proposed algorithms that can find
set covers with cost at most $1+\ln d$ times the optimal one, where $d$ is the
maximum cardinality of sets in $U$, while Chvátal [5] generalized their result
for the weighted version. Feige [8] showed that the set cover cannot be
approximated within a factor of $\ln n$. A review of several algorithms for
the SCP and comparison of their practical performances is given in [4].

There is much less work reported in the literature on quantum algorithms for
the SCP. Lucas [21] describes a QUBO formulation for the SCP, among other NP-
hard problems, but does not implement it on a quantum annealer. The number of
variables for his formulation is $m+n(\log m+1)$. For the related problem of
minimum vertex cover, Pelofske et al. [23] design a quantum annealing
algorithm that can deal with problems that are too large to fit onto the
quantum hardware by decomposing them into smaller subproblems. Zhang et al.
[32] propose quantum approximate optimization algorithm (QAOA), which uses a
model for quantum computing different from quantum annealing, for the minimum
vertex cover problem and apply it to graphs of ten vertices. However, the
minimum vertex cover problem is simpler than set cover in the sense that its
standard formulation involves a single equality constraint and no inequality
ones. Cao et al. [3] consider the set cover with pairs problem, introduced in
[10], and propose a QUBO formulation for it that uses $O(nm^{2})$ binary
variables. They run a quantum annealing simulator on instances requiring no
more than 19 variables, and also solve the QUBO problems using simulated
annealing on the same set of instances.

The augmented Lagrangian method (ALM) was introduced in 1969 by Hestenes and
Powell [11, 25] and has been widely studied in the field of optimization.
Different variants of the method have been applied to a wide range of
problems, including quadratic programming [27], nonlinear programming [6], and
convex optimization [15]. In quantum computing, ALM has been used in [26] for
solving the quantum-chemical ground-state energy problem on a gate-based
quantum computers. Yonaga et al. [31] use the alternating direction method of
multipliers, a variant of ALM, to solve the quadratic knapsack problem.
Djidjev [7] applies ALM for representing logical qubits in quantum annealers.

We also use higher order binary optimization (HUBO) problem formulations,
which have been the subject of intense research. Several works have focused on
quadratization, or reduction of the HUBO to a quadratic form. Kolmogorov and
Zabih [19] and Freedman and Drineas [9] show how to quadratize any monomial
with a negative coefficient, regardless of its degree, by introducing a single
auxiliary variable. Ishikawa [13] propose a method that results in a more
efficient quadratization for positive monomials, utilizing approximately half
the number of variables compared to previous techniques. Methods that don’t
introduce auxiliary variables but instead enumerate assignments to a small
subset of the variables in order to reduce the degree include using deductions
[28] and excludable local configurations [14]. Boros and Gruber [2] review the
previous work on quadratization and propose new techniques. HUBO formulations
have been used for quantum annealing by Pelofske et al. for Boolean tensor
networks [24], by Mato for molecule unfolding, and Jun for prime factorization
[17].

##  3 Methods

###  3.1 Using slack variables

The standard approach to deal with the inequalities of the SCP is using slack
variables in order to convert them to equalities and then use the quadratic
penalty method to incorporate the resulting equalities into the objective
function [21]. Specifically, we define a binary variable $x_{i}\in\\{0,1\\}$
for each set $S_{i}\in M$ that indicates whether the set is included in the
final solution or not. Then the objective function to minimize is

$Q_{A}=\sum_{i=1}^{m}x_{i}\mathit{wt}_{i}.$

If we define $\sigma_{i}=\\{j~{}|~{}S_{j}\ni\,i\\}$, the constraint that at
least one of the selected sets covers element $i$ is
$\sum_{j\in\sigma_{i}}x_{j}\geq 1$. We convert this into the equality

$\sum_{j\in\sigma_{i}}x_{j}=d_{i}+1,$

where $d_{i}$ is a new integer variable in $[0,m-1]$. We encode $d_{i}$ using
$k=\lfloor\log(m-1)+1\rfloor$ binary variables $x_{i,\alpha}$. Then the QUBO
encoding all the constraints is

$Q_{B}=\sum_{i=1}^{m}\bigg{(}\big{(}\sum_{j\in\sigma_{i}}x_{j}-\sum_{\alpha=0}^{k}2^{\alpha}x_{i,\alpha}-1\big{)}^{2}\bigg{)}.$
(2)

Finally, we combine $Q_{A}$ and $Q_{B}$ into a single QUBO

$Q=Q_{A}+\mu\,Q_{B},$

where $\mu$ is a constant satisfying $\mu>\max\\{\mathit{wt}_{i}\\}$ [21].
Although the penalty constant $\mu$ itself is usually not large, the issue is
with the coefficients $2^{\alpha}$, which may become as large as $m$. In the
next two sections, we describe the proposed new approaches.

###  3.2 Augmented Lagrangian version

####  3.2.1 The general method

The augmented Lagrangian method (ALM) for solving constrained problems
combines the penalty method, used in the previous subsection, with the method
of the Lagrangian multipliers. Specifically, in the case of inequalities, all
inequality constraints of type
$c_{i}(\bm{x})={\bm{a}_{i}}^{\intercal}\bm{x}-b_{i}\leq 0$, $i=1,\dots,n$, can
be included into the objective as an additive term

$\bm{\lambda}^{\intercal}\bm{c}(\bm{x})+\dfrac{\mu}{2}||\bm{c}(\bm{x})||^{2}=\sum_{i=1}^{n}\big{(}\lambda_{i}c_{i}(\bm{x})+\dfrac{\mu}{2}||c_{i}(\bm{x})||^{2}\big{)},$

where the coefficients $\lambda_{i}$ are called Lagrangian multipliers,
$\bm{\lambda}=\\{\lambda_{1},\dots,\lambda_{n}\\}$, and
$\bm{c}=\\{c_{1},\dots,c_{n}\\}$. Coefficients $\bm{\lambda}$ and $\mu$ are
estimated using an iterative procedure as described in Algorithm 1. Note that
a version with both equality and inequality constraints is possible.

Algorithm 1 Augmented Lagrangian method for inequalities

1:initial $\mu>0$, initial $\bm{\lambda}$, increase factor $\rho>1$

2:final $\mu$, $\bm{\lambda}$ for use in the augmented Lagrangian function

3:repeat

4:
$\bm{x}\leftarrow\arg\min_{\scaleto{\bm{x}}{4.6pt}\in\\{0,1\\}^{n}}\,\big{(}f(\bm{x})+\bm{\lambda}^{T}\,\bm{c}(\bm{x})+\dfrac{\mu}{2}\,||\bm{c}(\bm{x})||^{2}\,\big{)}$
$\triangleright$ Can solve on a quantum annealer

5: for $i=1,\dots,n$ do

6: if $c_{i}(\bm{x})>0$ then

7: $\lambda_{i}\leftarrow\lambda_{i}+\mu c_{i}(\bm{x})$

8: end if

9: end for

10: $\mu\leftarrow\rho\mu$

11:until $\bm{c}(\bm{x})\leq\bm{0}$ or iteration limit reached
$\triangleright$ Other stopping criteria possible

In the next subsection we apply the method to the SCP.

####  3.2.2 Applying the AL method to the SCP

First, we formulate the SCP problem as a 0$-$1 linear program with constraints
and then we apply the AL method to get rid of the inequalities. As in Section
3.1, we define a binary variable $x_{i}$ for each $i\in[1,m]$ such that
$x_{i}=1$, if subset $S_{i}$ is selected for the cover, or $x_{i}=0$,
otherwise. Then the SCP can be formulated as

$\displaystyle\underset{x_{i}}{\mbox{minimize}}\quad\sum_{i=1}^{m}x_{i}\mathit{wt}_{i}$
(3) $\displaystyle\mbox{subject to}\quad\sum_{j\in\sigma_{i}}x_{j}\geq
1,\;x_{j}\in\\{0,1\\},\;i=1,...,n.$ (4)

The corresponding augmented Lagrangian function, which is the new objective of
the minimization problem, is

$\mathrm{AL}(\bm{x})=\sum_{i=1}^{m}x_{i}\mathit{wt}_{i}+\sum_{i=1}^{n}\lambda_{i}\big{(}1-\sum_{j\in\sigma_{i}}x_{j}\big{)}+\frac{\mu}{2}\sum_{i=1}^{n}\big{(}1-\sum_{j\in\sigma_{i}}x_{j}\big{)}^{2}$
$=\sum_{i=1}^{m}x_{i}\mathit{wt}_{i}+\sum_{i=1}^{n}\Big{(}\lambda_{i}\big{(}1-\sum_{j\in\sigma_{i}}x_{j}\big{)}+\frac{\mu}{2}\big{(}1-\sum_{j\in\sigma_{i}}x_{j}\big{)}^{2}\Big{)}.$

Since

$\big{(}1-\sum_{j\in\sigma_{i}}x_{j}\big{)}^{2}=1+\sum_{j\in\sigma_{i}}x_{j}^{2}+2\\!\\!\\!\sum_{j<k\in\sigma_{i}}\\!\\!\\!x_{j}x_{k}-2\sum_{j\in\sigma_{i}}x_{j},$

and using that $x_{j}^{2}=x_{j}$ for $x_{j}\in\\{0,1\\}$, we get

$\mathrm{AL}(\bm{x})=\sum_{i=1}^{m}x_{i}\mathit{wt}_{i}+\sum_{i=1}^{n}\Big{(}-\lambda_{i}\sum_{j\in\sigma_{i}}x_{j}+\mu\\!\\!\\!\sum_{j<k\in\sigma_{i}}\\!\\!\\!x_{j}x_{k}-\frac{\mu}{2}\sum_{j\in\sigma_{i}}x_{j}\Big{)}+C,$
$=\sum_{i=1}^{m}x_{i}\mathit{wt}_{i}+\sum_{i=1}^{n}\Big{(}(-\lambda_{i}-\frac{\mu}{2})\sum_{j\in\sigma_{i}}x_{j}+\mu\\!\\!\\!\sum_{j<k\in\sigma_{i}}\\!\\!\\!x_{j}x_{k}\Big{)}+C,$
(5)

where $C$ is a constant, which can be ignored when solving the optimization
problem. From 5, we can get the coefficients $J_{ij}$ and $h_{i}$ of the QUBO
representation 1 and solve that QUBO on a quantum annealer, updating
parameters $\bm{\lambda}$ and $\mu$ at each iteration as specified in
Algorithm 1.

###  3.3 HUBO version

Higher-order binary optimization (HUBO) is a generalization of quadratic
unconstrained binary optimization (QUBO) to higher-order polynomials. Each
HUBO problem can be converted into a QUBO problem by defining auxiliary
variables that encode products of other variables in a way that leads to
decreasing the polynomial degree. For instance, if monomial $x_{1}x_{2}x_{3}$
is part of a HUBO, one can define a new variable $u=x_{1}x_{2}$ and replace
$x_{1}x_{2}x_{3}$ by $uv_{3}$. The constraint $u=x_{1}x_{2}$ can be enforced
using a well-known penalty quadratic function $x_{1}x_{2}-2(x_{1}+x_{2})u+3u$
[12] resulting into the QUBO

$ux_{3}+\mu(x_{1}x_{2}-2(x_{1}+x_{2})u+3u),$

which can replace $x_{1}x_{2}x_{3}$ for sufficiently large penalty $\mu$.
Applying this repeatedly can convert HUBO of any degree into a QUBO.

To formulate a HUBO version of the SCP, we start with the 0-1 linear program
3–4. The constraint $\sum_{j\in\sigma_{i}}x_{j}\geq 1$ means that, for at
least one $j\in\sigma_{i}$, $x_{j}=1$, which means that, for at least for one
$j\in\sigma_{i}$, $1-x_{j}=0$. Hence, the constraints 4 are equivalent to

$\prod_{j\in\sigma_{i}}(1-x_{j})=0,\;i=1,...,n.$

and

$\sum_{i=1}^{n}\prod_{j\in\sigma_{i}}(1-x_{j})=0.$ (6)

Define new binary variables $y_{i}=1-x_{i}$, $i=1,\dots,m$. Replacing $x_{j}$
in 6 and 3 with $y_{j}$ and combining them into a single function, we get the
HUBO formulation of the SCP

$\mbox{minimize}\bigg{(}m-\sum_{i=1}^{m}y_{i}\mathit{wt}_{i}+\mu\sum_{i=1}^{n}\prod_{j\in\sigma_{i}}y_{j}\bigg{)},$
(7)

where $\mu$ is a penalty coefficient. Clearly, the constant $m$ can be ignored
for the optimization.

It is easy to see that it is enough to choose $\mu>\max\\{\mathit{wt}_{i}\\}$.
Assume that $\mu>\max\\{\mathit{wt}_{i}\\}$ and, in a solution $\\{y_{i}\\}$
of 7, there is an element $i$ that is not covered, i.e.,
$\prod_{j\in\sigma_{i}}y_{j}=1$. Then there will exist a set $S_{i}$ that
contains $j$ and is not in the cover, i.e., $y_{i}=1$ (since by the SCP
definition $\bigcup_{i}S_{i}=U$). Adding $S_{i}$ to the cover, i.e., changing
$y_{i}$ to 0, will change the value of the objective function by
$\mathit{wt}_{i}-\mu<0$, which contradicts the assumption that $\\{y_{i}\\}$
is an optimal solution of 7.

##  4 Results

###  4.1 Implementation of the algorithms

For our experiments, we use the D-Wave Advantage_system4.1 quantum annealer,
which we call DWA hereafter, available through the Leap quantum cloud service.
The annealing parameters we use to control the annealing process are
$\texttt{num\\_reads}=1000$ for the number of samples returned per call to the
annealer, $\texttt{annealing\\_time}=100$, which sets the number of
microseconds for the annealing time, and
$\texttt{chain\\_break\\_method}=\texttt{MinimizeEnergy}$. We also use the
flux_biases parameter, which is used to help control some hardware biases.
Unless explicitly stated otherwise, all other parameters are set to their
default values.

We test the proposed new algorithms and compare them against the standard
approach on random instances of the SCP. To generate the test problems, our
generator for SCP instances takes as an input the number of sets $m$, the
number of elements $n$, and the coverage $c$, defined as the average number of
sets covering an element of $U=\\{1,\dots,n\\}$. The generator initially
creates $m$ empty sets and then randomly places elements in sets $S_{i}$ until
the following conditions are satisfied at completion: (i) each element of $U$
is contained in at least two sets, (ii) each set contains at least one
element, and (iii) $\sum_{i}|S_{i}|\geq mc$. Finally, a random weight
$\mathit{wt}_{i}\in[0,1]$ for $i=1,\dots,m$ is assigned to each set $S_{i}$.
We use different size parameters for $m$ in $\\{50,75,100,\dots,400\\}$ and,
for each $m$, we compute a set ${N}_{m}$ of three values for $n$ defined as
${N}_{m}=\\{\lceil{0.5m}\rceil,\,\lceil{0.75m}\rceil,\,m\\}$. For the value of
coverage we use $c=3$. For each type of experiment and combination of $m$ and
$n$, we generate three random instances and average the values of the measured
characteristic over the three instances. If the annealer returns an infeasible
solution, i.e., one that doesn’t correspond to a valid cover, we define as
cost of that solution the sum of all weights, which corresponds to the trivial
cover consisting of all sets.

Name | Inequalities method | Optimization method | Classical/quantum
---|---|---|---
SV_QA | slack variable | quantum annealing | quantum
AL_SA | augmented Lagrangian | simulated annealing | classical
AL_QA | augmented Lagrangian | quantum annealing | quantum
HUBO_SA | HUBO | simulated annealing | classical
HUBO_QA | HUBO | quantum annealing | quantum
Table 1: Implemented algorithms used for the experimental analysis.

The algorithms used in the experiments are the following (see Table 1).

SV_QA: based on the standard slack-variable formulation from Section 3.1 and
implementation on a quantum annealer.

AL_SA: uses the augmented Lagrangian representation as a model and the
simulated annealing method to do the optimization. Simulated annealing [29] is
a classical optimization method that uses a probabilistic approach to explore
the solution space in search for a global minimum. To escape local minima, it
uses gradual cooling from a high temperature to a low one, which changes the
probability of accepting moves increasing the objective function from higher
in the beginning to low at the end. We have opted for simulated annealing in
our experiments because it serves as a classical analogue of quantum annealing
being a general-purpose heuristic for global optimization. To carry out our
experiments, we utilized the implementation of simulated annealing that is
included in the Ocean D-Wave software. The ALM parameters from Algorithm 1 are
$0.5$ for the initial $\mu$, $0$ for the initial $\lambda$, and $1.1$ for the
increase factor $\rho$.

AL_QA: the quantum version of AL_SA,combining ALM with quantum annealing. We
use the same ALM parameters as AL_SA.

HUBO_SA: uses the HUBO formulation from Section 3.3. For solving the HUBO we
use the simulated annealing solver of Ocean. Note that, since simulated
annealing is not restricted to quadratic models, we don’t need a conversion to
a QUBO in this algorithm.

HUBO_QA: the quantum version of HUBO_SA as it uses the HUBO formulation of the
SCP. But in order to solve on DWA, we also need a HUBO-to-QUBO converter, for
which we use the one supplied by the Ocean software.

In the next three subsections, we will analyze the implementations of the
augmented Lagrangian and HUBO algorithms and compare all five algorithms with
respect to the qualities of the proposed solutions to SCP.

###  4.2 Augmented Lagrangian method iterations

Our focus here is on analyzing the iterations of ALM, specifically, examining
the decrease rate in the number of non-satisfied constraints as the iteration
number increases. Figure 1 shows the percentage of uncovered elements (non-
satisfied constraints), for AL_SA and AL_QA and for number of variables
$m\in\\{50,100,150\\}$. The number of ALM iterations is set to 10. We observe
that, in all cases, for smaller values for $n$ we see a better performance.
Although this may seem self-evident, we will see in the next subsection that
decreasing the value of $n$ makes the problem harder in other aspects.
Regarding the dependence on the value of $m$, we see for AL_QA that with
increasing the value of $m$ the percentage of uncovered elements also goes up.
We cannot see such a clear trend for AL_SA. Comparing AL_SA and AL_QA
implementations, we see that, with respect to this criterion, AL_QA has a
better performance.

Figure 1: Number of elements not covered, per iteration step, for
$m\in\\{50,100,150\\}$. Shown at the top is AL_SA and at the bottom it is
AL_QA. A “$\times$” symbol indicates the average iteration number at which the
best solution for the corresponding problem size combination has been
obtained.

Finally, we look at the value of the penalty factor $\mu$ from Algorithm 1.
The magnitude of $\mu$ is important, especially in the case of quantum
annealing, since it is used to scale up some coefficients, and large values of
$\mu$ can negatively affect the accuracy of the annealing, as discussed
earlier. While Figure 1 doesn’t directly show these values, they can be easily
calculated given the iteration number $i$, e.g., in our implementation,
$\mu(i)=0.5(1.1)^{i}$. Figure 1 shows the average iteration number where the
best solution was found. We observe that, in the case $m=50$, the best-
iteration numbers and the corresponding values of $\mu$ are lower compared to
$m=100$ and $m=150$, especially for algorithm AL_QA. However, we don’t see
significant difference when we compare $m=100$ and $m=150$. One possible
explanation is that the number of iterations, ten, may be not enough for some
instances with large values of $m$ and the best iteration number for $m=150$
may be greater than 10. But increasing the number of iterations also increases
the cost of the algorithm.

###  4.3 Number of variables and embeddability

The computational complexity of a classical optimization algorithm goes up
with the number of variables of the instance. For quantum algorithms, larger
number of variables usually means lower quality of the solution. But a large
number of variables may also means that the QUBO does not fit on the quantum
device, so the problem may not be solvable at all. Two key factors influence
the sizes of problems that can be solved on a quantum computing device: the
number of qubits and the connections between them. The DWA device has more
than 5000 qubits, with each qubit connected to no more than 15 other qubits.
The sparsity of connections means that, for most problems, several connected
qubits have to be combined into a single logical qubit to represents one
binary variable. Hence, the sizes of the problems that can be solved on DWA,
in terms of the number of binary variables, may be much smaller than the
number of available qubits. In this section, we analyze what sizes of problems
are solvable for each of the three quantum algorithms: SV_QA, AL_QA, and
HUBO_QA.

First, let us analyze the number of the variables. SV_QA uses a QUBO
formulation with $m+n(\log m+1)$ [21]. AL_QA also uses a QUBO formulation and
its number of binary variables is $m$. Finally, HUBO_QA uses a HUBO
formulation with $m$ variables. But in order to convert it to a quadratic
form, one should define a number of auxiliary variables to reduce the degree
of the monomials, so the final number of variables is higher. The number of
auxiliary variables depends on the implementation, and we don’t have a formula
for the one used by the Ocean software. But we can experimentally analyze them
by counting the number of final variables for each test instance.

Figure 2, bottom, shows the number of binary variables for HUBO_QA for
different SCP sizes. For each value of $m$, we observe that the number of
binary variables increases with decreasing the value of $n$. While this might
look counter-intuitive, it is based on the fact that the number of auxiliary
variables depends on the degrees of the monomials in the HUBO, as each
monomial of degree larger than two needs additional variables to quadratize
it. Furthermore, the degree of the $i$-th monomial is the size of the set
$\sigma_{i}$ (see 7), which is the number of all sets containing element $i$.
Smaller values of $n$ means more sets covering each element, on average and,
hence, higher degree monomials. We can also see that the number of variables
for HUBO_QA grows when $m$ is varied between 50 an 400 from 95.5 to 794.2, for
$n=m$, and from 133.6 to 1170.4, for $n=m/2$.

Figure 2: Comparison of the three quantum algorithms with respect to the
number of variables or couplers. On the top is SV_QA, in the middle it is
AL_QA, and at the bottom it is HUBO_QA. The colors show what fraction the
problem instances are embeddable into the DWA chip.

Figure 2, top and middle, displays the number of couplers (quadratic
coefficients) for SV_QA and AL_QA, respectively. We don’t plot the number of
QUBO variables for these algorithms as they can be calculated using explicit
formulas, which are $m+n(\log m+1)$ and $m$. The number of couplers is
important since, without all-to-all connectivity, larger number of couplers
usually means that more qubits are needed to represent a single QUBO variable,
which reduces the sizes of the embeddable problems. We observe a similar
pattern, which is that the number of couplers increases when $m$ goes up or
$n$ goes down. Specifically, when $n$ decreases, the average size of
$\sigma_{i}$ increases, while the numbers of couplers for $i$ goes up as
roughly $(|\sigma_{i}|+\log m)^{2}$ for SV_QA and $|\sigma_{i}|$ for AL_QA,
see 2 and 5. Specifically, for $m=275$ and $n=138$, the number of couplers for
SV_QA is 7794 and for AL_QA it is 2326.

Whether a particular instance of the SCP can be embedded in the DWA chip
depends on both the number of variables and the number of couplers, but it
also depends on the specific connection patterns, e.g., the locality of the
connections. Since the connection patterns are hard to quantify, the ultimate
criterion for evaluating and comparing the algorithms we use is the likelihood
for problems with given values of $n$ and $m$ to be embeddable. On Figure 2,
the colors indicate what portion of the ten random problems generated for each
specific combination of $n$ and $m$ could be successfully embedded. We can
observe that, for $m$ upto 150, all algorithms produce embeddable QUBOs 100%
of the time. At $m=175$, SV_QA could still embed all problems, except one
instance for $n=138$, and the other two algorithms can embed all variables
100% of the time. Despite that single infeasible instance, we take $m=175$ to
be the largest value for $m$ where all methods are able to produce solutions,
and in the next subsection we compare the methods with respect to their
accuracy for $m$ upto 175.

For values of $m$ between $200$ and $275$, all instances of AL_QA and HUBO_QA
are embeddable, while none of the SV_QA instances are. Although HUBO_QA has
more embeddable instances for $m\geq 300$, both methods have a 33%
embeddability rate for $m=400$, and for $m=425$, no instances are embeddable.
It is noteworthy that HUBO_QA performs just slightly worse than AL_QA with
respect to embeddability in DWA, despite having a substantially higher number
of QUBO variables (ranging from two to three times more, depending on $n$).

###  4.4 Solution quality comparison

Figure 3 shows the average costs of the covers computed by the five algorithms
discussed in this paper for different values of $m$ and $n$. To make the
comparisons easier, the weights on the sets for each value of $m$ have been
renormalized so that the bar height for HUBO_QA is one. Overall, in terms of
solution quality, HUBO_QA performs the best, while SV_QA performs the worst.
To provide a more detailed analysis, we compare the performance of quantum and
classical algorithms, followed by a comparison of the algorithms based on the
optimization method used..

Our analysis compares the performance of quantum and classical versions of our
proposed methods, specifically AL_QA vs. AL_SA and HUBO_QA vs. HUBO_SA. We
found that quantum annealing (QA) consistently outperforms simulated annealing
(SA), particularly in the most challenging case where $n=\lceil m/2\rceil$ and
when comparing AL_QA with AL_SA. On the other had, HUBO_SA is just slightly
worse than HUBO_QA. A possible explanation is that, while the HUBO_QA
performance is degraded with increasing the size of the problem due to the
QUBO-to-HUBO conversion, our implementation of HUBO_SA directly applies
simulated annealing to the HUBO representation and does not suffer from an
increased number of variables.

Figure 3: Comparison of the methods based on the set cover costs of computed
solutions. Each method and combination of $m\in\\{50,75,\dots,175\\}$ is
represented by a bar, with shading indicating a different value of $n$.

Finally, let us compare the algorithms with respect to the methods used to
implement the optimization constraints: slack variables (SV_QA), augmented
Lagrangian (AL_SA and AL_QA), and HUBO (HUBO_QA and HUBO_SA). The slack
variables approach finds worse solutions than the others. Except in one case,
$m=175$, where AL_SA is a bit worse than, in the other cases the solutions
found by SV_QA have cost about twice greater on average than the costs of the
solutions produced by the other methods. In contrast to the other methods, the
performance of SV_QA is relatively consistent across different values of $n$
when $m$ is fixed, resulting in little variation in the quality of the
solutions obtained. The HUBO method performed the best, with both the QA and
SA implementations finding better solutions than the other methods. ALM, in
particular its QA version, performed slightly worse, but seems to be a viable
alternative.

##  5 Conclusion

This paper focuses on the set cover problem and the challenge of implementing
multiple inequality constraints on a quantum annealer. We compare the standard
approach based on slack variables [21] with two new approaches based on the
augmented Lagrangian method (ALM) and higher-order binary optimization (HUBO),
respectively. Our experimental analysis shows that both new approaches
outperform the standard approach. The HUBO approach finds the highest quality
solutions and is easy to implement. However, unlike ALM, which is more general
and can be applied straightforwardly to other problems with inequality
constraints, the HUBO formulation relies on the specific structure of the set
cover problem and may not be applicable to other problems with inequality
constraints. Also, the HUBO approach is less scalable than ALM in terms of
embeddability in the quantum chip of D-Wave.

We also demonstrate that even with a large number of inequality constraints,
the augmented Lagrangian method may be a viable approach for solving
constrained problems on a quantum annealer. We perform experiments with
problems having up to $400$ constraints and find good quality solutions.
However, the quality of solutions tends to degrade slowly with increasing $n$
and $m$, but it is much more affected by the number of quadratic couplers of
the QUBO, which is determined by the ratio of $n/m$.

Our results could be applicable to solving other optimization problems with
constraints on quantum annealers. Future research directions could include
improving the implementation of the augmented Lagrangian optimization
procedure (Algorithm 1), for instance, by updating the stopping criterion or
the method for updating the values of $\mu$ and $\bm{\lambda}$. Additionally,
the conversion from HUBO to QUBO could be improved to produce fewer auxiliary
variables. The Ocean implementation we used is good, but better
implementations are possible, especially those that take into account the
specific structure of the HUBO.

## Acknowledgments

This work was supported by grant number KP-06-DB-11 of the Bulgarian National
Science Fund and by the Laboratory Directed Research and Development program
of Los Alamos National Laboratory under project 20210114ER. Los Alamos
National Laboratory is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of U.S. Department of Energy
(contract No. 89233218CNA000001).

## References

  * [1] F. Barahona.  On the computational complexity of Ising spin glass models.  Journal of Physics A: Mathematical and General, 15(10):3241, 1982\.
  * [2] Endre Boros and Aritanan Gruber.  On quadratization of pseudo-boolean functions.  arXiv preprint arXiv:1404.6538, 2014.
  * [3] Yudong Cao, Shuxian Jiang, Debbie Perouli, and Sabre Kais.  Solving set cover with pairs problem using quantum annealing.  Scientific reports, 6(1):1–15, 2016.
  * [4] Alberto Caprara, Paolo Toth, and Matteo Fischetti.  Algorithms for the set covering problem.  Annals of Operations Research, 98(1):353–371, 2000.
  * [5] Vasek Chvatal.  A greedy heuristic for the set-covering problem.  Mathematics of operations research, 4(3):233–235, 1979.
  * [6] G Di Pillo and L Grippo.  A new class of augmented lagrangians in nonlinear programming.  SIAM Journal on Control and Optimization, 17(5):618–628, 1979.
  * [7] Hristo N Djidjev.  Logical qubit implementation for quantum annealing: augmented lagrangian approach.  arXiv preprint arXiv:2301.12393, 2023.
  * [8] Uriel Feige.  A threshold of ln n for approximating set cover.  Journal of the ACM (JACM), 45(4):634–652, 1998.
  * [9] Daniel Freedman and Petros Drineas.  Energy minimization via graph cuts: Settling what is possible.  In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), volume 2, pages 939–946. IEEE, 2005.
  * [10] Refael Hassin and Danny Segev.  The set cover with pairs problem.  In International Conference on Foundations of Software Technology and Theoretical Computer Science, pages 164–176. Springer, 2005.
  * [11] Magnus R Hestenes.  Multiplier and gradient methods.  Journal of optimization theory and applications, 4(5):303–320, 1969\.
  * [12] D-Wave Inc.  Problem-solving handbook. reformulating a problem.  https://docs.dwavesys.com/docs/latest/handbook_reformulating.html.  Accessed: 2023-01-16.
  * [13] Hiroshi Ishikawa.  Transformation of general binary MRF minimization to the first-order case.  IEEE transactions on pattern analysis and machine intelligence, 33(6):1234–1249, 2011.
  * [14] Hiroshi Ishikawa.  Higher-order clique reduction without auxiliary variables.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1362–1369, 2014.
  * [15] Alfredo N Iusem.  Augmented lagrangian methods and proximal point methods for convex optimization.  Investigación Operativa, 8(11-49):7, 1999.
  * [16] David S. Johnson.  Approximation algorithms for combinatorial problems.  J. Comput. Syst. Sci., 9(3):256–278, dec 1974.
  * [17] Kyungtaek Jun.  HUBO and QUBO models for prime factorization.  arXiv preprint arXiv:2301.06738, 2023.
  * [18] Richard M. Karp.  Reducibility among combinatorial problems.  In Raymond E. Miller and James W. Thatcher, editors, Complexity of Computer Computations, The IBM Research Symposia Series, pages 85–103. Plenum Press, New York, 1972.
  * [19] Vladimir Kolmogorov and Ramin Zabin.  What energy functions can be minimized via graph cuts?  IEEE transactions on pattern analysis and machine intelligence, 26(2):147–159, 2004.
  * [20] László Lovász.  On the ratio of optimal integral and fractional covers.  Discrete mathematics, 13(4):383–390, 1975.
  * [21] A. Lucas.  Ising formulations of many NP problems.  Front Physics, 2(5), 2014.
  * [22] E. Pelofske, G. Hahn, and H. Djidjev.  Solving large Maximum Clique problems on a quantum annealer.  Proceedings of the First International Workshop on Quantum Technology and Optimization Problems QTOP’19, 2019.
  * [23] Elijah Pelofske, Georg Hahn, and Hristo Djidjev.  Solving large minimum vertex cover problems on a quantum annealer.  In Proceedings of the 16th ACM International Conference on Computing Frontiers, pages 76–84, 2019.
  * [24] Elijah Pelofske, Georg Hahn, Daniel O’Malley, Hristo N Djidjev, and Boian S Alexandrov.  Quantum annealing algorithms for boolean tensor networks.  Scientific Reports, 12(1):8539, 2022.
  * [25] Michael JD Powell.  A method for nonlinear constraints in minimization problems.  In Fletcher, editor, Optimization, pages 283–298. Academic Press, 1969.
  * [26] Robert Schade, Carsten Bauer, Konstantin Tamoev, Lukas Mazur, Christian Plessl, and Thomas D Kühne.  Parallel quantum chemistry on noisy intermediate-scale quantum computers.  Physical Review Research, 4(3):033160, 2022.
  * [27] Klaus Schittkowski.  On the convergence of a sequential quadratic programming method with an augmented Lagrangian line search function.  Mathematische Operationsforschung und Statistik. Series Optimization, 14(2):197–216, 1983.
  * [28] Richard Tanburn, Emile Okada, and Nike Dattani.  Reducing multi-qubit interactions in adiabatic quantum computation without adding auxiliary qubits. part 1: The" deduc-reduc" method and its application to quantum factorization of numbers.  arXiv preprint arXiv:1508.04816, 2015.
  * [29] Peter JM Van Laarhoven, Emile HL Aarts, Peter JM van Laarhoven, and Emile HL Aarts.  Simulated annealing.  Springer, 1987.
  * [30] Rao R Vemuganti.  Applications of set covering, set packing and set partitioning models: A survey.  Handbook of Combinatorial Optimization: Volume1–3, pages 573–746, 1998.
  * [31] Kouki Yonaga, Masamichi J Miyama, and Masayuki Ohzeki.  Solving inequality-constrained binary optimization problems on quantum annealer.  arXiv preprint arXiv:2012.06119, 2020.
  * [32] YJ Zhang, XD Mu, Xiao-Wen Liu, XY Wang, Xingyi Zhang, Kai Li, TY Wu, D Zhao, and Cheng Dong.  Applying the quantum approximate optimization algorithm to the minimum vertex cover problem.  Applied Soft Computing, 118:108554, 2022.