# Gaussian Process Inference Using Mini-batch Stochastic Gradient Descent:
Convergence Guarantees and Empirical Benefits

Hao Chen<EMAIL_ADDRESS>
Lili Zheng11footnotemark: 1<EMAIL_ADDRESS>
Department of Statistics
University of Wisconsin-Madison
1300 University Avenue
Madison, WI 53706, USA
Raed Al Kontar<EMAIL_ADDRESS>
Department of Industrial and Operations Engineering
University of Michigan
1891 IOE Building 1205, Beal Ave
Ann Arbor, MI 48109, USA
Garvesh Raskutti<EMAIL_ADDRESS>
Department of Statistics
University of Wisconsin-Madison
1300 University Avenue
Madison, WI 53706, USA Equal contribution.Corresponding author.

###### Abstract

Stochastic gradient descent (SGD) and its variants have established themselves
as the go-to algorithms for large-scale machine learning problems with
independent samples due to their generalization performance and intrinsic
computational advantage. However, the fact that the stochastic gradient is a
biased estimator of the full gradient with correlated samples has led to the
lack of theoretical understanding of how SGD behaves under correlated settings
and hindered its use in such cases. In this paper, we focus on hyperparameter
estimation for the Gaussian process (GP) and take a step forward towards
breaking the barrier by proving minibatch SGD converges to a critical point of
the full log-likelihood loss function, and recovers model hyperparameters with
rate $O(\frac{1}{K})$ for $K$ iterations, up to a statistical error term
depending on the minibatch size. Our theoretical guarantees hold provided that
the kernel functions exhibit exponential or polynomial eigendecay which is
satisfied by a wide range of kernels commonly used in GPs. Numerical studies
on both simulated and real datasets demonstrate that minibatch SGD has better
generalization over state-of-the-art GP methods while reducing the
computational burden and opening a new, previously unexplored, data size
regime for GPs.

Keywords: Stochastic Optimization, Gaussian Processes, Convergence Rate,
Scalability

##  1 Introduction

The Gaussian process (GP) has seen many success stories in various domains, be
it in optimization (Yue and Al Kontar, 2020; Snoek et al., 2012),
reinforcement learning (Srinivas et al., 2009; Krause and Ong, 2011), time
series analysis (Kontar et al., 2020; Álvarez and Lawrence, 2011), control
theory (Kocijan et al., 2004; Mesbah, 2016) and simulation meta-modeling (Zhou
et al., 2011; Qian and Wu, 2008). One can attribute such success to its
natural Bayesian interpretation, uncertainty quantification capability and
highly flexible model priors. Yet its main limitation is the $O(n^{3})$
computation and $O(n^{2})$ storage for $n$ training points (Rasmussen, 2003).
Indeed, as mentioned in Hensman et al. (2013), a traditional large dataset for
a GP is one with a few thousand data points and even those often require
approximation techniques.

As a result, in the past two decades, a large body of work on GPs tackled
approximate inference procedures to reduce the computational demands and
numerical instabilities (mainly due to the need for matrix inversions). This
push towards scalability dates back to the seminal paper by Quiñonero-Candela
and Rasmussen (2005) in 2005 which unified previous approximation methods into
a single probabilistic framework based on inducing points. Since then, many
new methods have also been introduced. Most notable are: variational inference
procedures that laid the theoretical foundation for the class of inducing
point methods (Damianou et al., 2016; Nguyen et al., 2014; Zhao and Sun, 2016;
Álvarez et al., 2010; Wilson et al., 2016), mixture of experts models
(Deisenroth and Ng, 2015; Tresp, 2000), covariance tapering (Furrer et al.,
2006; Kaufman et al., 2008) and kernel expansions (Le et al., 2013; Rahimi and
Recht, 2008; Yang et al., 2015). On the other hand, there has been a recent
push to utilize increasing computational power and GPU acceleration to solve
exact GPs. This recent literature includes distributed Cholesky factorizations
(Nguyen et al., 2019), preconditioned conjugate gradients (PCG) to solve
linear systems (Gardner et al., 2018) and kernel matrix partitioning to
perform all matrix-vector multiplications (Wang et al., 2019). Interestingly,
Wang et al. (2019) was able to fit a bit more than 1 million data points using
8 GPUs in a few days.

One possible solution to extend GPs far beyond what is currently possible is
through stochastic gradient decent (SGD) and its variants: _drawing $m<<n$
samples at each iteration and updating model parameters following the gradient
of the log-likelihood loss function on the $m$ subsamples_. Indeed, SGD, or
more generally the capability of inference via minibatches (possibly also with
second order information), has been a key propeller behind the success of deep
learning (LeCun et al., 2015) in its various forms and other objectives. The
caveat in GPs, however, is that, unlike empirical loss minimization, there
exists correlation across all samples where any finite collection of the
samples have a joint Gaussian distribution with covariance characterized by an
empirical kernel matrix. Hence the log-likelihood loss function is no longer
the sum of losses evaluated at each sample, translating to the stochastic
gradient being a biased estimator of the full gradient when taking expectation
with respect to the random sampling. The lack of theoretical backing and
understanding of how SGD behaves has long stood in the way of using SGD to
conduct inference in GPs (Hensman et al., 2013) and in most settings where
correlation amongst samples is high.

In this paper, we establish convergence guarantees of SGD for GPs for both the
full gradient and the model parameters. Interestingly, without convexity or
even Liptchitz conditions on the loss function, the structure of GP leads to
an optimization error term of $O(\frac{1}{K})$ for converging to a critical
point and recovering the true noise variance up to a statistical error that
vanishes as $m$ tends to $+\infty$, for both RBF kernels and Matérn kernels.
Our proof involves two key steps: first we concentrate the stochastic gradient
to its conditional expectation using a covering argument and then we show that
the latter satisfies a property similar to strong convexity by exploiting
eigenvalues of the empirical kernel matrix. The proof and key findings offer
standalone value beyond GPs and we hope they encourage researchers to further
investigate SGD in other correlated settings such as Lévy, Itô and Markov
processes.

Most importantly, our results open up a new data size regime to explore GPs.
We are able to train $n\approx 1.2\times 10^{6}$ data points using a single
CPU core in around 30 minutes. Recall, it took the most recent advancements in
exact GPs a couple of days using 8 GPUs to train when $n\approx 10^{6}$, and
$n$ is limited to approximately $10^{4}$ without GPU. We find that GPs
inferred using SGD offer remarkably better performance in various case studies
with different dataset sizes, noise levels and input dimensions.

###  1.1 Main Contributions

We establish convergence guarantees for the minibatch SGD algorithm when
training a GP under sampling with or without replacement and conduct numerical
experiments to validate and supplement our theoretical results. Our main
contributions are summarized as follows:

  * •

Convergence guarantees: For a large enough minibatch size $m$, minibatch SGD
converges to a critical point of the full log-likelihood loss function, and
recovers the true noise variance up to a statistical error depending on $m$,
when the kernel function exhibits exponential (RBF kernels) or polynomial
eigendecay (Matérn kernels). To be specific, the full gradient and the
estimation error of the noise variance evaluated at the $K$th iterate are
bounded by an optimization error term $O(\frac{1}{K})$ and a statistical error
term $O(m^{-\frac{1}{2}})$, see Theorems 3.1 to 3.4.

    * –

Proof techniques for statistical error: Since the stochastic gradient is
biased for estimating full gradients, we instead bound the difference from its
conditional expectation given the covariates in the corresponding minibatch,
uniformly over all possible parameter iterates. We use novel truncation and
covering arguments to prove the uniform error bound, in order to avoid the
dependence between past parameter iterates and the minibatch in the current
iteration. This contributes to the statistical error term
$O(m^{-\frac{1}{2}})$ in the convergence error bound.

    * –

Proof techniques for optimization error: To guarantee the $O(\frac{1}{K})$
optimization error bound, no convexity or even Liptchitz condition on the loss
function are assumed. Instead, we prove that the conditional expectation of
the loss function given covariates $\bm{\mathbf{X}}_{n}$ satisfies a relaxed
property of strong convexity (see e.g., Lemma 4), where the “curvature”
parameter is lower bounded by a constant regardless of minibatch size $m$.
This proof relies on careful analysis for bounding the eigenvalues of
empirical kernel matrices.

  * •

Numerical findings: Through benchmarking with state-of-the-art methods on
various datasets we show that SGD offers great value from both computational
and statistical perspectives. Computationally, we scale to dataset sizes
previously unexplored in GPs in a fraction of time needed for competing
methods. Meanwhile statistically, we find that SGD improves generalization in
GPs, specifically in large data settings.

###  1.2 Related Work

As mentioned earlier, there are several methods trying to tackle the
computational complexity of GPs. Those can be roughly split into the following
three categories, though it is by no means an exhaustive list (see the survey
in [1]).

  * •

Exact inference via matrix vector multiplications (MVM): This recent class of
literature has had the most success in scaling GPs. Initially such approaches
depended on a structured kernel matrix where data lies in a regularly spaced
grid (Saatçi, 2012; Wilson and Nickisch, 2015). Then with the help of GPU
acceleration, conjugate gradient and distributed Cholesky factorization, MVMs
were applied to more general settings (Wang et al., 2019; Gardner et al.,
2018; Ubaru et al., 2017). Such approaches have training complexity of
$O(n^{2})$ ($O(n\,\mbox{log}\,n)$ possible on spaced grids), yet amenable to
distributed computation and GPU acceleration.

  * •

Sparse approximate inference: This class of methods is based on a low rank
approximation of the empirical kernel matrix where
$\bm{\mathbf{K}}_{n}\approx\bm{\mathbf{K}}_{nz}\bm{\mathbf{K}}^{-1}_{zz}\bm{\mathbf{K}}_{zn}$
and $z$ denotes a set of inducing points with $\mbox{cardinality}(z)=n_{z}<<n$
(Kontar et al., 2018; Alvarez and Lawrence, 2009; Damianou et al., 2016; Zhao
and Sun, 2016; Snelson and Ghahramani, 2006). Their time complexity is mainly
$O(n_{z}^{2}n)$ which can be reduced to $O(n+cn_{z})$ for structured and
regularly spaced grids. Indeed, sparse GPs have gained increased attention
since variational inference (VI) laid the theoretical foundation of this class
of inducing points/kernel approximations (starting from the early work of
Titsias (2009)).

  * •

Stochastic variational inference (SVI): Following the work of Hoffman et al.
(2013), SVI was introduced to GPs in Hensman et al. (2013). The key idea is to
introduce a variational distribution over the inducing points so that the VI
framework is amenable to stochastic optimization. This leads to a complexity
of $O(n_{z}^{3})$ at each iteration (Hoang et al., 2015; Blei et al., 2017).
Unfortunately, recent results in Burt et al. (2019) show the need for at least
$O(\mbox{log}^{D}n)$ inducing points for Gaussian kernels, which implies a
superlinear growth with the input dimension $D$. Although many of the
aforementioned methods are proposed in the context of model inference
(prediction), the idea can often carry over to the model selection task, and
the complexity is multiplied by the number of iterations.

The theoretical analysis for SGD has also been extensively studied under
various assumptions (Nemirovski et al., 2009; Rakhlin et al., 2011; Frostig et
al., 2015; Bottou et al., 2018). In particular, in the context of empirical
risk minimization where the objective loss function is summed over $n$ data
points and stochastic gradients are calculated for i.i.d. sampled data at all
iterations, it is known that the expected squared error of SGD iterates at
iteration $K$ compared to the true minimizer is $O(1/K)$, with diminishing
step size for strongly convex objectives. Furthermore, some recent literature
(Hardt et al., 2016; Keskar et al., 2016) suggests that SGD has good
generalization power, which has encouraged more practitioners to apply SGD to
various application scenarios.

However, there are much fewer results when the stochastic gradients are biased
estimators for the full gradients, despite the fact that unbiased estimates
for the full gradient can be expensive or unavailable in certain cases. Some
examples include: learning graph neural networks (Chen et al., 2018),
distributed parallel optimization where sparsified stochastic gradient is
applied (Stich and Karimireddy, 2019) and performing model selection for GPs.
Homem-de Mello (2008); Chen and Luss (2018); Ajalloeian and Stich (2020) study
the stochastic gradient algorithms under non-i.i.d. sampling or when the
stochastic gradients are biased, and provide error bounds involving the bias
term, or convergence guarantees built on consistency assumptions. Our paper
does not make assumptions on the consistency of stochastic gradients or
convexity of the full loss function, but exploits the nature of GP loss
function and kernel matrices instead.

###  1.3 Organization

The paper is organized as follows: The problem setup is described in Section
2, and theoretical guarantees are provided in Section 3; Section 4 is devoted
to a proof outline, key lemmas and the proof of some main steps; Section 5
presents practical considerations for applying minibatch SGD on GPs, while
Section 6 includes our numerical results. We point out some open problems in
Section 7 and conclude in Section 8.

##  2 Problem Setup

#### Notation

Vectors and matrices are denoted by boldface letters, e.g.,
$\bm{\mathbf{K}}_{n}$, $\bm{\mathbf{\theta}}$, except for the full gradient
$\nabla\ell(\bm{\mathbf{\theta}})$ and stochastic gradient
$g(\bm{\mathbf{\theta}})$. For any vector $\bm{\mathbf{u}}\in\mathbb{R}^{p}$,
$u_{i}$ denotes its $i$th entry, and
$\|\bm{\mathbf{u}}\|_{2}=\left(\sum_{i=1}^{p}u_{i}^{2}\right)^{\frac{1}{2}}$
denotes its $\ell_{2}$ norm. For any square matrix $\bm{\mathbf{A}}$,
$\lambda_{i}(\bm{\mathbf{A}})$ denotes its $i$th largest eigenvalue. Also see
a table of important notations in Appendix A.

We consider the Gaussian process model

$\begin{split}f\sim\mathcal{GP}(\mu(\cdot),k(\cdot,\cdot)),&\quad\mathbf{x}_{1},\dots,\mathbf{x}_{n}\overset{\text{i.i.d.}}{\sim}\mathbb{P},\\\
y_{i}=f(\mathbf{x}_{i})+\epsilon_{i},&\quad\epsilon_{i}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\sigma_{\epsilon}^{2}),\quad
1\leq i\leq n,\end{split}$ (1)

where $\mathbf{x}_{i}\in\mathcal{X}\subset\mathbb{R}^{D}$ is the input,
$\mu(\cdot):\mathcal{X}\rightarrow\mathbb{R}$ is the prior mean function,
$k(\cdot,\cdot):\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ is the
prior covariance function, and $\epsilon_{i}$ is the observational noise with
variance $\sigma_{\epsilon}^{2}$. Without loss of generality, we consider
constant 0 mean function. In addition, the prior covariance function
$k(\cdot,\cdot)=\sigma_{f}^{2}k_{0}(\cdot,\cdot)$ involves a known kernel
function $k_{0}(\cdot,\cdot)$ and a signal variance parameter
$\sigma_{f}^{2}$. We observe data points
$\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{n}$ generated from (1) and organize them
into
$(\mathbf{X}_{n},\mathbf{y}_{n})=((\mathbf{x}_{1},\dotsc,\mathbf{x}_{n})^{\top},(y_{1},\dotsc,y_{n})^{\top})$,
from which we aim to learn the hyperparameters in order to predict outputs
from new inputs based on the posterior process.

Denote by
$\bm{\mathbf{\theta}}^{*}=(\sigma_{f}^{2},\sigma_{\epsilon}^{2})^{\top}\in\mathbb{R}^{2}$
the underlying hyperparameters to be determined, and for notational
convenience, we may also use $\theta_{1}^{*}$ to denote $\sigma_{f}^{2}$ and
$\theta_{2}^{*}$ to denote $\sigma_{\epsilon}^{2}$ in the following. One
direct approach to estimate $\bm{\mathbf{\theta}}^{*}$ is by applying gradient
descent to minimize the scaled negative log marginal likelihood function

$\begin{split}\ell(\bm{\mathbf{\theta}};\mathbf{X}_{n},\mathbf{y}_{n})&=-\frac{1}{n}\log
p(\mathbf{y}_{n}|\mathbf{X}_{n},\bm{\mathbf{\theta}})\\\
&=\frac{1}{2n}[\mathbf{y}_{n}^{\top}\mathbf{K}_{n}^{-1}(\bm{\mathbf{\theta}})\mathbf{y}_{n}+\log|\mathbf{K}_{n}(\bm{\mathbf{\theta}})|+n\log(2\pi)]\end{split}$
(2)

over $\bm{\mathbf{\theta}}\in(0,\infty)^{2}$, where
$\mathbf{K}_{n}(\bm{\mathbf{\theta}})=\theta_{1}\mathbf{K}_{f,n}+\theta_{2}\mathbf{I}_{n}\in\mathbb{R}^{n\times
n}$ is the marginal covariance matrix for noisy observations $\mathbf{y}_{n}$
given $\bm{\mathbf{X}}_{n}$ and $\mathbf{K}_{f,n}\in\mathbb{R}^{n\times n}$ is
the kernel matrix of $k_{0}$ evaluated at $\mathbf{X}_{n}$, i.e.
$(\bm{\mathbf{K}}_{f,n})_{i,j}=k_{0}(\bm{\mathbf{x}}_{i},\bm{\mathbf{x}}_{j})$.
For notational convenience we will omit
$\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}})$ to $\bm{\mathbf{K}}_{n}$ when
$\bm{\mathbf{\theta}}$ is clear from the context and denote
$\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}}^{*})$ by $\bm{\mathbf{K}}_{n}^{*}$.
In this case, the derivative of $\ell(\bm{\mathbf{\theta}})$ is of particular
interest to us where each of its entries takes the form

$\begin{split}\left(\nabla\ell(\bm{\mathbf{\theta}};\mathbf{X}_{n},\mathbf{y}_{n})\right)_{l}&=\frac{1}{2n}\left[-\mathbf{y}_{n}^{\top}\mathbf{K}_{n}^{-1}\frac{\partial\mathbf{K}_{n}}{\partial\theta_{l}}\mathbf{K}_{n}^{-1}\mathbf{y}_{n}+\text{tr}\left(\mathbf{K}_{n}^{-1}\frac{\partial\mathbf{K}_{n}}{\partial\theta_{l}}\right)\right]\\\
&=\frac{1}{2n}\text{tr}\left[(\mathbf{K}_{n}^{-1}(\mathbf{I}_{n}-\mathbf{y}_{n}\mathbf{y}_{n}^{T}\mathbf{K}_{n}^{-1})\frac{\partial\mathbf{K}_{n}}{\partial\theta_{l}}\right],\end{split}$
(3)

where $\theta_{l}$ is the $l$th element of $\bm{\mathbf{\theta}}$ and
$(\partial\mathbf{K}_{n}/\partial\theta_{l})_{jk}=\partial(\mathbf{K}_{n})_{jk}/\partial\theta_{l}$.
For notational convenience we will suppress
$\bm{\mathbf{X}}_{n},\bm{\mathbf{y}}_{n}$ and use
$\nabla\ell(\bm{\mathbf{\theta}})$ instead. Notice that the computation in (3)
is dominated by the calculation of $\mathbf{K}_{n}^{-1}$ which requires
$O(n^{3})$ time. In order to reduce the computational cost of training, we
consider the minibatch stochastic gradient descent approach to optimize (2).

###  2.1 Minibatch SGD algorithm

Let $\xi$ be a random subset of $\\{i\\}_{i=1}^{n}$ of size $m$, then
$\\{(\mathbf{x}_{i},y_{i})\\}_{i\in\xi}$ is the corresponding subset of data
points which we organize into $(\mathbf{X}_{\xi},\mathbf{y}_{\xi})$, where
$\mathbf{X}_{\xi}$ is the submatrix formed by the rows of $\mathbf{X}_{n}$ and
$\mathbf{y}_{\xi}$ is the subvector of $\mathbf{y}_{n}$, both indexed by
$\xi$. Define
$g(\bm{\mathbf{\theta}};\mathbf{X}_{\xi},\mathbf{y}_{\xi})\in\mathbb{R}^{2}$
as an approximation to $\nabla\ell(\theta;\mathbf{X}_{n},\mathbf{y}_{n})$ that
can be calculated from this subset, i.e.

$\left(g(\bm{\mathbf{\theta}};\mathbf{X}_{\xi},\mathbf{y}_{\xi})\right)_{l}=\frac{1}{2s_{l}(m)}\text{tr}\left[(\mathbf{K}_{\xi}^{-1}(\mathbf{I}_{m}-\mathbf{y}_{\xi}\mathbf{y}_{\xi}^{\top}\mathbf{K}_{\xi}^{-1})\frac{\partial\mathbf{K}_{\xi}}{\partial\theta_{l}}\right],\quad
1\leq l\leq 2,$ (4)

where $\mathbf{K}_{\xi}$ is the covariance matrix for $\mathbf{y}_{\xi}$ while
also the principle submatrix formed by the rows and columns of
$\mathbf{K}_{n}$ indexed by $\xi$. In the following we will also let
$\bm{\mathbf{K}}_{f,\xi}$ denote the $m\times m$ block of
$\bm{\mathbf{K}}_{f,n}$ indexed by $\xi$. A natural choice for $s_{l}(m)$ is
$m$, but we will see in Section 3 that if kernels $k_{0}$ have exponential
eigendecay, setting $s_{1}(m)\asymp\log m$ and $s_{2}(m)=m$ would lead to the
convergence of both $\theta^{(k)}_{1},\theta^{(k)}_{2}$ to the true
hyperparameters. Algorithm 1 summarizes the steps of minibatch SGD, where we
do not specify whether minibatches are sampled with or without replacement
since our convergence guarantees will hold true under both cases, if minibatch
size $m$ is large enough (details provided in Section 3). We consider
_diminishing step sizes_ : the step size at the $k$ th iteration is
$\alpha_{k}=\frac{\alpha_{1}}{k}$. It is noteworthy that the time complexity
of Algorithm 1 is $O(Km^{3})$ , compared to $O(Kn^{3})$ for running gradient
descent with $K$ iterations.

###  2.2 Sampling Methods

In Algorithm 1 we conduct uniform sampling for each minibatch, that is, any
subset of indices $\xi\subset[n]$ of size $m$ has the same probability of
being selected. An alternative to uniform sampling is to sample data points
that are close to each other, which we call nearby sampling. One particular
nearby sampling strategy is nearest neighbor search, where a minibatch
consists of a uniformly sampled data point and its $m-1$ nearest neighbors
within the data pool. We may construct a $k$-$d$ tree to conduct nearest
neighbor search, which finds the $m-1$ nearest neighbors for every data point
in a given dataset of size $n$ in $O(n\log n)$ time and $O(n)$ space. That is
to say, the time complexity for minibatch SGD with this nearby sampling method
(only line 3 in Algorithm 1 changes) is $O(Km^{3}+n\log n)$ for $K$
iterations.

Our main theoretical contribution is establishing convergence guarantees for
uniform sampling SGD in Algorithm 1, but in addition to that, we will also
provide some theoretical insights and numerical experiments for understanding
the effect of nearby sampling.

1 Input: $\bm{\mathbf{\theta}}^{(0)}\in\mathbb{R}^{2}$, initial step size
$\alpha_{1}>0$.

2 for  _$k=1,2,\dotsc,K$_ do

3 Randomly sample a subset of indices $\xi_{k}$ of size $m$;

4 Compute the stochastic gradient
$g(\bm{\mathbf{\theta}}^{(k-1)};\mathbf{X}_{\xi_{k}},\mathbf{y}_{\xi_{k}})$;

5 $\alpha_{k}\leftarrow\frac{\alpha_{1}}{k}$;

6
$\bm{\mathbf{\theta}}^{(k)}\leftarrow\bm{\mathbf{\theta}}^{(k-1)}-\alpha_{k}g(\bm{\mathbf{\theta}}^{(k-1)};\mathbf{X}_{\xi_{k}},\mathbf{y}_{\xi_{k}})$;

7

8 end for

Algorithm 1 Minibatch SGD with uniform sampling

##  3 Theoretical Guarantees

In this section, we present convergence guarantees for Algorithm 1, including
error bounds for $\theta_{l}^{(k)}-\theta_{l}^{*},l=1,2$ and
$\nabla\ell(\bm{\mathbf{\theta}}^{(k)})$. Two types of kernels are considered:
those with exponential eigendecay (Section 3.1) and those with polynomial
eigendecay (Section 3.2). For both types of kernels, the convergence of
$\theta^{(k)}_{2}$ to the true noise variance $\sigma_{\epsilon}^{2}$ and the
full gradient $\nabla\ell(\bm{\mathbf{\theta}}^{(k)})$ to $0$ is guaranteed.
In particular, for the former type of kernel, $\theta^{(k)}_{1}$ is also
guaranteed to converge to the truth $\sigma_{f}^{2}$ under appropriate choice
of $s_{1}(m)$.

Furthermore, to understand the faster convergence observed in numerical
studies when nearby sampling is applied rather than uniform sampling
(Algorithm 1), we provide theoretical insights into how nearby sampling
influences the curvature parameter in Section 3.3.

First we state the assumptions needed for our convergence guarantees.

######  Assumption 3.1 (Bounded iterates)

Both $\bm{\mathbf{\theta}}^{*}$ and $\bm{\mathbf{\theta}}^{(k)}$ for $0\leq
k\leq K$ lie in $[\theta_{\min},\theta_{\max}]^{2}$, where
$0<\theta_{\min}<\theta_{\max}$.

######  Remark 1 (Justification for Assumption 3.1)

The boundedness of parameter iterates is usually assumed in the literature of
theoretical analysis for SGD (e.g., Nemirovski et al., 2009). As will be
revealed by Theorem 3.1 and Theorem 3.4, $\bm{\mathbf{\theta}}^{(k)}$ is
guaranteed to be bounded within a region around $\bm{\mathbf{\theta}}^{*}$
that gets smaller as $k$ increases w.h.p., whenever the previous
$\bm{\mathbf{\theta}}^{(0)},\dots,\bm{\mathbf{\theta}}^{(k-1)}$ are bounded
within some $[\theta_{\min},\theta_{\max}]$. Hence, we only need to be careful
about the initial steps of Algorithm 1, ensuring that the parameters
$\theta^{(k)}_{j},j=1,2$ are always positive and bounded. Moreover, our
numerical experiments suggest that the iterate $\bm{\mathbf{\theta}}^{(k)}$ of
Algorithm 1 is always bounded as long as the initial step size is chosen
appropriately (see Figure 2 in Section 6).

######  Assumption 3.2 (Bounded stochastic gradient)

For all $0\leq k<K$,

$\|g(\bm{\mathbf{\theta}}^{(k)};\bm{\mathbf{X}}_{\xi_{k+1}},\bm{\mathbf{y}}_{\xi_{k+1}})\|_{2}\leq
G$

for some $G>0$.

######  Remark 2 (Justification for Assumption 3.2)

The boundedness assumption for stochastic gradients is also commonly seen in
the literature (e.g., Hazan and Kale, 2011). Furthermore, Assumption 3.1 can
imply Assumption 3.2 with high probability, under similar conditions to those
in Theorem 3.1 or Theorem 3.3. The key idea is the stochastic gradients will
be shown to be close to their conditional expectation
$\mathbb{E}(g(\bm{\mathbf{\theta}}^{(k)};\bm{\mathbf{X}}_{\xi_{k+1}},\bm{\mathbf{y}}_{\xi_{k+1}}|\bm{\mathbf{X}}_{\xi_{k+1}})$,
and the latter is bounded with high probability given Assumption 3.1 and the
eigendecay assumptions for kernels in Theorem 3.1 or Theorem 3.3. We defer the
detailed explanation on this in Appendix F.

###  3.1 Kernels with Exponential Eigendecay

The exponential eigendecay assumption is stated in detail as follows:

######  Assumption 3.3 (Exponential eigendecay)

Consider the kernel operator $\mathcal{K}:L^{2}(\mathbb{P})\rightarrow
L^{2}(\mathbb{P})$ that satisfies
$\mathcal{K}\phi(\cdot)=\int\phi(\bm{\mathbf{x}})k_{0}(\cdot,\bm{\mathbf{x}})d\mathbb{P}(\bm{\mathbf{x}})$,
where $\mathbb{P}$ is the probability measure of the input as defined in (1).
The eigenvalues of $\mathcal{K}$ are $\\{Ce^{-bj}\\}_{j=0}^{\infty}$, where
$b>0$ and $C\leq 1$ are regarded as constants.

The exponential eigendecay assumption is satisfied by the Radial Basis
Function (RBF) kernels
$k(\bm{\mathbf{x}},\bm{\mathbf{x}}^{\prime})=\exp\\{-\|\bm{\mathbf{x}}-\bm{\mathbf{x}}^{\prime}\|_{2}^{2}/(2l^{2})\\}$
when the probability distribution $\mathbb{P}$ of the input is Gaussian (see
Section 4.3.1 of Rasmussen (2003)), which is widely seen in the GP literature.
The specific decay rate $b$ depends on the lengthscale parameter $l$ of the
corresponding kernel $k_{0}$. The requirement $C\leq 1$ is only for
theoretical convenience, and it suffices to have bounded $C$. The following
theorem guarantees the convergence of the parameter iterates under the
aforementioned assumptions.

######  Theorem 3.1 (Convergence of parameter iterates, exponential
eigendecay)

Consider the output $\bm{\mathbf{\theta}}^{(K)}$ of Algorithm 1, the minibatch
SGD algorithm with diminishing step sizes. Under Assumptions 3.1 to 3.3, when
$m>C$ for some constant $C>0$, we have the following results under two
corresponding conditions on $s_{l}(m)$:

  1. 1.

If $s_{2}(m)=m$, initial step size $\alpha_{1}$ satisfies
$\frac{3}{2\gamma}\leq\alpha_{1}\leq\frac{2}{\gamma}$ where
$\gamma=\frac{1}{4\theta_{\max}^{2}}$, then for any
$0<\varepsilon<C\frac{\log\log m}{\log m}$, with probability at least
$1-CK\exp\\{-cm^{2\varepsilon}\\}$,

$(\theta^{(K)}_{2}-\theta^{*}_{2})^{2}\leq\frac{8G^{2}}{\gamma^{2}(K+1)}+Cm^{-\frac{1}{2}+\varepsilon}.$
(5)

  2. 2.

If in addition to $s_{2}(m)=m$, $s_{1}(m)$ is set as $\tau\log m$ where
$\tau>C$, $\frac{3}{2\gamma}\leq\alpha_{1}\leq\frac{2}{\gamma}$ where $\gamma$
depends on $\tau$, then for any $0<\varepsilon<\frac{1}{2}$, with probability
at least $1-CK\exp\\{-c(\log m)^{2\varepsilon}\\}$,

$(\theta^{(K)}_{1}-\theta^{*}_{1})^{2}+(\theta^{(K)}_{2}-\theta^{*}_{2})^{2}\leq\frac{8G^{2}}{\gamma^{2}(K+1)}+C(\log
m)^{-\frac{1}{2}+\varepsilon}.$ (6)

Here $c,C>0$ depend only on $\theta_{\min},\theta_{\max},b$.

###### Remark 3.1

Theorem 3.1 suggests that the noise variance parameter $\theta^{(K)}_{2}$ is
guaranteed to converge to the truth $\theta^{*}_{2}$, with the optimization
error term $O(\frac{1}{K})$ and the statistical error term
$O(m^{-\frac{1}{2}+\varepsilon})$ with high probability if $\varepsilon\log m$
is large, when the initial stepsize is appropriately chosen and $s_{2}(m)=m$.
Furthermore, if we let $s_{1}(m)=\tau\log m$, then Algorithm 1 achieves
convergence for both $\theta^{(K)}_{1}$ and $\theta^{(K)}_{2}$ with
statistical error $O((\log m)^{-\frac{1}{2}+\varepsilon})$.

###### Remark 3.2

The different rates of statistical errors for estimating $\theta^{*}_{1}$ and
$\theta^{*}_{2}$ come from the different eigenvalue structures between
$\bm{\mathbf{K}}_{f,\xi}$ (the $m\times m$ block of $\bm{\mathbf{K}}_{f,n}$
indexed by $\xi$) and $\bm{\mathbf{I}}_{m}$. One may also note that the
statistical errors depend on $m$ instead of $n$: this is due to the
correlation among $\bm{\mathbf{y}}_{\xi}$ from different minibatches $\xi$,
conditioning on $\bm{\mathbf{X}}$, which is different from the problems with
independent samples.

###### Remark 3.3

The choice $s_{1}(m)\asymp\log m$ is because

$s_{1}(m)g_{1}(\bm{\mathbf{\theta}})=\mathrm{tr}\left[\bm{\mathbf{K}}_{\xi}^{-1}(\bm{\mathbf{I}}_{m}-\bm{\mathbf{y}}_{\xi}\bm{\mathbf{y}}_{\xi}^{\top}\bm{\mathbf{K}}_{\xi}^{-1})\frac{\partial\bm{\mathbf{K}}_{\xi}}{\partial\theta_{1}}\right]\asymp\log
m,$ (7)

and thus this choice of $s_{1}(m)$ ensures that $g_{1}(\bm{\mathbf{\theta}})$
has the same scale as $g_{2}(\bm{\mathbf{\theta}})$ (constant scale).

###### Remark 3.4

For the second case where $s_{1}(m)=\tau\log m$, we need
$\tau>\frac{64\theta_{\max}^{4}}{b\theta_{\min}^{4}}$ and

$\gamma=\min\left\\{\frac{1}{32\tau
b\theta_{\max}^{2}},\frac{1}{4\theta_{\max}^{2}}-\frac{2\theta_{\max}^{2}}{\tau
b\theta_{\min}^{4}}\right\\}.$

###### Remark 3.5

The optimization error $O(\frac{1}{K})$ is credited to the structure of the GP
loss function, which satisfies a relaxation of strong convexity (see Lemma 4
in Section 4). $\gamma$ can be viewed a lower bound of an approximate
“curvature” of the loss function, in the sense of a relaxed convexity. We will
revisit this curvature term in Section 3.3 and illustrate the potential
improvement nearby sampling brings to the curvature.

Based on Theorem 3.1, we also derive the following convergence guarantee for
the full gradient.

######  Theorem 3.2 (Convergence of full gradient, exponential eigendecay)

Consider the output $\bm{\mathbf{\theta}}^{(K)}$ of Algorithm 1, the minibatch
SGD algorithm with diminishing step sizes. Under Assumptions 3.1 to 3.3, if
$\frac{3}{2\gamma}\leq\alpha_{1}\leq\frac{2}{\gamma}$ for
$\gamma=\frac{1}{4\theta_{\max}^{2}}$, $m>C$, $s_{2}(m)=m$, then for any
$0<\varepsilon<C\frac{\log\log m}{\log m}$, with probability at least
$1-CK\exp\\{-cm^{2\varepsilon}\\}$,

$\|\nabla\ell(\bm{\mathbf{\theta}}^{(K)})\|_{2}^{2}\leq
C\left[\frac{G^{2}}{K+1}+m^{-\frac{1}{2}+\varepsilon}\right],$ (8)

holds, where $c,C>0$ depend only on $\theta_{\min},\theta_{\max},b$.

Theorem 3.2 implies that running SGD for sufficiently many iterations with
large minibatch size leads to the convergence to a critical point of
$\ell(\bm{\mathbf{\theta}})$. As we will show in the proof sketch, the error
bound for the full gradient is dominated by
$(\theta^{(K)}_{2}-\theta^{*}_{2})^{2}$, the estimation error of the noise
variance, thus it scales the same as (5).

###  3.2 Kernels with Polynomial Eigendecay

Now we consider the kernels with polynomial eigendecay, which captures much
stronger correlation than the kernels with exponential eigendecay, and thus
broadens the applications of GP to a wider class of data sets. Due to this
reason, it is of both practical and theoretical interest to investigate how
SGD performs for this type of kernels. The polynomial eigendecay assumption is
stated in detail as follows:

######  Assumption 3.4 (Polynomial eigendecay)

Consider the kernel operator $\mathcal{K}:L^{2}(\mathbb{P})\rightarrow
L^{2}(\mathbb{P})$ that satisfies
$\mathcal{K}\phi(\cdot)=\int\phi(\bm{\mathbf{x}})k_{0}(\cdot,\bm{\mathbf{x}})d\mathbb{P}(\bm{\mathbf{x}})$,
where $\mathbb{P}$ is the probability measure of the input as defined in (1).
The eigenvalues of $\mathcal{K}$ are $\\{Cj^{-2b}\\}_{j=0}^{\infty}$, where
$b>\frac{\sqrt{21}+3}{4}$, and $C\leq 1$ are regarded as constants.

This assumption is satisfied by the Matérn kernels (see section 2.3 in Bach
(2017), Kanagawa et al. (2018)), another important kernel function class
widely used in GP:

$k_{\alpha,h}(\bm{\mathbf{x}},\bm{\mathbf{x}}^{\prime})=\frac{1}{2^{\alpha-1}\Gamma(\alpha)}\left(\frac{\sqrt{2\alpha\|\bm{\mathbf{x}}-\bm{\mathbf{x}}\|_{2}^{2}}}{h}\right)^{\alpha}B_{\alpha}\left(\frac{\sqrt{2\alpha\|\bm{\mathbf{x}}-\bm{\mathbf{x}}\|_{2}^{2}}}{h}\right),$
(9)

where $B_{\alpha}(\cdot)$ is the modified Bessel function of the second kind
of order $\alpha$, and larger $\alpha$ leads to faster decay rate $b>0$.

######  Theorem 3.3 (Convergence of parameter iterates, polynomial eigendecay
)

Consider the output $\bm{\mathbf{\theta}}^{(K)}$ of Algorithm 1, the minibatch
SGD algorithm with diminishing step sizes. Under Assumptions 3.1 to 3.2 and
Assumption 3.4, when $m>C$ for some constant $C>0$, $s_{2}(m)=m$,
$\frac{3}{2\gamma}\leq\alpha_{1}\leq\frac{2}{\gamma}$ where
$\gamma=\frac{1}{8\theta_{\max}^{2}}$, then for any
$\varepsilon\in(\max\\{0,f_{1}(b)\\},\frac{1}{2})$, with probability at least
$1-CKm^{-f_{2}(b)\left[\varepsilon-
f_{1}(b)\right]}-CK\exp\\{-cm^{2\varepsilon}\\}$,

$(\theta^{(K)}_{2}-\theta^{*}_{2})^{2}\leq\frac{8G^{2}}{\gamma^{2}(K+1)}+Cm^{-\frac{1}{2}+\varepsilon}.$
(10)

Here $c,C>0$ depend only on $\theta_{\min},\theta_{\max},b$, and
$f_{1}(b)=-\frac{2b^{2}-5b-3}{2b(2b-1)}$, $f_{2}(b)=\frac{4b(2b-1)}{4b+3}$.

######  Remark 3.6 (Comparison with Theorem 3.1: error bounds)

Compared with Theorem 3.1, Theorem 3.3 reflects the influence of stronger
correlations (slower eigendecay of kernels) on the convergence of SGD. More
specifically, $\varepsilon>f_{1}(b)$ is required in addition to
$0<\varepsilon<\frac{1}{2}$. That is to say, when $0<f_{1}(b)<\frac{1}{2}$
($\frac{\sqrt{21}+3}{4}<b<3$), the statistical error scales at least as
$m^{-(\frac{1}{2}-f_{1}(b))}>m^{-\frac{1}{2}}$, which decreases as $b$
increases on $(\frac{\sqrt{21}+3}{4},3)$; while if $f_{1}(b)\leq 0$ ($b\geq
3$) which means the correlation is not too strong, the statistical error still
scales roughly as $m^{-\frac{1}{2}}$, the same as the exponential eigendecay
case. Therefore, the slower eigendecay of kernels (stronger correlation
structure) may lead to a slower convergence of SGD; while the good news is
that, for moderately fast polynomial eigendecay ($b\geq 3$) we still have the
same rate as the exponential eigendecay case.

######  Remark 3.7 (Comparison with Theorem 3.1: probability terms)

Another difference between Theorem 3.1 and Theorem 3.3 lies in the probability
term. When $b<3$ and thus $f_{1}(b)>0$, $CKm^{-f_{2}(b)[\varepsilon-
f_{1}(b)]}$ always dominates $CK\exp\\{-cm^{2\varepsilon}\\}$ for
$\varepsilon>f_{1}(b)$, which is a lower probability for the error bound in
(10) compared to (5). This is also the price we need to pay when considering
kernels with slower eigendecay.

To extend the theoretical results for kernels with exponential eigendecay to
polynomial eigendecay, we develop novel upper and lower bounds for
$\sum_{j=1}^{n}\lambda_{j}^{l}(\theta_{1}^{(k)}\lambda_{j}+\theta^{(k)}_{2})^{-2},l=0,1,2$
where $\lambda_{j}$ is the $j$th largest eigenvalue of $K_{f,\xi}$, see Lemma
12. The proof for Lemma 12 requires careful analysis and different arguments
from the proof for Lemma 11 that is established for kernels with exponential
eigendecay, although they are both based on error bounds for eigenvalues of
empirical kernel matrices in Braun (2006).

We briefly explain the reason behind the different scalings of statistical
errors between Theorem 3.3 and the first part of Theorem 3.1 in the following.
In fact, the statistical error term for
$(\theta^{(K)}_{2}-\theta^{*}_{2})^{2}$ is composed of two parts: one is
$m^{-\frac{1}{2}+\varepsilon}$ for both types of kernels, another is caused by
the fact that $\theta^{(K)}_{1}$ may not be estimated well, and the error
induced by this fact depends on the eigendecay of the kernel, which scales as
$m^{-(\frac{1}{2}-f_{1}(b))}$ for kernels with polynomial eigendecay. While
for kernels with exponential eigendecay, this error term scales as $\frac{\log
m}{m}$ and is dominated by the first error term
$m^{-\frac{1}{2}+\varepsilon}$.

###### Remark 3.8

For kernels with polynomial eigendecay, we don’t have convergence guarantee
for $\theta^{(K)}_{1}$ (signal variance of the kernel with the slowest
eigendecay). This is due to that it is very hard to derive matched upper and
lower bounds for the stochastic gradient
$s_{1}(m)g_{1}(\bm{\mathbf{\theta}})=\mathrm{tr}\left[(\bm{\mathbf{K}}_{\xi}^{-1}(\bm{\mathbf{I}}_{m}-\bm{\mathbf{y}}_{\xi}\bm{\mathbf{y}}_{\xi}^{\top}\bm{\mathbf{K}}_{\xi}^{-1})\frac{\partial\bm{\mathbf{K}}_{\xi}}{\partial\theta_{1}}\right]$
(which scales as $\log m$ for exponential eigendecay), and thus we cannot
specify the choice for $s_{1}(m)$ in order to make
$g_{1}(\bm{\mathbf{\theta}})$ scales similarly from
$g_{2}(\bm{\mathbf{\theta}})$.

######  Theorem 3.4 (Convergence of full gradient, polynomial eigendecay)

Consider the output $\bm{\mathbf{\theta}}^{(K)}$ of Algorithm 1, the minibatch
SGD algorithm with diminishing step sizes. Under the same conditions as
Theorem 3.4, for any $\varepsilon\in(\max\\{0,f_{1}(b)\\},\frac{1}{2})$, with
probability at least $1-CK\left(m^{-f_{2}(b)\left[\varepsilon-
f_{1}(b)\right]}+\exp\\{-cm^{2\varepsilon}\\}\right)$,

$\|\nabla\ell(\bm{\mathbf{\theta}}^{(K)})\|_{2}^{2}\leq
C\left[\frac{G^{2}}{K+1}+m^{-\frac{1}{2}+\varepsilon}\right],$ (11)

holds, where $c,C>0$ depend only on $\theta_{\min},\theta_{\max},b$,
$f_{1}(b)$ and $f_{2}(b)$ are defined as in Theorem 3.3.

As mentioned after Theorem 3.2, the bound (11) for the full gradient scales
the same as the bound (10) for $(\theta^{(K)}_{2}-\theta^{*}_{2})^{2}$.

###  3.3 Effect of Nearby Sampling

Although our theoretical guarantees are all derived for uniform sampling, some
empirical evidence suggests that sampling nearby points for each minibatch can
lead to lower errors for learning the noise variance
$\theta_{2}^{*}=\sigma_{\epsilon}^{2}$ and sometimes improved prediction
performance (some comparisons are provided in Section 6.2). In this section,
we present some theoretical insights to understand why and how nearby sampling
helps, in terms of learning $\theta_{2}^{*}$. Note that we focus on
$\theta_{2}^{*}$ instead of $\theta_{1}^{*}$, since
$\theta^{(k)}_{2}-\theta_{2}^{*}$ dominates the convergence of the upper bound
for $\|\nabla\ell(\bm{\mathbf{\theta}}^{(k)}\|_{2}$ at the $k$th iteration.

In the following, we investigate the following hypothesis: the approximate
“curvature” for $\theta_{2}^{(k)}$ is improved when the points within each
minibatch are closer to each other, so that nearby sampling leads to larger
curvature and hence faster convergence. More specifically, some calculation
shows that at the $(k+1)$th iteration, the “approximate” curvature term w.r.t.
the noise variance is

$\gamma(\bm{\mathbf{\theta}}^{(k)}):=\frac{\partial\mathbb{E}(g_{2}(\bm{\mathbf{\theta}}^{(k)};\bm{\mathbf{X}}_{\xi_{k+1}},\bm{\mathbf{y}}_{\xi_{k+1}})|\bm{\mathbf{X}}_{\xi_{k+1}})}{\partial\theta_{2}}=\frac{1}{2m}\sum_{j=1}^{m}(\theta^{(k)}_{1}\lambda_{j}^{(k)}+\theta^{(k)}_{2})^{-2},$
(12)

where $\lambda_{j}^{(k)}$ is the $j$th largest eigenvalue of
$\bm{\mathbf{K}}_{f,\xi_{k+1}}$. As revealed in our proof of Theorem 3.1, 3.4,
our convergence rate results critically depends on a lower bound for (12).

Figure 1: Curvature $\gamma(\bm{\mathbf{\theta}})$, defined in (12) under
different minibatch sizes $m$ and sampling schemes. For each minibatch size
$m$ and sampling scheme, the mean value of $\gamma(\bm{\mathbf{\theta}})$ over
$50$ replicates is taken and the standard deviation is marked as the error
bar.

#### Numerical evidence:

Figure 1 demonstrates an example of the “approximate” curvature term (12) with
various minibatch sizes $m$ under the uniform sampling and nearby sampling
schemes, where nearby sampling leads to larger curvatures.

The detailed numerical experiments for generating Figure 1 is as follows. We
randomly generate a full data pool including $n=2048$ data points:
$x_{1},\dots,x_{n}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{N}(0,10^{2})$. For
each minibatch size $m$, we perform uniform sampling and nearby sampling for
$50$ replicates; for uniform sampling the minibatch $\xi$ of size $m$ is
sampled from $1,\dots,n$ uniformly at random; for nearby sampling, the first
index $\xi_{1}$ in the minibatch $\xi$ is sampled uniformly at random from
$1,\dots,n$, and the rest $m-1$ indices correspond to the nearest neighbors of
$\bm{\mathbf{x}}_{\xi_{1}}$. Then for each replicate, we calculate
$\gamma(\bm{\mathbf{\theta}})=\frac{1}{2m}\sum_{j=1}^{m}(\theta_{1}\lambda_{j}+\theta_{2})^{-2}$,
where $\bm{\mathbf{\theta}}=(4,1)^{\top}$ and $\lambda_{j}$ is the $j$th
eigenvalue of $\bm{\mathbf{K}}_{f,\xi}$, the kernel matrix formed by
$\bm{\mathbf{x}}_{i},i\in\xi$. The kernel function $k_{0}(\cdot)$ is set as
the RBF kernel with lengthscale $l=0.5$. We then take the mean of
$\gamma(\bm{\mathbf{\theta}})$ over the $50$ replicates.

#### Theoretical insights:

Now we provide more theoretical insights into the influence of nearby sampling
upon $\gamma(\bm{\mathbf{\theta}})$. For simplicity of analysis, we fix our
focus on the $D=1$ case, and let the kernel $k_{0}$ be the RBF kernel:
$k_{0}(x,x^{\prime})=\exp\\{-(x-x^{\prime})^{2}/(2l^{2})\\}$,
$x_{1},\dots,x_{n}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\sigma^{2})$.
First note that $\frac{\lambda_{j}^{(k)}}{m}$ converges to $\lambda_{j}^{*}$
(the $j$th largest eigenvalue of kernel $k_{0}(\cdot,\cdot)$) under uniform
sampling, when minibatch size $m$ tends to $\infty$. Define
$\widetilde{\gamma}$ by substituting $\lambda_{j}^{(k)}$ with
$m\lambda_{j}^{*}$ ($\lambda_{j}^{*}$ is the $j$th eigenvalue of kernel $k$)
in the definition (12) of $\gamma^{(k)}$:

$\widetilde{\gamma}(\bm{\mathbf{\theta}})=\frac{1}{m}\sum_{j=1}^{m}\left(\theta_{1}\lambda_{j}^{*}m+\theta_{2}\right)^{-2},$
(13)

then $\widetilde{\gamma}(\bm{\mathbf{\theta}}^{(k)})$ is a reasonable
approximation for $\gamma^{(k)}$ when $m$ is large enough under uniform
sampling. If we increase the length scale $l$ under uniform sampling, then it
is equivalent to decreasing the distance between points, and thus can serve as
an approximation to nearby sampling.

The following lemma illustrates how the length scale $l$ influences
$\widetilde{\gamma}(\bm{\mathbf{\theta}})$:

###### Lemma 3

For any $l_{0}>0$, there exists $m_{0}>0$ depending on
$\theta_{1},\theta_{2},\sigma,l_{0}$ such that as long as $m>m_{0}$,
$\widetilde{\gamma}(\bm{\mathbf{\theta}})$ is an increasing function of $l\geq
l_{0}$.

Lemma 3 suggests that, for large enough minibatches, larger length scale leads
to faster convergence for $\theta^{(k)}_{2}$ to
$\theta^{*}_{2}=\sigma_{\epsilon}^{2}$. Therefore, our hypothesis that “nearby
sampling may improve the curvature” is plausible. More numerical support will
be provided in Section 6.

###  3.4 Extension to Summation of Multiple Kernels

So far we have assumed the kernel function $k_{0}(\cdot,\cdot)$ to be known.
However, sometimes there might be several choices of potential kernels and it
is desired to learn which kernel is the most appropriate from the data instead
of manually picking one kernel. One possible approach is to let the covariance
function $k(\cdot,\cdot)$ in (1) be a linear combination of all potential $M$
kernels:

$k(\cdot,\cdot)=\sum_{l=1}^{M}\sigma_{f,l}^{2}k_{l}(\cdot,\cdot),$ (14)

and then learn the signal variances $\sigma_{f,l}^{2}$ associated with the
$M>1$ kernels (Rasmussen, 2003). The kernel selection problem then translates
to how well we can learn the signal variance parameters $\sigma_{f,l}^{2}$,
$l=1,\dots M$.

Under this extension of the classical GP model, let the hyperparameter

$\bm{\mathbf{\theta}}^{*}=(\sigma_{f,1}^{2},\dots,\sigma_{f,M}^{2},\sigma_{\epsilon}^{2})^{\top}\in\mathbb{R}^{M+1}.$

We can still write out the log-likelihood loss $\ell(\bm{\mathbf{\theta}})$ as
in (2), while the only difference lies in the formulation of
$\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}})$:

$\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}})=\sum_{l=1}^{M}\theta_{l}\bm{\mathbf{K}}_{f,n}^{(l)}+\theta_{M+1}\bm{\mathbf{I}}_{n},$

where $\bm{\mathbf{K}}_{f,n}^{(l)}$ is the kernel matrix of
$k_{l}(\cdot,\cdot)$ evaluated at $\bm{\mathbf{X}}_{n}$, i.e.,
$(\bm{\mathbf{K}}_{f,n}^{(l)})_{ij}=k_{l}(\bm{\mathbf{x}}_{i},\bm{\mathbf{x}}_{j})$.
Then it is straightforward to extend Algorithm 1 to this setting for learning
$\bm{\mathbf{\theta}}^{*}$.

Similar to Theorem 3.1 to Theorem 3.4, under some additional conditions on the
kernel matrices, we also have convergence guarantees for $\sigma_{f,1}^{2}$,
$\sigma_{\epsilon}^{2}$, and the full gradient when $k_{l}(\cdot,\cdot),1\leq
l\leq M$ have exponential eigendecay; for $\sigma_{\epsilon}^{2}$ and the full
gradient when $k_{l}(\cdot,\cdot),1\leq l\leq M$ have polynomial eigendecay.
The detailed theoretical results for this setting are included in Appendix B.

##  4 Proof Sketch

In this section, we present the proof sketch for the first part of Theorem 3.1
and Theorem 3.2 (kernels with exponential eigendecay). The proof of the second
part in Theorem 3.1, Theorem 3.3 and Theorem 3.4 follows similar ideas
although requiring more careful analysis. With a bit abuse of notation, we
will omit
$g(\bm{\mathbf{\theta}}^{(k)};\bm{\mathbf{X}}_{\xi_{k+1}},\bm{\mathbf{y}}_{\xi_{k+1}})$
to $g(\bm{\mathbf{\theta}}^{(k)})$ and denote its conditional expectation
$\mathbb{E}(g(\bm{\mathbf{\theta}}^{(k)})|\bm{\mathbf{X}}_{\xi_{k+1}})$ by
$g^{*}(\bm{\mathbf{\theta}}^{(k)})$. Similarly we define
$\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})=\mathbb{E}(\nabla\ell(\bm{\mathbf{\theta}}^{(k)})|\bm{\mathbf{X}}_{n})$.

Due to the bias in the stochastic gradient, we take the followings steps
instead of directly drawing the connection between
$g(\bm{\mathbf{\theta}}^{(k)})$ and $\nabla\ell(\bm{\mathbf{\theta}}^{(k)})$:

  * •

For proving the first part of Theorem 3.1:

    * –

We first show that the conditional expectation
$g^{*}(\bm{\mathbf{\theta}}^{(k)})$ of the stochastic gradient has a property
similar to strong convexity, see Lemma 4.

    * –

We then prove that $g(\bm{\mathbf{\theta}})$ is close to its conditional
expectation $g^{*}(\bm{\mathbf{\theta}})$ uniformly over all possible
$\bm{\mathbf{\theta}}$, and thus $g(\bm{\mathbf{\theta}}^{(k)})$ is close to
$g^{*}(\bm{\mathbf{\theta}}^{(k)})$. Applying Lemma 5 to each minibatch leads
to the desired result.

These two steps lead to the $O(\frac{1}{K})$ optimization error rate for
$(\theta^{(k)}_{2}-\theta^{*}_{2})^{2}$, and a statistical error rate
depending on $m$, as shown in Theorem 3.1.

  * •

For proving Theorem 3.2:

    * –

Lemma 5 suggests that $\nabla\ell(\bm{\mathbf{\theta}}^{(k)})$ is close to
$\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})$

    * –

The eigendecay of kernel matrices (see Lemma 6) ensures that
$\|\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})\|_{2}$ is controlled by
$(\theta^{(k)}_{2}-\theta^{*}_{2})^{2}$, which is upper bounded in Theorem
3.1.

These steps above provide us with the same error bound of
$\|\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})\|_{2}$ from that of
$(\theta^{(k)}_{2}-\theta^{*}_{2})^{2}$ in Theorem 3.1.

###  4.1 Key Lemmas

The following two lemmas are the key building blocks of the proof: one shows
that $g^{*}(\bm{\mathbf{\theta}}^{(k)})$ enjoys a property similar to strong
convexity, the other establishes a uniform bound for the statistical error
$\nabla\ell(\bm{\mathbf{\theta}})-\nabla\ell^{*}(\bm{\mathbf{\theta}})$ over
$\bm{\mathbf{\theta}}\in[\theta_{\min},\theta_{\max}]^{2}$, and thus also
bounds for $g(\bm{\mathbf{\theta}}^{(k)})-g^{*}(\bm{\mathbf{\theta}}^{(k)})$;

######  Lemma 4 (Strongly convex-like property of
$g^{*}(\bm{\mathbf{\theta}}^{(k)})$)

Under Assumptions 3.1 to 3.3, if $s_{2}(m)=m$, $m>C$, then with probability at
least $1-3Km^{-c}$, the following claim holds true for $0\leq k<K$:

$(\theta^{(k)}_{2}-\theta^{*}_{2})(g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{2}\geq\frac{1}{8\theta_{\max}^{2}}(\theta^{(k)}_{2}-\theta^{*}_{2})^{2}-\frac{C\log
m}{m},$ (15)

Here $C>0$ depends only on $\theta_{\min},\theta_{\max},b$.

Lemma 4 is a relaxation of strong convexity, but leads to similar convergence
guarantees from running SGD on strongly convex objectives. The approximate
“curvature" parameter, $\frac{1}{8\theta_{\max}^{2}}$ on the R.H.S of (15),
remains a constant regardless of how large $m$ is. To guarantee the constant
“curvature”, we establish novel upper and lower bounds on
$\sum_{j=1}^{m}\lambda_{j}^{l}(\theta^{(k)}_{1}\lambda_{j}+\theta^{(k)}_{2})^{-2},l=0,1,2$
with high probability when $m$ is large, where $\lambda_{j}$ is the $j$th
largest eigenvalue of $\bm{\mathbf{K}}_{f,\xi}$ (see Lemma 6). The proof of
Lemma 6 is based on established error bounds for the empirical eigenvalues in
Braun (2006) and the eigendecay of the kernel $k_{0}(\cdot,\cdot)$.

######  Lemma 5 (Uniform statistical error)

Under Assumptions 3.1 to 3.2, Assumption 3.3 or 3.4, for any
$0<\varepsilon<\frac{1}{2}$, $1\leq i\leq 2$, we have

$\mathbb{P}\left(\sup_{\bm{\mathbf{\theta}}\in[\theta_{\min},\theta_{\max}]^{2}}\left|(\nabla\ell(\bm{\mathbf{\theta}}))_{i}-(\nabla\ell^{*}(\bm{\mathbf{\theta}}))_{i}\right|>Cn^{-\frac{1}{2}+\varepsilon}\right)\leq
C\exp\\{-cn^{2\varepsilon}\\}.$ (16)

Here $c,C>0$ only depend on $\theta_{\min}$, $\theta_{\max},b$.

The major difficulty in the proof of Lemma 5 is to control the error term
uniformly over $\bm{\mathbf{\theta}}\in[\theta_{\min},\theta_{\max}]^{2}$. We
need a uniform error bound, since $g^{*}(\bm{\mathbf{\theta}}^{(k)})$ is no
longer the conditional expectation of $g(\bm{\mathbf{\theta}}^{(k)})$ if
conditioning on the past iterate $\bm{\mathbf{\theta}}^{(k)}$. Although the
set $[\theta_{\min},\theta_{\max}]^{2}$ has constant dimension, the kernel
matrix $\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}})\in\mathbb{R}^{n\times n}$ is
of high dimension and is determined by $\bm{\mathbf{\theta}}$ in a non-linear
way. Our solution is to explore the Taylor’s expansion of
$\nabla\ell(\bm{\mathbf{\theta}})-\nabla\ell^{*}(\bm{\mathbf{\theta}})$, then
use truncation and covering arguments.

###  4.2 Proof of the First Part of Theorem 3.1

Let
$\widehat{e}_{k}=(g(\bm{\mathbf{\theta}}^{(k)}))_{2}-(g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{2}$.
Due to Lemma 4 and Assumption 3.2, we have

$\begin{split}\left(\theta^{(k)}_{2}-\theta^{*}_{2}\right)^{2}=&\left(\theta^{(k-1)}_{2}-\theta^{*}_{2}\right)^{2}-2\alpha_{k}(\theta^{(k-1)}_{2}-\theta^{*}_{2})(g(\bm{\mathbf{\theta}}^{(k-1)}))_{2}\\\
&+\alpha_{k}^{2}(g(\bm{\mathbf{\theta}}^{(k-1)}))_{2}^{2}\\\
\leq&\left(\theta^{(k-1)}_{2}-\theta^{*}_{2}\right)^{2}(1-\alpha_{k}\gamma)+\alpha_{k}^{2}G^{2}\\\
&+2\alpha_{k}\left(\frac{C\log
m}{m}-(\theta^{(k-1)}_{2}-\theta^{*}_{2})\widehat{e}_{k-1}\right),\end{split}$
(17)

where $\gamma=\frac{1}{4\theta_{\max}^{2}}$. Recall that
$\frac{3}{2\gamma}\leq\alpha_{1}\leq\frac{2}{\gamma}$, and
$\alpha_{k}=\frac{\alpha_{1}}{k}$ for all $k\geq 1$. Now we prove the
following statement for $k\geq 1$ by induction:

$\left(\theta^{(k)}_{2}-\theta^{*}_{2}\right)^{2}\leq\frac{2\alpha_{1}^{2}G^{2}}{k+1}+\sum_{i=0}^{k-1}\eta_{k,i}\left(\frac{C\log
m}{m}-(\theta^{(k-1)}_{2}-\theta^{*}_{2})\widehat{e}_{k-1}\right),$ (18)

where $\eta_{k,i}=2\alpha_{i+1}\prod_{j=i+2}^{k}(1-\alpha_{j}\gamma)$. When
$k=1$, by (17) and the fact that $1-\alpha_{1}\gamma<0$,

$\left(\theta^{(1)}_{2}-\theta^{*}_{2}\right)^{2}\leq\alpha_{1}^{2}G^{2}+\eta_{1,0}\left(\frac{C\log
m}{m}-(\theta^{(0)}_{2}-\theta^{*}_{2})\widehat{e}_{0}\right).$ (19)

Assuming (18) holds for $k=l\geq 1$, then due to (17) and the fact that
$1-\alpha_{l+1}\gamma\geq 0$ for $l\geq 1$, we have

$\begin{split}&\left(\theta^{(l+1)}_{2}-\theta^{*}_{2}\right)^{2}\\\
\leq&\left(\frac{2\alpha_{1}^{2}G^{2}}{l+1}+\sum_{i=0}^{l-1}\eta_{l,i}\left(\frac{C\log
m}{m}-(\theta^{(i)}_{2}-\theta^{*}_{2})\widehat{e}_{i}\right)\right)(1-\alpha_{l+1}\gamma)+\alpha_{l+1}^{2}G^{2}\\\
&+2\alpha_{l+1}\left(\frac{C\log
m}{m}-(\theta^{(l)}_{2}-\theta^{*}_{2})\widehat{e}_{l}\right)\\\
\leq&\frac{2\alpha_{1}^{2}G^{2}(l+1-\alpha_{1}\gamma)}{(l+1)^{2}}+\frac{\alpha_{1}^{2}G^{2}}{(l+1)^{2}}+\sum_{i=0}^{l}\eta_{l+1,i}\left(\frac{C\log
m}{m}-(\theta^{(i)}_{2}-\theta^{*}_{2})\widehat{e}_{i}\right)\\\
\leq&\frac{2\alpha_{1}^{2}G^{2}}{l+2}+\sum_{i=0}^{l}\eta_{l+1,i}\left(\frac{C\log
m}{m}-(\theta^{(i)}_{2}-\theta^{*}_{2})\widehat{e}_{i}\right).\end{split}$
(20)

Here the last two lines are due to range of $\alpha_{1}$ and the definitions
of $\eta_{l,i}$. The next step is to bound
$\sum_{i=0}^{K-1}\eta_{K,i}\left(\frac{C\log
m}{m}-(\theta^{(i)}_{2}-\theta^{*}_{2})\widehat{e}_{i}\right)$. First we have

$\begin{split}&\left|\sum_{i=0}^{K-1}\eta_{K,i}\left(\frac{C\log
m}{m}-(\theta^{(i)}_{2}-\theta^{*}_{2})\widehat{e}_{i}\right)\right|\\\
\leq&\frac{2\alpha_{1}}{K}\sum_{i=0}^{K-1}\left|\theta^{(i)}_{2}-\theta^{*}_{2}\right|\widehat{e}_{i}|+\frac{C\alpha_{1}\log
m}{m}\\\ \leq&C\left(\max_{0\leq i\leq K-1}|\widehat{e}_{i}|+\frac{\log
m}{m}\right).\end{split}$ (21)

Note that the distribution of each minibatch
$\\{\bm{\mathbf{X}}_{\xi_{k+1}},\bm{\mathbf{y}}_{\xi_{k+1}}\\}_{i=1}^{m}$ is
the same as sampling $m$ independent $\bm{\mathbf{x}}_{i}$ from $\mathbb{P}$,
and then sampling
$\bm{\mathbf{y}}_{\xi_{k+1}}\sim\mathcal{N}(0,\bm{\mathbf{K}}^{*}_{\xi_{k+1}})$,
thus we can apply the Lemma 5 on each
$\widehat{e}_{i}=g(\bm{\mathbf{\theta}}^{(i)})_{2}-g^{*}(\bm{\mathbf{\theta}}^{(i)})_{2}$
and take a union bound over $0\leq i\leq K-1$:

$\mathbb{P}\left(\max_{0\leq i\leq
K-1}\left|\widehat{e}_{i}\right|>Cm^{-\frac{1}{2}+\varepsilon}\right)\leq
CK\exp\\{-cm^{2\varepsilon}\\},$ (22)

for any $\varepsilon>0$. Therefore,

$\left(\theta^{(k)}_{2}-\theta^{*}_{2}\right)^{2}\leq\frac{2\alpha_{1}^{2}G^{2}}{k+1}+Cm^{-\frac{1}{2}+\varepsilon},$
(23)

with probability at least

$1-CK\exp\\{-c\min\\{\log m,m^{2\varepsilon}\\}\geq
1-CK\exp\\{-cm^{2\varepsilon}\\},$ (24)

for any $0<\varepsilon<C\frac{\log\log m}{\log m}$, when $m>C$ for some $C>0$
depending on $\theta_{\min},\theta_{\max},b$.

###  4.3 Proof of Theorem 3.2

We start from bounding $\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})$, the
conditional expectation of $\nabla\ell(\bm{\mathbf{\theta}}^{(k)})$ given
$\bm{\mathbf{x}}_{1},\dots,\bm{\mathbf{x}}_{n}$, then control the statistical
error
$\nabla\ell(\bm{\mathbf{\theta}}^{(k)})-\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})$.
By the definition of $\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})$, for $1\leq
i\leq 2$,

$\begin{split}\left(\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})\right)_{i}=&\frac{1}{2n}\mathrm{tr}\left[\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}}^{(k)})^{-1}(\bm{\mathbf{I}}_{n}-\bm{\mathbf{K}}_{n}^{*}\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}}^{(k)})^{-1})\frac{\partial\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}}^{(k)})}{\partial\theta^{(k)}_{i}}\right]\\\
=&\frac{1}{2n}\sum_{j=1}^{n}\frac{(\theta^{(k)}_{1}-\theta^{*}_{1})\lambda_{j}^{1+\mathbbm{1}_{\\{i=1\\}}}+(\theta^{(k)}_{2}-\theta^{*}_{2})\lambda_{j}^{\mathbbm{1}_{\\{i=1\\}}}}{\left(\theta^{(k)}_{1}\lambda_{j}+\theta^{(k)}_{2}\right)^{2}},\end{split}$
(25)

where $\lambda_{j}$ is the $j$th largest eigenvalue of
$\bm{\mathbf{K}}_{f,n}$. The following lemma provides bounds for
$\sum_{j=1}^{n}\frac{\lambda_{j}^{l}}{\left(\theta^{(k)}_{1}\lambda_{j}+\theta^{(k)}_{2}\right)^{2}}$
for all $l=0,1,2$.

###### Lemma 6

Under Assumption 3.3, for any $\alpha>0$, if $n>C$ for $C>0$ depending on $b$,
then with probability at least $1-3n^{-\alpha}$,

$\begin{split}\text{ if }l=1\text{ or
}2,\sum_{j=1}^{n}\frac{\lambda_{j}^{l}}{\left(\theta_{1}\lambda_{j}+\theta_{2}\right)^{2}}\leq&\frac{2(2+\alpha)}{b\theta_{\min}^{2}}\log
n,\\\
\sum_{j=1}^{n}\frac{1}{\left(\theta_{1}\lambda_{j}+\theta_{2}\right)^{2}}\leq&\frac{n}{\theta_{\min}^{2}},\end{split}$
(26)

holds for any $\bm{\mathbf{\theta}}\in[\theta_{\min},\theta_{\max}]^{2}$.

We prove Lemma 6 by exploiting the error bounds for eigenvalues of empirical
kernel matrices from the population eigenvalues of the kernel operator. A
detailed version of Lemma 6 including also lower bounds for
$\frac{\lambda_{j}^{l}}{\left(\theta_{1}\lambda_{j}+\theta_{2}\right)^{2}}$ is
presented in the Appendix, which is a key result for proving Lemma 4.

For any constant $c>0$, apply Lemma 6 with $\alpha=c$, then (26) holds with
probability at least $1-3n^{-c}$, if $n>C$ for $C$ depending on $b$. Combining
this result and (25) together implies that

$\begin{split}\left|\left(\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})\right)_{1}\right|\leq&\frac{C\log
n}{n}\end{split}$ (27)

where $C>0$ depends on $\theta_{\min},\theta_{\max},b$. Meanwhile,

$\begin{split}\left|\left(\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})\right)_{2}\right|\leq
C\left(|\theta^{(k)}_{2}-\theta^{*}_{2}|+\frac{\log n}{n}\right).\end{split}$
(28)

Thus we have

$\begin{split}\|\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})\|_{2}^{2}\leq&C\left[\left(\frac{\log
n}{n}\right)^{2}+(\theta^{(k)}_{2}-\theta^{*}_{2})^{2}\right].\end{split}$
(29)

For bounding
$\nabla\ell(\bm{\mathbf{\theta}}^{(k)})-\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})$,
we can apply Lemma 5. By (16) and Theorem 3.1, for any
$0<\varepsilon<C\frac{\log\log m}{\log m}$, if $m>C$, then with probability at
least $1-CK\exp\\{-cm^{2\varepsilon}\\}$, we have

$\|\nabla\ell(\bm{\mathbf{\theta}}^{(K)})\|_{2}^{2}\leq
C\left[\frac{G^{2}}{K+1}+m^{-\frac{1}{2}+\varepsilon}\right],$ (30)

where $c,C>0$ depend only on $\theta_{\min},\theta_{\max},b$.

##  5 Practical Considerations for Applying SGD on GP

###  5.1 Sampling Scheme

As discussed in Section 2, one may consider both uniformly and nearby
sampling. In Section 3.3, we provide both numerical and theoretical evidence
that hint to the fact that sampling nearby points for each minibatch leads to
faster convergence specifically for the noise variance parameter
$\theta_{2}^{(k)}$ towards $\theta_{2}^{*}=\sigma_{\epsilon}^{2}$. Besides
that, our case studies in numerical experiments in Section 6 also support this
claim.

Below we highlight the Algorithm for nearby sampling which is a simple
extension of Algorithm 1

1 Input: $\bm{\mathbf{\theta}}^{(0)}\in\mathbb{R}^{2}$, initial step size
$\alpha_{1}>0$.

2 for  _$k=1,2,\dotsc,K$_ do

3 Sample a data point uniformly from the data pool, and then select its $m-1$
nearest neighbors, which forms $(\mathbf{X}_{\xi_{k}},\mathbf{y}_{\xi_{k}})$
of size $m$;

4 Compute the stochastic gradient
$g(\bm{\mathbf{\theta}}^{(k-1)};\mathbf{X}_{\xi^{(2)}_{k}},\mathbf{y}_{\xi^{(2)}_{k}})$;

5 $\alpha_{k}\leftarrow\frac{\alpha_{1}}{k}$;

6
$\bm{\mathbf{\theta}}^{(k)}\leftarrow\bm{\mathbf{\theta}}^{(k-1)}-\alpha_{k}g(\bm{\mathbf{\theta}}^{(k-1)};\mathbf{X}_{\xi_{k}},\mathbf{y}_{\xi_{k}})$;

7

8 end for

Algorithm 2 Minibatch SGD with nearby sampling

###  5.2 Optimizing other hyperparameters

In practice, we may also need to determine other hyperparameters of the kernel
function besides signal variance and noise variance. For example, when
considering the RBF kernel
$k(\mathbf{x},\mathbf{x}^{\prime})=\exp\\{-\sum_{j=1}^{d}\frac{(x_{j}-x^{\prime}_{j})^{2}}{2l_{j}^{2}}\\}$,
we need to estimate the lengthscale parameters $l_{j},1\leq j\leq d$; when
considering the Matèrn kernel (9), the hyperparameters $\alpha$ and $h$ are
also unknown and require estimation. Similar to Algorithm 1 and Algorithm 2,
we can update these parameters alongside the variance parameters using
minibatch SGD. Our numerical experiments in Section 6 suggest that nearby
sampling method may also be a good option for the lengthscale parameters.

###  5.3 Prediction

Although our main focus in this paper is estimating the hyperparameters, or
model selection, the last step when applying GP in real applications is always
prediction. Following the model selection process from which we obtain optimal
hyperparameters, various strategies can be applied to calculate the predictive
mean for $\mathbf{x}_{*}$ and the predictive covariance between
$\mathbf{x}_{*}$ and $\mathbf{x}_{*}^{\prime}$ using the well known predictive
equation below

$\mu_{\text{pred}}(\mathbf{x}_{*})=\mathbf{k}_{\mathbf{X}_{n}\mathbf{x}_{*}}^{\top}\mathbf{K}_{n}^{-1}\mathbf{y}_{n},\quad\quad
k_{\text{pred}}(\mathbf{x}_{*},\mathbf{x}_{*}^{\prime})=k(\mathbf{x}_{*},\mathbf{x}_{*}^{\prime})-\mathbf{k}_{\mathbf{X}_{n}\mathbf{x}_{*}}^{\top}\mathbf{K}_{n}^{-1}\mathbf{k}_{\mathbf{X}_{n}\mathbf{x}_{*}^{\prime}}\,,$
(31)

where
$\mathbf{k}_{\mathbf{X}_{n}\mathbf{x}_{*}}=(k(\mathbf{x}_{1},\mathbf{x}_{*})\dotsc,k(\mathbf{x}_{n},\mathbf{x}_{*}))^{\top}$.
The main computational cost of the predictive mean and the predictive
covariance come from $\mathbf{K}_{n}^{-1}$. In general, for $n<10^{4}$, they
can be computed via Cholesky decomposition; for $n<10^{5}$, preconditioned
conjugate gradient (PCG) (Gardner et al., 2018) can be applied for
acceleration; for $n<10^{6}$, PCG with partitioned kernel (Wang et al., 2019)
could provide further speed up, if distributed computational resources are
available. Another practical but less ideal strategy when predicting with
extremely large $n$ is to follow the same approach as nearby sampling and
utilize only $\tilde{n}$ nearest neighbors of $\mathbf{x}_{*}$ within the
observed data to solve (31), where $\tilde{n}<n$ is determined by the
available computational resource. Fortunately, prediction is a one-shot
process compared to the iterative training process.

##  6 Numerical Results

###  6.1 Numerical Illustration of Theory

In this section, we conduct simulation studies to verify our theoretical
results.

We consider $n=1,024$,
$\mathbf{x}_{i}\overset{i.i.d.}{\sim}\mathcal{N}(0,5^{2})$ and
$\mathbf{y}_{n}\sim\mathcal{N}(\mathbf{0},\sigma_{f}^{2}\mathbf{K}_{f,n}+\sigma_{\epsilon}^{2}\mathbf{I}_{n})$,
where $\mathbf{K}_{f,n}$ is an RBF kernel matrix with known lengthscale
$l=0.5$. The underlying true parameters are outputscale $\sigma_{f}^{2}=4$ and
noise variance $\sigma_{\epsilon}^{2}=1$. In each experiment, we perform 25
epochs of minibatch SGD with diminishing step sizes $\alpha_{k}=\alpha_{1}/k$.
We set scaling factors to $s_{1}(m)=3\log m$ for $\sigma_{f}^{2}$ and
$s_{2}(m)=m$ for $\sigma_{n}^{2}$. Each experiment is repeated 10 times with
independent data pools for different repetitions.

Figure 2: Illustration of the convergence of parameters under uniform
sampling. We consider $m=128$, and demonstrate three cases with varying
initial points where initial stepsizes are $\alpha_{1}=$ $9,9$ and $6$,
respectively. Lines in black denote the true parameters.

Fig. 2 shows the convergence of parameters under uniform sampling, varying
initializations and step sizes. All, the curves display $O(\frac{1}{K})$
convergence rates which are consistent with our results in Theorem 3.1.
Moreover, the locations where the updates of $\sigma_{\epsilon}^{2}$ converges
to, are significantly more concentrated around the truth compared to that of
$\sigma_{f}^{2}$, which is consistent with the $O((\log m)^{-\frac{1}{2}})$
statistical error for $\sigma_{f}^{2}$ and $O(m^{-\frac{1}{2}})$ statistical
error for $\sigma_{\epsilon}^{2}$, also stated in Theorem 3.1.

Figure 3: Comparison of the convergence of the full gradient under uniform
sampling with varying minibatch sizes. The mean of
$\|\nabla\ell(\bm{\mathbf{\theta}}^{(k)})\|_{2}^{2}$ is shown in blue and the
region within its one standard error over 10 repetitions is shown in light
blue, both under log scale. The three experiments share initial point
$\bm{\mathbf{\theta}}^{(0)}=(5.0,3.0)$ and inital step size $\alpha_{1}=9$.

Fig. 3 displays the effect of minibatch size on the convergence of the full
gradient. To start with, the curves flatten slower and converge to larger
values as minibatch size decreases, suggesting that a larger minibatch leads
to faster convergence of the full gradient, as well as a full gradient with
smaller statistical error. In addition, the convergence points of
$\log(\|\nabla\ell(\bm{\mathbf{\theta}}^{(k)})\|_{2}^{2})$ scale linearly with
minibatch size $m$, indicating a $O(m^{-\frac{1}{2}})$ statistical error for
$\|\nabla\ell(\bm{\mathbf{\theta}}^{(k)})\|_{2}^{2}$. The above observations
confirm our statements in Theorem 3.2.

Figure 4: Comparison of the convergence of parameters with varying minibatch
sizes. Lines in black denote the true parameters. The three experiments share
initial point $\bm{\mathbf{\theta}}^{(0)}=(5.0,3.0)$ and inital step size
$\alpha_{1}=9$.

We also investigate how minibatch size $m$ influences the convergence of
parameters, which is illustrated in Fig. 4. As highlighted in our theory, we
find that a larger mini-batch size results in faster convergence and smaller
statistical error (more concentrated curves) for the parameters. Here we note
that similar results hold for the Matérn kernel and hence we omit the figures.

###  6.2 Case Studies

In this section, we test our model’s performance on publicly available real
and simulated datasets. Our benchmarked models are: (i) Exact inference using
matrix vector multiplication denoted as EGP (Gardner et al., 2018; Wang et
al., 2019), (ii) Vecchia’s GP approximation denoted as Vecchia (Guinness,
2018; Katzfuss et al., 2020), (iii) sparse GP regression denoted as SGPR
(Titsias, 2009) and (iv) stochastic variational GP denoted as SVGP (Hensman et
al., 2013). Our stochastic gradient-based GP approach is denoted as sgGP.

All models are tested on real datasets from the UCI repository (Dua and Graff,
2017) and simulated datasets from the Virtual Library of Simulation
Experiments (Surjanovic and Bingham, ).The real datasets are: Bike, Energy,
PM2.5, Protein and Query. The simulated dataset are Levy, Greiwank and
Borehole. We also use two other simulated datasets from the Virtual Library
Simulation Experiments that represent real-life systems. The OTL circuit
models an output transformerless push-pull circuit while the Wing Weight
models a light aircraft wing.

Throughout all experiments, we consider constant zero prior mean function and
the scaled RBF covariance function with a separate lengthscale for each input
dimension. We run Adam to learn the signal variance, noise variance, and also
lengthscales, as an extension from our problem set-up stated in Section 2. We
conduct 10 independent trials on each dataset. In each trial, we randomly
split the dataset into 60% training set and 40% test set. In addition, the
training set is normalized to 0 mean and 1 standard deviation, and the test
set is scaled accordingly.

During model selection, the hyperparameters and variational parameters are
learned through minimizing the negative log marginal likelihood or its
surrogate. (i) For sgGP, we consider both uniform and nearest neighbor
sampling schemes, where we perform 100 epochs of Adam with minibatch size
$m=16$ and a learning rate of 0.01. (ii) For Vecchia, we order the observed
data following the maximum minimum distance (MMD) ordering (Guinness, 2018).
MMD ordering works by first selecting a center point, and then sequentially
selecting the next point to have maximum minimum distance to all previously
selected points. We let each observed response condition on its $m=16$ nearest
neighbors within its predecessors from the ordered set. We also carry out 100
iterations of the Fisher scoring (Guinness, 2021) algorithm. (iii) For other
methods, we follow the theoretical recommendations in Burt et al. (2019) and
the practical recommendations in Wang et al. (2019). Further, for EGP, we
perform 100 iterations of Adam with a learning rate of 0.1. For SGPR, we use
$m=512$ inducing points and carry out 100 iterations of Adam with a learning
rate of 0.1. For SVGP, we use $m=1,024$ inducing points and perform 100 epochs
of Adam with a minibatch size of $1,024$ and a learning rate of 0.01. To
ensure fairness of comparison, we do not perform any pretraining or fine-
tuning, and we let different methods share a common but randomly selected
starting point in each trial.

Regarding the prediction of sgGP, we adopt the PCG algorithm in EGP to
approximate (31). While for prediction in Vecchia, we order the inputs to be
predicted using MMD ordering and append them to the ordered observed inputs.
We set the conditioning-set size to $m=64$.

Model selection of sgGP is coded with R, where RANN package (Arya et al.,
2019) is used for finding nearest neighbors. Vecchia is coded using R, where
we utilize GpGp package (Guinness, 2018) to find ordered nearest neighbors and
implement Fisher’s scoring algorithm. The prediction of sgGP, together with
EGP, SGPR and SVGP are implemented through GPyTorch (Gardner et al., 2018).
Each experiment is performed on a single core of Intel Xeon E5-2680 v3 @
2.50GHz CPU. For simulated datasets, we manually inject noise to the response.
For query dataset, we constrain the learned noise to be at least 0.1 to
regularize the ill-conditioned kernel matrix. Due to memory limit, for
Borehole, OTL Circuit and Wing Weight datasets, we use PCG algorithm for
prediction in sgGP yet using only 60,000 nearest neighbors of each test point.

The results of our experiments are shown in Tables 3 \- 4. We start first by
analyzing Tables 3 and 3. Table 3 summarizes the prediction accuracy of all
benchmarked methods while Table 3 highlights the accuracy of the learned noise
variance $\hat{\sigma}_{n}^{2}$ on simulated datasets where we know the
underlying truth. Based on the results, one can derive many insights.

First, we find that sgGP equipped with nearest-neighbor sampling (sgGP (nn))
exhibits the best predictive performance among the various methods on datasets
with varying sizes, input dimensions, and noise levels. In addition, its
learned noise variance is significantly more accurate than all benchmarks.
Second, while sgGP (uni) can sometimes achieve good performance, it performs
poorly in comparison to sgGP (nn), Vecchia, and EGP. This supports our
numerical and theoretical evidence of the advantages of nearby sampling in
Section 3.3. Third, Vecchia does perform well overall in terms of prediction
performance, which is contrary to the finding in Jankowiak and Pleiss (2021).
Most likely, the heuristic MMD ordering we adopted offers significant
improvement in model approximation over the default coordinate-based ordering
(Guinness, 2018). However, Vecchia significantly underestimates the noise
level and subsequently yields lower prediction accuracy than sgGP (nn). Here,
it should be noted that the ordering of observations is crucial for the
quality of Vecchia’s approximation, and therefore, extensive effort towards
dataset-specific tuning may be required, yet there lacks heuristic guidance
and theoretical support for datasets of higher dimensions. Fourth, EGP
exhibits inferior prediction accuracy compared to sgGP (nn). This highlights
the ability of sgGP (nn) to learn parameters that generalize better as both
sgGP and EGP aim at exact inference. Yet, it should be noted that while EGP
tackles exact inference, it features many approximations within. Finally, we
find that SGPR and SVGP both do poorly overall and yield twice the prediction
errors of sgGP (nn) on datasets like Levy, PM 2.5 and Query. Also, SGPR and
SVGP (especially) tend to exaggerate the noise level (Bauer et al., 2016;
Jankowiak et al., 2020), as seen in Table 3. Similar to Wang et al. (2019),
this finding sheds light on the ability of exact GPs to significantly benefit
from the increase in the number of training points.

Table 3 summarizes the training time of all competing methods. The results
exhibit the overwhelming time advantage of sgGP in training, which
significantly scales with dataset size. Not only does sgGP achieve better
generalization, but it also does that in a fraction of the training time
needed for competing methods. This result is again confirmed by our test of
the application-driven simulated datasets of size $2\times 10^{6}$ in Table 4.
Remarkably, it takes around 30 minutes to perform model selection for OTL
Circuit dataset using a single core with R functions that are not designed for
fast execution. In addition, sgGP enjoys superior memory efficiency due to the
use of minibatches. These experiments justify that SGD does open up a new data
size regime for exploring GPs. Here we note that we are aware that EGP is
designed to leverage multiple GPU parallelization; however, much like regular
SGD, sgGP can be readily extended to a batch version where the gradient
estimate in each update is the average of $M$ gradient estimates from $M$
mini-batches. This allows sgGP to take advantage of parallel computing when
the hardware is available.

Table 1: We summarize the RMSE of sgGP and other GPs on benchmark datasets.
Here and elsewhere, we report the averages $\pm$ standard errors over $10$
dataset splits. Best results are in bold (lower is better). sgGP (unif)
utilizes uniform minibatches for training and sgGP (nn) utilizes nearest
neighbor minibatches for training. For query and borehole datasets, we are
unable to train with EGP due to memory limit.

RMSE Dataset Size $D$ sgGP (uni) sgGP (nn) Vecchia EGP SGPR SVGP Levy 10,000 4
$0.593\pm 0.002$ $\mathbf{0.264\pm 0.002}$ $0.309\pm 0.003$ $0.316\pm 0.004$
$0.564\pm 0.010$ $0.582\pm 0.013$ Griewank 10,000 6 $0.149\pm 0.003$ $0.070\pm
0.000$ $0.081\pm 0.007$ $\mathbf{0.064\pm 0.000}$ $0.132\pm 0.003$ $0.093\pm
0.005$ Bike 17,379 17 $0.227\pm 0.002$ $\mathbf{0.220\pm 0.002}$ $0.223\pm
0.002$ $0.228\pm 0.002$ $0.276\pm 0.004$ $0.250\pm 0.010$ Energy 19,735 27
$\mathbf{0.712\pm 0.005}$ $0.786\pm 0.008$ $0.738\pm 0.006$ $0.802\pm 0.007$
$0.843\pm 0.006$ $0.795\pm 0.005$ PM2.5 41,757 15 $0.573\pm 0.003$
$\mathbf{0.286\pm 0.002}$ $0.385\pm 0.004$ $0.287\pm 0.003$ $0.638\pm 0.005$
$0.540\pm 0.010$ Protein 45,730 9 $0.829\pm 0.002$ $0.659\pm 0.004$
$\mathbf{0.597\pm 0.001}$ $0.696\pm 0.004$ $0.715\pm 0.003$ $0.676\pm 0.004$
Query 100,000 4 $0.128\pm 0.000$ $0.027\pm 0.000$ $\mathbf{0.024\pm 0.000}$ ––
$0.058\pm 0.002$ $0.061\pm 0.000$ Borehole 1,000,000 8 $0.173\pm 0.000$
$\mathbf{0.172\pm 0.000}$ $0.174\pm 0.000$ –– $0.176\pm 0.000$ $0.173\pm
0.000$

Table 2: We summarize the accuracy of learned noise variance of sgGP and other
GPs on simulated datasets. Results follow the experiments in Table 3.

Learned Noise Variance/True Noise Variance
($\bm{\hat{\sigma}_{n}^{2}}/\bm{\sigma_{n}^{2}}$) Dataset Size $D$ sgGP (uni)
sgGP (nn) Vecchia EGP SGPR SVGP Levy 10,000 4 $9.95\pm 0.16$ $1.32\pm 0.05$
$0.12\pm 0.01$ $0.14\pm 0.01$ $2.32\pm 0.12$ $11.05\pm 0.45$ Griewank 10,000 6
$0.03\pm 0.00$ $1.38\pm 0.03$ $0.02\pm 0.01$ $0.48\pm 0.02$ $9.18\pm 0.49$
$22.77\pm 6.74$ Borehole 1,000,000 8 $0.96\pm 0.02$ $0.99\pm 0.02$ $0.00\pm
0.00$ –– $0.86\pm 0.08$ $1.97\pm 0.21$

Table 3: We summarize the training time of sgGP and other GPs on benchmark
dataset. Results follow the experiments in Table 3.

Training Time (min) Dataset Size $D$ sgGP (uni) sgGP (nn) Vecchia EGP SGPR
SVGP Levy 10,000 4 $0.41\pm 0.02$ $0.35\pm 0.03$ $3.15\pm 0.12$ $11.12\pm
0.71$ $3.55\pm 0.22$ $14.74\pm 0.69$ Griewank 10,000 6 $0.54\pm 0.02$ $0.48\pm
0.03$ $4.24\pm 0.09$ $13.37\pm 1.18$ $1.76\pm 0.12$ $14.60\pm 0.65$ Bike
17,379 17 $1.50\pm 0.12$ $1.55\pm 0.10$ $6.70\pm 0.47$ $29.48\pm 3.96$
$5.31\pm 2.05$ $25.26\pm 3.97$ Energy 19,735 27 $3.03\pm 0.02$ $2.58\pm 0.17$
$10.44\pm 1.29$ $53.25\pm 2.47$ $5.41\pm 0.73$ $25.09\pm 5.50$ PM2.5 41,757 15
$4.90\pm 0.33$ $4.24\pm 0.29$ $13.69\pm 0.70$ $372.88\pm 16.78$ $13.59\pm
2.30$ $52.46\pm 10.08$ Protein 45,730 9 $3.12\pm 0.01$ $2.63\pm 0.17$
$13.06\pm 0.12$ $453.40\pm 21.31$ $19.55\pm 1.66$ $55.27\pm 13.09$ Query
100,000 4 $4.73\pm 0.36$ $5.03\pm 0.36$ $30.86\pm 1.69$ –– $20.73\pm 1.63$
$124.73\pm 22.25$ Borehole 1,000,000 8 $54.82\pm 2.69$ $65.19\pm 2.98$
$235.74\pm 18.00$ –– $857.60\pm 76.02$ $1380.86\pm 11.32$

Table 4: We summarize the results of sgGP (nn) on simulated application-driven
datasets. We follow similar setups of the experiments in Table 3.

Dataset Size $D$ RMSE Training Time (min) Memory Usage (GB) OTL Circuit
2,000,000 6 $0.401\pm 0.000$ $33.43\pm 4.40$ $0.99\pm 0.00$ Wing Weight
2,000,000 10 $0.072\pm 0.004$ $78.78\pm 9.26$ $1.22\pm 0.00$

##  7 Open Problems

There still exist some open problems that are worth future investigations.

  1. 1.

The extension to convergence guarantees for learning the lengthscale parameter
in RBF kernel is an interesting but extremely challenging problem: our case
studies suggest that SGD may still be used for estimating the lengthscale in
practice, but the proof for both Lemma 4 and Lemma 5 presents additional
challenges if looking at the lengthscale. This is due to it being wrapped
within the exponential term as a denominator, which translates to different
eigenvectors for $\frac{\partial\bm{\mathbf{K}}_{\xi}}{\partial l}$,
$\bm{\mathbf{K}}_{\xi}$ and $\bm{\mathbf{K}}_{\xi}^{*}$. To see the difficulty
for proving Lemma 4, note that the curvature term for estimating $l$ involves

$\mathrm{tr}\left(\bm{\mathbf{K}}_{\xi}(\bm{\mathbf{\theta}}^{*})\bm{\mathbf{K}}_{\xi}(\bm{\mathbf{\theta}})^{-1}\frac{\partial\bm{\mathbf{K}}_{\xi}(\bm{\mathbf{\theta}})}{\partial
l}\bm{\mathbf{K}}_{\xi}(\bm{\mathbf{\theta}})^{-1}\frac{\partial\bm{\mathbf{K}}_{\xi}(\bm{\mathbf{\theta}})}{\partial
l}\bm{\mathbf{K}}_{\xi}(\bm{\mathbf{\theta}})^{-1}\right)$ (32)

and cannot be expressed as a function of the eigenvalues of kernel matrices
due to their different eigenvectors. It is also very hard to upper bound the
statistical error in Lemma 5 due to similar reasons.

  2. 2.

It would also be interesting to extend theoretical guarantees from $M=1$ to
$M>1$ without Assumption B.2. The technical challenge for this part is similar
to the previous point: without Assumption B.2, $K_{n}$ and $K_{n}^{*}$ can not
be simultaneously diagonalized, which hinders the proofs of Lemma 4 and Lemma
5.

  3. 3.

Another open problem is to establish convergence guarantees for running SGD
with nearby sampling and to explore different techniques for nearby sampling
upon the choice of the kernel.

##  8 Conclusion

In this paper, we provide theoretical guarantees for the minibatch SGD for the
model selection of Gaussian process (GP). In particular, we prove that the
iterates of SGD converge to the true hyperparameters and the critical point of
the full loss function, with rate $O(\frac{1}{K})$ up to a statistical error
term depending on minibatch size. Given the correlation structure of GPs, the
challenge lies in the bias of stochastic gradient when taking expectation
w.r.t. random sampling. Numerical studies support our theoretical results and
show that minibatch SGD has better performance than state-of-the-art methods
on various datasets while enjoying huge computational benefits.

##  A Table of Notations

Notations | Description
---|---
$n$ | number of data points in the full data set
$m$ | number of data points in a minibatch
$K$ | number of iterations of minibatch SGD
$G$ | An upper bound for $\|g(\bm{\mathbf{\theta}}^{(k)})\|_{2}$, specified in Assumption 3.2
$\sigma_{f}^{2}=\theta^{*}_{1}$ | true signal variance parameter
$\sigma_{\epsilon}^{2}=\theta^{*}_{2}$ | true noise variance parameter
$\bm{\mathbf{\theta}}^{(k)}$ | output of minibatch SGD at the $k$th iteration, as an estimate of $\bm{\mathbf{\theta}}^{*}$
$\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}})=\theta_{1}\bm{\mathbf{K}}_{f,n}+\theta_{2}\bm{\mathbf{I}}_{n}$ | covariance of $\bm{\mathbf{y}}_{n}$ given $\bm{\mathbf{X}}_{n}$, if the hyperparameter is $\bm{\mathbf{\theta}}$
$\bm{\mathbf{K}}_{\xi}(\bm{\mathbf{\theta}})$ | submatrix of $\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}})$ with rows and columns both indexed by $\xi$
$\bm{\mathbf{K}}_{f,n}$ | kernel matrix evaluated at $\bm{\mathbf{X}}_{n}$
$\nabla\ell(\bm{\mathbf{\theta}};\bm{\mathbf{X}}_{n},\bm{\mathbf{y}}_{n})$ or $\nabla\ell(\bm{\mathbf{\theta}})$ | full gradient evaluated at $\bm{\mathbf{\theta}}$ and full data $\bm{\mathbf{X}}_{n}$, $\bm{\mathbf{y}}_{n}$
$g(\bm{\mathbf{\theta}};\bm{\mathbf{X}}_{\xi},\bm{\mathbf{y}}_{\xi})$ | stochastic gradient evaluated at $\bm{\mathbf{\theta}}$ and minibatch $\bm{\mathbf{X}}_{\xi}$, $\bm{\mathbf{y}}_{\xi}$
$\alpha_{k}=\frac{\alpha_{1}}{k}$ | step size at the $k$th iteration
$g^{*}(\bm{\mathbf{\theta}}^{(k)};\bm{\mathbf{X}}_{\xi_{k+1}})$ or $g^{*}(\bm{\mathbf{\theta}}^{(k)})$ | conditional expectation of $g(\bm{\mathbf{\theta}}^{(k)};\bm{\mathbf{X}}_{\xi_{k+1}},\bm{\mathbf{y}}_{\xi_{k+1}})$
at the $k$th iteration given $\bm{\mathbf{X}}_{\xi_{k+1}}$
$\lambda_{j}^{(k)}$ | the $j$th largest eigenvalue of $\bm{\mathbf{K}}_{f,\xi_{k+1}}$
$\lambda_{j}^{*}$ | the $j$th largest eigenvalue of $\bm{\mathbf{K}}_{f,n}$
Table 5: Important notations used throughout the paper

##  B Theoretical Guarantees for Section 3.4

Before presenting the theoretical guarantees under this setting, we first
provide a formal definition for the considered minibatch SGD algorithm. With
sampled indices $\xi\subset[n]$, let the stochastic gradient
$g(\bm{\mathbf{\theta}};\mathbf{X}_{\xi},\mathbf{y}_{\xi})\in\mathbb{R}^{M+1}$
be defined as follows:

$\left(g(\bm{\mathbf{\theta}};\mathbf{X}_{\xi},\mathbf{y}_{\xi})\right)_{l}=\frac{1}{2s_{l}(m)}\text{tr}\left[(\mathbf{K}_{\xi}^{-1}(\mathbf{I}_{m}-\mathbf{y}_{\xi}\mathbf{y}_{\xi}^{\top}\mathbf{K}_{\xi}^{-1})\frac{\partial\mathbf{K}_{\xi}}{\partial\theta_{l}}\right],\quad
1\leq l\leq M+1,$ (33)

where $\mathbf{K}_{\xi}$ is the principle submatrix formed by the rows and
columns of $\mathbf{K}_{n}$ indexed by $\xi$. In the following we will also
let $\bm{\mathbf{K}}_{f,\xi}^{(l)}$ denote the $m\times m$ block of
$\bm{\mathbf{K}}_{f,n}^{(l)}$ indexed by $\xi$. Algorithm 3 summarizes the
steps of minibatch SGD.

1 Input: $\bm{\mathbf{\theta}}^{(0)}\in\mathbb{R}^{M+1}$, initial step size
$\alpha_{1}>0$.

2 for  _$k=1,2,\dotsc,K$_ do

3 Randomly sample a subset of indices $\xi_{k}$ of size $m$;

4 Compute the stochastic gradient
$g(\bm{\mathbf{\theta}}^{(k-1)};\mathbf{X}_{\xi_{k}},\mathbf{y}_{\xi_{k}})\in\mathbb{R}^{M+1}$;

5 $\alpha_{k}\leftarrow\frac{\alpha_{1}}{k}$;

6
$\bm{\mathbf{\theta}}^{(k)}\leftarrow\bm{\mathbf{\theta}}^{(k-1)}-\alpha_{k}g(\bm{\mathbf{\theta}}^{(k-1)};\mathbf{X}_{\xi_{k}},\mathbf{y}_{\xi_{k}})$;

7

8 end for

Algorithm 3 Minibatch SGD with uniform sampling when the covariance function
is the sum of multiple kernels

In the following, we present convergence guarantees for Algorithm 3 when
kernels exhibit exponential or polynomial eigendecay. The assumptions are
similar to the ones presented in Section 3.

######  Assumption B.1 (Bounded iterates)

Both $\bm{\mathbf{\theta}}^{*}$ and $\bm{\mathbf{\theta}}^{(k)}$ for $0\leq
k\leq K$ lie in $[\theta_{\min},\theta_{\max}]^{M+1}$, where
$0<\theta_{\min}<\theta_{\max}$.

###### Assumption B.2

For any $n>0$ and sample $\\{\bm{\mathbf{x}}_{i}\\}_{i=1}^{n}$, the kernel
matrices $\bm{\mathbf{K}}_{f,n}^{(1)},\dots,\bm{\mathbf{K}}_{f,n}^{(M)}$ share
the same eigenvectors.

######  Remark 7 (Explanation for Assumption B.2)

When extending the theoretical guarantees from $M=1$ to $M>1$, we find it
extremely challenging without Assumption B.2, which ensures that matrix
$\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}})$ and $\bm{\mathbf{K}}_{n}^{*}$ are
simultaneously diagonalizable, and thus facilitates the analysis for the
gradient. It remains an open question to establish theoretical results without
this assumption. We believe that if the eigenvectors of the kernel matrices
are “close” our results should still hold.

######  Assumption B.3 (Bounded stochastic gradient)

For all $0\leq k<K$,

$\|g(\bm{\mathbf{\theta}}^{(k)};\bm{\mathbf{X}}_{\xi_{k+1}},\bm{\mathbf{y}}_{\xi_{k+1}})\|_{2}\leq
G$

for some $G>0$.

###  B.1 Kernels with Exponential Eigendecay

######  Assumption B.4 (Exponential eigendecay)

For $1\leq i\leq M$, the eigenvalues of kernel function $k_{i}$ w.r.t.
probability measure $\mathbb{P}$ are $\\{C_{i}e^{-b_{i}j}\\}_{j=0}^{\infty}$,
where $0<b_{1}<b_{2}<\cdots<b_{M}$, and $C_{i}\leq 1$ are regarded as
constants.

######  Theorem B.1 (Convergence of parameter iterates, exponential
eigendecay)

Under Assumptions B.1 to B.4, when $m>C$ for some constant $C>0$, we have the
following results under two corresponding conditions on $s_{l}(m)$:

  1. 1.

If $s_{M+1}(m)=m$, initial step size $\alpha_{1}$ satisfies
$\frac{3}{2\gamma}\leq\alpha_{1}\leq\frac{2}{\gamma}$ where
$\gamma=\frac{1}{4\theta_{\max}^{2}}$, then for any
$0<\varepsilon<C\frac{\log\log m}{\log m}$, with probability at least
$1-CK\exp\\{-cm^{2\varepsilon}\\}$,

$(\theta^{(K)}_{M+1}-\theta^{*}_{M+1})^{2}\leq\frac{8G^{2}}{\gamma^{2}(K+1)}+Cm^{-\frac{1}{2}+\varepsilon}.$
(34)

  2. 2.

If in addition to $s_{M+1}(m)=m$, $s_{l}(m)$ is set as $\tau\log m$ for $1\leq
l\leq M$ where $\tau>C$, the eigendecay rates $b_{2}>2b_{1}$ when $M\geq 2$,
$\frac{3}{2\gamma}\leq\alpha_{1}\leq\frac{2}{\gamma}$ where $\gamma$ depends
on $\tau$, then for any $0<\varepsilon<\frac{1}{2}$, with probability at least
$1-CK\exp\\{-c(\log m)^{2\varepsilon}\\}$,

$(\theta^{(K)}_{1}-\theta^{*}_{1})^{2}+(\theta^{(K)}_{M+1}-\theta^{*}_{M+1})^{2}\leq\frac{8G^{2}}{\gamma^{2}(K+1)}+C(\log
m)^{-\frac{1}{2}+\varepsilon}.$ (35)

Here $c,C>0$ depend only on $M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$.

###### Remark B.1

Here we do not provide estimation error bounds for $\theta^{*}_{l}$ or
$\sigma_{f,l}^{2}$, for $2\leq l\leq M$, since they are associated with
kernels with faster eigendecay than $k_{1}$ and thus are not identifiable. The
technical condition on decay rates $b_{2}>2b_{1}$ is to ensure the convergence
of $\theta^{(K)}_{1}$ although $|\theta^{(K)}_{l}-\theta^{*}_{l}|$ does not
converge to $0$ for $2\leq l\leq M$.

###### Remark B.2

For the second case where $s_{l}(m)=\tau\log m$, $\tau$ and $\gamma$ need to
satisfy

$\tau>\frac{8b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{4}}{3b_{1}^{2}(b_{2}-2b_{1})\theta_{\min}^{4}},$
(36) $\gamma=\min\left\\{\frac{3(b_{2}-2b_{1})}{8\tau
b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}},\frac{1}{4\theta_{\max}^{2}}-\frac{2b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}}{3\tau
b_{1}^{2}(b_{2}-2b_{1})\theta_{\min}^{4}}\right\\}.$ (37)

######  Theorem B.2 (Convergence of full gradient, exponential eigendecay)

Under Assumptions B.1 to B.4, if
$\frac{3}{2\gamma}\leq\alpha_{1}\leq\frac{2}{\gamma}$ for
$\gamma=\frac{1}{4\theta_{\max}^{2}}$, $m>C$, $s_{M+1}(m)=m$, then for any
$0<\varepsilon<C\frac{\log\log m}{\log m}$, with probability at least
$1-CK\exp\\{-cm^{2\varepsilon}\\}$,

$\|\nabla\ell(\bm{\mathbf{\theta}}^{(K)})\|_{2}^{2}\leq
C\left[\frac{G^{2}}{K+1}+m^{-\frac{1}{2}+\varepsilon}\right],$ (38)

holds, where $c,C>0$ depend only on
$M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$.

###  B.2 Kernels with Polynomial Eigendecay

######  Assumption B.5 (Polynomial eigendecay)

For $1\leq l\leq M$, the eigenvalues of kernel function $k_{l}$ w.r.t.
probability measure $\mathbb{P}$ are $\\{C_{l}j^{-2b_{l}}\\}_{j=0}^{\infty}$,
where $\frac{\sqrt{21}+3}{4}<b_{1}<b_{2}<\cdots<b_{M}$, and $C_{l}\leq 1$ are
regarded as constants.

######  Theorem B.3 (Convergence of parameter iterates, polynomial eigendecay
)

Under Assumptions B.1 to B.3 and Assumption B.5, when $m>C$ for some constant
$C>0$, $s_{M+1}(m)=m$, $\frac{3}{2\gamma}\leq\alpha_{1}\leq\frac{2}{\gamma}$
where $\gamma=\frac{1}{8\theta_{\max}^{2}}$, then for any
$\varepsilon\in(\max\\{0,f_{1}(b_{1})\\},\frac{1}{2})$, with probability at
least $1-CKm^{-f_{2}(b_{1})\left[\varepsilon-
f_{1}(b_{1})\right]}-CK\exp\\{-cm^{2\varepsilon}\\}$,

$(\theta^{(K)}_{M+1}-\theta^{*}_{M+1})^{2}\leq\frac{8G^{2}}{\gamma^{2}(K+1)}+Cm^{-\frac{1}{2}+\varepsilon}.$
(39)

Here $c,C>0$ depend only on $M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$,
and $f_{1}(\cdot)$, $f_{2}(\cdot)$ are defined as in Theorem 3.3.

######  Theorem B.4 (Convergence of full gradient, polynomial eigendecay)

Under the same conditions as Theorem B.4, for any
$\varepsilon\in(\max\\{0,f_{1}(b_{1})\\},\frac{1}{2})$, with probability at
least $1-CK\left(m^{-f_{2}(b_{1})\left[\varepsilon-
f_{1}(b_{1})\right]}+\exp\\{-cm^{2\varepsilon}\\}\right)$,

$\|\nabla\ell(\bm{\mathbf{\theta}}^{(K)})\|_{2}^{2}\leq
C\left[\frac{G^{2}}{K+1}+m^{-\frac{1}{2}+\varepsilon}\right],$ (40)

holds, where $c,C>0$ depend only on
$M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$, $f_{1}(\cdot)$ and
$f_{2}(\cdot)$ are defined as in Theorem 3.3.

##  C Detailed Versions of Key Lemmas

First we present the detailed versions of the key lemmas (Lemma 4, 5 and 6)
discussed in Section 4, which are useful for proving the second part of
Theorem 3.1, Theorem 3.3 to Theorem B.4, and other supporting lemmas. We will
focus on the general case where the covariance function is a linear
combination of multiple kernels, introduced in Section 3.4. The model
introduced in Section 2 can be viewed a special case of this general model,
with $M=1$. Since we are considering the general setting in our proofs, the
notations are consistent with the ones introduced in Section 3.4: when $M=1$,
the only kernel function is referred to as $k_{1}(\cdot,\cdot)$ instead of
$k_{0}(\cdot,\cdot)$, we use $b_{1}$ to denote its eigendecay rate instead of
$b$.

######  Lemma 8 (Strongly convex-like property of
$g^{*}(\bm{\mathbf{\theta}}^{(k)})$, exponential eigendecay)

Under Assumptions 3.1 to 3.3,

  1. 1.

if $s_{M+1}(m)=m$, $m>C$ for some $C>0$, then with probability at least
$1-3MKm^{-c}$, the following claim holds true for $0\leq k<K$:

$\langle\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widetilde{g}_{k}^{*}\rangle\geq\frac{\gamma}{2}\|\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}-\varepsilon,$
(41)

where $\widetilde{\bm{\mathbf{\theta}}}^{(k)}=\theta^{(k)}_{M+1}$,
$\widetilde{\bm{\mathbf{\theta}}}^{*}=\theta^{*}_{M+1}$,
$\widetilde{g}_{k}^{*}=(g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{M+1}$,
$\gamma=\frac{1}{4\theta_{\max}^{2}}$, $\varepsilon=\frac{C\log m}{m}$;

  2. 2.

if $M\geq 2$, in addition to $s_{M+1}(m)=m$, we also have $s_{i}(m)=\tau\log
m$ for $1\leq i\leq M$, and $\tau$ satisfies (36), $b_{2}>2b_{1}$, then for
any
$0<\alpha<\min\\{\frac{2b_{1}+b_{2}}{2b_{1}},\frac{2b_{2}-4b_{1}}{14b_{1}+b_{2}}\\}$,
with probability at least $1-3MKm^{-\alpha}$, (41) holds for
$\widetilde{\bm{\mathbf{\theta}}}^{(k)}=(\theta^{(k)}_{1},\theta^{(k)}_{M+1})$,
$\widetilde{\bm{\mathbf{\theta}}}^{*}=(\theta^{*}_{1},\theta^{*}_{M+1})$,
$\widetilde{g}_{k}^{*}=((g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{1},(g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{M+1})^{\top}$,

$\gamma=\min\left\\{\frac{3(b_{2}-2b_{1})}{8\tau
b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}},\frac{1}{4\theta_{\max}^{2}}-\frac{2b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}}{3\tau
b_{1}^{2}(b_{2}-2b_{1})\theta_{\min}^{4}}\right\\};$ (42)

and $\varepsilon=C(\alpha+(\log m)^{-1})$;

  3. 3.

if $M=1$, in addition to $s_{M+1}(m)=m$, we also have $s_{1}(m)=\tau\log m$
where $\tau>\frac{64\theta_{\max}^{4}}{b_{1}\theta_{\min}^{4}}$, then with
probability at least $1-2Km^{-c}$, (41) holds for
$\widetilde{\bm{\mathbf{\theta}}}^{(k)}=\bm{\mathbf{\theta}}^{(k)}$,
$\widetilde{\bm{\mathbf{\theta}}}^{*}=\bm{\mathbf{\theta}}^{*}$,
$\widetilde{g}_{k}^{*}=g^{*}(\bm{\mathbf{\theta}})$,

$\gamma=\min\left\\{\frac{1}{32\tau
b_{1}\theta_{\max}^{2}},\frac{1}{4\theta_{\max}^{2}}-\frac{2\theta_{\max}^{2}}{\tau
b_{1}\theta_{\min}^{4}}\right\\}.$ (43)

and $\varepsilon=C\frac{\log m}{m}$.

Here $C>0$ depends only on $M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$.

######  Lemma 9 (Strongly convex-like property of
$g^{*}(\bm{\mathbf{\theta}}^{(k)})$, polynomial eigendecay)

If
$s_{M+1}(m)=m$, $m>C$ for some $C>0$, then for any
$0<\alpha<\frac{8b_{1}^{2}-12b_{1}-6}{4b_{1}+3}$, with probability at least
$1-MKm^{-\alpha}$, the following claim holds true for $0\leq k<K$:

$(g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{M+1}(\theta^{(k)}_{M+1}-\theta^{*}_{M+1})\geq\frac{\gamma}{2}(\theta^{(k)}_{M+1}-\theta_{M+1}^{*})^{2}-\epsilon,$
(44)

where $\gamma=\frac{1}{8\theta_{\max}^{2}}$,
$\epsilon=Cm^{-\frac{8b_{1}^{2}-12b_{1}-6-\alpha(4b_{1}+3)}{4b_{1}(2b_{1}-1)}}$.
Here $C>0$ depends only on $M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$.

Lemma 8 and Lemma 9 are detailed versions of Lemma 4.

######  Lemma 10 (Uniform statistical error)

For any $x>0$, $1\leq i\leq M+1$, we have

$\mathbb{P}\left(\sup_{\bm{\mathbf{\theta}}\in[\theta_{\min},\theta_{\max}]^{M+1}}\frac{n}{s_{i}(n)}\left|(\nabla\ell(\bm{\mathbf{\theta}}))_{i}-(\nabla\ell^{*}(\bm{\mathbf{\theta}}))_{i}\right|>Cx\right)\leq\delta(x).$
(45)

If Assumption 3.3 holds, $s_{i}(n)=\tau\log n$ for $\tau$ satisfying (36),
$n>C$ for some $C>0$, then

$\delta(x)\leq Cn^{-c}+C(\log x)^{2M+2}\exp\\{-c\log n\min\\{x^{2},x\\}\\}.$

If Assumption 3.3 or 3.4 hold, $s_{i}(n)=n$,

$\delta(x)\leq C(\log x)^{2M+2}\exp\\{-cn\min\\{x^{2},x\\}\\}.$

Here $c,C>0$ only depend on $M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$.

Lemma 10 is a detailed version of Lemma 5, including results for kernels with
exponential or polynomial eigendecay.

######  Lemma 11 (Detailed version of Lemma 6, exponential eigendecay)

Under Assumption 3.3, for any $\alpha>0$, if $n>C$ for $C>0$ depending on
$M,b_{1},\dots,b_{M}$, then with probability at least $1-3Mn^{-\alpha}$,

$\begin{split}\text{ if }l\wedge l^{\prime}\leq
M,\sum_{j=1}^{n}\frac{\lambda_{lj}\lambda_{l^{\prime}j}}{\left(\sum_{h=1}^{M+1}\theta_{h}\lambda_{hj}\right)^{2}}\leq&\frac{2(2+\alpha)}{b_{1}\theta_{\min}^{2}}\log
n,\\\ \frac{n-C(\alpha)\log
n}{4\theta_{\max}^{2}}\leq\sum_{j=1}^{n}\frac{\lambda_{M+1,j}^{2}}{\left(\sum_{h=1}^{M+1}\theta_{h}\lambda_{hj}\right)^{2}}\leq&\frac{n}{\theta_{\min}^{2}},\\\
\sum_{j=1}^{n}\frac{\lambda_{1j}\lambda_{M+1,j}}{\left(\sum_{h=1}^{M+1}\theta_{h}\lambda_{hj}\right)^{2}}\leq&\frac{5+2\alpha}{7b_{1}\theta_{\min}^{2}}\log
n,\end{split}$ (46)

holds for any $\bm{\mathbf{\theta}}\in[\theta_{\min},\theta_{\max}]^{M+1}$,
where $C(\alpha)>0$ depends only on $\alpha,b_{1}$. Furthermore,

  * •

if $M=1$, then for any $0<\alpha,\epsilon<1$, with probability at least
$1-2n^{-\alpha}$, in addition to (46) we have

$\sum_{j=1}^{n}\frac{\lambda_{1j}^{2}}{\left(\sum_{h=1}^{2}\theta_{h}\lambda_{hj}\right)^{2}}\geq\frac{\epsilon\log
n}{8b\theta_{\max}^{2}};$ (47)

  * •

if $b_{2}>2b_{1}$ holds, then for any $\frac{2b_{1}}{b_{2}}<\epsilon<1$,
$0<\alpha<\min\left\\{\epsilon,\frac{2\epsilon b_{2}-4b_{1}}{6b_{1}+\epsilon
b_{2}}\right\\}$, with probability at least $1-3Mn^{-\alpha}$, in addition to
(46) we have

$\begin{split}\sum_{j=1}^{n}\frac{\lambda_{1j}^{2}}{\left(\sum_{h=1}^{M+1}\theta_{h}\lambda_{hj}\right)^{2}}\geq&\frac{\epsilon(b_{2}-2b_{1})}{2b_{1}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}}\log
n,\\\ \text{for $1<i\leq M$, $l\in S_{i}$,
}\sum_{j=1}^{n}\frac{\lambda_{lj}\lambda_{ij}}{\left(\sum_{h=1}^{M+1}\theta_{h}\lambda_{hj}\right)^{2}}\leq&\frac{(6b_{1}+b_{2})\alpha\log
n}{2b_{1}(4b_{l}-b_{2}-6b_{1})\theta_{\min}^{2}}+\frac{C(\epsilon)}{\theta_{\min}^{2}};\end{split}$
(48)

as long as $n>C(\epsilon)$, where $C(\epsilon)>0$ depends on
$M,b_{1},\dots,b_{M}$ and $\epsilon$. Here $S_{i}=\\{1,i,i+1,\dots,M+1\\}$.

######  Lemma 12 (Detailed version of Lemma 6, polynomial eigendecay)

Under Assumption 3.4, for any
$0<\alpha<\frac{8b_{1}^{2}-12b_{1}-6}{4b_{1}+3}$, with probability at least
$1-Mn^{-\alpha}$, for any
$\bm{\mathbf{\theta}}\in[\theta_{\min},\theta_{\max}]^{M}$,

$\begin{split}\text{ if }l\leq l^{\prime}\leq
M,&\sum_{j=1}^{n}\frac{\lambda_{lj}\lambda_{l^{\prime}j}}{\left(\sum_{h=1}^{M+1}\theta_{h}\lambda_{hj}\right)^{2}}\leq
n^{\frac{(2+\alpha)(4b_{l}+3)}{4b_{l}(2b_{l}-1)}}\left(\frac{1}{\theta_{\min}^{2}}+\frac{a_{l}(4b_{l}+3)}{\theta_{\min}^{2}(4b_{l}^{2}-6b_{l}-3)}\right),\\\
&\frac{n-M\max_{l}a_{l}n^{\frac{(2+\alpha)(4b_{1}+3)}{4b_{1}(2b_{1}-1)}}}{4\theta_{\max}^{2}}\leq\sum_{j=1}^{n}\frac{\lambda_{M+1,j}^{2}}{\left(\sum_{h=1}^{M+1}\theta_{h}\lambda_{hj}\right)^{2}}\leq\frac{n}{\theta_{\min}^{2}},\end{split}$
(49)

where
$a_{l}=2\sqrt{2C_{l}}+\sqrt{\frac{2C_{l}}{2b_{l}-1}}+\frac{C_{l}}{2b_{l}-1}$.

##  D Proofs of Theorem 3.1, 3.3 and 3.4

Proof [Proof of Theorem 3.1] First we apply Lemma 8 under both cases of
$s_{i}(m)$: for the first case ($s_{M+1}(m)=m$) discussed in Lemma 8, define
$\widetilde{g}(\bm{\mathbf{\theta}}^{(k)})=(g(\bm{\mathbf{\theta}}^{(k)}))_{M+1}$,
and for the second case ($s_{i}(m)=\tau\log m$ and $s_{M+1}(m)=m$), define
$\widetilde{g}(\bm{\mathbf{\theta}}^{(k)})=((g(\bm{\mathbf{\theta}}^{(k)}))_{1},(g(\bm{\mathbf{\theta}}^{(k)})_{M+1})^{\top}$.
Then let
$\widehat{\bm{\mathbf{e}}}_{k}=\widetilde{g}(\bm{\mathbf{\theta}}^{(k)})-\widetilde{g}_{k}^{*}$.
Due to Lemma 8 and Assumption 3.2, we have

$\begin{split}\|\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}=&\|\widetilde{\bm{\mathbf{\theta}}}^{(k-1)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}-2\alpha_{k}\langle\widetilde{\bm{\mathbf{\theta}}}^{(k-1)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widetilde{g}(\bm{\mathbf{\theta}}^{(k-1)})\rangle+\alpha_{k}^{2}\|\widetilde{g}(\bm{\mathbf{\theta}}^{(k-1)})\|_{2}^{2}\\\
\leq&\|\widetilde{\bm{\mathbf{\theta}}}^{(k-1)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}(1-\alpha_{k}\gamma)+\alpha_{k}^{2}G^{2}+2\alpha_{k}\left(\varepsilon-\langle\widetilde{\bm{\mathbf{\theta}}}^{(k-1)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widehat{\bm{\mathbf{e}}}_{k-1}\rangle\right).\end{split}$
(50)

Recall that $\frac{3}{2\gamma}\leq\alpha_{1}\leq\frac{2}{\gamma}$, and
$\alpha_{k}=\frac{\alpha_{1}}{k}$ for all $k\geq 1$. Now we prove the
following statement for $k\geq 1$ by induction:

$\|\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}\leq\frac{2\alpha_{1}^{2}G^{2}}{k+1}+\sum_{i=0}^{k-1}\eta_{k,i}\left(\varepsilon-\langle\widetilde{\bm{\mathbf{\theta}}}^{(i)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widehat{\bm{\mathbf{e}}}_{i}\rangle\right),$
(51)

where $\eta_{k,i}=2\alpha_{i+1}\prod_{j=i+2}^{k}(1-\alpha_{j}\gamma)$. When
$k=1$, by (50) and the fact that $1-\alpha_{1}\gamma<0$,

$\|\widetilde{\bm{\mathbf{\theta}}}^{(1)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}\leq\alpha_{1}^{2}G^{2}+\eta_{1,0}\left(\varepsilon-\langle\widetilde{\bm{\mathbf{\theta}}}^{(0)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widehat{\bm{\mathbf{e}}}_{0}\rangle\right).$
(52)

Assuming (51) holds for $k=l\geq 1$, then due to (50) and the fact that
$1-\alpha_{l+1}\gamma\geq 0$ for $l\geq 1$, we have

$\begin{split}&\|\widetilde{\bm{\mathbf{\theta}}}^{(l+1)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}\\\
\leq&\left(\frac{2\alpha_{1}^{2}G^{2}}{l+1}+\sum_{i=0}^{l-1}\eta_{l,i}\left(\varepsilon-\langle\widetilde{\bm{\mathbf{\theta}}}^{(i)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widehat{\bm{\mathbf{e}}}_{i}\rangle\right)\right)(1-\alpha_{l+1}\gamma)+\alpha_{l+1}^{2}G^{2}\\\
&+2\alpha_{l+1}\left(\varepsilon-\langle\widetilde{\bm{\mathbf{\theta}}}^{(l)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widehat{\bm{\mathbf{e}}}_{l}\rangle\right)\\\
\leq&\frac{2\alpha_{1}^{2}G^{2}(l+1-\alpha_{1}\gamma)}{(l+1)^{2}}+\frac{\alpha_{1}^{2}G^{2}}{(l+1)^{2}}+\sum_{i=0}^{l}\eta_{l+1,i}\left(\varepsilon-\langle\widetilde{\bm{\mathbf{\theta}}}^{(i)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widehat{\bm{\mathbf{e}}}_{i}\rangle\right)\\\
\leq&\frac{2\alpha_{1}^{2}G^{2}}{l+2}+\sum_{i=0}^{l}\eta_{l+1,i}\left(\varepsilon-\langle\widetilde{\bm{\mathbf{\theta}}}^{(i)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widehat{\bm{\mathbf{e}}}_{i}\rangle\right).\end{split}$
(53)

Here the last two lines are due to range of $\alpha_{1}$ and the definitions
of $\eta_{l,i}$. The next step is to bound
$\sum_{i=0}^{K-1}\eta_{K,i}\left(\varepsilon-\langle\widetilde{\bm{\mathbf{\theta}}}^{(i)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widehat{\bm{\mathbf{e}}}_{i}\rangle\right)$.
First we have

$\begin{split}&\left|\sum_{i=0}^{K-1}\eta_{K,i}\left(\varepsilon-\langle\widetilde{\bm{\mathbf{\theta}}}^{(i)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widehat{\bm{\mathbf{e}}}_{i}\rangle\right)\right|\\\
\leq&\frac{2\alpha_{1}}{K}\sum_{i=0}^{K-1}\|\widetilde{\bm{\mathbf{\theta}}}^{(i)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}\|\widehat{\bm{\mathbf{e}}}_{i}\|_{2}+2\alpha_{1}\varepsilon\\\
\leq&C\left(\max_{0\leq i\leq
K-1}\|\widehat{\bm{\mathbf{e}}}_{i}\|_{2}+\varepsilon\right).\end{split}$ (54)

Note that the distribution of each minibatch
$\\{\bm{\mathbf{X}}_{\xi_{k+1}},\bm{\mathbf{y}}_{\xi_{k+1}}\\}_{i=1}^{m}$ is
the same as sampling $m$ independent $\bm{\mathbf{x}}_{i}$ from $\mathbb{P}$,
and then sampling
$\bm{\mathbf{y}}_{\xi_{k+1}}\sim\mathcal{N}(0,\bm{\mathbf{K}}^{*}_{\xi_{k+1}})$,
thus we can apply Lemma 10 on $\widetilde{g}(\bm{\mathbf{\theta}}^{(k)})$ and
$\widetilde{g}_{k}^{*}$. Combining Lemma 8, Lemma 10 and (51) leads to the
following conclusion.

  1. 1.

If $s_{M+1}(m)=m$, $m>C$, then for any $0<\varepsilon<\frac{1}{2}$, with
probability at least $1-CKm^{-c}-CK\exp\\{-cm^{2\varepsilon}\\}$,

$(\theta^{(K)}_{M+1}-\theta^{*}_{M+1})^{2}\leq\frac{8G^{2}}{\gamma^{2}(K+1)}+Cm^{-\frac{1}{2}+\varepsilon},$
(55)

where $\gamma=\frac{1}{4\theta_{\max}^{2}}$. Let $\varepsilon<C\frac{\log\log
m}{\log m}$, then $K\exp\\{-cm^{2\varepsilon}\\}\geq CKm^{-c}$, thus the
probability term is $1-CK\exp\\{-cm^{2\varepsilon}\\}$.

  2. 2.

If $M=1$, $s_{1}(m)=\tau\log m$, $s_{2}(m)=m$, $m>C$, then for any
$0<\varepsilon<\frac{1}{2}$, with probability at least

$1-CK\exp\\{-c(\log m)^{2\varepsilon}\\},$

we have

$(\theta^{(K)}_{1}-\theta^{*}_{1})^{2}+(\theta^{(K)}_{2}-\theta^{*}_{2})^{2}\leq\frac{8G^{2}}{\gamma^{2}(K+1)}+C(\log
m)^{-\frac{1}{2}+\varepsilon},$ (56)

where $\gamma=\min\left\\{\frac{1}{32\tau
b\theta_{\max}^{2}},\frac{1}{4\theta_{\max}^{2}}-\frac{2\theta_{\max}^{2}}{\tau
b\theta_{\min}^{4}}\right\\}$.

  3. 3.

If $s_{i}(m)=\tau\log m$ for $1\leq i\leq M$, $s_{M+1}(m)=m$, $m>C$,
$b_{2}>2b_{1}$, then for any
$0<\alpha<\min\\{\frac{2b_{1}+b_{2}}{2b_{1}},\frac{2b_{2}-4b_{1}}{14b_{1}+b_{2}}\\}$,$0<\varepsilon<\frac{1}{2}$,
with probability at least

$1-CKm^{-\alpha}-CK\exp\\{-c(\log m)^{2\varepsilon}\\},$

we have

$(\theta^{(K)}_{1}-\theta^{*}_{1})^{2}+(\theta^{(K)}_{M+1}-\theta^{*}_{M+1})^{2}\leq\frac{8G^{2}}{\gamma^{2}(K+1)}+C(\log
m)^{-\frac{1}{2}+\varepsilon}+C\alpha,$ (57)

where $\gamma$ is defined in (43). Let $c(\log
m)^{-1+2\varepsilon)}\alpha<C(\log m)^{-\frac{1}{2}+\varepsilon}$, then
$Km^{-\alpha}\geq\exp\\{-c(\log m)^{2\varepsilon}\\}$ and thus the probability
term can be written as $1-CK\exp\\{-c(\log m)^{2\varepsilon}\\}$ and the error
bound is $\frac{8G^{2}}{\gamma^{2}(K+1)}+C(\log
m)^{-\frac{1}{2}+\varepsilon}$.

Here $c,C>0$ depend only on $M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$.

Proof [Proof of Theorem 3.3] Define
$\widehat{e}_{k}=(g(\bm{\mathbf{\theta}}^{(k)}))_{M+1}-(g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{M+1}$.
Following similar arguments from the proof of Theorem 3.1 and applying Lemma
9, one can show that

$\begin{split}(\theta^{(k)}_{M+1}-\theta^{*}_{M+1})^{2}\leq\frac{2\alpha_{1}^{2}G^{2}}{k+1}+\sum_{i=0}^{k-1}\eta_{k,i}\left(\epsilon-(\theta^{(i)}_{M+1}-\theta^{*}_{M+1})\widehat{e}_{i}\right),\end{split}$
(58)

where $\eta_{k,i}=2\alpha_{i+1}\prod_{j=i+2}^{k}(1-\alpha_{j}\gamma)$,
$\gamma=\frac{1}{8\theta_{\max}^{2}}$ and
$\varepsilon=Cm^{-\frac{8b_{1}^{2}-12b_{1}-6-\alpha(4b_{1}+3)}{4b_{1}(2b_{1}-1)}}$.
Also note that

$\begin{split}&\left|\sum_{i=0}^{K-1}\eta_{K,i}\left(\epsilon-(\theta^{(i)}_{M+1}-\theta^{*}_{M+1})\widehat{e}_{i}\right)\right|\\\
\leq&2\alpha_{1}(\theta_{\max}-\theta_{\min})\max_{0\leq
i<K}\widehat{e}_{i}|+2\alpha_{1}\epsilon.\end{split}$ (59)

Similarly from the proof of Theorem 3.2, we can apply Lemma 10 on
$(g(\bm{\mathbf{\theta}}^{(k)}))_{M+1}$ and
$(g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{M+1}$. Therefore, combining (58) and
Lemma 10 leads to the following result: If $s_{M+1}(m)=m$, $m>C$, then for any
$0<\alpha<\frac{8b_{1}^{2}-12b_{1}-6}{4b_{1}+3}$, $0<\varepsilon<\frac{1}{2}$,
with probability at least $1-MKm^{-\alpha}-CK\exp\\{-cm^{2\varepsilon}\\}$,

$(\theta^{(K)}_{M+1}-\theta^{*}_{M+1})^{2}\leq\frac{8G^{2}}{\gamma^{2}(K+1)}+Cm^{-\frac{8b_{1}^{2}-12b_{1}-6-\alpha(4b_{1}+3)}{4b_{1}(2b_{1}-1)}}+Cm^{-\frac{1}{2}+\varepsilon},$
(60)

where $\gamma=\frac{1}{8\theta_{\max}^{2}}$. Here $c,C>0$ depend only on
$M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$. Let

$\frac{8b_{1}^{2}-12b_{1}-6-\alpha(4b_{1}+3)}{4b_{1}(2b_{1}-1)}=\frac{1}{2}-\varepsilon,$

we arrive at the final conclusion.

Proof [Proof of Theorem 3.4] Similarly from the proof of Theorem 3.2, we
utilize (25) and let $\lambda_{lj}$ be the $j$th largest eigenvalue of
$\bm{\mathbf{K}}_{f,n}^{(l)}$, $\lambda_{M+1,j}=1$. By (25) and Lemma 12, for
any $0<\alpha<\frac{8b_{1}^{2}-12b_{1}-6}{4b_{1}+3}$, with probability at
least $1-Mn^{-\alpha}$,

$\begin{split}\left|\left(\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})\right)_{i}\right|\leq&Cn^{-\frac{8b_{1}^{2}-12b_{1}-6-\alpha(4b_{1}+3)}{4b_{1}(2b_{1}-1)}},\end{split}$
(61)

for $1\leq i\leq M$, where $C>0$ depends on
$M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$. Meanwhile,

$\begin{split}\left|\left(\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})\right)_{M+1}\right|\leq
C\left(|\theta^{(k)}_{M+1}-\theta^{*}_{M+1}|+n^{-\frac{8b_{1}^{2}-12b_{1}-6-\alpha(4b_{1}+3)}{4b_{1}(2b_{1}-1)}}\right).\end{split}$
(62)

Thus we have

$\begin{split}\|\nabla\ell^{*}(\bm{\mathbf{\theta}}^{(k)})\|_{2}^{2}\leq&C\left[n^{-\frac{8b_{1}^{2}-12b_{1}-6-\alpha(4b_{1}+3)}{2b_{1}(2b_{1}-1)}}+(\theta^{(k)}_{M+1}-\theta^{*}_{M+1})^{2}\right].\end{split}$
(63)

By (45), Theorem 3.3 and Lemma 10, for any
$\varepsilon\in(\max\\{0,f_{1}(b_{1})\\},\frac{1}{2})$, if $m>C$, then with
probability at least $1-CK\left(m^{-f_{2}(b_{1})\left[\varepsilon-
f_{1}(b_{1})\right]}+\exp\\{-cm^{2\varepsilon}\\}\right)$, we have

$\|\nabla\ell(\bm{\mathbf{\theta}}^{(K)})\|_{2}^{2}\leq
C\left[\frac{G^{2}}{K+1}+m^{-\frac{1}{2}+\varepsilon}\right],$ (64)

where $c,C>0$ depend only on
$M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$.

##  E Proofs of Supporting Lemmas

Proof [proof of Lemma 8] Let $\lambda^{(k)}_{lj}$ be the $j$th eigenvalue of
$\bm{\mathbf{K}}_{f,\xi_{k+1}}^{(l)}$ for $1\leq l\leq M$, and
$\lambda^{(k)}_{M+1,j}=1$, then by the definition of
$g^{*}(\bm{\mathbf{\theta}}^{(k)})$, we have

$\begin{split}(g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{1}=&\frac{1}{2s_{1}(m)}\mathrm{tr}\Bigg{[}\bm{\mathbf{K}}_{\xi_{k+1}}(\bm{\mathbf{\theta}}^{(k)})^{-1}\left(\bm{\mathbf{I}}_{m}-\bm{\mathbf{K}}_{\xi_{k+1}}(\bm{\mathbf{\theta}}^{*})\bm{\mathbf{K}}_{\xi_{k+1}}(\bm{\mathbf{\theta}}^{(k)})^{-1}\right)\bm{\mathbf{K}}^{(1)}_{f,\xi_{k+1}}\Bigg{]}\\\
=&\frac{1}{2s_{1}(m)}\mathrm{tr}\Bigg{[}\bm{\mathbf{K}}_{\xi_{k+1}}(\bm{\mathbf{\theta}}^{(k)})^{-1}\left(\sum_{l=1}^{M}(\theta^{(k)}_{l}-\theta_{l}^{*})\bm{\mathbf{K}}^{(l)}_{f,\xi_{k+1}}+(\theta^{(k)}_{M+1}-\theta^{*}_{M+1})\bm{\mathbf{I}}_{m}\right)\\\
&\bm{\mathbf{K}}_{\xi_{k+1}}(\bm{\mathbf{\theta}}^{(k)})^{-1}\bm{\mathbf{K}}^{(1)}_{f,\xi_{k+1}}\Bigg{]}\\\
=&\frac{1}{2s_{1}(m)}\sum_{l=1}^{M+1}(\theta^{(k)}_{l}-\theta_{l}^{*})\sum_{j=1}^{m}\frac{\lambda_{lj}^{(k)}\lambda_{1j}^{(k)}}{\left(\sum_{l=1}^{M+1}\theta^{(k)}_{l}\lambda_{lj}^{(k)}\right)^{2}},\end{split}$
(65)

and

$\begin{split}(g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{M+1}=&\frac{1}{2m}\mathrm{tr}\Bigg{[}\bm{\mathbf{K}}_{\xi_{k+1}}(\bm{\mathbf{\theta}}^{(k)})^{-1}\left(\bm{\mathbf{I}}_{m}-\bm{\mathbf{K}}_{\xi_{k+1}}(\bm{\mathbf{\theta}}^{*})\bm{\mathbf{K}}_{\xi_{k+1}}(\bm{\mathbf{\theta}}^{(k)})^{-1}\right)\Bigg{]}\\\
=&\frac{1}{2m}\sum_{l=1}^{M+1}(\theta^{(k)}_{l}-\theta_{l}^{*})\sum_{j=1}^{m}\frac{\lambda_{lj}^{(k)}}{\left(\sum_{l=1}^{M+1}\theta^{(k)}_{l}\lambda_{lj}^{(k)}\right)^{2}}.\end{split}$
(66)

We prove Lemma 8 under two cases separately.

  1. 1.

$s_{i}(m)=\tau\log m$ for $1\leq i\leq M$, $s_{M+1}(m)=m$,
$\widetilde{\bm{\mathbf{\theta}}}^{(k)}=(\theta^{(k)}_{1},\theta^{(k)}_{M+1})^{\top}$,
$\widetilde{\bm{\mathbf{\theta}}}^{*}=(\theta^{*}_{1},\theta^{*}_{M+1})^{\top}$
and
$\widetilde{g}_{k}^{*}=((g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{1},(g^{*}(\bm{\mathbf{\theta}}^{(k)})_{M+1})^{\top}$.
Under this case, we can write
$\langle\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widetilde{g}_{k}^{*}\rangle$
as

$\begin{split}\langle\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widetilde{g}_{k}^{*}\rangle=(\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*})^{\top}\bm{\mathbf{A}}(\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*})+\varepsilon,\end{split}$

where each entry $A_{ij}$ of $\bm{\mathbf{A}}\in\mathbb{R}^{2\times 2}$ is
defined as follows:

$\begin{split}A_{11}=&\frac{1}{2\tau\log
m}\sum_{j=1}^{m}\frac{\lambda_{1j}^{(k)2}}{\left(\sum_{l=1}^{M+1}\theta^{(k)}_{l}\lambda_{lj}^{(k)}\right)^{2}},\\\
A_{12}=&\frac{1}{2\tau\log
m}\sum_{j=1}^{m}\frac{\lambda_{1j}^{(k)}}{\left(\sum_{l=1}^{M+1}\theta^{(k)}_{l}\lambda_{lj}^{(k)}\right)^{2}},\\\
A_{21}=&\frac{1}{2m}\sum_{j=1}^{m}\frac{\lambda_{1j}^{(k)}}{\left(\sum_{l=1}^{M+1}\theta^{(k)}_{l}\lambda_{lj}^{(k)}\right)^{2}},\\\
A_{22}=&\frac{1}{2m}\sum_{j=1}^{m}\frac{1}{\left(\sum_{l=1}^{M+1}\theta^{(k)}_{l}\lambda_{lj}^{(k)}\right)^{2}},\end{split}$

and

$\begin{split}\varepsilon=&\frac{1}{2\tau\log
m}\sum_{l=2}^{M}(\theta^{(k)}_{l}-\theta_{l}^{*})(\theta_{1}^{(k)}-\theta_{1}^{*})\sum_{j=1}^{m}\frac{\lambda_{lj}^{(k)}\lambda_{1j}^{(k)}}{\left(\sum_{l=1}^{M+1}\theta^{(k)}_{l}\lambda_{lj}^{(k)}\right)^{2}}\\\
&+\frac{1}{2m}\sum_{l=2}^{M}(\theta^{(k)}_{l}-\theta_{l}^{*})(\theta_{M+1}^{(k)}-\theta_{M+1}^{*})\sum_{j=1}^{m}\frac{\lambda_{lj}^{(k)}}{\left(\sum_{l=1}^{M+1}\theta^{(k)}_{l}\lambda_{lj}^{(k)}\right)^{2}},\end{split}$
(67)

for $M\geq 2$ and $\varepsilon=0$ for $M=1$. Note that the distribution of
each minibatch $\bm{\mathbf{X}}_{\xi_{k+1}}$ can be seen as $m$ independent
samples from $\mathbb{P}$, thus we can still apply Lemma 11, but substituting
$n$ by $m$.

     * •

When $M=1$, apply (47) in Lemma 11 with $\epsilon=\frac{1}{2}$, then for any
$0<\alpha<1$, with probability at least $1-2m^{-\alpha}$, we have

$\begin{split}A_{11}\geq\frac{1}{32\tau b_{1}\theta_{\max}^{2}},&\quad
A_{22}\geq\frac{1}{8\theta_{\max}^{2}}\left(1-\frac{C\log m}{m}\right),\\\
A_{12}\leq\frac{1}{2\tau b_{1}\theta_{\min}^{2}},&\quad A_{21}\leq\frac{\log
m}{2b_{1}\theta_{\min}^{2}m}.\end{split}$ (68)

Also note that for any $\omega>0$,

$\begin{split}&(\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*})^{\top}\bm{\mathbf{A}}(\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*})\\\
\geq&\left(A_{11}-\frac{(A_{12}+A_{21})\omega}{2}\right)(\theta^{(k)}_{1}-\theta_{1}^{*})^{2}+\left(A_{22}-\frac{(A_{12}+A_{21})}{2\omega}\right)(\theta^{(k)}_{2}-\theta_{2}^{*})^{2}\end{split}$
(69)

Let $\omega=\frac{\theta_{\min}^{2}}{16\theta_{\max}^{2}}$, then by (68) and
(69), one can show that

$\begin{split}&\langle\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widetilde{g}_{k}^{*}\rangle\\\
\geq&\frac{1}{64\tau
b\theta_{\max}^{2}}(\theta^{(k)}_{1}-\theta_{1}^{*})^{2}+\left(\frac{1}{8\theta_{\max}^{2}}-\frac{4\theta_{\max}^{2}}{\tau
b\theta_{\min}^{4}}\right)(\theta^{(k)}_{2}-\theta_{2}^{*})^{2}-C\frac{\log
m}{m}\\\
\geq&\frac{\gamma}{2}\|\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}-C\frac{\log
m}{m},\end{split}$ (70)

where

$\gamma=\min\left\\{\frac{1}{32\tau
b\theta_{\max}^{2}},\frac{1}{4\theta_{\max}^{2}}-\frac{2\theta_{\max}^{2}}{\tau
b\theta_{\min}^{4}}\right\\},$ (71)

and $C>0$ depends on $\theta_{\min},\theta_{\max},b_{1}$. It is guaranteed
that $\gamma>0$ Since we have assumed

$\tau>\frac{64\theta_{\max}^{4}}{b\theta_{\min}^{4}}.$

Therefore, if $m>C$, for any $0<\alpha<1$, with probability $1-2m^{-\alpha}$,
the following claims holds true:

$\langle\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widetilde{g}_{k}^{*}\rangle\geq\frac{\gamma}{2}\|\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}-\varepsilon,$
(72)

where $\varepsilon=C\frac{\log m}{m}$ for some constant $C>0$ depending on
$\theta_{\min},\theta_{\max},b_{1}$.

     * •

When $M\geq 2$ and $b_{2}>2b_{1}$, apply (48) in Lemma 11 with
$\epsilon=\frac{2b_{1}+b_{2}}{2b_{2}}$, then for any
$0<\alpha<\min\\{\frac{2b_{1}+b_{2}}{2b_{2}},\frac{2b_{2}-4b_{1}}{14b_{1}+b_{2}}\\}$,
with probability at least $1-3Mm^{-\alpha}$, we have

$\begin{split}|\varepsilon|\leq
CM(\theta_{\max}-\theta_{\min})^{2}\left(\frac{\alpha+1/\log
m}{\tau\theta_{\min}^{2}}+\frac{\log m}{m}\right),\end{split}$ (73)
$\begin{split}A_{11}\geq\frac{(2b_{1}+b_{2})(b_{2}-2b_{1})}{8\tau
b_{1}b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}},&\quad
A_{22}\geq\frac{1}{8\theta_{\max}^{2}}\left(1-\frac{C\log m}{m}\right),\\\
A_{12}\leq\frac{1}{2\tau b_{1}\theta_{\min}^{2}},&\quad A_{21}\leq\frac{\log
m}{2b_{1}\theta_{\min}^{2}m}.\end{split}$ (74)

Also note that for any $\omega>0$,

$\begin{split}&(\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*})^{\top}\bm{\mathbf{A}}(\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*})\\\
\geq&\left(A_{11}-\frac{(A_{12}+A_{21})\omega}{2}\right)(\theta^{(k)}_{1}-\theta_{1}^{*})^{2}+\left(A_{22}-\frac{(A_{12}+A_{21})}{2\omega}\right)(\theta^{(k)}_{M+1}-\theta_{M+1}^{*})^{2}\end{split}$
(75)

Let

$\omega=\frac{(2b_{1}+b_{2})(b_{2}-2b_{1})\theta_{\min}^{2}}{4b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}},$

then by (73), (74) and (75), one can show that

$\begin{split}&\langle\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widetilde{g}_{k}^{*}\rangle\\\
\geq&\frac{(2b_{1}+b_{2})(b_{2}-2b_{1})}{16\tau
b_{1}b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}}(\theta^{(k)}_{1}-\theta_{1}^{*})^{2}\\\
&+\left(\frac{1}{8\theta_{\max}^{2}}-\frac{b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}}{\tau
b_{1}(2b_{1}+b_{2})(b_{2}-2b_{1})\theta_{\min}^{4}}\right)(\theta^{(k)}_{M+1}-\theta_{M+1}^{*})^{2}\\\
&-\frac{CM(\theta_{\max}-\theta_{\min})^{2}}{\tau\theta_{\min}^{2}}(\alpha+(\log
m)^{-1})\\\
&-\frac{C(M+1)^{2}\theta_{\max}^{4}(\theta_{\max}-\theta_{\min})^{2}}{\theta_{\min}^{4}}\frac{\log
m}{m}\\\
\geq&\frac{\gamma}{2}\|\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}-C(\alpha+(\log
m)^{-1}),\end{split}$ (76)

where

$\gamma=\min\left\\{\frac{3(b_{2}-2b_{1})}{8\tau
b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}},\frac{1}{4\theta_{\max}^{2}}-\frac{2b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{2}}{3\tau
b_{1}^{2}(b_{2}-2b_{1})\theta_{\min}^{4}}\right\\},$ (77)

and $C>0$ depends on $M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$. $C$
does not depend on $\tau$ since we have assumed

$\tau>\frac{8b_{2}(b_{2}-b_{1})(M+1)^{2}\theta_{\max}^{4}}{3b_{1}^{2}(b_{2}-2b_{1})\theta_{\min}^{4}},$

which implies $\gamma>0$. Therefore, if $m>C$, for any
$0<\alpha<\min\\{\frac{2b_{1}+b_{2}}{2b_{2}},\frac{2b_{2}-4b_{1}}{14b_{1}+b_{2}}\\}$,
with probability $1-3Mm^{-\alpha}$, the following claims holds true:

$\langle\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widetilde{g}_{k}^{*}\rangle\geq\frac{\gamma}{2}\|\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*}\|_{2}^{2}-\varepsilon,$
(78)

where $\varepsilon=C(\alpha+(\log m)^{-1})$ for some constant $C>0$ depending
on $M,\theta_{\min},\theta_{\max},b_{1},\dots,b_{M}$.

  2. 2.

$s_{M+1}(m)=m$, $\widetilde{\bm{\mathbf{\theta}}}^{(k)}=\theta^{(k)}_{M+1}$,
$\widetilde{\bm{\mathbf{\theta}}}^{*}=\theta^{*}_{M+1}$ and
$\widetilde{g}_{k}^{*}=(g^{*}(\bm{\mathbf{\theta}}^{(k)}))_{M+1}$
Under this case, we can apply (46) in Lemma 11. Following similar arguments
from the first case, one can show that with probability at least $1-3Mm^{-c}$,

$\langle\widetilde{\bm{\mathbf{\theta}}}^{(k)}-\widetilde{\bm{\mathbf{\theta}}}^{*},\widetilde{g}_{k}^{*}\rangle\geq\frac{\gamma}{2}(\theta^{(k)}_{M+1}-\theta_{M+1}^{*})^{2}-\varepsilon,$
(79)

where $\gamma=\frac{1}{4\theta_{\max}^{2}}$, $\varepsilon=\frac{C\log m}{m}$,
if $m>C$. Here $C>0$ depends only on
$\theta_{\min},\theta_{\max},M,b_{1},\dots,b_{M}$.

Proof [proof of Lemma 10] Without loss of generality, we start from bounding
$\left(\nabla\ell(\bm{\mathbf{\theta}})\right)_{i}-\left(\nabla\ell^{*}(\bm{\mathbf{\theta}})\right)_{i}$
for an arbitrary $1\leq i\leq M+1$. By the definition of
$\nabla\ell(\bm{\mathbf{\theta}})$ and $\nabla\ell^{*}(\bm{\mathbf{\theta}})$,
we have

$\begin{split}&\left(\nabla\ell(\bm{\mathbf{\theta}})\right)_{i}-\left(\nabla\ell^{*}(\bm{\mathbf{\theta}})\right)_{i}\\\
=&-\frac{1}{2n}\left[\bm{\mathbf{y}}_{n}^{\top}\bm{\mathbf{K}}_{n}^{-1}(\bm{\mathbf{\theta}})\bm{\mathbf{K}}_{f,n}^{(i)}\bm{\mathbf{K}}_{n}^{-1}(\bm{\mathbf{\theta}})\bm{\mathbf{y}}_{n}-\mathrm{tr}\left(\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}})^{-1}\bm{\mathbf{K}}_{f,n}^{(i)}\bm{\mathbf{K}}_{n}(\bm{\mathbf{\theta}})^{-1}\bm{\mathbf{K}}_{n}^{*}\right)\right]\\\
=&-\left((\bm{\mathbf{K}}_{n}^{*})^{-\frac{1}{2}}\bm{\mathbf{y}}_{n}\right)^{\top}\bm{\mathbf{A}}(\bm{\mathbf{\theta}})\left((\bm{\mathbf{K}}_{n}^{*})^{-\frac{1}{2}}\bm{\mathbf{y}}_{n}\right)+\mathrm{tr}(\bm{\mathbf{A}}(\bm{\mathbf{\theta}})),\end{split}$
(80)

where
$\bm{\mathbf{A}}(\bm{\mathbf{\theta}})=\frac{1}{2n}\bm{\mathbf{K}}_{n}^{*\frac{1}{2}}\bm{\mathbf{K}}_{n}^{-1}(\bm{\mathbf{\theta}})\bm{\mathbf{K}}_{f,n}^{(i)}\bm{\mathbf{K}}_{n}^{-1}(\bm{\mathbf{\theta}})\bm{\mathbf{K}}_{n}^{*\frac{1}{2}}$.
Since we have assumed that the eigenvectors of
$\bm{\mathbf{K}}_{f,n}^{(1)},\dots,\bm{\mathbf{K}}_{f,n}^{(M)}$ are all the
same in Assumption B.2, we can write
$\bm{\mathbf{K}}_{f,n}^{(j)}=\bm{\mathbf{P}}^{\top}\bm{\mathbf{\Lambda}}_{j}\bm{\mathbf{P}}$
for all $j$, where $\bm{\mathbf{P}}$ is an orthogonal matrix and
$\bm{\mathbf{\Lambda}}_{j}$ is a diagonal matrix consisting of the eigenvalues
of $\bm{\mathbf{K}}_{f,n}^{(j)}$. Also let
$\bm{\mathbf{\Lambda}}_{M+1}=\bm{\mathbf{I}}_{n}$, then we have

$\bm{\mathbf{A}}(\bm{\mathbf{\theta}}^{(k)})=\bm{\mathbf{P}}^{\top}\frac{1}{2n}\left(\sum_{l=1}^{M+1}\theta_{l}^{*}\bm{\mathbf{\Lambda}}_{l}\right)\left(\sum_{l=1}^{M+1}\theta_{l}\bm{\mathbf{\Lambda}}_{l}\right)^{-2}\bm{\mathbf{\Lambda}}_{i}\bm{\mathbf{P}}.$
(81)

Let
$\bm{\mathbf{z}}_{n}=\bm{\mathbf{P}}(\bm{\mathbf{K}}_{n}^{*})^{-\frac{1}{2}}\bm{\mathbf{y}}_{n}$,
and
$\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}})=\left(\sum_{l=1}^{M+1}\theta_{l}^{*}\bm{\mathbf{\Lambda}}_{l}\right)\left(\sum_{l=1}^{M+1}\theta_{l}\bm{\mathbf{\Lambda}}_{l}\right)^{-2}\bm{\mathbf{\Lambda}}_{i}$,
where $\theta_{l}$ is the $l$th entry of $\bm{\mathbf{\theta}}$, then our goal
is to derive a bound for

$\begin{split}&\sup_{\bm{\mathbf{\theta}}\in[\theta_{\min},\theta_{\max}]^{M+1}}\frac{1}{2n}\left|\bm{\mathbf{z}}_{n}^{\top}\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}})\bm{\mathbf{z}}_{n}-\mathrm{tr}(\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}}))\right|.\end{split}$

We claim that there exists an $\varepsilon$-net
$\\{\bm{\mathbf{\theta}}_{\varepsilon}^{(1)},\dots,\bm{\mathbf{\theta}}_{\varepsilon}^{(N)}\\}$
of $[\theta_{\min},\theta_{\max}]^{M+1}$ under $\|\cdot\|_{\infty}$, with size
$N=(1+\frac{(\theta_{\max}-\theta_{\min})}{\varepsilon})^{M+1}$. That is to
say, for any $\bm{\mathbf{\theta}}\in[\theta_{\min},\theta_{\max}]^{M+1}$,
$\exists\widetilde{\bm{\mathbf{\theta}}}\in\\{\bm{\mathbf{\theta}}_{\varepsilon}^{(1)},\dots,\bm{\mathbf{\theta}}_{\varepsilon}^{(N)}\\}$
such that
$\bm{\mathbf{\Delta}}=\bm{\mathbf{\theta}}-\widetilde{\bm{\mathbf{\theta}}}$
satisfies$\|\bm{\mathbf{\Delta}}\|_{\infty}\leq\varepsilon$. The following
proof of this claim is very similar to the proof of Lemma 5.2 in Vershynin
(2010).

Define
$\bm{\mathbf{\theta}}_{c}=(\frac{\theta_{\min}+\theta_{\max}}{2},\dots,\frac{\theta_{\min}+\theta_{\max}}{2})\in\mathbb{R}^{M+1}$,
then an alternative way to represent $[\theta_{\min},\theta_{\max}]^{M+1}$ is
$\bm{\mathbf{\theta}}_{c}+B_{\infty}(\frac{\theta_{\max}-\theta_{\min}}{2})$.
Let
$\\{\bm{\mathbf{\theta}}^{(1)}_{\varepsilon},\dots,\bm{\mathbf{\theta}}^{(N)}_{\varepsilon}\\}$
be a maximal $\varepsilon$-separated subset of
$\bm{\mathbf{\theta}}_{c}+B_{\infty}(\frac{\theta_{\max}-\theta_{\min}}{2})$
(not the iterates of the SGD algorithm), which means that it is an
$\varepsilon$-net of
$\bm{\mathbf{\theta}}_{c}+B_{\infty}(\frac{\theta_{\max}-\theta_{\min}}{2})$,
and $\forall 1\leq i\neq j\leq N$,
$\|\bm{\mathbf{\theta}}_{\varepsilon}^{(i)}-\bm{\mathbf{\theta}}_{\varepsilon}^{(j)}\|_{\infty}\geq\varepsilon$.
Consider the $\ell_{\infty}$ balls with centers
$\\{\bm{\mathbf{\theta}}_{\varepsilon}^{(i)}\\}_{i=1}^{N}$ and radius
$\frac{\varepsilon}{2}$, then these balls are disjoint and are subsets of
$\bm{\mathbf{\theta}}_{c}+B_{\infty}(\frac{\theta_{\max}-\theta_{\min}+\varepsilon}{2})$.
Thus the sum of volumes of these balls is bounded by that of
$\bm{\mathbf{\theta}}_{c}+B_{\infty}(\frac{\theta_{\max}-\theta_{\min}+\varepsilon}{2})$,
which finishes the proof of

$N\leq\left(1+\frac{\theta_{\max}-\theta_{\min}}{\varepsilon}\right)^{M+1},$
(82)

In the following we linearize
$\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}})=\bm{\mathbf{\Lambda}}(\widetilde{\bm{\mathbf{\theta}}}+\bm{\mathbf{\Delta}})$
based on the Taylor series expression of each of its diagonal entries, so that
the upper bound for
$\left|\bm{\mathbf{z}}_{n}^{\top}\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}})\bm{\mathbf{z}}_{n}-\mathrm{tr}(\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}}))\right|$
can be implied by some bounds related to $\widetilde{\bm{\mathbf{\theta}}}$.
For any $1\leq j\leq m$, denote the $j$th diagonal entry of
$\bm{\mathbf{\Lambda}}_{l}$ by $\lambda_{lj}$ which is independent of
$\bm{\mathbf{\theta}}$, then the $j$th diagonal entry of
$\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}})$ can be written as follows:

$\begin{split}\bm{\mathbf{\Lambda}}_{jj}(\bm{\mathbf{\theta}})=\frac{\sum_{l=1}^{M+1}\lambda_{lj}\lambda_{ij}\theta_{l}^{*}}{\left(\sum_{l=1}^{M+1}\lambda_{lj}\theta_{l}\right)^{2}}.\end{split}$
(83)

Meanwhile, let $\Delta_{l}$ and $\widetilde{\theta}_{l}$ be the $l$th entry of
$\bm{\mathbf{\Delta}}$ and $\bm{\mathbf{\theta}}$, then one can show that

$\begin{split}\frac{1}{\left(\sum_{l=1}^{M+1}\lambda_{lj}\theta_{l}\right)^{2}}=&\frac{1}{\left(\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l}+\sum_{l=1}^{M+1}\lambda_{lj}\Delta_{l}\right)^{2}}\\\
=&\left(\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l}\right)^{-2}\left(1+\frac{\sum_{l=1}^{M+1}\lambda_{lj}\Delta_{l}}{\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l}}\right)^{-2}\\\
=&\left(\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l}\right)^{-2}\sum_{h=0}^{H-1}\frac{h+1}{\left(-\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l}\right)^{h}}\left(\sum_{l=1}^{M+1}\lambda_{lj}\Delta_{l}\right)^{h}\\\
&+\frac{H+1}{(1+\xi)^{H+2}\left(-\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l}\right)^{H+2}}\left(\sum_{l=1}^{M+1}\lambda_{lj}\Delta_{l}\right)^{H}\\\
=&\left(\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l}\right)^{-2}\left(\sum_{h_{1}+\dots+h_{M+1}\leq
H-1}\alpha^{(j)}_{h_{1},\dots,h_{M+1}}\prod_{l=1}^{M+1}\Delta_{l}^{h_{l}}+\text{RES}_{H}^{(j)}(\bm{\mathbf{\theta}})\right),\end{split}$
(84)

where the third equality holds if
$\left|\sum_{l=1}^{M+1}\lambda_{lj}\Delta_{l}\right|<\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l}$,
which is implied by $\|\bm{\mathbf{\Delta}}\|_{\infty}\leq\theta_{\min}$, and
we will choose $\varepsilon$ small enough to satisfy this. Here $\xi$ lies
between $0$ and
$\frac{\sum_{l=1}^{M+1}\lambda_{lj}\Delta_{l}}{\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l}}$,

$\begin{split}\alpha_{h_{1},\dots,h_{M+1}}^{(j)}=&\frac{(\sum_{l=1}^{M+1}h_{l}+1)!\prod_{l=1}^{M+1}\lambda_{lj}^{h_{l}}}{h_{1}!\cdots
h_{M+1}!(-\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l})^{\sum_{l=1}^{M+1}h_{l}}},\\\
\text{RES}_{H}^{(j)}(\bm{\mathbf{\theta}})=&\sum_{h_{1}+\cdots+h_{M+1}=H}\frac{(H+1)!\prod_{l=1}^{M+1}\lambda_{lj}^{h_{l}}\Delta_{l}^{h_{l}}}{h_{1}!\cdots
h_{M+1}!(1+\xi)^{H+2}(-\sum_{l=1}^{M+1}\lambda_{lj}\widetilde{\theta}_{l})^{H}}.\end{split}$
(85)

The quantities above satisfy

$\begin{split}|\alpha_{h_{1},\dots,h_{M+1}}^{(j)}|\leq\left(\sum_{l=1}^{M+1}h_{l}+1\right)\left(\frac{M+1}{\theta_{\min}}\right)^{\sum_{l=1}^{M+1}h_{l}},|\text{RES}_{H}^{(j)}(\bm{\mathbf{\theta}})|\leq&(H+1)\left(\frac{\varepsilon(M+1)}{\theta_{\min}}\right)^{H},\end{split}$
(86)

since $\sum_{h_{1}+\cdots+h_{M+1}=h}\frac{h!}{h_{1}!\cdots
h_{M+1}!}=(M+1)^{h}$. Define the following diagonal matrices:
$\bm{\mathbf{\Lambda}}^{(h_{1},\dots,h_{M+1})}(\widetilde{\bm{\mathbf{\theta}}}),\bm{\mathbf{\Lambda}}^{(H)}(\bm{\mathbf{\theta}})\in\mathbb{R}^{n\times
n}$ are with diagonal entries

$\bm{\mathbf{\Lambda}}^{(h_{1},\dots,h_{M+1})}_{jj}(\widetilde{\bm{\mathbf{\theta}}})=\alpha_{h_{1},\dots,h_{M+1}}^{(j)}\bm{\mathbf{\Lambda}}_{jj}(\widetilde{\bm{\mathbf{\theta}}}),\bm{\mathbf{\Lambda}}^{(H)}_{jj}(\bm{\mathbf{\theta}})=\text{RES}_{H}^{(j)}(\bm{\mathbf{\theta}})\bm{\mathbf{\Lambda}}_{jj}(\widetilde{\bm{\mathbf{\theta}}}).$
(87)

Then we can write

$\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}})=\sum_{h_{1}+\cdots+h_{M+1}\leq
H-1}\bm{\mathbf{\Lambda}}^{(h_{1},\dots,h_{M+1})}(\widetilde{\bm{\mathbf{\theta}}})\prod_{l=1}^{M+1}\Delta_{l}^{h_{l}}+\bm{\mathbf{\Lambda}}^{(H)}(\bm{\mathbf{\theta}}),$

and thus

$\begin{split}&\left|\bm{\mathbf{z}}_{n}^{\top}\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}})\bm{\mathbf{z}}_{n}-\mathrm{tr}(\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}}))\right|\\\
\leq&\max_{1\leq k\leq N}\sum_{h_{1}+\cdots+h_{M+1}\leq
H-1}\varepsilon^{\sum_{l=1}^{M+1}h_{l}}\left|\bm{\mathbf{z}}_{n}^{\top}\bm{\mathbf{\Lambda}}^{(h_{1},\dots,h_{M+1})}(\bm{\mathbf{\theta}}_{\varepsilon}^{(k)})\bm{\mathbf{z}}_{n}-\mathrm{tr}(\bm{\mathbf{\Lambda}}^{(h_{1},\dots,h_{M+1})}(\bm{\mathbf{\theta}}_{\varepsilon}^{(k)}))\right|\\\
&+\left|\bm{\mathbf{z}}_{n}^{\top}\bm{\mathbf{\Lambda}}^{(H)}(\bm{\mathbf{\theta}})\bm{\mathbf{z}}_{n}-\mathrm{tr}(\bm{\mathbf{\Lambda}}^{(H)}(\bm{\mathbf{\theta}}))\right|.\end{split}$
(88)

In order to provide an upper bound for the first term above, we first bound

$\left|\bm{\mathbf{z}}_{n}^{\top}\bm{\mathbf{\Lambda}}^{(h_{1},\dots,h_{M+1})}(\bm{\mathbf{\theta}}_{\varepsilon}^{(k)})\bm{\mathbf{z}}_{n}-\mathrm{tr}(\bm{\mathbf{\Lambda}}^{(h_{1},\dots,h_{M+1})}(\bm{\mathbf{\theta}}_{\varepsilon}^{(k)}))\right|$

for an arbitrary $k$. First note that for any $1\leq k\leq N$,

$\|\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}}_{\varepsilon}^{(k)})\|_{2}=\max_{j}\frac{\sum_{l=1}^{M+1}\lambda_{lj}\lambda_{ij}\theta_{l}^{*}}{\left(\sum_{l=1}^{M+1}\lambda_{lj}\theta^{(k)}_{l}\right)^{2}}\leq\frac{\theta_{\max}}{\theta_{\min}^{2}}.$
(89)

While for
$\|\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}}_{\varepsilon}^{(k)})\|_{F}^{2}$,
one can show that

$\begin{split}\|\bm{\mathbf{\Lambda}}(\bm{\mathbf{\theta}}_{\varepsilon}^{(k)})\|_{F}^{2}\leq&\theta_{\max}^{2}\sum_{j=1}^{n}\frac{(\sum_{l=1}^{M+1}\lambda_{lj}\lambda_{ij})^{2}}{(\sum_{l=1}^{M+1}\lambda_{lj}\theta^{(k)}_{l})^{4}}\\\
\leq&C\sum_{j=1}^{n}\frac{\sum_{l=1}^{M+1}\lambda_{lj}\lambda_{ij}}{(\sum_{l=1}^{M+1}\lambda_{lj}\theta_{l})^{2}}.\end{split}$
(90)

Let