Supervised learning with quantum enhanced feature spaces

References

Supplementary Information: Supervised learning with quantum enhanced feature spaces

Consider a classification task on a set $C=\{1,2\ldots c\}$ of $c$ classes (labels) in a supervised learning scenario. In such settings, we are given a training set $T$ and a test set $S$, both are assumed to be labeled by a map $m:T\cup S\rightarrow C$ unknown to the programmer. Both sets $S$ and $T$ are provided to the programmer, but the programmer only receives the labels of the training set. So, formally, the programmer has only access to a restriction $m_{|T}$ of the indexing map $m$:

Description of the Algorithm

We consider two different learning schemes. The first is referred to as “Quantum variational classification”, the second is referred to “Quantum kernel estimation”. Both schemes construct a separating hyperplane in the state space of $n$ qubits. The classical data is mapped to this space with $\mbox{dim}=4^{n}$ using a unitary circuit family starting from the reference state $|\,0\rangle\langle 0\,|^{n}$.

For our first classification approach we design a variational algorithm which exploits the large dimensional Hilbert space of our quantum processor to find an optimal cutting hyperplane in a similar vein as Support Vector Machines (SVM) do. The algorithm consists of two main parts: a training stage and a classification stage. For the training stage, a set of labeled data points are provided, on which the algorithm is performed. For the classification stage, we take a different set of data points and run the optimized classifying circuit on them without any label input. Then we compare the label of each data point to the output of the classifier to obtain a success ratio for the data set. For both the training and the classification stages, the quantum circuit that implements the algorithm comprises three main parts: the encoding of the feature map, the variational optimization and the measurement, Fig S1. The training phase consists of these steps.

The classification can be applied when the training phase is complete. The optimal parameters are used to decide the correct label for new input data. Again, the same circuit is applied as in Fig S1, however, this time the parameters are fixed and and the outcomes are combined to determine the label which is reported as output of the classifier.

Quantum kernel estimation

For the second classification protocol, we restrict ourselves to the binary label case, with $C=\{+1,-1\}$. In this protocol we only use the quantum computer to estimate the $|T|\times|T|$ kernel matrix $K(\vec{x}_{i},\vec{x}_{j})=|\left\langle\,\Phi(\vec{x}_{i})\,|\,\Phi(\vec{x}_{j})\,\right\rangle|^{2}$. For all pairs of points $\vec{x}_{i},\vec{x}_{j}\in T$ in the the training data, we sample the overlap to obtain the matrix entry in the kernel. This output probability can be estimated from the circuit depicted in Fig. S5.b. by sampling the output distribution with $R$ shots and only taking the $0^{n}$ count. After the kernel matrix for the full training data has been constructed we use the conventional (classical) support vector machine classifier. The optimal hyperplane is constructed by solving the dual problem $L_{D}$ in eqn. (6), which is completely specified after we have been given the labels $y_{i}$ and have estimated the kernel $K(\vec{x}_{i},\vec{x}_{j})$. The solution of the optimization problem is given in terms of the support vectors $N_{S}$ for which $\alpha_{i}>0$.

The Relationship of variational quantum classifiers to support vector machines

The references vapnik2013nature ; burges1998tutorial provide a detailed introduction to the construction of support vector machines for pattern recognition. Support vector machines are an important tool to construct classifiers for tasks in supervised learning. We will show that the variational circuit classifier bears many similarities to a classical non-linear support vector machine.

The task in constructing a linear support vector machine (SVM) in this scenario is the following. One is looking for a hyperplane that separates the data, with the largest possible distance between the two separated sets. The perpendicular distance between the plane and two points with different labels is called a margin and such points are referred to as ‘support vectors’. This means that we want to maximize the margin by minimizing $||{\bf w}||$, or equivalently $||{\bf w}||^{2}$ subject to the constraints as given in eqn. (2), for all data points in the training set $T$. The corresponding cost function can be written as:

where $\alpha_{i}\geq 0$ are Lagrange multipliers chosen to ensure the constraints are satisfied.

These slack variables are then used to modify the objective function by $1/2\|{\bf w}\|^{2}\rightarrow 1/2\|{\bf w}\|^{2}+C(\sum_{i}\xi_{i})^{r}+\sum_{i}\mu_{i}\xi$. When we choose $r\geq 1$ the optimization problem remains convex and a dual can be constructed. In particular, for $r=1$, neither the $\xi_{i}$ or their Lagrange multipliers $\mu_{i}$ appear in the dual Lagrangian.

It is very helpful to consider the dual of the original primal problem $L_{P}$ in eqn. (4). The primal problem is a convex, quadratic programming problem, for which the Wolfe dual cost function $L_{D}$ for the Lagrange multipliers can be readily derived by variation with respect to ${\bf w}$ and $b$. The dual optimization problem is

The variables of the primal are given in terms of the dual variables by

and the bias $b$ can be computed from the Karush-Kuhn-Tucker (KKT) conditions when the corresponding Lagrange multiplier does not vanish. The optimal variables satisfy the KKT conditions and play an important role in the understanding of the SVM. They are given for primal as

Note that the condition eqn. (12) ensures that either the optimal $\alpha_{i}=0$ or the corresponding constraint eqn. (12) is tight. This is a property referred to as complementary slackness, and indicates that only the vectors for which the constraint is tight give rise to non-zero $\alpha_{i}>0$. These vectors are referred to as the support vectors and we will write $N_{S}$ for their index set. The classifier in the dual picture is given by substituting ${\bf w}$ from eqn. (7) and $b$ into the classifier eqn. (3). The bias $b$ is obtained for any $i\in N_{S}$ from the equality in eqn. (12).

The method can be generalized to the case when the decision function does depend non-linearly on the data by using a trick from boser1992training and introducing a high-dimensional, non-linear feature map. The data is mapped via

It is important to note that it is in fact not necessary to construct the mapped data $\Phi(\vec{x}_{i})$ in $\mathcal{H}$ explicitly. Both the training data, as well as the new data to be classified enters only through inner products, in both the optimization problem for training, c.f. eqn. (6), as well as in the classifier, eqn. (3). Hence, we can construct the SVM for arbitrarily high dimensional feature maps $(\Phi,\mathcal{H})$, if we can efficiently evaluate the inner products $\Phi(\vec{x}_{i})\circ\Phi(\vec{x}_{j})$ and $\Phi(\vec{x}_{i})\circ\Phi(\vec{s})$, for $\vec{x}_{i}\in T$ and $\vec{s}\in S$. In particular, if we can find a kernel $K(\vec{x},\vec{y})=\Phi(\vec{x})\circ\Phi(\vec{y})$ that satisfies Mercer’s condition (which ensures that the kernel is positive semi-definite and can be interpreted as matrix of inner products) boser1992training ; vapnik2013nature , we can construct a classifier by setting

Variational circuit classifiers:

where we have defined the diagonal operator

In classification tasks we assign c.f. eqn. (34), the label with the highest empirical weight of the distribution $p_{y}$. We ask whether the outcome $+1$ is more likely than $-1$, or vice versa. That is, we ask, whether $p_{+1}>p_{-1}-b$ or whether the converse is true. This of course depends on the sign of the expectation value $\langle\Phi(\vec{x})\,|W^{\dagger}(\theta,\varphi)\;{\bf f}\;W(\theta,\varphi)|\,\Phi(\vec{x})\rangle$ for the data point $\vec{x}$.

To understand how this relates to the SVM in greater detail, we need to choose an orthonormal operator basis, such as for example the Pauli group

This means that both the measurement operator $W^{\dagger}(\theta)\,{\bf f}\,W(\theta)$ in the $W$-rotated frame as well as the state $|\,\Phi(\vec{x})\rangle\langle\Phi(\vec{x})\,|$ can be expanded in terms of the operator basis with only real coefficients as

This expression is identical to the conventional SVM classifier, c.f. eqn. (3), after the feature map has been applied. However, in the experiment we only have access to the probabilities $p_{y}$ through estimation. Furthermore, the $w_{\alpha}(\theta)$ are constrained to stem from the observable ${\bf f}$ measured in the rotated frame.

This means, that the correct feature space, where a linearly separating hyperplane is constructed is in fact the quantum state space of density matrices, and not the Hilbert space $\mathcal{H}_{2}^{n}$ itself. This is reasonable, since the physical states in $\mathcal{H}_{2}^{n}$ are only defined up to a global phase $|\,\psi\rangle\sim e^{i\eta}|\,\psi\rangle$. The equivalence of states up to a global phase would make it impossible to find a separating hyperplane, since both $|\,\psi\rangle$ and $-|\,\psi\rangle$ give rise to the same physical state but can lie on either side of a separating plane.

Encoding of the data using a suitable feature map

The action of the map can be understood by a unitary circuit family denoted by ${\cal U}_{\Phi(\vec{x})}$ that is applied to some reference state, e.g. $|\,0\rangle^{n}$. The resulting state is given by $|\,\Phi(\vec{x})\rangle={\cal U}_{\Phi(\vec{x})}|\,0\rangle^{n}$. The state in the feature space should depend non-linearly on the data. Let us discuss proposals for possible feature maps

There are many choices for the feature map $\Phi$. Let us first discuss what would happen if we were to choose a feature map that corresponds to a product input state. We assume a feature map, comprised of single qubit rotations $U(\varphi)\in\text{SU}(2)$ on every qubit on the quantum circuit. The angles for every qubit can be chosen as a non-linear function $\varphi:\vec{x}\rightarrow(0,2\pi]^{2}\times[0,\pi]$ into the space of Euler angles for the individual qubits, so that the full feature map can be implemented as:

One example for such an implementation is the unitary implementation of the feature map used in the context of the classical classifiers by Stoudenmire and Schwab stoudenmire2016supervised based on tensor networks. There each qubit encodes a single component $x_{i}$ of $\vec{x}\in^{n}$ so that $n$ qubits are used. The resulting state that that is prepared is then

Non-trivial feature map with entanglement

There are many choices of feature maps, that do not suffer from the malaise of the aforementioned product state feature maps. To obtain an quantum advantage we would like these maps to give rise to a kernel $K(\vec{x},\vec{y})=|\left\langle\,\Phi(\vec{x})\,|\,\Phi(\vec{y})\,\right\rangle|^{2}$ that is computationally hard to estimate up to an additive polynomially small error by classical means. Otherwise the map is immediately amenable to classical analysis and we are guaranteed to have lost any conceivable quantum advantage.

Let us therefore turn to a family of feature maps, c.f. Fig S2 for which we conjecture that it is hard to estimate the overlap $|\left\langle\,\Phi(\vec{x})\,|\,\Phi(\vec{y})\,\right\rangle|^{2}$ on a classical computer. We define the family of feature map circuit as follows

which leaves $|V|+|E|$, real parameters to encode the data. In particular, we know that we have at least $|V|=n$ real numbers to encode the data. Furthermore, depending on the connectivity of the interactions, we have $|E|\leq n(n-1)/2$ further parameters that can be used to encode more data or nonlinear relations of the initial data points.

This feature map encodes both the actual function $\Phi_{\vec{x}}(z)$ of the diagonal phases, as well as the corresponding Fourier-Walsh transform $\hat{\Phi}_{\vec{x}}(p)$ at $z,p\in\{0,1\}^{n}$

Since we have that the variance of the random variable is bounded by $1$ since $|\Phi_{\vec{x}}(z)|^{2}=1$ , we get an additive error that scales as $\mathcal{O}(\epsilon)$. This means for a single layer, the kernel can be estimated classically.

Quantum variational classification

Following the structure of the feature map circuit, we construct the classifier part of the variational algorithm by appending layers of single-qubit unitaries and entangling gates Fig. S3.a. Each subsequent layer, or depth, contains an additional set of entanglers across all the qubits used for the algorithm. We use a coherently controllable quantum mechanical system, such as for example the superconducting chip with $n$ transmon qubits to prepare a short depth quantum circuit $W(\vec{\theta})$. In the experiment here, comprising $n=2$ qubits, one controlled-phase gate is added per depth. The single-qubit unitaries used in the classifier are limited to $Y$ and $Z$ rotations to simplify the number of parameters to be handled by the classical optimizer. Our use of controlled-phase, rather than CNOT, gates for the entanglers is justified by our aim at increased generality in our software. Using controlled-phase gates does not require to particularize this part of the algorithm for different systems topologies. A specific entangling map for a given device can then be used by our compiler to translate each controlled-phase gate into the CNOTs available in our system.

The general circuit is comprised of the following sequence of single qubit and multi-qubit gates:

We apply a circuit of $l$ repeated entanglers as depicted in Fig S3.b and interleave them with layers comprised of local single qubit rotations:

This short-depth circuit can generate any unitary if sufficiently many layers $d$ are applied. The circuit Fig. S3.a can be understood as a bang-bang controlled evolution of an Ising model Hamiltonian $H_{0}=\sum_{(ij)\in E}J_{ij}Z_{i}Z_{j}$,c.f. Fig. S3.b, interspersed with single qubit control pulses in $SU(2)$ on every qubit. It is known that this set of drift steps together with all single control pulses are universal, so we have that any state can be prepared this way with sufficient circuit depth d2007introduction . For a general unitary gate sequence the entangling unitary has to be effectively generated from cross resonance gates, by applying single qubit control pulses.

Choosing the cost-function for the circuit optimization

The central goal is to find an optimal classifying circuit $W(\vec{\theta})$, c..f. eqn. (31) that separates the data sets with different labels. Since we can re-run the same classifying multiple times ($R$ shots), we may consider a ‘winner takes all’ scenario, where we assign the label according to the outcome with the largest probability. We choose a cost function, so that the optimization procedure minimizes the probability of assigning the wrong label after having constructed the distribution after $R$ shots.

There are multiple ways of performing a multi-label classification. We only need to modify the final measurement $M$, to correspond to multiple partitions. This can be achieved by multiple strategies. For example one could choose to measure again in the computational basis, i.e. the basis in which Pauli $Z$ are diagonal and then constructing classical labels form the measured samples, such as a labeling the outcome $z\in\{0,1\}^{n}$ according to a function $f:\{0,1\}^{n}\rightarrow\{1,\ldots,c\}$. The resulting $\{M_{y}\}_{y=1,\ldots,c}$ is therefore diagonal in the computational basis. Alternatively one could construct a commuting measurement akin to the syndrome check measurement for quantum stabilizers. For this approach we choose a set $\{g_{i}\}_{i=1\ldots\lceil{\log_{2}(c)}\rceil}$ of Pauli matrices $g_{i}\in\mathcal{P}_{N}$ that are commuting $[g_{i},g_{j}]=0$. The resulting measurement that would need to be performed is similar to that of an error correcting scheme. The measurement operators are given by $M_{y}=2^{-1}\left(1-\prod_{i=1}^{\lceil{\log_{2}(c)}\rceil}g_{i}^{y^{i}}\right)$. Here $y^{i}$ denotes the $i$’th bit in the binary representation of $y$. In either case, the decision rule that assigns the labels can be written as

This corresponds to taking $R$ shots in order to estimate the largest outcome probability from the outcome statistics of the measurement $M_{y}$ for $y=1,\ldots,c$. Labelling $T_{c}$ the subset of samples $T$ labelled with $c$, the overall expected misclassification rate is given by:

Binary label classification

Assume the programmer classifies into labels $y\in\{-1,1\}$ by taking $R$ shots for a single datapoint. She obtains an empirical estimates of probability of the datum being labeled by a label $y$

After $R=r_{y}+r_{-y}$ shots and a prior bias $b$, she misclassifies into a label $y$ if

The probability of her misclassifying a $y$ sample according to the argmax rule is hence estimated by

Assuming large $R$, computing this exactly may be difficult. Setting $Rp_{y}=a,Rp_{y}(1-p_{y})=\beta^{2}$ and $\lceil\left(\frac{1+yb}{2}\right)R\rceil=\gamma$, we can approximate the binomial CDF as an error function:

See Fig. S4. The error function can be consequently approximated with a sigmoid

as an estimate for misclassifying a sample $s$. The cost function to optimize is then given by using this in eqn. (35).

For multiple labels, one tries to optimize

We consider the case of three labels. For $R$ samples with frequencies $\{n_{0},n_{1},n_{2}\}$, drawn independently from the output probability distribution, the probability of misclassifying a sample $s\in T_{0}$ by argmax is given by

where the last inequality is derived as follows

Hence setting $\gamma=\frac{N+|n_{1}-n_{2}|}{3}$, it follows that

This however still depends on $n_{1},n_{2}$, which can’t be simply eliminated. Additionally, for a general $k$-label case, there is no simple analytic solution for $\gamma$. For this reason, we therefore try to estimate the above probability by simply taking $\gamma=\max_{c^{\prime}}\left(\{n_{c^{\prime}}\}_{c^{\prime}/c}\right)$. So for $k$-label case, the cost function terms are given by

Quantum kernel estimation

For the second classification method we only use the quantum computer to estimate the kernel $K(\vec{x}_{i},\vec{x}_{j})=|\left\langle\,\Phi(\vec{x}_{i})\,|\,\Phi(\vec{x}_{j})\,\right\rangle|^{2}$ for all the labeled training data $\vec{x}_{j}\in T$. Then we use the classical optimization problem as outlined again in eqn. (6) to obtain the optimal Lagrange multipliers $\alpha_{i}$ and support vectors $N_{S}$ can be obtained. From this the classifier can be constructed, c.f .eqn. (14). To apply the classifier to a new datum $\vec{s}\in S$ the kernel $K(\vec{x}_{i},\vec{s})$ between $\vec{s}$ and the support vectors in $i\in N_{S}$ has to be estimated. We discuss two methods to estimate this overlap for our setting.

The usual method of estimating the fidelity between two states is by using the swap test buhrman2001quantum . This circuit, however, is not a short depth circuit on a quantum computing architecture with geometrically local gates. It would require a sequence of controlled SWAP, also known as Fredkin gates, all conditioned on the state of the same ancilla qubit. A very nice protocol was recently developed in cincio2018learning . The authors have learned multiple ways of optimizing the conventional swap test. If only the value of the fidelity is needed, as is the case for our algorithm, the authors propose, c.f. cincio2018learning section III.C, a circuit that is constant depth if pairs of CNOT gates can be executed in parallel. This circuit, c.f. Fig S5.a, evaluates the expectation value $\langle\psi\,|\langle\phi\,|\textsf{SWAP}|\,\psi\rangle|\,\phi\rangle$ directly. The action of the algorithm can be understood as follows:

The SWAP gate is both a unitary gate and a hermitian observable $\textsf{SWAP}^{\dagger}=\textsf{SWAP}$ with eigenvalues $\pm 1$. The expectation value on two product states as we said given by $\langle\psi\,|\langle\phi\,|\textsf{SWAP}|\,\psi\rangle|\,\phi\rangle=|\left\langle\,\phi\,|\,\psi\,\right\rangle|^{2}$. Now the gate can be decomposed in to a product of two qubit swap gates $\textsf{SWAP}=\prod_{k=1}^{n}\textsf{SWAP}_{s_{k}t_{k}}$ that act all in parallel. To evaluate the expectation value one has to diagonalize the full gate. This is achieved by diagonalizing the individual two-qubit swap gates by observing that $\textsf{SWAP}_{ij}=\textsf{CNOT}_{i\rightarrow j}\textsf{CNOT}_{j\rightarrow i}\textsf{CNOT}_{i\rightarrow j}$. Furthermore using the circuit identity $\textsf{CNOT}_{j\rightarrow i}=H_{j}\textsf{CZ}_{ji}H_{j}$, we see that $\textsf{SWAP}_{ij}$ is diagonalized by $\textsf{CNOT}_{j\rightarrow i}H_{j}$ and has eigenvalue $(-1)^{x_{i}x_{j}}$. For the full circuit Fig S5.a one first applies a transversal set of CNOT gates across both registers followed by a single layer of Hadamard gates $H$ on the top register. Then the output is sampled and the average of the boolean function

is reported. The output bits on the top register are labeled by $s\in\{0,1\}^{n}$, while $t\in\{0,1\}^{n}$ is the output string on the lower register.

Note that the optimization problem, eqn. (6) is only concave, when the matrix $K\geq 0$ is positive semi-definite. It can happen, that the shot noise and other errors in the experiment lead to a $\hat{K}$ that is no longer positive semi-definite. We have indeed observed this multiple times in the experiment. A possible way of dealing with this problem is a method developed in smolin2012efficient , where an optimization problem is solved to find the closest positive semi-definite $K$-matrix in trace norm to $\hat{K}$ consistent with the constraint. In our experiments however, we have found this not to be necessary and the performance has been almost optimal without performing this method.

Device parameters

Our device is fabricated on a 720-$\mu$m-thick Si substrate. A single optical lithography step is used to define all CPW structures and the qubit capacitors with Nb. The Josephson junctions are patterned via electron beam lithography and made by double-angle deposition of Al.

The dispersive readout signals are amplified by Josephson Parametric Converters Bergeal2010 (JPC). Both the quantum processor and the JPC amplifiers are thermally anchored to the mixing chamber plate of a dilution refrigerator.

The two qubit fundamental transition frequencies are $\omega_{i}/2\pi=\{5.2760(4),5.2122(3)\}$ GHz, with anharmonicities $\Delta_{i}/2\pi=\{-330.3,-331.9\}$ MHz, where $i\in\{0,1\}$. The readout resonator frequencies used are $\omega_{Ri}/2\pi=\{6.530553,6.481651\}$ GHz, while the CPW bus resonator connecting $Q_{0}$ and $Q_{1}$ was unmeasured and designed to be 7.0 GHz. The dispersive shifts and line widths of the readout resonators are measured to be $2\chi_{i}/2\pi=\{-1.06,-1.02\}$ MHz and $\kappa_{i}/2\pi=\{661,681\}$ kHz, respectively.

The two qubit lifetimes and coherences were measured intermittently throughout our experiments. The observed mean values were $T_{1(i)}=\{55,38\}$, $T^{*}_{2(i)}=\{16,17\}$, $T^{\rm{echo}}_{2(i)}=\{43,46\}$ $\mu$s with $i\in\{0,1\}$

Gate characterization

Our experiments use calibrated $X-$rotations ($X_{\pi}$ and $X_{\pi/2}$) as single-qubit unitary primitives. $Y-$rotations are attained by appropiate adjustment of the pulse phases, whereas $Z-$rotations are implemented via frame changes in software mckay2017efficient . A time buffer is added at the end of each physical pulse to mitigate effects from reflections and evanescent waves in the cryogenic control lines and components.

We use two sets of gate times in order to perform the Richardson extrapolation of the noise in our system temme2017error . For the first set of gate times we use 83 ns for single-qubit gates and 333 ns for each cross-resonance pulse. The buffer after each physical pulse is 6.5 ns. The single-qubit gates are gaussian shaped pulses with $\sigma=20.75$ ns. The cross-resonance gates are flat with turn-on and -off gaussian shapes of $\sigma=10$ ns. Our implementation of a CNOT has a duration of two single-qubit pulses and two cross-resonance pulses, giving a total of 858 ns for the first set of gate times, including buffers. For our second set of gate times we use the times of the first set but stretched by a factor of 1.5, including the pulses $\sigma$s and the buffers. This gives a total CNOT time of 1.287 $\mu$s and single-qubit gates of $\sim 125$ ns.

We experimentally verified our single- and two-qubit unitaries by Randomized Benchmarking (RB) Gambetta2012 ; Corcoles2013 . The following table shows the RB results for all single-qubit gates used in our experiments, including individual and simultaneous RB.

Our two-qubit unitarias are CNOTs constructed from echo cross-resonance sequences Corcoles2013 ; Sheldon2016 . Each of the two cross-resonance pulses in a CNOT has durations of 333 and 500 ns for the two different gate lengths used in our experiments. For our two-qubit RB we obtain a CNOT error of $0.0373\pm.0015$ ($0.0636\pm.0021$) for the 333 (500) ns cross-resonance pulse.

Readout correction

Our readout assigned fidelity was $\sim 95\%$ for both qubits.

For each experiment, we run 4 ($2^{2}$) calibration sequences preparing our two qubits in their joint computational states. We gather statistics of these calibrations and create a measurement matrix $A_{ij}=P(|i\rangle||j\rangle)$ where $P(|n\rangle||m\rangle)$ is the probability of measuring state $|m\rangle$ having prepared state $|n\rangle$. We then correct the observed outcomes of our experiments by inverting this matrix and multiplying our output probability distributions by this inverse.

Support vectors

Here we show the support vectors and $\alpha_{i}$ as calculated for each of the three datasets studied from their $K$ matrices.

Error mitigation for kernel estimation

The experimental estimation of the kernel matrices shown in Fig. S7 and in Fig. 4 in the main text involves running the experiments at different gate lengths and extrapolating the expectation value of the observable of interest to its zero-noise value. While this technique can be extremely powerful in scenarios where the noise is invariant under time rescaling, it is particularly sensitive to measurement sampling noise. In many cases it is the experimental readout assignment fidelity that determines the bound on how precisely the observable can be estimated.

Even though for our Sets I and II we attain 100 % classification success over 10 randomly drawn test sets in each case, we can quantify how close our experimentally determined separating hyperplane is to the ideal.

The optimal hyperplane for a given training set can be expressed as the linear combination $\sum_{i}\alpha_{i}y_{i}\vec{x}_{i}={\bf w}$ (eqn. 7), where $\bf{w}$ is a vector orthogonal to the optimal separating hyperplane and $\vec{x}_{i}$ are the support vectors. We can therefore quantify the distance between the experimentally obtained hyperplane and the ideal hyperplane by calculating the inner product $\langle{\bf{w}},{\bf{w}_{\textrm{ideal}}}\rangle=\sum_{i\in N_{S}}\sum_{j\in N_{S}^{\prime}}y_{i}y_{j}\alpha_{i}^{*}\alpha_{j}|\langle\vec{x}_{i}\vec{x}_{j}\rangle|^{2}/||\bf{w}||||\bf{w}_{\textrm{ideal}}||$ where $N_{S}$ and $N_{S}^{\prime}$ are the sets of experimentally obtained and ideal support vectors, respectively.

In Fig. S8 we show the inner products between the ideal and all experimental hyperplanes, including the two sets of gate times used throughout our experiments, $c1$ and $c1.5$, as well as the error-mitigated hyperplanes.

For Sets I and II, which classify at 100 % success, it is clear that error mitigation improves our results very significantly. This is not the case for Set III, which classifies at 94.75 % success. In fact, for Set III we see that error-mitigation worsens the hyperplane, as the results are closer to the ideal for the unmitigated experiments than both Sets I and II.

A look at the calibration data taken along the direct kernel estimation experiments for each set, we see that the readout assignment fidelities of $Q_{0}$($Q_{1}$) are 96.56% (96.31%) for Set I, 95.90% (96.36%) for Set II, and 93.99% (95.47%) for Set III. The slightly worse readout fidelities for Set III could partially explain the worse classification results for this set, but other aspects of the protocol might also contribute to this, like for example some gates operating on a somewhat non-linear regime after the calibrations in this set.

Another symptom for the degree of classification success in a given set can be observed by looking at the combined weight of negative eigenvalues in the kernel matrix, $\sum_{k_{i}<0}|k_{i}|$, with $k_{i}$ the eigenvalues of the kernel. We obtain 1.40, 1.27 and 2.41 for Sets I, II, and III, respectively.