\section{\label{sec:intro}Introduction}

Supervised learning is one of the subareas of machine learning \cite{inner_product,schuld2015introduction,biamonte2017quantum} that consists of techniques to learn to classify new data taking as example a training set. More specifically, the computer is given a training set $X$, consisting on $n$ pairs of point and label, $(x,y)$. With the information, the computer is supposed to extract or infer the conditional probability distributions $p(y|x)$ and use it to classify new points $x$. This paradigm is in contrast with unsupervised learning, that like in the case of clustering, attempts to find structure to a set of points without labels; and reinforcement learning \cite{paparo2014quantum}, where an agent has to figure out the best policy or action for each situation it may face.

An important kind of supervised learning is what is usually called \textit{active learning} \cite{melnikov2018active}. To introduce this concept suppose that we have a supervised learning algorithm, with its corresponding training set. However, instead of directly trying to predict the label of new points, we give the classifier the option to pose us interesting questions in order to reduce the uncertainty in $p(y|x)$. In this setting, the algorithm will add to its training set new points in areas where it has a lot of uncertainty. 

To explain the concept better, let us give an example. Suppose we have an image classifier used for a self-driving car, that has to to distinguish between cars, pedestrians... Internally, images are decomposed into pixels that can be characterised by their combination of red, blue and green. Thus, an image can be expressed as an array of 3 dimensional vectors, or as a large vector if we flatten the array. Any image is then a vector which can be identified with a point of a high dimensional space containing all images.
The set of images of different classes forms a single manifold $\mathcal{M}$ embedded in that highly dimensional space. In particular, the dimension of the space is $3n_p$, for $n_p$ the number of pixels. The key idea of an active learning algorithm is one that is able to understand in what kind of images it has the most uncertainty, and request additional examples to be labeled and added to its training set. 

Another important concept in the context of supervised learning is that of \textit{adversarial attacks} or \textit{adversarial examples}, the name given to a phenomenon where a trained and accurate (usually a neural network) classifier, can be mislead into wrong classification, by producing a carefully chosen and slightly modified version of one point that is classified well. Adversarial examples were discovered quite recently \cite{Intriguing_properties,goodfellow2014explaining}, and have received a lot of attention leading to models robust to particular attacks \cite{Madry2018,sinha2017certifying,wong2017provable,raghunathan2018certified}. 

After the discovery of such adversarial examples, considerable effort has been put in explaining why they happen and how they can be avoided. One of the given theoretical reasons links their existence to a high codimension, the difference in dimension in the structure of the classes we are trying to separate with respect to the highly dimensional space in which they are embedded \cite{Geometry_of_adversarial_examples}. Intuitively this means that the dimension of the space of images of cars, for example, is much lower than the entire space of all possible images, of dimension $3n_p$ as said earlier; since we need to impose a lot of constraints for an image to be indeed the image of a car, even if such constraints are not easily definable. 

This suggest a strategy, proposed in theorem 5 of [12], that states that if we are able to cover our classes with a sufficiently fine sampling of the classes, our algorithm will be provably robust against these adversarial examples. However, how fine this sampling is depends crucially on the minimum distance $r_p$ between classes, and since we do not fully know the classes, we also do not know this minimum distance with precision. Overestimating $r_p$ will result in not being able to use theorem 5 of [12]; whereas underestimating it will mean oversampling the classes, with the associated cost.

In this article we present a quantum active learning algorithm that allows for fast sampling of the most informative points that could be added to the training set $X_\mathcal{L}$ to find out this minimum distance $r_p$. Here, the concept of `informativeness' will refer to an expected gain of information, in the sense of improving the estimate of $r_p$. Our aim is to sample points to be added to the training set with the highest possible informativeness. It will be defined as the product of the probability that a given point is in a class, and the inverse of the margin of the Quantum Support Vector Machine, which measures the amount of information gained if the point were really in that class.

The quantum algorithm will allow us to perform this sampling very efficiently, in polylogarithmic time in the dimension of the space $m$ and the number of already classified points, $n$, and polynomial cost in all other variables. To achieve such complexity we will need to use qRAMs \cite{QRAM_bucket_brigade} and techniques from \cite{QSVM} and \cite{QLSAchilds}. By comparison, linear algebra operations will require polynomial cost in $n$ and $m$, when performed using classical computing.

The origin of this advantage is that in order to calculate such `informativeness' we need to solve the Quantum Support Vector Machine, but we will not need to read out the solution. Rather, we will operate quantumly with the output to calculate the scalar value of `informativeness', which we define formally in the next section. Notice that most quantum linear algebra algorithms are better suited for the cases where one does not have to read out the solution entirely, but rather calculate some expected value. This is partially the intuition that motivated our research. 

We also hope that the framing of this article will highlight to the quantum machine learning community, often focused mostly on obtaining quantum advantages, the necessity of designing systems that are not only efficient and capable, but also robust and reliable.

The structure of the remaining sections of the article is the following. In the following subsection \ref{sec:related} we introduce related work on the topics of adversarial examples, active learning, and quantum machine learning; we also briefly mention the differences between adversarial examples and generative adversarial networks. In section \ref{sec:active_learning} we introduce a bit more of background on theorem 5 of \cite{Geometry_of_adversarial_examples}, and explain how to model the problem of finding the minimum distance between classes as an active learning problem. In particular we introduce the important concept of `informativeness'. Finally, some background on Support Vector Machines is reviewed in \ref{sec:background_SVM}. 

 Section \ref{sec:Main algorithm} constitutes the main corpus of the text. To perform the active learning algorithm, in \ref{sec:Probability} we explain how to obtain $P_c(\vec{x})$, the probability that an arbitrary point of the space is in a given class, which is one of the main components of the `informativeness' of such point. The other main component is described in the next two subsections. In section \ref{sec:QSVM} we review the main results that we will be using from our reference \cite{QSVM}, and in \ref{sec:norm} we use the result of the previous section to calculate $|\vec{w}|_{n+1}$, the second main component of the `informativeness' of $\vec{x}$. Subsection \ref{sec:Target} explains the strategy to select a point that with high probability improves the current estimate of the SVM. Section \ref{sec:results} reviews the main results, and section \ref{sec:Complexity} is dedicated to the calculation of the complexity of our algorithm. Finally, in section \ref{sec:Conclusion} we explain our conclusions. In the appendices we include definitions, theorems from other articles that we use, and some technical results.


\subsection{\label{sec:related}Related work.}

Adversarial examples are a danger for any classifier that needs to be robust to perturbation. For instance, adversarial examples can be dangerous when an autonomous car has to recognise traffic signals.
Since adversarial examples were discovered \cite{Intriguing_properties}, there has been lots of work to explain why they happen \cite{goodfellow2014explaining} and also to obtain provably robust models \cite{Madry2018}. In particular some of the most promising ideas to avoid them are related to adversarial training: training against those adversarial examples before the actual adversary has time to pose them to the classifier. This is for instance the model explained in \cite{Madry2018}, where they use projected gradient descent to minimize the maximum expected loss
\begin{equation}
    \min_\theta \left(\mathbb{E}_{(x,y) \sim \mathcal{D}} \left[\max_{\epsilon \in \mathcal{S}} L(\theta,x+\epsilon,y)\right]\right),
    \label{madry minimization}
\end{equation}
where $\mathcal{D}$ is the initial population of points $x$, with true label $y$, and perturbations $\epsilon$ can be taken from a small set $\mathcal{S}$. $L$ is the loss function: the function that measures the difference between the predicted and actual classification, with adjustable parameters $\theta$, and $\mathbb{E}$ indicates expected value. The maximisation represents the work of the adversary, whereas the minimization represents the work to make the classifier robust. This setup certainly works, as long as one can minimize \eqref{madry minimization}, and efforts have been put forward to perform this adversarial training more efficiently \cite{qin2019adversarial}. However, as pointed out in \cite{Geometry_of_adversarial_examples}, in order to work perfectly it would require an exponential number of adversarial examples in the dimension of the problem added to the training set. Thus, additional strategies are worth exploring.  

It is also worth noticing that adversarial attacks are gaining attention in the quantum community. Recently, it has been indicated that this phenomenon is also present in the case of quantum classifiers \cite{quantum_adversarial}, where the dimension plays a very important role: the higher the dimension the easier to carry out those adversarial attacks. Some experimental work on this line is done in \cite{lu2019quantum}.

Our work is additionally strongly related to several forms of \textit{active learning} algorithms. Active learning algorithms are those where the classifier can ask for new points to be classified and added to its training data base. This field can be divided in two main branches\cite{wang2015active}: \textit{query synthesis}, where new examples are created, and \textit{sampling}. The later is subdivided in \textit{stream-based sampling} and \textit{pool-based}. In stream-based sampling one selects one item at a time and decides whether it is worth the cost of classifying it \cite{atlas1990stream}. In pool-based sampling examples are sampled from a large pool of unlabelled data \cite{wang2015active}. As we will see, our algorithm can be used both as a pool-based sampling or as query synthesis, depending on the points used.
The most typical strategies to select the samples are \textit{uncertainty sampling}, that selects points with maximum uncertainty about which class they belong to \cite{lewis1994uncertainty}, or \textit{Query-by-committee}, where the space of classifiers that agree with the data is halved sequentially \cite{melville2004commitee}. Finally there is also the strategy of using \textit{expected error reduction} \cite{roy2001error} that selects those points that on average, weighted according to probability, make the loss function as low as possible. This last procedure is similar to what we are using, except that instead of minimizing the loss function, we select points that on expectation would achieve the margin of the SVM to be as low as possible.

Finally, let us make a brief introduction to quantum algorithms used for quantum machine learning, since the quantum SVM that we have used through the text is just one of the earlier works on this field. The possibility of applying quantum techniques in the machine learning and artificial intelligence domain have attracted in the last decade much interest, and a review of this efforts is \cite{dunjko2018machine}. For example, following work on quantum SVM, various kernel methods have been developed to improve Support Vector Machines, such as \cite{schuld2019quantum} or \cite{perez2020data}. Most of this efforts have been focused on approaches that can be implemented NISQ devices. For instance, the later reference, \cite{perez2020data}, was specially focused on constructing a universal quantum classifier with minimal amount of quantum resources, using a technique called data re-uploading.

Other machine learning techniques that have been quantized include Bayesian learning \cite{zhao2019bayesian}, where a polynomial speedup is expected thanks to the use of linear algebra techniques such as variations over HHL.

Finally, although considerable effort has been put on quantizing machine learning algorithms such as Support Vector Machines and kernel methods, given the success of the Deep Learning (aka neural networks) techniques in tackling many problems that include image recognition and text translation, significant effort has be put in quantizing them. One example of this is \cite{farhi2018classification}, or the recent publication of TensorFlow Quantum \cite{tensorflow_quantum}, the TensorFlow library to work with quantum neural networks. One particular example of a Deep Learning technique that has received much attention both in the classical and quantum machine learning communities are the generative adversarial networks \cite{zoufal2019quantum}. This networks are composed of two parts, a generator and a discriminator that compete, and a distribution $\mathcal{D}$. The generator tries to generate a distribution resembling as much as possible $\mathcal{D}$, and the discriminator tries to differentiate between the actual $\mathcal{D}$, and the one generated by the generator. The game ends when the generator has fully learned $\mathcal{D}$, and precisely this is the reason why adversarial generative networks have been devised as a method to prepare quantum states efficiently.

However, we would like to emphasise that although somewhat related to adversarial examples via adversarial training, the setup in the later case is different. For example, in the setting of adversarial examples our classifier does not need to know how to generate any state of the distribution, but only how to classify an input in one of the classes. In such respect, there is no generative device in adversarial examples, except perhaps the adversary who is trying to fool the classifier, whereas in the generative adversarial training setup there is a single distribution, and therefore a single class. Perhaps more importantly, adversarial examples seem to be a general and somewhat unfortunate feature of many different kinds of classifiers, whereas generative adversarial networks are a particular, although useful, technique.


\section{\label{sec:active_learning}Modelling the active learning problem}


As mentioned in the introduction, the problem we want to solve in this article is the following: 
\begin{Problem}
\textbf{Informative Points Sampling Problem} \label{Informative sampling problem}
Suppose we are given a set of points $S=(\vec{x}_i, \vec{c}_i)$, $i \in \{1,..,n\}$, where $\vec{x}_i \in \mathbb{R}^m$ are the $n$ data points already classified, and $\vec{c}_i$ are label vectors that indicate probability of membership to the possible classes where the point might be classified. Suppose also that these memberships sum up to, at most, 1 for each point, and for simplicity assume that there are just two classes. The problem is: what are the most informative points to be added to $S$ in order to learn the Support Vector Machine, and in particular the minimum distance between classes $r_p$, with better precision?
\end{Problem}

This problem makes reference to two key concepts: Support Vector Machine and `informativeness', that we introduce now. Although the technical definition of Support Vector Machine (SVM) is stated in the section \ref{sec:background_SVM}, it is basically the hyperplane of separation between two classes that maximizes the margin.

In our particular problem, we are interested in finding the minimum distance between two classes, in order to provide the required covering of the two classes that avoid having adversarial attacks. However, it might be the case that both classes have some points in common, e.g. $(\vec{c}_i)_j, (\vec{c}_i)_k> 0$, where $i$ indicates the point and $j$ and $k$ two classes. In order to avoid this, we establish that a point $\vec{x}_i$ is in class $j$ if $(\vec{c}_i)_{j} > 0.8$, and as $\sum_j (\vec{c}_i)_j = 1 $ for any point $\vec{x}_i$, this will avoid any overlap between classes. This value is arbitrary but should be strictly greater than $0.5$. This way we clearly separate the two sets, and the minimum distance between classes is greater than $0$. With this condition one wants to make classes clearly separated, since otherwise it does not make sense the concept of adversarial example, which will be central to our discussion. If there are more than two classes, then one can generalise the previous setting making comparison on all possible pairs of classes, which scales quadratically with the number of classes.

Now let us look at figure \ref{fig:SVM}. We can see that the initially calculated SVM has a greater margin than that of the SVM we would have if we were able to perfectly know both classes. If we used such margin in theorem \ref{Teorema 5}, we would not find a cover that avoids adversarial examples. Thus, we are interested in minimizing the maximum margin, which is given by the true SVM if we knew the actual classes and not just some samples. 

\begin{figure}
    \centering
    \includegraphics[scale=.75]{SVM.pdf}
    \caption{Example of a SVM, characterised by the equation $\vec{w}\cdot \vec{x}-b=0$. The points depicted are the initial training set, with colors according to their classification, and those in black are those not belonging to any class according to our criteria that for a point $\vec{x}_i$ to be in class $j$, $c_{ij}>0.8$. Initially we have the SVM in black, but we want to obtain the green one, that is more accurate. The dotted lines are the ones that we cannot see initially and must find out. The equation of the SVM is the one depicted and the margin $1/|\vec{w}|$ is chosen to be equal to $r_p$. This is achieved when we make the two parallel hyperplanes that indicate the margin to fulfill equations $\vec{w}\cdot \vec{x}-b=\pm 1$. Since making the SVM more precise implies making the margin smaller, we want to find an SVM that maximizes $|\vec{w}|$.}
    \label{fig:SVM}
\end{figure}

Our problem is however in contrast with usual pool-based sampling (see section \ref{sec:related} for more information). In their problem the algorithm usually seeks to classify new points near the SVM. This is not the case for our problem, since we want points that have a membership to a class greater than $0.8$. Points that have no clear membership to a single class are not so interesting and we do not include them in our data set. 

Now we turn to the concept of `informativeness'. How then do we measure how interesting could be to classify an arbitrary point? A good heuristics for our problem is that we are trying to find points that with high probability would decrease the margin $1/|\vec{w}|$ a lot, to get a better estimate of $r_p$. This is because there are two competing conditions on $r_p$. On the one hand, in order to fulfill the condition of the problem, we cannot take $r_p$ larger than it really is, as that would make us choose a cover of the classes that does not fulfill the conditions of the theorem 5 of \cite{Geometry_of_adversarial_examples}. On the other, the smaller $r_p$ is, the more expensive it is to establish the cover of the classes.

The previous paragraph suggests measuring the `informativeness' of a point by dividing the problem in two parts: we first find the probability that a given point $\vec{x}_{n+1}$ is in class $c$, and then multiply this probability by the inverse of the margin distance of the updated SVM, taking into account this would-be newly classified point, $|\vec{w}_{n+1}|$.

The `informativeness' will thus measure an expected value of information gain if a given point was in a given class. As such, it will be the product of a probability of $\vec {x}$ being in class $c$, $P_c(\vec{x})$, times the information we gain, that we will measure as $|\vec{w}|$. The reason for this choice is that we want to add points to training set, that give a more precise account of the classes. That is, we wish to minimize the margin of size $1/|\vec{w}|$, and therefore we want to maximize $|\vec{w}|$.

\begin{Definition}
\textit{`Informativeness'} measures the expected value of information we get by adding a point $\vec{x}_i$ to the training set. We will therefore define it as the product
\begin{equation}
    P_c(\vec{x}_i) \cdot |\vec{w}_{\vec{x}_i}|,
\end{equation}
where $P_c(\vec{x}_i)$ is the probability that $\vec{x}_i$ is in class $c$, and $1/|\vec{w}_{\vec{x}_i}|$ is the size of the margin of the resulting SVM if $\vec{x}_i$ really were in such class.
\end{Definition}

\section{\label{sec:background_SVM}Background on Support Vector Machines}

Since we will rely on them somewhat heavily, in order to carry out those results it is useful to remember some results from \cite{QSVM}, that explains how a quantum Support Vector Machine algorithm works. Let us first introduce the definition of a Support Vector Machine.

\begin{Definition} \textbf{Support Vector Machine (SVM)}
Let $\mathcal{M}$ be a $m$-dimensional manifold with two submanifolds or classes and a set of points already classified $\{\vec{x}_i,y_i\}$, $\vec{x}_i$ being the point and $y_i$ the label. Then a \textit{Support Vector Machine} is a hyperplane in the manifold separating both classes such that minimum distance between $\{\vec{x}_i\}$ and the hyperplane is as large as possible. If the SVM is linear it may be described by equation
\begin{equation}
    \vec{w}\cdot \vec{x} - b = 0,
\end{equation}
and the size of the margin is $1/|\vec{w}|$, which would be equal to $r_p/2$ if the SVM was perfect. One may also define a non-linear SVM using a non-linear kernel for the dot product.

We highlight that the name Support Vector Machine (SVM) will refer to both the algorithm and the separation hyperplane, the decision boundary.
\end{Definition}

The authors of \cite{QSVM} assume that each point $\vec{x}_i$ is labeled with a single class $y_i$, and there are only two classes $y_i = \pm 1$. Given the pairs of data and label $(\vec{x}_i,y_i)_{i\in \{1,...,n\}}$, the authors state that calculating the SVM is equivalent, in the dual formulation, to maximizing over the multipliers $\vec{\alpha}$ of the Lagrangian
\begin{equation}
    L(\vec{\alpha})=\sum_{i=1}^n \alpha_i y_i -\frac{1}{2}\sum_{i,k=1}^n \alpha_i K_{ik} \alpha_k,\label{lagrangian}
\end{equation}
with constraints $\sum_i \alpha_i = 0$ and $\alpha_i y_i \geq 0$ $\forall i$. After this, the result of the classifier is given, for a new point $\vec{x}$ as
\begin{equation}
    y(\vec{x})=\text{sgn}\left(\sum_{i=1}^n\alpha_i \sum_{l,k=1}^m K_{l,k}x_{i,l} x_k+b\right),
\end{equation}
with $K_{l,k}=\vec{x}_l \cdot \vec{x}_k$, and may be rewritten as 
\begin{equation}
    y(\vec{x})= \text{sgn} \left(K(\vec{w},\vec{x})+b\right),\quad \vec{w}:= \sum_i \alpha_i \vec{x}_i. \label{w definition}
\end{equation}
Notice that the matrix $K_{l,k}$ is a kernel matrix that defines a dot product. Since the margins are of length at least 1, this means that for the training data
\begin{equation}
    y_i(\vec{w}\cdot \vec{x}_i+b)\geq 1. \label{hard margin equation}
\end{equation}
 Since there are only two classes, $y_i=\pm 1$, and so $y_i^2=1$. We will later say, with a bit of abuse of notation, that a given point is in class $c$ to mean any of the two possible values of $y_i$. Additionally, we can transform the previous inequality results into equality adding slack variables $e_i$, 
\begin{equation}
    (\vec{w}\cdot \vec{x}_i+b)=y_i-y_i e_i,
\end{equation}
such that if we allow $e_i>0$, we will be allowing for a soft margin, that is, some points may not fulfill \eqref{hard margin equation}. In such case we will add to Lagrangian \eqref{lagrangian} a term $(\gamma/2) \sum_j e_j^2$, where $\gamma$ is specified by the user, to penalise any violation of \eqref{hard margin equation}.
Then, minimising the Lagrangian is equivalent to solving the following least square approximation \cite{QSVM},
\begin{equation}
    F
    \begin{pmatrix}
    b\\
    \vec{\alpha}
    \end{pmatrix}
=
\begin{pmatrix}
0 & \vec{1}^T\\
\vec{1} & K+\gamma^{-1}\mathbf{1} 
\end{pmatrix}
    \begin{pmatrix}
    b\\
    \vec{\alpha}
    \end{pmatrix}
    =    \begin{pmatrix}
    0\\
    \vec{y}
    \end{pmatrix}.
    \label{system of equations}
\end{equation}
Here is where \cite{QSVM} uses quantum linear algebra techniques to solve \eqref{system of equations}, as we will see later, and obtain the result $\ket{b,\vec{\alpha}}$.
To classify a new point $\ket{\vec{x}}$, one constructs, making use of qRAMs,
\begin{equation}
    \ket{\Tilde{u}}:=\frac{1}{\sqrt{N_{\Tilde{u}}}}
    \left(b\ket{0}\ket{0}+\sum_{k=1}^n\alpha_k|\vec{x}_k|\ket{k}\ket{\vec{x}_k}\right)
    \label{u tilde}
\end{equation}
and 
\begin{equation}
    \ket{\Tilde{x}}:=\frac{1}{\sqrt{N_{\Tilde{x}}}}
    \left(\ket{0}\ket{0}+\sum_{k=1}^n|\vec{x}|\ket{k}\ket{\vec{x}}\right).
    \label{x tilde}
\end{equation}
Remember that qRAMs perform $\sum_k\beta_k\ket{k}\rightarrow \sum_k\beta_k\ket{k}\ket{\vec{x}_k}$, for $\beta_k$ arbitrary amplitudes. Using the qRAM proposed in \cite{QRAM_bucket_brigade} only $O(\log n)$ gates are activated, although $O(n)$ should be present in the circuit.

If the previous is possible, one should also be able to prepare the state $\ket{\psi}= 1/\sqrt{2}(\ket{0}\ket{\Tilde{u}}+\ket{1}\ket{\Tilde{x}})$ and measures the probability of the ancilla being in state $\ket{-}=1/\sqrt{2}(\ket{0}-\ket{1})$,
\begin{equation}
    P=\frac{1}{2}(1-\braket{\Tilde{u}|\Tilde{x}}). \label{probability of 0}
\end{equation}
This is equivalent to performing a Hadamard gate over the first ancilla register and measuring the probability of obtaining $\ket{1}$, and is called Swap Test \cite{SupervisedQML}. If $P\geq 1/2$, $y(\vec{x})=1$, otherwise $y(\vec{x})=-1$.


\section{\label{sec:Main algorithm} Main algorithm}

As we have seen, our aim is to sample points to be added to the training set, $\vec{x}_{n+1}$, with the highest possible informativeness, defined as $P_c(\vec{x}_{n+1})|\vec{w}_{n+1}|$, where the first term indicates the probability that a given point is in a class, and the second measures how much that would improve the classifier.

Thus, we will employ the strategy explained in the Algorithm \ref{Main algorithm} to sample from the $1/C$ most relevant points.

\begin{figure}
\begin{algorithm}[H]
\caption{Active learning against adversarial examples.}\label{Main algorithm}
\begin{algorithmic}[1]
\Procedure{Active learning against adversarial examples}{}
\State\hskip-\ALG@thistlm \textbf{To find a point $\vec{x}$ that is in the $1-1/C$ quantile of `informativeness' with probability $1-e^{-\beta}$, iterate $O(C\beta)$ times:}
\State Uniformly at random sample a point $\vec{x}_{n+1}$ in the space.
\State Evaluate, as explained in section \ref{sec:Probability}, the probability that a point is in a given class $c$. We get $P_c(\vec{x}_{n+1})$. \label{step probability}
\State Solve the linear system of equations \eqref{system of equations} using the procedure of \cite{QSVM}, to obtain the state $\ket{b,\vec{\alpha}}_{n+1}$. It is explained in section \ref{sec:QSVM}. \label{step QSVM}
\State Perform \eqref{second matrix system}, using the Chebyshev approach taken from \cite{QLSAchilds}, to calculate $\ket{\vec{w}_{n+1}}$. Then use Amplitude Estimation to estimate $|\vec{w}_{n+1}|$. The procedure to calculate $\ket{\vec{w}_{n+1}}$ is described in section \ref{sec:norm}, and its norm in appendix \ref{sec:Norm_appendix}. \label{step norm}
\State Calculate the informativeness $P_c(\vec{x}_{n+1})\cdot |\vec{w}_{n+1}|$ of point $\vec{x}_{n+1}$. If it is higher than the previous best `informativeness', save the pair $(\vec{x}_{n+1},P_c(\vec{x}_{n+1})\cdot |\vec{w}_{n+1}|)$ to memory substituting the previous best pair of values.
\State After the $O(\beta C)$ iterations, output the saved best pair.
\EndProcedure
\end{algorithmic}
\end{algorithm}
\end{figure}

\subsection{\label{sec:Probability}Calculating $P_c(\vec{x}_{n+1})$}

The first thing we should care about is calculating $P_c(\vec{x}_{n+1})$, the probability that point $\vec{x}_{n+1}$ is in class $c$. 

The simplest way to calculate the probability would be to solve the SVM, calculate the distance from the point $\vec{x}_{n+1}$ to the decision boundary, and apply an activation function that converts the distance to a probability.
Solving the SVM can be done using \cite{QSVM} and reading each entry using Amplitude Estimation. It would return the vector $(b,\vec{\alpha})$, which can be used to create $\vec{w}=\sum_j \alpha_j \vec{x}_j$, and therefore the SVM.

A second, more elegant solution, is the following. In \cite{QSVM}, the authors propose a method to estimate, using the Swap Test and the output state $\ket{b, \vec{\alpha}}$ , to which class does a given point $\ket{\vec{x}}$ belong. They calculate that the success probability of measuring a $\ket{-}$ in the ancilla is
\begin{equation}
    P=\frac{1}{2}(1-\braket{\Tilde{u}|\Tilde{x}}),
\end{equation}
for $\ket{\Tilde{u}}$ and $\ket{\Tilde{x}}$ defined as in \eqref{u tilde} and \eqref{x tilde}.
The interesting thing to notice is that if $P>1/2$ the classification is in one class and if $P<1/2$ it is in the other. We modify this protocol slightly so that we perform Amplitude Estimation on this result instead of repeatedly measuring the expected value, which improves the complexity from $O(\epsilon^{-2})$ to $O(\epsilon^{-1})$. This allows us to obtain the amplitude $A_{\vec{x}_{n+1}} = \sqrt{P}$. Then, one may use $P$ as $P_c(\vec{x}_{n+1})$, or apply whatever activation function we consider appropriate to shape $P_c(\vec{x}_{n+1})$, like a sigmoid for example.

\subsection{\label{sec:QSVM} Quantum Support Vector Machine}

Now we focus on the main subroutine that our algorithm uses, which is mostly based in \cite{QSVM}, and that outputs $\ket{b,\vec{\alpha}}$, the solution to \eqref{system of equations}. This subroutine will be key to calculating $|\vec{w}_{n+1}|$, as one can use \eqref{w definition} to calculate $\vec{w}$, whose norm we will estimate later. In this section we will review how \cite{QSVM} performs Hamiltonian Simulation, and its role in the HHL algorithm \cite{HHL}, as well as some variations of that algorithm. We will also show the low rank approximation needed to achieve the exponential speedup, and the role of the condition number as a normalization factor in the estimation of $||(b,\vec{\alpha})||$. Finally, we will review the complexity of this calculation, which relies on the use of qRAMs to achieve an exponential speedup.

Thus, the first step is solving \eqref{system of equations}, which \cite{QSVM} does by using HHL algorithm \cite{HHL} with a different Hamiltonian simulation. This is due to the fact that in general the kernel matrix $K$ is dense, and the Hamiltonian simulation of \cite{HHL} is better suited for the sparse case.

Rather, the authors of \cite{QSVM} use a technique they had previously developed in \cite{QPCA} that allows for efficient Hamiltonian simulation of dense low-rank matrices. Let us explain how to do it. Suppose we have a given quantum state $\sigma$, and we are able to prepare $T$ copies of a density matrix $\rho$. One would like to perform the Hamiltonian simulation
\begin{equation}
    e^{-i t \rho}\sigma e^{i t \rho}. \label{matrix application}
\end{equation}
Repeated application of 
\begin{equation}
    \mathrm{tr}_1(e^{-i \Delta t S}\rho\otimes\sigma e^{i \Delta t S}) = \sigma -i \Delta t [\rho, \sigma] + O(\Delta t^2), \label{low rank hamiltonian simulation}
\end{equation}
$S$ the Swap operator, approximates
\begin{equation}
    e^{-i T\Delta t \rho}\sigma e^{i T\Delta t \rho},
\end{equation}
$T$ indicating the number of times we apply \eqref{low rank hamiltonian simulation}.
The authors of \cite{QPCA} observe that to simulate \eqref{matrix application} with precision $\epsilon^{-1}$ one has to repeat the process $ O(t^2\epsilon^{-1})$ times, each taking time $O(\log m)$ since that is the time it takes to prepare $\rho$ using a qRAM.

In our case $\rho$ is the matrix $F$, that will be represented by operator $\hat{F} = (J + K + \gamma^{-1} \mathbf{1})/ \mathrm{tr} F$, where
\begin{equation}
    J = \begin{pmatrix}
    0 & \vec{1}^T\\
    \vec{1} & 0
    \end{pmatrix}.
\end{equation}
Then, $e^{i\Delta t \hat{F}}= e^{i\Delta t J/\mathrm{tr} F}e^{i\Delta t K/\mathrm{tr} F}e^{i\Delta t \gamma^{-1} \mathbf{1}/\mathrm{tr} F}$. Simulating $J$ is described in \cite{childs2010relationship}, and the matrix $\gamma^{-1} \mathbf{1}$ is also easy. An important insight of \cite{QSVM} is how to prepare a state with density matrix $K/\mathrm{tr} F$. If we were able to prepare $K/\mathrm{tr} K$ then one would only need to correct for a scalar term $\mathrm{tr} K/ \mathrm{tr} F = O(1)$ in the simulation time. In the appendix A of \cite{QSVM} it is explained how to estimate $\mathrm{tr} K$, which can be used to estimate $\mathrm{tr} F = \mathrm{tr} K + \gamma \mathrm{tr}  \mathbf{1}$.

So, let us show how \cite{QSVM} prepares a density matrix state $K/\mathrm{tr} K$. To do so, prepare, using a qRAM:
\begin{equation}
    \ket{\chi} = \frac{1}{\sqrt{N_\chi}}\sum_i |\vec{x}_i| \ket{i}\ket{\vec{x}_i},
\end{equation}
Tracing out the second register prepares
\begin{equation}
\begin{split}
    \mathrm{tr}_2 (\ket{\chi}\bra{\chi}) &= \frac{1}{N_\chi} \sum_{i,j=1}^n\braket{\vec{x}_i|\vec{x}_j}|\vec{x}_i||\vec{x}_j|\ket{i}\bra{j}\\
    &= K/\mathrm{tr} K.
\end{split}
\end{equation}

The previous method allows to prepare the density matrix $K/\mathrm{tr} K$ needed to implement the Hamiltonian simulation of \eqref{low rank hamiltonian simulation}, in time complexity $O(\epsilon^{-1}t^2 k)$. Therefore, we can perform the Hamiltonian simulation necessary to implement the HHL algorithm efficiently.

The HHL algorithm consists on the following steps. First, one formally decomposes the state $\ket{0,\vec{y}}$ of \eqref{system of equations} in eigenvectors $\ket{0,\vec{y}} = \sum \beta_j \ket{u_j}$ of the Hamiltonian defined by matrix $F$ in \eqref{system of equations}. Next, one Phase Estimates the eigenvalues, using a Hamiltonian simulation algorithm such as the one explained above, obtaining $\sum \beta_j \ket{u_j}\ket{\lambda_j}$. Finally, one performs the controlled rotations
\begin{equation}
\begin{split}
    &\sum_j \beta_j \ket{u_j}\ket{\lambda_j}\ket{0}\rightarrow \\
    &\sum_j \beta_j \ket{u_j}\ket{\lambda_j}\left(\frac{1}{\lambda_j\kappa}\ket{1}+\sqrt{1-\frac{1}{\lambda_j^2\kappa^2}}\ket{0}\right),
\end{split} \label{HHL}
\end{equation}
and postselects in ancilla state $\ket{1}$ to obtain the solution to the system of equations \cite{HHL}.

The stated complexity of the original HHL algorithm is $O(\kappa^2 \epsilon^{-1}\text{poly}\log n)$, $\kappa$ the condition number. The complexity of the originally used Hamiltonian simulation algorithm is $O(d\epsilon^{-1}\text{poly}\log n)$, where $d$ the sparsity, for the sparse oracle access model \cite{HHL}. The complexity of the condition number has been lowered to $O(\kappa)$, as shown using the technique of Variable Time Amplitude Estimation in \cite{Ambainis}, which is optimal in this parameter, also proved in that article. On the other hand, relying on more efficient Hamiltonian simulation techniques, \cite{QLSAchilds} proposed a variation of the HHL algorithm that does not require Amplitude Estimation thus reducing the complexity from $O(\epsilon^{-1})$ to $O(\text{poly}\log \epsilon^{-1})$. Unfortunately, the Hamiltonian simulation described above still has complexity $O(\epsilon^{-1})$, making it impossible to reduce the complexity of the overall method further.

There is an alternative to the Hamiltonian simulation method that we have described. In \cite{Dense_hamiltonian_simulation}, a technique is introduced that allows to simulate a dense Hamiltonian with complexity $O(\sqrt{n}\text{poly}\log \epsilon^{-1})$, for Hamiltonians of any rank, and it can be used with the techniques from \cite{QLSAchilds}, to reduce the complexity on the precision from linear to polylogarithmic. The caveat is that it requires to use a special quantum-accessible data structure that is explained in the appendix \ref{sec:definitions} and plays the role of a Quantum Read-Only Memory, and has the complexity stated above, polynomial in $n$.

Coming back to our main discussion, once we have the state $\ket{b,\vec{\alpha}}$, we would like to recover the norm of vector $\vec{w}$. As a first step, we would like to obtain the norm $||(b,\vec{\alpha})||$. Suppose we are trying to solve $Ax=b$, where the largest eigenvalue of $A$ is less or equal to 1; else see the end of appendix \ref{sec:Norm_appendix} for a minor correction. As explained in appendix A of \cite{FEM}, we can calculate the norm of the solution $||x||$ using
\begin{equation}
    ||x||= \kappa ||b|| \sqrt{p_1}, \label{norm solution HHL}
\end{equation} 
where $\sqrt{p_1}$ represents the amplitude of the postselection ancilla of HHL algorithm being in the correct state, usually $\ket{1}$. This ancilla is used to perform the non-unitary part of the algorithm via a measurement. The reason for the previous equation \eqref{norm solution HHL} is because the acceptance probability of HHL scales as
\begin{equation}
    p_1 = \frac{||A^{-1}b||}{||b||\kappa^2}.
\end{equation}
Estimating such amplitude $\sqrt{p_1}$ thus requires of Amplitude Estimation \cite{Amplitude_estimation}, with cost $O(\epsilon^{-1})$. 

Now let us turn to the condition number of the system of equations that appears in the SVM we are solving, the condition number of the matrix $F$ in \eqref{system of equations}. Recall that the condition number is defined as
\begin{equation}
    \kappa= \frac{\sigma_{\max}}{\sigma_{\min}},
\end{equation}
where $\sigma_{\min}$ and $\sigma_{\max}$ are the minimum and maximum singular values respectively. The condition number will be important to correct the norm of the solution to the system of equations, as can be seen from \eqref{norm solution HHL}.
However, calculating it with the Quantum Singular Value estimation technique from \cite{kerenidis2017recommendation} might be too expensive. 

On the other hand, the authors of the quantum SVM article, \cite{QSVM}, propose that in the case where the kernel matrix has $O(1)$ eigenvalues of size $O(1)$, and $O(n)$ eigenvalues with values $O(1/n)$ as it is in our case, we can choose a condition number $\kappa_{\mathrm{eff}}=O(1)$ such that in the end we will get an additional error of order $O(1/\sqrt{n})$, in addition to $\epsilon$.  

This low rank approximation means that the algorithm in \cite{QSVM} only takes into account the eigenvalues $\lambda_i$ that are $\epsilon_K\leq \lambda_i \leq 1$. The main idea of the low rank approximation is filtering out those eigenvalues $\lambda< \kappa_{\mathrm{eff}}^{-1}$, so the final rotation of the HHL algorithm, indicated in \eqref{HHL}, is performed only if $\lambda_j > \kappa_{\mathrm{eff}}^{-1}$ and imposes $\kappa = \kappa_{\mathrm{eff}} = O(1)$.

In appendix C of the supplementary material of \cite{QSVM} the authors show that $||K-K_{\mathrm{eff}}|| = \sqrt{\sum_{\lambda_i = O(1/n)}\lambda_i^2}$, with $K_{\mathrm{eff}}$ the `filtered' low-rank approximation of $K$, $K_{\mathrm{eff}} = \sum_{\lambda_i \geq \kappa_{\mathrm{eff}}^{-1}}\lambda_i \ket{u_i}\bra{u_i}$. Since there are $O(n)$ eigenvalues of size $O(1/n)$, the induced error is of order $O(n^{-1/2})$. Then, one can use theorem 1 in the appendix of the HHL algorithm, \cite{HHL}, to show that, if $\ket{b,\vec{\alpha}}$ is the exact solution of the system of equation; and $\ket{b,\vec{\alpha}}_{\mathrm{eff}}$ the approximate one after postselecting on the subspace spanned by the eigenvalues $\lambda_j \geq \kappa_{\mathrm{eff}}^{-1}$ and the ancilla in state $\ket{1}$, then $||\ket{b,\vec{\alpha}}_{\mathrm{eff}}-\ket{b,\vec{\alpha}}|| = O(\epsilon +n^{-1/2})$. Thus, for relatively large $n$ the induced error with this approximation is small.

Overall, this implies that instead of \eqref{norm solution HHL}, we will have 
\begin{equation}
    ||x||= \kappa_{\mathrm{eff}} ||b|| \sqrt{p_1}, \label{norm solution QSVM}
\end{equation} 
In such case notice that we have imposed an effective condition number for all candidate points $\vec{x}_{n+1}$: $\kappa_{\vec{x}_{n+1}}=\kappa_{\mathrm{eff}}$, so we no longer have to care about the condition number: it will be the same scale factor for all points $\vec{x}_{n+1}$, and thus without relevance. The same will happen also for $||b||$ in \eqref{norm solution QSVM} since $b = (0, y_1,...,y_{n+1})^T$, with $y_{n+1}$ the same for all candidate points.
Since we suppose $n$ is large enough, the additional additive error we introduce will be small, of order $O(n^{-1/2})$ according to the appendix C in the supplementary material of \cite{QSVM}. Thus, the complexity of the algorithm is $O(\epsilon^{-3}\kappa_{\mathrm{eff}}^3\log(mn))$, where $m$ is the dimension of the space. 

To finish this section, let us recall how to calculate the complexity of the algorithm of \cite{QSVM}. They claim that \eqref{low rank hamiltonian simulation} has error $\epsilon = \Tilde{O}(||\hat{F}||^2\Delta t^2 )$, $||\hat{F}||$ the Frobenius norm of $\hat{F}$, and it is repeated $T$ periods. Thus, $\epsilon =\Tilde{O}(||\hat{F}||^2 \Delta t^2 T)= \Tilde{O}(||\hat{F}||^2 t^2/T)$, for $\Delta t = t/T$. This implies a simulation cost $T = \Tilde{O}(t^2 ||\hat{F}||^2 \epsilon^{-1})$, when implementing \eqref{low rank hamiltonian simulation} is taken at unit cost.

On the other hand, HHL requires phase estimating the eigenvalues, so the relative error of $\lambda^{-1}$ can be indicated as $\epsilon = O(1/\lambda t)\leq O(\kappa_{\mathrm{eff}}/t)$. Therefore $t = O(\kappa_{\mathrm{eff}} \epsilon^{-1})$, and substituting in the previous paragraph, $T = O(||\hat{F}||^2 \kappa_{\mathrm{eff}}^2 \epsilon^{-3} \log(mn))$. An additional $\kappa_{\mathrm{eff}}$ is needed to perform postselection on the result. Thus, the overall complexity of the algorithm is $O(o ||\hat{F}||^2  \kappa_{\mathrm{eff}}^3\epsilon^{-3}\log(mn))$, for a kernel of order $o$.

If instead of using the basic HHL technique we had decided to use more advanced ones, such as the Fourier approach described in \cite{QLSAchilds} (see appendix \ref{sec:definitions}), the complexity would have been $\Tilde{O}(T\alpha)$, for $T = \Tilde{O}(||\hat{F}||^2t^2\epsilon^{-1})$, $t = \Tilde{O}(\kappa)$ and $\alpha = \Tilde{O}(\kappa)$, ignoring polylogarithmic factors \cite{QLSAchilds}. Overall, the complexity would be $\Tilde{O}(||\hat{F}||^2\kappa^3\epsilon^{-1})$.


\subsection{\label{sec:norm}  The norm of $|\vec{w}_{n+1}|$}

The aim of the previous section was to obtain $\ket{b,\vec{\alpha}}$ and its norm. The aim of this one is to calculate $\vec{w}$ from that result. We will see that one may do that using a Fourier expansion, at cost $O(\epsilon^{-1})$; or better, use of Chebyshev series as in \cite{QLSAchilds}, at cost polylogarithmic in all variables, provided the use of qRAMs. In appendix \ref{sec:Norm_appendix} we explain how to calculate the norm of $|\vec{w}_{n+1}|$ from $\vec{w}$, which only requires performing Amplitude Estimation once, and normalizing with some factors.

The first, trivial, idea we had to calculate $\ket{\vec{w}_{n+1}}$ is reading all the entries of $\ket{b,\vec{\alpha}}$ using Amplitude Estimation. Then one may use classical computing to compute $\vec{w}= \sum_{i=1}^{n+1} \alpha_i \vec{x}_i $, and finally its norm $|\vec{w}_{n+1}|= \sqrt{\sum_{i=1}^m (\sum_{j=1}^{n+1} x_{ij,x_{n+1}} \alpha_{j,x_{n+1}})^2}$. Amplitude Estimation does not immediately recover the sign of each $\alpha_i$, but finding it is not complicated either. Once we know $|\alpha_i|$ and $|\alpha_j|$ we can prepare the state $C_{ij} (|\alpha_i| \ket{j} + |\alpha_j| \ket{i})$ and perform a Swap Test with the solution vector $\ket{b, \vec{\alpha}}$. If the relative sign of entries $i$ and $j$ is the same we will get a nonzero result proportional to $2C_{ij} \alpha_i \alpha_j$, but if the relative sign is opposite, the dot product will cancel out \cite{QIP}. Establishing the relative sign of two entries is enough, since the norm will not care about the global sign.

However, the procedure described above is time consuming, as one needs to prepare the solution of the system of equations \eqref{system of equations} $O(n+1)$ times in order to read all the entries of the solution and calculate their relative sign. Also, calculating $\vec{w}$ takes complexity $O(nm)$ using classical computing, so, instead of that we will try to prepare $\ket{\vec{w}_{n+1}}$ from $\ket{b,\vec{\alpha}}$, and estimate the change in the norm of the result. Notice in the first place that the Quantum SVM article \cite{QSVM} is able to classify a point without first calculating $\ket{\vec{w}_{n+1}}$, as indicated in \eqref{probability of 0} and subsequent paragraph. Also, it is not clear how to prepare $\ket{\vec{w}_{n+1}}$ from the state in \eqref{u tilde}. 

Instead, since $\vec{w}= \sum_i \alpha_i \vec{x}_i$, this suggest multiplying the solution vector from the Quantum SVM, $\ket{b, \vec{\alpha}}$ by a matrix operator whose entries are $x_{i,j}$. Explicitly, 
\begin{equation}
\begin{split}
A\ket{b,\vec{\alpha}} = 
\begin{pmatrix}
0 & x_{1,1} & ... & x_{1,m}\\
\vdots &    &     & \vdots \\
0 & x_{n+1,1} & ... & x_{n+1,m}
\end{pmatrix}
\begin{pmatrix}
b\\
\alpha_1\\
\vdots \\
\alpha_{n+1}
\end{pmatrix}\\
= 
\begin{pmatrix}
\sum_{i=1}^m  \alpha_j x_{1,j}\\
\vdots \\
\sum_{i=1}^m  \alpha_j x_{n+1,j}\\
\end{pmatrix}
= \vec{w}.
\end{split} \label{second matrix system}
\end{equation}
Performing this operation classically, requires $O(nm)$ operations. Can we perform this operation in a quantum way?

Our strategy to calculate the previous matrix-vector product will take inspiration from the algorithms presented in \cite{QLSAchilds} to simulate $A^{-1}$. 
The fact that in our case $A$ is not necessarily square should not pose any problem, since one can write
\begin{equation}
\begin{blockarray}{ccc}
 & m+1 &  n+1 \\
\begin{block}{c(cc)}
     m+1& \mathbf{0} & A^T   \\
     n+1 & A & \mathbf{0} \\
\end{block}
\end{blockarray}
    \quad     
    \begin{blockarray}{c}
    \\
    \begin{block}{(c)}
    b\\
    \vec{\alpha}\\
    \vec{0}_{n+1}\\
    \end{block}
    \end{blockarray}\quad = \quad
    \begin{blockarray}{c}
    \\
    \begin{block}{(c)}
    \vec{0}_{m+1}\\
    \vec{w}\\
    \end{block}
    \end{blockarray}\quad.
    \label{new Matrix system}
\end{equation}
Let us call $M$ the matrix from this previous equation. The first thing we have to do is to attach a register to $\ket{b,\vec{\alpha}}$, that we set to $\ket{0}$ for the entry $b$ (that is, the original register is in state $\ket{0}$) and $\ket{1}$ otherwise. This will allow us to substitute any of the $0$s in the matrix in \eqref{second matrix system}, but avoiding interference of the entry of $b$ with those of $\vec{\alpha}$. The aim of this is to avoid $M$ being ill-conditioned, although it is only a technical detail. Then all the singular values fall in $[1/\kappa_M, 1]$, and we can rewrite
\begin{equation}
    \begin{pmatrix}
    b\\
    \vec{\alpha}\\
    \vec{0}_{n+1}
    \end{pmatrix}
    \rightarrow
    \begin{pmatrix}
    b\\
    \vec{0}_{m}\\
    \vec{0}_{n+1}
    \end{pmatrix}
    \otimes \ket{0}
    +
        \begin{pmatrix}
    0\\
    \vec{\alpha}\\
    \vec{0}_{n+1}
    \end{pmatrix}\otimes\ket{1}.
\end{equation}
Notice that by linearity of quantum mechanics the unitary we were going to apply to the the left-hand side is applied individually to each of the two terms in the decomposition. Later we will only care about the case when the additional ancilla we have added is in state $\ket{1}$.

With this remark, the question becomes twofold: how to decompose $M$ into a Linear Combination of Unitaries (LCU), as it is done with $A^{-1}$ in \cite{QLSAchilds}, and how to simulate those linear operators efficiently. Notice that now $e^{-iM}$ is a unitary operator.

The first question is how to decompose $M$ in a linear combination of unitaries. If we choose, as it is frequently done, unitary operators of the form $e^{-iMt}$, we can write $M = \sum_i \beta_i e^{-iMt_i}$. Then, as both sides of the equality are diagonal in the same basis, one can consider this is equivalent to writing $x = \sum_i \beta_i e^{-ixt_i}$. The first idea is then to decompose $x$ as a Fourier series. Specifically, decomposing $x$ gives us
\begin{equation}
\begin{split}
    x &= \sum_{t=1}^\infty \frac{(-1)^t}{t} 2\sin(tx)\\
     &= i\sum_{t=1}^\infty  \frac{(-1)^t}{t} (e^{-itx}-e^{itx}).
\end{split}
\end{equation}
Since we cannot perform the entire summation, the question is how to truncate the series maintaining a given precision $\epsilon^{-1}$. Call $S_{\bar{t}}(x)$ the truncated series up to $\bar{t}$. Since all the singular values will be in $[1/\kappa_M,1]$, with $\kappa_M$ the condition number of $M$, then, we want to study the series in the range $[0,1]$. Notice that $f(x)=x$ is Lipschitz-continuous with constant $ 1$, so can use a result from \cite{jackson1930theory} to bound the error incur when truncating the series. According to Corollary I in \cite{jackson1930theory}, we can say that $|f(x)-S_{\bar{t}}(x)|\leq (a \log{\bar{t}})/\bar{t}$, where $a$ is a constant. Thus, if we want $|x-S_{\bar{t}}(x)|\leq\epsilon$, it is enough to choose $\bar{t}/\log \bar{t} = O(\epsilon^{-1})$. So the complexity is almost linear in $\epsilon^{-1}$, as we have to Hamiltonian-simulate $M$ during time $\bar{t}$. Furthermore, if we chose a Hamiltonian simulation algorithm, as $M$ is dense, the cost will be relatively large, $\Tilde{O}(\sqrt{n+m})$ if we were to use \cite{Dense_hamiltonian_simulation} for instance. In any case we would not obtain an exponential advantage, since Hamiltonian simulation of dense Hamiltonians is not expected to have complexity $\log n$ in general \cite{childs2009limitations}. Also notice that we cannot use the trick from the previous section, since that requires $M$ to be a density matrix and $\mathrm{tr} M = 0 \Rightarrow M$ is not a density matrix.

A second, better approach, also inspired by \cite{QLSAchilds}, is to decompose the matrix in Chebyshev polynomials, $x = \sum \gamma_i \mathcal{T}_{i}(x)$. In fact, this decomposition in polynomials is quite simple since the first Chebyshev polynomial is
\begin{equation}
    \mathcal{T}_{1}(x) = x,
\end{equation}
so the expansion will have a single term and will be exact, which is very positive. Let us now look at the cost of simulating a Chebyshev polynomial for a $\bar{n}\times \bar{n}$ $d$-sparse matrix $M$. One does that via Quantum Walks. Define a quantum walk as set of states $\{\ket{\psi_j}\in \mathbb{C}^{2\bar{n}}\otimes\mathbb{C}^{2\bar{n}}; j\in [\bar{n}]\}$ \cite{QLSAchilds}:
\begin{equation}
\begin{split}
    &\ket{\psi_j}:= \ket{j}\otimes \\
    &\frac{1}{d} \sum _{k\in [\bar{n}]; A_{jk}\neq 0}\left(\sqrt{A^*_{jk}}\ket{k}+\sqrt{1-|A_{jk}|}\ket{k+\bar{n}}\right).
\end{split}
\end{equation}
Then, one defines the isometry
\begin{equation}
    T= \sum_{j\in [\bar{n}]}\ket{\psi_j}\bra{j},
\end{equation}
and $S$ the Swap operator $S:\ket{j,k}\rightarrow\ket{k,j}$. The Walk operator is $W := S(2TT^\dagger - \mathbf{1})$. With these definitions, and using Lemma 15 of \cite{QLSAchilds} enunciated in the appendix \ref{sec:definitions}, one can prove that for any state $\ket{\psi}$, substituting $T$ by its corresponding unitary $T_U$ \cite{QLSAchilds}
\begin{equation}
\begin{split}
    &T_U^\dagger W T_U \ket{0^{\lceil \log 2\bar{n}\rceil +1}}\ket{\psi}\\& = \ket{0^{\lceil \log 2\bar{n}\rceil +1}} \mathcal{T}_{1}(\overline{M})\ket{\psi} + \ket{\Phi^{\perp}}\\
    &=\ket{0^{\lceil \log 2\bar{n}\rceil +1}} \overline{M}\ket{\psi} + \ket{\Phi^{\perp}},\label{Quantum walk}
\end{split}
\end{equation}
for $\overline{M}= M/d$, and $\ket{\Phi^{\perp}}$ such that $(\ket{0^{\lceil \log 2\bar{n}\rceil +1}}\bra{0^{\lceil \log 2\bar{n}\rceil +1}}\otimes \mathbf{1})\ket{\Phi^{\perp}}=0$. 

Let us first study the cost of applying $T$. According to Lemma 10 in \cite{berry2015hamiltonian}, to implement $T$ one must apply $\log d $ Hadamard gates, followed by $O(1)$ calls to an oracle $O_F$ that outputs the column index of the $l$-th non-zero element of a row $j$ of the matrix, $O_F:\ket{j,l}\rightarrow\ket{j, f(l,j)}$; and the sparse-access oracle $O_M$, that outputs entry $M_{jk}$ on input $(j,k)$. For simplicity and without loss of generality let us assume that in \eqref{new Matrix system} all entries of $A$ are non-zero. If such is the case, then the cost of the oracles is not high either. The first oracle outputs $O_F:\ket{j,l}\rightarrow\ket{j, l}$ or $O_F:\ket{j,l}\rightarrow\ket{j, l+n+1}$ depending on whether we are on the last $n+1$ rows or in the first $m+1$, respectively. The second oracle performs $O_M:\ket{j,k}\ket{z}\rightarrow\ket{j,k}\ket{z\oplus M_{jk}}$, which is nevertheless no more restrictive than the qRAM model that we were using to solve the quantum SVM.

What is the cost of implementing $W$ then? In the proof of their theorem 4, \cite{QLSAchilds} cites Lemma 10 from \cite{berry2015hamiltonian} to explicitly state that if we want to simulate $W$ with error $\epsilon'$, its complexity and the complexity of \eqref{Quantum walk} is $O(\log \bar{n} + \log^{2.5}(\kappa d / \epsilon'))$. This is not surprising since the cost only depends on $T$, its inverse $T^\dagger$, and a Swap gate. Notice that our matrix is dense, so $d = O (\bar{n})=O(n+m)$, but still the cost of this process is polylogarithmic in all variables! In fact this method seems applicable to any matrix $M$, sparse or not, whose entries are accessible through a sparse-access oracle in the form of a qRAM. One may even say that the qRAM is hiding the cost, since it requires time to prepare the data and also has $O(n)$ quantum gates even if it only uses a small number of them. Notice also that the linear dependence on $d$ and $\kappa$ for the Quantum Linear System Algorithm presented in \cite{QLSAchilds} does not come from implementing $W$ but because one has to apply it $O(d\kappa)$ times in the series expansion of $A^{-1}$. In our case we only have to apply $W$ once, so we do not have such linear complexity term.


Finally, we measure the amplitude of $\ket{1}$ of the ancilla that we attached to the state to mark the entry of $b$ after \eqref{second matrix system}. The amplitude of such state is, after taking into account a normalization factor for $M$ and $||(b,\vec{\alpha})||$, precisely $|\vec{w}_{n+1}|=\sqrt{\sum_{i=1}^m (\sum_{j=1}^{n+1} x_{ij,x_{n+1}} \alpha_{j,x_{n+1}})^2}$, what we were looking for. We can measure that amplitude using Amplitude Estimation.

In appendix \ref{sec:Norm_appendix} we explain in greater detail how to correct the value extracted with Amplitude Estimation, to find the norm $|\vec{w}_{n+1}|$ provided either the Linear Combination of Unitaries with the Fourier series, or the Chebyshev polynomial approach. The cost of Amplitude Estimation is $O(\epsilon^{-1})$.
Having calculated $|P_c(\vec{x}_{n+1})|$ and $|\vec{w}_{n+1}|$, the informativeness of a given point is just the product of both.

\subsection{\label{sec:Target} Finding the target point $\vec{x}_{n+1}$}

Now that we know how to calculate the `informativeness' of a given point, let us explain why the Algorithm \ref{Main algorithm} finds, with probability $1-e^{-\beta}$ a point that is in the quantile $1-1/C$ of informativeness. The arguments presented here are the same to those in appendix \ref{sec:classical_strategies}. 

If we sample uniformly at random, and calculate the informativeness of $c$ points, the probability that none of those points is in the quantile $1-1/C$ is clearly 
\begin{equation}
    p = \left(1-\frac{1}{C}\right)^c.
\end{equation}
We want to make such probability $p< e^{-\beta}$, so that at least one is in the $1-1/C$ quantile. That means
\begin{equation}
    c > \frac{-\beta}{\log \left(1-1\frac{1}{C}\right)}.
\end{equation}
Since 
\begin{equation}
    \lim_{C\rightarrow\infty} \frac{\frac{-1}{\log \left(1-\frac{1}{C}\right)}}{C} = 1,
\end{equation}
it is clear then that the number of sampled points should be $c = O(\beta C)$. This explains why our procedure is efficient in these variables. More information on alternative strategies can be found in the appendix \ref{sec:classical_strategies}.

\section{\label{sec:results}Main results}

In this article we propose a theoretical framework that allows us to think of active learning as sampling the most promising new points to be classified, so that the minimum distance between classes can be found, and theorem \ref{Teorema 5} used. We also propose a quantum algorithm that would allow us to perform the sampling efficiently with a exponential advantage over what is achieved classically.

In the quantum algorithm we do not use neural networks but rather a SVM \cite{SVM}, although it might be possible to use the general strategy in other setups using quantum neural networks.

Let us explain further the complexity comparison with the classical algorithms. One may argue that a good quantum strategy to solve Problem \ref{Informative sampling problem} would be to use Amplitude Estimation and bisection search to find a threshold for the quantile $1-1/C$, for example the $1\%$ most relevant points, in which case $C=100$. Then, one may use Amplitude Amplification to find one of the points in such quantile. This achieves an exponential advantage over the classical case when one wants to find such points using classical computing with certainty. In the usual case of Amplitude Amplification and Amplitude Estimation we have an oracle that tells us whether an element is marked or not, and this yields a quadratic advantage with respect to the classical case. However, in our situation, being a good point or not depends on its relative `informativeness' with respect to other points. This means that classically, in order to assert with certainty that a given point is within the top $1/C=1\%$ quantile, one should first calculate the informativeness of $0.99N$ points, out of $N$ points. This is clearly prohibitive, since $N = O(l^m)$, where $m$ is the dimension of the space, and $l$ its discretisation. Notice for example that for a $n_0\times n_0$ image, the dimension of the image would be $m = 3 n_0^2$, due to the three colours or channels needed to define each pixel. 

In contrast, if we want to solve this problem using Amplitude Estimation we do not incur in such cost. What we do is find, using bisection and Amplitude Estimation, an informativeness threshold above which there are only $1/C$ of the most informative points. Once that is the case, we can mark those points and use Amplitude Amplification to find them \cite{Amplitude_estimation}. The cost will then be $O(C\epsilon^{-1}\beta)$ for a given precision $\epsilon$ in the threshold, and success probability $1-e^{-\beta}$, and crucially $C$ independent of $N$. Amplitude Amplification also bears a cost $O(\sqrt{C})$.

On the other hand, one can also find probabilistic classical strategies that output a point that is in the $1-1/C$ quantile with probability exponentially high, $1-e^{-\beta}$, and the complexity would be polynomial in parameters $C$ and $\beta$. Since our quantum strategy would also have some exponentially small probability of failure, and the complexity would also be polynomial in those parameters, the advantage would be unclear. We discuss this point in depth in the appendix \ref{sec:classical_strategies}, for a problem that generalises our own, and find that this advantage would be quadratic in some parameters.

However, the algorithm that we present here achieves an exponential advantage over its classical counterpart. The speedup is due to a faster calculation of the `informativeness' of a given point, thanks to the use of qRAMs, the algorithm for quantum SVM \cite{QSVM}, and an approach to perform matrix-vector products using Chebyshev polynomials taking advantage of techniques developed in \cite{berry2015hamiltonian} and \cite{QLSAchilds}; in conclusion, a better use of quantum linear algebra techniques. Our algorithm here has overall complexity $\Tilde{O}(C\epsilon^{-1}\beta(o||F||^2\kappa_{\mathrm{eff}}^{-3}\epsilon^{-3}+ \mathrm{poly}\log (\kappa_M\epsilon^{-1}(n+m))))$, where $||F||$ is the Frobenius norm of the matrix in \eqref{system of equations}, $o$ the order of the kernel, $\kappa_{\mathrm{eff}}=O(1)$ an effective condition number that we choose, $\epsilon$ the precision on the `informativeness' of a given point, $1-e^{-\beta}$ the success probability with which one wants to be sure that the chosen point is in the $1-1/C$ quantile, and $\kappa_M$ the condition number of matrix M appearing in \eqref{new Matrix system}. 

In contrast, a classical algorithm would have polynomial complexity in $n$ and $m$ due for instance to matrix-vector product \eqref{new Matrix system}, or the final calculation of the norm of $|\vec{w}|$, central to our discussion.

In the next subsection we lay out the general strategy of our paper to solve this problem.



\section{\label{sec:Complexity} Complexity}

In this section we want to give a calculation of the complexity of the proposed quantum algorithm. 

The first step is calculating the probability that a given point of the high dimensional space is in a certain class, $P_c(\vec{x}_{n+1})$. This implies solving a quantum SVM \cite{QSVM}, with complexity $O(o ||F||^2\epsilon^{-3} \kappa_{\mathrm{eff}}^{-3} \log (mn))$, where $o$ is the order of the kernel, $||F||$ is the Frobenius norm of the matrix that appears in \eqref{system of equations}, $m$ the dimension of the space, $n$ the number of already classified points, and $\kappa_{\mathrm{eff}}$ is an effective condition number which can be taken $O(1)$ for this problem. 

Reading out the probability using Amplitude Estimation in the Swap Test costs an additional multiplicative $O(\epsilon^{-1})$, so the overall complexity of this step is $O(o ||F||^2 \epsilon^{-4} \kappa_{\mathrm{eff}}^{-3} \log (mn))$.  Since we are assuming that we have scaled the matrix to obtain the largest eigenvalue equal to 1, the procedure of \cite{QSVM} implies that the smallest one that we take into consideration is $\lambda^{-1} = \kappa_{\mathrm{eff}}$. The low-rank approximation induces an additional error of $O(n^{-1/2})$, as it is explained in the supplementary material of the original article \cite{QSVM}. 

Next step is calculating the Quantum SVM for the system with the added point, again with complexity $O(o ||F||^2 \epsilon^{-3} \kappa_{\mathrm{eff}}^{-3} \log (mn))$.
Then, preparing $\ket{\vec{w}_{n+1}}$ has complexity $O(\text{poly}\log (\kappa_M \epsilon^{-1}(n+m)))$, as reported in the corresponding section. Calculating $|\vec{w}_{n+1}|$ also have an additional multiplicative complexity of $\Tilde{O}(\epsilon^{-1})$ due to the Amplitude Estimation, given the procedure in appendix \ref{sec:Norm_appendix}.

Finally, if we are interested in finding a point in the $1-1/C$ quantile of `informativeness' with probability $1-e^{-\beta}$, then the procedure explained in section \ref{sec:Target} implies iterating the procedure over $O(C\beta)$ points and selecting the best.

Thus, in general, the complexity will be $\Tilde{O}(C\beta\epsilon^{-1}(o||F||^2 \kappa_{\mathrm{eff}}^{-3}\epsilon^{-3}+ \text{poly}\log (\kappa_M\epsilon^{-1}(n+m))))$, with error $\Tilde{O}(\epsilon + n^{-1/2})$. The $n^{-1/2}$ comes from the low-rank approximation of \cite{QSVM}, as it is explained on its appendix C.

\section{\label{sec:Conclusion}Conclusion}
In the previous sections we have seen that provided the use of qRAMs, it could be possible to establish a polynomial-complexity sampling procedure which lets us know what are the most promising points to be added to the training set in order to get a better approximation of the minimum distance between classes, $r_p$. Recall that knowing $r_p$ is a requisite to applying a $\delta$-cover and, using theorem \ref{Teorema 5}, avoid adversarial examples. This protocol is heuristic, which means that in order to check its actual performance we have to run it in a realistic quantum computer.

In this paper we have presented a procedure that allows to solve this problem efficiently, in time polylogarithmic in the dimension of the space $m$ and number of already classified points, $n$; and polynomial in $C$, $\beta$ and the precision $\epsilon^{-1}$ of the informativeness estimate.

There are nevertheless two important shortcomings of our algorithm. The first is that, since we rely heavily on reference \cite{QSVM}, and they use qRAMs to achieve their speedup, we also need qRAMs. The second is that there are some minor parts of the general algorithm that will necessarily have polynomial complexity in the dimension of the space $m$. Those are: sampling initial points, and loading their value in the qRAM. We strongly believe that this should nevertheless be no problem. The sampling procedure for instance will hardly be inefficient and can be easily done while the previous point is being processed. The loading in the qRAM is perhaps more expensive, but may be done while the point is asked to be classified by a human. Therefore we believe that the complexity result is still very valid in practice.

Finally, we want to recall that our work relies heavily on reference \cite{QSVM}. At the time that article was written, there was no known method to solve low-rank linear systems of equations in polylogarithmic time. However, \cite{chia2018quantum,arrazola2019quantum} recently proposed a method that achieves precisely that, complexity $O(||F||^6 k^6 \kappa^{16} \epsilon^{-6})$, $k$ the rank. It relies on the techniques developed by Ewin Tang \cite{ewin-tang}, which found a classical algorithm taking inspiration from \cite{Q_recommendation_sys}. Therefore, one could also solve the linear system of equations that appear in this case in polylogarithmic time in the dimension $n$. In the published version of this article it is argued that the algorithm could probably easily be dequantized, following such techniques, since the most complicated part is the solution of the system of equations. Here, we have added appendix \ref{sec:dequantization} to explain how to do it.

In conclusion, in this article we have framed a possible solution to adversarial examples using active learning, and a sampling methodology to find the most important points that could be added to the training set. We have also presented an algorithm that provides an exponential advantage over its classical counterpart, provided the use of qRAMs.

\section*{Acknowledgements}
We would like to thank Santiago Varona for useful comments on the manuscript, as well to Jaime Sevilla, Nikolas Bernaola and Javier Prieto for pointing us to useful statistic results for Appendix \ref{sec:classical_strategies}.
We acknowledge financial support from the Spanish MINECO grants MINECO/FEDER Projects FIS 2017-91460-EXP, PGC2018-099169-B-I00 FIS-2018 and from CAM/FEDER Project No. S2018/TCS-4342 (QUITEMAD-CM). The research of M.A.M.-D. has been partially supported by the U.S. Army Research Office through Grant No.  W911NF-14-1-0103. P. A. M. C. thanks the support of a FPU MECD Grant. 

\bibliographystyle{ieeetr}
