Title: QUANTUM SEMI-SUPERVISED KERNEL LEARNING

Abstract: Quantum machine learning methods have the potential to facilitate learning using extremely large datasets. While the availability of data for training machine learning models is steadily increasing, oftentimes it is much easier to collect feature vectors that to obtain the corresponding labels. One of the approaches for addressing this issue is to use semi-supervised learning, which leverages not only the labeled samples, but also unlabeled feature vectors. Here, we present a quantum machine learning algorithm for training Semi-Supervised Kernel Support Vector Machines. The algorithm uses recent advances in quantum sample-based Hamiltonian simulation to extend the existing Quantum LS-SVM algorithm to handle the semi-supervised term in the loss, while maintaining the same quantum speedup as the Quantum LS-SVM.

1. INTRODUCTION:

Data sets used for training machine learning models are becoming increasingly large, leading to continued interest in fast methods for solving large-scale classification problems. One of the approaches being explored is training the predictive model using a quantum algorithm that has access to the training set stored in quantum-accessible memory. In parallel to research on efficient architectures for quantum memory (Blencowe, 2010) , work on quantum machine learning algorithms and on quantum learning theory is under way (see for example Refs. (Biamonte et al., 2017; Dunjko & Briegel, 2018; Schuld & Petruccione, 2018) and (Arunachalam & de Wolf, 2017) for review). An early example of this approach is Quantum LS-SVM (Rebentrost et al., 2014a) , which achieves exponential speedup compared to classical LS-SVM algorithm. Quantum LS-SVM uses quadratic least-squares loss and squared-L 2 regularizer, and the optimization problem can be solved using the seminal HHL (Harrow et al., 2009) algorithm for solving quantum linear systems of equations. While progress has been made in quantum algorithms for supervised learning, it has been recently advocated that the focus should shift to unsupervised and semi-supervised setting (Perdomo-Ortiz et al., 2018) .

In many domains, the most laborious part of assembling a training set is the collection of sample labels. Thus, in many scenarios, in addition to the labeled training set of size m we have access to many more feature vectors with missing labels. One way of utilizing these additional data points to improve the classification model is through semi-supervised learning. In semi-supervised learning, we are given m observations x 1 , ..., x m drawn from the marginal distribution p(x), where the l (l m) first data points come with labels y 1 , ..., y l drawn from conditional distribution p(y|x). Semi-supervised learning algorithms exploit the underlying distribution of the data to improve classification accuracy on unseen samples. In the approach considered here, the training samples are connected by a graph that captures their similarity.

Here, we introduce a quantum algorithm for semi-supervised training of a kernel support vector machine classification model. We start with the existing Quantum LS-SVM (Rebentrost et al., 2014a) , and use techniques from sample-based Hamiltonian simulation (Kimmel et al., 2017) to add a semisupervised term based on Laplacian SVM (Melacci & Belkin, 2011) . As is standard in quantum machine learning (Li et al., 2019) , the algorithm accesses training points and the adjacency matrix of the graph connecting samples via a quantum oracle. We show that, with respect to the oracle, the proposed algorithm achieves the same quantum speedup as LS-SVM, that is, adding the semisupervised term does not lead to increased computational complexity.

2. PRELIMINARIES:

2.1 SEMI-SUPERVISED LEAST-SQUARES KERNEL SUPPORT VECTOR MACHINES.

Consider a problem where we are aiming to find predictors h(x) : X → R that are functions from a RKHS defined by a kernel K. In Semi-Supervised LS-SVMs in RKHS, we are looking for a function h ∈ H that minimizes

min h∈H,b∈R γ 2 l i=1 (y i -(h(x i ) + b)) 2 + 1 2 ||h|| 2 H + 1 2 ||∇h|| 2 E ,

where γ is a user define constant allowing for adjusting the regularization strength. The last term captures the squared norm of the graph gradient on the graph G that contains all training samples as vertices, and expresses similarity between samples through, possibly edges

G u,v 1 2 ||∇h|| 2 E = 1 2 u∼v G u,v ( hu -hv ) 2 = hT L h,

where hu is the function value h(x i ) for the vertex u corresponding to training point x i , and L is the

combinatorial graph Laplacian matrix such that L[i, j] = D j -G i,j .

The Representer Theorem states that if H is RKHS defined by kernel K : X × X → R, then the solution minimizing the problem above is achieved for a function h that uses only the representers of the training points, that is, a function of the form h

(x) = m j=1 c j K xj (x) = m j=1 c j K(x j , x

). Thus, we can translate the problem from RKHS into a constrained quadratic optimization problem over finite, real vectors

min c,ξ,b γ 2 m i=1 ξ 2 i + 1 2 c T Kc + 1 2 c T KLKc s.t. 1 -y i   b + m j=1 c j K[i, j]   = ξ i

where l ≤ m is the number of training points with labels (these are grouped at the beginning of the training set), and h = Kc, since function h is defined using representers K xi . The semi-supervised term, the squared norm of the graph gradient of h, 1/2||∇h|| 2 E , penalizes large changes of function h over edges of graph G. In defining the kernel K and the Laplacian L and in the two regularization terms we use all m samples. On the other hand, in calculating the empirical quadratic loss we only use the first l samples.

The solution to the Semi-Supervised LS-SVMs is given by solving the following (m + 1) × (m + 1) system of linear equations

EQUATION

where y = (y 1 , ..., y m ) T , 1 = (1, ..., 1) T , 1 is identity matrix, K is kernel matrix, L is the graph Laplacian matrix, γ is a hyperparameter and α = (α 1 , ..., α m ) T is the vector of Lagrangian multipliers.

2.2 QUANTUM COMPUTING AND QUANTUM LS-SVM Quantum computers are devices which perform computing according to the laws of quantum mechanics, a mathematical framework for describing physical theories, in language of linear algebra.

Quantum Systems. Any isolated, closed quantum physical system can be fully described by a unit-norm vector in a complex Hilbert space appropriate for that system; in quantum computing, the space is always finite-dimensional, C d . In quantum mechanics and quantum computing, Dirac notation for linear algebra is commonly used. In Dirac notation, a vector x ∈ C d and its complex conjugate x T , which represents a functional C d → R, are denoted by |x (called ket) and x| (called bra), respectively. We call {|e i } d i=1 the computational basis, where |e i = (0, ..., 1, ...0) T with exactly one 1 entry in the i-th position. Any |v = (v 1 , ..., v d ) T can be written as |v

= d i=1 v i |e i ; coefficient v i ∈ C are called probability amplitudes. Any unit vector |x ∈ C d describes a d-level quantum state. Such a system is called a pure state. An inner product of |x 1 , |x 2 ∈ C d is written as x 1 |x 2 . A two-level quantum state |ψ = α|0 + β|1 , where |0 = (1, 0) T , |1 = (0, 1) T and α, β ∈ C, |α| 2 + |β| 2 ,

is called a quantum bit, or qubit for short. When both α and β are nonzero, we say |ψ is in a superposition of the computational basis |0 and |1 ; the two superposition states

|+ = 1 √ 2 (|0 + |1 ) and |-= 1 √ 2 (|0 -|1

) are very common in quantum computing.

A composite quantum state of two distinct quantum systems |x 1 ∈ C d1 and |x 2 ∈ C d2 is described as tensor product of quantum states |x 1 ⊗ |x 2 ∈ C d1 ⊗ C d2 . Thus, a state of an n-qubit system is a vector in the tensor product space 

C 2 ⊗n = C 2 ⊗ C 2 ⊗ ... ⊗ C 2 ,

+ U|u 2 = U (|u 1 + |u 2 )

, applying a unitary to a superposition of states has the same effect as applying it separately to element of the superposition. In quantum circuit model unitary transformations are referred to as quantum gates -for example, one of the most common gates, the single-qubit Hadamard gate, is a unitary operator represented in the computational basis by the matrix

EQUATION

Note that H|0 = |+ and H|1 = |-.

Quantum Input Model. Quantum computation typically starts from all qubits in |0 state. To perform computation, access to input data is needed. In quantum computing, input is typically given by a unitary operator that transforms the initial state into the desired input state for the computation -such unitaries are commonly referred to as oracles, and the computational complexity of quantum algorithms is typically measured with access to an oracle as the unit. For problems involving large amounts of input data, such as for quantum machine learning algorithms, an oracle that abstracts random access memory is often assumed. Quantum random access memory (qRAM) uses log N qubits to address any quantum superposition of N memory cell which may contains either quantum or classical information. For example, qRAM allows accessing classical data entries x j i in quantum superposition by a transformation

1 √ mp m i=1 p j=1 |i, j |0...0 qRAM ---→ 1 √ mp m i=1 p j=1 |i, j |x j i ,

where |x j i is a binary representation up to a given precision. Several approaches for creating quantum RAM are being considered (Giovannetti et al., 2008; Arunachalam et al., 2015; Biamonte et al., 2017) , but it is still an open challenge, and subtle differences in qRAM architecture may erase any gains in computational complexity of a quantum algorithm Aaronson (2015) .

Quantum Linear Systems of Equations. Given an input matrix A ∈ C n×n and a vector b ∈ C n , the goal of linear system of equations problem is finding x ∈ C n such that Ax = b. When A is Hermitian and full rank, the unique solution is x = A -1 b. If A is not a full rank matrix then A -1 is replaced by the Moore-Penrose pseudo-inverse. HHL algorithm introduced an analogous problem in quantum setting: assuming an efficient algorithm for preparing b as a quantum state b = n i=1 b i |i using log n + 1 qubits, the algorithm applies quantum subroutines of phase estimation, controlled rotation, and inverse of phase estimation to obtain the state

EQUATION

)

Intuitively and at the risk of over-simplifying, HHL algorithm works as follows: if A has spectral decomposition A = n i=1 λ i v i v T i (where λ i and v i are corresponding eigenvalues and eigenstates of A), then A -1 maps λ i v i → 1 λ i v i . The vector b also can be written as the linear combination of the A's eigenvectors v i as b = n i=1 β i v i (we are not required to compute

β i ). Then A -1 b = n i=1 β i 1 λ i v i .

In general A and A -1 are not unitary (unless all A's eigenvalues have unit magnitude), therefore we are not able to apply A -1 directly on |b . However, since U = e iA = n i=1 e iλi v i v T i is unitary and has the same eigenvectors as A and A -1 , one can implement U and powers of U on a quantum computer by Hamiltonian simulation techniques; clearly for any expected speed-up, one need to enact e iA efficiently. The HHL algorithm uses the phase estimation subroutine to estimate an approximation of λ i up to a small error. The Next step computes a conditional rotation on the approximated value of λ i and an auxiliary qubit |0 and outputs

1 λ i |0 + 1 - 1 λ 2 i |1

. The last step involves the inverse of phase estimation and quantum measurement for getting rid of garbage qubits and outputs our desired state

|x = A -1 |b = n i=1 β i 1 λ i v i .

Density Operators. Density operator formalism is an alternative formulation for quantum mechanics that allows probabilistic mixtures of pure states, more generally referred to as mixed states.

A mixed state that describes an ensemble {p i , |ψ i } is written as

EQUATION

where k i=1 p i = 1 forms a probability distribution and ρ is called density operator, which in a finite-dimensional system, in computational basis, is a semi-definite positive matrix with T r(ρ) = 1.

A unitary operator U maps a quantum state expressed as a density operator ρ to UρU † , where U † is the complex conjugate of the operator U.

Partial Trace of Composite Quantum System. Consider a two-part quantum system in a state described by tensor product of two density operators ρ ⊗ σ. A partial trace, tracing out the second part of the quantum system, is defined as the linear operator that leaves the first part of the system in a state Tr 2 (ρ ⊗ σ) = ρ tr (σ), where Tr (σ) is the trace of the matrix σ.

To obtain Kernel matrix K as a density matrix, quantum LS-SVM (Rebentrost et al., 2014b) relies on partial trace, and on a quantum oracle that can convert, in superposition, each data point {x i } m i=1 , x i ∈ R p to a quantum state |x i = 1 xi p t=1 (x i ) t |t , where (x i ) t refers to the tth feature value in data point x i and assuming the oracle is given x i and y i . Vector of the labels is given in the same fashion as |y = 1 y m i=1 y i |i . For preparation the normalized Kernel matrix K = 1 tr(K) K where K = X T X, we need to prepare a quantum state combining all data points in quantum

EQUATION

The approach used above to construct density matrix corresponding to linear kernel matrix can be extended to polynomial kernels (Rebentrost et al., 2014b) .

LMR Technique for Density Operator Exponentiation. In HHL-based quantum machine learning algorithms , including in the method proposed here, matrix A for the Hamiltonian simulation within the HHL algorithm is based on data. For example, A can contain the kernel matrix K captured in the quantum system as a density matrix. Then, one need to be able to efficiently compute e -iK∆t , where K is scaled by the trace of kernel matrix. Since K is not sparse, a strategy similar to (Lloyd et al., 2014) is adapted for the exponentiation of a non-sparse density matrix:

EQUATION

where S = i,j |i j| ⊗ |j i| is the swap operator and the facts Tr 1 {S(K ⊗ σ)} = Kσ and Tr 1 {(K ⊗ σ)S} = σK are used. The equation ( 6) summarizes the LMR technique: approximating e -iK∆t σe iK∆t up to error O(∆t 2 ) is equivalent to simulating a swap operator S, applying it to the state K ⊗ σ and discarding the first system by taking partial trace operation. Since the swap operator is sparse, its simulation is efficient. Therefore the LMR trick provides an efficient way to approximate exponentiation of a non-sparse density matrix.

Quantum LS-SVM. Quantum LS-SVM (Rebentrost et al., 2014b) uses partial trace to construct density operator corresponding to the kernel matrix K. Once the kernel matrix K becomes available as a density operator, the quantum LS-SVM proceeds by applying the HHL algorithm for solving the system of linear equations associated with LS-LSVM, using the LMR technique for performing the density operator exponentiation e -iK∆t where the density matrix K encodes the kernel matrix.

3 QUANTUM SEMI-SUPERVISED LEAST SQUARE SVM.

Semi-Supervised Least Square SVM involves solving the following system of linear equations

b α = 0 1 T 1 K + KLK + γ -1 1 -1 0 y = A -1 0 y (7)

In quantum setting the task is to generate |b, α = Â-1 |0, y , where the normalized Â = A T r(A) .

The linear system differs from the one in LS-SVM in that instead of K, we have K + KLK. While this difference is of little significance for classical solvers, in quantum systems we cannot just multiply and then add the matrices and then apply quantum LS-SVM -we are limited by the unitary nature of quantum transformations.

In order to obtain the solution to the quantum Semi-Supervised Least Square SVM, we will use the following steps. First, we will read in the graph information to obtain scaled graph Laplacian matrix as a density operator. Next, we will use polynomial Hermitian exponentiation for computing the matrix inverse (K + KLK) -1 .

3.1. QUANTUM INPUT MODEL FOR THE GRAPH LAPLACIAN:

In the semi-supervised model used here, we assume that we have information on the similarity of the training samples, in a form of graph G that uses n edges to connect similar training samples, represented as m vertices. We assume that for each sample, G contains its k most similar other samples, that is, the degree of each vertex is d. To have the graph available as a quantum density operator, we observe that the graph Laplacian L is the Gram matrix of the rows of the m × n graph incidence matrix G I , L = G I G T I . We assume oracle access to the graph adjacency list, allowing us to construct, in superposition, states corresponding to rows of the graph incidence matrix

G I |v i = 1 √ d n t=1 G I [i, t]|t .

That is, state |v i has probability amplitude 1 √ d for each edge |t incident with vertex i, and null probability amplitude for all other edges. In superposition, we prepare a quantum state combining rows of the incidence matrix for all vertices, to obtain

|G I = 1 √ md m i=1 |i ⊗ |v i

The graph Laplacian matrix L, composed of inner products of the rows of G I , is obtained by discarding the second part of the system,

EQUATION

3.2. POLYNOMIAL HERMITIAN EXPONENTIATION FOR SEMI SUPERVISED LEARNING:

For computing the matrix inverse (K + KLK) -1 on a quantum computer that runs our quantum machine algorithm and HHL algorithm as a subroutine, we need to efficiently compute e -i(K+KLK)∆t σe i(K+KLK)∆t . For this purpose we adapt the generalized LMR technique for simulating Hermitian polynomials proposed in (Kimmel et al., 2017) to the specific case of e -i(K+KLK)∆t σe i(K+KLK)∆t . Simulation of e -iK∆t follows from the original LMR algorithm, and therefore we focus here only on simulation e -iKLK∆t . The final dynamics (K + KLK) -1 can be obtained by sampling from the two separate output states for e -iKLK∆t and e -iK∆t .

Simulating e iKLK∆t e iKLK∆t e iKLK∆t . Let D(H) denote the space of density operators associated with state space H. Let K † , K, L ∈ D(H) be the density operators associated with the kernel matrix and the Laplacian, respectively. We will need two separate systems with the kernel matrix K, to distinguish between them we will denote the first as K † and the second as K; since K is real and symmetric, these are indeed equal. The kernel and Laplacian matrices K † , K, L are not sparse therefore we adapt the generalized LMR technique for simulating Hermitian polynomials for our specific case B = K † LK.

For adapting the generalized LMR technique to our problem we need to generate a quantum state

EQUATION

where

B = ρ -ρ = 1 2 K † LK + 1 2 KLK † = KLK and S := |0 0| ⊗ S + |1 1| ⊗ (-S

) is a controlled partial swap in the forward (+S) and backward direction (-S) in time, and

e -iS ∆ = |0 0| ⊗ e -iS∆ + |1 1| ⊗ e iS∆ .

Therefore with one copy of ρ , we obtain the simulation of e -iB∆ up to error O(∆ 2 ). If we choose the time slice ∆ = δ/t and repeating the above procedure for t 2 /δ times, we are able to simulate e -iBt up to error O(δ) using n = O(t 2 /δ) copies of ρ . The analysis of the steps preformed by the circuit depicted in Fig. 1 is as follows. Let P be the cyclic permutation of three copies of H A that operates as P |j 1 , j 2 , j 3 = |j 3 , j 1 , j 2 . In operator form it can be written as 

Generating ρ = |0 0| ⊗ ρ + |1 1| ⊗ ρ ρ = |0 0| ⊗ ρ + |1 1| ⊗ ρ ρ = |0 0| ⊗ ρ + |1 1| ⊗ ρ .

EQUATION

3.3. QUANTUM SEMI-SUPERVISED LS-SVM ALGORITHM AND ITS COMPLEXITY:

The quantum LS-SVM in (Rebentrost et al., 2014b) offers exponential speedup O(log mp) over the classical time complexity for solving SVM as a quadratic problem, which requires time O(log( -1 )poly(p, m)), where is the desired error. The exponential speedup in p occurs as the result of fast quantum computing of kernel matrix, and relies on the existence of efficient oracle access to data. The speedup on m is due to applying quantum matrix inversion for solving LS-SVM, which is inherently due to fast algorithm for exponentiation of a resulting non-sparse matrix. Our algorithm introduces two additional steps: preparing the Laplacian density matrix, and Hamiltonian simulation for KLK. The first step involves oracle access to a sparse graph adjacency list representation, which is at least as efficient as the oracle access to non-sparse data points. The Hamiltonian simulation involves simulating a sparse conditional partial swap operator, which results an efficient strategy for applying e -iKLK∆t in time Õ(log(m)∆t), where the notation Õ hides more slowly growing factors in the simulation (Berry et al., 2007) .

3.4. COMPARISON WITH ALTERNATIVE APPROACHES:

Considerable effort has been devoted into designing fast classical algorithms for training SVMs. The decomposition-based methods such as SMO (Platt, 1998) are able to efficiently manage problems with large number of features p, but their computational complexities are super-linear in m. Other training strategies (Suykens & Vandewalle, 1999; Fung & Mangasarian, 2005; Keerthi & DeCoste, 2005) are linear in m but scale quadratically in p in the worst case. The Pegasos algorithm (Shalev-Shwartz et al., 2011) for non-linear kernel improves the complexity to Õ (m/(λ )), where λ, and are the regularization parameter of SVM and the error of the solution, respectively.

Beyond the classical realm, three quantum algorithms for training linear models have been proposed, the quantum LS-SVM that involves L 2 regularizer (Rebentrost et al., 2014a) , a recently proposed Quantum Sparse SVM which is limited to a linear kernel (Arodz & Saeedi, 2019) , and a quantum training algorithm that solves a maximin problem resulting from a maximum -not average -loss over the training set (Li et al., 2019) .

:

n -1 i=0 α i |i , where i is expressed using its binary representation; for example for n = 4, we have|2 = |0010 = |0 ⊗ |0 ⊗ |1 ⊗ |0 .

Figures and Tables:

Figure 1 shows the quantum circuit for creating ρ = |0 0| ⊗ ρ + |1 1| ⊗ ρ such that T r(ρ + ρ ) = 1 and B = ρ -ρ = KLK.

Figure 1: Quantum circuit for creating ρ = |0 0| ⊗ ρ + |1 1| ⊗ ρ . The circuit is to be read from left-to-right. Each wire at left shows its corresponding input state. The vertical rectangle denotes the cyclic permutation operator P on K, L, K defined in (10). H is the Hadamard gate, and the waste bins show partial trace. The measurement on the first quantum state is in computational basis.

10) Algorithm 1 Quantum Semi-Supervised LS-SVM Input: The datapoint set {x 1 , ...x l , ...x m } with the first l data points labeled and the rest unlabeled, y = (y 1 , ..., y l ) and the graph G Output: The classifier |α, b = A -1 |y 1: Quantum data preparation. Encode classical data points into quantum data points using quantum oracles O x : {x 1 , ...x l , ...x m } → |X = x i |x i and O x : y → |y . 2: Quantum Laplacian preparation. Prepare quantum density matrix using oracle access to G. 3: Matrix inversion. Compute the matrix inversion |α, b = A -1 |y via HHL algorithm. A quantum circuit for the HHL algorithm has three main steps: 4: Phase estimation, including efficient Hamiltonian simulation involving KLK (Section 3Classification. Based on Swap test algorithm, same as in Quantum LS-SVM.

and is written asTransforming and Measuring Quantum States. Quantum operations manipulate quantum states in order to obtain some desired final state. Two types of manipulation of a quantum system are allowed by laws of physics: unitary operators and measurements. Quantum measurement, if done in the computational basis, stochastically transforms the state of the system into one of the computa-

Table content: <table/>

Peer review (iclr2020):
This paper developes a quantum algorithm for kernel-based support vector machine working in a semi-supervised learning setting. The motivation is to utilise the significant advantage of quantum computation to train machine learning models on large-scale datasets efficiently. This paper reviews the existing work on using quantum computing for least-squares svm (via solving quantum linear systems of equations) and then extends it to deal with kernel svm in a semi-supervised setting. Strengths: This is an interesting emerging research topic that has its significance. Also, this paper prodives a nice tutorial on the key ideas of quantum machine learning and provides detailed derivations and analysis on the proposed algorithm. Weaknesses: The novelty of this work seems to be incremental. It largely extends the existing algorithms such as HHL and LMR. Minor issues: 1. Can any experimental study or applications be demonstrated, or only theoretical comptuational complexity can be compared? 2. In Section 1.1, L is defined as L = G_I G^T_I. In this definition, whether and how the edge weights are considered? Please clarify.