# Parallel hybrid quantum-classical machine learning for kernelized time-
series classification

Jack S. Baker1 Gilchan Park2 Kwangmin Yu2<EMAIL_ADDRESS>Ara Ghukasyan1 Oktay
Goktas1 Santosh Kumar Radha1<EMAIL_ADDRESS>1Agnostiq Inc., 325 Front St
W, Toronto, ON M5V 2Y1  2Computational Science Initiative, Brookhaven National
Laboratory, Upton, New York 11973, USA

###### Abstract

Supervised time-series classification garners widespread interest because of
its applicability throughout a broad application domain including finance,
astronomy, biosensors, and many others. In this work, we tackle this problem
with hybrid quantum-classical machine learning, deducing pairwise temporal
relationships between time-series instances using a time-series Hamiltonian
kernel (TSHK). A TSHK is constructed with a sum of inner products generated by
quantum states evolved using a parameterized time evolution operator. This sum
is then optimally weighted using techniques derived from multiple kernel
learning. Because we treat the kernel weighting step as a differentiable
convex optimization problem, our method can be regarded as an end-to-end
learnable hybrid quantum-classical-convex neural network, or QCC-net, whose
output is a data set-generalized kernel function suitable for use in any
kernelized machine learning technique such as the support vector machine
(SVM). Using our TSHK as input to a SVM, we classify univariate and
multivariate time-series using quantum circuit simulators and demonstrate the
efficient parallel deployment of the algorithm to 127-qubit superconducting
quantum processors using quantum multi-programming.

Quantum multi-programming

††preprint: APS/123-QED

##  I Introduction

Processes and systems which produce observable characteristics evolving with
time are present in topics as disparate as finance, sensor technologies,
medicine, astronomy and many others. As a result, new techniques for time-
series analysis have become among the most sought after in machine learning
(ML) where we seek to learn temporal trends and correlations from time-series
data to perform classification [1], anomaly detection [2, 3], regression [4],
forecasting [5, 6] and to generate synthetic time-series instances [7].
Focusing on classification, classical ML has provided a zoo of algorithms
where the present state-of-the-art is centered around deep learning. Popular
approaches include recurrent neural networks like long-short-term-memory [8]
(and its variations [9, 10]), gated recurrent units and, more recently,
transformer networks [11, 12, 13].

While popular for static data (i.e. data without time dependence), kernelized
methods like the well-known support vector machine (SVM) [14] see limited use
in time-series analysis. Although time-series kernel methods have been
proposed [15, 16, 17, 18] and used as analytical tools for interpreting deep
learning models [19], none are tailored to encode temporal trends by way of an
explicit and learnable time-dependent inner product space. Within classical
ML, it is not clear how such a space can be constructed in a non-trivial
manner. In this work, we show that time evolution as generated by a
parameterized Hamiltonian operator within quantum mechanics is a natural
approach for achieving such an objective. Indeed, by combining inner products
evolved to different points of time, we achieve a quantum kernel function
adapted for time-series data which we call the Time-Series Hamiltonian Kernel
(TSHK). Furthermore, we discuss situations where training a TSHK is useful and
provide theoretical tools to analyze the time dependence of the learnt inner
product spaces.

The construction of the TSHK is tied to a field known as Multiple Kernel
Learning (MKL). In MKL, multiple kernel functions are combined using various
tactics [20, 21] to give rise to a combined kernel that is by some measure
superior. The TSHK is a combined kernel built using a weighted linear
combination of quantum kernel functions each defined at different instances of
time $t$. These weights are chosen to maximize the separation between labeled
data classes as was proposed in the well-known EasyMKL [21] algorithm.
Furthermore, as was originally suggested in [22], the variational parameters
of quantum kernels and the kernel weights can be obtained simultaneously by
training an end-to-end learnable Quantum-Classical-Convex neural network (QCC-
net). The network producing the TSHK can then be regarded as a temporally-
aware variation of the original QCC-net.

Currently in the noisy intermediate-scale quantum (NISQ) era [23], various
sources of noise prevent the implementation of deep quantum circuits on real
hardware. In light of this, we must ensure that our quantum circuits
implementing the TSHK are sufficiently shallow. Conventionally, however, the
required time evolution operators are implemented using a potentially high-
depth Lie-Suzuki-Trotter expansion [24]. Fortunately, built upon the
eigendecomposition approach originally used in the variational fast forwarding
algorithm [25], a family of time evolution operators can be implemented with
low depth variational circuits. The space of accessible time evolution
operators can controlled by modifying the depth and structure of the circuit
Ansätz implementing the eigendecomposition. Several works [26, 27, 28, 29, 30]
have now appeared leveraging this approach to: learn constraint-enforcing
mixing operators in the Quantum Approximate Optimization Algorithm [26],
generate synthetic time-series instances [27] and detect anomalous behaviour
in time-series data [28]. The generalization bounds of this approach have also
been studied [29, 30].

Although the eigendecomposition approach does allow for a shallow circuit
implementation of the TSHK, the number of computations required to compute the
kernel matrices scales with the square of the size of the training data set
and linearly with the length of the time series. Fortunately, the computations
of kernel matrix elements are independent of one another allowing them to be
computed in parallel. Indeed, a new approach to efficiently utilize
contemporary NISQ computers is by overlapping multiple quantum circuits [31].
This method, called Quantum Multi-Programming (QMP), utilizes NISQ devices by
executing multiple quantum circuits concurrently. In this work, we use QMP to
compute the TSHK in parallel. Compared to serial execution, we achieve a
notable speed-up (at least $35$ times) with QMP without loss of accuracy using
two 127-qubit IBM quantum computers. In light of this large speed-up, this
work demonstrates the practical utility of parallelism in hybrid quantum-
classical ML workflows.

The rest of this work is organized as follows: in Sec. II we describe how we
achieve a time-dependent inner product space and show how a weighted sum of
inner products can be used to construct the TSHK (Sec. II.1). Also within this
Section, we show how the TSHK can be trained using a QCC-net (Sec. II.2) and
show how the resulting TSHK can be plugged in to a SVM to perform time-series
classifications (Sec II.3). In Sec. III, we describe in detail what is meant
by time-dependence in the trained kernel function and show how one can
efficiently probe this property. In Sec. IV, we demonstrate SVM classification
using the TSHK for a synthetic multivariate time-series (Sec. IV.1) and for a
real univariate time-series: the well-known gun-point data set [32] (Sec.
IV.2). Our experiments using QMP on 127 qubit superconducting transmon chips
are contained within Sec. V. The QMP approach is described in detail in Sec.
V.1 and pertinent QMP design considerations are explained in Sec. V.2. The
results of our QMP experiments are shown and discussed in Sec. V.3. Finally,
in Sec. VI, we conclude and remark on future directions in hybrid quantum-
classical ML exploiting parallelism.

##  II Algorithm description

Figure 1: Optimizing the combined quantum kernel for a simplified sinusoid
versus cosinusoid classification problem with only one training example from
class $y_{i}=1$ [$-\sin(t)$] and class $y_{i}=-1$ [$-\cos(t)$]. (i) Given some
initial $\bm{\theta}=\bm{\theta}_{\text{init}}$, calculate the kernels of Eq.
5 at each time $t_{l}$ and calculate the equally weighted Gram matrix
$K_{e}(\bm{\theta})=\sum_{t\in T}K_{t}(\bm{\theta})$. Unitary operator
diagrams are color-coded in the Figure and because there is only one unique
element (excluding the diagonal which is unity by definition) only one unitary
diagram is shown for each kernel computation. (ii) Solve the convex
optimization problem with respect to the Lagrangian dual variables $\bm{\phi}$
using a cone program. (iii) Calculate the quantum gradients
$\partial\mathcal{L}_{\phi_{\text{min}}}(\bm{\theta})/\partial\bm{\theta}$ by
differentiating through the cone program. The surface plot shows a convex
optimization landscape spanned by $\phi_{1}$ and $\phi_{2}$. (iv) Use the
calculated
$\partial\mathcal{L}_{\phi_{\text{min}}}(\bm{\theta})/\partial\bm{\theta}$ to
step forward a derivative-requiring optimization routine (i.e. gradient
descent, Adam etc. ). If the termination criteria (a fixed number of loss
iterations, convergence tolerance for $\mathcal{L}_{\phi_{\text{min}}}$ etc.)
are met with this new step, advance to (v) else, go back to (i) and enter new
loop iteration. Note the negative sign of
$\mathcal{L}_{\phi_{\text{min}}}(\bm{\theta})$ in part (iv) of the figure.
This is because in practice we maximize the function by minimizing its
negative value. (v) Extract the optimal kernel weights $\bm{\eta}^{\star}$
using Eq. 11 to yield the combined kernel function (the TSHK) from Eq. 6. This
kernel function can now be used to calculate kernel matrices
$K(\bm{\theta}^{\star},\bm{\eta}^{\star})$ for use in any kernelized learning
task including SVM.

###  II.1 Time-dependent inner product space

We begin with a formal definition of the type of data we are working with:
time-series. A time-series $\bm{x}$ is a sequence of $p\in\mathbb{Z}^{+}$
observations from a process/system arranged in chronological order

$\bm{x}:=(\bm{x}_{t}:t\in T),\quad T:=(t_{l}:1\leq l\leq p)$ (1)

where
$\bm{x}_{t}=[{}^{1}x_{t},{}^{2}x_{t},\ldots{}^{d}x_{t}]\in\mathbb{R}^{d}$,
$d\in\mathbb{Z}^{+}$, $l\in\mathbb{Z}^{+}$ and time $t\in\mathbb{R}^{+}$. In
what follows, we devise a technique for obtaining a trainable kernel function
adapted for time-series data of the form defined above by (i) formulating a
time-dependent quantum inner product space and (ii) using methods from the
classical MKL literature [21, 33] to find an optimally weighted sum of inner
products from different points in time.

The first step towards achieving our goal is to inject time dependence into a
quantum state. We achieve this through operating on the $n$-qubit state of
zeroes $|0\rangle^{\otimes n}$ using a parameterized unitary time evolution
operator $e^{-iH(\bm{\beta},\bm{\gamma})t^{\prime}}$ where
$t^{\prime}\in\mathbb{R}^{+}$ followed by a parameterized unitary matrix
$U(\bm{x}_{t},\bm{\alpha})$ used to embed $\bm{x}_{t}$ within the quantum
state

$|\bm{x}_{t},t^{\prime},\bm{\theta}\rangle:=U(\bm{x}_{t},\bm{\alpha})e^{-iH(\bm{\beta},\bm{\gamma})t^{\prime}}|0\rangle^{\otimes
n}$ (2)

for real-valued set of free parameter vectors
$\bm{\theta}=\\{\bm{\alpha},\bm{\beta},\bm{\gamma}\\}$ where the embedding
unitary matrix $U(\bm{x}_{t},\bm{\alpha})$ is permitted to be any quantum
feature map [34, 35] implemented using a layered parameterized quantum
circuit. It should be stressed that $U(\bm{x}_{t},\bm{\alpha})$ is a static
map. That is, despite the subscript $t$ bared by the argument $\bm{x}_{t}$,
the resulting unitary matrix is time-independent. Explicitly, should we have
$\bm{x}_{t_{1}}=\bm{x}_{t_{2}}$ , $t_{1}\neq t_{2}$ the unitary matrices
$U(\bm{x}_{t_{1}},\bm{\alpha})$ and $U(\bm{x}_{t_{2}},\bm{\alpha})$ are
indistinguishable. Time-dependence only comes as a result of the time
evolution operator $e^{-iH(\bm{\beta},\bm{\gamma})t^{\prime}}$ which we write
as an eigendecomposition

$V_{t^{\prime}}(\bm{\beta},\bm{\gamma}):=W^{\dagger}(\bm{\beta})D(\bm{\gamma},t^{\prime})W(\bm{\beta})=e^{-iH(\bm{\beta},\bm{\gamma})t^{\prime}}$
(3)

for parameterized unitary matrix of eigenvectors $W(\bm{\beta})$ and
parameterized time-encoded diagonal unitary matrix
$D(\bm{\gamma},t^{\prime})$. The above equality with
$e^{-iH(\bm{\beta},\bm{\gamma})t^{\prime}}$ holds following Stone’s theorem
for strongly continuous one-parameter unitary groups [36]. The strength of
this approach lies in implementing each unitary matrix in the
eigendecomposition as a layered parameterized quantum circuit allowing one to
control the search space of unitary time evolution operators. Should a small
number of layers prove sufficient to optimize the forthcoming loss function
(Eq. 8), our algorithm becomes suitable for present generation NISQ computers.
For notational short-hand, it is convenient to write a combined unitary matrix
as a product of the static and time-dependent parts

$\Omega(\bm{x}_{t},t^{\prime},\bm{\theta}):=U(\bm{x}_{t},\bm{\alpha})V_{t^{\prime}}(\bm{\beta},\bm{\gamma}).$
(4)

We can now begin to formally derive a time-dependent inner product. To do so,
we first must regard the density matrix
$\rho_{t^{\prime}}(\bm{x}_{t},\bm{\theta})=|\bm{x}_{t},t^{\prime},\bm{\theta}\rangle\langle\bm{x}_{t},t^{\prime},\bm{\theta}|$
as the unambiguous feature vector embedding the classical data point
$\bm{x}_{t}$ (regarding $|\bm{x}_{t},t^{\prime},\bm{\theta}\rangle$ as the
feature vector is ambiguous as the space $\mathbb{C}^{2^{n}}$ is only
physically defined up to a global phase [35]) . We now define a kernel by
setting $t^{\prime}\rightarrow t$ and taking the Frobenius inner product
$\langle\rho_{t}(\bm{x}^{\prime}_{t},\bm{\theta}),\rho_{t}(\bm{x}_{t},\bm{\theta})\rangle$.
Equivalently,

$\kappa_{t}(\bm{x}_{t},\bm{x}^{\prime}_{t},\bm{\theta}):=\text{Tr}[\rho_{t}(\bm{x}^{\prime}_{t},\bm{\theta})\rho_{t}(\bm{x}_{t},\bm{\theta})]\\\
=|\langle\bm{x}^{\prime}_{t},t,\bm{\theta}|\bm{x}_{t},t,\bm{\theta}\rangle|^{2}$
(5)

where $\bm{x}^{\prime}_{t}$ is an element of another time-series instance
$\bm{x}^{\prime}$ defined analogously to Eq. 1. Eq. 5 can be estimated
practically by preparing the state
$\Omega^{\dagger}(\bm{x}^{\prime}_{t},t,\bm{\theta})\Omega(\bm{x}_{t},t,\bm{\theta})|0\rangle^{\otimes
n}$ and taking the expectation value of the projector on the state of zeros
$P_{0}=|0\rangle^{\otimes n}\langle 0|^{\otimes n}$. This requires repeated
preparation and measurement of
$\Omega^{\dagger}(\bm{x}_{t},t,\bm{\theta})\Omega(\bm{x}_{t},t,\bm{\theta})|0\rangle^{\otimes
n}$ in the computational basis to estimate the probability of measuring the
$n$-bit string of zeros $00\ldots 0$ .

Eq. 5 fulfils the symmetry requirement of a valid kernel function. That is, we
have
$\kappa_{t}(\bm{x}_{t},\bm{x}^{\prime}_{t},\bm{\theta})=\kappa_{t}(\bm{x}^{\prime}_{t},\bm{x}_{t},\bm{\theta})$
for any choice of the two arguments $\bm{x}_{t}$,
$\bm{x}^{\prime}_{t}\in\mathbb{R^{d}}$ or the time argument $t$. Indeed, the
achievement of this symmetry motivates the order in which unitary operations
are applied in Eq. 2. If we were to switch the order in which we apply of
$U(\bm{x}_{t},\bm{\alpha})$ and $e^{-iH(\bm{\beta},\bm{\gamma})t^{\prime}}$ to
$|0\rangle^{\otimes n}$ we yield the state
$|\bm{x}_{t},t^{\prime},\bm{\theta}\rangle_{G}:=e^{-iH(\bm{\beta},\bm{\gamma})t^{\prime}}U(\bm{x}_{t},\bm{\alpha})|0\rangle^{\otimes
n}$. Computing the inner product of this state at $t^{\prime}=t$ with another
state at $t^{\prime}=t^{\prime}$ and multiplying the result by its own complex
conjugate, we yield the quantity $\mathbb{G}:=|{}_{G}\langle
x_{t},\bm{\alpha}|e^{-iH(t-t^{\prime})}|x^{\prime}_{t},\bm{\alpha}\rangle_{G}|^{2}$
where we have taken
$|x_{t},\bm{\alpha}\rangle_{G}:=U(\bm{x}_{t},\bm{\alpha})|0\rangle^{\otimes
n}$. We see immediately that for $t\neq t^{\prime}$ this expression is not a
valid kernel as the symmetry requirement detailed above is violated. It is
only at $t=t^{\prime}$ that symmetry is recovered but at the expense of losing
time dependence entirely as the time evolution operator becomes the identity.
It is worth pointing out that while a time-dependent kernel function is not
achievable in this way, the result is otherwise important. That is,
interestingly, $\mathbb{G}$ is equivalent to the product of the Green’s
function propagator (hence the subscript $G$ in the above) from quantum field
theory [37] with its complex conjugate
$G_{\bm{\theta}}(\bm{x}_{t},t;\bm{x}^{\prime}_{t},t^{\prime})G^{\star}_{\bm{\theta}}(\bm{x}_{t},t;\bm{x}^{\prime}_{t},t^{\prime})$.
Using parameterized Green’s function propagators could prove an interesting
avenue for non-kernelized time-series analysis in QML but is beyond the scope
of this work.

Returning now to the kernel function defined in Eq. 5, we compute it $\forall
t\in T$ and combine the $p$-terms in a weighted sum to achieve a combined
kernel function measuring the similarity between a pair of time-series
instances $(\bm{x},\bm{x}^{\prime})$

$\kappa(\bm{x},\bm{x}^{\prime},\bm{\theta},\bm{\eta}):=\sum_{t\in
T}\eta_{t}\cdot\kappa_{t}(\bm{x}_{t},\bm{x}^{\prime}_{t},\bm{\theta})$ (6)

for kernel coefficient vector $\bm{\eta}:=[\eta_{t}|t\in T]\in\mathbb{R}^{p}$
subject to the constraint $\sum_{t\in T}\eta_{t}=1$. We call Eq. 6 the TSHK.
The motivation for such a form of kernel function is to allow the measured
similarity between a pair of time-series $(\bm{x},\bm{x}^{\prime})$ to be
influenced by differently weighted contributions at distinct points in time.
Should the resulting kernel be used in a classification algorithm like SVM,
this allows for “smoking gun” values of $t$ to discriminate between different
classes of time-series. In the next Section, we present a strategy for
optimally setting the kernel parameters $\bm{\theta}$ and $\bm{\eta}$ where
the exact sense of optimality is to be defined.

Before we advance to the next Section, however, now is a good time to discuss
what structures of time series data may be learnt well by the kernel function
of Eq. 6 and for which we may expect performance to be no better than time
independent kernels. It is clear from the form of Eq. 6 and the time-dependent
kernels which construct it (Eq. 5) that great emphasis is placed on the notion
of a “shared time axis”; a significant inductive bias of our approach. That
is, we expect a learnable TSHK when each instance of the training time series
data begins at a common initial condition in time. Examples of such data
include the temperature in different locations for each day of the year (where
we may classify the southern hemisphere versus the northern hemisphere), or
traffic (the number of detected vehicles) on roadways given at some frequency
on different days (where we may classify weekday versus weekend traffic). In
absence of such a shared time axis, we may not expect a learnable TSHK. For
example, if we are given “windowed data” (sub-sequences taken from a larger
time series) the initial condition is not guaranteed to be shared among data
instances. One example of this is the classification of different genres of
music given small audio clips from a longer song. It should be noted that even
in this case, it could be possible increase the learnability of the TSHK by
preprocessing time series using a time-axis alignment protocol like Dynamic
Time Warping [38]. The success of such an approach would have to be assessed
on a case-by-case basis.

###  II.2 Training the kernel function

We now present a strategy for training the kernel function of Eq. 6 in a
supervised learning setting. Throughout this Section, Fig. 1 should be used as
a companion and will be referenced throughout. To begin, we assume to have
access to training set $X$ containing $N_{X}\in\mathbb{Z}^{+}$ time-series
instances $\bm{x}^{(i)}$ such that

$X:=\\{\bm{x}^{(i)}:i\in\mathbb{Z}^{+}_{\leq N_{X}}\\}$ (7)

where each time-series instance $\bm{x}^{(i)}$ has a corresponding binary
class label $y_{i}\in\\{1,-1\\}$. Using this labeled data, we define a
suitable loss function to be optimized.

For optimizing the kernel coefficients $\bm{\eta}$, it was proposed for the
EasyMKL algorithm [21] that the objective function from the Kernel
Optimization of the Margin Distribution (KOMD) [39] approach should be used.
The KOMD objective represents the separation between positive ($y_{i}=1$) and
negative ($y_{i}=-1$) training examples where the kernel weights are
represented implicitly with the Lagrangian dual variables
$\bm{\phi}:=[\phi_{1},\phi_{2},\ldots\phi_{N_{X}}]\in\mathbb{R}^{N_{X}}$. That
is, $\bm{\eta}$ is derivable using $\bm{\phi}$ as is shown later in Eq. 11
(see [21] for explicit steps in transforming the optimization problem in
$\bm{\eta}$ to that of $\bm{\phi}$ and vice-versa). Recalling that our kernel
function is also parameterized by the unitary operator parameters
$\bm{\theta}$, the separation between classes as defined in [21] is modified
to be written as a function of both parameter sets $\bm{\theta}$ and
$\bm{\phi}$:

$\mathcal{L}(\bm{\theta},\bm{\phi}):=(1-\lambda)\bm{\phi}^{T}\hat{Y}K_{e}(\bm{\theta})\hat{Y}\bm{\phi}+\lambda||\bm{\phi}||^{2}_{2}$
(8)

where $K_{e}(\bm{\theta})=\sum_{t\in T}K_{t}(\bm{\theta})$ is an equally
weighted sum of Gram matrices with elements
$K^{ij}_{t}(\bm{\theta})=\kappa_{t}(\bm{x}_{t}^{(i)},\bm{x}_{t}^{(j)},\bm{\theta})$
from Eq. 5, $\hat{Y}=\text{diag}(y_{1},y_{2},\ldots,y_{N_{x}})$ and
$\lambda\in[0,1]$ is a real-valued penalty scaling factor hyperparameter
intended to favor low variance solutions for $\bm{\phi}$ with increasing
$\lambda$. As shown in Fig. 1(i), $K_{e}(\bm{\theta})$ is initially computed
with some $\bm{\theta}_{\text{init}}$ (chosen randomly or otherwise) which is
used to prime the first step of the max-min optimization problem

$\mathcal{L}^{\star}:=\text{max}_{\bm{\theta}}[\text{min}_{\bm{\phi}}\mathcal{L}(\bm{\theta},\bm{\phi})]$
(9)

where $\mathcal{L}^{\star}$ is the optimal loss and the optimal arguments are
$\bm{\theta}^{\star}$ and $\bm{\phi}^{\star}$, respectively. Because the
internal minimization with respect to $\bm{\phi}$ is a convex objective, as
shown in Fig. 1(ii), for any given fixed
$\bm{\theta}=\bm{\theta}_{\text{fixed}}$, we can efficiently and accurately
find $\text{min}_{\bm{\phi}}[\mathcal{L}(\bm{\theta},\bm{\phi})]$ using a
convex solver such as a cone program [40]. Because we can reasonably assume
that the global minimum of a convex problem will be found, we can re-frame the
problem as

$\mathcal{L}^{\star}=\text{max}_{\bm{\theta}}[\mathcal{L}_{\bm{\phi}_{\text{min}}}(\bm{\theta})]$
(10)

where
$\mathcal{L}_{\bm{\phi}_{\text{min}}}(\bm{\theta}):=\text{min}_{\bm{\phi}}[\mathcal{L}(\bm{\theta},\bm{\phi})]$
and the subscript $\bm{\phi}_{\text{min}}$ indicates that we are working at
the $\bm{\phi}$ which minimizes $\mathcal{L}(\bm{\theta},\bm{\phi})$ with
respect to $\bm{\phi}$ for any given $\bm{\theta}$.

The maximization with respect to $\bm{\theta}$ can be performed using any
general purpose optimization algorithm. Importantly, the gradient
$\partial\mathcal{L}_{\bm{\phi}_{\text{min}}}(\bm{\theta})/\partial\bm{\theta}$
is obtainable analytically through an automatic differentiation engine [41]
because (i) partial gradients of the kernel function with respect to the
quantum circuit parameters $\bm{\theta}$ are known using techniques including
the parameter shift rule and its variants [42, 43, 41] and (ii) the partial
gradients of the minimum of the convex optimization problem with respect to
the quantum circuit parameters are known because the cone program is
differentiable [44, 45] [Fig. 1(iii)]. This allows one to update the quantum
parameters $\bm{\theta}$ [Fig. 1(iv)] using one of a whole host of gradient-
requiring optimizers available from the wider ML literature [46]. Indeed, our
method can be interpreted as a hybrid neural network with quantum layers and
classical convex optimization layers [44]: a Quantum-Classical-Convex network
(QCC-net).

After meeting the termination criteria of the optimization (which could be a
fixed number of iterations or tolerance on changes in the value of Eq. 8 when
the parameters are updated, for example), we arrive at an approximation of
$\bm{\theta}^{\star}$ and $\bm{\phi}^{\star}$ where the latter is given as
$\text{argmin}_{\bm{\phi}}[\mathcal{L}(\bm{\theta}^{\star},\bm{\phi})]$. We
use both optimal parameter vectors to extract the optimal kernel weights
$\eta_{t}^{*}$

$\eta_{t}^{*}=\frac{\bm{\phi}^{\star
T}\hat{Y}K_{t}(\bm{\theta}^{\star})\hat{Y}\bm{\phi^{\star}}}{\sum_{t\in
T}\bm{\phi}^{\star
T}\hat{Y}K_{t}(\bm{\theta}^{\star})\hat{Y}\bm{\phi^{\star}}}.$ (11)

Now, in Eq. 6 we set $\bm{\theta}\rightarrow\bm{\theta}^{\star}$ and
$\bm{\eta}\rightarrow\bm{\eta}^{\star}$ to obtain our goal: a trained TSHK
[Fig. 1(v)].

It should also be noted that during the training cycle, it is often
appropriate to work with mini-batches of the training data. This both reduces
the computational effort required in training and is known to improve
generalization performance [47] in classical ML algorithms. Formally, a mini-
batched iteration of the training cycle involves sampling, at random, a
training mini-batch $X_{\text{batch}}\subset X$ with cardinality
$|X_{\text{batch}}|=N_{\text{batch}}$. During this cycle, Gram matrices
$K_{t}(\bm{\theta})$ and $K_{e}(\bm{\theta})$ and the diagonal label matrix
$\hat{Y}$ have the reduced dimension
$\mathbb{R}^{N_{\text{batch}}}\times\mathbb{R}^{N_{\text{batch}}}$ and the
Lagrangian parameter vector $\bm{\phi}$ reduces its dimension to
$\mathbb{R}^{N_{\text{batch}}}$. It should be clear now that the wall-time for
a single mini-batch iteration using a fixed number of shots to evaluate the
quantum kernel elements and a fixed number of iterations to solve the convex
problem scales quadratically with $\mathbb{R}^{N_{\text{batch}}}$. When these
batches are sampled, we must ensure that each mini-batch contains at least one
time-series instance with label $y_{i}=1$ and another with $y_{i}=-1$. Without
this, we cannot define a separation between classes and therefore the loss
function of Eq. 8 is undefined.

Before moving on, we must consider the trainability of the TSHK. Among the key
ingredients of the TSHK are of course quantum kernels. It is known that for
static quantum kernels, quantum circuit gradients vanish exponentially with
the number of qubits when classical data is embedded using Ansätze without
sufficient inductive biases [48]. This is one of the manifestations of the
barren plateaus phenomenon [49]. Within our QCC-net framework, because quantum
circuit gradients are obtained by differentiating through convex optimization
layers, it is unclear to what extent QCC-nets suffer from barren plateaus.
Indeed, not all quantum neural network architectures are found with this
problem [50]. An interesting topic of future study would be a theoretical
analysis of gradient scaling in QCC-nets and how to use convex optimization
layers to impart different inductive biases.

###  II.3 Using the trained kernel function: classification

Now equipped with an optimal kernel function, we can use it in any kernelized
ML technique like kernelized SVM or kernel ridge regression. In this work we
focus on classification using the kernelized soft margin variation of SVM
which from now on-wards we refer to as “SVM”. While a full presentation of the
formalism describing SVM is beyond the scope of this work (and is discussed
elsewhere [14]), we now provide a high level explanation. In binary
classification, SVM separates data points belonging to two different classes
by learning a decision boundary, also known as a hyperplane, from labelled
training data. It allows for some misclassifications by introducing a soft
margin, where the appetite for misclassfication is controlled by a single
hyperparameter $C\in\mathbb{R}$. A kernel function is used to map the data
into a higher-dimensional space where the hyperplane can separate the classes
more easily.

For our purposes, SVM can be regarded as a decision function for an unseen
data point $\bm{x}$ given four items:
$\kappa(\bm{x},\bm{x}^{\prime}),X,\bm{y}$ and $C$ where each, in order, is a
valid kernel function, a set of training data, the vector of binary class
labels for that training data and the misclassification appetite
hyperparameter. That is, we have

$\mathcal{D}(\bm{x}|\kappa(\bm{x},\bm{x}^{\prime}),X,\bm{y},C):=\mathcal{D}^{\star}(\bm{x})\in\mathbb{R}.$
(12)

Binary class predictions can now be made using

$y_{\text{pred}}(\bm{x})=\text{sgn}[\mathcal{D}^{\star}(\bm{x})]\in\\{-1,1\\}.$
(13)

The most natural way to use this framework to perform time-series
classifications is simply to set $\kappa(\bm{x},\bm{x}^{\prime})$ to a trained
kernel function of the form of Eq. 6 at the optimal parameters
$\bm{\theta}^{\star}$ and $\bm{\eta}^{\star}$. Another way of doing so is the
break the problem down into $p$-many time-dependent SVMs. We then have

$\mathcal{D}_{t}(\bm{x}_{t}|\kappa_{t}(\bm{x}_{t},\bm{x}^{\prime}_{t},\bm{\theta}^{\star}),X_{t},\bm{y},C):=\mathcal{D}_{t}^{\star}(\bm{x}_{t})\in\mathbb{R}$
(14)

where $X_{t}=\\{\bm{x}^{1}_{t},\bm{x}^{2}_{t},\ldots\bm{x}^{N_{X}}_{t}\\}$.
This function is then used to make a prediction at each $t$ with
$y^{t}_{\text{pred}}(\bm{x}_{t})=\text{sgn}[\mathcal{D}_{t}^{\star}(\bm{x}_{t})]$.
The final class prediction for the entire time-series $\bm{x}$ is then given
by a weighted majority vote from each $t$

$y^{\text{vote}}_{\text{pred}}(\bm{x})=\text{sgn}\left[\sum_{t\in
T}\eta_{t}\cdot y^{t}_{\text{pred}}(\bm{x}_{t})\right]\in\\{-1,1\\}.$ (15)

While it is not guaranteed that $y_{\text{pred}}(\bm{x})$ will provide more
accurate predictions than $y^{\text{vote}}_{\text{pred}}(\bm{x})$ in general,
we note that because our training procedure defined in Sec. II.2 refines a
combined kernel [used to define $y_{\text{pred}}(\bm{x})$], it is intutive to
expect this. Indeed, we find that $y_{\text{{pred}}}(\bm{x})$ is most
effective for the data sets treated in Sec. IV.1 , IV.2 and V.3, so we choose
to make predictions using Eq. 13. However, analyzing
$\mathcal{D}_{t}^{\star}(\bm{x}_{t})$ and $y^{t}_{\text{pred}}(\bm{x}_{t})$
can provide valuable insights into the time-dependent nature of a given TSHK
function which are used in this work.

##  III Probing time dependence

After having achieved a reasonable approximation of $\bm{\theta}^{\star}$
using a classical optimization routine, a relevant quantity to probe is the
time dependence of the inner product space. This is an important question to
ask because while usage of the TSHK allows for a time dependent inner product
space, it does not mean that strong time dependence was achieved when
minimizing the cost function of Eq. 8 which will be strictly dependent on the
data set and the choice of parameters/circuit Ansätze used to construct the
QCC-net. As an extreme case, there can exist some $\bm{\beta}$ and
$\bm{\gamma}$ such that
$e^{-iH(\bm{\beta},\bm{\gamma})t}\approx\mathbb{I}\quad\forall t$ thus the
resulting TSHK is approximately time-independent. Naively, one may initially
think to examine the value of
$\kappa_{t}(\bm{x}_{t},\bm{x}^{\prime}_{t},\bm{\theta}^{\star})$ as a function
of $t$ where $\bm{x}_{t}$ and $\bm{x}^{\prime}_{t}$ are drawn from a testing
data set (or iterate over an entire testing set to calculate an entire kernel
matrix). While there is nothing wrong with this approach, we must note that in
the ideal case from training, time dependence in this sense will disappear.
That is, at each point in time $t\in\mathcal{T}$ the kernel function was
trained to embed classical data with label $y_{i}=1$ in a separate region of
Hilbert space to those with label $y_{i}=-1$. In the ideal case, these two
regions do not overlap at all and hence
$\kappa_{t}(\bm{x}_{t},\bm{x}^{\prime}_{t},\bm{\theta}^{\star})=0$ $\forall t$
when $\bm{x}_{t}$ and $\bm{x}^{\prime}_{t}$ belong to different classes.
Before the kernel function is trained, finite values of
$\kappa_{t}(\bm{x}_{t},\bm{x}^{\prime}_{t},\bm{\theta})$ are of course
permitted under the same conditions for the input. The above discussion is
best explained alongside Fig. 2. Intended for illustrative purposes only, Fig.
2 compresses high-dimensional vectors $|v\rangle\in\mathbb{C}^{2^{n}}$ into
two dimensional regions on the surface of a sphere. Differently colored
regions represent the space of feature vectors
$|\bm{x}_{t},t,\bm{\theta}\rangle$ accessible in different domains defined by
the classical data $\bm{x}_{t}$ (see the figure caption for more information).
Looking at Fig. 2(a), we can see that at any $t$, the time-dependent Hilbert
space regions $\mathcal{H}_{t}^{y_{i}=1}$ do not overlap with
$\mathcal{H}_{t}^{y_{i}=-1}$ when the kernel function is ideally trained
(regions marked “T”). For an untrained kernel (or at least not ideally trained
kernel; regions marked “U”), overlap between these regions are permitted.

To provide a more intuitive notion of time dependence, we must ask another
question: “given a general classical data vector
$\bm{\chi}=[{}^{1}\chi,{}^{2}\chi,\ldots{}^{d}\chi]\in\mathbb{R}^{d}$, does
the quantum state embedding of $\bm{\chi}$ given by Eq. 2 when we set
$\bm{x}_{t}\rightarrow\bm{\chi}$ change as $t\rightarrow t+\delta t$?”.

Figure 2: Time-evolving Hilbert space illustrations. Spheres are divided into
$Q$-many petals each indexed by $q\in\mathbb{Z}^{\neq 0}$. In each petal, the
grey area represents the vector space spanning $\mathbb{C}^{2^{n}}$ such that
each point on a petal represents a vector $|v\rangle\in\mathbb{C}^{2^{n}}$.

A full rotation of the sphere corresponds to the period $\mathcal{T}$ of time
evolution operator $e^{-iH(\bm{\beta,\gamma})t}$. We note that a period
$\mathcal{T}$ formally exists when the spectrum of $H$ is free of irrational
eigenvalues. The $q=0$ case is omitted to retain visual clarity in (b) where
colored regions would now (by definition of Eq. 16 and 17) maximally overlap.
(a) Hilbert space regions accessible by embedded classical data belonging to
class 1 or -1 where the kernel function has been Trained (T) or is untrained
(or at least not ideally trained; U) at different times $t=t_{q}$. (b) Hilbert
space regions accessible by the entire classical data input space
$\bm{\chi}\in\mathbb{R}^{d}$ at time $t$ and time $t+q\Delta t$ where $\Delta
t$ is a fixed time step.

The quantity of interest to probe this time dependence is the overlap of the
$t$ and $t+\delta t$ quantum embedding spaces integrated over all $\bm{\chi}$.
That is

$f_{\bm{\theta}}(\delta
t):=\int_{{}^{d}\chi_{\text{lo}}}^{{}^{d}\chi_{\text{hi}}}\ldots\int_{{}^{1}\chi_{\text{lo}}}^{{}^{1}\chi_{\text{hi}}}|\langle\bm{\chi},t,\bm{\theta}|\bm{\chi},t+\delta
t,\bm{\theta}\rangle|^{2}d^{1}\chi\ldots d^{d}\chi$ (16)

where the limits ${}^{m}\chi_{lo}$ and ${}^{m}\chi_{hi}$ are the minimum and
maximum values permitted for classical data embedded into the chosen quantum
feature space, respectively. Many quantum feature maps are periodic such that
${}^{m}\chi\in[{}^{m}\chi_{\text{lo}},{}^{m}\chi_{\text{hi}})$ where
${}^{m}\chi_{\text{lo}}=0$ and ${}^{m}\chi_{\text{hi}}=2\pi$. Should either
limit be unbound, limits should be set to enclose the zone of interest (i.e.
some region enclosing all training and testing examples) to the data set.
After simplifying the integrand using Eq. 2 (namely,
$U^{\dagger}(\bm{\chi},\bm{\alpha})U(\bm{\chi},\bm{\alpha})=\mathbb{I}$ and
cancelling $t$ in the exponent of the time evolution operator) we find

$f_{\bm{\beta},\bm{\gamma}}(\delta t):=\langle a_{\delta
t},\bm{\beta},\bm{\gamma}|P_{0}|a_{\delta
t},\bm{\beta},\bm{\gamma}\rangle\prod_{j=1}^{d}(^{j}\chi_{hi}-{}^{j}\chi_{lo})$
(17)

where $|a_{\delta
t},\bm{\beta},\bm{\gamma}\rangle:=e^{-iH(\bm{\beta},\bm{\gamma})\delta
t}|0\rangle^{\otimes n}$ and $P_{0}=|0\rangle^{\otimes n}\langle 0|^{\otimes
n}$. Noting also that the result of the product over the index $j$ in Eq. 17
is a real constant scaling factor $k$, it is useful to set $k=1$ choosing to
work with the quantity $\mathcal{F}_{\bm{\beta},\bm{\gamma}}(\delta t)=\langle
a_{\delta t},\bm{\beta},\bm{\gamma}|P_{0}|a_{\delta
t},\bm{\beta},\bm{\gamma}\rangle\in[0,1]$. Examining
$\mathcal{F}_{\bm{\beta},\bm{\gamma}}(\delta t)$, it is clear that we need
only consider the time evolution operator acting on the zero state with no
consideration of classical data points. Indeed, much like Eq. 5,
$\mathcal{F}_{\bm{\beta},\bm{\gamma}}(\delta t)$ can be efficiently computed
by estimating the probability of measuring the bit-string of all zeros in the
computational basis when a quantum computer is prepared in the state
$|a_{\delta t},\bm{\beta},\bm{\gamma}\rangle$.

The meaning of Eq. 16 (or, equivalently, Eq. 17) is demonstrated visually in
Fig. 2(b). The sphere is constructed in the same manner as Fig. 2(a) but
colored regions now represent the Hilbert space region accessible by the
embedding of $\bm{\chi}$ at time $t$ ($\mathcal{H}_{t}$; red regions) and at
time $t+q\Delta t$ ($\mathcal{H}_{t+\Delta t}$; yellow regions) where
$q\in\mathbb{Z^{+}}$ and $\Delta t\in\mathbb{R}_{>0}$ are fixed constants. The
squared modulus of overlapping regions of red and yellow on Fig. 2. (b) can be
interpreted as either Eq. 16 or 17.

##  IV Demonstrations with quantum circuit simulators

###  IV.1 Building intuition: a synthetic example

Figure 3: Quantum circuit simulator experiments with the moons2circles data
set. (a) The evolution of the time-dependent class prediction
$y_{\text{vote}}^{t}$ as moons (far left panel) gradually become circles (far
right panel). Regions predicting $y_{i}=1$ are shown in dark red and regions
predicting $y_{i}=-1$ are shown in dark blue. The Testing data set is overlaid
and are colored by their truth values as given in the legend. The figure is
generated using the best model obtained for the QAOA-3-SEL-3 runs described in
the main text. The upper panel displays results with time dependence (wtd) in
the inner product space while the lower panel displays results without time
dependence in the inner product space (w/otd) (b-e) The normalized loss
$\mathcal{L}_{\bm{\phi_{\text{min}}}}(\bm{\theta})/N^{2}_{\text{batch}}$ as a
function of mini-batch iterations for the $R_{y}$-SEL-1, $R_{y}$-SEL-3,
QAOA-1-SEL-1 and QAOA-3-SEL-3 quantum circuits, respectively. The mean value
of 50 independent runs is shown in red and the probability density over these
runs is shaded in grey. Recall that according to Eq. 10 we are solving a
maximization problem, so we expect the loss in increase with the number of
mini-batch iterations. (f-i) For the same quantum circuits and order defined
in (b-e), the time-resolved embedding overlap
$\mathcal{F}_{\bm{\beta}^{\star},\bm{\gamma}^{\star}}(\delta t)$ is plotted
using the best model obtained using each circuit. Green markers are inserted
at $\mathcal{F}_{\bm{\beta}^{\star},\bm{\gamma}^{\star}}(t_{l})$ to show
overlaps of embedding spaces where classical data was seen in training.

As a first example of our time-series classification algorithm, we deal with a
didactic case to build intuition. Our goal is to create a synthetic time-
series classification problem where (i) data is not linearly separable (thus
kernelized SVM is required) and (ii) the shape and position of the decision
boundary must evolve with time in order to successfully classify data at
different values of $t$. Accordingly, we generate a training set
$X_{\text{syn}}$ with $|X_{\text{syn}}|=100$ two-dimensional ($d=2$) time-
series instances. We construct the beginning and endpoints of each time series
using the well-known moons and circles data sets, respectively. Intermediate
points in the series are generated using a 10-step (i.e. $p=10$) linear
interpolation between the start and end points. The result is a gradual
transition from moons to circles which we call the moons2circles data set.
Testing data from moons2circles is shown as markers overlaid on the panels of
Fig. 3(a).

We use two different varieties of $n=3$ qubit quantum circuit structures to
implement the quantum components of the QCC-net . We denote first variety
$R_{y}$-SEL-$g$ and the second QAOA-$g$-SEL-$g$. We implement $g=1$ and $3$
varieties where $g$ denotes the number of repeated layers of the layered
Ansätze. $R_{y}$ means a fixed embedding circuit (i.e. $\bm{\alpha}$ is
dropped) of $U(\bm{x}_{t})=R_{y}({}^{1}x_{t})\otimes R_{y}({}^{2}x_{t})\otimes
I$ is used and QAOA means an embedding circuit inspired by the Quantum
Alternating Operator Ansätz [51, 52, 53] is used to implement
$U(\bm{x}_{t},\bm{\alpha})$. SEL means a Strongly Entangling Layers Ansätz
[54] is used to implement $W(\bm{\beta})$. For all circuits,
$D(\bm{\gamma},t)$ is implemented using a circuit representing a truncated
n-local Walsh-operator expansion as described in [55]. The kernel function of
Eq. 6 is then trained using the process described in Sec. II.2 for 250 mini-
batch iterations with $N_{\text{batch}}=4$ using the Adam optimizer [56] with
an initial learning rate of 0.05. For solving the convex problem (i.e, finding
the Lagrangian dual variables $\bm{\phi}$) we use the splitting conic solver
[57, 58] which interfaces cvxpyayers [44] for gradient computations [45] which
are passed to the automatic differentiation in Pennylane [41]. The training is
repeated 50 times where each run is seeded with initial parameters
$\bm{\theta}_{\text{init}}$ drawn uniformly at random from the range
$[-\pi,\pi]$.

The results of training are shown in Fig. 3(b-e). On average we can see that
QAOA-3-SEL-3 produces the largest loss values both initially (i.e. at the
first mini-batch iteration) and after training for the 250 mini-batch
iterations. While it is intuitive that the loss after training should be
largest for QAOA-3-SEL-3 simply because it has the largest number of quantum
circuit parameters $\bm{\theta}$, this does not explain the large initial
loss. Indeed, at the first mini-batch iteration, we measure the separation of
binary classes for optimally combined random kernels. This means that combined
random kernels become better separators as the depth and number of parameters
increases for the moons2circles data set. This effect is so pronounced that
untrained QAOA-3-SEL-3 models, on average, have larger losses than trained
$R_{y}$-SEL-1 and $R_{y}$-SEL-3 models. Examining the density heat maps [grey
shaded areas on Fig. 3(b-e)], we can see that training is highly stochastic,
much of which is owed to the small mini-batch size. We do see this
stochasticity reduce for the QAOA-$g$-SEL-$g$ models where it is more obvious
that an envelope of high density forms around the mean curve showed in dark
red.

Using the best (as determined by the largest value of the loss evaluated using
the entire testing data set) trained TSHK with the QAOA-3-SEL-3 Ansätz, we
train a SVM with $C=100$ at each of the 10 time steps and compute
$y_{\text{vote}}^{t}(\bm{x}_{t})$ on a fine $100\times 100$ grid at each $t$.
These results are shown on the upper panel of Fig. 3(a) as heat-maps at each
point in time. Regions colored in dark red/dark blue show areas where
$y_{\text{vote}}^{t}=1$/ $y_{\text{vote}}^{t}=-1$ The result is clear: the
decision function is time-evolving and adapts to separate data belonging to
the two different classes at every time step. Although this is a simple
example, we note that the best TSHK for any circuit structure variation
achieves 100% classification accuracy when passed to a SVM. This is despite
the data at intermediate time steps not being perfectly separated (data of
different classes overlap). Upon inspection of the kernel weights, we see that
the beginning and end-points are most heavily weighted which is intuitive
since these points see no overlap between classes. This observation exposes
the power of our algorithm: time-series can be classified based on “smoking-
gun” time points as these points are able to be heavily weighted by the convex
optimizer.

The effects of including a time-dependent inner product space become even
clearer should we observe the performance on the moons2circles data set in its
absence. This is shown in the lower panel of Fig. 3(a) as achieved by fixing
$t=1$ (an arbitrary choice) in the time evolution operator. Other settings are
identical to those which produce the results of the upper panel of Fig. 3(a).
While there remains some change in the position of the decision boundary as we
advance in time, it is certainly much less malleable than in the case with the
time dependent inner product. Now unable to use a different kernel function at
each $t$, training finds a single kernel function which strikes a balance
between being able to classify moons-like and circles-like data, which, when
inserted into a SVM, results in decision boundaries with a striped character.
Moreover, classification accuracy with this model slips down to 78%. From
these investigations, we conclude that using the TSHK is most effective for
data sets where (i) the correct notion of similarity between pairs of data
points changes with time and (ii) there exists a shared time axis in the data
as was described in Sec. II.1.

In Fig. 3(f-g) we probe the time dependence of the best trained models using
the approach presented in Sec. III. With the exception of $R_{y}$-SEL-3 [Fig.
3(g)] $\mathcal{F}_{\bm{\beta}^{\star},\bm{\gamma}^{\star}}(\delta t)$ is
observed to be highly oscillatory. Importantly, apart from QAOA-1-SEL-1 [Fig.
3(h)], $\mathcal{F}_{\bm{\beta}^{\star},\bm{\gamma}^{\star}}(\delta t)$ is a
broadly decreasing function intuitively suggesting that the quantum embedding
space when data is more moons-like becomes different as the data becomes more
circles-like. Indeed, even for QAOA-SEL-1, the last point (shown with green
markers) of all models shows that
$\mathcal{F}_{\bm{\beta}^{\star},\bm{\gamma}^{\star}}(\delta t)$ commensurate
with the time difference between exactly moons [the left-most panel of Fig.
3(a)] and exactly circles [the right-most panel of Fig. 3(a)] is very small
which means that the quantum embedding space used for classifying moons is
strongly different than that used for classifying circles.

###  IV.2 Univariate time-series: the gun-point data set

In this section, we demonstrate the performance of our algorithm for
classification on a real and popular benchmark: the gun-point data set [32].
Although this data set has been described in detail elsewhere, we summarize
briefly here. Two participants, with their hands beginning at their sides, are
asked to either (i) point their right finger towards a target or (ii) point a
small firearm (a hand gun) at a target and bring their hand back down to the
starting position. During this process, the forward motion on their hand is
tracked with a sensor at regular intervals and each finger-point or gun-point
instance is recorded as a univariate ($d=1$) time-series instance. Labelled
gun-point instances belong to class 1 ($y_{i}=1$) and finger-point instances
belong to class 2 ($y_{i}=-1$). The training set $X_{tr}$ has size
$|X_{tr}|=50$ and the testing set $X_{te}$ has size $|X_{te}|=150$. Each
instance has 150 time stamps ($p=150$).

Figure 4: Insights and performance metrics for trained TSHKs passed to SVMs on
the gun-point data set. (a) Testing data set time-series instances reported to
be the most class 1 (gun-point; orange dotted and dashed line) and most class
2 (finger point; blue dotted and dashed line) as determined by the value of
$\mathcal{D}^{\star}(\bm{x})$. (b) The kernels weights $\eta_{t}$ at each time
for the $n=2$ (gold line) and $n=8$ (blue line). (c) $F_{1}$, $A_{B}$, and RA
scores (left axis) and the normalized optimal loss
$\mathcal{L}^{\star}/N_{batch}^{2}$ (right axis) for the $n=2,4,6$ and $8$
models.

Using the QAOA-3-SEL-3 circuit structure described in Sec. IV.1, we train four
models with $n=2,4,6$ and $8$ qubits. Each model is selected from 20 training
attempts each using 500 mini-batch iterations with $N_{\text{batch}}=4$ with
$\bm{\theta}_{\text{init}}$ drawn uniformly at random from $[-\pi,\pi]$. Each
model is selected based on the largest value of the loss function evaluated on
the testing data set. Before the SVMs are trained, given $\bm{\theta}^{\star}$
each optimal model is passed through the convex optimization step a final time
using the entire training data set (i.e, $N_{\text{batch}}=N_{X}$) to refine
the kernel weights $\bm{\eta}$. We now train $C=100$ SVM decision functions
$\mathcal{D}^{\star}(\bm{x})$ which we use to classify each time-series
instance in $X_{\text{test}}$. With these predictions, we calculate the
$F_{1}$ score, the balanced accuracy score ($A_{B}$) and the receiver
operating characteristic area under the curve score (ROC AUC; abbreviated
further to RA for the rest of this work).

Fig. 4 shows the results of these experiments. In Fig. 4(a), The time-series
instances which were identified as the most class 1 and the most class 2 as
determined by having the highest and lowest values of
$\mathcal{D}^{\star}(\bm{x})$ are shown. Directly below on the same time axes,
in Fig. 4(b), we have the time-resolved kernel weights $\eta_{t}$ for the
$n=2$ and $n=8$ models. It is clear that different points in time are weighted
very differently to others. In particular, peaks near the regions where the
finger/gun is being lifted and peaks near where the finger/gun is being put
back down again are seen. This is in stark contrast to the beginning, end and
middle points of the series where finger/gun is approximately stationary. This
points to the key discriminating factors between the two cases being a
reaction time difference and/or other characteristic differences during the
time stamps where the finger/gun is in motion. It should also be noted that
there is a remarkable symmetry in the profile of $\eta_{t}$ about the mid-
point of time (where the finger/gun is stationary); finger/gun instances can
be discriminated equally well from motion on the way up to pointing at the
target as they can on the way down. All of the models produce a similar
profile for $\eta_{t}$, which, for $n=2$ and $n=8$ are shown in Fig. 4(a). It
can be seen that increases in the qubit size for the model lead only to small
changes in the weighting of different time stamps in the classification.
Looking now at Fig. 4(c), we see that the balanced accuracy score is already
at 97.2% for the $n=2$ model, and raises a small amount to $98.0\%$ by the
time we reach $n=8$. Similar small gains are seen for $F_{1}$ scores and the
ROC AUC score. Remarkably, these scores mean that our hybrid quantum-classical
algorithm is competitive with purely classical algorithms for this data set.
Specifically, when compared against 9 state-of-the-art deep learning
approaches [1], our algorithm is beaten in terms of accuracy scores by only
two of them: fully convolutional neural networks (FCN; $A_{B}=1.000$) and
residual neural networks (ResNet; $A_{B}=0.991$) which, in comparison, require
far greater computational resources.

Figure 5: Discerning the impact of different mechanisms defining the TSHK as
applied to SVM classification of the gun-point data set. Different approaches
include the full TSHK method (red bars), the otherwise full method with time
dependence in the inner product space forbidden (a single static quantum
kernel is used at each $t$; grey bars), the otherwise full method but with
randomly chosen quantum circuit parameters $\bm{\theta}$ (dark blue bars) and
an approach using random $\bm{\theta}$, equal kernel weights ($\eta_{t}=1/p$
$\forall t$) and a time-dependent inner product space.
$\mathcal{A}_{train/test}=1-1/p^{2}\sum_{i}^{p}\sum_{j}^{p}|K_{train/test,ij}-K_{truth,ij}|$
is the kernel alignment of the training/testing data sets, $K_{train/test,ij}$
are the elements of combined kernels for the training/testing datasets and
$K_{truth,ij}=y_{i}y_{j}$ are the elements of the truth matrix. (a)
$N_{batch}=4$ and $\lambda=0.1$. (b) $N_{batch}=50$ (the entire training data
set) and $\lambda=0.8$. An expanded version of this figure is available in the
Supplemental Material [59].

To quantify which components of our approach are most significant in our
experiments and what training hyperparameters yield the best results, we
conduct a further study focusing on the $n=2$ qubit case. A summary of the
results are shown in Fig 5 and an expanded version of Fig. 5 is available in
the Supplemental Material [59]. In general, we find that training with small
batch sizes and small $\lambda$ is the most effective (Fig. 5(a)). This effect
is enhanced when paired with the full TSHK formalism (red bars) which, for
these hyperparameters, yields the largest $A_{B}$ and $F_{1}$ scores as well
as the largest training kernel alignment $\mathcal{A}_{train}$ (see caption of
Fig. 5 for a definition) of all approaches taken in this study. This peak
performance must, however, be viewed in the context of the other bars on Fig.
5(a). Examining the case where the time dependent inner product is forbidden
(grey bars), the full approach offers only a 1% improvement in $A_{B}$ and
$F_{1}$. Indeed, examining the general level of scores for all other bars on
both Fig. 5(a) and 5(b), we can see that any of the taken approaches can score
in the mid-90% range. It must therefore be concluded that most of the
performance we achieve on this data set is gained from the classical portions
of the algorithm (mostly the SVM itself) correctly classifying many instances
using any of the combined kernels we supply it. Only a few extra points are
gained by using the full TSHK formalism. Given how closely performant the top
classical deep learning approaches are for this data set [1], it is, however,
not unreasonable to suggest that these gains are indeed significant.
Furthermore, we view the results presented Fig. 4 to be sensitive to the
choice of circuit Ansätz. It is unclear what choice should be made when
presented with classical data and is a topic of current research to figure out
how to encode inductive biases into classical data encoding Ansätze [60].

Lastly, we remark that while scores only increase modestly as the number of
qubits increases [Fig. 4(c)], the loss does gradually get larger. Recalling
that this loss represents the separation between positive and negative
classes, this means that, within the parameters of our study, higher qubit
models are able to find a time-dependent quantum embedding space where
classical data is better separated than smaller qubit models. This larger
separation is known to give rise to higher generalization performance for SVM
[61].

##  V Parallel execution on quantum hardware

###  V.1 Quantum multi-programming

To more efficiently use NISQ hardware, it is possible to run multiple quantum
circuits concurrently. Known as Quantum Multi-Programming (QMP), this
technique allows NISQ devices to execute several quantum circuits
concurrently, even if they differ in structure or complexity. In this Section,
we show that computation of TSHKs can be significantly sped up using QMP.

The primary motivation for QMP is driven by the fact that the number of qubits
present in NISQ devices is often significantly higher than their Quantum
Volume (QV). Specifically, when compared to trapped ion devices,
superconducting quantum computers have restricted qubit connectivity and lower
QV. Taking this further, following IBM’s quantum volume $V_{Q}$ we have
$\displaystyle\log_{2}V_{Q}=\text{argmax}_{m}\left\\{\min\left[m,d(m)\right]\right\\}$
where $d(m)$ is the depth of a model circuit with $m$-qubits [62]. One depth
of the model circuit is composed of a random permutation of the qubits
involved in the test, followed by random two-qubit gates. For instance, the QV
of the ibm_washington device (used in this work) is $64$, with $127$ qubits
available. This means the model circuit with $6$-qubits runs reliably at six
depths of the model circuit on average on ibm_washington. Other
superconducting quantum computers share similar properties such as limited
connectivity between qubits, and many more qubits than $\log_{2}V_{Q}$.
Therefore, it is crucial to use the capacity of modern superconducting quantum
devices more efficiently, and QMP has been proposed as a method to achieve
this goal. Therefore, QMP fills the gap between the relatively many qubits and
the relatively low QV of NISQ devices by executing multiple quantum circuits
concurrently, enhancing the throughput and utilization of NISQ devices.

However, implementing QMP on NISQ devices has several issues to address
because QMP accompanies unfavorable impacts on the whole system such as
measurement timing of the concurrent circuits [31] and crosstalk between
different circuits [63]. Efficiently mapping qubits between logical and
physical states, as well as task scheduling [64, 65, 66], have been studied
alongside the aforementioned issues. Despite preliminary investigations into
QMP, its integration with quantum algorithms, and especially QML, remains
underdeveloped. Nonetheless, there has been recent noteworthy progress in
applying QMP to Grover’s search algorithm [67], resulting in improved success
probabilities compared to previous attempts.

###  V.2 QMP implementation considerations

Before proceeding to calculate TSHKs with QMP, we discuss how the typical
issues arising out of QMP. Starting with crosstalk, this effect is studied in
detail by P. Murali et al [68]. They suggested several rules to mitigate
crosstalk between qubits. The relation between the number of physical buffers
and the error rate was studied by Ohkura et al [63]. They introduced a
physical buffer, which is the number of idle qubits between quantum circuits.
They tested a different number of controlled-X (CX) gates with a different
number of physical buffers. Ohkura et al. [63] concluded that only $1$
physical buffer is sufficient until $30$ CX gates. Therefore, we use one
physical buffer for our implementations.

One important issue of QMP is measurement timing. When quantum circuits having
different circuit depths run concurrently, the measurement timing influences
the shortest circuit. This effect was studied in Ref. [31, 69]. In these
studies, they suggested delaying shorter circuits to align the measurement
timing. In our case, because the circuits implementing elements of the TSHK
each have the same depth, we do not need to make any alterations to
measurement timings.

Another crucial issue is efficient physical qubit mapping. Systematic qubit
mapping algorithms have been developed [31, 64, 66, 65]. Even though these
qubit mapping algorithms are also important for quantum circuits executed in
serial (i.e. without QMP), efficient qubit mapping becomes more challenging
and crucial in QMP because practical QMP circuits use more qubits than serial
circuits. In our particular circumstance where each quantum thread is composed
of two logical qubits, the mapping problem is solved by hand by mapping each
thread to two physical qubits with two-qubit gate connectivity subject to the
constraint that at least one buffer qubit separates individual threads. The
QMP layout can be seen in Fig. 6 for ibm_washington and ibm_sherbrooke.

Finally, we mention partial measurement in QMP. After the whole measurement of
a QMP circuit, we need to extract the measurements from each parallelized
circuit. This is called partial measurement and is implemented by post-
processing after the whole measurement since the partial measurement
probability is the sum of all other uninvolved qubits measurement probability.
This is detailed in Ref. [67] and in the Supplemental Material [59]. Since the
example code in Appendix A in Ref. [67] iterates over the $2^{n}$ where $n$ is
the number of the measured qubits, the computing time increases exponentially
as more qubits are measured. Hence, we re-implement the partial measurement
code, which extracts the partial measurement out of the whole measurement, to
depend on the number of shots regardless of the number of the measured qubits.
Our new implementation is given in the Supplemental Material [59].

Figure 6: The mapping between logical and physical qubits for the QMP circuits
used in our experiments on (a)ibm_washington (Eagle r1 processor) and (b)
ibm_sherbrooke (Eagle r3 processor). Both machines have 127 physical qubits
and have a similar (but not identical) qubit connectivity graph.
ibm_washington has a QV of 64 while ibm_sherbrooke has a QV of 32. On
ibm_washington, the native gates are CX, ID, RZ, SX and X and ibm_sherbrooke
uses the same single qubit gates but ECR gates are used for 2-qubit
interactions. Idle qubits are shown as blue circles and active qubits are
black circles containing integers indexing the logical qubits.

###  V.3 QMP experiments with the gun-point data set

We conduct QMP experiments using the gun-point data set described in Sec. IV.2
where the time dimension is decimated by $3\times$ ($p$ is reduced from 150 to
50) to reduce the total QPU run time to a reasonable level. An optimal
simulated (ideal) model for the decimated data set is first obtained using the
same method described in Sec. IV.2. Using the $\bm{\theta}^{\star}$ derived
from the simulated model, we generate the quantum circuits used to calculate
the matrix elements of the kernel matrices required to train an SVM. These
circuits are deployed on two 127-qubit superconducting QPUs: ibm_washington
(64 QV) and ibm_sherbrooke (32 QV). The experiments utilize the maximum number
of active qubits which can be used in a single QMP run. For both devices,
including the buffer qubit constraint, this number is 70. Fig. 6 shows the
qubit layout on ibm_washington [Fig. 6(a)] and ibm_sherbrooke [Fig. 6(b)] . In
addition to mapping considerations already discussed in Sec. V.2, we manually
choose the qubits with low error ratios (read assignment error, etc.).
Considering each circuit has two qubits, we can run 35 circuits in parallel
with a single QMP run. The Trial Reduction Factor (TRF) [31] describes the
ratio of the number of trials (shots) that executes individually (baseline) to
the number of trials (shots) of the QMP circuit. This factor represents the
efficiency of the QMP in terms of trials (shots). In our case, the TRF is
$35$. The QMP used $8192$ shot numbers to measure $35$ circuits whereas non-
QMP implementations need $8192\times 35$ shots to measure $35$ circuits.

Recalling the total size of the training and testing data sets in the
decimated gun-point data set, a naive serial QPU implementation (without QMP)
requires $436,250$ individual calls to a quantum device to calculate all of
the relevant kernel matrices to train and test an SVM. By applying QMP, the
total number of calls to the QPU is significantly reduced by a factor of 35 to
$12,465$ because the QMP circuit parallelizes $35$ quantum circuits. As a
preliminary step, we first verified that our QMP implementation produced
consistent outputs with serial execution by comparing the output of both using
a quantum circuit simulator. The outputs for serial and parallel are in
perfect agreement. The details of these tests, and others, are given in the
Supplemental Material [59]. Importantly, in practice, QMP reduces the queuing
time of the circuit. Because of the much-reduced number of calls to the QPU,
there are far fewer circuits in the queuing system of an individual device.
The advantage of this is two-fold as (i) the wait time is reduced for the user
and other users and (ii) schedulers often use a “fair-share” allowance on the
number of jobs submitted by a given user. The fewer number of circuits sent
means that job priority levels administered by the scheduler remain high.
Finally, we note that because of our choice to use two separate 127 qubit
machines, circuits were sent to both in parallel thus, two queues were being
occupied at all times which itself leads to a practical speedup.

Figure 7: Simulator versus device results on the sub-sampled gun-point data
set. (a) The $F_{1}$, $A_{B}$ (balanced accuracy) and RA (ROC AUC) scores for
the simulator (red) and device (grey; QPU $\times 2$) results. (b) The
percentage error in the kernel weights from the device results vs the
simulator given as
$\bm{\eta}_{error}=100\cdot(\bm{\eta}_{device}-\bm{\eta}_{ideal})/\bm{\eta}_{ideal}$
plotted as a function of time index $l$. (c-d) The SVM decision function
$\mathcal{D}^{\star}(\bm{x})$ for each of the 150 testing time-series using
the simulated kernel matrices (c) or the device kernel matrices (d). Points
are colored by their truth label. Blue points above the black dotted line are
correctly classified while those below it are incorrectly classified. The
reverse is true for gold points.

The results of our QMP experiments are shown in Fig. 7. Fig. 7(a) compares
$F_{1}$, $A_{B}$ and RA scores for the simulated results (red bars) and the
QMP runs executed on ibm_washington and ibm_sherbrooke (QPU $\times 2$; grey
bars). Compared with the simulator results from Sec. IV.2, we see that MKL SVM
on decimated gun-point data set gives rise to lower scores for all three
metrics. This is intuitive as the algorithm is being shown fewer data points
(100 fewer time points) in training. When run on $2\times$ QPU, scores fall
$\approx$ 10-20% across the board as a result of all of the various types of
hardware noise. We should note that besides the use of buffer qubits and
selection of low-error qubits on the machines, no other optimizations/error
mitigation techniques were used so the results presented in Fig. 7 can be
considered a baseline. In reality, more highly optimized transpiler passes
could be used alongside a whole host of dedicated error mitigation techniques
including (but not limited to) Pauli twirling, dynamical decoupling and
various flavors of read-out mitigation. Fig. 7(b) shows the error in the
kernel weights $\bm{\eta}$ as a function of time (i.e. $\eta_{t}$ is shown).
Errors exceeding 25% occur at low and high $t$ while errors are very small at
intermediate $t$. Interestingly (although not shown), areas where the weights
themselves are larger (at intermediate $t$) have the smallest errors which
means despite a significant amount of noise entering our kernel matrix
computations, the kernel weight values are quite robust. Fig. 7(c) and 7(d)
show the SVM decision function $\mathcal{D}^{\star}(\bm{x})$ (discussed in
Sec. II.3) evaluated using all 150 of the testing time-series instances using
where the TSHK was either computed with an ideal quantum circuit simulator
[Fig. 7(c)] or calculated on real devices using QMP [Fig. 7 (d)]. Recall that
$\text{sgn}[\mathcal{D}^{\star}(\bm{x})]$ distinguishes the predicted binary
class label such that points above the black dotted line at
$\mathcal{D}^{\star}(\bm{x})=0$ are predicted as class 1 while those below it
are predicted as class -1. Because each point is colored by its truth label
$y_{i}$, it is possible to see from Fig. 7(c) and 7(d) which points are
correctly and incorrectly classified. Looking first at the simulator results
[Fig. 7(c) ] it can be seen that many points sit very close to
$\mathcal{D}^{\star}(\bm{x})=0$ meaning their classification is a close call.
This finding is important as it provides some insight into why performance
metrics decrease for the QPU results. That is, we can see from Fig. 7(d) that
those close call points which were sitting close to
$\mathcal{D}^{\star}(\bm{x})=0$ for the simulated results have been smeared
across the decision boundary by a non-trivial propagation of error incurred by
(i) training the SVM with a noisy (and only approximately positive semi-
definite) Gram matrix and (ii) making predictions using already noisy trained
SVM model using a similarly noisy kernel function. We observe also that those
points which were sitting further from the decision boundary in the simulator
results now even further from it.

Following the approach taken in Ref. [70], errors incurred by the processes
described in points (i) and (ii) above can be mitigated against. Using the
Tikhonov regularization scheme [70], Gram matrices $K$ can regularized (be
made positive semi-definite) with a simple correction:

$K^{\text{reg}}=\begin{cases}K-\epsilon_{\text{min}}I&\text{if
}\epsilon_{\text{min}}<0\\\ K&\text{otherwise}\end{cases}$ (18)

for regularized gram matrix $K^{\text{reg}}$, smallest eigenvalue in the
spectrum of $K$, $\epsilon_{\text{min}}$ and identity matrix $I$. We apply the
correction of Eq. 18 to each of the noisy time-dependent Gram matrices $K_{t}$
and use the regularized output to find the optimal kernel weights
$\eta_{t}^{\star}$ given in Eq. 11. Finally, we build a regularized combined
kernel using Eq. 6. We now use the regularized combined kernel as input to a
SVM to classify the gun-point data set. We find that all considered metrics
are improved: $A_{B}=0.81$, $F_{1}=0.82$ and $RA=0.90$. This shows the success
of even simple regularization techniques and supports the need for future
works assessing the best possible way to use kernel methods on noisy QPUs.

##  VI Conclusions and Perspectives

In this work, we proposed and demonstrated a hybrid quantum-classical
supervised ML algorithm for time-series classification. The algorithm works by
utilizing time-dependent inner product spaces as generated by a common time-
evolution operator and combining the inner products using classical MKL
techniques. Together, the components of the algorithm create a temporally-
aware QCC-net able to tailor kernel functions for arbitrary time-series data.
When these tailored kernel functions were used with kernelized SVM, we found
classification performance to be comparable with purely classical methods for
the well-known gun-point data set, although, a great deal of the performance
stemmed from the classical portions of the algorithm. Through conducting
experiments on the synthetic moons2circles data set, we verified that the TSHK
is most effective in situations where the correct notion of similarity between
pairs of data points evolves with time and there exists a shared time axis in
the data as was described in Sec. II.1.

Furthermore, we developed a method to study the resulting time-dependence
present within a trained time-series kernel function. We must stress that
algorithms trained to explicitly learn time-dependent inner product spaces are
novel even within the classical ML literature and indeed, we do not know of an
existing classical approach that can precisely emulate the method proposed in
this work.

We also found that the total run time of our method can be greatly reduced by
utilizing the state-of-the-art method, QMP. That is, because the computation
of quantum kernel matrix elements is an embarrassingly parallel problem, many
matrix elements can be computed in parallel by distributing multiple quantum
threads across large QPUs. We demonstrated this parallelism using two
127-qubit superconducting quantum processors to perform SVM classification on
a decimated version of the gun-point data set. With this method, we were able
to reduce the total number of required QPU calls by $35\times$, which,
assuming constant job queuing times for both 127-qubit QPUs, leads to a
$70\times$ speed-up compared to serial execution on a single QPU. We have
therefore demonstrated, for the first time, that QMP is a valuable tool for
significantly speeding up QML algorithms on present and near future NISQ
machines. With recent advancements in low-latency CPU-QPU interactions by
Qiskit Runtime and others, in the near future, we expect to see the emergence
of real-time training of error-mitigated QML where the quantum components of
the algorithm are accelerated by QMP and the classical parts are accelerated
using conventional parallel processing paradigms like MPI and OpenMP. This
work, therefore, makes strides along the path toward the much-anticipated
intersection of quantum computation and classical high performance computing.

## Acknowledgments

We acknowledge the use of IBM Quantum services for this work. The views
expressed are those of the authors and do not reflect the official policy or
position of IBM or the IBM Quantum team.

## Data and code availability

The data used to produce Fig. 3, 4 and 7 can be found at
https://doi.org/10.5281/zenodo.7996534 alongside a code example implementing a
QCC-net for the Sine vs Cosine classification shown in Fig. 1.

## Competing interests

The authors declare no competing personal or financial interests.

## Funding

This research used quantum computing resources of the Oak Ridge Leadership
Computing Facility, which is a DOE Office of Science User Facility supported
under Contract DE-AC05-00OR22725. This research used resources of the National
Energy Research Scientific Computing Center, a DOE Office of Science User
Facility supported by the Office of Science of the U.S. Department of Energy
under Contract No. DE-AC02-05CH11231 using NERSC award DDR-ERCAP0024165.

## Author’s contributions

Jack S. Baker made the figures, conducted the numerical experiments on quantum
circuit simulators, wrote the manuscript and helped develop the theoretical
method. Gilchan Park conducted the QMP experiments and contributed to the
manuscript. Kwangmin Yu oversaw the QML experiments, contributed to the
manuscript and the supplemental information. Ara Ghukasayan helped to develop
the theoretical method and reviewed the manuscript. Oktay Goktas oversaw the
numerical experiments and reviewed the manuscript. Santosh Kumar Radha seeded
the idea of using time evolution operators to create a time dependent quantum
kernel function and reviewed the manuscript.

## References

  * Fawaz  _et al._ [2019] H. I. Fawaz, G. Forestier, J. Weber, L. Idoumghar, and P.-A. Muller, Deep learning for time series classification: a review, Data Mining and Knowledge Discovery 33, 917 (2019).
  * Blázquez-García  _et al._ [2021] A. Blázquez-García, A. Conde, U. Mori, and J. A. Lozano, A review on outlier/anomaly detection in time series data, ACM Computing Surveys 54, 1 (2021).
  * Choi  _et al._ [2021] K. Choi, J. Yi, C. Park, and S. Yoon, Deep learning for anomaly detection in time-series data: Review, analysis, and guidelines, IEEE Access 9, 120043 (2021).
  * Clark  _et al._ [2020] S. Clark, R. J. Hyndman, D. Pagendam, and L. M. Ryan, Modern strategies for time series regression, International Statistical Review 88, 10.1111/insr.12432 (2020).
  * Deb  _et al._ [2017] C. Deb, F. Zhang, J. Yang, S. E. Lee, and K. W. Shah, A review on time series forecasting techniques for building energy consumption, Renewable and Sustainable Energy Reviews 74, 902 (2017).
  * Torres  _et al._ [2021] J. F. Torres, D. Hadjout, A. Sebaa, F. Martínez-Álvarez, and A. Troncoso, Deep learning for time series forecasting: A survey, Big Data 9, 3 (2021).
  * Zhang  _et al._ [2018] C. Zhang, S. R. Kuppannagari, R. Kannan, and V. K. Prasanna, Generative adversarial network for synthetic time series data generation in smart grids, in  _2018 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm)_ (IEEE, 2018).
  * Hochreiter and Schmidhuber [1997] S. Hochreiter and J. Schmidhuber, Long short-term memory, Neural Computation 9, 1735 (1997).
  * Melis  _et al._ [2019] G. Melis, T. Kočiský, and P. Blunsom, Mogrifier lstm (2019).
  * Nguyen  _et al._ [2020] T. Nguyen, R. Baraniuk, A. Bertozzi, S. Osher, and B. Wang, Momentumrnn: Integrating momentum into recurrent neural networks, in  _Advances in Neural Information Processing Systems_ , Vol. 33, edited by H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin (Curran Associates, Inc., 2020) pp. 1924–1936.
  * Vaswani  _et al._ [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin, Attention is all you need, in  _Advances in Neural Information Processing Systems_ , Vol. 30, edited by I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Curran Associates, Inc., 2017).
  * Yang  _et al._ [2021] C.-H. H. Yang, Y.-Y. Tsai, and P.-Y. Chen, Voice2series: Reprogramming acoustic models for time series classification, in  _International Conference on Machine Learning_ (PMLR, 2021) pp. 11808–11819.
  * Zerveas  _et al._ [2021] G. Zerveas, S. Jayaraman, D. Patel, A. Bhamidipaty, and C. Eickhoff, A transformer-based framework for multivariate time series representation learning, in  _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery; Data Mining_ (ACM, 2021).
  * Cortes and Vapnik [1995] C. Cortes and V. Vapnik, Support-vector networks, Machine learning 20, 273 (1995).
  * Badiane  _et al._ [2018] M. Badiane, M. O’Reilly, and P. Cunningham, Kernel methods for time series classification and regression., in  _AICS_ (2018) pp. 54–65.
  * Bailly [2018] A. Bailly,  _Time series classification algorithms with applications in remote sensing_ , Ph.D. thesis, Université Rennes 2 (2018).
  * Fábregues de los Santos [2017] L. Fábregues de los Santos,  _Forecasting financial time series using multiple Kernel Learning_ , Master’s thesis, Universitat Politècnica de Catalunya (2017).
  * Rüping [2001] S. Rüping,  _SVM kernels for time series analysis_ , Tech. Rep. (Technical report, 2001).
  * Tino [2020] P. Tino, Dynamical systems as temporal feature spaces, The Journal of Machine Learning Research 21, 1649 (2020).
  * Gönen and Alpaydın [2011a] M. Gönen and E. Alpaydın, Multiple kernel learning algorithms, The Journal of Machine Learning Research 12, 2211 (2011a).
  * Aiolli and Donini [2015] F. Aiolli and M. Donini, EasyMKL: a scalable multiple kernel learning algorithm, Neurocomputing 169, 215 (2015).
  * Ghukasyan  _et al._ [2023] A. Ghukasyan, J. S. Baker, O. Goktas, J. Carrasquilla, and S. K. Radha, Quantum-classical multiple kernel learning (2023), arXiv:2305.17707 [quant-ph] .
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Noisy intermediate-scale quantum algorithms, Rev. Mod. Phys. 94, 015004 (2022).
  * Wiebe  _et al._ [2010] N. Wiebe, D. Berry, P. Høyer, and B. C. Sanders, Higher order decompositions of ordered operator exponentials, Journal of Physics A: Mathematical and Theoretical 43, 065203 (2010).
  * Cîrstoiu  _et al._ [2020] C. Cîrstoiu, Z. Holmes, J. Iosue, L. Cincio, P. J. Coles, and A. Sornborger, Variational fast forwarding for quantum simulation beyond the coherence time, npj Quantum Information 6, 10.1038/s41534-020-00302-0 (2020).
  * Radha [2021] S. K. Radha, Quantum constraint learning for quantum approximate optimization algorithm (2021).
  * Horowitz  _et al._ [2022] H. Horowitz, P. Rao, and S. K. Radha, A quantum generative model for multi-dimensional time series using hamiltonian learning (2022).
  * Baker  _et al._ [2022] J. S. Baker, H. Horowitz, S. K. Radha, S. Fernandes, C. Jones, N. Noorani, V. Skavysh, P. Lamontangne, and B. C. Sanders, Quantum variational rewinding for time series anomaly detection (2022).
  * Gibbs  _et al._ [2022] J. Gibbs, Z. Holmes, M. C. Caro, N. Ezzell, H.-Y. Huang, L. Cincio, A. T. Sornborger, and P. J. Coles, Dynamical simulation via quantum machine learning with provable generalization (2022).
  * Caro  _et al._ [2022] M. C. Caro, H.-Y. Huang, N. Ezzell, J. Gibbs, A. T. Sornborger, L. Cincio, P. J. Coles, and Z. Holmes, Out-of-distribution generalization for learning quantum dynamics (2022).
  * Das  _et al._ [2019] P. Das, S. S. Tannu, P. J. Nair, and M. Qureshi, A case for multi-programming quantum computers, in  _Proceedings of the 52nd Annual IEEE/ACM International Symposium on Microarchitecture_ (2019) pp. 291–303.
  * Ratanamahatana and Keogh [2005] C. A. Ratanamahatana and E. Keogh, Three myths about dynamic time warping data mining, in  _Proceedings of the 2005 SIAM international conference on data mining_ (SIAM, 2005) pp. 506–510.
  * Gönen and Alpaydın [2011b] M. Gönen and E. Alpaydın, Multiple kernel learning algorithms, The Journal of Machine Learning Research 12, 2211 (2011b).
  * Schuld and Killoran [2019] M. Schuld and N. Killoran, Quantum machine learning in feature hilbert spaces, Physical Review Letters 122, 10.1103/physrevlett.122.040504 (2019).
  * Havlíček  _et al._ [2019] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Supervised learning with quantum-enhanced feature spaces, Nature 567, 209 (2019).
  * Stone [1932] M. H. Stone, On one-parameter unitary groups in hilbert space, The Annals of Mathematics 33, 643 (1932).
  * Bjorken and Drell [1965] J. D. Bjorken and S. D. Drell,  _Relativistic quantum fields_ (McGraw-Hill, 1965).
  * Sakoe and Chiba [1978] H. Sakoe and S. Chiba, Dynamic programming algorithm optimization for spoken word recognition, IEEE Transactions on Acoustics, Speech, and Signal Processing 26, 43 (1978).
  * Aiolli  _et al._ [2008] F. Aiolli, G. D. S. Martino, and A. Sperduti, A kernel method for the optimization of the margin distribution, in  _Artificial Neural Networks - ICANN 2008_ (Springer Berlin Heidelberg, 2008) pp. 305–314.
  * O’Donoghue  _et al._ [2016a] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd, Conic optimization via operator splitting and homogeneous self-dual embedding, Journal of Optimization Theory and Applications 169, 1042 (2016a).
  * Bergholm  _et al._ [2018] V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, S. Ahmed, V. Ajith, M. S. Alam, G. Alonso-Linaje, B. AkashNarayanan, A. Asadi, J. M. Arrazola, U. Azad, S. Banning, C. Blank, T. R. Bromley, B. A. Cordier, J. Ceroni, A. Delgado, O. Di Matteo, A. Dusko, T. Garg, D. Guala, A. Hayes, R. Hill, A. Ijaz, T. Isacsson, D. Ittah, S. Jahangiri, P. Jain, E. Jiang, A. Khandelwal, K. Kottmann, R. A. Lang, C. Lee, T. Loke, A. Lowe, K. McKiernan, J. J. Meyer, J. A. Montañez-Barrera, R. Moyard, Z. Niu, L. J. O’Riordan, S. Oud, A. Panigrahi, C.-Y. Park, D. Polatajko, N. Quesada, C. Roberts, N. Sá, I. Schoch, B. Shi, S. Shu, S. Sim, A. Singh, I. Strandberg, J. Soni, A. Száva, S. Thabet, R. A. Vargas-Hernández, T. Vincent, N. Vitucci, M. Weber, D. Wierichs, R. Wiersema, M. Willmann, V. Wong, S. Zhang, and N. Killoran, Pennylane: Automatic differentiation of hybrid quantum-classical computations (2018).
  * Mitarai  _et al._ [2018] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Quantum circuit learning, Physical Review A 98, 10.1103/physreva.98.032309 (2018).
  * Schuld  _et al._ [2019] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, and N. Killoran, Evaluating analytic gradients on quantum hardware, Physical Review A 99, 10.1103/physreva.99.032331 (2019).
  * Agrawal  _et al._ [2019a] A. Agrawal, B. Amos, S. Barratt, S. Boyd, S. Diamond, and J. Z. Kolter, Differentiable convex optimization layers, in  _Advances in Neural Information Processing Systems_ , Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).
  * Agrawal  _et al._ [2019b] A. Agrawal, S. Barratt, S. Boyd, E. Busseti, and W. M. Moursi, Differentiating through a cone program (2019b).
  * Bottou  _et al._ [2018] L. Bottou, F. E. Curtis, and J. Nocedal, Optimization methods for large-scale machine learning, SIAM Review 60, 223 (2018).
  * Keskar  _et al._ [2017] N. S. Keskar, D. Mudigere, J. Nocedal, M. Smelyanskiy, and P. T. P. Tang, On large-batch training for deep learning: Generalization gap and sharp minima (2017), arXiv:1609.04836 [cs.LG] .
  * Kübler  _et al._ [2021] J. M. Kübler, S. Buchholz, and B. Schölkopf, The inductive bias of quantum kernels (2021), arXiv:2106.03747 [quant-ph] .
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 10.1038/s41467-018-07090-4 (2018).
  * Pesah  _et al._ [2021] A. Pesah, M. Cerezo, S. Wang, T. Volkoff, A. T. Sornborger, and P. J. Coles, Absence of barren plateaus in quantum convolutional neural networks, Phys. Rev. X 11, 041011 (2021).
  * Lloyd  _et al._ [2020] S. Lloyd, M. Schuld, A. Ijaz, J. Izaac, and N. Killoran, Quantum embeddings for machine learning (2020), arXiv:2001.03622 [quant-ph] .
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm, arXiv preprint arXiv:1411.4028 (2014).
  * Hadfield  _et al._ [2019] S. Hadfield, Z. Wang, B. O’Gorman, E. G. Rieffel, D. Venturelli, and R. Biswas, From the quantum approximate optimization algorithm to a quantum alternating operator ansatz, Algorithms 12 (2019).
  * Schuld  _et al._ [2020] M. Schuld, A. Bocharov, K. M. Svore, and N. Wiebe, Circuit-centric quantum classifiers, Phys. Rev. A 101, 032308 (2020).
  * Welch  _et al._ [2014] J. Welch, D. Greenbaum, S. Mostame, and A. Aspuru-Guzik, Efficient quantum circuits for diagonal unitaries without ancillas, New Journal of Physics 16, 033040 (2014).
  * Kingma and Ba [2017] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2017), arXiv:1412.6980 [cs.LG] .
  * O’Donoghue  _et al._ [2016b] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd, Conic optimization via operator splitting and homogeneous self-dual embedding, Journal of Optimization Theory and Applications 169, 1042 (2016b).
  * O’Donoghue [2021] B. O’Donoghue, Operator splitting for a homogeneous embedding of the linear complementarity problem, SIAM Journal on Optimization 31, 1999 (2021).
  * [59] See Supplemental Material [publisher url] for prelimary studies on the impact of different components of the QCC-net as tested on the gun-point data set referenced in the main article and for further explanation of how QMP was implemented.
  * Bowles  _et al._ [2023] J. Bowles, V. J. Wright, M. Farkas, N. Killoran, and M. Schuld, Contextuality and inductive bias in quantum machine learning (2023), arXiv:2302.01365 [quant-ph] .
  * Hastie  _et al._ [2009] T. Hastie, R. Tibshirani, J. H. Friedman, and J. H. Friedman,  _The elements of statistical learning: data mining, inference, and prediction_ , Vol. 2 (Springer, 2009) p. 134.
  * Cross  _et al._ [2019] A. W. Cross, L. S. Bishop, S. Sheldon, P. D. Nation, and J. M. Gambetta, Validating quantum computers using randomized model circuits, Physical Review A 100, 032328 (2019).
  * Ohkura  _et al._ [2022] Y. Ohkura, T. Satoh, and R. Van Meter, Simultaneous execution of quantum circuits on current and near-future NISQ systems, IEEE Transactions on Quantum Engineering (2022).
  * Liu and Dou [2021] L. Liu and X. Dou, QuCloud: A new qubit mapping mechanism for multi-programming quantum computing in cloud environment, in  _2021 IEEE International Symposium on High-Performance Computer Architecture (HPCA)_ (IEEE, 2021) pp. 167–178.
  * Niu and Todri-Sanial [2021] S. Niu and A. Todri-Sanial, Enabling multi-programming mechanism for quantum computing in the NISQ era, arXiv preprint arXiv:2102.05321 (2021).
  * Niu and Todri-Sanial [2022] S. Niu and A. Todri-Sanial, How parallel circuit execution can be useful for NISQ computing?, in  _2022 Design, Automation & Test in Europe Conference & Exhibition (DATE)_ (IEEE, 2022) pp. 1065–1070.
  * Park  _et al._ [2023] G. Park, K. Zhang, K. Yu, and V. Korepin, Quantum multi-programming for Grover’s search, Quantum Information Processing 22, 54 (2023).
  * Murali  _et al._ [2020] P. Murali, D. C. McKay, M. Martonosi, and A. Javadi-Abhari, Software mitigation of crosstalk on noisy intermediate-scale quantum computers, in  _Proceedings of the Twenty-Fifth International Conference on Architectural Support for Programming Languages and Operating Systems_ (2020) pp. 1001–1016.
  * Ohkura [2021] Y. Ohkura, Crosstalk-aware NISQ Multi-programming, Faculty Policy Manage., Keio Univ., Tokyo, Japan (2021).
  * Hubregtsen  _et al._ [2022] T. Hubregtsen, D. Wierichs, E. Gil-Fuster, P.-J. H. S. Derks, P. K. Faehrmann, and J. J. Meyer, Training quantum embedding kernels on near-term quantum computers, Phys. Rev. A 106, 042431 (2022).
  * Yirka [2021] B. Yirka, IBM announces development of 127-qubit quantum processor (2021).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Moll  _et al._ [2018] N. Moll, P. Barkoutsos, L. S. Bishop, J. M. Chow, A. Cross, D. J. Egger, S. Filipp, A. Fuhrer, J. M. Gambetta, M. Ganzhorn, _et al._ , Quantum optimization using variational algorithms on near-term quantum devices, Quantum Science and Technology 3, 030503 (2018).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Variational quantum algorithms, Nature Reviews Physics 3, 625 (2021).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, A variational eigenvalue solver on a photonic quantum processor, Nature communications 5, 1 (2014).
  * Aaronson and Rall [2020] S. Aaronson and P. Rall, Quantum approximate counting, simplified, Symposium on Simplicity in Algorithms , 24–32 (2020).
  * Suzuki  _et al._ [2020] Y. Suzuki, S. Uno, R. Raymond, T. Tanaka, T. Onodera, and N. Yamamoto, Amplitude estimation without phase estimation, Quantum Information Processing 19, 75 (2020).
  * Yu  _et al._ [2020] K. Yu, H. Lim, and P. Rao, Practical numerical integration on NISQ devices, in  _Quantum Information Science, Sensing, and Computation XII_ , Vol. 11391 (SPIE, 2020) p. 1139106.
  * Grinko  _et al._ [2021] D. Grinko, J. Gacon, C. Zoufal, and S. Woerner, Iterative quantum amplitude estimation, npj Quantum Information 7, 1 (2021).
  * Rao  _et al._ [2020] P. Rao, K. Yu, H. Lim, D. Jin, and D. Choi, Quantum amplitude estimation algorithms on IBM quantum devices, in  _Quantum Communications and Quantum Imaging XVIII_ , Vol. 11507 (SPIE, 2020) pp. 49–60.
  * Zhang and Korepin [2020] K. Zhang and V. E. Korepin, Depth optimization of quantum search algorithms beyond Grover’s algorithm, Physical Review A 101, 032346 (2020).
  * Zhang  _et al._ [2021] K. Zhang, P. Rao, K. Yu, H. Lim, and V. Korepin, Implementation of efficient quantum search algorithms on NISQ computers, Quantum Information Processing 20, 1 (2021).
  * Zhang  _et al._ [2022] K. Zhang, K. Yu, and V. Korepin, Quantum search on noisy intermediate-scale quantum devices, Europhysics Letters 140, 18002 (2022).
  * Lubinski  _et al._ [2021] T. Lubinski, S. Johri, P. Varosy, J. Coleman, L. Zhao, J. Necaise, C. H. Baldwin, K. Mayer, and T. Proctor, Application-oriented performance benchmarks for quantum computing, arXiv preprint arXiv:2110.03137 (2021).

Supplemental Material: Parallel hybrid quantum-classical machine learning for
kernelized time-series classification

## Overview

This document contains supplemental material pertaining to the article
”Massively parallel hybrid quantum-classical machine learning for kernelized
time-series classification”. We begin by presenting a preliminary study
regarding the gun-point data set, examining the effects of hyperparameter
tuning on model performance and the effects of turning on or off optimization
over the parameters of the quantum kernels and the coefficients of the kernels
in the multiple kernel combination. We then provide a didactic example for how
quantum multiprocessing is conducted with our study and show agreement between
our method and standard serial quantum computations.

## Hyperparameter tuning: gun-point data set

While the benefits of using our temporally aware QCC-net are clear in the
context of the synthetic data set discussed in the main article, realistic
data requires more careful assessment to see if (i) the allowance for time
dependent inner product spaces or (ii) the optimization over kernel parameters
and coefficients has a marked effect when it comes to performance. We
conducted such a study and the results are shown in Fig. S1.

The study is structured as follows. Using the gun-point data set described in
the main text, we performed an SVM classification on the testing time-series
instances given (1) a value of the penalty scaling hyperparameter
$\lambda\in\\{0.1,0.3,0.8\\}$, (2) a value for the batch size
$N_{batch}\in\\{4,8,16,50\\}$ and (3) a method for obtaining a combined kernel
function used in the SVM classification. Looking further into point (3), we
look at four cases: (i) The case where the variational parameters
$\bm{\theta}$ and the kernel coefficients $\bm{\eta}$ are optimized and the
inner product space is permitted to be time-dependent (red bars on Fig. S1);
the full form of the QCC net proposed in the main article. (ii) The same
optimization as (i) but time-dependence is forbidden in the inner product
space so time-dependence is only retained classically (i.e. through the time-
dependent kernel weights; grey bars on Fig. S1). (iii) Kernels are generated
with a random value of $\bm{\theta}$, $\bm{\eta}$ is optimized and time-
dependence is allowed in the inner product space (dark blue bars on Fig. S1).
(iv) Kernels are generated with a random value of $\bm{\theta}$, the kernel
weights are $\bm{\eta}$ are fixed as $1/p$ and time-dependence is allowed in
the inner product space (gold bars on Fig. S1). After each run, we record the
balanced accuracy score $A_{B}$, the $F_{1}$ scores, the ROC AUC score (RA in
shorthand) and the kernel alignments $\mathcal{A}_{train/test}$ given by
$1-1/N\sum_{i}\sum_{j}|K_{train/test,ij}-K_{truth,ij}|$ where N is the total
number of matrix elements, $K_{train/test,ij}$ are the elements of the kernel
matrix obtained using the given combination of (1), (2) and (3) on either the
training or testing data set and $K_{truth,ij}$ are the truth matrix elements
given by $y_{i}y_{j}$ for labels $y_{i},y_{j}\in\\{1,-1\\}$. Each result is
given as the average of five random initializations of the parameters
$\bm{\theta}$. Those runs where $\bm{\theta}$ is optimized are trained for 500
mini-batch iterations.

We see that of all the runs, the scenario achieving the best overall
performance in all of the considered metrics is low $\lambda$, low $N_{batch}$
and the full optimization over all parameters in a time-dependent QCC-net.
Specifically, best performance is seen at $\lambda=0.1$ and $N_{batch}=4$ thus
motivating the settings used in the main article. It should be noted that
gains are small over the other methods and for intermediate $\lambda$ and
$N_{batch}$ are often seen at all. There are also some cases where the full
time-aware QCC-net is slightly outperformed by the version with the time-
independent inner product space (grey bars on Fig. S1). It should be noted
also that despite all of these variations in $A_{B}$, $F_{1}$ and $RA$, we
only see very slight variations in $\mathcal{A}_{train}$ and
$\mathcal{A}_{test}$. Indeed, this can be rationalized by examining Fig. 6(c)
and 6(d) from the main article by examining how many time-series instances sit
so close to the decision boundary. It is quite possible that very marginal
changes in the kernel alignment bring close call time-series instances over
the decision boundary to becoming correctly classified.

Figure S1: A study on the performance of the time-resolved QCC-net for
different values of the $\lambda$ and batch size $N_{batch}$ hyperparameters.
The study also includes the effects of removing time dependence from the inner
product space and from turning on or off various parts of the optimization
procedure.

## Quantum Multiprocessing

In the past few years, companies such as Google, IBM, IonQ, Quantinuum, and
Rigetti have achieved remarkable advancements in the development of quantum
computers. A number of companies have now launched quantum computers with
around 50 qubits, with IBM leading the way with a 127-qubit quantum computer -
the largest number of qubits available at the time of writing [71]. Despite
the fact that the accuracy of quantum computers is dependent on the
implementation technique of the Noisy-Intermediate Scale Quantum (NISQ) [72]
hardware, such as superconducting or trapped ion, regardless of their
architecture, they are prone to various noises, errors, and decoherence.
Contemporary quantum devices are restricted in their ability to produce
dependable results for practical computing problems. Despite vendors releasing
more advanced quantum computers with higher Quantum Volume [73, 62], which is
a means of quantifying a quantum device’s computational power, these devices
are still a long way from achieving quantum supremacy for practical problems.
Moreover, the attainment of fault-tolerant quantum computers is currently
unfeasible and may remain so for several decades. Consequently, it is
imperative to maximize the utilization of NISQ devices to gain a quantum
advantage on such devices.

To optimize the utilization of NISQ devices, various methods have been
proposed. Variational Quantum Algorithms (VQAs) [74], including the Quantum
Approximate Optimization Algorithm (QAOA) [52] and the Variational Quantum
Eigensolver (VQE) [75], have gained significant attention as potential
approaches to achieve quantum advantage on Noisy Intermediate-Scale Quantum
(NISQ) devices. An alternative strategy involves enhancing classical-quantum
algorithms, such as Quantum Amplitude Estimation (QAE) [76, 77, 78, 79, 80]
and Grover’s search [81, 82, 83], so that they can be implemented on NISQ
devices. These methods are aimed at formulating quantum algorithms that are
compatible with NISQ devices, but face challenges with respect to long quantum
circuit depth. Although these algorithms offer better precision and fidelity
due to reduced circuit depth, the use of classical computers in the process
weakens the potential for quantum speedup.

Figure S2: Quantum circuit diagram of the QMP. This diagram shows the quantum
circuit diagram of the QMP overlapping four independent circuits.

(a) Whole measurement

(b) Partial measurement

Figure S3: The measurement histogram of the circuit of Fig. S2 with $10000$
trials (shots). Figure (a) shows the result of the whole measurement of the
circuit in Fig. S2. The result of ‘circuit 1’ is the partial measurement of
qubit $q_{0}$ and $q_{1}$. This partial measurement in Fig. (b) is extracted
by the partial measurement algorithm in Appendix A in Ref. [67] based on the
whole measurement shown in Fig. (a). Figure S4: The measurement histogram of
each circuit (2 qubits) in the circuit diagram in Fig. S2. This histogram is
the results of each circuit consisting of the QMP circuit in Fig. S2. This
histogram is presented to compare with the QMP result shown in Fig. S3 (b).

###  VI.1 Example

In this section, we explain the QMP with an example from the Qiskit simulator.
In the bit ordering convention, the rightmost and the uppermost qubits are the
least significant bit (LSB) in the bra-ket notation and the quantum circuit,
respectively.

Figure S2 shows a QMP quantum circuit diagram of four quantum circuits for the
training for the kernelized time-series classification with the gun-point
data. The QMP circuit parallelizes four independent quantum circuits. The
first quantum circuit (‘circuit 1’) is placed on the qubit $q_{0}$ and
$q_{1}$, and the second quantum circuit (‘circuit 2’) is placed on the qubit
$q_{2}$ and $q_{3}$, respectively.

From the sampling result of $10,000$ trials (shots) shown in Fig. S3 (a), we
extract the partial measurement of each circuit. We reimplement the partial
measurement code based on the code in Appendix A in Ref. [67]. The partial
measurement results of each circuit of the QMP are shown in Fig. S3 (b). The
results of ‘circuit 1’ and ‘circuit 2’ are the partial measurements of qubit
$q_{0}$, $q_{1}$ and $q_{2}$, $q_{3}$, respectively.

The Trial Reduction Factor (TRF) [31] describes the ratio of the number of
trials (shots) that executes individually (baseline) to the number of trials
(shots) of the QMP circuit. This factor represents the efficiency of the QMP
in terms of trials (shots). In this example, the TRF is $4$ because each
circuit has $10,000$ shots and the total shot number is $40,000$ when the
circuits are executed in serial.

###  VI.2 Result Fidelity

Result fidelity is a benchmark metric introduced in Ref. [84]. It quantifies
the distance between the output probabilistic distribution $P_{\text{out}}$
with the ideal $P_{\text{ideal}}$ with the normalization in terms of the
uniform distribution $P_{\text{uni}}$. The definition is

$F(P_{\text{out}},P_{\text{ideal}})=\max\\{F_{raw}(P_{\text{out}},P_{\text{ideal}}),0\\}$
(S1)

where

$F_{raw}(P_{\text{out}},P_{\text{ideal}})=\frac{f(P_{\text{out}},P_{\text{ideal}})-f(P_{\text{uni}},P_{\text{ideal}})}{1-f(P_{\text{uni}},P_{\text{ideal}})}$
(S2)

with the distance

$f(P_{1},P_{2})=\left(\sum_{j}\sqrt{P_{1}(j)P_{2}(j)}\right)^{2}.$ (S3)

Note that we have $F(P_{\text{ideal}},P_{\text{ideal}})=1$,
$F(P_{\text{uni}},P_{\text{ideal}})=0$, and
$F(P_{\text{out}},P_{\text{ideal}})\leq 1$. The result fidelity quantifies the
degree of degradation due to the noises. Hence, we quantify the difference
between the QMP output and the outputs of each circuit (2 qubit) execution
(comparison between Fig. S3 (b) and Fig. S4). We set the outputs of each
circuit execution as the ideal result distribution and compare the QMP
outputs. In this setting, we have $0.99999$, $0.99889$, $0.99986$, and
$0.99979$ result fidelity of ‘circuit 1’, ‘circuit 2’, ‘circuit 3’, and
‘circuit 4’, respectively. When we consider the finite sampling number of the
distribution, this result fidelity shows the QMP result is consistent with the
independent circuit execution result.

The gun-point data consists of 50 training and 50 test time-series instances.
The training data (matrices) for the SVM has symmetrical $50\times 50$
dimensions, and the test data has $150\times 50$ dimensions. The number of
circuits to compute for each training sample is $1275=(50\times 51)/2$, and
$7500=150\times 50$ for each test sample. By applying the QMP, the number of
device runs is significantly reduced to 37 and 215 times for each training and
test matrix. As a preliminary step, we first verified that the QMP produced
consistent outputs with the conventional method (non-QMP) by comparing both
performances on an IBM simulator and a device with 35 test circuits. Due to
the limit of the maximum number of qubits on a simulator, we used 12 circuits,
whereas all 35 circuits were run on a device for the QMP at a time. Table 1
demonstrates that QMP and non-QMP generated almost identical outcomes with
over 0.99 fidelity on both simulation (NQ_S vs Q_S) and real quantum machine
(NQ_D vs Q_D) environments. The fidelity between a simulator and a device for
the QMP (Q_S vs Q_D) also showed over 0.95, which is even higher than the non-
QMP score, 0.91 (NQ_S vs NQ_D). Since both the simulator and device comparison
between QMP and non-QMP show more than 0.99 fidelity, we conclude that the QMP
implementation does not have accuracy degradation with respect to the
canonical (non-QMP) implementation in the experiments.

Table 1:  Result fidelity [84] scores between Non-QMP (NQ), QMP (Q) on the simulator (S) and device (D).  | NQ_S | Q_S | NQ_D | Q_D
---|---|---|---|---
NQ_S | NA | 0.9997 | 0.9134 | 0.91769
Q_S | 0.9997 | NA | 0.9446 | 0.9508
NQ_D | 0.9134 | 0.9446 | NA | 0.9909
Q_D | 0.9177 | 0.9508 | 0.9909 | NA

###  VI.3 Partial Measurement Code

⬇

1import numpy as np

2

3def get_partial_measurement(count, least, n):

4 """

5 parameters:

6 - count: return of result.get_counts()

7 the total number of counts is the number shots (e.g., 8192)

8 - least: the least significant bit index of the partial measurement

9 - n : the number of bits for partial measurement

10 """

11 cnt = {}

12 for i in range(2**n):

13 key = f’{{0:0>{n}b}}’.format(i)

14 cnt[key] = 0

15

16 for m, c in count.items(): # measurement keys, counts

17 start_idx = len(m) - (least+n)

18 end_idx = len(m) - least

19

20 if key == m[start_idx:end_idx]:

21 cnt[key] += c

22 return cnt

Listing 1: A sample code for partial measurement extraction from the whole
measurement

## References

  * Fawaz  _et al._ [2019] H. I. Fawaz, G. Forestier, J. Weber, L. Idoumghar, and P.-A. Muller, Deep learning for time series classification: a review, Data Mining and Knowledge Discovery 33, 917 (2019).
  * Blázquez-García  _et al._ [2021] A. Blázquez-García, A. Conde, U. Mori, and J. A. Lozano, A review on outlier/anomaly detection in time series data, ACM Computing Surveys 54, 1 (2021).
  * Choi  _et al._ [2021] K. Choi, J. Yi, C. Park, and S. Yoon, Deep learning for anomaly detection in time-series data: Review, analysis, and guidelines, IEEE Access 9, 120043 (2021).
  * Clark  _et al._ [2020] S. Clark, R. J. Hyndman, D. Pagendam, and L. M. Ryan, Modern strategies for time series regression, International Statistical Review 88, 10.1111/insr.12432 (2020).
  * Deb  _et al._ [2017] C. Deb, F. Zhang, J. Yang, S. E. Lee, and K. W. Shah, A review on time series forecasting techniques for building energy consumption, Renewable and Sustainable Energy Reviews 74, 902 (2017).
  * Torres  _et al._ [2021] J. F. Torres, D. Hadjout, A. Sebaa, F. Martínez-Álvarez, and A. Troncoso, Deep learning for time series forecasting: A survey, Big Data 9, 3 (2021).
  * Zhang  _et al._ [2018] C. Zhang, S. R. Kuppannagari, R. Kannan, and V. K. Prasanna, Generative adversarial network for synthetic time series data generation in smart grids, in  _2018 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm)_ (IEEE, 2018).
  * Hochreiter and Schmidhuber [1997] S. Hochreiter and J. Schmidhuber, Long short-term memory, Neural Computation 9, 1735 (1997).
  * Melis  _et al._ [2019] G. Melis, T. Kočiský, and P. Blunsom, Mogrifier lstm (2019).
  * Nguyen  _et al._ [2020] T. Nguyen, R. Baraniuk, A. Bertozzi, S. Osher, and B. Wang, Momentumrnn: Integrating momentum into recurrent neural networks, in  _Advances in Neural Information Processing Systems_ , Vol. 33, edited by H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin (Curran Associates, Inc., 2020) pp. 1924–1936.
  * Vaswani  _et al._ [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin, Attention is all you need, in  _Advances in Neural Information Processing Systems_ , Vol. 30, edited by I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Curran Associates, Inc., 2017).
  * Yang  _et al._ [2021] C.-H. H. Yang, Y.-Y. Tsai, and P.-Y. Chen, Voice2series: Reprogramming acoustic models for time series classification, in  _International Conference on Machine Learning_ (PMLR, 2021) pp. 11808–11819.
  * Zerveas  _et al._ [2021] G. Zerveas, S. Jayaraman, D. Patel, A. Bhamidipaty, and C. Eickhoff, A transformer-based framework for multivariate time series representation learning, in  _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery; Data Mining_ (ACM, 2021).
  * Cortes and Vapnik [1995] C. Cortes and V. Vapnik, Support-vector networks, Machine learning 20, 273 (1995).
  * Badiane  _et al._ [2018] M. Badiane, M. O’Reilly, and P. Cunningham, Kernel methods for time series classification and regression., in  _AICS_ (2018) pp. 54–65.
  * Bailly [2018] A. Bailly,  _Time series classification algorithms with applications in remote sensing_ , Ph.D. thesis, Université Rennes 2 (2018).
  * Fábregues de los Santos [2017] L. Fábregues de los Santos,  _Forecasting financial time series using multiple Kernel Learning_ , Master’s thesis, Universitat Politècnica de Catalunya (2017).
  * Rüping [2001] S. Rüping,  _SVM kernels for time series analysis_ , Tech. Rep. (Technical report, 2001).
  * Tino [2020] P. Tino, Dynamical systems as temporal feature spaces, The Journal of Machine Learning Research 21, 1649 (2020).
  * Gönen and Alpaydın [2011a] M. Gönen and E. Alpaydın, Multiple kernel learning algorithms, The Journal of Machine Learning Research 12, 2211 (2011a).
  * Aiolli and Donini [2015] F. Aiolli and M. Donini, EasyMKL: a scalable multiple kernel learning algorithm, Neurocomputing 169, 215 (2015).
  * Ghukasyan  _et al._ [2023] A. Ghukasyan, J. S. Baker, O. Goktas, J. Carrasquilla, and S. K. Radha, Quantum-classical multiple kernel learning (2023), arXiv:2305.17707 [quant-ph] .
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Noisy intermediate-scale quantum algorithms, Rev. Mod. Phys. 94, 015004 (2022).
  * Wiebe  _et al._ [2010] N. Wiebe, D. Berry, P. Høyer, and B. C. Sanders, Higher order decompositions of ordered operator exponentials, Journal of Physics A: Mathematical and Theoretical 43, 065203 (2010).
  * Cîrstoiu  _et al._ [2020] C. Cîrstoiu, Z. Holmes, J. Iosue, L. Cincio, P. J. Coles, and A. Sornborger, Variational fast forwarding for quantum simulation beyond the coherence time, npj Quantum Information 6, 10.1038/s41534-020-00302-0 (2020).
  * Radha [2021] S. K. Radha, Quantum constraint learning for quantum approximate optimization algorithm (2021).
  * Horowitz  _et al._ [2022] H. Horowitz, P. Rao, and S. K. Radha, A quantum generative model for multi-dimensional time series using hamiltonian learning (2022).
  * Baker  _et al._ [2022] J. S. Baker, H. Horowitz, S. K. Radha, S. Fernandes, C. Jones, N. Noorani, V. Skavysh, P. Lamontangne, and B. C. Sanders, Quantum variational rewinding for time series anomaly detection (2022).
  * Gibbs  _et al._ [2022] J. Gibbs, Z. Holmes, M. C. Caro, N. Ezzell, H.-Y. Huang, L. Cincio, A. T. Sornborger, and P. J. Coles, Dynamical simulation via quantum machine learning with provable generalization (2022).
  * Caro  _et al._ [2022] M. C. Caro, H.-Y. Huang, N. Ezzell, J. Gibbs, A. T. Sornborger, L. Cincio, P. J. Coles, and Z. Holmes, Out-of-distribution generalization for learning quantum dynamics (2022).
  * Das  _et al._ [2019] P. Das, S. S. Tannu, P. J. Nair, and M. Qureshi, A case for multi-programming quantum computers, in  _Proceedings of the 52nd Annual IEEE/ACM International Symposium on Microarchitecture_ (2019) pp. 291–303.
  * Ratanamahatana and Keogh [2005] C. A. Ratanamahatana and E. Keogh, Three myths about dynamic time warping data mining, in  _Proceedings of the 2005 SIAM international conference on data mining_ (SIAM, 2005) pp. 506–510.
  * Gönen and Alpaydın [2011b] M. Gönen and E. Alpaydın, Multiple kernel learning algorithms, The Journal of Machine Learning Research 12, 2211 (2011b).
  * Schuld and Killoran [2019] M. Schuld and N. Killoran, Quantum machine learning in feature hilbert spaces, Physical Review Letters 122, 10.1103/physrevlett.122.040504 (2019).
  * Havlíček  _et al._ [2019] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Supervised learning with quantum-enhanced feature spaces, Nature 567, 209 (2019).
  * Stone [1932] M. H. Stone, On one-parameter unitary groups in hilbert space, The Annals of Mathematics 33, 643 (1932).
  * Bjorken and Drell [1965] J. D. Bjorken and S. D. Drell,  _Relativistic quantum fields_ (McGraw-Hill, 1965).
  * Sakoe and Chiba [1978] H. Sakoe and S. Chiba, Dynamic programming algorithm optimization for spoken word recognition, IEEE Transactions on Acoustics, Speech, and Signal Processing 26, 43 (1978).
  * Aiolli  _et al._ [2008] F. Aiolli, G. D. S. Martino, and A. Sperduti, A kernel method for the optimization of the margin distribution, in  _Artificial Neural Networks - ICANN 2008_ (Springer Berlin Heidelberg, 2008) pp. 305–314.
  * O’Donoghue  _et al._ [2016a] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd, Conic optimization via operator splitting and homogeneous self-dual embedding, Journal of Optimization Theory and Applications 169, 1042 (2016a).
  * Bergholm  _et al._ [2018] V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, S. Ahmed, V. Ajith, M. S. Alam, G. Alonso-Linaje, B. AkashNarayanan, A. Asadi, J. M. Arrazola, U. Azad, S. Banning, C. Blank, T. R. Bromley, B. A. Cordier, J. Ceroni, A. Delgado, O. Di Matteo, A. Dusko, T. Garg, D. Guala, A. Hayes, R. Hill, A. Ijaz, T. Isacsson, D. Ittah, S. Jahangiri, P. Jain, E. Jiang, A. Khandelwal, K. Kottmann, R. A. Lang, C. Lee, T. Loke, A. Lowe, K. McKiernan, J. J. Meyer, J. A. Montañez-Barrera, R. Moyard, Z. Niu, L. J. O’Riordan, S. Oud, A. Panigrahi, C.-Y. Park, D. Polatajko, N. Quesada, C. Roberts, N. Sá, I. Schoch, B. Shi, S. Shu, S. Sim, A. Singh, I. Strandberg, J. Soni, A. Száva, S. Thabet, R. A. Vargas-Hernández, T. Vincent, N. Vitucci, M. Weber, D. Wierichs, R. Wiersema, M. Willmann, V. Wong, S. Zhang, and N. Killoran, Pennylane: Automatic differentiation of hybrid quantum-classical computations (2018).
  * Mitarai  _et al._ [2018] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Quantum circuit learning, Physical Review A 98, 10.1103/physreva.98.032309 (2018).
  * Schuld  _et al._ [2019] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, and N. Killoran, Evaluating analytic gradients on quantum hardware, Physical Review A 99, 10.1103/physreva.99.032331 (2019).
  * Agrawal  _et al._ [2019a] A. Agrawal, B. Amos, S. Barratt, S. Boyd, S. Diamond, and J. Z. Kolter, Differentiable convex optimization layers, in  _Advances in Neural Information Processing Systems_ , Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).
  * Agrawal  _et al._ [2019b] A. Agrawal, S. Barratt, S. Boyd, E. Busseti, and W. M. Moursi, Differentiating through a cone program (2019b).
  * Bottou  _et al._ [2018] L. Bottou, F. E. Curtis, and J. Nocedal, Optimization methods for large-scale machine learning, SIAM Review 60, 223 (2018).
  * Keskar  _et al._ [2017] N. S. Keskar, D. Mudigere, J. Nocedal, M. Smelyanskiy, and P. T. P. Tang, On large-batch training for deep learning: Generalization gap and sharp minima (2017), arXiv:1609.04836 [cs.LG] .
  * Kübler  _et al._ [2021] J. M. Kübler, S. Buchholz, and B. Schölkopf, The inductive bias of quantum kernels (2021), arXiv:2106.03747 [quant-ph] .
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 10.1038/s41467-018-07090-4 (2018).
  * Pesah  _et al._ [2021] A. Pesah, M. Cerezo, S. Wang, T. Volkoff, A. T. Sornborger, and P. J. Coles, Absence of barren plateaus in quantum convolutional neural networks, Phys. Rev. X 11, 041011 (2021).
  * Lloyd  _et al._ [2020] S. Lloyd, M. Schuld, A. Ijaz, J. Izaac, and N. Killoran, Quantum embeddings for machine learning (2020), arXiv:2001.03622 [quant-ph] .
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm, arXiv preprint arXiv:1411.4028 (2014).
  * Hadfield  _et al._ [2019] S. Hadfield, Z. Wang, B. O’Gorman, E. G. Rieffel, D. Venturelli, and R. Biswas, From the quantum approximate optimization algorithm to a quantum alternating operator ansatz, Algorithms 12 (2019).
  * Schuld  _et al._ [2020] M. Schuld, A. Bocharov, K. M. Svore, and N. Wiebe, Circuit-centric quantum classifiers, Phys. Rev. A 101, 032308 (2020).
  * Welch  _et al._ [2014] J. Welch, D. Greenbaum, S. Mostame, and A. Aspuru-Guzik, Efficient quantum circuits for diagonal unitaries without ancillas, New Journal of Physics 16, 033040 (2014).
  * Kingma and Ba [2017] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2017), arXiv:1412.6980 [cs.LG] .
  * O’Donoghue  _et al._ [2016b] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd, Conic optimization via operator splitting and homogeneous self-dual embedding, Journal of Optimization Theory and Applications 169, 1042 (2016b).
  * O’Donoghue [2021] B. O’Donoghue, Operator splitting for a homogeneous embedding of the linear complementarity problem, SIAM Journal on Optimization 31, 1999 (2021).
  * [59] See Supplemental Material [publisher url] for prelimary studies on the impact of different components of the QCC-net as tested on the gun-point data set referenced in the main article and for further explanation of how QMP was implemented.
  * Bowles  _et al._ [2023] J. Bowles, V. J. Wright, M. Farkas, N. Killoran, and M. Schuld, Contextuality and inductive bias in quantum machine learning (2023), arXiv:2302.01365 [quant-ph] .
  * Hastie  _et al._ [2009] T. Hastie, R. Tibshirani, J. H. Friedman, and J. H. Friedman,  _The elements of statistical learning: data mining, inference, and prediction_ , Vol. 2 (Springer, 2009) p. 134.
  * Cross  _et al._ [2019] A. W. Cross, L. S. Bishop, S. Sheldon, P. D. Nation, and J. M. Gambetta, Validating quantum computers using randomized model circuits, Physical Review A 100, 032328 (2019).
  * Ohkura  _et al._ [2022] Y. Ohkura, T. Satoh, and R. Van Meter, Simultaneous execution of quantum circuits on current and near-future NISQ systems, IEEE Transactions on Quantum Engineering (2022).
  * Liu and Dou [2021] L. Liu and X. Dou, QuCloud: A new qubit mapping mechanism for multi-programming quantum computing in cloud environment, in  _2021 IEEE International Symposium on High-Performance Computer Architecture (HPCA)_ (IEEE, 2021) pp. 167–178.
  * Niu and Todri-Sanial [2021] S. Niu and A. Todri-Sanial, Enabling multi-programming mechanism for quantum computing in the NISQ era, arXiv preprint arXiv:2102.05321 (2021).
  * Niu and Todri-Sanial [2022] S. Niu and A. Todri-Sanial, How parallel circuit execution can be useful for NISQ computing?, in  _2022 Design, Automation & Test in Europe Conference & Exhibition (DATE)_ (IEEE, 2022) pp. 1065–1070.
  * Park  _et al._ [2023] G. Park, K. Zhang, K. Yu, and V. Korepin, Quantum multi-programming for Grover’s search, Quantum Information Processing 22, 54 (2023).
  * Murali  _et al._ [2020] P. Murali, D. C. McKay, M. Martonosi, and A. Javadi-Abhari, Software mitigation of crosstalk on noisy intermediate-scale quantum computers, in  _Proceedings of the Twenty-Fifth International Conference on Architectural Support for Programming Languages and Operating Systems_ (2020) pp. 1001–1016.
  * Ohkura [2021] Y. Ohkura, Crosstalk-aware NISQ Multi-programming, Faculty Policy Manage., Keio Univ., Tokyo, Japan (2021).
  * Hubregtsen  _et al._ [2022] T. Hubregtsen, D. Wierichs, E. Gil-Fuster, P.-J. H. S. Derks, P. K. Faehrmann, and J. J. Meyer, Training quantum embedding kernels on near-term quantum computers, Phys. Rev. A 106, 042431 (2022).
  * Yirka [2021] B. Yirka, IBM announces development of 127-qubit quantum processor (2021).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Moll  _et al._ [2018] N. Moll, P. Barkoutsos, L. S. Bishop, J. M. Chow, A. Cross, D. J. Egger, S. Filipp, A. Fuhrer, J. M. Gambetta, M. Ganzhorn, _et al._ , Quantum optimization using variational algorithms on near-term quantum devices, Quantum Science and Technology 3, 030503 (2018).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Variational quantum algorithms, Nature Reviews Physics 3, 625 (2021).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, A variational eigenvalue solver on a photonic quantum processor, Nature communications 5, 1 (2014).
  * Aaronson and Rall [2020] S. Aaronson and P. Rall, Quantum approximate counting, simplified, Symposium on Simplicity in Algorithms , 24–32 (2020).
  * Suzuki  _et al._ [2020] Y. Suzuki, S. Uno, R. Raymond, T. Tanaka, T. Onodera, and N. Yamamoto, Amplitude estimation without phase estimation, Quantum Information Processing 19, 75 (2020).
  * Yu  _et al._ [2020] K. Yu, H. Lim, and P. Rao, Practical numerical integration on NISQ devices, in  _Quantum Information Science, Sensing, and Computation XII_ , Vol. 11391 (SPIE, 2020) p. 1139106.
  * Grinko  _et al._ [2021] D. Grinko, J. Gacon, C. Zoufal, and S. Woerner, Iterative quantum amplitude estimation, npj Quantum Information 7, 1 (2021).
  * Rao  _et al._ [2020] P. Rao, K. Yu, H. Lim, D. Jin, and D. Choi, Quantum amplitude estimation algorithms on IBM quantum devices, in  _Quantum Communications and Quantum Imaging XVIII_ , Vol. 11507 (SPIE, 2020) pp. 49–60.
  * Zhang and Korepin [2020] K. Zhang and V. E. Korepin, Depth optimization of quantum search algorithms beyond Grover’s algorithm, Physical Review A 101, 032346 (2020).
  * Zhang  _et al._ [2021] K. Zhang, P. Rao, K. Yu, H. Lim, and V. Korepin, Implementation of efficient quantum search algorithms on NISQ computers, Quantum Information Processing 20, 1 (2021).
  * Zhang  _et al._ [2022] K. Zhang, K. Yu, and V. Korepin, Quantum search on noisy intermediate-scale quantum devices, Europhysics Letters 140, 18002 (2022).
  * Lubinski  _et al._ [2021] T. Lubinski, S. Johri, P. Varosy, J. Coleman, L. Zhao, J. Necaise, C. H. Baldwin, K. Mayer, and T. Proctor, Application-oriented performance benchmarks for quantum computing, arXiv preprint arXiv:2110.03137 (2021).