A quantum extension of SVM-perf for training nonlinear SVMs in almost linear time

Training a soft-margin $\ell_{1}$-SVM with parameter $C>0$ corresponds to solving the following optimization problem:

Note that, following [1], we divide $\sum_{i}\xi_{i}$ by $m$ to capture how $C$ scales with the training set size. The trivial case $\Phi(\mathbf{x})=\mathbf{x}$ corresponds to a linear SVM, i.e. a separating hyperplane is sought in the original input space. When one considers feature maps $\Phi(\mathbf{x})$ in a high dimensional space, it is more practical to consider the dual optimization problem, which is expressed in terms of inner products, and hence the kernel trick can be employed.

This is a convex quadratic program with box constraints, for which many classical solvers are available, and which requires time polynomial in $m$ to solve. For instance, using the barrier method [30] a solution can be found to within $\varepsilon_{b}$ in time $O(m^{4}\log(m/\varepsilon_{b}))$. Indeed, even the computation of the kernel matrix $K$ takes time $\Theta(m^{2})$, so obtaining subquadratic training times via direct evaluation of $K$ is not possible.

Joachims [1] showed that an efficient approximation algorithm - with running time $O(m)$ - for linear SVMs could be obtained by considering a slightly different but related model known as a structural SVM [31], which makes use of linear combinations of label-weighted feature vectors:

With this notation, the structural SVM primal and dual optimization problems are:

Whereas the original SVM problem OP 1 is defined by $m$ constraints and $m$ slack variables $\xi_{i}$, the structural SVM OP 3 has only one slack variable $\xi$ but $2^{m}$ constraints, corresponding to each possible binary vector $\mathbf{c}\in\{0,1\}^{m}$. In spite of these differences, the solutions to the two problems are equivalent in the following sense.

While elegant, Joachims’ algorithm can achieve $O(m)$ scaling only for linear SVMs — as it requires explicitly computing a set of vectors $\{\Psi_{\mathbf{c}}\}$ and their inner products — or to shift-invariant or low-rank kernels where sampling methods can be employed. For high dimensional feature maps $\Phi$ not corresponding to shift invariant kernels, computing $\Psi_{\mathbf{c}}$ classically is inefficient. We propose instead to embed the feature mapped vectors $\Phi(\mathbf{x})$ and linear combinations $\Psi_{\mathbf{c}}$ in the amplitudes of quantum states, and compute the required inner products efficiently using a quantum computer.

In Section 3 we will formally introduce the concept of a quantum feature map. For now it is sufficient to view this as a quantum circuit which, in time $T_{\Phi}$, realizes a feature map $\Phi:\mathbb{R}^{d}\rightarrow\mathcal{H}$, with maximum norm $\max_{\mathbf{x}}\left\|\Phi(\mathbf{x})\right\|=R$, by mapping the classical data into the state of a multi-qubit system.

Our first main result is a quantum algorithm with running time linear in $m$ that generates an approximately optimal solution for the structural SVM problem. By Theorem 1, this is equivalent to solving the original soft-margin $\ell_{1}$-SVM.

Here and in what follows, the tilde big-O notation hides polylogarithmic terms. In the Simulation section we show that, in practice, the running time of the algorithm can be significantly faster than the theoretical upper-bound. The solution $\hat{\alpha}$ is a $t_{\max}$-sparse vector of total dimension $2^{m}$. Once it has been found, a new data point $\mathbf{x}$ can be classified according to

where $y_{pred}$ is the predicted label of $\mathbf{x}$. This is a sum of $O(mt_{\max})$ inner products in feature space, which classical methods require time $O(mt_{\max})$ to evaluate in general. Our second result is a quantum algorithm for carrying out this classification with running time independent of $m$.

On the surface, the structural SVM problems OP 3 and OP 4 look more complicated to solve than the original SVM problems OP 1 and OP 2. However, it turns out that the solution $\alpha^{*}$ to OP 4 is highly sparse and, consequently, the structural SVM admits an efficient algorithm. Joachims’ original procedure is presented in Algorithm 1.

The main idea behind Algorithm 1 is to iteratively solve successively more constrained versions of problem OP 3. That is, a working set of indices $\mathcal{W}\subseteq\{0,1\}^{m}$ is maintained such that, at each iteration, the solution $(\mathbf{w},\xi)$ is only required to satisfy the constraints $\frac{1}{m}\sum_{i=1}^{m}c_{i}-\left\langle\mathbf{w},\Psi_{c}\right\rangle\leq\xi$ for $\mathbf{c}\in\mathcal{W}$. The inner for loop then finds a new index $\mathbf{c}^{*}$ which corresponds to the maximally violated constraint in OP 3, and this index is added to the working set. The algorithm proceeds until no constraint is violated by more than $\epsilon$. It can be shown that each iteration must improve the value of the dual objective by a constant amount, from which it follows that the algorithm terminates in a number of rounds independent of $m$.

In terms of time cost, each iteration $t$ of the algorithm involves solving the restricted optimization problem

which is done in practice by solving the corresponding dual problem, i.e. the same as OP 4 but with summations over $\mathbf{c}\in\mathcal{W}$ instead of over all $\mathbf{c}\in\{0,1\}$. This involves computing

• •

  $O(t^{2})$ matrix elements $J_{\mathbf{c}\mathbf{c}^{\prime}}=\left\langle\Psi_{\mathbf{c}},\Psi_{\mathbf{c}^{\prime}}\right\rangle$

• •

  $m$ inner products $\left\langle\mathbf{w},\mathbf{x}_{i}\right\rangle=\left\langle\sum_{\mathbf{c}}\alpha_{\mathbf{c}}\Psi_{\mathbf{c}},\mathbf{x}_{i}\right\rangle$

where $\alpha_{\mathbf{c}}$ is the solution to the dual of the optimization problem in the body of Algorithm 1. In the case of linear SVMs, $\Psi_{\mathbf{c}}=\frac{1}{m}\sum_{i=1}^{m}c_{i}y_{i}\mathbf{x}_{i}$ and $\left\langle\Psi_{\mathbf{c}},\Psi_{\mathbf{c}^{\prime}}\right\rangle$ can each be explicitly computed in time $O(m)$. The cost of computing matrix $J$ is $O(mt^{2})$, and subsequently solving the dual takes time polynomial in $t$. As Joachims showed that $t\leq\max\left\{\frac{2}{\epsilon},\frac{8C\max_{i}\left\|\mathbf{x}_{i}\right\|^{2}}{\epsilon^{2}}\right\}$, and since $\left\langle\mathbf{w},\mathbf{x}_{i}\right\rangle$ can be computed in time $O(d)$, the entire algorithm therefore has running time linear in $m$. Note that Joachims considered the special case of $s$-sparse data vectors, for which $\left\langle\mathbf{w},\mathbf{x}_{i}\right\rangle$ can be computed in time $O(s)$ rather than $O(d)$. In what follows we will not consider any sparsity restrictions.

For nonlinear SVMs, the feature maps $\Phi(\mathbf{x}_{i})$ may be of very large dimension, which precludes explicitly computing $\Psi_{\mathbf{c}}=\frac{1}{m}\sum_{i=1}^{m}c_{i}y_{i}\Phi(\mathbf{x}_{i})$ and directly evaluating $\left\langle\Psi_{\mathbf{c}},\Psi_{\mathbf{c}^{\prime}}\right\rangle$ as is done in SVM-perf. Instead, one must compute $J_{\mathbf{c}\mathbf{c}^{\prime}}=\left\langle\Psi_{\mathbf{c}},\Psi_{\mathbf{c}^{\prime}}\right\rangle=\frac{1}{m^{2}}\sum_{i,j=1}^{m}c_{i}c_{j}y_{i}y_{j}\left\langle\Phi(\mathbf{x}_{i}),\Phi(\mathbf{x}_{j})\right\rangle$ as a sum of $O(m^{2})$ inner products $\left\langle\Phi(\mathbf{x}_{i}),\Phi(\mathbf{x}_{j})\right\rangle$, which are then each evaluated using the kernel trick. This rules out the possibility of an $O(m)$ algorithm, at least using methods that rely on the kernel trick to evaluate each $\left\langle\Phi(\mathbf{x}_{i}),\Phi(\mathbf{x}_{j})\right\rangle$. Noting that $\mathbf{w}=\sum_{\mathbf{c}}\alpha_{\mathbf{c}}\Psi_{\mathbf{c}}$ , the inner products $\left\langle\mathbf{w},\Phi(\mathbf{x}_{i})\right\rangle$ are similarly expensive to compute directly classically if the dimension of the feature map is large.

We now show how quantum computing can be used to efficiently approximate the inner products $\left\langle\Psi_{\mathbf{c}},\Psi_{\mathbf{c}^{\prime}}\right\rangle$ and $\left\langle\sum_{\mathbf{c}}\alpha_{\mathbf{c}}\Psi_{\mathbf{c}},\Phi(\mathbf{x}_{i})\right\rangle$, where high dimensional $\Psi_{\mathbf{c}}$ can be implemented by a quantum circuit using only a number of qubits logarithmic in the dimension. We first assume that the data vectors $\mathbf{x}_{i}\in\mathbb{R}^{d}$ are encoded in the state of an $O(d)$-qubit register $\left|\mathbf{x}_{i}\right\rangle$ via some suitable encoding scheme, e.g. given an integer $k$, $\mathbf{x}_{i}$ could be encoded in $n=kd$ qubits by approximating each of the $d$ values of $\mathbf{x}_{i}$ by a length $k$ bit string which is then encoded in a computational basis state of $n$ qubits. Once encoded, a quantum feature map encodes this information in larger space in the following way:

Note that the states $\left|\Phi(\mathbf{x})\right\rangle$ are not necessarily orthogonal. Implementing such a quantum feature map could be done, for instance, through a controlled parameterized quantum circuit.

We also define the quantum state analogy of $\Psi_{\mathbf{c}}$ from Definition 1:

Let real vectors $x,y\in\mathbb{R}^{d}$ have corresponding normalized quantum states $\left|x\right\rangle=\frac{1}{\left\|x\right\|}\sum_{i=1}^{d}x_{i}\left|i\right\rangle$ and $\left|y\right\rangle=\frac{1}{\left\|y\right\|}\sum_{i=1}^{d}y_{i}\left|i\right\rangle$. The following result shows how the inner product $\left\langle x,y\right\rangle=\langle x|y\rangle\left\|x\right\|\left\|y\right\|$ can be estimated efficiently on a quantum computer.

Thus, if one can efficiently create quantum states $\left|\Psi_{\mathbf{c}}\right\rangle$ and estimate the norms $\left\|\Psi_{\mathbf{c}}\right\|$, then the corresponding $J_{\mathbf{c}\mathbf{c}^{\prime}}=\left\langle\Psi_{\mathbf{c}},\Psi_{\mathbf{c}^{\prime}}\right\rangle=\left\|\Psi_{\mathbf{c}}\right\|\left\|\Psi_{\mathbf{c}^{\prime}}\right\|\langle\Psi_{\mathbf{c}}|\Psi_{\mathbf{c}^{\prime}}\rangle$ can be approximated efficiently. In this section we show that this is possible with a quantum random access memory (qRAM), which is a device that allows classical data to be queried efficiently in superposition. That is, if $x\in\mathbb{R}^{d}$ is stored in qRAM, then a query to the qRAM implements the unitary $\sum_{j}\alpha_{j}\left|j\right\rangle\left|0\right\rangle\rightarrow\sum_{j}\alpha_{j}\left|j\right\rangle\left|x_{j}\right\rangle$. If the elements $x_{j}$ of $x$ arrive as a stream of entries $(j,x_{j})$ in some arbitrary order, then $x$ can be stored in a particular data structure [33] in time $\tilde{O}(d)$ and, once stored, $\left|x\right\rangle=\frac{1}{\left\|x\right\|}\sum_{j}x_{j}\left|j\right\rangle$ can be created in time polylogarithmic in $d$. Note that when we refer to real-valued data being stored in qRAM, it is implied that the information is stored as a binary representation of the data, so that it may be loaded into a qubit register.

A similar result applies to estimating inner products of the form $\sum_{\mathbf{c}\alpha_{c}}\left\langle\Psi_{\mathbf{c}},y_{i}\Phi(\mathbf{x}_{i})\right\rangle$.

Proofs of Theorems 4 and 5 are given in Appendix A.

As is standard in SVM theory, the solution $\hat{\alpha}$ from Algorithm 2 can be used to classify a new data point $\mathbf{x}$ according to

where $y_{pred}$ is the predicted label of $\mathbf{x}$. From Theorem 5, and noting that $\left|\mathcal{W}\right|\leq t_{\max}$, we obtain the following result:

Taking the sign of the output then completes the classification.

To test our algorithm we need to choose both a data set as well as a quantum feature map. The general question of what constitutes a good quantum feature map, especially for classifying classical data sets, is an open problem and beyond the scope of this investigation. However, if the data is generated from a quantum problem, then physical intuition may guide our choice of feature map. We therefore consider the following toy example which is nonetheless instructive. Let $H_{N}$ be the Hamiltonian of a generalized Ising Hamiltonian on $N$ spins

where $\vec{J},\vec{\Delta},\vec{\Gamma}$ are vectors of real parameters to be chosen, and $Z_{j},X_{j}$ are Pauli $Z$ and $X$ operators acting on the $j$-th qubit in the chain, respectively. We generate a data set by randomly selecting $m$ points $(\vec{J},\vec{\Delta},\vec{\Gamma})$ and labelling them according to whether the expectation value of the operator $M=\frac{1}{N}\left(\sum_{j}Z_{j}\right)^{2}$ with respect to the ground state of $H_{N}(\vec{J},\vec{\Delta},\vec{\Gamma})$ satisfies

for some cut-off value $\mu_{0}$, i.e. the points are labelled depending on whether the average total magnetism squared is above or below $\mu_{0}$. In our simulations we consider a special case of (3) where $J_{j}=J\cos\frac{k_{J}\pi(j-1)}{N}$, $\Delta_{j}=\Delta\sin\frac{k_{\Delta}\pi j}{N}$ and $\Gamma_{j}=\Gamma$, where $J,k_{J},\Delta,k_{\Delta},\Gamma$ are real. Examples of data sets $S_{N,m}$ corresponding to such a Hamiltonian, whose parameters we notate by

can be found in Fig 1.

For quantum feature map we choose

where $\left|\psi_{GS}\right\rangle$ is the ground state of (3) and, as it is a normalized state, has corresponding value of $R=1$. We compute such feature maps classically by explicitly diagonalizing $H_{N}$. In a real implementation of our algorithm on a quantum computer, such a feature map would be implemented by a controlled unitary for generating the (approximate) ground state of $H_{N}$, which could be done by a variety of methods e.g. by digitized adiabatic evolution or methods based on imaginary time evolution [34, 35], with running time $T_{\phi}$ dependent on the degree of accuracy required. The choice of (5) is motivated by noting that condition (4) is equivalent to determining the sign of $\left\langle W,\Psi\right\rangle$, where $W$ is a vector which depends only on $\mu_{0}$, and not on the choice of parameters in $H_{N}$ (see Appendix E). By construction, $W$ defines a separating hyperplane for the data, so the chosen quantum feature map separate the data in feature space. As the Hamiltonian is real, it has a set of real eigenvectors and hence $\left|\Psi\right\rangle$ can be defined to have real amplitudes, as required.

We first evaluate the performance of our algorithm on data sets $S_{N,m}$ for $N=6$ and increasing orders of $m$ from $10^{2}$ to $10^{5}$.

• •

  For each value of $m$, a data set $S_{6,m}(\mu_{0},J,k_{J},\Delta,k_{\Delta},\Gamma)$ was generated for points $(J,\Delta)$ sampled uniformly at random in the range $[-2,2]^{2}$.

• •

  The values of $\mu_{0},k_{j},k_{\Delta},\Gamma$ were fixed and chosen to give roughly balanced data, i.e. the ratio of $+1$ to $-1$ labels is no more than 70:30 in favor of either label.

• •

  Each set of $m$ data points was divided into training and test sets in the ratio 70:30, and training was performed according to Algorithm 2 with parameters $(C,\epsilon,\delta,t_{\max})=(10^{4},10^{-2},10^{-1},50)$.

• •

  These values of $C$ and $\epsilon$ were selected to give classification accuracy competitive with classical SVM algorithms utilizing standard Gaussian radial basis function (RBF) kernels, with hyperparameters trained using a subset of the training set of size $20\%$ used for hold-out validation. Note that the quantum feature maps do not have any similar tunable parameters, and a modification of (5), for instance to include a tunable weighting between the two parts of the superposition, could be introduced to further improve performance.

• •

  The quantum $IP_{\epsilon,\delta}$ inner product estimations in the algorithm were approximated by adding Gaussian random noise to the true inner product, such that the resulting inner product was within $\epsilon$ of the true value with probability at least $1-\delta$. Classically simulating quantum $IP_{\epsilon,\delta}$ inner product estimation with inner products distributed according to the actual quantum procedures underlying Theorems 4 and 5 was too computationally intensive in general to perform. However, these were tested on small data sets and quantum feature vectors, and found to behave very similarly to adding Gaussian random noise. This is consistent with the results of the numerical simulations in [32].

Note that the values of $C,\epsilon,\delta$ chosen correspond to $\max\left\{\frac{4}{\epsilon},\frac{16CR^{2}}{\epsilon^{2}}\right\}>10^{9}$. This is an upper-bound on the number of iterations $t_{\max}$ needed for the algorithm to converge to a good solution. However, we find empirically that $t_{\max}=50$ is sufficient for the algorithm to terminate with a good solution across the range of $m$ we consider.

The results are shown in Table 1. We find that (i) with these choices of $C,\epsilon,\delta,t_{\max}$ our algorithm has high classification accuracy, competitive with standard classical SVM algorithms utilizing RBF kernels with optimized hyperparameters. (ii) $\Psi_{\min}$ is of the order $10^{-2}$ in these cases, and does scale noticeably over the range of $m$ from $10^{2}$ to $10^{5}$. If $\Psi_{\min}$ were to decrease polynomially (or worse, exponentially) in $m$ then this would be a severe limitation of our algorithm. Fortunately this does not appear to be the case.

We further investigate the behaviour of $\Psi_{\min}$ by generating data sets $S_{N,m}$ for fixed $m=1000$ and $N$ ranging from $4$ to $8$. For each $N$, we generate $100$ random data sets $S_{N,m}(\mu_{0},J,k_{J},\Delta,k_{\Delta},\Gamma)$, where each data set consists of $1000$ points $(J,\Delta)$ sampled uniformly at random in the range $[-2,2]^{2}$, and random values of $\mu_{o},k_{J},k_{\Delta},\Gamma$ chosen to give roughly balanced data sets as before. Unlike before, we do not divide the data into training and test sets. Instead, we perform training on the entire data set, and record the value of $\Psi_{\min}$ in each instance. The results are given in Table 2 and show that across this range of $N$ (i) the average value $\bar{\Psi}_{\min}$ is of order $10^{-2}$ (ii) the spread around this average is fairly tight, and the minimum value of $\Psi_{\min}$ in any single instance is of order $10^{-3}$. These support the results of the first experiment, and indicate that the value of $\Psi_{\min}$ may not adversely affect the running time of the algorithm in practice.