# Quantum Sparse Coding

Yaniv Romano Equal contribution. Departments of Electrical and Computer
Engineering and of Computer Science, Technion IIT, Israel  LightSolver LTD.,
Tel Aviv, Israel  Harel Primack††footnotemark:  LightSolver LTD., Tel Aviv,
Israel  Talya Vaknin††footnotemark:  LightSolver LTD., Tel Aviv, Israel  Idan
Meirzada  LightSolver LTD., Tel Aviv, Israel  Ilan Karpas  LightSolver LTD.,
Tel Aviv, Israel  Dov Furman  LightSolver LTD., Tel Aviv, Israel  Chene
Tradonsky  LightSolver LTD., Tel Aviv, Israel  Ruti Ben Shlomi  LightSolver
LTD., Tel Aviv, Israel

###### Abstract

The ultimate goal of any sparse coding method is to accurately recover from a
few noisy linear measurements, an unknown sparse vector. Unfortunately, this
estimation problem is NP-hard in general, and it is therefore always
approached with an approximation method, such as lasso or orthogonal matching
pursuit, thus trading off accuracy for less computational complexity. In this
paper, we develop a quantum-inspired algorithm for sparse coding, with the
premise that the emergence of quantum computers and Ising machines can
potentially lead to more accurate estimations compared to classical
approximation methods. To this end, we formulate the most general sparse
coding problem as a quadratic unconstrained binary optimization (QUBO) task,
which can be efficiently minimized using quantum technology. To derive at a
QUBO model that is also efficient in terms of the number of spins (space
complexity), we separate our analysis into three different scenarios. These
are defined by the number of bits required to express the underlying sparse
vector: binary, 2-bit, and a general fixed-point representation. We conduct
numerical experiments with simulated data on LightSolver’s quantum-inspired
digital platform to verify the correctness of our QUBO formulation and to
demonstrate its advantage over baseline methods.

Keywords: Quantum Annealing, Sparse Regularization, Compressed Sensing,
Simulated Annealing, Sparse Coding, Sparse Pursuit, Feature Selection.

##  1 Introduction

A ubiquitous problem in machine learning, statistics, and signal processing is
to accurately estimate an unknown sparse vector from a few noisy linear
measurements. This estimation problem, which we refer to as _sparse coding_ ,
is at the heart of the field of compressed sensing, revealing that under
sparsity assumptions it is possible to successfully recover a signal that
sampled significantly below the Nyquist rate [1, 2]. This, in turn, led to a
dramatic increase in magnetic resonance imaging (MRI) scanning session speed
[3]. Another exciting application that also builds on the sparsity assumption
is unsupervised representation learning, i.e., given high-dimensional input
data, such as an image, finding a low-dimensional representation that captures
the intrinsic underlying structure in the input [4, 5, 6]. These
representations are often used in image restoration tasks to effectively
remove noise (denoising) [7, 8], fill-in missing pixels (inpainting) [9, 10,
11], and to achieve high quality digital zoom (super-resolution) [10, 12, 13,
14]. Sparsity also plays a key role in linear regression when given a large
pool of features, to form a predictive rule that estimates an unknown response
using a smaller, _interpretable_ subset of features that manifests the
strongest effects [15, 16, 17, 18].

To formalize the sparse coding problem, which is central for tackling the
aforementioned applications, we consider the following linear model:

${\bm{b}}={\bm{A}}{\bm{x}}+{\bm{v}},$ (1)

where ${\bm{A}}$ is a matrix of size ${M\times N}$, the vector ${\bm{x}}$ is
of length $N$, and ${\bm{v}}$ is a noise vector of length ${M}$. In this
paper, we focus on a challenging setting in which $M\ll N$, where a crucial
assumption we make is that the vector ${\bm{x}}$ is $k$-sparse, i.e., it
contains only $k$ non-zero elements with $k\ll N$ [2, 1, 19]. In the context
of sparse linear regression, the vector ${\bm{b}}$ is called a _response_ ,
the columns of ${\bm{A}}$ are the _features_ , and the rows correspond to the
_observations_. In the context of unsupervised representation learning, the
matrix ${\bm{A}}$ is called a _dictionary_ , whose columns are referred to as
basis elements. The assumption in this literature is that the underlying
signal can be represented as a _sparse_ linear combination of a few columns
taken from the dictionary ${\bm{A}}$, providing in a form of a compression
effect [4, 5, 6].

Suppose we are given a matrix ${\bm{A}}$ and a noisy vector ${\bm{b}}$, and
our goal is to find a sparse solution $\hat{{\bm{x}}}$ to the following
optimization problem:

$\displaystyle\hat{{\bm{x}}}=\arg\min_{{\bm{x}}}\|{\bm{A}}{\bm{x}}-{\bm{b}}\|_{2}^{2}\
\ \text{subject to}\ \|{\bm{x}}\|_{0}=k,$ (2)

where $\|{\bm{x}}\|_{0}$ is the $L_{0}$ pseudo-norm, counting the number of
non-zeros in ${\bm{x}}$. Unfortunately, the above is a non-convex optimization
problem, and finding a sparse solution that minimizes (2) is NP-hard in
general [20, 21]: intuitively, to solve (2), we should sweep over all possible
combinations of the $k$-sparse vectors and choose the one that minimizes the
squared error term. Problem (2) is often expressed in the following Lagrangian
form:

$\displaystyle\hat{{\bm{x}}}=\arg\min_{{\bm{x}}}\|{\bm{A}}{\bm{x}}-{\bm{b}}\|_{2}^{2}+\lambda\|{\bm{x}}\|_{0},$
(3)

where a proper choice of $\lambda$—controlling the strength of the sparsity
penalty (larger $\lambda$ leads to sparser $\hat{{\bm{x}}}$)—may render (2)
and (3) to yield the same solution. The desire to find a solution for the non-
convex sparse coding problem has led to practical approximation algorithms,
which are supported by theoretical guarantees. For example, the well-known
lasso [15] and basis pursuit [17] algorithms suggest replacing the non-convex
$L_{0}$ with the convex $L_{1}$ norm in (3). This choice is attractive since
one can utilize the convexity of the modified objective, and minimize it
efficiently [22, 23, 24]. Another popular algorithm is the orthogonal matching
pursuit (OMP) [25, 26, 27]—a method that greedily picks one non-zero at a time
in lieu of exhaustively sweeping over all the possible solutions. Loosely
speaking, following [19, 28], these methods perform particularly well when (i)
the sparsity level $k$ of the underlying ${\bm{x}}$ is “small enough”; and
(ii) the columns of ${\bm{A}}$ are of low coherence, implying that any subset
of $k$ columns is almost an orthonormal submatrix [29]. In fact, under certain
assumptions on $k$, the correlation structure of ${\bm{A}}$, and the noise
level, lasso and OMP can provably find the correct locations of the non-zero
elements in the unknown ${\bm{x}}$, and, in turn, obtain an estimate
$\hat{{\bm{x}}}$ that is close to ${\bm{x}}$, e.g., in the Euclidean distance
[27, 28, 1, 2, 19].

The NP-hardness nature of the sparse coding problem brings us to the
fascinating and ever-evolving field of quantum computing. The rise of quantum
or quantum-inspired technologies unlock new opportunities to solve more
accurately NP-hard problems than classic computing machines. Notable examples
include the traveling salesman problem [30], protein folding [31], and graph
clustering [32]. Many of these problems can be cast as a quadratic
unconstrained binary optimization (QUBO) problem [33], a formulation that is
of particular interest to us as well. Furthermore, there has been growing
activity in leveraging quantum technologies to tackle machine learning
problems, such as linear regression and classification [34], unsupervised
clustering [34, 35], and even for improving deep generative models [36].

###  1.1 Our contribution

Due to the computational difficulty of minimizing the original problem (3), it
is always being replaced by one of its surrogates, e.g., lasso and OMP.
However, in many cases, the conditions for guaranteeing the validity of the
surrogates are either unknown or not satisfied. This paper “goes back to
basics” and introduces novel quantum machine learning algorithms, directly
minimizing the non-convex problem (3) with the ultimate goal of obtaining more
accurate solutions than standard approximation algorithms. Concretely, we show
how to formulate (3) as a QUBO problem, which, in turn, can be minimized
efficiently using cutting-edge solving platforms, such as quantum computers,
Ising machines [37], or well-known heuristics including Simulated Annealing
and Tabu Search [38, 39]. Our derivation, presented in Section 3, covers (i) a
simple case where each entry in the unknown sparse vector ${\bm{x}}$ is binary
(1-bit representation); (ii) a more flexible setting, using a 2-bit
representation per element in ${\bm{x}}$; and (iii) the most general
formulation where ${\bm{x}}$ is represented using any arbitrary fixed-point
precision. We also analyze the complexity of the proposed formulation with
respect to the number of qubits (or spins) required to solve the QUBO problem.
Naturally, the binary and 2-bit cases are highly efficient in terms of space
complexity—we show how to form the QUBO problem without introducing any
auxiliary (ancilla) spins. By contrast, the most general formulation in which
${\bm{x}}$ can have fractional, non-integer entries comes at the cost of
increasing the space complexity. Yet, our formulation is highly efficient as
it requires only one ancilla spin for each entry in ${\bm{x}}$, as we
carefully argue in Section 3.2. Numerical experiments, given in Section 4,
demonstrate the superiority of our QUBO-based solutions over the classic and
most commonly used approximation algorithms—lasso and OMP. Specifically, we
minimize the proposed QUBO problem via the quantum-inspired computing
framework of LightSolver [37], and show how our method tends to produce more
accurate estimations of the unknown ${\bm{x}}$ in regimes where the sample
size is small or when the number of non-zero elements to be estimated is
relatively large.

###  1.2 Related work

The general QUBO problem can be written as [33]

$\min_{{\bm{q}}\in\mathbb{B}^{D}}{\bm{q}}^{\top}{\bm{W}}{\bm{q}},$ (4)

where $\mathbb{B}=\\{0,1\\}$ a binary set and ${\bm{W}}\in\mathbb{R}^{D\times
D}$ is a real-valued QUBO matrix. The ability to efficiently minimize (4)
using quantum computers paves the way to the design of powerful solvers for
fundamental machine learning problems. Importantly, QUBO optimization becomes
appealing when the runtime of quantum-based algorithms is smaller than that of
classic computers. Motivated by this capability, [40] presented a QUBO
formulation for solving a linear system of equations
${\bm{A}}{\bm{x}}={\bm{b}}$. The authors of [34] broadened this scope even
further and showed how to fit least squares regressors, support vector machine
(SVM) classifiers, and balanced k-means clustering algorithms, by translating
these problems into a QUBO scheme, too. Notably, the above algorithms were
shown to have improved time and space complexities compared to classic SVM and
balanced k-means clustering algorithms, and similar complexity to that of
classic least squares solvers. Our work enriches the quantum machine learning
toolbox by introducing a novel QUBO formulation for NP-hard sparse coding
problems. Stressing this point, since we express the sparsity penalty in
general QUBO terms, we could augment it to the QUBO SVM method of [34] and
tackle the sparse SVM problem [41] instead. We believe such an extension could
be of great interest by itself, and we leave it for future work.

Recently, [42, 43, 44] presented a QUBO formulation for solving a special case
of (3)—also referred to as binary compressive sensing—for which the unknown
${\bm{x}}$ is sparse and _binary_. These quantum-based solutions were shown to
yield improved performance compared to the classic lasso method, an
observation that is also corroborated by our experiments. Our paper builds
upon the above line of work and extends the QUBO formulation to _non-binary_
${\bm{x}}$, represented in fixed-point arithmetic. A biological problem that
is tightly connected to binary sparse signal recovery is the estimation of the
transcription factor in DNA binding. By making a linear modeling assumption,
[45] expressed the gene sequence as a QUBO problem with _binary_ sparse
regularization. The authors compared their quantum-based estimation to
existing baseline methods (including lasso) and reported that the most
significant improvement is obtained in small sample size regimes. This
observation is also supported by our experiments, although conducted in a
different context as well as in the challenging setting where ${\bm{x}}$ is a
non-binary vector. The method reported in [46] is the closet to ours, as it
presents a QUBO formulation for the recovery of sparse signals represented in
fixed-point arithmetic. Our formulation, however, is much more efficient with
respect to the number of spins required to derive the QUBO model. Let $P$ be
the number of bits defining the resolution of the fixed point representation.
To present (3) as a QUBO problem for $P\geq 3$, the approach presented in [46]
introduces $P-2$ auxiliary spins per variable while ours uses only one! We
conclude with [47], which “translated to QUBO” several feature selection
methods. Instead of minimizing (3), the authors of [47] assigned each feature
with an information-theoretic measure of importance—e.g., its correlation with
the target variable—and selected the ones that deemed to be the most
important. This approximation allows handling a non-binary ${\bm{x}}$ by
assigning a single binary spin per feature, indicating whether the feature
should be selected or not, akin to binary compressive sensing problems.
Importantly, the formulations of the feature selection methods in [47] are
substantially different than ours—we aim to tackle the original, ubiquitous
sparse coding problem (3) in its most general form, rather than an alternative
formulation of it.

##  2 A motivating example

Before describing the proposed QUBO formulations, we pause to present a
motivating, small-scale example that demonstrates the advantage of an
exhaustive method for solving the sparse coding problem compared to classical
approximation methods. To this end, we generate a matrix ${\bm{A}}$ with
$L_{2}$-normalized columns, i.e., $\|{\bm{a}}_{i}\|_{2}=1$ for all $1\leq
i\leq N$, where ${\bm{a}}_{i}$ is the $i$th column of ${\bm{A}}$. We design
the matrix ${\bm{A}}$ such that its _mutual coherence_ , defined as [48, 28,
19]

$\mu({\bm{A}})=\max_{i\neq j}|{\bm{a}}_{i}^{\top}{\bm{a}}_{j}|,$ (5)

is minimized.111A detailed procedure for generating matrices with low
coherence is given in Appendix A. We make this design choice since low mutual
coherence is crucial for the success of lasso and OMP in recovering the
unknown ${\bm{x}}$ [19]. Put simply, we wish to show that an algorithm that
exhaustively search for a solution for (3) is highly valuable, even when
“playing on the turf” of lasso and OMP.

(a) Cardinality $k=3$, noise $\sigma=0.1$.

(b) Cardinality $k=3$, rows $M=8$.

(c) Rows $M=8$, noise $\sigma=0.1$.

Figure 1: Comparison between the solutions obtained by Algorithm 1, which
exhaustively minimizes (2), and two classic approximation algorithms—lasso and
OMP. (a) Estimation error as a function the number of rows. (b) Estimation
error as a function the the noise level $\sigma$. (c) Estimation error as a
function of the number non-zeros $k$ in the true ${\bm{x}}$. Each point in the
graphs is an average over 30 independent realizations of ${\bm{A}},{\bm{x}}$,
and ${\bm{v}}$; the length of ${\bm{x}}$ is fixed and equals to $N=16$.

Given such a matrix ${\bm{A}}$, we generate a vector
$\bar{{\bm{b}}}={\bm{A}}{\bm{x}}$ by drawing an ${\bm{x}}$ with $k$ randomly
chosen non-zero entries that are equal to $1$, and then forming its noisy
version ${\bm{b}}=\bar{{\bm{b}}}+{\bm{v}}$, where the noise component
${\bm{v}}$ is sampled from a Normal distribution with mean zero and standard
deviation (STD) $\sigma$. Now, given ${\bm{A}}$ and the noisy ${\bm{b}}$, we
attempt to reconstruct ${\bm{x}}$ via lasso, OMP, and an exhaustive algorithm
for minimizing (2), detailed in Algorithm 1. Due to the high computational
complexity of Algorithm 1, we could run this experiment only for a matrix
${\bm{A}}$ with a small the number of columns, which we set to $N=16$. The
hyper-parameter of each method (the $L_{1}$ regularization strength in lasso,
and the cardinality level in OMP and in Algorithm 1) is tuned by an “oracle”
that has access to the true ${\bm{x}}$, optimizing each of the following two
error metrics. The first metric evaluates the reconstruction error, defined as
$\|{\bm{x}}-\hat{{\bm{x}}}\|_{2}/\|{\bm{x}}\|_{2}$. The second measures the
support recovery error, defined as $\|\mathcal{S}-\hat{\mathcal{S}}\|_{0}$,
where $\mathcal{S}\in\mathbb{B}^{N}$ is a binary vector that contains the
value $\mathcal{S}_{i}=1$ in the $i$th index if ${\bm{x}}_{i}\neq 0$, or the
value $\mathcal{S}_{i}=0$ otherwise. The vector
$\hat{\mathcal{S}}\in\mathbb{B}^{N}$ is defined similarly, but with respect to
the estimated vector $\hat{{\bm{x}}}$. In plain words, the support recovery
error is the number of non-zero elements in the unknown ${\bm{x}}$ and
estimated $\hat{{\bm{x}}}$ that do not coincide in their locations.

Figure 1 compares the performance of each method as a function of the sample
size $M$ (left), noise level $\sigma$ (middle), and cardinality $k$ (right).
Following that figure, we can see that the exhaustive method for minimizing
(2) outperforms lasso and OMP by a large gain, especially (i) in small sample
size regimes, (ii) for small noise levels, and (iii) for higher cardinality
values. This toy experiment, which explicitly shows how powerful is the
original $L_{0}$ regularization compared to its surrogates, elucidates our
strong desire to leverage quantum computing technologies to better approximate
the NP-hard sparse coding problem, and further scale-up the dimension of
${\bm{A}}$ to have a much larger number of columns.

Algorithm 1 An exhaustive solver for (2).

A matrix ${\bm{A}}$ of size $M\times N$, a noisy vector ${\bm{b}}$ of length
$M$, and a sparsity level $k$.

A $k$-sparse vector $\hat{{\bm{x}}}$ minimizing (2).

$S\leftarrow$ a set containing all the possible binary support vectors with
$k$ non-zero elements $\triangleright$ $N\choose k$ vectors

$\hat{{\bm{x}}}\leftarrow\bm{0}$, where $\bm{0}$ is an $N$-dimensional vector
of zeros.

$R\leftarrow\|{\bm{b}}\|_{2}^{2}$. $\triangleright$ Squared error of the zero-
vector solution

for all ${\bm{s}}$ in $S$ do

${\bm{A}}_{{\bm{s}}}\leftarrow$ a matrix of size $M\times k$ containing $k$
columns from ${\bm{A}}$, specified by the support vector ${\bm{s}}$.

$\hat{{\bm{z}}}_{{\bm{s}}}\leftarrow({\bm{A}}_{{\bm{s}}}^{\top}{\bm{A}}_{{\bm{s}}})^{-1}{\bm{A}}_{{\bm{s}}}^{\top}{\bm{b}}$.
$\triangleright$ Solution for the least squares problem:
$\min_{{\bm{z}}}\|{\bm{A}}_{{\bm{s}}}{\bm{z}}-{\bm{b}}\|_{2}^{2}$

${r}\leftarrow\|{\bm{A}}_{{\bm{s}}}\hat{{\bm{z}}}_{{\bm{s}}}-{\bm{b}}\|_{2}^{2}$.

if ${r}<R$ then $\triangleright$ The current solution is better that the best
we found thus far

$R\leftarrow r$.

$\hat{{\bm{x}}}_{{\bm{s}}}\leftarrow\hat{{\bm{z}}}_{s}$, where the off-support
elements that are not specified by ${\bm{s}}$ are set to zero.

end if

end for

##  3 Quantum solvers for sparse coding problems

The problem of solving a system of linear equations with minimum errors can be
reformulated as minimizing a quadratic function of binary variables, i.e., as
a QUBO problem. The derivation presented in this section decomposes the
objective function (3) into the squared error term
$\|{\bm{A}}{\bm{x}}-{\bm{b}}\|_{2}^{2}$ and the sparsity penalty term
$\|{\bm{x}}\|_{0}$, forming the QUBO matrix for each component, denoted by
${\bm{W}}^{L_{2}}$ and ${\bm{W}}^{L_{0}}$, respectively. At a high level, the
QUBO problem that we formalize takes the following form:

$\min_{{\bm{q}}\in\mathbb{B}^{D}}{{\bm{q}}}^{\top}({{\bm{W}}}^{L_{2}}+\lambda{{\bm{W}}}^{L_{0}}){\bm{q}},$
(6)

where ${\bm{q}}$ is a binary vector of spins. These spins express the unknown
elements in the vector ${\bm{x}}$, using fixed-point arithmetic as follows:

$x_{i}=c_{i}^{\text{min}}+d_{i}\sum_{p=1}^{P}q_{i,p}2^{p-1},\quad P\geq 1.$
(7)

Above, $c_{i}^{\text{min}},d_{i}$ and $P$ are pre-specified constants:
$c_{i}^{\text{min}}$ is the minimal value that can be expressed, $d_{i}$ is a
scaling factor, and $P$ is the number of bits allocated; the variables

$q_{i,p}\in\mathbb{B},\ \ 1\leq i\leq N,\ 1\leq p\leq P$

are the spins, stored in the $NP$ dimensional vector

${\bm{q}}=[q_{1,1},\dots,q_{1,P},q_{2,1},\dots,q_{2,P},\dots,q_{N,1},\dots,q_{N,P}]^{\top}.$
(8)

The above is the most general representation for ${\bm{x}}$. For example,
consider the case where $c_{i}^{\text{min}}=-1$, $d_{i}=1$, and $P=2$. Here,
$x_{i}$ can get one of the values in the set $\\{-1,0,1,2\\}$, corresponding
to the four possible combinations of $q_{i,p}\in\mathbb{B},\ 1\leq p\leq 2$.
Importantly, throughout this paper, we require that the value $0$ is always
contained in the possible values that $x_{i}$ can get; otherwise this feature
would always be considered a non-zero. We shall note that the analysis of the
most general fixed-point representation for ${\bm{x}}$ with $P\geq 3$ requires
us to expand ${\bm{q}}$ by adding one auxiliary spin per entry $x_{i}$ and
revising (6) accordingly; we discuss this in detail in Section 3.2.

###  3.1 QUBO matrix for the squared error term

In this section, we focus only on the minimization of the squared error term
$\|{\bm{A}}{\bm{x}}-{\bm{b}}\|_{2}^{2}$ with respect to ${\bm{x}}$, and show
how to formulate it as a QUBO problem that minimizes the same objective,
however, with respect to the spins ${\bm{q}}$. To do so, we first expand the
$L_{2}$ norm and rewrite it as follows:

$H_{L_{2}}:=\|{\bm{A}}{\bm{x}}-{\bm{b}}\|_{2}^{2}=\sum_{i=1}^{N}\sum_{j=1}^{N}x_{i}x_{j}W_{i,j}^{\text{base,1}}+\sum_{i=1}^{N}x_{i}W_{i,i}^{\text{base,2}}+h^{\text{base}},$
(9)

where ${\bm{W}}^{\text{base},1}$ is a matrix of size $N\times N$ with entries

$\displaystyle W_{i,j}^{\text{base},1}:=\sum_{m=1}^{M}A_{m,i}A_{m,j},\ \ 1\leq
i,j\leq N;$ (10)

the diagonal matrix ${\bm{W}}^{\text{base},2}$ is also of size $N\times N$
with diagonal elements

$\displaystyle W_{i,i}^{\text{base},2}:=-2\sum_{m=1}^{M}A_{m,i}b_{m},\ \ 1\leq
i\leq N;$ (11)

and $h^{\text{base}}$ is a constant scalar, defined as

$\displaystyle h^{\text{base}}:=\sum_{m=1}^{M}b_{m}^{2}.$ (12)

Next, we plug into (9) the explicit fixed-point expression for $x_{i}$ defined
in (7), expressing $H_{L_{2}}$ using the binary spins $q_{i,p}\in\mathbb{B}$
in lieu of $x_{i}$. After applying basic algebraic manipulations, we get the
following expression:

$H_{L_{2}}=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{s=1}^{P}\sum_{p=1}^{P}q_{i,s}q_{j,p}W_{s+P(i-1),p+P(j-1)}^{L_{2},1}+\sum_{i=1}^{N}\sum_{p=1}^{P}q_{i,p}W_{p+P(i-1),p+P(i-1)}^{L_{2},2}+h^{L_{2}}.$
(13)

Above, ${\bm{W}}^{L_{2},1}$ is a matrix of size $NP\times NP$ whose entries
are given by

$\displaystyle
W_{s+P(i-1),p+P(j-1)}^{L_{2},1}:=2^{s+p-2}W_{i,j}^{\text{base},1}d_{i}d_{j},\
\ 1\leq i,j\leq N,\ \ 1\leq s,p\leq P.$ (14)

In addition, ${\bm{W}}^{L_{2},2}$ is an $NP\times NP$ diagonal matrix whose
diagonal elements are expressed as

$\displaystyle
W_{p+P(i-1),p+P(i-1)}^{L_{2},2}:=2^{p-1}d_{i}\left(W_{i,i}^{\text{base},2}+2\sum_{j=1}^{N}c_{j}^{\text{min}}W_{i,j}^{\text{base},1}\right),\
1\leq i\leq N,\ \ \ \ 1\leq p\leq P,$ (15)

and the constant $h^{L_{2}}$ is formulated as

$\displaystyle
h^{L_{2}}:=\sum_{i=1}^{N}\sum_{j=1}^{N}W_{i,j}^{\text{base},1}c_{i}^{\text{min}}c_{j}^{\text{min}}+\sum_{i=1}^{N}W_{i,i}^{\text{base},2}c_{i}^{\text{min}}+h^{\text{base}}.$
(16)

Now, we make two observations. First, since $q_{i,p}$ is binary we have
$q_{i,p}q_{i,p}=q_{i,p}^{2}=q_{i,p}$, and therefore we can add the diagonal
elements of ${\bm{W}}^{L_{2},1}_{ip,ip}$ to those of
${\bm{W}}^{L_{2},2}_{ip,ip}$. Second, the minimization of $H_{L_{2}}$ with
respect to ${\bm{q}}$ is not affected by $h^{L_{2}}$ since the latter is a
constant. These two observations complete our derivation of the QUBO matrix
for the squared error term in (3), concluding that the matrix
${\bm{W}}^{L_{2}}$ in (6) is given by

${\bm{W}}^{L_{2}}={\bm{W}}^{L_{2},1}+{\bm{W}}^{L_{2},2}.$ (17)

###  3.2 QUBO matrix for the cardinality term

We turn to develop an explicit expression for the QUBO matrix that corresponds
to the sparsity penalty. The idea is to compute the cardinality of the
estimated ${\bm{x}}$ under the fixed-point representation, implemented via the
binary vector ${\bm{q}}$. Since we allow $x_{i}$ to have negative values, we
use the following strategy to facilitate the representation of the zero
elements in ${\bm{x}}$. Denote the combination of the binary elements that
leads to $x_{i}=0$ by $c_{i,p}^{0}$, satisfying

$x_{i}=0=c_{i}^{\text{min}}+d_{i}\sum_{p=1}^{P}c_{i,p}^{0}2^{p-1}.$ (18)

Recall our example where we set $c_{i}^{\text{min}}=-1$, $d_{i}=1$, and $P=2$,
in which $x_{i}\in\\{-1,0,1,2\\}$. In this case, the value $x_{i}=0$ is
obtained for $c^{0}_{i,1}=1$ and $c^{0}_{i,2}=0$. Importantly, the binary
elements $c^{0}_{i,p}$ are known constants, as these are derived from the pre-
determined constants $c_{i}^{\text{min}}$, $d_{i}$, and $P$. Given the
constants $c^{0}_{i,p}$, we define

$y_{i,p}=\begin{cases}1-q_{i,p},&c_{i,p}^{0}=0,\\\
q_{i,p},&c_{i,p}^{0}=1,\end{cases}$ (19)

where the binary variables $y_{i,p}\in\mathbb{B}$ can be thought of as
“transformed spins”. Importantly, under the above formulation, we have

$\displaystyle y_{i,p}=1\ \ \forall 1\leq p\leq P\ \ \xleftrightarrow{\ }\
c^{0}_{i,p}=q_{i,p}\ \ \forall 1\leq p\leq P\ \ \xleftrightarrow{\ }\
x_{i}=0,$ (20)

where $\leftrightarrow$ stands for “if and only if”. In words, the above
implies that $x_{i}$ is equal to zero if and only if all the transformed spins
$y_{i,p}$, $1\leq p\leq P$ are ‘active’ and equal to $1$. Compactly, this can
be expressed as

$z_{i}=y_{i,1}y_{i,2}\dots y_{i,P},$ (21)

where $z_{i}=1$ if and only if $x_{i}=0$. Now, we can compute how many non-
zeros we have in ${\bm{x}}$ by

$H_{L_{0}}:=\|{\bm{x}}\|_{0}=\sum_{i=1}^{N}(1-z_{i}).$ (22)

Therefore, the formulation of the transformed spins $y_{i,p}$ allows us to
express the sparsity penalty in QUBO terms, for the most general fixed-point
representation. Unfortunately, in this general case, the variables $z_{i}$ are
not quadratic in the spins $q_{i,p}$, and therefore (22) can not be used in
its present form to construct the matrix ${\bm{W}}^{L_{0}}$ in (6). Yet, we
will show how to overcome this challenge using auxiliary ancilla spins. Before
doing so, however, we pause to discuss two important special cases. The first
is a case where ${\bm{x}}$ is binary with $P=1$ and $d_{i}=1$ in (7), and the
second extends the former to the 2-bit case, with $P=2$. In contrast to the
$P\geq 3$ case, both the binary and 2-bit representation do not require the
use of ancilla spins.

#### Binary representation: $\mathbold{P}$=1

In this simple setting, $x_{i}=q_{i}$ by construction, as $q_{i}$ is the
binary spin that corresponds to the $i$th entry in the binary ${\bm{x}}$.
Therefore, the number of non-zeros in the estimated ${\bm{x}}$ is nothing but
the sum over all the original spins, which can be written as follows:

$H_{L0}^{\text{binary}}=\|{\bm{x}}\|_{0}=\sum_{i=1}^{N}q_{i}=\sum_{i=1}^{N}q_{i}q_{i}={\bm{q}}^{\top}{\bm{W}}^{L_{0}}{\bm{q}}.$
(23)

Above, the third equality holds since $q_{i}\in\mathbb{B}$, and the forth
equality holds for the identity matrix ${\bm{W}}^{L_{0}}={\bm{I}}$ of size
$N\times N$.

#### 2-bit representation: $\mathbold{P}$=2

In contrast to the binary case, in this setting the cardinality of ${\bm{x}}$
is not equal anymore to the sum the original spins $q_{i,p}$, $1\leq i\leq N$,
$1\leq p\leq 2$. Yet, following (22), we can express the sparsity penalty
using the two transformed spins $y_{i,1}$ and $y_{i,2}$, where this
transformation is pre-defined via the constants $c^{0}_{i,p}\in\mathbb{B}$,
$1\leq p\leq 2$. Below, we show that for the special case of $P=2$, the
expression in (22) is a quadratic function of the original spins $q_{i,p}$,
which is crucial to yield the QUBO matrix ${\bm{W}}^{L_{0}}$. The idea is to
map the transformed spins $y_{i,p}$ back to the original space of $q_{i,p}$,
for all $i,p$, and then rewrite the resulting expression in a form of
${\bm{q}}^{\top}{\bm{W}}^{L_{0}}{\bm{q}}$. Concretely, consider our running
example, in which $c_{i}^{\text{min}}=-1$, $d_{i}=1$, and $P=2$. Here, $x_{i}$
can get negative values, however we know that $x_{i}=0$ for the constants
$c^{0}_{i,1}=1$ and $c^{0}_{i,2}=0$. Now, imagine we start with a $2N\times
2N$ matrix ${\bm{W}}^{L_{0}}$ of zeros, and we wish to fill-in the entries in
${\bm{W}}^{L_{0}}$ that correspond to $x_{i}$. According to (19), since
$c_{i,1}^{0}=1,\ c_{i,2}^{0}=0$, we have

$y_{i,1}=q_{i,1}\ \text{and}\ y_{i,2}=1-q_{i,2},$ (24)

leading to $1-z_{i}$ that is quadratic in the spins $q_{i,p}$, since

$1-z_{i}=1-y_{i,1}y_{i,2}=1-q_{i,1}(1-q_{i,2})=1-q_{i,1}q_{i,1}+q_{i,1}q_{i,2}.$
(25)

The above implies that we should set the corresponding entries in the QUBO
matrix as follows:

$W^{L_{0}}_{1+2(i-1),1+2(i-1)}=-1;\
W^{L_{0}}_{1+2(i-1),2+2(i-1)}=\frac{1}{2};\ \text{and}\
W^{L_{0}}_{2+2(i-1),1+2(i-1)}=\frac{1}{2}.$ (26)

Note that we ignore the leading constant ‘1’ in the right-hand-side of (25)
since constants do not affect the solution of the underlying optimization
problem. The remaining elements in the matrix can be determined by repeating a
similar set of steps for all $1\leq i\leq N$, however with possibly different
values of $c^{0}_{i,1},c^{0}_{i,2}$. Notice that there are four cases for
$c^{0}_{i,p}\in\mathbb{B}$, $1\leq p\leq 2$ that should be considered, while
above we studied only one of these for which $c^{0}_{i,1}=1$ and
$c^{0}_{i,2}=0$. For completeness, we present below the three remaining cases.
For $c_{i,1}^{0}=0,\ c_{i,2}^{0}=0$, we follow (19) and get
$y_{i,1}=1-q_{i,1}\ \text{and}\ y_{i,2}=1-q_{i,2}$, leading to
$1-z_{i}=q_{i,1}+q_{i,2}-q_{i,1}q_{i,2}$. Analogously, when $c_{i,1}^{0}=0,\
c_{i,2}^{0}=1$ we have $y_{i,1}=1-q_{i,1}\ \text{and}\ y_{i,2}=q_{i,2}$,
leading to $1-z_{i}=1-q_{i,2}+q_{i,1}q_{i,2}$. The last case deals with
$q_{i,1}^{0}=1,\ q_{i,2}^{0}=1$; here, $y_{i,1}=q_{i,1}\ \text{and}\
y_{i,2}=q_{i,2}$, and thus $1-z_{i}=1-q_{i,1}q_{i,2}$. It is immediate to
explicitly express ${\bm{W}}^{L_{0}}$ for each case as in (26); we omit this
in the interest of space.

#### The general case: $\mathbold{P\geq}$ 3

This is the most difficult case to address, as a naive extension of (21) to
the $P$-bit case would result in $H_{L_{0}}$ in (22) that has high-order
interactions between the original spins, breaking the bilinear structure of
the QUBO problem. As a way out, we introduce ancilla spins, which allow us to
express $H_{L_{0}}$ for the most general representation of ${\bm{x}}$ in a
quadratic form that perfectly fits the QUBO structure. We take inspiration
from [49, 50] and offer a solution that is extremely efficient with respect to
the number of ancilla spins: we add only one auxiliary spin per feature, which
we view as the minimal number of spins that one can hope for, in such a
general case.

Turning to the details, denote the ancilla spin for $x_{i}$ by
$s_{i}\in\mathbb{B}$ and define the function

$F(y_{i,1},y_{i,2},\dots,y_{i,P},s_{i}):=s_{i}\cdot\left(y_{i,1}+y_{i,2}+\dots+y_{i,P}-(P-1)\right),$
(27)

which gets as input all the transformed spins $y_{i,p}$ as well as the ancilla
spin $s_{i}$ of $x_{i}$. Now, we invoke a beautiful result presented in [49]
and [50, Section 4.2], stating that

$-z_{i}=-y_{i,1}y_{i,2}\dots
y_{i,P}=\min_{s_{i}}-F(y_{i,1},y_{i,2},\dots,y_{i,P},s_{i}).$ (28)

In plain words, the minimal value of the function
$-F(y_{i,1},y_{i,2},\dots,y_{i,P},s_{i})$, optimized with respect to the
ancilla spin $s_{i}$, is in fact equal to $-z_{i}$. Since the function $F$ in
(27) is a sum of $s_{i}y_{i,p}$, $1\leq p\leq P$, we can harness it to form a
cost function that is equivalent to $H_{L_{0}}$ in (22), but involves only
bilinear spin terms. Specifically, we define our cardinality Hamiltonian as

$\displaystyle H_{L_{0}}=\|{\bm{x}}\|_{0}=\sum_{i=1}^{N}(1-z_{i})$
$\displaystyle=\sum_{i=1}^{N}(1-F(y_{i,1},y_{i,2},\dots,y_{i,P}))$ (29)
$\displaystyle=\sum_{i=1}^{N}\left(1-s_{i}y_{i,1}-s_{i}y_{i,2}-\dots-
s_{i}y_{i,P}+s_{i}\left(P-1\right)\right).$ (30)

Armed with the above cost function, we are able to form the structure and
content of the matrix ${\bm{W}}^{L_{0}}$. One technical challenge is that
$y_{i,p}$ are the transformed versions of the original spins $q_{i,p}$, and so
we must substitute the spins $q_{i,p}$ directly into (29). This step is
similar in spirit to the 2-bit case that we have already discussed in depth
earlier. By doing so, we obtain an expression that is quadratic with respect
to the $NP$ original spins $q_{i,p}$, as well as with the additional $N$
ancilla spins $s_{i}$. Concretely, we define the vector
$\tilde{{\bm{q}}}=[{\bm{q}}\ ;\ {\bm{s}}]^{\top}$ that contains a total of
$N(P+1)$ spins; the first $NP$ elements are the original spins ${\bm{q}}$, and
the rest are the $N$ ancilla spins

${\bm{s}}=[s_{1},s_{2},\dots,s_{N}].$

Consequently, the matrix ${\bm{W}}^{L_{0}}$ is of size $N(P+1)\times N(P+1)$,
where we use the convention that the first $NP$ rows (resp. columns)
correspond to the spins $q_{i,p}$ and the rest $N$ rows (resp. columns)
correspond to the ancilla spins $s_{i}$. Here, the QUBO problem defined with
the vector of spins $\tilde{{\bm{q}}}$ is given by

$\min_{\tilde{{\bm{q}}}\in\mathbb{B}^{d}}\tilde{{\bm{q}}}^{\top}(\widetilde{{\bm{W}}}^{L_{2}}+\lambda{\bm{W}}^{L_{0}})\tilde{{\bm{q}}},$

where $\widetilde{{\bm{W}}}^{L_{2}}$ is a matrix of size $N(P+1)\times
N(P+1)$, obtained by padding the matrix ${\bm{W}}^{L_{2}}$ in (17) with zeros
in all the entries that correspond to the ancilla spins.

To better understand how to construct ${\bm{W}}^{L_{0}}$, it may be best to
consider an example with $P=4$. Observe that each of the four transformed
spins $y_{i,1},y_{i,2},y_{i,3},y_{i,4}$ can get one of two possible
expressions, determined by the constants $c_{i,1}^{0}$ and $c_{i,2}^{0}$,
$c_{i,3}^{0}$, $c_{i,4}^{0}$, respectively. This creates 16 different explicit
forms for the term $F(y_{i,1},y_{i,2},y_{i,3},y_{i,4},s_{i})$ as a function of
$c_{i,1}^{0},c_{i,2}^{0},c_{i,3}^{0},c_{i,4}^{0}$. For instance, following
(19), the choice of $c_{i,1}^{0}=1,c_{i,2}^{0}=0,c_{i,3}^{0}=0,c_{i,4}^{0}=1$
implies that $y_{i,1}=q_{i,1}\ y_{i,2}=1-q_{i,2},\ y_{i,3}=1-q_{i,3},\
y_{i,4}=q_{i,4}$. Plugging the latter into the explicit expression of $F$
results in

$\displaystyle 1-z_{i}=1-F(y_{i,1},y_{i,2},y_{i,3},y_{i,4},s_{i})$
$\displaystyle=1-F(q_{i,1},1-q_{i,2},1-q_{i,3},q_{i,4},s_{i,1})$ (31)
$\displaystyle=1-s_{i}(q_{i,1}+(1-q_{i,2})+(1-q_{i,3})+q_{i,4}-(4-1))$ (32)
$\displaystyle=1-s_{i}(q_{i,1}-q_{i,2}-q_{i,3}+q_{i,4}-1).$ (33)

Observe that the minimal value of (33) is 0, obtained for the combination
$q_{i,1}=c_{i,1}^{0}=1,q_{i,2}=c_{i,2}^{0}=0,q_{i,3}=c_{i,3}^{0}=0$,
$q_{i,4}=c_{i,4}^{0}=1$, and $s_{i}=1$. Now, since the leading constant ‘1’ in
(33) does not affect the combination of spins that minimizes this expression,
we can set the entries in ${\bm{W}}^{L_{0}}$ that correspond to (33) as
follows:

$W^{L_{0}}_{1+4(i-1),N+i}=-1,\ W^{L_{0}}_{2+4(i-1),N+i}=1,\
W^{L_{0}}_{3+4(i-1),N+i}=1,\ W^{L_{0}}_{4+4(i-1),N+i}=-1,\ \text{and}\
W^{L_{0}}_{N+i,N+i}=+1.$

A similar analysis can be done for the remaining 15 cases, which we do not
provide here in the interest of space.

In closing, this section introduced a novel decomposition of the sparse coding
problem (3) as a QUBO model. To formulate an efficient solution with respect
to the number of spins, our analysis is separated into three cases: binary
with $P=1$, 2-bit with $P=2$, and a general one with $P\geq 3$, where the
total number of spins varies with $P$. For $P=1$ we have $N$ spins, for $P=2$
we have $2N$ spins, and for the general $P\geq 3$ case we have $N(P+1)$ spins
in total.

##  4 Experiments

(a) Cardinality $k=30$, noise $\sigma=0.1$.

(b) Cardinality $k=30$, rows $M=80$.

(c) Rows $M=80$, noise $\sigma=0.1$.

Figure 2: Results obtained by LightSolver’s digital simulator, minimizing the
proposed QUBO problem for a binary ${\bm{x}}$ of length $N=160$. Each point in
the graphs is an average over 20 independent realizations of
${\bm{A}},{\bm{x}}$, and ${\bm{v}}$. Other details are as in Figure 1.

##### Binary sparse vectors.

We return to the experiment from Section 2, but now increase significantly the
dimensions of the problem. We scale-up the dimensions of ${\bm{A}}$ and
${\bm{x}}$ by a factor of 10, focusing on a sparse _binary_ vector ${\bm{x}}$
($P=1$) of dimension $N=160$. The matrix ${\bm{A}}$ is designed to have low
mutual coherence according to (5), generated as described in Appendix A. Since
the vector ${\bm{x}}$ is binary, we use the QUBO formulation with $P=1$, where
the total number of spins is equal to $N=160$. Armed with ${\bm{W}}^{L_{2}}$
and ${\bm{W}}^{L_{0}}$, we minimize (6) using a quantum-inspired annealer
implemented by LightSolver’s digital simulator. We then compare the quality of
our quantum-inspired estimation $\hat{{\bm{x}}}$ to lasso and OMP, serving as
baseline methods. Naturally, we cannot include Algorithm 1 as a baseline,
since it is infeasible to run such an exhaustive optimization method with the
dimensions of the problem that we study here. Similarly to Section 2, the
hyper-parameter $\lambda$ of our method, as well as the hyper-parameters of
lasso and OMP, are chosen by an “orcale” that has access to ${\bm{x}}$,
seeking the parameter that achieves the smallest possible error metric under
study, for each method. Figure 2 presents the reconstruction and support
recovery errors as a function of $M$ (Figure 2(a)), $\sigma$ (Figure 2(b)),
and $k$ (Figure 2(c)). Observe how the quantum-inspired estimations are highly
accurate, portraying nearly perfect recovery in all cases. In particular, the
quantum-based estimations are superior than those obtained by the widely-used
lasso and OMP, all across the board. This dramatic improvement can be
explained as follows. First, the annealer accurately minimizes the NP-hard
QUBO problem of interest. Second, lasso and OMP do not leverage the prior
knowledge that ${\bm{x}}$ is binary, in contrast to the annealer that utilizes
this knowledge, encapsulated apriori in the QUBO matrices.

##### 2-bit sparse vectors.

We now turn to expand the above experiment, by considering a more flexible
representation of ${\bm{x}}$ with $P=2$. Recall that our QUBO formulation for
the 2-bit setting requires $2N$ spins in total, in contrast to the binary case
that requires only $N$ spins. Therefore, we consider a smaller matrix
${\bm{A}}$ with $N=80$ columns, keeping the total number of spins identical
across the two experiments. The 2-bit representation of ${\bm{x}}$ follows
(7), with $c_{i}^{\text{min}}=0$ and $d_{i}=1$; consequently, each _non-zero_
entry is randomly sampled from the discrete set $\\{1,2,3\\}$. Figure 3
displays the reconstruction and support recovery errors as a function of the
sample-size $M$, while fixing the cardinality ($k=10$) and noise
($\sigma=0.1$) levels. As portrayed, the annealer achieves impressive
estimation errors, outperforming OMP and lasso especially in the most
challenging setting in which the sample size is relatively small.

Figure 3: Results obtained by LightSolver’s digital simulator, minimizing the
proposed QUBO problem for a 2-bit ${\bm{x}}$ of length $N=80$. We fixed the
cardinality to $k=10$ and set the noise level to $\sigma=0.1$. Error metrics
are evaluated over $20$ independent realizations of ${\bm{A}},{\bm{x}}$ and
${\bm{v}}$. Other details are as in Figure 1.

##  5 Conclusions and future directions

In this work, we formulate the ubiquitous sparse coding task as a quadratic
unconstrained binary optimization (QUBO) problem, which, potentially, can be
minimized efficiently using quantum computers and Ising machines. Numerical
experiments demonstrate the superiority of our quantum-based estimations
compared to widely-used, classic sparse approximation algorithms—lasso and
OMP. In particular, we report significant improvements in challenging regimes
where (i) the sample-size is relatively small, (ii) the noise-level is
relatively large, and (iii) the cardinality is relatively large. While our
QUBO formulation can handle the most general fixed-point representation of the
underlying sparse vector, we believe it would be of great interest to extend
the proposed framework to a floating point representation. Such a formulation
may reduce the number of spins required to express the solution with a
sufficient precision, which is of paramount importance given the challenge of
quantum/Ising platforms to increase the number of spins. We also believe it
would be exciting to expand our work beyond regression settings, and provide
QUBO formulations to sparse-regularized classification algorithms.

In a broader view, there is a growing interest in using quantum technologies
in biological sciences [51], and the methods presented in this paper may
provide performance improvements in applications where sparse linear
regression plays a role. As an example, in genome-wide association studies
(GWAS) scientists are interested in accurately identifying which of the
thousands of single-nucleotide polymorphisms (SNPs) are associated with high
cholesterol, or any other phenotype of interest [52, 53]. The ability to find
a sparse, interpretable set of ‘important’ SNPs can be used to improve medical
treatment as well as to expand our understanding of the human genome. In these
applications, it is common for the number of observations (e.g., individuals
participating in a study) to be smaller than the number of features (SNPs), as
biotech development allows a massive collection of genetic variants for each
individual [54]. The desire to form an interpretable predictive rule together
with the high-dimensional nature of the data—the number of SNPs (columns of
$A$) is larger than the observations (rows of $A$)—may explain the popularity
of sparse linear regression methods in this field. Although conducted on
simulated data, the experiments presented in this paper indicate that our
quantum-based approach thrives in ultra high-dimensional regimes, and
therefore we believe it would be exciting to explore the potential use of our
quantum sparse coding method in this domain.

## References

  * [1] Emmanuel J Candes, Justin K Romberg, and Terence Tao.  Stable signal recovery from incomplete and inaccurate measurements.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 59(8):1207–1223, 2006.
  * [2] David L Donoho.  Compressed sensing.  IEEE Transactions on information theory, 52(4):1289–1306, 2006\.
  * [3] Michael Lustig, David Donoho, and John M Pauly.  Sparse MRI: The application of compressed sensing for rapid MR imaging.  Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, 58(6):1182–1195, 2007\.
  * [4] Kjersti Engan, Sven Ole Aase, and J Hakon Husoy.  Method of optimal directions for frame design.  In International Conference on Acoustics, Speech, and Signal Processing. Proceedings, volume 5, pages 2443–2446. IEEE, 1999.
  * [5] Michal Aharon, Michael Elad, and Alfred Bruckstein.  K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation.  IEEE Transactions on signal processing, 54(11):4311–4322, 2006\.
  * [6] Julien Mairal, Francis Bach, Jean Ponce, and Guillermo Sapiro.  Online dictionary learning for sparse coding.  In Proceedings of the 26th annual international conference on machine learning, pages 689–696, 2009.
  * [7] Michael Elad and Michal Aharon.  Image denoising via sparse and redundant representations over learned dictionaries.  IEEE Transactions on Image processing, 15(12):3736–3745, 2006.
  * [8] Weisheng Dong, Xin Li, Lei Zhang, and Guangming Shi.  Sparsity-based image denoising via dictionary learning and structural clustering.  In CVPR 2011, pages 457–464. IEEE, 2011.
  * [9] Bin Shen, Wei Hu, Yimin Zhang, and Yu-Jin Zhang.  Image inpainting via sparse representation.  In 2009 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 697–700. IEEE, 2009.
  * [10] Mohamed-Jalal Fadili, J-L Starck, and Fionn Murtagh.  Inpainting and zooming using sparse representations.  The Computer Journal, 52(1):64–79, 2009.
  * [11] Jeremias Sulam, Boaz Ophir, Michael Zibulevsky, and Michael Elad.  Trainlets: Dictionary learning in high dimensions.  IEEE Transactions on Signal Processing, 64(12):3180–3193, 2016\.
  * [12] Jianchao Yang, John Wright, Thomas S Huang, and Yi Ma.  Image super-resolution via sparse representation.  IEEE transactions on image processing, 19(11):2861–2873, 2010.
  * [13] Roman Zeyde, Michael Elad, and Matan Protter.  On single image scale-up using sparse-representations.  In International conference on curves and surfaces, pages 711–730. Springer, 2010.
  * [14] Yaniv Romano, Matan Protter, and Michael Elad.  Single image interpolation via adaptive nonlocal sparsity-based modeling.  IEEE Transactions on Image Processing, 23(7):3085–3098, 2014.
  * [15] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.
  * [16] Shaobing Chen and David Donoho.  Basis pursuit.  In Proceedings of 1994 28th Asilomar Conference on Signals, Systems and Computers, volume 1, pages 41–44. IEEE, 1994.
  * [17] Scott Shaobing Chen, David L Donoho, and Michael A Saunders.  Atomic decomposition by basis pursuit.  SIAM review, 43(1):129–159, 2001.
  * [18] Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani.  Least angle regression.  The Annals of statistics, 32(2):407–499, 2004.
  * [19] Michael Elad.  Sparse and redundant representations: from theory to applications in signal and image processing, volume 2.  Springer, 2010.
  * [20] Geoff Davis, Stephane Mallat, and Marco Avellaneda.  Adaptive greedy approximations.  Constructive approximation, 13(1):57–98, 1997.
  * [21] Balas Kausik Natarajan.  Sparse approximate solutions to linear systems.  SIAM journal on computing, 24(2):227–234, 1995.
  * [22] Ingrid Daubechies, Michel Defrise, and Christine De Mol.  An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 57(11):1413–1457, 2004.
  * [23] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM journal on imaging sciences, 2(1):183–202, 2009.
  * [24] Patrick L Combettes and Valérie R Wajs.  Signal recovery by proximal forward-backward splitting.  Multiscale modeling & simulation, 4(4):1168–1200, 2005.
  * [25] Stéphane G Mallat and Zhifeng Zhang.  Matching pursuits with time-frequency dictionaries.  IEEE Transactions on signal processing, 41(12):3397–3415, 1993\.
  * [26] Yagyensh Chandra Pati, Ramin Rezaiifar, and Perinkulam Sambamurthy Krishnaprasad.  Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition.  In Proceedings of 27th Asilomar conference on signals, systems and computers, pages 40–44. IEEE, 1993.
  * [27] Joel A Tropp.  Greed is good: Algorithmic results for sparse approximation.  IEEE Transactions on Information theory, 50(10):2231–2242, 2004\.
  * [28] Joel A Tropp.  Just relax: Convex programming methods for identifying sparse signals in noise.  IEEE transactions on information theory, 52(3):1030–1051, 2006\.
  * [29] Emmanuel J Candes and Terence Tao.  Decoding by linear programming.  IEEE transactions on information theory, 51(12):4203–4215, 2005\.
  * [30] Andrew Lucas.  Ising formulations of many NP problems.  Frontiers in physics, page 5, 2014.
  * [31] Ken A Dill and Justin L MacCallum.  The protein-folding problem, 50 years on.  science, 338(6110):1042–1046, 2012.
  * [32] Satu Elisa Schaeffer.  Graph clustering.  Computer science review, 1(1):27–64, 2007.
  * [33] Robert Patton, Catherine Schuman, Thomas Potok, et al.  Efficiently embedding QUBO problems on adiabatic quantum computers.  Quantum Information Processing, 18(4):1–31, 2019.
  * [34] Prasanna Date, Davis Arthur, and Lauren Pusey-Nazzaro.  QUBO formulations for training machine learning models.  Scientific reports, 11(1):1–10, 2021.
  * [35] Johannes S Otterbach, Riccardo Manenti, Nasser Alidoust, A Bestwick, M Block, B Bloom, S Caldwell, N Didier, E Schuyler Fried, S Hong, et al.  Unsupervised machine learning on a hybrid quantum computer.  arXiv preprint arXiv:1712.05771, 2017.
  * [36] Nicholas Gao, Max Wilson, Thomas Vandal, Walter Vinci, Ramakrishna Nemani, and Eleanor Rieffel.  High-dimensional similarity search with quantum-assisted variational autoencoder.  In Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining, pages 956–964, 2020.
  * [37] Idan Meirzada, Assaf Kalinski, Dov Furman, Tsafrir Armon, Talya Vaknin, Harel Primack, Chene Tradonsky, and Ruti Ben-Shlomi.  Lightsolver–a new quantum-inspired solver cracks the 3-regular 3-xorsat challenge.  arXiv preprint arXiv:2207.09517, 2022.
  * [38] Fred Glover.  Tabu search—part I.  ORSA Journal on computing, 1(3):190–206, 1989.
  * [39] Fred Glover.  Tabu search—part II.  ORSA Journal on computing, 2(1):4–32, 1990.
  * [40] Kyungtaek Jun and Kwangmin Yu.  QUBO formulations for system of linear equations.  arXiv preprint arXiv:2106.10819, 2021.
  * [41] Paul S Bradley and Olvi L Mangasarian.  Feature selection via concave minimization and support vector machines.  In ICML, volume 98, pages 82–90, 1998.
  * [42] Robert S Wezeman, Irina Chiscop, Laura Anitori, and Wim van Rossum.  Quantum-classical solution methods for binary compressive sensing problems.  In International Conference on Computational Science, pages 107–121. Springer, 2022.
  * [43] Ramin Ayanzadeh, Seyedahmad Mousavi, Milton Halem, and Tim Finin.  Quantum annealing based binary compressive sensing with matrix uncertainty.  arXiv preprint arXiv:1901.00088, 2019.
  * [44] Ramin Ayanzadeh, Milton Halem, and Tim Finin.  An ensemble approach for compressive sensing with quantum annealers.  In IEEE International Geoscience and Remote Sensing Symposium, pages 3517–3520. IEEE, 2020.
  * [45] Richard Y Li, Rosa Di Felice, Remo Rohs, and Daniel A Lidar.  Quantum annealing versus classical machine learning applied to a simplified computational biology problem.  NPJ quantum information, 4(1):1–10, 2018.
  * [46] Naoki Ide and Masayuki Ohzeki.  Sparse signal reconstruction with QUBO formulation in l0-regularized linear regression.  arXiv preprint arXiv:2202.00452, 2022.
  * [47] Maurizio Ferrari Dacrema, Fabio Moroni, Riccardo Nembrini, Nicola Ferro, Guglielmo Faggioli, and Paolo Cremonesi.  Towards feature selection for ranking and classification exploiting quantum annealers.  In Proceedings of the 45th International ACM SIGIR Conference on Research and Development in Information Retrieval, pages 2814–2824, 2022.
  * [48] David L Donoho, Michael Elad, and Vladimir N Temlyakov.  Stable recovery of sparse overcomplete representations in the presence of noise.  IEEE Transactions on information theory, 52(1):6–18, 2005.
  * [49] 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.
  * [50] 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, 2010.
  * [51] Prashant S Emani, Jonathan Warrell, Alan Anticevic, Stefan Bekiranov, Michael Gandal, Michael J McConnell, Guillermo Sapiro, Alán Aspuru-Guzik, Justin T Baker, Matteo Bastiani, et al.  Quantum computing at the frontiers of biological sciences.  Nature Methods, 18(7):701–709, 2021.
  * [52] M Graziano Usai, Mike E Goddard, and Ben J Hayes.  LASSO with cross-validation for genomic selection.  Genetics research, 91(6):427–436, 2009.
  * [53] Patrik Waldmann, Gábor Mészáros, Birgit Gredler, Christian Fuerst, and Johann Sölkner.  Evaluation of the lasso and the elastic net in genome-wide association studies.  Frontiers in genetics, 4:270, 2013.
  * [54] Matteo Sesia, Stephen Bates, Emmanuel Candès, Jonathan Marchini, and Chiara Sabatti.  False discovery rate control in genome-wide association studies with population structure.  Proceedings of the National Academy of Sciences, 118(40):e2105841118, 2021.

##  Appendix A Generating a rectangular matrix with low mutual coherence

Algorithm 2 presents a heuristic iterative procedure for generating a matrix
${\bm{A}}$ whose mutual coherence $\mu({\bm{A}})$ is minimized. The idea
behind this procedure is to design an ${\bm{A}}$ such that the gram matrix
${\bm{A}}^{\top}{\bm{A}}\in\mathbb{R}^{N\times N}$ will be as close as
possible to the identity matrix.

Algorithm 2 Heuristic gradient decent algorithm for generating ${\bm{A}}$ with
low mutual coherence.

Number of rows $M$, number of columns $N$, step size $\alpha$, number of
iterations $K$.

A matrix ${\bm{A}}$ with normalized columns such that $\mu({\bm{A}})$ defined
in (5) is minimized.

Initialize ${\bm{A}}_{1}$ with entries drawn from a Standard Normal
distribution.

Normalize ${\bm{A}}_{1}$ such each column $a_{i}$, $1\leq i\leq N$, has a unit
$L_{2}$ norm.

for $k=1,\dots,K$ do

${\bm{A}}_{k+1}\leftarrow
A_{k}-\alpha\cdot\nabla_{A_{k}}\|A_{k}^{\top}A_{k}-I_{N\times N}\|_{F}^{2}$,
where $I_{N\times N}$ is the $N\times N$ identity matrix.

Normalize ${\bm{A}}_{k+1}$ to have columns with a unit $L_{2}$ norm.

end for