# Quantum Matching Pursuit: a Quantum Algorithm for Sparse Representations

Armando Bellante  Politecnico di Milano, DEIB, Via Ponzio 34/5 – Building 20,
Milan 20133, Italy<EMAIL_ADDRESS>Stefano Zanero  Politecnico di
Milano, DEIB, Via Ponzio 34/5 – Building 20, Milan 20133, Italy.

(August 30, 2024)

###### Abstract

Representing signals with sparse vectors has a wide range of applications that
range from image and video coding to shape representation and health
monitoring. In many applications with real-time requirements, or that deal
with high-dimensional signals, the computational complexity of the encoder
that finds the sparse representation plays an important role. Quantum
computing has recently shown promising speed-ups in many representation
learning tasks. In this work, we propose a quantum version of the well-known
matching pursuit algorithm. Assuming the availability of a fault-tolerant
quantum random access memory, our quantum matching pursuit lowers the
complexity of its classical counterpart of a polynomial factor, at the cost of
some error in the computation of the inner products, enabling the computation
of sparse representation of high-dimensional signals. Besides proving the
computational complexity of our new algorithm, we provide numerical
experiments that show that its error is negligible in practice. This work
opens the path to further research on quantum algorithms for finding sparse
representations, showing suitable quantum computing applications in signal
processing.

††preprint: APS/123-QED

##  I Introduction

Finding a sparse representation is the problem of representing a dense signal
as a linear combination of a few unit vectors, also referred to as _atoms_.
Usually, the set of atoms is larger than the space where the signal lies, as
over-complete sets of atoms enable sparser representations [1]. Once a set of
atoms, or _dictionary_ , is fixed, the sparse representation of the signal is
the set of coefficients of their linear combination.

Signals of the same type are likely to be represented sparsely over the same
dictionary. For instance, the widely-used JPEG algorithm exploits the fact
that images are sparse with respect to the Discrete Cosine Transform basis to
perform compression [2]. Finding sparse representations is a subject of
interest in many fields, and its applications range from data compression to
denoising and anomaly detection [3, 4].

When these applications have real-time requirements, or deal with high-
dimensional signals, the computational cost of finding the representation is
crucial. Unfortunately, finding the sparsest representation that approximates
the signal is an NP-hard problem, and is intractable in practice. For this
reason, researchers have developed a series of greedy algorithms that, through
local optimizations, find approximate solutions in an acceptable running time.

In recent years, the effectiveness of quantum computing in representation
learning has become increasingly evident. Recent research has proven
computational advantages for algorithms such as principal component analysis
[5], slow feature analysis [6], and spectral clustering [7].

In this work, we propose an end-to-end quantum algorithm for learning sparse
representations using a matching pursuit approach. We develop a quantum
version of the well-known matching pursuit algorithm [8], reducing its
computational cost of a polynomial factor. While there are some known speed-
ups in the case of specific analytical dictionaries [9], to our knowledge,
there are no classical algorithms that compare with our run-time over a
general dictionary.

Besides thoroughly analyzing the running time and error of our novel
algorithm, we describe a suitable quantum processing framework that can be
used as a starting point to construct other quantum pursuit algorithms.

The remainder of the paper is organized as follows. In Section II, we discuss
previous work that relates to ours. Section III describes the classical
algorithm for matching pursuit, introducing the necessary notation for both
the quantum and the classical versions. In Section IV, we briefly introduce
the concept of quantum computation and some subroutines that will serve as a
basis for the novel quantum algorithm. Section V presents the novel quantum
matching pursuit algorithm, providing a thorough run-time and error analysis.
Finally, in Section VI, we run numerical experiments to show that the quantum
matching pursuit can find representations that are as sparse as those of its
classical counterpart.

##  II Related work

The matching pursuit algorithm was first introduced by Mallat and Zhang [8].
The original version of the algorithm has a running time of $O(knm)$, where
$k$ is the number of optimization iterations, $n$ is the length of the signal,
and $m$ is the number of atoms in the dictionary.

Among the attempts to speed up the matching pursuit algorithm, the closest to
ours is the one of Krstulovic and Gribonval [9]. They exploited particular
properties of some analytic (non-learned) dictionaries, like the multiscale
time-frequency Gabor dictionary [8], to reduce the run-time of matching
pursuit to $O(kn\log(n))$. However, their algorithm is still slow on non-
analytical dictionaries.

While some previous works suggest the use of matching pursuit to simulate the
dynamics of quantum mechanical processes [10, 11, 12, 13], to our knowledge
there is no previous work that discusses quantum speed-ups for finding sparse
representations of signals over large dictionaries.

##  III Classical Matching Pursuit

###  III.1 Notation

We denote matrices using capital letters and use lowercase letters for vectors
and scalars. Given a matrix $A$, its $i^{th}$ row and column are denoted by
$a_{i,\cdot}$ and $a_{i}$ respectively. The component identified by the
$i^{th}$ row and the $j^{th}$ column is denoted $a_{ij}$. We write the
$j^{th}$ element of a vector $u$ as $u_{j}$.

Let $x\in\mathbb{C}^{n}$ be a unit vector. Using Dirac’s notation, we use
$\ket{x}$ to represent it as a column vector and $\bra{x}$ to denote its
complex-conjugate row vector. We use $\braket{a_{i},b_{j}}$ to denote the
inner product between two vectors $a_{i},b_{j}$.

The notation $\left\lVert\cdot\right\rVert_{2}$ indicates the Euclidean norm
of a vector. The pseudo-norm $\left\lVert\cdot\right\rVert_{0}$ is the number
of non-zero components of a vector. The symbol
$\left\lVert\cdot\right\rVert_{F}$ indicates the Frobenius norm of a matrix.
For a matrix $A\in\mathbb{R}^{n\times m}$, the Frobenius norm is defined as
$\left\lVert A\right\rVert_{F}=\sqrt{\sum_{i}^{n}\sum_{j}^{m}a_{ij}^{2}}$.

When stating the complexity of an algorithm, the $\widetilde{O}()$ notation
omits poly-logarithmic terms in the input data size (e.g., if an algorithm
uses a matrix $A\in\mathbb{R}^{n\times m}$, $\widetilde{O}(1)\equiv
O(\text{polylog}(nm))$).

###  III.2 Problem statement

We can represent a signal as a vector $s\in\mathbb{R}^{n}$. We use $d_{j}$ to
denote the $j^{th}$ atom over which we search the sparse representation. Each
atom is a unit vector, meaning that for every $j$ we have $\left\lVert
d_{j}\right\rVert_{2}=1$. A dictionary is a matrix $D\in\mathbb{R}^{n\times
m}$ whose columns are the atoms $d_{j}$, for $j\in\\{0,\dots,m-1\\}$. In most
of the interesting cases, the dictionary is over-complete (i.e. $m>n$).

Formally, given a signal $s\in\mathbb{R}^{n}$ and a dictionary
$D\in\mathbb{R}^{n\times m}$, the problem of finding a sparse representation
$x\in\mathbb{R}^{m}$ of the signal is known as $\mathcal{P}_{0}^{\epsilon}$.

######  Definition 1 (Problem $\mathcal{P}_{0}^{\epsilon}$).

Given $s\in\mathbb{R}^{n},D\in\mathbb{R}^{n\times m}$, and
$\epsilon\in\mathbb{R}^{+}$, problem $\mathcal{P}_{0}^{\epsilon}$ is defined
as:

$\displaystyle\operatorname*{arg\,min}_{x}\left\lVert
x\right\rVert_{0}\text{such that}\left\lVert Dx-s\right\rVert_{2}\leq\epsilon$
(1)

Finding the exact solution to problem $\mathcal{P}_{0}^{\epsilon}$ is an NP-
hard task [14]. While quantum computers are not expected to solve NP-hard
problems in polynomial time (indeed, it is widely believed that
$\textsc{NP}\nsubseteq\textsc{BQP}$ [15]), they can still provide speed-ups of
practical use on greedy algorithms that compute approximate solutions. In this
paper we propose a quantum version of the matching pursuit algorithm, a greedy
approach to approximately solve the $\mathcal{P}_{0}^{\epsilon}$ problem in
polynomial time.

###  III.3 Algorithm

The strategy behind the matching pursuit algorithm is to face the problem
through subsequent optimization steps. Starting from an empty solution
$x=0^{\otimes m}$, the matching pursuit searches for the atom that best
reduces the difference between the representation $Dx$ and the signal at each
iteration, updating the solution iteratively. We now discuss the matching
pursuit algorithm in detail.

As an initialization step, we create a residual vector and an empty solution

$\displaystyle r=s,$ (2) $\displaystyle x=0^{\otimes m}.$ (3)

Since we are at the beginning of the algorithm and
$s-Dx=s$, the residual is set equal to the signal.

Once the initialization is complete, the algorithm searches for the closest
atom to the residual by computing

$\displaystyle j^{*}=\operatorname*{arg\,min}_{j}\left\lVert
r-z_{j}d_{j}\right\rVert_{2},$ (4)

where $z_{j}\in\mathbb{R}$ is the best scaling factor for the atom $d_{j}$

$\displaystyle z_{j}=\operatorname*{arg\,min}_{z}\left\lVert
r-zd_{j}\right\rVert.$ (5)

It is possible to show that $z_{j}=\braket{r,d_{j}}$ [8], from which we can
derive the following equivalence:

$\displaystyle\left\lVert r-z_{j}d_{j}\right\rVert_{2}^{2}=\left\lVert
r\right\rVert_{2}^{2}-|\braket{r,d_{j}}|^{2}$ (6)

Because of Equation 6, finding the best atom (Equation 4) is equivalent to
searching for the maximum absolute value of the inner products between the
current residual and the atoms

$\displaystyle j^{*}=\operatorname*{arg\,max}_{j}|\braket{r,d_{j}}|.$ (7)

This step is known as the _sweep stage_ , as we need to iterate over all the
atoms in the dictionary to compute the inner products and choose the best one.

After selecting the best atom, both the solution and the residual get updated

$\displaystyle x_{j^{*}}=x_{j^{*}}+z_{j^{*}},$ (8) $\displaystyle
r=r+z_{j^{*}}d_{j^{*}}.$ (9)

Updating the residual makes it so that the algorithm does not consider the
part of the signal that has been modeled so far. At each iteration, the
residual is $r=s-Dx$.

The algorithm continuously searches for the best approximating atom and
performs the updates until the following stopping condition is met

$\displaystyle\left\lVert x\right\rVert_{0}>L\text{ or }\left\lVert
r\right\rVert_{2}\leq\epsilon,$ (10)

for a sparsity threshold $L\in\mathbb{N}^{+}$ and an error reconstruction
tolerance $\epsilon\in\mathbb{R}^{+}$. We remark that, at each iteration, the
norm of the residual $\left\lVert r\right\rVert_{2}=\left\lVert
s-Dx\right\rVert_{2}$ expresses how well our solution approximates the
original signal.

We summarize this procedure in Algorithm 1.

Algorithm 1 Matching pursuit

1:Initialize $r=s$, $x=0^{\otimes m}$.

2:while not ($\left\lVert x\right\rVert_{0}>L$ or $\left\lVert
r\right\rVert_{2}\leq\epsilon$) do

3: for all $j\in[m]$ do

4: Compute $\braket{d_{j},r}$

5: end for

6: Select $j^{*}=argmax(\left|\braket{d_{j},r}\right|)$

7: Assign $z=\braket{d_{j^{*}},r}$

8: Update the solution $x_{j}=x_{j}+z$

9: Update the residual $r=r-zd_{j^{*}}$

10:end while

11:Output $x$.

###  III.4 Computational complexity

The analysis of the run-time of the algorithm proceeds as follows. The
initialization step is linear in the length of the residual $O(n)$. The
computation of the sweep stage (steps 3-5) is the bottleneck of this
algorithm. Indeed, the algorithm computes $m$ inner products of vectors of
length $n$, which needs time $O(nm)$. Selecting the best atom has a negligible
cost, as it can be done during the computation of the inner products without
significant overhead. Updating the solution is $O(1)$ and the residual’s
update is bounded by $O(n)$.

The complexity of the sweep stage dominates all the other complexities in the
loop. Therefore, assuming that the matching pursuit converges after $k$
iteration, its asymptotic computational complexity scales as

$\displaystyle O(knm).$ (11)

##  IV Quantum Computing background

###  IV.1 Quantum computation

Just like a bit is the fundamental unit of information in classical computing,
a qubit is the fundamental information unit in quantum computing. A qubit is a
mathematical representation of a quantum mechanical object and can be
described as an $\ell_{2}$ normalized vector of $\mathbb{C}^{2}$. The state of
a $n$-qubit system (a register of a quantum computer) is the tensor product of
single qubits: a unitary vector $\ket{x}\in H^{\otimes
n}\simeq\mathbb{C}^{2^{n}}$. In other words, with $\ket{i}\in H^{\otimes n}$
we denote a quantum register that contains the binary expansion of number $i$.
Its corresponding complex vector is a vector of length $2^{n}$, full of
zeroes, with the $i^{th}$ element equal to one. For instance,

$\displaystyle\ket{3}\in H^{\otimes 2}=\ket{1}\ket{1}=\begin{bmatrix}0\\\ 0\\\
0\\\ 1\\\ \end{bmatrix}.$ (12)

Given a basis $\\{\ket{i}\\}_{0}^{n-1}$ for $H^{\otimes\log_{2}(n)}$, with
$\log_{2}(n)$ qubits, we can describe a quantum state
$\ket{\psi}=\sum_{i}^{n}\alpha_{i}\ket{i}$ with
$\sum_{i}^{n}|\alpha_{i}|^{2}=1$. The values $\alpha_{i}\in\mathbb{C}$ are
called amplitudes of the quantum states $\ket{i}$ for the state $\ket{\psi}$.

The evolution of a quantum system is described by unitary matrices $U$.
Unitary matrices are norm-preserving and thus can be used as a suitable
mathematical description of pure quantum evolutions. Any quantum algorithm
that does not perform measurements can be represented by a unitary matrix.

Quantum states can be measured, but measurements alter the state. In this
work, measurements are performed with respect to the computational basis
$\\{\ket{i}\\}_{0}^{n-1}$ of $H^{\otimes\log_{2}(n)}$. This means that if we
measure a quantum register $\ket{\psi}=\sum_{i}^{n}\alpha_{i}\ket{i}$, it can
collapse to any state $\ket{i}$, each with probability $|\alpha_{i}|^{2}$. It
is important to recall that no quantum algorithm can create a copy of a
generic quantum state. Therefore, to measure a state multiple times, it is
necessary to create it again from scratch every time.

For a deeper introduction to the subject, we encourage the reader to consult
Nielsen and Chuang [16].

###  IV.2 Quantum random access memory

A QRAM, or quantum ram, is a device that: given a list of $n$ bit strings
$x_{i}\in\\{0,1\\}^{m}$ of length $m$, performs the following mapping in time
$O(\text{polylog}(n))$:

$\displaystyle\sum_{i}^{n}\alpha_{i}\ket{i}\ket{0}\mapsto\sum_{i}^{n}\alpha_{i}\ket{i}\ket{x_{i}},$
(13)

where $\alpha_{i}\geq 0$ and $\ket{x_{i}}\in H^{\otimes m}$ is the state of
the computational basis that corresponds to the bit string $x_{i}$.

The reader can think of it as a quantum equivalent to the classical ram. In a
classical ram, we can store $n$ values and query any of those in time $O(1)$,
considering it can access the $m$ bits in parallel. The main difference with a
classical ram is that a quantum ram needs to perform queries in superposition.

As explained in the next section, our quantum algorithm, like many previous
ones, assumes the availability of such a device to encode scalars, matrices,
and vectors in quantum states efficiently.

Building a fault-tolerant, hardware-efficient QRAM is not an easy task. One of
the most promising proposals to structure the quantum random access memory is
the bucket-brigade architecture. First presented in Giovannetti  _et al._
[17], this architecture is composed of $O(n)$ gates, while its circuit is only
$O(\log(n))$ deep. Current error analysis research claim that algorithms that
query the bucket-brigade QRAM a super-polylogarithmic number of times (e.g.,
Grover’s search, and consequently ours) likely require the bucket-brigade QRAM
to be error corrected [18, 19]. This requires additional hardware and created
skepticism on the effective speed-up of this class of algorithms in this input
model [18]. Recently, Hann  _et al._ [19] have shown that the bucket-brigade
architecture is highly resilient to generic errors and that its architectural
advantage persists even in case of error correction, contrary to what was
previously believed. At the same time, a recent work has shown how to build a
fault-tolerant bucket-brigade QRAM by parallelizing Clifford + T gates [20],
at the cost of using $O(n)$ ancillary qubits.

In this work, we will perform our analysis assuming access to a fault-tolerant
QRAM, capable of performing queries in time
$O(\text{polylog}(n))\sim\widetilde{O}(1)$, with $n$ the number of entries
stored in the QRAM. With this assumption in mind, we can explain our data
access model.

###  IV.3 Quantum data access

We can encode a scalar $a\in\mathbb{R}$ in a quantum register $\ket{a}$ using
its binary encoding and retrieve it by measuring the register in the
computational basis, just as discussed in the previous section. To encode a
scalar, we need as many qubits as the classical bits that store it. For more
detailed information on how to encode a real number in a quantum register, we
suggest reading Nielsen and Chuang [16, Chapter 5, Sections 5.1, 5.2] to see
how phase angles are encoded in a state, and Rebentrost  _et al._ [21] for a
more formal definition of a quantum arithmetic model with fixed point
precision.

On the other hand, the components of a vector $a\in\mathbb{R}^{m}$ can be
encoded with fewer qubits than classical bits, using the amplitudes of a
quantum state. When measured, the quantum state collapses to an index of the
vector with probability proportional to the magnitude of the indexed
component. We call this representation state-vector.

######  Definition 2 (State-vector).

Given a vector $x\in\mathbb{R}^{n}$, the corresponding state-vector is the
following quantum state $\ket{x}=\frac{1}{\left\lVert
x\right\rVert_{2}}\sum_{i}^{n}x_{i}\ket{i}$, which is encoded in $\lceil\log
n\rceil$ qubits.

We say to have quantum access to a classical vector $x\in\mathbb{R}^{n}$ if we
have access to a unitary operator that performs the mapping
$U_{x}:\ket{0}\mapsto\ket{x}$ in time $O(log(n))$. Similarly, we define the
concept of quantum access to a matrix.

######  Definition 3 (Quantum access to a matrix).

We say to have quantum access to a matrix $A\in\mathbb{R}^{n\times m}$ if we
can perform the following mappings in time $O(\text{polylog}(nm))$:
•
$U:\ket{i}\ket{0}\mapsto\ket{i}\ket{a_{i,\cdot}}=\ket{i}\frac{1}{||a_{i,\cdot}||}\sum_{j}^{m}a_{ij}\ket{j}$,
for $i\in\mathbb{R}^{n}$;  •
$V:\ket{0}\mapsto\frac{1}{\|A\|_{F}}\sum_{i}^{n}\|a_{i,\cdot}\|\ket{i}$.

By combining the two unitaries above, it is possible to represent a matrix in
a quantum state

$\displaystyle\ket{A}=U(V\otimes\mathbb{I})\ket{0}\ket{0}$
$\displaystyle=\frac{1}{\left\lVert
A\right\rVert_{F}}\sum_{i}^{n}\sum_{j}^{m}a_{ij}\ket{i}\ket{j}$ (14)
$\displaystyle=\frac{1}{\left\lVert A\right\rVert_{F}}\sum_{i}^{n}\left\lVert
a_{i,\cdot}\right\rVert\ket{i}\ket{a_{i,\cdot}}$ (15)

using two registers of $\lceil\log(n)\rceil+\lceil\log(m)\rceil$ qubits.

The appendix of Kerenidis and Prakash [22] shows in detail how to construct a
classical data structure that enables efficient computation of the unitary
matrices that give quantum access to vectors and matrices. This classical data
structure can be created in $O(nm\log^{2}(nm))$ for a $n\times m$ matrix and
in $O(n\log(n))$ for a vector of length $n$. Assuming the ability to perform
quantum queries to the entries of this data structure in superposition (i.e.,
assuming the availability of a QRAM that stores the entries of these trees),
the authors show how to provide quantum access to a vector/matrix in time
$\widetilde{O}(1)$.

In practice, in this input model, at the cost of some classical pre-
processing, it is possible to encode vectors and matrices in quantum states
using a small number of qubits in time $\widetilde{O}(1)$.

###  IV.4 Relevant quantum subroutines

We now introduce two quantum subroutines that are particularly important to
our work. By slightly modifying these two results, we introduce two new
corollaries that are more suitable to our needs.

######  Lemma 4 (Inner product estimation [23]).

Let there be quantum access to the matrices $V\in\mathbb{R}^{n\times m}$ and
$C\in\mathbb{R}^{k\times m}$, through the unitaries
$U_{v}:\ket{i}\ket{0}\mapsto\ket{i}\ket{v_{i,\cdot}}$ and
$U_{c}:\ket{j}\ket{0}\mapsto\ket{j}\ket{c_{j,\cdot}}$, that run in time $T$.
Then, for any $\delta>0$ and $\epsilon>0$, there exists a quantum algorithm
that computes
$\ket{i}\ket{j}\ket{0}\mapsto\ket{i}\ket{j}\ket{\overline{\braket{v_{i,\cdot},c_{j,\cdot}}}}$,
such that
$|\overline{\braket{v_{i,\cdot},c_{j,\cdot}}}-\braket{v_{i,\cdot},c_{j,\cdot}}|\leq\epsilon$,
with probability greater than $1-2\delta$ in time
$\widetilde{O}\left(\frac{T\log(1/\delta)}{\epsilon}\right)$.

In our quantum matching pursuit we will need to perform inner products between
the columns of a matrix and a vector. It is possible to use this result to
prepare a state that stores the inner products between the columns of a matrix
$A$ and a column vector $x$. Indeed, we only need quantum access to the
matrix’s transpose $A^{T}$ and to the vector $x$ via unitaries $U_{A}$ and
$U_{x}$. If we have that, the two unitaries that prepare the states trivially
become $U_{v}=U_{A^{T}}$ and $U_{c}=(\mathbb{I}\otimes U_{x})$, and we can
ignore the second register. We also stress that having quantum access to a
matrix is equivalent to having access to its transpose. Indeed, if we swap the
first and the second register of Equation 14, we have a quantum state that
represents the transposed matrix. Finally, using consistent phase estimation
from Ta-Shma [24], it is possible to modify the algorithm above so that the
error across several runs is consistent.

######  Corollary 5 (Matrix-vector product estimation).

Let there be quantum access to a matrix $A\in\mathbb{R}^{n\times m}$ and to a
unit vector $x\in\mathbb{R}^{m}$, through the unitaries
$U_{A}:\ket{0}\ket{0}\mapsto\frac{1}{\left\lVert
A\right\rVert_{F}}\sum_{i}^{n}\ket{i}\ket{a_{i,\cdot}}$ and
$U_{x}:\ket{0}\mapsto\frac{1}{\left\lVert
x\right\rVert_{2}}\sum_{i}^{m}x_{i}\ket{i}$, that run in time less than $T$.
Then, for any $\delta>0$ and $\epsilon>0$, there exists a quantum algorithm
that computes $\ket{0}\ket{0}\mapsto\frac{1}{\left\lVert
A\right\rVert_{F}}\sum_{i}^{n}\ket{i}\ket{\overline{\braket{a_{i,\cdot},x}}}$,
such that
$|\overline{\braket{a_{i,\cdot},x}}-\braket{a_{i,\cdot},x}|\leq\epsilon$
consistently across multiple runs, with probability greater than $1-2\delta$
in time $\widetilde{O}\left(\frac{T\log(m/\delta)}{\epsilon}\right)$.

To prove the bound on the success probability of the corollary above, we can
exploit the union bound, also known as Boole’s inequality

$\displaystyle P(\cup_{i}^{m}p_{f}(i))\leq\sum_{i}^{m}p_{f}(i).$ (16)

The union bound states that, given a set of likely events, the probability
that any one of them happens is lower than the sum of their individual
probabilities.

The run-time overhead of Lemma 4 to bound the failure probability of one inner
product with $p_{f}(i)\leq 2\delta$ is $O(\log(1/\delta))$. Bounding the
failure probability of one product by $p_{f}(i)\leq 2\delta^{\prime}$, the
probability $P_{f}$ that one of the $m$ inner products fail can be bounded by

$\displaystyle P_{f}\leq\sum_{i}^{m}p_{f}(i)\leq 2m\delta^{\prime}.$ (17)

Choosing $\delta^{\prime}=\frac{\delta}{m}$, we have that the algorithm
succeeds with probability $1-2\delta$ with a run-time overhead of
$O\left(\log\left(\frac{m}{\delta}\right)\right)$.

Finally, if $x$ is not a unit vector, we can multiply the result by
$\left\lVert x\right\rVert_{2}$ to get an estimate within $\epsilon\left\lVert
x\right\rVert_{2}$ error.

The second algorithm that we introduce enables searching for the minimum value
of an unsorted array quadratically faster than what we can do classically.

######  Lemma 6 (Finding the minimum [25]).

Let there be quantum access to a vector $u\in[0,1]^{N}$ via the operation
$\ket{j}\ket{0}\to\ket{j}\ket{u_{j}}$ in time $T$. Then, we can find the
minimum $u_{\min}=\min_{j\in[N]}u_{j}$ and its index
$j_{\min}=\operatorname*{arg\,min}_{j\in[N]}u_{j}$ with success probability
$1-\delta$ in time $O\left(T\sqrt{N}\log\left(\frac{1}{\delta}\right)\right)$.

This algorithm is built around the famous Grover’s search algorithm [26].
Grover’s algorithm takes advantage of an oracle to mark the elements of the
superposition that satisfy the search conditions. An oracle is a function
$f:\mathbb{R}\mapsto\\{0,1\\}$. The Finding the minimum routine uses Grover’s
search many times, using oracles of the type

$\displaystyle f_{i}(j)=\begin{cases}1&\mbox{if }u_{j}<u_{i}\\\
0&\mbox{otherwise}\end{cases}.$ (18)

By using the same algorithm described in Durr and Hoyer [25] with oracles

$\displaystyle f^{\prime}_{i}(j)=\begin{cases}1&\mbox{if }|u_{j}|>|u_{i}|\\\
0&\mbox{otherwise}\end{cases},$ (19)

it is possible to find the maximum absolute value of an array in the same
running time.

######  Corollary 7 (Finding the maximum absolute value).

Let there be quantum access to a vector $u\in[0,1]^{N}$ via the operation
$\ket{j}\ket{0}\to\ket{j}\ket{u_{j}}$ in time $T$. Then, we can find the
maximum absolute value $u_{\max}=\max_{j\in[N]}|u_{j}|$ and its index
$j_{\min}=\operatorname*{arg\,min}_{j\in[N]}u_{j}$ with success probability
$1-\delta$ in time $O\left(T\sqrt{N}\log\left(\frac{1}{\delta}\right)\right)$.

##  V Quantum matching pursuit

Now that we have introduced the necessary quantum background, we show how to
construct the quantum matching pursuit algorithm. The proposed algorithm
closely follows its classical counterpart, but it takes advantage of the
quantum routines presented in Section IV.4. The main idea is to exploit the
quantum regime to speed up the sweep stage.

###  V.1 Data access

First, we need quantum access to the dictionary $D$. We recall that, without
loss of generality, we can consider the columns of $D$ to be $\ell_{2}$
normalized vectors, meaning that the Frobenius norm of $D$ is $\left\lVert
D\right\rVert_{F}=\sqrt{m}$. Recall that, with a pre-processing time of
$O(nm\log^{2}(nm))$, we can create quantum access to the following quantum
state in time $\widetilde{O}(1)$

$\displaystyle\ket{D}=\frac{1}{\sqrt{m}}\sum_{j}^{m}\ket{j}\ket{d_{j}},$ (20)

where $d_{j}$ is the $j^{th}$ column of $D$ and $\ket{d_{j}}$ is a state-
vector. This pre-processing cost needs to be paid only once and allows
applying the quantum matching pursuit on any signal $s\in\mathbb{R}^{n}$ that
is sparse over $D$. For this reason, we will not include this cost in our run-
time analysis.

Figure 1: The tree structure that enables efficient quantum access to the
vector of the residuals. Each node stores the sum of squares of the leaves
that descend from that node.

Similarly, we can create quantum access to the residual $r\in\mathbb{R}^{n}$
by using a tree data structure. For the sake of clarity, we report the data
structure in Figure 1. This data structure can be created in time
$O(n\log(n))$ for a vector with $n$ components and enables quantum access to
the following quantum state in time $\widetilde{O}(1)$ [22, 27]

$\displaystyle\ket{r}=\frac{1}{\left\lVert
r\right\rVert_{2}}\sum_{i}^{n}r_{i}\ket{i}.$ (21)

Preparing access to $\ket{r}$ and $\ket{D}$ means implementing the circuits
described by the unitaries discussed in Section IV.3. Finally, both the signal
$s\in\mathbb{R}^{n}$ and its sparse representation $x\in\mathbb{R}^{m}$ are
represented as classical arrays, as they will not be encoded in quantum
states. To further lower the memory complexity, it is possible to exploit the
sparseness of $x$ and store it in a data structure which uses $O(\left\lVert
x\right\rVert_{0})$ space (e.g., a hash table). This will not affect the
overall time complexity as long as the data structure has a constant time
insertion/update cost.

###  V.2 Quantum algorithm

We start by initiating the residual structure with the signal components.
Since we have quantum access to the residual and to the dictionary as in
Equations 21 and 20, we can use the matrix-vector product estimation procedure
from Corollary 5 to produce the state

$\displaystyle\ket{\varphi}=\frac{1}{\sqrt{m}}\sum_{j}^{m}\ket{j}\ket{\overline{z}_{j}}.$
(22)

We recall that $z_{j}=\braket{d_{j},r}$ and that the procedure computes
$\overline{z}_{j}$ such that

$\displaystyle|\overline{z}_{j}-z_{j}|\leq\xi\left\lVert r\right\rVert_{2}$
(23)

where $\xi\in\mathbb{R}^{+}$ is a parameter of arbitrary choice and
$\left\lVert r\right\rVert_{2}$ is the residual’s $\ell_{2}$ norm at the
current iteration.

Once we have a quantum register $\ket{\varphi}$ with a superposition of all
the inner products, we can perform the Finding the maximum absolute value
routine from Corollary 7 to find $j^{*}$ and $z_{j^{*}}$. With the best $j$
and $z_{j}$, we can proceed to update the solution $x\in\mathbb{R}^{m}$ and
the residual $r\in R^{n}$.

Just like in the classical algorithm, we repeat this procedure until the norm
of the residual is lower than a threshold $\epsilon$ or the solution is such
that $\left\lVert x\right\rVert_{0}>L$ for a threshold $L\in\mathbb{N}^{+}$.

Algorithm 2 concisely summarizes the quantum matching pursuit procedure.

Algorithm 2 Quantum matching pursuit

1:Initialize $r=s$, $x=0^{\otimes m}$.

2:while not ($\left\lVert x\right\rVert_{0}>L$ or $\left\lVert
r\right\rVert_{2}\leq\epsilon$) do

3: Prepare $\ket{r}=\frac{1}{\left\lVert
r\right\rVert_{2}}\sum_{i}^{n}r_{i}\ket{i}$

4: Prepare $\ket{D}=\frac{1}{\sqrt{m}}\sum_{j}^{m}\ket{j}\ket{d_{j}}$

5: Use Inner product estimation to create
$\ket{\varphi}=\frac{1}{\sqrt{m}}\sum_{j}^{m}\ket{j}\ket{\overline{z}_{j}}$,
where $z_{j}=\braket{d_{j},r}$ and $|\overline{z}_{j}-z_{j}|\leq\xi\left\lVert
r\right\rVert_{2}$

6: Apply Finding the maximum absolute value to $\ket{\varphi}$ to obtain
$j^{*}$ and $z_{j^{*}}$

7: Update the solution $x_{j^{*}}=x_{j^{*}}+z_{j^{*}}$

8: Update quantum access to the residual $r=r-z_{j^{*}}d_{j^{*}}$

9:end while

10:Output $x$.

It is important to remark that the error in the inner products can introduce
convergence issues and could lead to worse solutions than the ones provided by
the classical algorithm. We have identified two alternative versions of the
quantum matching pursuit algorithm:

  1. 1.

_Single-error_ : The error only affects the computation of $j^{*}$. We re-
compute $z_{j^{*}}$ classically after the finding procedure outputs the atom’s
index;

  2. 2.

_Double-error_ : The error affects both the computation of $j^{*}$ and
$z_{j^{*}}$, which are retrieved quantumly by the finding procedure.

While the asymptotic complexity of the algorithm does not change between the
two versions, the second one is slightly faster to compute. In Section VI.1,
we provide experimental evidence that the two versions of the algorithm do not
significantly affect the performance of quantum matching pursuit in practice,
supporting the use of the second version. We also provide numerical
experiments to check whether the quantum algorithm can find representations of
the same quality as its classical counterpart.

###  V.3 Success probability

Differently from its classical counterpart, the quantum algorithm that we
propose has a probability of failing. In this section, we discuss the success
probability of our algorithm, showing that we can arbitrarily bound its
failure probability with little running time overhead.

The probability of failure is due to the computation of the inner products at
step 5 and the search at step 6. In particular, Corollaries 5 and 7 fail with
a probability smaller than $2\delta_{1}^{\prime}$ and $\delta_{2}^{\prime}$,
with a run-time overhead of $O(\log(\frac{m}{\delta_{1}^{\prime}}))$ and
$O(\log(\frac{1}{\delta_{2}^{\prime}}))$ respectively.

Let us denote the probability of failure of the $i^{th}$ loop iteration by
$p_{f}(i)$ and consider
$\delta_{1}^{\prime}=\delta_{2}^{\prime}=\delta^{\prime}$. The probability of
success of the $i^{th}$ loop iteration is:

$\displaystyle
1-p_{f}(i)\geq(1-2\delta^{\prime})(1-\delta^{\prime})=1-3\delta^{\prime}+2\delta^{\prime
2}\geq 1-3\delta^{\prime}$ (24)

with a run-time overhead of
$O(\log(\frac{m}{\delta^{\prime}})\log(\frac{1}{\delta^{\prime}}))=O(\log(\frac{m+1}{\delta^{\prime}}))$.
So, the probability of failure of the $i^{th}$ loop iteration is bounded by
$3\delta^{\prime}$.

Given that the loop is executed $k$ times, by union bound, we can bound the
total probability of failure as

$\displaystyle p(\cup_{i}^{k}p(i))\leq\sum_{i}^{k}p(i)\leq 3k\delta^{\prime}.$
(25)

It follows that the success probability of our algorithm is
$1-3k\delta^{\prime}$ with a run-time overhead of
$O(\log(\frac{m+1}{\delta^{\prime}}))$.

By choosing $\delta^{\prime}=\frac{\delta}{3k}$, we can affirm that our
algorithm succeeds with probability greater than $1-\delta$ with a run-time
overhead of $O(\log(\frac{3k(m+1)}{\delta}))\sim O(\log(\frac{3km}{\delta}))$.

Using the same proof technique, we could choose two different probabilities of
failure $\delta_{1}^{\prime}$ and $\delta_{2}^{\prime}$ and further reduce the
overhead by a constant factor. For instance, by choosing
$\delta_{1}^{\prime}=\frac{\delta}{2k}$ and
$\delta_{2}^{\prime}=\frac{\delta}{2km}$, we can bound the run-time overhead
with $O(\frac{2km}{\delta})$. Finally, note that even if the exact value of
$k$ is not known beforehand, we can find suitable $\delta_{1}^{\prime}$ and
$\delta_{2}^{\prime}$ by considering $k\sim O(L)$.

###  V.4 Running time

Finally, we provide a thorough analysis of the run-time of the new proposed
algorithm, proving its computational complexity. As already said, we assume
that the dictionary has been stored in an appropriate data structure, and we
analyze the time to compute the sparse representation of a signal over that
dictionary.

The first step of the algorithm consists in initializing the residual. This
can be done by building the data structure in Figure 1, and it requires time
$O(n\log(n))$. Since the residual and the dictionary are stored in adequate
data structures, the preparation of states $\ket{D}$ and $\ket{r}$ at steps 3
and 4 is $\widetilde{O}(1)$, as discussed in Section IV.3.

Given that the cost of preparing the two quantum states is $\widetilde{O}(1)$,
from Corollary 5, we know that the cost of step 5 is
$\widetilde{O}(\frac{1}{\xi})$ and that it encodes in the register values
$\overline{z}_{j}$, with error $|\overline{z}_{j}-z_{j}|\leq\xi\left\lVert
r\right\rVert_{2}$. Therefore, we can prepare the state

$\displaystyle\ket{\varphi}=\frac{1}{\sqrt{m}}\sum_{j}^{m}\ket{j}\ket{\overline{z}_{j}}$
(26)

in time $\widetilde{O}(\frac{1}{\xi})$.

We can ignore the run-time terms that depend on $\delta$ as they are related
to the success probability studied in the previous section. We will include
this overhead in the run-time at the end of this section.

Recall that the state above is the superposition of $m$ inner products, among
which we need to identify the one with the greatest maximum absolute value.
Using Corollary 7, we can find the value $z_{j^{*}}$ and its index $j^{*}$ in
time $O(\frac{\sqrt{m}}{\xi})$.

Once these values are computed, we need to update the solution and the
residual. To update the solution, we need to add $\overline{z}_{j_{*}}$ to one
component of the solution vector. We can do it in time $O(1)$. On the other
hand, updating the residual is more demanding, as we need to modify the tree
data structure. Each element update costs $O(\log(n))$, as we need to update
one leaf and all its parent nodes. Since the update is
$r=r+z_{j^{*}}d_{j^{*}}$, we need to update $\left\lVert
d_{j^{*}}\right\rVert_{0}$ elements of the residual. This step costs
$O(\left\lVert d_{j^{*}}\right\rVert_{0}\log(n))$.

Let us denote by $d_{j^{*}}^{(i)}$ the best atom at iteration $i$. The cost of
the $i^{th}$ iteration of the loop, results in

$\displaystyle\widetilde{O}\left(\|d_{j^{*}}^{(i)}\|_{0}\log(n)+\frac{\sqrt{m}}{\xi}\right).$
(27)

By assuming that we perform $k$ iterations of the loop, considering the
initial cost of initiating the residual, and introducing the term that
accounts for the success probability, the total complexity of our algorithm is

$\displaystyle\widetilde{O}\left(n\log
n+\sum_{i}^{k}\|d_{j*}^{(i)}\|_{0}\log(n)+k\frac{\sqrt{m}}{\xi}\log\left(\frac{3km}{\delta}\right)\right).$
(28)

If the atoms of our dictionary are not sparse, we can write the run-time in a
more compact way by observing that $\forall j$, $\left\lVert
d_{j}\right\rVert_{0}\leq n$

$\displaystyle\widetilde{O}\left(kn\log
n+k\frac{\sqrt{m}}{\xi}\log\left(\frac{3km}{\delta}\right)\right).$ (29)

We can concisely summarise our result in the following theorem.

######  Theorem 8 (Quantum matching pursuit).

Let there be quantum access to a dictionary $D\in\mathbb{R}^{n\times m}$ and
to a vector $s\in\mathbb{R}^{n}$. Let $\xi,\delta\in\mathbb{R}_{>0}$ be
precision parameters. There exists a quantum algorithm that simulates the
matching pursuit algorithm, with error
$\left|z_{j}^{(i)}-\overline{z}_{j}^{(i)}\right|\leq\xi\left\lVert
r^{(i)}\right\rVert_{2}$ on the inner products estimation at each $i^{th}$
iteration, in running time $\widetilde{O}\left(kn\log
n+k\frac{\sqrt{m}}{\xi}\log\left(\frac{3km}{\delta}\right)\right)$, where $k$
is the total number of iterations. The algorithm succeeds with probability
greater than $1-\delta$.

##  VI Numerical experiments

The quantum algorithm that we propose executes the exact steps of its
classical counterpart, but it introduces some random error along the
computation, possibly affecting the quality of the representation. We present
numerical simulations on synthetic data to better study the run-time advantage
and the convergence properties of the quantum matching pursuit. We study
whether the _Single-error_ and _Double-error_ versions of the algorithm,
introduced in Section V.2, actually affect the algorithm’s correctness. Our
numerical simulations are performed on a classical computer by introducing
artificial errors $\xi\in[-0.01,0.01]$ in the computation of the inner
products.

###  VI.1 Single/double-error and representation quality

We run experiments to test whether the two implementations of the quantum
matching pursuit algorithm produce different results and whether the
representation quality is different from the one obtained by the classical
matching pursuit.

To do so, we proceed by creating $100$ batches, each of which composed of:

  * •

a random dictionary $D\in\mathbb{R}^{100\times 512}$, with unit columns

  * •

$100$ sparse vectors $x\in\mathbb{R}^{m}$, with $\|x\|_{0}=17$

  * •

$100$ signals of the form $s=Dx+\epsilon$, where $\epsilon$ is a vector
containing Gaussian truncated noise.

The dictionary, the signals and the sparse vectors are provided by Scikit-
Learn’s [28] _make_sparse_coded_signal_ function. We add the $\epsilon$ noise
artificially.

We simulate the two versions of the quantum matching pursuit algorithm and the
classical one for each batch, to compare the sparse representations on the
same data. In order to assess the convergence properties, we do not set any
threshold $L$ that limits the sparsity of the solution.

We observe that all three versions represent the signals with $18$ components
on average. For both the _Single-error_ and the _Double-error_ version, we
compute, for each batch, the average sparsity of the solution for that batch,
divided by the one obtained using the classical version. The average value for
this metric is respectively $1.0048$ and $1.0060$, indicating that the quantum
and classical algorithms compute solutions with similar sparsity, and also
that the two quantum versions do not differ much.

To better prove the latter point, we run a statistical test. First, we check
whether the metric values obtained for each quantum algorithm version are
normally distributed, using a Shapiro-Wilk normality test.

The resulting $p$-values are respectively $0.37490$ for the _Double-error_
version and $0.07854$ for the _Single-error_ one. Since the _Single-error_
results are not normally distributed, and the experiments have been conducted
on the same batches of data, we run a Wilcoxon signed-rank test to determine
whether the performances are different. The test outcome is a $p$-value of
$0.34790$, which does not allow us to reject the hypothesis that the two
algorithms perform similarly.

###  VI.2 Run-time simulation

While it is true that the classical complexity of the matching pursuit, in the
general case, scales as $O(knm)$, and that the quantum version scales as
$O(kn\log(n)+k\frac{\sqrt{m}}{\xi}\log(\frac{3km}{\delta})\text{polylog}(nm))$,
it is legitimate to wonder whether the number of iterations $k$ is the same
for both the classical and the quantum version.

To analyze the run-time of the algorithms, we compare the classical matching
pursuit with the _Double-error_ quantum matching pursuit on the same dataset.

We generate datasets in the same way discussed in the previous section, with
the difference that we study the run-time as the length $n$ of the signal
increases. For each batch, the number of atoms $m$ is set to twice the length
of the signal, while solution is five times sparser than the original signal.

Experimentally, we see that there is not a significant difference in the
number of iterations required to converge. Indeed, the two algorithms produce
solutions of the same sparsity. Figure 2 illustrates the run-times of the
classical and quantum algorithms, considering a probability of failure
$\delta=0.01$. Each point represents the average number of operations required
to find the sparse representation of the signal, while the red bars show the
standard deviation.

We notice that the quantum algorithm is not advantageous until a certain
number of signal components. While the exact signal length for quantum
advantage depends on the characteristics of the problem (e.g., number of
atoms, expected sparsity), this experiment shows that the quantum matching
algorithm provides a substantial speed-up over its classical counterpart on
high-dimensional signals.

Figure 2:  Running times for classical and quantum matching pursuit. Each
point is the average number of operations required to find the sparse
representation of a signal $s\in\mathbb{R}^{n}$. We report the standard
deviation with red bars.

##  VII Conclusions

In this work, we have presented a novel quantum algorithm to find sparse
representations of signals in the QRAM input model. The new algorithm is a
quantum version of the classical matching pursuit algorithm. We have
leveraged, modified, and combined two quantum subroutines to speed up the
sweep stage of the matching pursuit. The result is a new routine with a
polynomial advantage over its classical counterpart while retrieving solutions
of the same quality.

This work shows that quantum computing can impact learning sparse
representation, and opens the path to further research on other quantum
pursuit algorithms.

###### Acknowledgements.

The authors would like to thank Dr Alessandro Luongo for his valuable advice,
and professors Ferruccio Resta and Donatella Sciuto for their support. They
are particularly grateful to Prof. Giacomo Boracchi for his inspiring lectures
on sparse representations, and to Prof. Timothy J. Sluckin for his precious
feedback on the first draft of this manuscript.

## References

  * Rauhut  _et al._ [2008] H. Rauhut, K. Schnass, and P. Vandergheynst, Compressed sensing and redundant dictionaries, IEEE Transactions on Information Theory 54, 2210 (2008).
  * Pennebaker and Mitchell [1992] W. B. Pennebaker and J. L. Mitchell,  _JPEG: Still image data compression standard_ (Springer Science & Business Media, 1992).
  * Elad and Aharon [2006] M. Elad and M. Aharon, Image denoising via sparse and redundant representations over learned dictionaries, IEEE Transactions on Image processing 15, 3736 (2006).
  * Adler  _et al._ [2015] A. Adler, M. Elad, Y. Hel-Or, and E. Rivlin, Sparse coding with anomaly detection, Journal of Signal Processing Systems 79, 179 (2015).
  * Bellante  _et al._ [2021] A. Bellante, A. Luongo, and S. Zanero, Quantum algorithms for data representation and analysis, arXiv preprint arXiv:2104.08987 (2021).
  * Kerenidis and Luongo [2020] I. Kerenidis and A. Luongo, Classification of the mnist data set with quantum slow feature analysis, Physical Review A 101, 062327 (2020).
  * Kerenidis and Landman [2021] I. Kerenidis and J. Landman, Quantum spectral clustering, Physical Review A 103, 042415 (2021).
  * Mallat and Zhang [1993] S. G. Mallat and Z. Zhang, Matching pursuits with time-frequency dictionaries, IEEE Transactions on signal processing 41, 3397 (1993).
  * Krstulovic and Gribonval [2006] S. Krstulovic and R. Gribonval, Mptk: Matching pursuit made tractable, in  _2006 IEEE International Conference on Acoustics Speech and Signal Processing Proceedings_ , Vol. 3 (IEEE, 2006) pp. III–III.
  * Wu and Batista [2003] Y. Wu and V. S. Batista, Matching-pursuit for simulations of quantum processes, The Journal of Chemical Physics 118, 6720 (2003).
  * Wu and Batista [2004] Y. Wu and V. S. Batista, Quantum tunneling dynamics in multidimensional systems: A matching-pursuit description, The Journal of Chemical Physics 121, 1676 (2004).
  * Chen and Batista [2006] X. Chen and V. S. Batista, Matching-pursuit/split-operator-fourier-transform simulations of excited-state nonadiabatic quantum dynamics in pyrazine, The Journal of Chemical Physics 125, 124313 (2006).
  * Wu  _et al._ [2005] Y. Wu, M. F. Herman, and V. S. Batista, Matching-pursuit/ split-operator fourier-transform simulations of nonadiabatic quantum dynamics, The Journal of Chemical Physics 122, 114114 (2005).
  * Natarajan [1995] B. K. Natarajan, Sparse approximate solutions to linear systems, SIAM journal on computing 24, 227 (1995).
  * Bennett  _et al._ [1997] C. H. Bennett, E. Bernstein, G. Brassard, and U. Vazirani, Strengths and weaknesses of quantum computing, SIAM journal on Computing 26, 1510 (1997).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, 2010).
  * Giovannetti  _et al._ [2008] V. Giovannetti, S. Lloyd, and L. Maccone, Architectures for a quantum random access memory, Physical Review A 78, 052310 (2008).
  * Arunachalam  _et al._ [2015] S. Arunachalam, V. Gheorghiu, T. Jochym-O’Connor, M. Mosca, and P. V. Srinivasan, On the robustness of bucket brigade quantum ram, New Journal of Physics 17, 123010 (2015).
  * Hann  _et al._ [2021] C. T. Hann, G. Lee, S.M. Girvin, and L. Jiang, Resilience of quantum random access memory to generic noise, PRX Quantum 2, 020311 (2021).
  * Paler  _et al._ [2020] A. Paler, O. Oumarou, and R. Basmadjian, Parallelizing the queries in a bucket-brigade quantum random access memory, Physical Review A 102, 032608 (2020).
  * Rebentrost  _et al._ [2021] P. Rebentrost, M. Santha, and S. Yang, Quantum alphatron, arXiv preprint arXiv:2108.11670 (2021).
  * Kerenidis and Prakash [2017] I. Kerenidis and A. Prakash, Quantum recommendation systems, in  _8th Innovations in Theoretical Computer Science Conference (ITCS 2017)_ (Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2017).
  * Kerenidis  _et al._ [2019] I. Kerenidis, J. Landman, A. Luongo, and A. Prakash, q-means: A quantum algorithm for unsupervised machine learning, in  _NeurIPS 2019_ (2019).
  * Ta-Shma [2013] A. Ta-Shma, Inverting well conditioned matrices in quantum logspace, in  _Proceedings of the forty-fifth annual ACM symposium on Theory of computing_ (2013) pp. 881–890.
  * Durr and Hoyer [1996] C. Durr and P. Hoyer, A quantum algorithm for finding the minimum, arXiv preprint quant-ph/9607014 (1996).
  * Grover [1996] L. K. Grover, A fast quantum mechanical algorithm for database search, in  _Proceedings of the twenty-eighth annual ACM symposium on Theory of computing_ (1996) pp. 212–219.
  * Grover and Rudolph [2002] L. Grover and T. Rudolph, Creating superpositions that correspond to efficiently integrable probability distributions, arXiv preprint quant-ph/0208112 (2002).
  * Pedregosa  _et al._ [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, _et al._ , Scikit-learn: Machine learning in python, the Journal of Machine Learning Research 12, 2825 (2011).