# Quark: A Gradient-Free Quantum Learning
Framework for Classification Tasks

Zhihao Zhang∗
Carnegie Mellon University
<EMAIL_ADDRESS>&Zhuoming Chen∗
Tsinghua University
<EMAIL_ADDRESS>&Heyang Huang
Microsoft
<EMAIL_ADDRESS>&Zhihao Jia
Carnegie Mellon University
<EMAIL_ADDRESS>
###### Abstract

As more practical and scalable quantum computers emerge, much attention has
been focused on realizing quantum supremacy in machine learning. Existing
quantum ML methods either (1) embed a classical model into a target
Hamiltonian to enable quantum optimization or (2) represent a quantum model
using variational quantum circuits and apply classical gradient-based
optimization. The former method leverages the power of quantum optimization
but only supports simple ML models, while the latter provides flexibility in
model design but relies on gradient calculation, resulting in barren plateau
(i.e., gradient vanishing) and frequent classical-quantum interactions. To
address the limitations of existing quantum ML methods, we introduce Quark, a
gradient-free quantum learning framework that optimizes quantum ML models
using quantum optimization. Quark does not rely on gradient computation and
therefore avoids barren plateau and frequent classical-quantum interactions.
In addition, Quark can support more general ML models than prior quantum ML
methods and achieves a dataset-size-independent optimization complexity.
Theoretically, we prove that Quark can outperform classical gradient-based
methods by reducing model query complexity for highly non-convex problems;
empirically, evaluations on the Edge Detection and Tiny-MNIST tasks show that
Quark can support complex ML models and significantly reduce the number of
measurements needed for discovering near-optimal weights for these tasks.

**footnotetext: These authors contributed equally to this work

##  1 Introduction

Quantum computing provides a new computational paradigm to achieve exponential
speedups over classical counterparts for various tasks, such as cryptography
(Shor, 1994), scientific simulation (Tazhigulov et al., 2022), and data
analytics (Arute et al., 2019). A key advantage of quantum computing is its
ability to entangle multiple quantum bits, called qubits, allowing $n$ qubits
to encode a $2^{n}$-dimensional vector, while encoding this vector in
classical computing requires $2^{n}$ bits.

Inspired by this potential, recent work (Jaderberg et al., 2022; Macaluso et
al., 2020b; Torta et al., 2021; Kapoor et al., 2016; Bauer et al., 2020; Farhi
& Neven, 2018a; Schuld et al., 2014; Cong et al., 2019b) has focused on
realizing quantum speedups over classical algorithms in the field of
supervised learning. Existing quantum ML work can be divided into two
categories: classical model with quantum optimization (CMQO) and quantum model
with classical optimization (QMCO). First, CMQO methods embed a classical ML
model jointly with the optimization problem into a target Hamiltonian and
optimize the model using quantum adiabatic evolution (QAE) (Finnila et al.,
1994) or quantum approximate optimization algorithm (QAOA) (Farhi et al.,
2014; Torta et al., 2021). As the transition between a classical model and the
target Hamiltonian only applies to low-order polynomial activations (see
Figure 2), CMQO methods do not support ML models with non-linear activations
that cannot be represented in low-order polynomial (e.g., ReLU). Second, QMCO
methods optimize variational quantum models 111They are also known as
variational quantum circuits (VQC)-based models in the quantum literature. by
iteratively performing gradient descent using classical optimizers. QMCO
methods are fundamentally limited by barren plateau (i.e., gradient vanishing
(McClean et al., 2018)) and the high cost of frequent quantum-classical
interactions.

Figure 1: An overview of the Quark optimization framework. Each horizontal
line indicates a qubit, and each box on these lines represents one or multiple
quantum gates applied on these qubits. Quark’s optimization pipeline includes
three stages: (1) $\Psi_{0}$ preparation, which initializes
$R_{W},R_{D},R_{O}$ and performs model’s forward processing $U_{M}$, (2)
amplitude amplification using a Grover-based algorithm, and (3) weights
measurement. Quark uses $K$-parallel datasets (KPD) to maximize the
probability of observing highly accurate weights in each measurement. Darker
bars in the probability plots denote weights with higher accuracies.

To address the limitations of existing quantum ML methods, we introduce Quark,
a gradient-free quantum learning framework for classification tasks that
optimizes quantum models with quantum optimization (QMQO). Figure 1 shows an
overview of Quark. A key idea behind Quark is entangling the weight
222Throughout the paper, we use the term weight to refer to the set of all
trainable parameters of a model. of an ML model (i.e., the $R_{W}$ register in
Figure 1) and the encoded dataset (i.e., the $R_{D}$ register in Figure 1) in
a quantum state, where model weights that achieve optimal classification
accuracy on the training dataset can be observed with the highest
probabilities in a measurement. Therefore, users can obtain highly accurate
model weights by directly measuring the updated $R_{W}$ weight register. To
maximize the probability of observing optimal weights, we introduce two key
techniques.

Amplitude amplification. Quark uses a Grover-based mechanism to iteratively
update the probability distribution of weights based on their training
accuracies. As a result, the probability of observing weights with higher
accuracy increases after each Grover iteration, as shown in Figure 1.

$K$-parallel datasets (KPD). Applying amplitude amplification on one dataset
results in a linear amplification scenario where the measuring probability of
each weight is proportional to its training accuracy $J(w_{i})$. We further
introduce $K$-parallel datasets, a technique to enable exponential
amplification. Specifically, by entangling $k$ identical training datasets
with model weights in parallel (using $k\times R_{D}$, as shown in Figure 1),
the probability of observing weight $w_{i}$ in a measurement is proportional
to $J(w_{i})^{k}$. Therefore, as $k$ increases, the optimized probability
distribution of weights gradually converges to the optimal weights.

Compared with CMQO methods, Quark provides more flexibility in model design by
composing models directly on quantum circuits and therefore supports a broader
range of ML models. Compared with QMCO methods, Quark does not require
gradient calculation and therefore does not suffer from barren plateau. Quark
avoids frequent classical-quantum interactions by realizing both model design
and optimization fully on quantum. Besides, by using basis encoding for the
training dataset, Quark supports non-linear operations (e.g., ReLU) in its
model architecture, and the optimization complexity is independent of the
training dataset size.

Theoretically, we compare model query complexity 333Number of model forward
being called between Quark and gradient-based methods on a balanced $C$-way
classification task, and prove that Quark can outperform gradient-based
methods by reducing model query complexity for highly non-convex problems. In
addition, we prove that using $K$-parallel datasets can further reduce model
query complexity under certain circumstances.

Simulations on two tasks (i.e., Edge Detection and Tiny-MNIST) show that Quark
supports complex ML models, which can include quantum convolution, pooling,
fully connected, and ReLU layers. In addition, Quark can significantly reduce
the number of measurements needed for discovering a near-optimal weight by
applying amplitude amplification and KPD.

##### Contributions.

This paper makes the following contributions:

  * •

We propose Quark, a gradient-free quantum learning framework for
classification tasks that optimizes quantum ML models with quantum
optimization. Quark avoids barren plateau and frequent classical-quantum
interactions, supports more general models than prior quantum ML frameworks,
and achieves a dataset-size-independent optimization complexity.

  * •

Theoretically, we prove that Quark can outperform gradient-based methods by
reducing model query complexity for highly non-convex problems and that using
KPD can further reduce model query complexity.

  * •

Empirically, we show that Quark can support complex ML models and
significantly reduce the number of measurements needed for discovering a near-
optimal weight for the Edge Detection and Tiny-MNIST tasks.

##  2 Related Work

Figure 2: Comparison between CMQO, QMCO, and Quark (QMQO).

Figure 2 compares Quark with existing quantum ML approaches.

###  2.1 Classical model with quantum optimization

Existing CMQO methods aim at solving classical ML problems with quantum
optimization techniques by leveraging the advantage of quantum parallelism
(Nielsen & Chuang, 2002). Based on the well-established algorithmic foundation
in quantum annealing (Finnila et al., 1994; Kadowaki & Nishimori, 1998; Brooke
et al., 1999; Santoro et al., 2002; Santoro & Tosatti, 2006) and adiabatic
quantum computing (Farhi et al., 2001; Albash & Lidar, 2018), prior work
(Denil & De Freitas, 2011; Dumoulin et al., 2014; Adachi & Henderson, 2015)
attempts for the quantum restricted Boltzmann machine (RBM) by formulating RBM
as an Ising model (Cipra, 1987). Inspired by the quantum approximate
optimization algorithm (QAOA) (Farhi et al., 2014), Torta et al. (2021) embeds
a single binary perceptron layer into a target Hamiltonian to search for
optimal weights. However, CMQO methods are limited by the locality restriction
of the target Hamiltonian and can only embed models with low-order polynomial
activations (e.g., square). This limitation prevents CMQO methods from
supporting practical deep learning architectures, which generally contain non-
polynomial activations such as ReLU and sigmoid.

Similar to Quark, Kapoor et al. (2016) also uses Grover’s algorithm to find a
hyperplane that can perfectly separate the training dataset. However, this
method only applies to an idealistic setup where a hyperplane with perfect
classification exists in its search space. Besides, the method cannot adapt to
generic model architectures other than single-layer perceptrons.

###  2.2 Quantum model with classical optimization

Motivated by the recent advances in variational quantum algorithms (VQAs)
(Cerezo et al., 2021), QMCO methods use variational quantum circuits (VQC)
(Benedetti et al., 2019) to represent the trainable parameters of an ML model.
Havlíček et al. (2019); Schuld & Killoran (2019) use VQC as a variational
feature map to reproduce linear support vector machines (SVM) and kernel
methods on quantum circuits, which can outperform classical counterparts under
certain circumstances (Liu et al., 2021).

Besides conventional ML methods, recent work has also explored the feasibility
of classical neural networks on quantum circuits (Massoli et al., 2022). Farhi
& Neven (2018b); Macaluso et al. (2020a); Killoran et al. (2019) use VQC as
building blocks for their quantum perceptron models with a classical gradient-
based optimizer. Quantum dissipative neural network (Beer et al., 2020) (QDNN)
and quantum convolutional neural network (Cong et al., 2019a) (QCNN), on the
other hand, move a step forward towards more complicated neural architectures.
QDNN enlarges its model space by applying unitary operators on both the input
and output qubits, while QCNN uses a measurement-controlled operation to
enable non-linear operations. However, McClean et al. (2018) shows that the
barren plateau phenomenon commonly exists in VQC-based methods, where
gradients vanish exponentially with the model size. Though Beer et al. (2020)
claims to design a VQC model immune to barren plateau, Sharma et al. (2022)
contradicts such claim with analytical proof. Though Du et al. (2021) uses a
Grover algorithm as part of their method, they still require VQC as their
model building blocks that require gradients update.

Besides, due to amplitude-based data encoding, VQC-based methods in general
suffer from a linear dependency with respect to dataset size in terms of model
query complexity during training. Another drawback for amplitude encoding is
that due to the unitary constraint of quantum transformations, non-linear
operations are hard to implement for VQC-based methods. In contrast, our
method uses basis encoding that concerns only qubits state transformation
rather than amplitude transformation, which enables more efficient model query
complexity and more general non-linear transformations.

##  3 preliminaries

###  3.1 notations

Let $\mathcal{D}=\\{(x_{i},y_{i})\\}_{i\in N}$ denote a training dataset,
where $x_{i}\in\\{0,1\\}^{d_{x}}$ is the binarized feature vector associated
with the $i$-th sample, and $y_{i}\in\\{0,1\\}^{d_{y}}$ is its label. Let
$\hat{y}=f(w,\hat{x}):\\{0,1\\}^{d_{w}}\times\\{0,1\\}^{d_{x}}\rightarrow\\{0,1\\}^{d_{y}}$
denote our model parameterized by $w\in\\{0,1\\}^{d_{w}}$. Given an objective
$l(\hat{y},y)$, our goal is to find a near-optimal $w^{*}$ that
minimizes/maximizes the overall objective
$J(w)=\frac{1}{N}\sum_{(x_{i},y_{i})\in\mathcal{D}}l(f(w,x_{i}),y_{i})$. We
focus on classification tasks and use the objective
$l(\hat{y},y)=\mathbbm{1}(\hat{y}=y)$. We use $|\cdot|$ to denote the
cardinality of a set and absolute value of a scalar, and use $\|\cdot\|_{2}$
and $\|\cdot\|_{\infty}$ to denote the $L_{2}$-norm and infinity norm of a
vector. Finally, we use $\otimes$ to denote tensor product, and use $\neg$,
$\oplus$, $\wedge$, and $\vee$ to denote NEGATE, XOR, AND, and OR in logical
expressions, respectively.

###  3.2 quantum basics

A bit in the quantum regime, called a qubit, is represented by a super-
position of $|0\rangle$ and $|1\rangle$, which is formally defined as
$|z\rangle=\alpha|0\rangle+e^{i\phi}\beta|1\rangle$, where $\alpha$ and
$\beta$ are the amplitudes, and $e^{i\phi}$ is the relative phase.
Furthermore, the rule also enforces $\langle z|z\rangle=1$ where $\langle z|$
is the conjugate transpose of $|z\rangle$. In an n-qubit system, a quantum
state is represented as a superposition of $2^{n}$ basis states. For
computational simplicity, we will be using the basis states that are spanned
by $\\{|0\rangle,|1\rangle\\}^{n}$ and denote the superposition of a n-qubits
state as $|\Psi\rangle=\sum_{i=0}^{2^{n}-1}\frac{1}{\sqrt{2^{n}}}|i\rangle$
where $|i\rangle$ is the corresponding computational basis. We thus use
$|w_{i}\rangle\in\\{|0\rangle,|1\rangle\\}^{d_{w}}$ to denote weight basis
state, $|x_{j},y_{j}\rangle$ for the entangled data basis state where
$|x_{j}\rangle\in\\{|0\rangle,|1\rangle\\}^{d_{x}}$ is the feature and
$|y_{j}\rangle\in\\{|0\rangle,|1\rangle\\}^{d_{y}}$ the label.

##  4 Methodology

Figure 3: Toy illustration of Quark’s insight, where the dataset
$\\{(x_{0}=0,y_{0}=0),(x_{1}=1,y_{1}=1)\\}$ is constructed by oracle function
$g(x)=0\oplus x$. We define our learner model as $o=f(x,w)=w\oplus x$. The
colored qubit strings stands for the states activated by model forward (cyan
for $|w\rangle=0$ and yellow for $|w\rangle=1$). violet bars stand for
solution states measuring probability, gray stands for non-solution states
measuring probability.

To realize quantum supremacy, we circumvent the overhead induced by frequent
classical-quantum interactions and the gradient calculation step in the VQC-
based methods, which may result in barren plateau as the circuit depth
increases. In addition, as VQC-based methods use the amplitude encoding scheme
for data encoding, non-linear operations are hard to implement, while their
training time can linearly depend on the training dataset size.

To this end, we designed Quark in a QMQO fashion, which (1) achieves
gradients-free optimization over the entire weight distribution, (2) uses
basis encoding to achieve potential speedup as sample size scales up, and (3)
enables a flexible quantum model design that can easily incorporate non-
linearities. Figure 1 shows an overview of the quantum circuit design in
Quark.

As we are using basis encoding over the entire data distribution, we
initialize the data register $R_{D}$ as
$|D\rangle=\sum_{(x_{i},y_{i})\in\mathcal{D}}\frac{1}{\sqrt{|\mathcal{D}|}}|x_{i},y_{i}\rangle$.
To initialize a model’s weights register $R_{W}$, we construct a uniform state
by applying the Hadamard gate on each weight qubit:
$|W\rangle=\sum_{i=0}^{2^{d_{w}}-1}\frac{1}{\sqrt{2^{d_{w}}}}|w_{i}\rangle$.
In addition, Quark includes an auxiliary register $R_{O}$ for storing
intermediate results and model’s output, which is initialized as
$|O\rangle=|0\rangle^{d_{o}}$. Therefore, the initial state is represented as
$|W\rangle\otimes|D\rangle\otimes|O\rangle$. By encoding the quantum model
architecture as a unitary matrix $U_{M}$, the model’s forward processing is
defined as:

$\displaystyle\displaystyle
U_{M}\left(|W\rangle\otimes|D\rangle\otimes|O\rangle\right)=\sum_{i}\frac{1}{2^{\frac{d_{w}}{2}}}|w_{i}\rangle\sum_{j}\frac{1}{\sqrt{|\mathcal{D}|}}|x_{j},y_{j}\rangle|f(w_{i},x_{j})\rangle$
(1)

where $f(w_{i},x_{j})$ is the model’s output for weight $w_{i}$ and sample
$x_{j}$ 444We omit the model’s intermediate results $|O\rangle$ for
simplicity.. Given the above expression for the model forward processing, the
key insight for Quark is to update the probability of observing weight $w$ in
a measurement based on its training objective $J(w)$. This is achieved through
Grover’s algorithm by defining the solution state space
$\mathcal{S}=\\{|w_{i}\rangle|x_{j},y_{j}\rangle|f(w_{i},x_{j})\rangle\mid
y_{j}=f(w_{i},x_{j})\\}$. Section A.1 includes an introduction to the original
Grover’s algorithm. After amplitude amplification, we observe a model weight
by measuring the weight register $R_{W}$.

To further illustrate our insight, Figure 3 shows a toy example, where the
underline oracle function that generates the training data is $g(x)=0\oplus
x$. Using one qubit for $x$, the training dataset is constructed as
$\\{(x_{0}=0,y_{0}=0),(x_{1}=1,y_{1}=1)\\}$. By defining the learner model as
$o=f(x,w)=w\oplus x$, where $w$ is the trainable parameter, Quark includes
four qubits (i.e., $|w\rangle$, $|x\rangle$, $|y\rangle$, and $|o\rangle$).
For initial state preparation, we simply go through the model forward circuit
with a uniform weight initialization that result in an uniform amplitude
state, as shown on the left of Figure 3. As no optimization happens at this
stage, the probabilities for observing $|w\rangle=|0\rangle$ and
$|w\rangle=|1\rangle$ in a measurement are equivalent. After amplitude
amplification, Quark reaches an optimized amplitude state (shown on the right
of Figure 3), where the probability of observing the optimal weight
$|w\rangle=|0\rangle$ in a measurement is much higher.

For the rest of this section, Section 4.1 introduces Quark’s on-circuit model
design philosophy, Section 4.2 describes the Grover-based method for amplitude
amplification, and Section 4.3 introduce K-Parallel Dataset (KPD), a technique
that enables exponential amplification to further improve the optimization
algorithm.

###  4.1 Model Design

(a) Convolution

(b) MaxPool

Figure 4: Logical gates illustration of quantum module design. (a)
Convolution: for a $1\times 3$ Conv kernel with weights
$W=[W_{1},W_{2},W_{3}]$ and input $X_{i}=[X_{i1},X_{i2},X_{i3}]$, the output
is given by $O_{i}=(\neg(W_{1}\oplus X_{i1}))\wedge(\neg(W_{2}\oplus
X_{i2}))\wedge(\neg(W_{3}\oplus X_{i3}))$. (b) MaxPooling: for a input
$X_{i}=[X_{i1},X_{i2},X_{i3}]$ the max pooling output is given by
$O_{i}=X_{i1}\vee X_{i2}\vee X_{i3}$.

The model circuit $U_{M}$ can be designed through arbitrary combinations of
base gates (CNOT, Rotational gates, Hadamard gate, etc.), which can then
entangle $|W\rangle,|D\rangle$ to get a trainable model.

Leveraging the flexibility of basis encoding, Quark is capable of realizing
modules used in classical deep models. Figure 4 illustrates two Quark modules
(Conv, MaxPool) used in our experiments. Besides Convolution and MaxPooling,
Quark can support more diverse operations given enough qubits. We include the
demonstration of the Fully Connective and ReLU modules in Section A.4.

During training, we encode the entire data distribution through basis encoding
so we can apply model forward simultaneously for all data samples. During
inference, we encode each sample as a single state $|\hat{x}\rangle$.
Prediction can be obtained by measurements over the output register after the
model forward step
$U_{M}|w^{*}\rangle|\hat{x},0\rangle|0\rangle=|w^{*}\rangle|\hat{x},0\rangle|f(w^{*},\hat{x})\rangle$,
where $w^{*}$ is the measured optimal weight.

###  4.2 Amplitude Amplification

To increase the probability of measuring the weights with the highest
classification accuracy, we use Grover’s algorithm to amplify the amplitude
for any state
$|w_{i}\rangle|x_{j},y_{j}\rangle|f(w_{i},x_{j})\rangle\in\mathcal{S}$, thus
the measuring probability of optimal weights would be amplified the most. In
this section, we will formally define the unitary operators we need for
Grover’s update, followed by the complexity analysis of the optimization
algorithm.

For conventional Grover’s algorithm, we need two reflection unitaries, namely
$U_{s^{\top}},U_{\Psi_{0}}$, where $U_{s^{\top}}$ is the reflection against
the non-solution sub-space and $U_{\Psi_{0}}$ is the reflection against the
initial state. In our case, by defining solution state

$|s\rangle=\sum_{|w_{i}\rangle|x_{j},y_{j}\rangle|f(w_{i},x_{j})\rangle\in\mathcal{S}}\frac{1}{\sqrt{|\mathcal{S}|}}|w_{i}\rangle|x_{j},y_{j}\rangle|f(w_{i},x_{j})\rangle$

and initial state

$\Psi_{0}=\sum_{i,j}\frac{1}{\sqrt{2^{d_{w}}|\mathcal{D}|}}|w_{i}\rangle|x_{j},y_{j}\rangle|f(w_{i},x_{j})\rangle$

, the Grover operator can be easily built from
$U_{\text{comb}}=U_{\Psi_{0}}U_{s^{\top}}$ where
$U_{s^{\top}}=\mathbb{I}-2|s\rangle\langle
s|,U_{\Psi_{0}}=2|\Psi_{0}\rangle\langle\Psi_{0}|-\mathbb{I}$.

Now that we have the Grover operator well-defined, we can formally do analysis
on the algorithm in terms of model query complexity.

###### Theorem 1

By defining
$\alpha=\frac{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})d\delta}{\int_{w_{i}\in\mathcal{W}}J(w_{i})d\delta}$
as the probability to measure an $\epsilon$-optimal solution after the
Grover’s update, the model query complexity for getting an $\epsilon$-optimal
solution $w_{i}\in\mathcal{W_{\epsilon}}$ on a balanced $C$-way classification
task is $O(\frac{\sqrt{C}}{\alpha})$.

$J(w_{i})$ is the objective value (accuracy) for weights $w_{i}$,
$\mathcal{W}$ is the complete weight space, and $\mathcal{W}_{\epsilon}$ is
the $\epsilon$-optimal weight subspace. The proof of Theorem 1 is included in
Section A.3.1, it follows from the analysis that we need $O(\sqrt{C})$ Grover
iterations per measurement and
$O(\frac{\int_{w_{i}\in\mathcal{W}}J(w_{i})d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})d\delta})$
measurements in expectation for sampling an $\epsilon$-optimal weight.

Notice that the complexity of our method does not depend on sample size $N$.
As we are sampling solutions from an optimized distribution for a general non-
convex problem, our method does depend on an
$O(\frac{\int_{w_{i}\in\mathcal{W}}J(w_{i})d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})d\delta})$
term. However, to achieve an $\epsilon$-optimal solution on a general non-
convex problem, the worst case scenario for VQC-based methods with gradient-
based optimizers needs
$O(\frac{\int_{w_{i}\in\mathcal{W}}d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}d\delta})$
iterations to sample initial points lie within convex regions that contain
$\epsilon$-optimal solutions. With per iteration gradient evaluation cost that
is in the order of $O(N)$, this gives an overall complexity of
$O(\frac{\int_{w_{i}\in\mathcal{W}}d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}d\delta}N)$
for VQC-based methods. Since
$O(\frac{\int_{w_{i}\in\mathcal{W}}J(w_{i})d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})d\delta})<O(\frac{\int_{w_{i}\in\mathcal{W}}d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}d\delta})$
and $\sqrt{C}\ll N$ in practice. Our method provides a speed-up for finding an
$\epsilon$-optimal solution in the general non-convex setup.

In addition, as Quark does not require gradients calculation, it will not
suffer from the notorious problem of barren plateau. For more idealistic case
of convex problems, Garg et al. (2020) has given a proof of no quantum speed-
ups can be obtained in this case with no exception of our method.

###  4.3 K-Parallel Datasets (KPD)

A naive implementation can only achieve $p(w_{i})\propto J(w_{i})$, which
requires more measurements for sampling $\epsilon$-optimal solution. To this
end, we further extend our method to achieve $p(w_{i})\propto J(w_{i})^{k}$
through $K$-Parallel Dataset (KPD), as shown in Figure 1 when $k>1$. The
intuition is similar to simulated annealing, by updating weights distribution
proportional to $J(w_{i})^{k}$, the probability mass would gradually converge
to the global optimal solutions.

We do so by concatenating the same dataset $K$ times for increasing solution
states in the order of $K$ as:

$U_{M}\left(|W\rangle\bigotimes_{k=0}^{K-1}\left(|D_{k}\rangle\otimes|O_{k}\rangle\right)\right)=\sum_{i}\frac{1}{2^{\frac{d_{w}}{2}}}|w_{i}\rangle\prod_{k=0}^{K-1}\left(\sum_{j}\frac{1}{\sqrt{|\mathcal{D}_{k}|}}|x_{j},y_{j}\rangle|f(w_{i},x_{j})\rangle\right)$

the solution set is now defined as:

$\mathcal{S}_{K}:\\{|w_{i}\rangle\bigotimes_{k=0}^{K-1}|x_{j_{k}},y_{j_{k}}\rangle|f(w_{i},x_{j_{k}})\rangle\mid\bigwedge_{k=0}^{K-1}\left(y_{j_{k}}=f(w_{i},x_{j_{k}})\right)\\}$

We then update $|s\rangle,|\Psi_{0}\rangle,U_{comb}$ accordingly as in Section
A.2. Now the ratio of solution states between different weights can grow
exponentially in terms of $K$.

###### Theorem 2

By defining
$\beta=\frac{\int_{w_{i}\in\mathcal{W_{\epsilon}}}d\delta}{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}d\delta}$
as the volume ratio between $\epsilon$-optimal weight subspace
$|\mathcal{W}_{\epsilon}|$ and non-$\epsilon$-optimal weight subspace
$|\mathcal{W}-\mathcal{W}_{\epsilon}|$, the model query complexity for getting
an $\epsilon$-optimal solution $w_{i}\in\mathcal{W_{\epsilon}}$ on a balanced
$C$-way classification task with k-parallel dataset is
$O((1+\beta^{k-1}(\frac{1}{\alpha}-1)^{k})kC^{\frac{k}{2}})$.

The proof of Theorem 2 is included in Section A.3.2. Here, the trade-off is
that the number of Grover iterations needed per measurement grows to
$O\left(C^{\frac{k}{2}}\right)$. As we also need to apply the model forward
for each dataset, the KPD Grover complexity in terms of model query complexity
is now $O(kC^{\frac{k}{2}})$. On the other hand, as we are converging to the
global solution, the number of measurements needed can be reduced from
$O(\frac{1}{\alpha})$ to $O(1+\beta^{k-1}(\frac{1}{\alpha}-1)^{k})$ as
$\beta<\alpha$. Thus the optimal value for $k$ depends on the specification of
$\alpha,\beta,C$.

###### Theorem 3

Given $\alpha,\beta,C$, the optimal value of $k$ is
$k^{*}=\lfloor\log_{\frac{\alpha}{\beta}}\frac{1}{\alpha}\rfloor+1$ if
$\frac{\alpha}{\beta}\geq{m}^{{\frac{1}{m-1}}}\sqrt{C}$ where
$m=\lfloor\log_{\frac{\alpha}{\beta}}\frac{1}{\alpha}\rfloor+1$ else
$k^{*}=1$.

Thus for cases where
${m}^{{\frac{1}{m-1}}}\sqrt{C}\leq\frac{\alpha}{\beta}\ll\frac{1}{\alpha}$, we
have our optimal $k>1$ which proves the effectiveness of KPD under certain
circumstances. We include the proof of Theorem 3 in Section A.3.3.

##  5 experiments

In this section, we empirically verify the effectiveness of Quark on two
tasks, namely Edge Detection and Tiny-MNIST. Algorithm 1 formally states the
pipeline of Quark 555Preprocessing() and Evaluate() are included in Section
A.5.1. Instead of leveraging the approximated iteration number from
theoretical derivations for Grover’s update, we use a more precise but
efficient estimation in practice to achieve a more accurate amplitude
amplification effect. We include the detail of Grover iteration number
estimation in Section A.5.1. Notice that most of the simulation results are
obtained numerically since we have no access to large scale quantum devices
that fit our setup. However, we do use Qiskit Aer (Anis & et al., 2021) to
verify the reproducibility of our numerical simulation results on models that
are applicable, results are included in Section A.6.

Input: Data Oracle: $U_{D}$; Model Oracle: $U_{M}$; Objective Oracle: $U_{L}$;
Number of parallel dataset: $k$; Measurement budget: $m$; Weights buffer:
$\mathcal{B}$

Output: Optimized Model Weight: $w^{*}$

$g,U_{\text{comb}},|\Psi_{0}\rangle=$ Preprocessing$(U_{D},U_{M},U_{L},k)$ ;

// get Grover iteration $g$, Grover operator $U_{\text{comb}}$ and initial
state $|\Psi_{0}\rangle$

for  _i = 0; $i<m$; ++i_ do

for  _j = 0; $j<g$; ++j_ do

$|\Psi_{j+1}\rangle=U_{\text{comb}}|\Psi_{j}\rangle$ ;

// Grover update

$w_{i}=\textbf{Measure}(|\Psi_{g}\rangle)$ ;

// Weights measurement on $R_{W}$

$\mathcal{B}.\text{add}(w_{i})$;

return $\arg\max_{w\in\mathcal{B}}(\textbf{Evaluate}(w))$

Algorithm 1 Quark’s optimization pipeline.

###  5.1 Edge Detection

For Edge Detection, our goal is to identify if a $3\times 3$ binary matrix has
1) both vertical and horizontal lines 2) vertical lines only 3) horizontal
lines only 4) no lines, where a line is defined by three consecutive ’1’s in a
row or column. Thus the task is a 4-way classification task with 512
instances. We split the 512 samples into a training set with 400 randomly
selected instances and a test set with the rest. We use a Quantum
Convolutional Model that consists of several $1\times 3$ convolution kernels
with Maxpooling modules as described in Figure 4 for this task. Results are
demonstrated in Figure 5 (normalized accuracy:
$\hat{J}(w_{i})=\frac{J(w_{i})}{\sum_{i}J(w_{i})}$), from which we can clearly
see a linear relationship between model accuracy and optimized weights
distribution using 1-PD. For 4-PD, we can observe that weights distribution
further concentrates on the global optimal solutions.

(a) Initial distribution

(b) Amplitude ampl. w/ 1-PD

(c) Amplitude ampl. w/ 4-PD

Figure 5: Amplitude ampl. + KPD on Edge Detection. normalized accuracy:
$\hat{J}(w_{i})=\frac{J(w_{i})}{\sum_{i}J(w_{i})}$.

Indeed, as we concatenate more training datasets, the probability mass will
gradually converge to the optimal solutions, thus reducing the number of
measurements we need significantly. We include the results for 2-PD and 3-PD
in Section A.5.3.

###  5.2 Tiny-MNIST

For Tiny-MNIST, we down-sample original $28\times 28$ images from MNIST to
$3\times 3$ and remove duplicated samples. Due to device limitation, we only
consider classes $1,2,7$, which makes this task a 3-way classification task.
We use Weighted Mask modules (details are included in Section A.4) with
MaxPooling modules to compose our model, which can achieve $86.08\%$ training
accuracy and $82.61\%$ testing accuracy by the best model in the search space.
Whereas a well-optimized classical single layer perceptron model can achieve
$85\%$ in training accuracy and $80\%$ test accuracy.

Similarly, as shown in Figure 6 (normalized accuracy:
$\hat{J}(w_{i})=\frac{J(w_{i})}{\sum_{i}J(w_{i})}$), a linear relationship
between model accuracy and optimized weights distribution can be observed
using 1-PD, while 4-PD can further increase the probability for sampling
optimal solutions. We also include the results of 2-PD and 3-PD on Tiny-MNIST
in Section A.5.4.

Besides distribution evolution, we also statistically demonstrate the
relationship between measured top-1 model’s test accuracy and measurements
budget using uniform random sampling, 1-PD, and 4-PD respectively in Figure 7.
In order to achieve equally well-performed model, 4-PD only requires $\sim 30$
shots on Edge Detection task for a mean test accuracy $>98\%$ while uniform
random sampling needs $\sim 900$ ($\sim 30\times$ more) shots with much higher
variance. Similar trend can be observed on Tiny-MNIST where 4-PD only requires
$\sim 1000$ shots for a mean test accuracy $>76.5\%$ while uniform random
sampling needs $\sim 20000$ ($\sim 20\times$ more) shots with higher variance.
We include the training ones in Section A.5.5.

(a) Initial Distribution

(b) Amplitude ampl. w/ 1-PD

(c) Amplitude ampl. w/ 4-PD

Figure 6: Amplitude ampl. + KPD on Tiny-MNIST. normalized accuracy:
$\hat{J}(w_{i})=\frac{J(w_{i})}{\sum_{i}J(w_{i})}$.

(a) Edge Detection

(b) Tiny-MNIST

Figure 7: The mean $\pm$ std of the test accuracy of the best discovered model
with different measurement budget (in shots). URS shows the result where the
weights are uniformly random sampled.

##  6 conclusion

In this paper, we propose a new quantum learning framework Quark that does not
involve gradients calculation and operates in a fully-quantum fashion.
Acknowledging the notorious problem of barren plateaus from VQC based methods,
Quark shed some lights on circumventing this phenomenon through a gradient-
free optimization pipeline. Quark also enables a more general set of module
design due to basis encoding, so that non-linear operations can be easily
implemented. Theoretically, we present some evidences in terms of model query
complexity for Quark to demonstrate trade-offs between VQC based methods and
Quark. Empirically, we have verified the effectiveness of Quark through
numerical simulations on Edge Detection and Tiny-MNIST.

## References

  * Adachi & Henderson (2015) Steven H Adachi and Maxwell P Henderson.  Application of quantum annealing to training of deep neural networks.  _arXiv preprint arXiv:1510.06356_ , 2015.
  * Albash & Lidar (2018) Tameem Albash and Daniel A Lidar.  Adiabatic quantum computation.  _Reviews of Modern Physics_ , 90(1):015002, 2018\.
  * Anis & et al. (2021) Md Sajid Anis and et al.  Qiskit: An open-source framework for quantum computing, 2021.
  * Arute et al. (2019) Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C Bardin, Rami Barends, Rupak Biswas, Sergio Boixo, Fernando GSL Brandao, David A Buell, et al.  Quantum supremacy using a programmable superconducting processor.  _Nature_ , 574(7779):505–510, 2019.
  * Bauer et al. (2020) Bela Bauer, Sergey Bravyi, Mario Motta, and Garnet Kin-Lic Chan.  Quantum algorithms for quantum chemistry and quantum materials science.  _Chemical Reviews_ , 120(22):12685–12717, 2020\.
  * Beer et al. (2020) Kerstin Beer, Dmytro Bondarenko, Terry Farrelly, Tobias J. Osborne, Robert Salzmann, Daniel Scheiermann, and Ramona Wolf.  Training deep quantum neural networks.  _Nature Communications_ , 11, 2020.
  * Benedetti et al. (2019) Marcello Benedetti, Erika Lloyd, Stefan Sack, and Mattia Fiorentini.  Parameterized quantum circuits as machine learning models.  _Quantum Science and Technology_ , 4(4):043001, 2019.
  * Brooke et al. (1999) J Brooke, David Bitko, Rosenbaum, and Gabriel Aeppli.  Quantum annealing of a disordered magnet.  _Science_ , 284(5415):779–781, 1999.
  * Cerezo et al. (2021) Marco Cerezo, Andrew Arrasmith, Ryan Babbush, Simon Benjamin, Suguro Endo, Keisuke Fujii, Jarrod Ryan McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick Coles.  Variational quantum algorithms.  _Nature Reviews Physics_ , 2021.  URL https://www.nature.com/articles/s42254-021-00348-9.
  * Cipra (1987) Barry A Cipra.  An introduction to the ising model.  _The American Mathematical Monthly_ , 94(10):937–959, 1987.
  * Cong et al. (2019a) Iris Cong, Soonwon Choi, and Mikhail D. Lukin.  Quantum convolutional neural networks.  _Nature Physics_ , 15(12):1273–1278, aug 2019a.  doi: 10.1038/s41567-019-0648-8.  URL https://doi.org/10.1038%2Fs41567-019-0648-8.
  * Cong et al. (2019b) Iris Cong, Soonwon Choi, and Mikhail D Lukin.  Quantum convolutional neural networks.  _Nature Physics_ , 15(12):1273–1278, 2019b.
  * Denil & De Freitas (2011) Misha Denil and Nando De Freitas.  Toward the implementation of a quantum rbm.  2011\.
  * Du et al. (2021) Yuxuan Du, Min-Hsiu Hsieh, Tongliang Liu, and Dacheng Tao.  A grover-search based quantum learning scheme for classification.  _New Journal of Physics_ , 23(2):023020, feb 2021\.  doi: 10.1088/1367-2630/abdefa.  URL https://doi.org/10.1088/1367-2630/abdefa.
  * Dumoulin et al. (2014) Vincent Dumoulin, Ian Goodfellow, Aaron Courville, and Yoshua Bengio.  On the challenges of physical implementations of rbms.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 28, 2014.
  * Farhi & Neven (2018a) Edward Farhi and Hartmut Neven.  Classification with quantum neural networks on near term processors.  _arXiv preprint arXiv:1802.06002_ , 2018a.
  * Farhi & Neven (2018b) Edward Farhi and Hartmut Neven.  Classification with quantum neural networks on near term processors, 2018b.  URL https://arxiv.org/abs/1802.06002.
  * Farhi et al. (2001) Edward Farhi, Jeffrey Goldstone, Sam Gutmann, Joshua Lapan, Andrew Lundgren, and Daniel Preda.  A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem.  _Science_ , 292(5516):472–475, 2001.
  * Farhi et al. (2014) Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  _arXiv preprint arXiv:1411.4028_ , 2014.
  * Finnila et al. (1994) Aleta Berk Finnila, MA Gomez, C Sebenik, Catherine Stenson, and Jimmie D Doll.  Quantum annealing: A new method for minimizing multidimensional functions.  _Chemical physics letters_ , 219(5-6):343–348, 1994.
  * Garg et al. (2020) Ankit Garg, Robin Kothari, Praneeth Netrapalli, and Suhail Sherif.  No quantum speedup over gradient descent for non-smooth convex optimization.  _arXiv preprint arXiv:2010.01801_ , 2020.
  * Havlíček et al. (2019) Vojtěch Havlíček, Antonio D. Córcoles, Kristan Temme, Aram W. Harrow, Abhinav Kandala, Jerry M. Chow, and Jay M. Gambetta.  Supervised learning with quantum-enhanced feature spaces.  _Nature_ , 567(7747):209–212, Mar 2019.  ISSN 1476-4687.  doi: 10.1038/s41586-019-0980-2.  URL https://doi.org/10.1038/s41586-019-0980-2.
  * Jaderberg et al. (2022) Ben Jaderberg, Lewis W Anderson, Weidi Xie, Samuel Albanie, Martin Kiffner, and Dieter Jaksch.  Quantum self-supervised learning.  _Quantum Science and Technology_ , 7(3):035005, 2022.
  * Kadowaki & Nishimori (1998) Tadashi Kadowaki and Hidetoshi Nishimori.  Quantum annealing in the transverse ising model.  _Physical Review E_ , 58(5):5355, 1998.
  * Kapoor et al. (2016) Ashish Kapoor, Nathan Wiebe, and Krysta Svore.  Quantum perceptron models.  _Advances in neural information processing systems_ , 29, 2016.
  * Killoran et al. (2019) Nathan Killoran, Thomas R Bromley, Juan Miguel Arrazola, Maria Schuld, Nicolás Quesada, and Seth Lloyd.  Continuous-variable quantum neural networks.  _Physical Review Research_ , 1(3):033063, 2019\.
  * Liu et al. (2021) Yunchao Liu, Srinivasan Arunachalam, and Kristan Temme.  A rigorous and robust quantum speed-up in supervised machine learning.  _Nature Physics_ , 17(9):1013–1017, Sep 2021\.  ISSN 1745-2481.  doi: 10.1038/s41567-021-01287-z.  URL https://doi.org/10.1038/s41567-021-01287-z.
  * Macaluso et al. (2020a) Antonio Macaluso, Luca Clissa, Stefano Lodi, and Claudio Sartori.  A variational algorithm for quantum neural networks.  _Computational Science – ICCS 2020_ , 12142:591 – 604, 2020a.
  * Macaluso et al. (2020b) Antonio Macaluso, Luca Clissa, Stefano Lodi, and Claudio Sartori.  A variational algorithm for quantum neural networks.  In _International Conference on Computational Science_ , pp. 591–604. Springer, 2020b.
  * Massoli et al. (2022) Fabio Valerio Massoli, Lucia Vadicamo, Giuseppe Amato, and Fabrizio Falchi.  A leap among quantum computing and quantum neural networks: A survey.  _ACM Comput. Surv._ , mar 2022.  ISSN 0360-0300.  doi: 10.1145/3529756.  URL https://doi.org/10.1145/3529756.  Just Accepted.
  * McClean et al. (2018) Jarrod R McClean, Sergio Boixo, Vadim N Smelyanskiy, Ryan Babbush, and Hartmut Neven.  Barren plateaus in quantum neural network training landscapes.  _Nature communications_ , 9(1):1–6, 2018.
  * Nielsen & Chuang (2002) Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information, 2002.
  * Santoro & Tosatti (2006) Giuseppe E Santoro and Erio Tosatti.  Optimization using quantum mechanics: quantum annealing through adiabatic evolution.  _Journal of Physics A: Mathematical and General_ , 39(36):R393, 2006.
  * Santoro et al. (2002) Giuseppe E Santoro, Roman Martonák, Erio Tosatti, and Roberto Car.  Theory of quantum annealing of an ising spin glass.  _Science_ , 295(5564):2427–2430, 2002.
  * Schuld & Killoran (2019) Maria Schuld and Nathan Killoran.  Quantum machine learning in feature hilbert spaces.  _Phys. Rev. Lett._ , 122:040504, Feb 2019.  doi: 10.1103/PhysRevLett.122.040504.  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.040504.
  * Schuld et al. (2014) Maria Schuld, Ilya Sinayskiy, and Francesco Petruccione.  The quest for a quantum neural network.  _Quantum Information Processing_ , 13(11):2567–2586, 2014.
  * Sharma et al. (2022) Kunal Sharma, M. Cerezo, Lukasz Cincio, and Patrick J. Coles.  Trainability of dissipative perceptron-based quantum neural networks.  _Physical Review Letters_ , 128(18), may 2022.  doi: 10.1103/physrevlett.128.180505.  URL https://doi.org/10.1103%2Fphysrevlett.128.180505.
  * Shor (1994) Peter W Shor.  Algorithms for quantum computation: discrete logarithms and factoring.  In _Proceedings 35th annual symposium on foundations of computer science_ , pp. 124–134. Ieee, 1994.
  * Tazhigulov et al. (2022) Ruslan Tazhigulov, Shi-Ning Sun, Reza Haghshenas, Huanchen Zhai, Adrian Tan, Nicholas Rubin, Ryan Babbush, Austin Minnich, and Garnet Kin-Lic Chan.  Simulating challenging correlated molecules and materials on the sycamore quantum processor.  _arXiv:2203.15291_ , 2022.  URL https://arxiv.org/abs/2203.15291.
  * Torta et al. (2021) Pietro Torta, Glen B Mbeng, Carlo Baldassi, Riccardo Zecchina, and Giuseppe E Santoro.  Quantum approximate optimization algorithm applied to the binary perceptron.  _arXiv preprint arXiv:2112.10219_ , 2021.

##  Appendix A Appendix

###  A.1 Grover’s algorithm

A well-known algorithm for amplitude amplification is Grover’s algorithm. The
original Grover’s algorithm is trying to search specific states that satisfy
some properties which are called solution states. Instead of enumerating over
all possible states to find the solution states that lie in the solution set
$\mathcal{S}$, Grover’s algorithm tries to amplify the solution states’
amplitudes using two reflection unitary matrices $U_{s^{\top}},U_{\Psi_{0}}$.
Let $|\Psi_{0}\rangle$ denote the initial state of all qubits and
$|s\rangle=\sum_{s_{i}\in\mathcal{S}}\frac{1}{\sqrt{|\mathcal{S}|}}|s_{i}\rangle$
represent the basis state spanned by all solution states. $U_{s^{\top}}$ and
$U_{\Psi_{0}}$ are constructed as follows:

$\displaystyle U_{s^{\top}}$ $\displaystyle=\mathbb{I}-2|s\rangle\langle s|$
(2) $\displaystyle U_{\Psi_{0}}$
$\displaystyle=2|\Psi_{0}\rangle\langle\Psi_{0}|-\mathbb{I}$ (3)

Geometrically, $U_{s^{\top}}$ is the reflection operator over
$|s^{\top}\rangle=\sum_{s_{i}\notin\mathcal{S}}\frac{1}{\sqrt{|\\{|0\rangle,|1\rangle\\}^{n}\backslash\mathcal{S}|}}|s_{i}\rangle$,
which is an state orthogonal to the solution space. Similarly, $U_{\Psi_{0}}$
is the reflection operator over $|\Psi_{0}\rangle$. Given
$\theta=\arcsin{(\langle\Psi_{0}|s^{\top}\rangle)}$, $|\Psi_{0}\rangle$ can be
expressed as:

$|\Psi_{0}\rangle=\cos{\theta}|s^{\top}\rangle+\sin{\theta}|s\rangle$

The combination of the two reflection unitary matrices
$U_{\text{comb}}=U_{\Psi_{0}}U_{s^{\top}}$ is equivalent to a rotation of
$2\theta$ on the plane spanned by $|s^{\top}\rangle$ and $|\Psi_{0}\rangle$.
Therefore, applying the combination of the two reflection matrices $k$ times
gives:

$\displaystyle|\Psi_{k}\rangle$
$\displaystyle=\left(\prod_{i=1}^{k}U_{\text{comb}}\right)|\Psi_{0}\rangle$
(4)
$\displaystyle=\cos{((2k+1)\theta)}|s^{\top}\rangle+\sin{((2k+1)\theta)}|s\rangle$
(5)

As for most practical setups the solutions are always sparse and existed, we
have $0<\theta\ll\frac{\pi}{3}$. To maximize the amplitude for $|s\rangle$,
$k$ should be in the order of $O(\frac{1}{\theta})$.

###  A.2 Definition

The $|s\rangle,|\Psi_{0}\rangle$ in KPD are updated as:

$|s\rangle=\sum_{|w_{i}\rangle\bigotimes_{k=0}^{K-1}|x_{j_{k}},y_{j_{k}}\rangle|f(w_{i},x_{j_{k}})\rangle\in\mathcal{S}_{K}}\frac{1}{\sqrt{|\mathcal{S}_{K}|}}|w_{i}\rangle\bigotimes_{k=0}^{K-1}|x_{j_{k}},y_{j_{k}}\rangle|f(w_{i},x_{j_{k}})\rangle$

and

$\Psi_{0}=\sum_{i,j_{0},\cdots,j_{K-1}}\frac{1}{\sqrt{2^{d_{w}}|\mathcal{D}|^{K}}}|w_{i}\rangle\bigotimes_{k=0}^{K-1}|x_{j_{k}},y_{j_{k}}\rangle|f(w_{i},x_{j_{k}})\rangle$

,again the Grover operator can be easily built from
$U_{\text{comb}}=U_{s^{\top}}U_{\Psi_{0}}$ where
$U_{s^{\top}}=\mathbb{I}-2|s\rangle\langle
s|,U_{\Psi_{0}}=2|\Psi_{0}\rangle\langle\Psi_{0}|-\mathbb{I}$.

###  A.3 Proof

####  A.3.1 Theorem 1

Given the probability for sampling a solution state is $p$, then the Grover
iterations to achieve the maximum amplitude amplification effect is
$\sqrt{\frac{1}{p}}$. As in our case, the probability for sampling a solution
state is given by
$\frac{\int_{w_{i}\in\mathcal{W}}J(w_{i})d\delta}{\int_{w_{i}\in\mathcal{W}}d\delta}$,
then suppose we are doing a balanced $C$-way classification task the Grover
iterations we need is:

$\displaystyle\sqrt{\frac{\int_{w_{i}\in\mathcal{W}}d\delta}{\int_{w_{i}\in\mathcal{W}}J(w_{i})d\delta}}$
$\displaystyle=$
$\displaystyle\sqrt{\frac{1}{\mathbb{E}_{w_{i}\sim\mathcal{W}}\left[J(w_{i})\right]}}$
(6) $\displaystyle=$ $\displaystyle\sqrt{\frac{1}{\frac{1}{C}}}$ (7)
$\displaystyle=$ $\displaystyle\sqrt{C}$ (8)

$(9)\rightarrow(10)$ is due to for a balanced $C$-way classification, we
should expect the accuracy for a random model to be the same as a random guess
$\frac{1}{C}$

As we have defined
$\alpha=\frac{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})d\delta}{\int_{w_{i}\in\mathcal{W}}J(w_{i})d\delta}$
as the probability to measure a $\epsilon$-optimal solution, thus in order to
sample a $\epsilon$-optimal solution by measurements, it takes
$O(\frac{1}{\alpha})$ trials. Which gives us an overall complexity of
$O(\frac{\sqrt{C}}{\alpha})$

####  A.3.2 Theorem 2

Given we are doing a balanced $C$-way classification task for $k$ parallel
dataset, the probability for sampling a solution state is now
$\frac{\int_{w_{i}\in\mathcal{W}}J(w_{i})^{k}d\delta}{\int_{w_{i}\in\mathcal{W}}1^{k}d\delta}=\mathbb{E}_{w_{i}\sim\mathcal{W}}[J(w_{i})^{k}]$.
Since the objective function we defined is non-negative $J(w_{i})\geq 0$, thus
we have $J(w_{i})^{k}$ to be a convex function in $J(w_{i})$. As expectation
operator preserve convexity we have
$\mathbb{E}[J(w_{i})^{k}]\geq\mathbb{E}[J(w_{i})]^{k}$, we have:

$\displaystyle\sqrt{\frac{\int_{w_{i}\in\mathcal{W}}1^{k}d\delta}{\int_{w_{i}\in\mathcal{W}}J(w_{i})^{k}d\delta}}$
$\displaystyle=$
$\displaystyle\sqrt{\frac{1}{\mathbb{E}_{w_{i}\sim\mathcal{W}}[J(w_{i})^{k}]}}$
(9) $\displaystyle\leq$
$\displaystyle\sqrt{\frac{1}{\mathbb{E}_{w_{i}\sim\mathcal{W}}[J(w_{i})]^{k}}}$
(10) $\displaystyle=$ $\displaystyle\sqrt{C^{k}}$ (11)

Thus the Grover iteration needed is upper bounded by $C^{\frac{k}{2}}$

For measurement, the expected iterations we need is
$\frac{\int_{w_{i}\in\mathcal{W}}J(w_{i})^{k}d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})^{k}d\delta}$
which can be approximated as:

$\displaystyle\frac{\int_{w_{i}\in\mathcal{W}}J(w_{i})^{k}d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})^{k}d\delta}$
$\displaystyle=$
$\displaystyle\frac{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})^{k}d\delta+\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}J(w_{i})^{k}d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})^{k}d\delta}$
(12) $\displaystyle=$ $\displaystyle
1+\frac{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}J(w_{i})^{k}d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})^{k}d\delta}$
(13) $\displaystyle=$ $\displaystyle
1+\frac{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}J(w_{i})^{k}d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})^{k}d\delta}\frac{\int_{w_{i}\in\mathcal{W_{\epsilon}}}1^{k}d\delta}{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}1^{k}d\delta}\frac{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}1^{k}d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}1^{k}d\delta}$
(14) $\displaystyle=$ $\displaystyle
1+\frac{\frac{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}J(w_{i})^{k}d\delta}{{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}1^{k}d\delta}}}{\frac{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})^{k}d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}1^{k}d\delta}}\frac{1}{\beta}$
(15) $\displaystyle=$ $\displaystyle
1+\frac{\mathbb{E}_{w_{i}\in\mathcal{W/W_{\epsilon}}}[J(w_{i})^{k}]}{\mathbb{E}_{w_{i}\in\mathcal{W_{\epsilon}}}[J(w_{i})^{k}]}\frac{1}{\beta}$
(16) $\displaystyle\approx$ $\displaystyle
1+\frac{\mathbb{E}_{w_{i}\in\mathcal{W/W_{\epsilon}}}[J(w_{i})]^{k}}{\mathbb{E}_{w_{i}\in\mathcal{W_{\epsilon}}}[J(w_{i})]^{k}}\frac{1}{\beta}$
(17) $\displaystyle=$ $\displaystyle
1+(\frac{\mathbb{E}_{w_{i}\in\mathcal{W/W_{\epsilon}}}[J(w_{i})]}{\mathbb{E}_{w_{i}\in\mathcal{W_{\epsilon}}}[J(w_{i})]})^{k}\frac{1}{\beta}$
(18) $\displaystyle=$ $\displaystyle
1+(\frac{\frac{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}J(w_{i})d\delta}{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}1d\delta}}{\frac{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}1d\delta}})^{k}\frac{1}{\beta}$
(19) $\displaystyle=$ $\displaystyle
1+(\frac{\int_{w_{i}\in\mathcal{W_{\epsilon}}}1d\delta}{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}1d\delta}\frac{\int_{w_{i}\in\mathcal{W/W_{\epsilon}}}J(w_{i})d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})d\delta})^{k}\frac{1}{\beta}$
(20) $\displaystyle=$ $\displaystyle
1+(\beta\frac{\int_{w_{i}\in\mathcal{W}}J(w_{i})d\delta-\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})d\delta}{\int_{w_{i}\in\mathcal{W_{\epsilon}}}J(w_{i})d\delta})^{k}\frac{1}{\beta}$
(21) $\displaystyle=$ $\displaystyle
1+(\beta(\frac{1}{\alpha}-1))^{k}\frac{1}{\beta}$ (22) $\displaystyle=$
$\displaystyle 1+\beta^{k-1}(\frac{1}{\alpha}-1)^{k}$ (23)

$(18)\rightarrow(19)$ is assuming
$\mathbb{E}_{w_{i}\in\mathcal{W}/\mathcal{W}_{\epsilon}}[J(w_{i})]>\text{Var}_{w_{i}\in\mathcal{W}/\mathcal{W}_{\epsilon}}[J(w_{i})]^{\frac{1}{2}}$
which is usually the case in practice since $J(w_{i})$ itself is upper bounded
by $[0,1-O(\epsilon)]$ within the non-$\epsilon$ optimal region and $J(w_{i})$
can be assumed to be near-uniformly distributed within
$\mathcal{W}/\mathcal{W}_{\epsilon}$ across range $[0,1-O(\epsilon)]$, which
is a general assumption. Thus
$\frac{\mathbb{E}_{w_{i}\in\mathcal{W/W_{\epsilon}}}[J(w_{i})^{k}]}{\mathbb{E}_{w_{i}\in\mathcal{W_{\epsilon}}}[J(w_{i})^{k}]}$
can be dominated by
$\frac{\mathbb{E}_{w_{i}\in\mathcal{W/W_{\epsilon}}}[J(w_{i})]^{k}}{\mathbb{E}_{w_{i}\in\mathcal{W_{\epsilon}}}[J(w_{i})^{k}]}\leq\frac{\mathbb{E}_{w_{i}\in\mathcal{W/W_{\epsilon}}}[J(w_{i})]^{k}}{\mathbb{E}_{w_{i}\in\mathcal{W_{\epsilon}}}[J(w_{i})]^{k}}$.

####  A.3.3 Theorem 3

As the trade-off is between measurements needed versus Grover iterations per
measurement, we can list their rate of change to make the comparison. For
$k$-parallel dataset, Grover iterations is increased by
$kC^{\frac{k-1}{2}}=\frac{kC^{\frac{k}{2}}}{C^{\frac{1}{2}}}$. For
measurements needed, the exact ratio should be
$\frac{\frac{1}{\alpha}}{1+\beta^{k-1}(\frac{1}{\alpha}-1)^{k}}$, however we
can further approximate this ratio as
$\frac{\frac{1}{\alpha}}{\frac{\beta}{\alpha}^{k-1}\frac{1}{\alpha}}=(\frac{\alpha}{\beta})^{k-1}$
given $\beta^{k-1}(\frac{1}{\alpha}-1)^{k}>1$ which is equivalent to
$k\leq\lfloor\log_{\frac{\alpha}{\beta}}\frac{1}{\alpha}\rfloor+1$. Thus the
optimal $k$ should satisfy $\frac{\alpha}{\beta}\geq
k^{\frac{1}{k-1}}C^{\frac{1}{2}}$ as well as
$k\leq\lfloor\log_{\frac{\alpha}{\beta}}\frac{1}{\alpha}\rfloor+1$. As
$k^{\frac{1}{k-1}}C^{\frac{1}{2}},k\geq 2$ is monotonically decreasing with
respect to $k$ thus as long as $\frac{\alpha}{\beta}\geq
m^{\frac{1}{m-1}}C^{\frac{1}{2}}$ for
$m=\lfloor\log_{\frac{\alpha}{\beta}}\frac{1}{\alpha}\rfloor+1$ we could have
our optimal $k=m=\lfloor\log_{\frac{\alpha}{\beta}}\frac{1}{\alpha}\rfloor+1$,
otherwise $k=1$.

###  A.4 Module Design Details

(a) Fully connected layer

(b) ReLU

Figure 8: Representing a fully connected layer and ReLU activation as quantum
circuits in the Quark optimization framework. (a) Fully connective layer: for
a $2\times 3$ FC layer with weights $W_{ij}$ and input
$X=[X_{1},X_{2},X_{3}]$, the output is given by $O_{i}=(W_{i1}\wedge
X_{1})+(W_{i2}\wedge X_{2})+(W_{i3}\wedge X_{3})$. (b) ReLU: for an input
$X_{i}=[S,X_{1},X_{2}]$ where $S$ is the sign qubit, the ReLU output is given
by $X^{\prime}=(X\oplus X)\vee(\neg S\wedge X)$.

In addition to the convolution and maxpooling layers shown in Figure 4(a) and
Figure 4(b), Quark can also incorporate other commonly used tensor algebra
operators. Figure 8 demonstrates how to represent a fully connected and a ReLU
layer as quantum circuits in Quark.

Due to the limitation of the number of qubits, our parameterization can be
viewed as a binary model over a bounded weight space $\mathcal{W}$. In our
experiments, for Edge Detection task, the model output is formally defined as:

$\displaystyle O_{0}=(\bigvee_{i=0}^{2}(\bigwedge_{j=0}^{2}W_{0,j}\oplus
X_{i,j}))\oplus W_{0,3}$ (24) $\displaystyle
O_{1}=(\bigvee_{i=0}^{2}(\bigwedge_{j=0}^{2}W_{1,j}\oplus X_{j,i}))\oplus
W_{1,3}$ (25)

We use $|O_{0}O_{1}\rangle=|00\rangle,|01\rangle,|10\rangle,|11\rangle$ to
express the 4 different predictions respectively.
For Tiny MNIST task, the model output is formally defined as:

$\displaystyle
O_{k}=(\bigvee_{i=0}^{2}(\bigvee_{j=0}^{2}W_{k,3i+j}X_{i,j}))\oplus
W_{k,9}\quad k=0,1$ (26)

We use $|O_{0}O_{1}\rangle=|10\rangle,|11\rangle$ to express the prediction of
Number 1 and use $|O_{0}O_{1}\rangle=|01\rangle,|00\rangle$ to express the
prediction of Number 2 and Number 7 respectively.

###  A.5 Experiments

####  A.5.1 Quark Process

(a) $QCA$

(b) $QCB$

Figure 9: $QCA$ is the quantum circuit for calculating the average accuracy
for different parameters, which is used in computing the number of Grover
iterations as shown in Algorithm 2. $QCB$ shows the quantum circuit for
computing the accuracy of a specific weight, which is used in Algorithm 3.

We use $\textbf{Preprocessing}()$ to:

  * •

Calculate number of Grover iterations we need (The procedure is illustrated in
Figure 9(a) and Algorithm 2).

  * •

Construct Grover Operator $U_{comb}$.

  * •

Uniformly initialize $R_{W}$ with Hadamard gates.

  * •

Initialize $k\times R_{D}$ with $U_{D}$s to encode k identical training
dataset in parallel.

In Algorithm 2, we may find that $\theta$ is close to $\frac{\pi}{4}$, making
$G$ extremely large or non-existent. However, we can introduce extra samples
to dataset to resolve the issue. These auxiliary samples are identified by an
additional qubit to classify them as non-solution states.

In function $\textbf{Evaluate}()$, we evaluate the objective value of a
specific weight $w$. The procedure is illustrated in Figure 9(b) and Algorithm
3.

Input: Data Encoder Oracle: $U_{D}$; Model Forward Oracle: $U_{M}$; Objective
Function Oracle: $U_{L}$; number of parallel dataset: $k$; Shots to estimate
accuracy: $s$

Output: Grover Iterations Number: $g$

Initialize $\text{Obj}=0$;

$\text{QC}=\text{QCA}(U_{D},U_{M},U_{L},k)$ ;

// Construct the quantum circuit shown in Figure 9(a)

for  _i = 0; $i<s$; ++i_ do

$\text{Obj}=\text{Obj}+\text{QC.measure}(R_{O}[-1])$

$\text{Obj}=\text{Obj}/S$;

$\theta=\arcsin(\text{Obj})$;

$g=[\frac{m\pi+\frac{\pi}{2}-\theta}{2\theta}]$, $m\in\mathcal{N}$;

return $g$;

Algorithm 2 Get Grover Iterations Number

Input: Data Encoder Oracle: $U_{D}$; Model Forward Oracle: $U_{M}$; Loss
Function Oracle: $U_{L}$; number of parallel dataset: $k$; Shots to estimate
accuracy: $s$; Weight: $w$

Output: $J(w)$

Initialize $J=0$;

$\text{QC}=\text{QCB}(U_{D},U_{M},U_{L},k,w)$ ;

// Construct the quantum circuit shown in Figure 9(b)

for  _i = 0; $i<s$; ++i_ do

$J=J+\text{QC.measure}(R_{O}[-1])$

$J=J/s$;

return $J$;

Algorithm 3 Evaluate

####  A.5.2 Tiny-MNIST Experimental Setup

To construct Tiny-MNIST dataset, we select images with label 1, 2, 7 form the
original MNIST training set and downsample them to 3x3 images with
binarization applied to form $D_{1}$. As samples with same representations in
$D_{1}$ cannot share different labels, we use majority voting to decide labels
for duplicate samples within $D_{1}$ to form our final dataset. We apply same
procedure for both training and testing datasets. This gives us a Tiny-MNIST
dataset with 79 instances for training and 46 instances for testing. We use
Tiny-MNIST for evaluating both classical methods and Quark. The settings of
classical methods are shown in Table 1.

Architecture |  $9\times 3$ FC + Softmax
---|---
Loss Function | CrossEntropyLoss
Optimizer | Adadelta
learning rate | 1.0
batch size | 8
epochs | 100
seed | 0,1,2,3,4,5,6,7,8,9
Table 1:  Basic Settings of Classical Method

####  A.5.3 Edge Detection Results

(a) 2-PD

(b) 3-PD

Figure 10: The effect of amplitude amplification for Edge Detection.

####  A.5.4 Tiny-MNIST Results

(a) 2-PD

(b) 3-PD

Figure 11: The effect of amplitude amplification for Tiny-MNIST.

####  A.5.5 Shots & Training Accuracy

(a) Edge Detection

(b) Tiny-MNIST

Figure 12: mean $\pm$ std of the best model’s training accuracy measured
within a given measurement budget.

###  A.6 QUANTUM SIMULATION

To further verify our framework, we use Qiskit Aer (Anis & et al., 2021) to
simulate the process of solving a simplified Edge Detection task with Quark.
The goal is to identify whether a 3$\times$3 binary matrix has horizontal
lines. The task is a binary classification task with 512 instances. We
randomly select 400 of them as the training set and the rest as the test set.
The result of simulation is shown in Table 2. Each data point is acquired with
20 runs. The corresponding result of numerical simulation is shown in Table 3.

shots | 1 | 2 | 4 | 8 | 16 | 32 | 64 | 128
---|---|---|---|---|---|---|---|---
train | 59.14% | 67.9% | 68.36% | 75.33% | 85.04% | 85.79% | 91.93% | 94.14%
std | 12.21% | 10.34% | 9.05% | 11.26% | 9.27% | 9.01% | 6.18% | 4.30%
test | 58.17% | 65.04% | 65.54% | 75.63% | 84.69% | 84.55% | 90.18% | 92.41%
std | 13.54% | 10.35% | 9.60% | 11.32% | 10.09% | 10.59% | 7.33% | 5.77%
Table 2: Relationships between training $\&$ test accuracy and shots. (Qiskit Simulation)  shots | 1 | 2 | 4 | 8 | 16 | 32 | 64 | 128
---|---|---|---|---|---|---|---|---
train | 60.74% | 66.20% | 71.52% | 77.49% | 83.82% | 89.59% | 93.93% | 96.88%
std | 10.72% | 10.26% | 10.96% | 11.12% | 10.09% | 7.87% | 5.53% | 3.86%
test | 60.03% | 64.99% | 70.11% | 76.18% | 82.83% | 88.90% | 93.22% | 95.96%
std | 11.17% | 11.07% | 11.97% | 12.28% | 11.23% | 8.72% | 6.11% | 4.95%
Table 3: Relationships between training $\&$ test accuracy and shots.
(Numerical Simulation)

Figure 13: The quantum circuit for simplified edge detection.


# Long-Term Autonomous Ocean Monitoring with Streaming Samples

Weizhe Chen, Lantao Liu
School of Informatics, Computing, and Engineering
Indiana University, Bloomington, IN 47408, USA.
E-mail: {chenweiz<EMAIL_ADDRESS>
###### Abstract

In the autonomous ocean monitoring task, the sampling robot moves in the
environment and accumulates data continuously. The widely adopted spatial
modeling method — standard Gaussian process (GP) regression — becomes
inadequate in processing the growing sensing data of a large size. To overcome
the computational challenge, this paper presents an environmental modeling
framework using a sparse variant of GP called streaming sparse GP (SSGP). The
SSGP is able to handle streaming data in an online and incremental manner, and
is therefore suitable for long-term autonomous environmental monitoring. The
SSGP summarizes the collected data using a small set of pseudo data points
that best represent the whole dataset, and updates the hyperparameters and
pseudo point locations in a streaming fashion, leading to high-quality
approximation of the underlying environmental model with significantly reduced
computational cost and memory demand.

##  I Introduction and Related Work

To autonomously monitor our aquatic environments such as the oceans,
intelligent robotic platforms such as the autonomous underwater vehicles
(AUVs) and unmanned surface vessels (USVs) have been increasingly utilized in
scientific information gathering missions due to the attractive mobility,
flexibility, and adaptivity of these platforms [12]. Fig. 1 illustrates the
interaction loop between a sampling robot and the environment. Typically, the
robot needs to first compute and plan a sampling path following which
environmental samples can be collected. The path is computed based on some
learned environmental model, which is usually a probabilistic model of
targeted environmental attributes such as the non-uniform salinity of the
ocean. Then, the robot samples the environment and uses the newly obtained
samples to update its estimated model, which in turn influences the sampling
path computation in the next round.

Figure 1:  A typical loop between the sampling robot and the environment. The
example here is to estimate the bathymetry of ocean floor.

The sampling path computation for environmental sensing and monitoring is
usually formulated as the informative path planning, the objective of which is
to maximize the “informativeness” of collected samples. Such informative path
navigates the robot to collect samples from locations of greater importance
for better estimating the underlying environmental model. Oftentimes, the
informativeness of a new observation (sampling) is quantified by its entropy,
variance reduction of prediction, or mutual information between the sampled
and un-sampled spaces. Representative informative planning methods include,
e.g., recursive-greedy algorithms that exploit the submodularity property of
mutual information [15, 18, 26, 32], dynamic programming algorithms that
select a set of waypoints with maximum information [8, 21, 22, 27], sampling-
based motion planning algorithms [16, 19], Bayesian optimization [2, 20, 24],
evolutionary methods [29], temporal difference reinforcement learning [10],
and Monte Carlo tree search [1, 3, 9, 28].

Most existing informative planners have been built upon probabilistic
regression models. This also requires the environmental modeling component to
be efficient and accurate. Gaussian process (GP) regression has become the de-
facto standard for spatial modeling due to its well-calibrated uncertainty
estimate, modeling flexibility, and robustness to overfitting. Nevertheless,
there are still many challenges to be overcome. For instance, one of the big
challenges is the limited scalability to large dataset. The time complexity
and the storage complexity of a vanilla GP regression with $N$ collected data
samples are $\mathcal{O}(N^{3})$ and $\mathcal{O}(N^{2})$, respectively.
Another challenge is the lack of a principled method to deal with streaming
data. As the robot collects data sequentially, we would like to update
prediction and (hyper) parameters of our model in a real-time and incremental
fashion. However, most existing work tackled the problem by initializing the
parameters through a pilot survey or prior data, and keeping them fixed during
the sampling process [4, 13, 14, 17]. Unfortunately, fixed parameters will
undoubtedly limit the adaptivity of the sampling robot.

This paper presents a sparse variant of GP which is able to handle streaming
data in an online and incremental manner, and is therefore suitable for long-
term autonomous environmental monitoring. Specifically, the predictive
distribution, hyperparameters, and the locations of pseudo points (a small set
of points summarizing the training data) are updated in real-time, leading to
an accurate estimation of the environmental state with reasonable
computational cost and memory demand. This work relates to sparse online GP
(SOGP) [11, 21] and sparse pseudo-inputs GP (SPGP) [27, 33]. The key
difference between this work and our former SOGP-based environmental sampling
work [21, 23] is that, instead of employing crafted heuristics for optimizing
and selecting a set of representative sparse samples, we present a principled
framework for learning hyperparameters and optimizing the representative
pseudo inputs. Also, the SPGP framework treats the pseudo inputs as additional
kernel hyperparameters, which might lead to overfitting when we jointly
optimize hyperparameters and pseudo inputs [34]. Moreover, in SPGP, there is
no discrepancy measure between the approximate model and the exact model. In
contrast, the SSGP [5] models the sparse pseudo points as variational
parameters which are optimized by minimizing the Kullback-Leibler (KL)
divergence between the approximate GP and the exact posterior GP (see Section
III). Such a mechanism prevents the learning algorithm from overfitting and
leads to high-quality approximation.

##  II Background

In this section, we briefly review GP regression and a variational framework
for sparse GP regression.

###  II-A Gaussian Process Regression

Assume we have $N$ training inputs (i.e. sampling locations) and the
corresponding real-value outputs (namely, observations)
$\\{\bm{x}_{n},y_{n}\\}^{N}_{n=1}$. In the vanilla GP regression, we assume
$y_{n}=f(\bm{x}_{n})+\epsilon_{n}$, where $f$ is an unknown function and the
observation noise is drawn from a Gaussian distribution
$\epsilon\sim\mathcal{N}(0,\sigma^{2}_{y})$. For notational simplicity, $f$ is
typically assumed to be drawn from a zero-mean GP prior
$f\sim\mathcal{GP}(\bm{0},k(\cdot,\cdot|\bm{\theta}))$, where
$k(\cdot,\cdot|\bm{\theta})$ is the covariance function with a set of
hyperparameters $\bm{\theta}$. After taking the observed data into account,
our prior GP can be updated to a posterior GP which is specified by a
posterior mean function and a posterior covariance function: 111We have
aggregated the observations into a vector $\bm{y}=\\{y_{n}\\}^{N}_{n=1}$.

$\displaystyle m_{\bm{y}}(\bm{x})$
$\displaystyle=K_{\bm{x}N}(K_{NN}+\sigma^{2}_{y}I)^{-1}\bm{y}$ $\displaystyle
k_{\bm{y}}(\bm{x},\bm{x}^{\prime})$
$\displaystyle=k(\bm{x},\bm{x}^{\prime})-K_{\bm{x}N}(K_{NN}+\sigma^{2}_{y}I)^{-1}K_{N\bm{x}^{\prime}}.$

Here, $K_{\bm{x}n}$ is a $1\times n$ covariance vector between the test input
$\bm{x}$ and the $N$ training inputs, and $K_{N\bm{x}}=K_{\bm{x}N}^{\top}$.
$K_{NN}$ is an $N\times N$ covariance matrix on the training inputs, and $I$
is an $N\times N$ identity matrix. The posterior GP depends on the values of
the hyperparameters which can be optimized by maximizing the log marginal
likelihood given by

$\log
p(\bm{y})=-\frac{1}{2}\bm{y}^{\top}K^{-1}\bm{y}-\frac{1}{2}\log|K|-\frac{n}{2}\log(2\pi).$

Here we use the shorthand $K=K_{NN}+\sigma^{2}_{y}I$.

However, the standard GP suffers from poor scalability because it requires
matrix inverse operations with time complexity of $\mathcal{O}(N^{3})$ and
memory (storage space) complexity of $\mathcal{O}(N^{2})$. This computational
challenge has led to many sparse approximation paradigms [6, 30]. Next, we
will introduce a variational framework for the sparse GP.

###  II-B Variational Sparse Gaussian Process

Variational sparse GP (VSGP) approximates the intractable posterior
$p(f|\bm{y},\bm{\theta})$ through an approximate posterior $q(f)$. We measure
the discrepancy between the approximate posterior and the true posterior using
Kullback-Leibler (KL) divergence. Then we can minimize this quantity so that
the approximate posterior is as “close” as possible to the exact posterior.
Since $p(f|\bm{y},\bm{\theta})$ is unknown, computing the KL divergence
$\mathbb{KL}[q(f)||p(f|\bm{y},\bm{\theta})]$ is infeasible. Fortunately, using
Bayes’ theorem, we have the following important equation:

$\mathbb{KL}[q(f)||p(f|\bm{y},\bm{\theta})]=\log
p(\bm{y}|\bm{\theta})-\mathbb{ELBO}[q(f),\bm{\theta}],$ (1)

where $\mathbb{ELBO}[q(f),\bm{\theta}]=\int
q(f)\log\frac{p(\bm{y},f|\bm{\theta})}{q(f)}\mathop{}\\!\mathrm{d}f$ is the
Evidence Lower BOund (ELBO). Since $\log p(\bm{y}|\bm{\theta})$ does not
depend on $q(f)$, maximizing the ELBO w.r.t. $q(f)$ is equivalent to
minimizing the KL divergence, which implies that the approximate posterior
gets closer to the true posterior. Furthermore, given the fact that the KL
divergence is non-negative, ELBO lower bounds the marginal likelihood $\log
p(\bm{y}|\bm{\theta})$ so that it can be used for learning the hyperparameters
$\bm{\theta}$.

We assume that the approximate posterior has the form
$q(f)=p(f_{\neq\bm{u}}|\bm{u},\bm{\theta})q(\bm{u})$, where $\bm{u}$ \- the
pseudo points - is a small subset of $f$, $q(\bm{u})$ is a variational
distribution over $\bm{u}$, and $p(f_{\neq\bm{u}}|\bm{u},\bm{\theta})$ is the
conditional prior of $f_{\neq\bm{u}}$. This assumption induces cancellation of
the uncountbly infinite parts in the equation and provides a computationally
tractable lower bound:

$\displaystyle\mathbb{ELBO}=\int
q(f)\log\frac{p(\bm{y}|f,\bm{\theta})p(\bm{u}|\bm{\theta})\cancel{p(f_{\neq\bm{u}}|\bm{u},\bm{\theta})}}{\cancel{p(f_{\neq\bm{u}}|\bm{u},\bm{u})}q(\bm{u})}\mathop{}\\!\mathrm{d}f.$

The closed-form expression for the optimal variational distribution can be
obtained by maximizing $\mathbb{ELBO}[q(\bm{u}),\bm{\theta}]$ w.r.t.
$q(\bm{u})$. Hyperparameters can also be computed by maximizing
$\mathbb{ELBO}[q(\bm{u}),\bm{\theta}]$ w.r.t. $\bm{\theta}$.

This framework assumes that the data arrives in batches. In such streaming
setting, we need to append newly acquired data to a continuously-growing
dataset, and then re-train the model with the combined dataset. In the
following section, we shall discuss how to update the approximate posterior
and hyperparameters by integrating the information from old approximation and
the new data in a mathematically sound way.

##  III Streaming Sparse Gaussian Process

Our goal is to derive the new approximate posterior and marginal likelihood
using the old approximation and the new data. Formally, let
$q_{\text{old}}(f)$ be the approximate posterior obtained at the previous
step. According to Bayes’ rule, we have

$\displaystyle q_{\text{old}}(f)$ $\displaystyle\approx
p(f|\bm{y}_{\text{old}})=\frac{p(f|\bm{\theta}_{\text{old}})p(\bm{y}_{\text{old}}|f)}{Z_{1}(\bm{\theta}_{\text{old}})}$
(2) $\displaystyle p(f|\bm{y}_{\text{old}},\bm{y}_{\text{new}})$
$\displaystyle=\frac{p(f|\bm{\theta}_{\text{new}})p(\bm{y}_{\text{old}}|f)p(\bm{y}_{\text{new}}|f)}{Z_{2}(\bm{\theta}_{\text{new}})},$
(3)

where $Z_{1}(\bm{\theta}_{\text{old}}),Z_{2}(\bm{\theta}_{\text{new}})$ are
the normalizing constants. Since the new posterior should only rely on the new
data $\bm{y}_{\text{new}}$ and the old approximation, we rearrange Eq. (2) and
substitute $p(\bm{y}_{\text{old}}|f)$ into Eq. (3), which yields

$\hat{p}(f|\bm{y}_{\text{old}},\bm{y}_{\text{new}})=\frac{Z_{1}(\bm{\theta}_{\text{old}})}{Z_{2}(\bm{\theta}_{\text{new}})}p(f|\bm{\theta}_{\text{new}})p(\bm{y}_{\text{new}}|f)\frac{q_{\text{old}}(f)}{p(f|\bm{\theta}_{\text{old}})}.$

This new posterior fuses the old approximation $q_{\text{old}}(f)$, the new
likelihood $p(\bm{y}_{\text{new}}|f)$, and our priors. However, we cannot use
this as the new approximate posterior,
$q_{\text{new}}(f)=\hat{p}(f|\bm{y}_{\text{old}},\bm{y}_{\text{new}})$,
because this recovers exact GP regression and it is intractable [5].
Therefore, we consider approximation by minimizing the KL divergence between
$q_{\text{new}}(f)$ and $p(\bm{y}_{\text{new}}|f)$.

Let $\bm{a}=f(\bm{z}_{\text{old}})$ and $\bm{b}=f(\bm{z}_{\text{new}})$ be the
function values at the pseudo-inputs before and after seeing new data. Note
that the number of pseudo points $M_{\bm{a}}=|\bm{a}|$ and
$M_{\bm{b}}=|\bm{b}|$ are not necessarily the same, and the new pseudo inputs
might be different from the old ones. This is required when new regions of
input space are gradually explored in the environmental monitoring scenario.
The forms of the approximate posteriors are assumed to be
$q_{\text{old}}(f)=p(f_{\neq\bm{a}}|\bm{a},\bm{\theta}_{\text{old}})q_{\text{old}}(\bm{a})$
and
$q_{\text{new}}(f)=p(f_{\neq\bm{b}}|\bm{b},\bm{\theta}_{\text{new}})q_{\text{new}}(\bm{b})$
as that of VSGP.

$\displaystyle\underbrace{\text{KL}\Big{[}q_{\text{new}}(f)\Big{|}\Big{|}\hat{p}(f|\bm{y}_{\text{old}},\bm{y}_{\text{new}})\Big{]}}_{\text{non-
negative}}=\underbrace{\log\frac{Z_{2}(\bm{\theta}_{\text{new}})}{Z_{1}(\bm{\theta}_{\text{old}})}}_{\text{constant}}$
$\displaystyle-\underbrace{\int
q_{\text{new}}(f)\left[\log\frac{p(\bm{b}|\bm{\theta}_{\text{new}})q_{\text{old}}(\bm{a})p(\bm{y}_{\text{new}}|f))}{p(\bm{a}|\bm{\theta}_{\text{old}})q_{\text{new}}(\bm{b})}\right]\operatorname{d}\\!{f}}_{\mathbb{ELBO}(q(\bm{b}),\bm{\theta}_{\text{new}})}.$

Since the constant term does not depend on $q(\bm{b})$, maximizing
$\mathbb{ELBO}(q(\bm{b}),\bm{\theta}_{\text{new}})$ w.r.t. $q(\bm{b})$
guarantees that $q_{\text{new}}(f)$ gets closer to
$\hat{p}(f|\bm{y}_{\text{old}},\bm{y}_{\text{new}})$222We omitted the
subscript of $q_{\text{new}}(\bm{b})$. Setting the derivative of
$\mathbb{ELBO}(q(\bm{b}),\bm{\theta}_{\text{new}})$ w.r.t. $q(\bm{b})$ equal
to $0$ gives us the optimal new approximate distribution $q(\bm{b})^{*}$.
Also, $\mathbb{ELBO}(q(\bm{b}),\bm{\theta}_{\text{new}})$ lower bounds the
online log marginal likelihood333$Z_{2}/Z1\approx
p(\bm{y}_{\text{new}}|\bm{y}_{\text{old}})$ since the KL divergence is non-
negative. This lower bound can be used to learn the hyperparameters and
optimize pseudo inputs. This gives us a principled framework for deploying
sparse GP in the streaming setting, providing online hyperparameter learning
and pseudo-input optimization. See [5] for more details.

##  IV Experiments

Our experimental evaluations aim at answering the following questions about
SSGP:

  * (Q1)

Is SSGP able to achieve competitive accuracy with improved computational
complexity and memory usage?

  * (Q2)

Can it effectively characterize the environment by learning hyperparameters?

  * (Q3)

How does the number of pseudo points influence the accuracy and efficiency?

###  IV-A Experiment Setup

(Setting) We have conducted extensive evaluations with both synthetic and
real-world data, and compared the SSGP with the following baseline methods:

  1. 1.

standard GP regression (GPR) with the whole collected dataset [31],

  2. 2.

standard GP regression with the most recent $500$444We keep the last $500$
samples so that the runtime of GPR roughly matches that of SSGP (GPR500)
training data points,

  3. 3.

variational sparse GP (VSGP) with the whole dataset [34],

  4. 4.

sparse pseudo-inputs GP (SPGP) with the whole dataset [27, 33].

The covariance function used throughout the experiments is the squared
exponential (SE) kernel with learned lengthscales $\ell_{d}$ (as a
hyperparameter):

$k(\bm{x},\bm{x}^{\prime})=\sigma^{2}_{f}\exp\left[-\frac{1}{2}\sum_{d=1}^{D}\left(\frac{x_{d}-x^{\prime}_{d}}{\ell_{d}}\right)^{2}\right],$

where the number of input dimensions $D$ is $2$ in our case. The
hyperparameters and pseudo inputs are initialized with the same values, and
then optimized using L-BFGS-B with the same stopping criteria. All the methods
were implemented on GPflow [25] and run on a standard desktop with a 3.6GHz
Intel i7 processor and 16GB of RAM.

(Data) Both the synthetic data and the real-world data are simulated on a
$100\text{m}\times 100\text{m}$ grid map (discrete scalar field) with the grid
resolution of 1m in each dimension, thus the test set contains $10000$ data
points. The synthetic environment is drawn from a two-dimensional GP with
hyperparameters $\\{\sigma^{2}_{y}:0.01$, $\sigma^{2}_{f}:1$, $\ell_{1}:0.3$,
$\ell_{2}:0.7\\}$. We use the real field Sea Surface Temperature (SST) data
provided by National Oceanic and Atmospheric Administration (NOAA). We adopt a
naïve lawnmower sampling method to best demonstrate the idea and reduce the
impact of the sampling mechanism on the modeling and learning results. The
sampling path allows us to collect $4312$ training data points which are
observed sequentially along the path. Specifically, the sampling robot follows
the planned lawnmower path and gathers a small batch of data before updating
the GP model. In this way we can split the entire lawnmower path into $98$
small batches, each of which contains $44$ samples. All the GP models of
different batches will update their predictions, hyperparameters, and the
pseudo inputs after receiving each batch of data. The optimized
hyperparameters and pseudo inputs estimated from the previous step will be
used as initialization for the next round.

Figure 2: Evaluation results on synthetic data Figure 3: Evaluation results
with real-world data

(Metrics) To measure the learning performance, we use the Root Mean Squared
Error (RMSE) and the average Negative Log Predictive Density (NLPD) on the
test data:

$\displaystyle\text{RMSE}=$
$\displaystyle\sqrt{\frac{1}{T}\sum_{t=1}^{T}(f_{t}-m_{t})^{2}},$
$\displaystyle NLPD=$ $\displaystyle\frac{1}{T}\sum_{t=1}^{T}-\log
p(y_{t}|\bm{x}_{t}),$

where $T$ is the number of the test data points, $f_{t}$ and $y_{t}$ are the
underlying function value and its corresponding noisy observation,
respectively; $m_{t}$ is the posterior mean value, and $p(y_{t}|\bm{x}_{t})$
is the probability density of the predictive distribution. RMSE only takes a
point prediction into account while NLPD penalizes over-confident predictions
and under-confident ones. We also show the numbers of on-board data including
training data points and pseudo points to compare the memory usage. We are
revealing that, the training runtime, prediction runtime, and overall runtime
demonstrate the computational efficiency of the SSGP method.

Figure 4: The learned hyperparameters with the ground truth in orange. Figure
5: Relationship between the number of pseudo points and the modeling accuracy
and efficiency. The values of $\alpha$ are set to be $0.5,1,2,4$ in
$M=\lceil\alpha\log^{2}N+1\rceil$.

###  IV-B Results

We first evaluate the algorithms’ accuracy, practical runtimes, and memory
usage (Q1). Fig. 2 shows the results of performance evaluation. Although
GPR500 has the fastest runtime, it throws away useful historical information,
leading to low accuracy in terms of RMSE and NLPD. The SSGP, GPR, VSGP, and
SPGP have similar RMSEs555We moved RMSE and NLPD of SPGP down a bit for better
visualization., but the computational cost of SSGP is much less than the other
three. The overall runtime of SSGP is slightly more than that of GPR500.
(Note, the runtimes of GPR are plotted in log scale.)

The number of on-board data is a reflection of the memory usage. The SSGP only
needs to store the current batch of data and a small set of pseudo points,
hence it entails the lowest memory consumption. This feature is extremely
important for long-term autonomous environmental monitoring because the robot
has limited computation and memory resources while the data size continuously
grows.

As noted by [5], there is a discrepancy between the NLPD of SSGP and that of
the other three as more and more data have been collected. This phenomenon is
more obvious on the real-world data as shown in Fig. 3. In the next section,
we will show that using more pseudo points can mitigate this issue without
increasing much computational cost and memory consumption.

To investigate the effectiveness of environmental characterization by learning
hyperparameters (Q2), we visualize the learned hyperparameters of all the
algorithms together with the ground truth (Fig. 4). The SSGP is able to
capture the underlying amplitude ($\sigma^{2}_{f}$) and lengthscale along
x-axis, though it underestimates the lengthscale along y-axis. Through the
lengthscales of SSGP, we can also infer that the synthetic data varies rapidly
along the x direction. Some snapshots of the modeling process of both datasets
are shown in Fig. 6 and Fig. 7.

We then investigate the relationship between the number of pseudo points and
the accuracy of the SSGP (Q3). We summarize the training data using $M$ pseudo
points, reducing the computational cost to $\mathcal{O}(NM^{2})$. However, as
seen in Fig. 2 and Fig. 3, the quality of approximation will be limited if the
number of pseudo points is not sufficient. (See the increased performance
discrepancy between the SSGP and other methods along with the growth of data
points where the number of pseudo points keeps unchanged.) This implies that
the number of pseudo points $M$ implicitly depends on the number of training
data $N$. For example, if such dependence is linear, i.e., $M$ scales linearly
with $N$, then the computational complexity is still $\mathcal{O}(N^{3})$.
Fortunately, a recent work [7] shows that $M=\mathcal{O}(\log^{D}N)$ suffices
for regression with normally distributed inputs in $D$-dimensions with the SE
kernel. To better demonstrate the trends, we ran SSGP with
$M=\lceil\alpha\log^{2}N+1\rceil$ pseudo inputs and a set of differing
$\alpha$ values to control the number $M$, i.e., $\alpha={0.5,1,2,4}$. We then
plotted the RMSE and NLPD of GPR for comparison. As shown in Fig. 5, when
$\alpha$ is $0.5$ or $1$, there are significant gaps between the RMSE and NLPD
of SSGP and those of GPR. When $\alpha=2$, the performance of SSGP matches to
that of GPR whilst the computational times grow with a small margin. When
$\alpha=4$, it brings a slight improvement in precision at the cost of more
runtime.

##  V Conclusions

This paper presents an environmental model learning framework using streaming
sparse Gaussian process (SSGP). The SSGP updates the predictive distribution,
the hyperparameters, and the pseudo points which summarize the historical
data, in a streaming manner. We have evaluated the performances of SSGP by
comparing it against other baseline methods. Our evaluations show that the
SSGP produces competitive prediction accuracy with dramatically reduced
computational cost and memory demand, making it suitable for long-term
autonomous environmental monitoring. We have also empirically investigated how
the number of pseudo points influences the learning accuracy and efficiency.
Our experimental result reveals that $\lceil 2\log^{2}N+1\rceil$ pseudo points
are sufficient for regression with $N$ $2$D inputs and SE ARD kernel.

## References

  * [1] A. Arora, P. M. Furlong, R. Fitch, S. Sukkarieh, and T. Fong, “Multi-modal active perception for information gathering in science missions,” _Autonomous Robots_ , pp. 1–27, 2018.
  * [2] S. Bai, J. Wang, F. Chen, and B. Englot, “Information-theoretic exploration with bayesian optimization,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2016, pp. 1816–1822.
  * [3] G. Best, O. M. Cliff, T. Patten, R. R. Mettu, and R. Fitch, “Dec-mcts: Decentralized planning for multi-robot active perception,” _The International Journal of Robotics Research_ , vol. 38, no. 2-3, pp. 316–337, 2019\.
  * [4] J. Binney, A. Krause, and G. S. Sukhatme, “Informative path planning for an autonomous underwater vehicle,” in _2010 IEEE International Conference on Robotics and Automation_. IEEE, 2010, pp. 4791–4796.
  * [5] T. D. Bui, C. Nguyen, and R. E. Turner, “Streaming sparse gaussian process approximations,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 3299–3307.
  * [6] T. D. Bui, J. Yan, and R. E. Turner, “A unifying framework for gaussian process pseudo-point approximations using power expectation propagation,” _The Journal of Machine Learning Research_ , vol. 18, no. 1, pp. 3649–3720, 2017.
  * [7] D. Burt, C. E. Rasmussen, and M. Van Der Wilk, “Rates of convergence for sparse variational gaussian process regression,” in _International Conference on Machine Learning_ , 2019, pp. 862–871.
  * [8] N. Cao, K. H. Low, and J. M. Dolan, “Multi-robot informative path planning for active sensing of environmental phenomena: A tale of two algorithms,” in _Proceedings of the 2013 international conference on Autonomous agents and multi-agent systems_. International Foundation for Autonomous Agents and Multiagent Systems, 2013, pp. 7–14.
  * [9] W. Chen and L. Liu, “Pareto monte carlo tree search for multi-objective informative planning,” in _Proceedings of Robotics: Science and Systems_ , FreiburgimBreisgau, Germany, June 2019.
  * [10] J. J. Chung, N. R. Lawrance, and S. Sukkarieh, “Learning to soar: Resource-constrained exploration in reinforcement learning,” _The International Journal of Robotics Research_ , vol. 34, no. 2, pp. 158–172, 2015\.
  * [11] L. Csató and M. Opper, “Sparse on-line gaussian processes,” _Neural computation_ , vol. 14, no. 3, pp. 641–668, 2002.
  * [12] M. Dunbabin and L. Marques, “Robots for environmental monitoring: Significant advancements and applications,” _IEEE Robotics & Automation Magazine_, vol. 19, no. 1, pp. 24–39, 2012.
  * [13] G. Flaspohler, N. Roy, and Y. Girdhar, “Near-optimal irrevocable sample selection for periodic data streams with applications to marine robotics,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 1–8.
  * [14] G. Hitz, A. Gotovos, M.-É. Garneau, C. Pradalier, A. Krause, R. Y. Siegwart _et al._ , “Fully autonomous focused exploration for robotic environmental monitoring,” in _2014 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2014, pp. 2658–2664.
  * [15] G. A. Hollinger, B. Englot, F. S. Hover, U. Mitra, and G. S. Sukhatme, “Active planning for underwater inspection and the benefit of adaptivity,” _The International Journal of Robotics Research_ , vol. 32, no. 1, pp. 3–18, 2013.
  * [16] G. A. Hollinger and G. S. Sukhatme, “Sampling-based robotic information gathering algorithms,” _The International Journal of Robotics Research_ , vol. 33, no. 9, pp. 1271–1287, 2014.
  * [17] S. Kemna, H. Heiarsson, and G. S. Sukhatme, “On-board adaptive informative sampling for auvs: a feasibility study,” in _OCEANS 2018 MTS/IEEE Charleston_. IEEE, 2018, pp. 1–10.
  * [18] A. Krause, A. Singh, and C. Guestrin, “Near-optimal sensor placements in gaussian processes: Theory, efficient algorithms and empirical studies,” _Journal of Machine Learning Research_ , vol. 9, no. Feb, pp. 235–284, 2008\.
  * [19] X. Lan and M. Schwager, “Rapidly exploring random cycles: Persistent estimation of spatiotemporal fields with multiple sensing robots,” _IEEE Transactions on Robotics_ , vol. 32, no. 5, pp. 1230–1244, 2016.
  * [20] C. K. Ling, K. H. Low, and P. Jaillet, “Gaussian process planning with lipschitz continuous reward functions: Towards unifying bayesian optimization, active learning, and beyond,” in _Thirtieth AAAI Conference on Artificial Intelligence_ , 2016.
  * [21] K.-C. Ma, L. Liu, H. K. Heidarsson, and G. S. Sukhatme, “Data-driven learning and planning for environmental sampling,” _Journal of Field Robotics_ , vol. 35, no. 5, pp. 643–661, 2018.
  * [22] K.-C. Ma, L. Liu, and G. S. Sukhatme, “An information-driven and disturbance-aware planning method for long-term ocean monitoring,” in _Intelligent Robots and Systems (IROS), 2016 IEEE/RSJ International Conference on_. IEEE, 2016, pp. 2102–2108.
  * [23] ——, “Informative planning and online learning with sparse gaussian processes,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 4292–4298.
  * [24] R. Marchant and F. Ramos, “Bayesian optimisation for intelligent environmental monitoring,” in _2012 IEEE/RSJ international conference on intelligent robots and systems_. IEEE, 2012, pp. 2242–2249.
  * [25] D. G. Matthews, G. Alexander, M. Van Der Wilk, T. Nickson, K. Fujii, A. Boukouvalas, P. León-Villagrá, Z. Ghahramani, and J. Hensman, “Gpflow: A gaussian process library using tensorflow,” _The Journal of Machine Learning Research_ , vol. 18, no. 1, pp. 1299–1304, 2017.
  * [26] A. Meliou, A. Krause, C. Guestrin, and J. M. Hellerstein, “Nonmyopic informative path planning in spatio-temporal models,” in _AAAI_ , vol. 10, no. 4, 2007, pp. 16–7.
  * [27] R. Mishra, M. Chitre, and S. Swarup, “Online informative path planning using sparse gaussian processes,” in _2018 OCEANS-MTS/IEEE Kobe Techno-Oceans (OTO)_. IEEE, 2018, pp. 1–5.
  * [28] P. Morere, R. Marchant, and F. Ramos, “Sequential bayesian optimization as a pomdp for environment monitoring with uavs,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 6381–6388.
  * [29] M. Popović, G. Hitz, J. Nieto, I. Sa, R. Siegwart, and E. Galceran, “Online informative path planning for active classification using uavs,” in _2017 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2017, pp. 5753–5758.
  * [30] J. Quiñonero-Candela and C. E. Rasmussen, “A unifying view of sparse approximate gaussian process regression,” _Journal of Machine Learning Research_ , vol. 6, no. Dec, pp. 1939–1959, 2005.
  * [31] C. E. Rasmussen, “Gaussian processes in machine learning,” in _Summer School on Machine Learning_. Springer, 2003, pp. 63–71.
  * [32] A. Singh, A. Krause, C. Guestrin, and W. J. Kaiser, “Efficient informative sensing using multiple robots,” _Journal of Artificial Intelligence Research_ , vol. 34, pp. 707–755, 2009.
  * [33] E. Snelson and Z. Ghahramani, “Sparse gaussian processes using pseudo-inputs,” in _Advances in neural information processing systems_ , 2006, pp. 1257–1264.
  * [34] M. Titsias, “Variational learning of inducing variables in sparse gaussian processes,” in _Artificial Intelligence and Statistics_ , 2009, pp. 567–574.

Figure 6: Snapshots of simulation on the synthetic data. The first column
denotes the environmental ground truth and the robot’s lawnmower sampling
path. The 2nd to the 6th columns are the predictions of GPR, GPR500, SSGP,
VSGP, and SPGP. White croses are the pseudo inputs.  Figure 7: Snapshots of
simulation on the Sea Surface Temperature (SST) data from National Oceanic and
Atmospheric Administration (NOAA).