# A Sieve Stochastic Gradient Descent Estimator for Online Nonparametric
Regression in Sobolev ellipsoids

Tianyu Zhang  Noah Simon

###### Abstract

The goal of regression is to recover an unknown underlying function that best
links a set of predictors to an outcome from noisy observations. In
nonparametric regression, one assumes that the regression function belongs to
a pre-specified infinite-dimensional function space (the hypothesis space). In
the online setting, when the observations come in a stream, it is
computationally-preferable to iteratively update an estimate rather than
refitting an entire model repeatedly. Inspired by nonparametric sieve
estimation and stochastic approximation methods, we propose a sieve stochastic
gradient descent estimator (Sieve-SGD) when the hypothesis space is a Sobolev
ellipsoid. We show that Sieve-SGD has rate-optimal mean squared error (MSE)
under a set of simple and direct conditions. The proposed estimator can be
constructed with a low computational (time and space) expense: We also
formally show that Sieve-SGD requires almost minimal memory usage among all
statistically rate-optimal estimators.

##  1 Introduction

It is commonly of interest to understand the association between a number of
features (or predictors) and a quantitative outcome. To this end, one often
estimates an underlying regression function that best links these two
quantities from noisy observations. More formally, suppose we obtain $n$
samples, $\left(X_{i},Y_{i}\right)$, where
$X_{i}\in\mathcal{X}\subset\mathbb{R}^{p}$ denotes a $p$-vector of features
from the $i$-th sample we observe, and $Y_{i}\in\mathbb{R}$ denotes the $i$-th
outcome. Further suppose that each pair $\left(X_{i},Y_{i}\right)$ is
independently and identically distributed (i.i.d.) from a fixed but unknown
distribution $\rho$ over
$\mathcal{X}\times\mathbb{R}\subset\mathbb{R}^{p}\times\mathbb{R}$. A common
target of estimation is the conditional mean $f_{\rho}(X):=E_{\rho}[Y|X]$.
Under extremely mild conditions, this conditional mean is the optimal function
for predicting $Y$ from $X$ with regard to mean squared error. More formally,

$f_{\rho}=\operatorname{argmin}_{f\in
L^{2}_{\rho_{X}}}E_{\rho}\left[\left(Y-f(X)\right)^{2}\right],$ (1)

where $L^{2}_{\rho_{X}}$ is the collection of all $\rho_{X}$-mean square
integrable functions and $\rho_{X}$ is the marginal distribution of $X$. Our
goal is to estimate $f_{\rho}$ from our collection of observed data.

In order to make a tractable estimation of $f_{\rho}$ from data, we need to
make additional assumptions on its smoothness/structure: The entire
$L^{2}_{\rho_{X}}$ space is too big to search within [4, 35]. We often
formally assume that $f_{\rho}$ belongs to a pre-specified function space
$\mathcal{F}\subsetneqq L^{2}_{\rho_{X}}$. This $\mathcal{F}$ is known as the
_hypothesis space_ of the regression problem.

If $\mathcal{F}$ can be indexed by a finite-dimensional parameter set
$\Theta\subset\mathbb{R}^{d}$, $d\in\mathbb{N}^{+}$, we refer to $\mathcal{F}$
as a _parametric function space_ or a _parametric class_. One common
parametric class is
$\mathcal{F}=\\{X^{\top}\beta\,|\,\beta\in\mathbb{R}^{d}\\}$, the set of all
linear functions of $X$. Parametric classes can impose overly restrictive
assumptions on the form of the regression function that may not be realistic
in practice. As such, it has become popular to assume less restrictive
structure: It is common to define the hypothesis space based on constraints on
derivatives, monotonicity, or other shape-related properties. Such an
$\mathcal{F}$ is most naturally written as an infinite-dimensional subset of
$L^{2}_{\rho_{X}}$. Commonly used examples of $\mathcal{F}$ in the statistics
community include Hölder balls, Sobolev spaces [26, 44, 61], reproducing
kernel Hilbert spaces (RKHS) [5, 13] and Besov spaces [28]. These are known as
_nonparametric function spaces_ , as they cannot naturally be parametrized
using a finite length vector. The Sobolev ellipsoid, in particular, is a
simple and useful abstraction of many important function spaces [61].
Therefore, we focus on them exclusively as the hypothesis spaces in this
paper.

In this paper, we propose an estimator for _online_ nonparametric regression.
In online estimation, the data are seen sequentially, one sample at a time.
After each sample is observed, our estimate of $f_{\rho}$ must be updated, as
a prediction may be required at any point in time before all the available
samples are processed. In an online problem with $n$ observations, we must
sequentially construct $n$ estimates. This is in contrast to the classical
batch learning setting where we collect all the data initially and perform
estimation only once. In the online setting, it is generally computationally
infeasible to repeatedly refit the whole model from scratch for each new
observation. Thus, online algorithms are generally carefully developed to
permit more tractable _updates_ after each new observation [17, 34].

An ideal estimator in online settings should be: i) statistically rate-
optimal, i.e. achieve the minimax-rate for estimating $f_{\rho}$ over
$\mathcal{F}$; and ii) computationally inexpensive to construct/update. In
this paper, we present such an online nonparametric estimator for use when the
hypothesis space is a Sobolev ellipsoid, which we term the _Sieve Stochastic
Gradient Descent estimator (Sieve-SGD)_. This method can be thought of as an
online version of the classical projection estimator [58], where the latter is
a specific example of sieve estimators [25, 51]. We use the more general term
“sieve” in naming our method to emphasize its nonparametric nature and avoid
confusion with the term “stochastic projection” [59]. We will show that Sieve-
SGD can achieve rate-optimal estimation error for $\mathcal{F}$ a Sobolev
ellipsoid and asymptotically uses minimal memory (up to a log factor) among
all rate-optimal estimators. In addition, our estimator has the same
computational cost (up to a constant) as merely examining each allocated
memory location every time a new sample $X_{i}$ is collected. This intimates
that in scenarios when our estimator has near optimal space complexity, it may
also have near optimal time complexity (though formal investigation of lower
bounds for time complexity in this problem is beyond the scope of the current
manuscript).

The structure of our paper continues as follows. In Section 2 we briefly cover
classical results for batch, nonparametric estimation in Sobolev ellipsoids,
focusing on projection estimators (which motivate our method). In Section 3 we
return to the online setting and explore intuition for how one might combine
projection estimation and stochastic gradient descent (SGD) [7]. The latter is
a well-studied method that has been applied fruitfully to online parametric
regression problems. This will help motivate our proposed method, which, as we
will see, can be thought of as an SGD estimator with a parameter space of
increasing dimension. In Section 4 we discuss existing nonparametric SGD
estimators, and identify some notable drawbacks of current methods. In Section
5, we introduce the formal construction of Sieve-SGD and analyze its
computational expense. From there, we show that our estimator has a
dramatically smaller “dimension” than existing methods and discuss how this
helps to reduce the computational expense. In Section 6, we give a theoretical
analysis of the statistical properties of Sieve-SGD. In constructing our
estimator, we need to decide how quickly to grow the dimension it projects
onto. Under minimal assumptions, we characterize the required growth rate and
learning rate for our estimator to be statistically and computationally (near)
optimal. We will also investigate under what conditions such an optimality
result is adaptive/insensitive to our choice of the “dimension-specific
learning rate”. Section 7 provides simulation studies to illustrate our
theoretical results. Finally, in Section 8, we have some further discussion of
Sieve-SGD and possible future research directions.

Notation: In this paper, we use $C$ to denote a generic constant that does not
depend on sample size $n$ (The value of $C$ may be different in different
parts of the manuscript). Additionally the notation $a_{n}=\Theta(b_{n})$
means $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$. The function $\lfloor x\rfloor$
maps $x$ to the largest integer smaller than $x$. For a vector
$x\in\mathbb{R}^{p}$, $x^{(i)}$ is the $i$-th component of $x$. The notation
$x\vee y$ (resp. $x\wedge y$) is shorthand for $\max\\{x,y\\}$ (resp.
$\min\\{x,y\\}$). The $\|\cdot\|_{\infty}$ norm of a continuous function $f$
is defined as $\|f\|_{\infty}=\sup_{x\in\mathcal{X}}|f(x)|$, where
$\mathcal{X}$ is the domain of $f$.

##  2 Batch Learning and the Projection Estimator

In this section we consider estimation in the classical batch setting where
our estimate is constructed once after all $n$ samples are observed. We will
begin by formally introducing a Sobolev ellipsoid: This is the hypothesis
space we will use throughout this manuscript. This will be followed by
presenting the classical projection estimator [58].

Consider a user-specified measure $\nu$ whose support contains $\mathcal{X}$,
and the corresponding square-integrable function space $L^{2}_{\nu}$. In many
interesting cases $\nu$ can be simply taken as Lebesgue measure over
$\mathcal{X}$ but it is not necessary in the general form of our theory. To
define a Sobolev ellipsoid in $L^{2}_{\nu}$, suppose we have a complete
orthonormal basis $\\{\psi_{j},j=1,2,...\\}\subset L^{2}_{\nu}$ of
$L^{2}_{\nu}$ [30]. This means

  1. i)

For any $f\in L^{2}_{\nu}$, there exists a unique sequence
$(\theta_{j})_{j=1}^{\infty}\in\ell^{2}$ such that

$\lim_{N\rightarrow\infty}\int\left|f(z)-\sum_{j=1}^{N}\theta_{j}\psi_{j}(z)\right|^{2}d\nu(z)=0\quad\text{(completeness)}$
(2)

where $\ell^{2}$ is the space of square convergent series.

  2. ii)

$\\{\psi_{j}\\}$ is an orthonormal system:

$\int\psi_{i}(z)\psi_{j}(z)d\nu(z)=\delta_{ij}\quad\text{(orthonormality)}$
(3)

where $\delta_{ij}$ is the Kronecker delta.

We define the _Sobolev ellipsoid_ $W(s,Q,\\{\psi_{j}\\})$ as:

$W\left(s,Q,\\{\psi_{j}\\}\right):=\left\\{\left.f=\sum_{j=1}^{\infty}\theta_{j}\psi_{j}\
\right|\ \sum_{j=1}^{\infty}\left(\theta_{j}j^{s}\right)^{2}\leq
Q^{2}\right\\}$ (4)

We refer to $(\theta_{j})_{j=1}^{\infty}$ as the (general) Fourier
coefficients of a function $f$. Throughout this manuscript, we assume the
measure $\nu$, basis functions $\psi_{j}$ and the regularity parameter $s$ are
all known. When it is clear which $\psi_{j}$ we are using, we will denote a
Sobolev ellipsoid simply by $W(s,Q)$. We may also use the further simplified
notation $W(s)$ because the diameter $Q$ usually plays a secondary role in our
theoretical analysis and the proposed method is adaptive to it. Intuitively,
by saying a function $f$ belongs to a Sobolev ellipsoid, we are requiring its
coefficients $\\{\theta_{j}\\}$ to converge to zero faster than $j^{-(s+1/2)}$
(if not, the sum $\sum_{j=1}^{\infty}\left(\theta_{j}j^{s}\right)^{2}$ would
diverge to infinity). The larger $s$ is, the faster the decay of $\theta_{j}$
will be, and thus the stronger our assumption is.

Sobolev ellipsoids are popular spaces to study for two reasons: 1) They impose
a useful structure for theories and computations, especially as a basic
example of hypothesis spaces with finite metric entropy; and 2) Many natural
spaces of regular functions are Sobolev ellipsoids. For example, if
$\mathcal{X}=[0,1]$ with $\nu$ as Lebesgue measure, then for any $s>0$, the
periodic Sobolev space

$\mathcal{F}=\left\\{f\in
L^{2}_{\nu}\mid\int\left(f^{(s)}(x)\right)^{2}dx<Q^{2},f^{(k)}(0)=f^{(k)}(1),k=0,1,..,s-1\right\\}$
(5)

can be written as a Sobolev ellipsoid, using an orthogonal basis of
trigonometric functions [61, Chapter 2]. More generally, for any RKHS
$\left(\mathcal{H},\langle\cdot,\cdot\rangle_{\mathcal{H}}\right)$, it is
possible to find a set of $\psi_{j}$ such that
$W(s,Q,\\{\psi_{j}\\})=\\{f\in\mathcal{H}\ |\ \|f\|_{\mathcal{H}}\leq Q\\}$ ,
i.e. a ball in an RKHS is a Sobolev ellipsoid (see [15, 55]).

In everything that follows we will assume that $f_{\rho}$, our target of
estimation, lives in a known Sobolev ellipsoid
$W\left(s,Q,\\{\psi_{j}\\}\right)$; with $\\{\psi_{j}\\}$ specified, and
orthonormal w.r.t. a specified measure $\nu$ (not necessarily equal to
$\rho_{X}$); and $s$ known (we allow $Q$ to be unknown).

The _Projection Estimator_ is a classical estimator naturally associated with
a Sobolev ellipsoid. We can treat it as a special case of general sieve
estimation [25, Chapter 10]: The estimates can be characterized by a sequence
of finite dimensional linear spaces of increasing dimension (the dimension
increases with sample size). For any given $f\in W(s,Q)$, the magnitude of its
Fourier coefficients must asymptotically decrease with $j$ fast enough. Thus,
it might be sensible to consider an estimator that discards the basis
functions far into the tail. This is precisely what the projection estimator
does. More formally, for a user-specified truncation level $J_{n}$, the
projection estimator is given by

$\hat{f}_{n,J_{n}}=\sum_{j=1}^{J_{n}}\hat{\theta}_{j}\psi_{j}$ (6)

where $\hat{\theta}=(\hat{\theta}_{1},..,\hat{\theta}_{J_{n}})^{\top}$ is the
solution of the least square problem:

$\min_{\theta\in\mathbb{R}^{J_{n}}}\sum_{i=1}^{n}\left(Y_{i}-\sum_{j=1}^{J_{n}}\theta_{j}\psi_{j}(X_{i})\right)^{2}$
(7)

It has been shown (e.g. [58], Theorem 1.9) that when we choose
$J_{n}=\Theta(n^{\frac{1}{2s+1}})$, the projection estimator is a rate-optimal
estimator over $W(s,Q)$, i.e.

$\limsup_{n\rightarrow\infty}\sup_{f_{\rho}\in
W(s,Q)}E\left[\left\|\hat{f}_{n,J_{n}}-f_{\rho}\right\|_{2}^{2}\right]=O(n^{-\frac{2s}{2s+1}})$
(8)

This result is usually shown in the literature for $X_{i}$ equally spaced, or
drawn from a uniform distribution. But in our theoretical analysis (Section
6), we allow $\rho_{X}$ to be a much more general distribution.

Sieve-SGD is inspired by this (batch) projection estimator. The key here is
that the number of basis functions we need to use can be dramatically smaller
than sample size, and their analytical forms do not depend on the data
(usually reproducing kernel methods use basis functions “centered” at the
feature vectors $X_{i}$). This possibility has been rarely explored [68] by
existing nonparametric online estimation research.

##  3 Online Learning and Stochastic Approximation

We now move to the online learning setting where observations are collected
sequentially from a data stream, and an estimate of our function is required
after each sample. Such an infinite data stream may really exist, for example,
with simulated samples as in reinforcement learning. Or the stream may serve
as an abstraction used with large-scale data sets where it is not favorable to
handle all the samples at once. It is generally computationally prohibitive to
use a method developed for the “batch” setting and completely refit it after
each observation. Instead methods that iteratively update are preferred. For
example, fitting a single projection estimator (solving (7)) with $n$
observations using $J_{n}=n^{\frac{1}{2s+1}}$ requires computation of
$\Theta(n^{1+\frac{2}{2s+1}})$. Refitting a projection estimator (from
scratch) after each observation $i=1,\ldots,n$ with $J_{i}=\lfloor
i^{\frac{1}{2s+1}}\rfloor$ would require an accumulated computation of
$\sum_{i=1}^{n}i^{1+\frac{2}{2s+1}}=\Theta(n^{2+\frac{2}{2s+1}})$. This scales
worse than quadratically in $n$. Our goal in the online nonparametric setting
is to find a statistically rate-optimal estimator whose computation scales
only slightly worse than linearly in $n$.

Online learning has been thoroughly studied for parametric $\mathcal{F}$. Many
proposed methods are based on the concept of _stochastic approximation_ [34].
One of the most popular methods in stochastic approximation is _Stochastic
Gradient Descent_ (SGD) [7]. In the parametric setting, SGD gives a
statistically rate-optimal estimator $\hat{f}_{n}$ whose population mean-
square-error $E\|\hat{f}_{n}-f_{\rho}\|_{L^{2}_{\rho_{X}}}^{2}$ is of order
$O(\frac{1}{n})$ [2, 3, 21]. Both vanilla SGD and its variants have been
applied to general convex loss functions and are shown to be statistically
rate-optimal under mild conditions [17].

###  3.1 Parametric SGD

To motivate stochastic optimization in the nonparametric setting, we first
give more details on SGD for parametric classes. Here we consider a specific
class of functions
$\mathcal{F}=\\{f=\sum_{j=1}^{d}\beta^{(j)}\psi_{j},\beta\in\mathbb{R}^{d}\\}$
for a set of pre-specified basis functions
$\psi_{j}:\mathbb{R}^{p}~{}\rightarrow~{}\mathbb{R},j=1,\ldots,d$. We use this
example to illustrate the principle of (parametric) SGD. Solving
$\operatorname{argmin}_{f\in\mathcal{F}}E\left[(Y-f(X))^{2}\right]$ reduces to
solving

$\min_{\beta\in\mathbb{R}^{d}}\ell(\beta):=\min_{\beta\in\mathbb{R}^{d}}E\left[\left(Y-\sum_{j=1}^{d}\beta^{(j)}\psi_{j}(X)\right)^{2}\right]$
(9)

We assume the minimizer of $\ell(\beta)$ exists and denote it as $\beta^{*}$.

If we knew the true joint distribution $\rho$ of $(X,Y)$ (which never happens
in practice), then equation (9) is just a numerical optimization problem which
does not involve data. We could use gradient descent to solve it. The gradient
of $\ell$ at any point $\beta$ is

$\displaystyle\nabla\ell(\beta)=-2E\left[\left(Y-\sum_{j=1}^{d}\beta^{(j)}\psi_{j}(X)\right)\left(\psi_{1}(X),...,\psi_{d}(X)\right)^{\top}\right]$
(10)

Thus, the gradient descent updating rule one could use is:

$\displaystyle\hat{\beta}_{0}$ $\displaystyle=0$ (11)
$\displaystyle\hat{\beta}_{n}$
$\displaystyle=\hat{\beta}_{n-1}-\gamma_{n}\nabla\ell(\hat{\beta}_{n-1})$

where $\\{\gamma_{n}\\}$ is a pre-specified sequence of step-sizes (or
learning rate) and $\hat{\beta}_{n}\in\mathbb{R}^{d}$ is the sequence of
approximations of $\beta^{*}$.

In practice, we do not know the joint distribution $\rho$: we must use data to
estimate $\beta^{*}$. In the framework of SGD, this is done by using the data
to get unbiased estimates of the gradients and substituting the estimates into
our updating rule (11). In particular we note that
$\widehat{\nabla\ell(\beta)}:=-2\left(Y_{i}-\sum_{j=1}^{d}\beta^{(j)}\psi_{j}(X_{i})\right)(\psi_{1}(X_{i}),...,\psi_{d}(X_{i}))^{\top}$
is an unbiased estimator of the gradient $\nabla\ell(\beta)$ based on one
sample. This results in the SGD updating rule.

$\displaystyle\hat{\beta}_{0}$ $\displaystyle=0$ (12)
$\displaystyle\hat{\beta}_{n}$
$\displaystyle=\hat{\beta}_{n-1}-\gamma_{n}\widehat{\nabla\ell(\hat{\beta}_{n-1})}$
$\displaystyle=\hat{\beta}_{n-1}+2\gamma_{n}\left(Y_{n}-\sum_{j=1}^{d}\hat{\beta}_{n-1}^{(j)}\psi_{j}(X_{n})\right)(\psi_{1}(X_{n}),...,\psi_{d}(X_{n}))^{\top}$

So our estimator $\hat{f}_{n}$ of $f_{\rho}$ has the following functional
update rule, derived from (12):

$\hat{f}_{n}=\hat{f}_{n-1}+2\gamma_{n}\left(Y_{n}-\hat{f}_{n-1}(X_{n})\right)\sum_{j=1}^{d}\psi_{j}(X_{n})\psi_{j}.$
(13)

Here we have shifted to considering our estimator $\hat{f}_{n}$ as a function,
rather than thinking about $\hat{\beta}_{n}$ a vector of coefficients. This
will be important in the nonparametric setting.

###  3.2 From parametric SGD to nonparametric SGD

In this subsection we discuss the intuition in moving from SGD in a finite
dimensional parametric space to an infinite dimensional space.

We assume $f_{\rho}\in W(s,Q,\\{\psi_{j}\\})\subset L^{2}_{\nu}$. Since
$\psi_{j}$ is a complete basis of $L^{2}_{\nu}$, we can always find an
expansion of $f_{\rho}$ w.r.t. $\\{\psi_{j}\\}$:

$f=\sum_{j=1}^{\infty}\theta_{j}\psi_{j}.$ (14)

In Subsection 3.1, we already discussed the SGD updating rule for a
$d$-dimensional model $f(X)=\sum_{j=1}^{d}\beta^{(j)}\psi_{j}(X)$. In the
nonparametric scenario, the number of basis function is increased from $d$ to
infinity: This causes problems if care is not taken.

One might naturally consider applying a direct analog to the finite-
dimensional SGD rule (13) here (we omit the constant 2):

$\hat{f}_{n}=\hat{f}_{n-1}+\gamma_{n}\left(Y_{n}-\hat{f}_{n-1}(X_{n})\right)\sum_{j=1}^{\infty}\psi_{j}(X_{n})\psi_{j}.$
(15)

Unfortunately we run into a severe problem: The series
$\sum_{j=1}^{\infty}\psi_{j}(X_{n})\psi_{j}$ does not converge even if all
$\psi_{j}$ are bounded (it is direct to check when $X_{n}=0$ and $\psi_{j}$
are trigonometric functions). However, as we assume $f_{\rho}\in W(s)$, we
know that those higher order components, $\psi_{j}$, $j\gg 1$ should have very
small coefficients. Thus, one natural solution is to use a different step size
per component, that decreases as $j$ increases. By doing “less fitting” for
larger $j$, we can stabilize our update (smaller variance), and yet might
still appropriately fit the overall regression function. In particular one
might modify (15) to

$\hat{f}_{n}=\hat{f}_{n-1}+\gamma_{n}\left(Y_{n}-\hat{f}_{n-1}(X_{n})\right)\sum_{j=1}^{\infty}t_{j}\psi_{j}(X_{n})\psi_{j},$
(16)

where the component-specific (or dimension-specific) learning rate $t_{j}>0$
are monotonically decreasing with $j$. For $t_{j}$ decreasing fast enough and
uniformly bounded $\psi_{j}$, the function series
$\sum_{j=1}^{\infty}t_{j}\psi_{j}(X_{n})\psi_{j}$ is absolutely convergent.
Now (16) becomes a sensible nonparametric SGD updating rule when the
hypothesis space is a Sobolev ellipsoid. In addition, sometimes
$\sum_{j=1}^{\infty}t_{j}\psi_{j}(X_{n})\psi_{j}$ actually has a simply
characterized closed form (in particular, for many RKHS). In such cases, (16)
results in a relatively straightforward algorithm. More specifically, one can
show that when $t_{j}=j^{-2s}$ and $\gamma_{n}=\Theta(n^{-\frac{1}{2s+1}})$,
the average

$\bar{f}_{n}:=\frac{1}{n}\sum_{i=1}^{n}\hat{f}_{i}$ (17)

is a rate-optimal estimator of $f_{\rho}\in W(s)$. This was recently proposed
(though motivated quite differently) in the context of RKHS hypothesis spaces
[16]. The authors there engage directly with the _kernel function_ for the
RKHS (though their updating rule is equivalent to eq (16)). This will be
discussed in more detail in Section 4. Our work engages and extends these
ideas (in combination with sieve estimation) to form a statistically rate-
optimal online estimator with greatly reduced computational and memory
complexity.

##  4 Related work

Nonparametric online learning is a relatively new area. A few remarkable
functional stochastic approximation algorithms have been proposed in the last
two decades [9, 16, 40, 57, 65]. The key ideas in that body of work are
intimately related to those mentioned in Section 3.2, however, they engage
those ideas from a different direction: They assume that the hypothesis
function space $\mathcal{F}$ is an RKHS, and then leverage the kernel in that
space. In particular, the RKHS structure makes it possible to take the
gradient of the evaluation functional $L_{x}(f):=f(x)$, with respect to the
RKHS inner product $\langle\cdot,\cdot\rangle_{K}$, i.e.

$L_{x}(f+\epsilon g)=f(x)+\epsilon g(x)=L_{x}(f)+\epsilon\langle
g,K_{x}\rangle_{K}.$ (18)

Thus, $K_{x}(\cdot):=K(x,\cdot)\in\mathcal{F}$ is the gradient of functional
$L_{x}$ at $f$. However, one cannot do this in the general $L^{2}_{\rho_{X}}$
space where the evaluation functional is no longer a bounded operator.

Thus when $\mathcal{F}$ is an RKHS associated with kernel $K$, there is a
simple nonparametric SGD updating rule for minimizing $E[(Y-f(X))^{2}]$ over
$\mathcal{F}$:

$\displaystyle\hat{f}_{0}$ $\displaystyle=0$ (19) $\displaystyle\hat{f}_{n}$
$\displaystyle=\hat{f}_{n-1}+\gamma_{n}\left(Y_{n}-\hat{f}_{n-1}(X_{n})\right)K(X_{n},\cdot)$

Here, because the gradient is taken with respect to the RKHS inner product, we
do not have the issue encountered in (15) where our series representation of
the “gradient” actually did not converge. In fact, by working with the RKHS
inner-product, we implicitly carry out the proposal of Section 3.2 and
decrease the component-specific learning rate of higher order terms. More
specifically, we usually have the Mercer expansion of the kernel function:

$K(x,z)=\sum_{j=1}^{\infty}t_{j}\psi_{j}(x)\psi_{j}(z),$ (20)

with respect to an orthonormal basis $\\{\psi_{j}\\}$ of $L^{2}_{\nu}$. For
many common RKHS, we have $t_{j}=\Theta(j^{-u})$ for some $u>1$ [20, Appendix
A]. Thus, (19) corresponds precisely to the previously discussed update (16).
Most popular RKHS have a kernel $K(x,z)$ with a closed form representation,
and thus, rather than having to store an infinite number of coefficients,
after $n$ steps the estimate from (19) would take the form of a weighted
linear combination of $n$ kernel functions [16]:

$\hat{f}_{n}=\sum_{i=1}^{n}b_{i}K(X_{i},\cdot).$ (21)

Although such estimators (with one more Polyak averaging step (17)) have been
shown to give rate-optimal MSE [16], updating them with a new observation
$(X_{n+1},Y_{n+1})$ usually involves evaluating $n$ kernel functions at
$X_{n+1}$, with computational expense of order $\Theta(n)$. This is in
contrast with the constant update cost of $\Theta(d)$ in parametric SGD, where
$d$ is the dimension of the parameter. Thus, the time expense of nonparametric
kernel SGD will accumulate at order $\Theta(n^{2})$. Also, one is required to
store the $n$ feature-values $\\{X_{i}\\}_{i=1}^{n}$ to evaluate the estimator
which results in $\Theta(n)$ space expense. This relatively large time and
space complexity indicates that those kernel-based SGD estimators are not
ideal as methods that are nominally designed to deal with large data sets.

There has been several works in the literature aiming at improving the
computational aspect of kernel SGD methods [52, 37, 33]. These methods select
a subset of the $n$ kernel functions centered at the feature vectors and use
them as basis functions to construct estimators (which is also related to
Nyström projection). These works emphasize the application aspect of the
proposed methods. Either the statistical performance or the computational
expense is guaranteed to be optimal. Also, the theoretical analysis in these
works typically require the noise variable to have extremely light tails.

There has also been recent work [9, 40] aimed at improving kernel SGD
algorithms by leveraging approximate second order information (SGD only uses
the first order information). The estimator in [40] is shown to give rate-
optimal MSE and have better (theoretical) computational efficiency than the
vanilla kernel SGD mentioned above. However, these algorithms are usually
dramatically more complicated and have a myriad of hyper parameters that need
to be tuned.

There is another branch of research also called ”online nonparametric
regression” that engages with a different setting [22, 47]. They do not aim to
minimize the (population) generalization error directly. Their definition of
regret is based on comparing a running average of prediction error and the
empirical risk minimizer’s training error. While this is an interesting area
of research, and might be used to engage with population generalization error,
it is less directly applicable as training error is only useful as a means to
getting good generalization error.

##  5 Online Learning and the Projection Estimator: Sieve-SGD

In this section, we combine ideas from the projection estimator (in the batch
learning setting), and stochastic gradient descent to develop an estimator
that is suitable for online nonparametric regression. The estimator we will
propose achieves the minimax rate for MSE over a Sobolev ellipsoid, and is
much more computationally efficient than standard kernel SGD methods.

As a reminder, the kernel SGD estimator based on (19) has minimax rate optimal
MSE. When $\sum_{j=1}^{\infty}t_{j}\psi_{i}(s)\psi_{j}(t)$ has an available
closed form, it requires $\Theta(n)$ memory and $\Theta(n^{2})$ computation
for sequentially processing $n$ observations. We aim to improve over this and
furthermore to propose an effective estimator in cases where
$\sum_{j=1}^{\infty}t_{j}\psi_{i}(s)\psi_{j}(t)$ has no closed form.

Motivated by the projection estimator, we opt to use truncated series in the
updating rule, modifying (16) (or equivalently (19)) to get

$\hat{f}_{n}=\hat{f}_{n-1}+\gamma_{n}\left(Y_{n}-\hat{f}_{n-1}(X_{n})\right)\sum_{j=1}^{J_{n}}t_{j}\psi_{j}(X_{n})\psi_{j}$
(22)

Here $J_{n}$ is an increasing sequence of integers that grows as we collect
more observations. When $J_{n}$ is larger, the updating rule (22) is closer to
our original form (16); however, a smaller $J_{n}$ is desirable because it
results in a lower computational expense. Part of our task is identifying a
“minimal” $J_{n}$ that still maintains favorable statistical properties.

Unlike for classical projection estimators, when $t_{j}$ is properly selected,
this truncation level, $J_{n}$ (so long as it is suitably large) does not
impact the bias-variance tradeoff (to first order). By using a truncated,
rather than infinite, series, we slightly reduce variance and add some minor
higher order bias. For $J_{n}$ sufficiently large, the first order terms for
bias and variance are determined by the sequence $\\{t_{j}\\}$ (not the
truncation level). This is akin to using a truncated basis for penalized
regression in the batch learning setting. For example, in [27] and [64,
Section 5.2], the authors propose to estimate $f_{\rho}$ by solving a
penalized regression spline problem, where they use a reduced spline basis for
improved computation (rather than including a knot at every point). The
bias/variance trade-off there is controlled via the penalty: They are careful
to include enough basis elements so that the use of a reduced basis only
contributes a second order term to the bias. We will also show that, when
$J_{n}$ is properly selected, so long as the component-specific learning rate
$t_{j}$ is not too large (controlling the variance in the dynamic of SGD) or
too small (controlling the bias term), Sieve-SGD can always achieve near
optimal performance with very low computational expense. In this setting it is
the truncation level, rather than $t_{j}$, balancing the trade-off between
bias and variance, which is analogous to the batch projection estimator.

We will next give details of our proposal. For this proposal we are assuming
that $f_{\rho}\in W(s,Q,\\{\psi_{j}\\})\subset L_{\nu}^{2}$, and that $s$ is
known. Based on this, we choose our component-specific step-sizes as
$t_{j}=j^{-2\omega}$ (for some $1/2<\omega\leq s$). We also define

$K_{x,J_{n}}(\cdot)=K^{\omega}_{x,J_{n}}(\cdot)=\sum_{j=1}^{J_{n}}j^{-2\omega}\psi_{j}(x)\psi_{j}(\cdot).$
(23)

In addition to simplifying exposition, this notation relates our method to
(20). The function $K_{x,J_{n}}(\cdot)$ can be seen as a truncated
approximation of the kernel function
$K(x,\cdot)=K^{\omega}(x,\cdot)=\sum_{j=1}^{\infty}j^{-2\omega}\psi_{j}(x)\psi_{j}(\cdot)$
that drops all the $\psi_{j}$ with index $j>J_{n}$.

###  5.1 Sieve Stochastic Gradient Descent

We now explicitly give our Sieve Stochastic Gradient Descent algorithm (Sieve-
SGD) for estimation of $f_{\rho}$ in a Sobolev ellipsoid
$W(s,Q,\\{\psi_{j}\\})$.

Let $J_{n}=\lfloor n^{\alpha}\rfloor$ for some specified $\alpha>0$ and
$\omega\in(\frac{1}{2},s]$. The parameter $\alpha$ is usually taken between
$\frac{1}{2s+1}$ and 1. We use $\gamma_{i}$ to denote the step size (learning
rate) of the $i$-th update and typically choose
$\gamma_{i}=\Theta(i^{-\frac{1}{2s+1}})$.

Proposed Algorithm: Sieve Stochastic Gradient Descent (Sieve-SGD)

Set $\alpha,\omega>0$, step size $\\{\gamma_{i}\\}$ and basis functions
$\\{\psi_{j}\\}$. Initialize $\bar{f}_{0}=\hat{f}_{0}=0$.

For $i=1,2,...$ :

  1. 1.

Calculate $J_{i}=\lfloor i^{\alpha}\rfloor$, collect data pair
$(X_{i},Y_{i})$.

  2. 2.

Update $\hat{f}_{i}$:

$\displaystyle\hat{f}_{i}$
$\displaystyle=\hat{f}_{i-1}+\gamma_{i}\left(Y_{i}-\hat{f}_{i-1}(X_{i})\right)\sum_{j=1}^{J_{i}}j^{-2\omega}\psi_{j}(X_{i})\psi_{j}$
(24)
$\displaystyle=\hat{f}_{i-1}+\gamma_{i}\left(Y_{i}-\hat{f}_{i-1}(X_{i})\right)K_{X_{i},J_{i}}$

  3. 3.

Polyak averaging: Update $\bar{f}_{i}$ by

$\displaystyle\bar{f}_{i}$
$\displaystyle=\frac{1}{i+1}\sum_{k=0}^{i}\hat{f}_{k}$ (25)
$\displaystyle\left(=\frac{i}{i+1}\bar{f}_{i-1}+\frac{1}{i+1}\hat{f}_{i}\right)$

We refer to the function $\bar{f}_{i}$ as the Sieve-SGD estimate of
$f_{\rho}$. We will later show that $\bar{f}_{i}$ has rate-optimal MSE for
estimating any $f_{\rho}\in W(s)$. Here we use the language of “updating a
function”, but in practice one would update the coefficient vector
corresponding to the functions $\\{\psi_{j}\\}_{j=1}^{J_{n}}$. In Appendix A
we attach a presentation of the algorithm that works directly with the
coefficients. This estimator is quite simple, though it does require apriori
selection/knowledge of $\\{\psi_{j}\\}$ and $s$ (which can be done using a
left-out validation set in practice). Unfortunately showing its favorable
statistical properties (in Section 6) is somewhat more complex!

###  5.2 Computational expense

After examining the updating rule above, one can see that $\hat{f}_{i}$ has
the form:

$\hat{f}_{i}(x)=\sum_{j=1}^{J_{i}}b_{j}\psi_{j}(x)$ (26)

This requires storing the coefficients $\\{b_{j}\\}_{j=1}^{J_{i}}$ in memory.
The main computational burden of each update step is calculating
$\hat{f}_{i-1}(X_{i})$ and $K_{X_{i},J_{i}}$. Both require evaluating $J_{i}$
basis functions at $X_{i}$. Thus, the computational expense of the “Update
$\hat{f}_{i}$” step above is of order $J_{i}=\Theta(i^{\alpha})$, when we take
evaluating one basis function at one point as $O(1)$. And the total expense of
processing $n$ samples is of order $\Theta\left(n^{1+\alpha}\right)$. The
space expense is of the same order $\Theta(i^{\alpha})$: We need only store
coefficients of $J_{i}$ basis functions. In Section 6.4 we will show that,
under mild conditions, this memory complexity is near optimal among all
estimators with rate-optimal MSE.

This compares favorably with standard kernel SGD (21) which uses $i$ basis
functions at step $i$: Our estimator uses fewer when $\alpha<1$; as we will
show later, $\alpha$ can be taken as small as $\frac{1}{2s+1}$ which is a
substantial improvement. In practice, the parameter $\alpha$ can either be
selected based on our assumptions about $s$ (belief on the smoothness of
$f_{\rho}$) or heuristically tuned for empirical performance.

###  5.3 General Convex loss

Although the main focus of this paper is regression with squared-error loss,
our algorithm has a straightforward extension to general convex loss. Suppose
we want to minimize the population loss

$E\left[\ell(Y,f(X))\right]$ (27)

over all functions $f\in W(s,Q,\\{\psi_{j}\\})$ and the loss function
$\ell(Y,\cdot)$ is convex for each $Y$. In this case, we need only modify step
$2$ of the Sieve-SGD estimator in Section 5.1. Given loss $\ell(\cdot,\cdot)$,
the updating rule for $\hat{f}_{i}$ takes the general form:

  1. 2’)

Update $\hat{f}_{i}$:

$\hat{f}_{i}=\hat{f}_{i-1}+\left.\gamma_{i}\frac{\partial}{\partial
v}\ell\left(u,v\right)\right|_{\left(Y_{i},\hat{f}_{i-1}(X_{i})\right)}K_{X_{i},J_{i}}$
(28)

For example, with $Y=\\{1,-1\\}$ considering nonparametric logistic
regression, the loss function one would use is
$\ell(Y,f(X))=\log(1+\exp(-Yf(X)))$. In this case, we have

$\left.\frac{\partial}{\partial
v}\ell\left(u,v\right)\right|_{\left(Y_{i},\hat{f}_{i-1}(X_{i})\right)}=\left(1+\exp(Y_{i}\hat{f}_{i-1}(X_{i}))\right)^{-1}Y_{i}\in\mathbb{R}$
(29)

Theoretical guarantees for Sieve-SGD using general convex loss are beyond the
scope of this paper. However, in Section 7 we provide simulated experiments
that show the empirical performance of Sieve-SGD for nonparametric logistic
regression. These empirical results intimate that perhaps similar theoretical
guarantees to those shown for squared-error-loss hold in a more general
setting.

###  5.4 Choice of Basis Functions & Multivariate Problems

In practice, there are many available choices of univariate $\psi_{j}$ that in
general lead to interesting (Sobolev-type) hypothesis spaces. For example,

$\psi_{1}(x)=1,\quad\psi_{j}=\sqrt{2}\cos((j-1)\pi x),\text{ for }j\geq 2.$
(30)

This set of basis functions are the “eigenfunctions” of Sobolev spaces over
$[0,1]$ (Appendix A.2 in [45]), which means they are orthogonal w.r.t to the
Lebesgue inner product and the Sobolev inner product simultaneously. The
corresponding Sobolev ellipsoid does not impose periodicity assumptions of
$f_{\rho}$ and is very convenient to use in practice. Among many other
choices, we can also use algebraic polynomials, or a combination of algebraic
polynomial and (periodic) Fourier basis [19].

In most applications, the covariate $X_{i}$’s take value in $\mathbb{R}^{p}$
where $p>1$. In some situations, there are some “canonical” choices of basis
function $\psi(x):\mathbb{R}^{p}\rightarrow\mathbb{R}$ that people might use
for identifying their (multivariate) Sobolev ellipsoid. For example, when
considering estimating a function on a sphere $\mathbb{S}^{2}$, $\psi_{j}$
could be taken as the orthonormal spherical harmonics ([31], [42]).

In many situations, the basis function $\psi_{j}$ can conveniently be taken as
a tensor product of a one-dimensional complete basis, and Sieve-SGD can be
directly applied in this multivariate setting. If we are using a univariate
Sobolev ellipsoid to represent a ball in an RKHS, then the ellipsoid defined
by the tensor product basis will correspond to a ball in the RKHS spanned by
the tensor product kernel (though care needs to be taken with the ordering of
the basis vectors). Some technical details and numerical examples on this can
be found in Appendix B and the reference therein. In all of these cases, our
theoretical results will hold (so long as the function $f_{\rho}$ belongs to
the specified space).

A common alternative approach in multivariate problems is to impose some
additional structure on the hypothesis space to make estimation more
tractable. This is particularly true when the feature dimension $p$ is large.
One popular model is the nonparametric additive model [54, 29, 67], which is
thought to effectively balance model flexibility and interpretability. For
$x\in\mathbb{R}^{p}$, we might consider assuming/imposing an additive
structure on the regression function:

$f_{\rho}(x)=\sum_{k=1}^{p}f_{\rho,k}\left(x^{(k)}\right)$ (31)

where each of the component functions $f_{\rho,k}$ belong to a Sobolev
ellipsoid $W_{k}(s_{k},Q_{k},\\{\psi_{jk}\\})$. For ease of exposition, in
(31), we assume $E[Y]=0$ to avoid the need for a common intercept term. For a
more complete version with common intercept, see Appendix B. For a fixed
dimension $p$, when all $W_{k}=W^{*}$ (for some Sobolev ellipsoid $W^{*}$),
the minimax rate for estimating such an additive model is identical (up to a
multiplicative constant $p$) to the minimax rate in the analogous one-
dimension nonparametric regression problem with the same hypothesis space
$W^{*}$ [48, 54]. For the additive model (31), the updating rule (24) of
Sieve-SGD could be replaced by:

$\hat{f}_{i}=\hat{f}_{i-1}+\gamma_{i}\left(Y_{i}-\sum_{k=1}^{p}\hat{f}_{i-1,k}\left(X_{i}^{(k)}\right)\right)\sum_{k=1}^{p}\sum_{j=1}^{J_{ik}}j^{-2\omega_{k}}\psi_{jk}\left(X_{i}^{(k)}\right)\psi_{jk}$
(32)

here $J_{ik}$ is the truncation level of $k$-th dimension when the sample size
$=i$ and $\hat{f}_{i-1,k}$ is the estimate of $f_{\rho,k}$. Most of the theory
that we develop in Section 6 could apply here.

##  6 Generalization Guarantees of Sieve-SGD

In this section, we show Sieve-SGD achieves the minimax rate for nonparametric
estimation in Sobolev ellipsoids under mild assumptions. We also show that
Sieve-SGD has near minimal memory complexity among all estimators that are
minimax rate-optimal for estimation in a Sobolev ellipsoid. The conditions on
the hyperparameters can be used as theoretical guidance when applying Sieve-
SGD to real data problems.

###  6.1 Model Assumptions

We begin by listing the conditions we will require in our proof. They reflect
different aspects of the problem: independent observations (A1), distribution
of $X$ (A2), the hypothesis space assumed for $f_{\rho}$ (A3) and tail
behaviour of the noise (A4). These conditions ensure the MSE rate-optimality
of Sieve-SGD.

  1. A1

(i.i.d. data) The data points
$(X_{n},Y_{n})_{n\in\mathbb{N}}\in\mathcal{X}\times\mathbb{R}$ are
independently, identically sampled from a distribution $\rho(X,Y)$.

  2. A2

(feature distribution) Let $\nu$ be a user-specified measure that is strictly
positive on $\mathcal{X}$. Assume the distribution of feature $X$, $\rho_{X}$,
is absolutely continuous w.r.t. $\nu$. Let $p_{X}=d\rho_{X}/d\nu$ denote its
Radon–Nikodym derivative. We assume for some $u,\ell$ such that
$0<\ell<u<\infty$:

$\ell\leq p_{X}(x)\leq u\quad\text{for all }x\in\mathcal{X}$

  3. A3

(Sobolev ellipsoid) Let $\\{\psi_{j}\\}_{j=1}^{\infty}$ be a set of uniformly
bounded ($\|\psi_{j}\|_{\infty}\leq M$), continuous, orthonormal basis of
$L^{2}_{\nu}$. We assume the regression function $f_{\rho}$ falls in a Sobolev
ellipsoid, with basis functions given by $\\{\psi_{j}\\}$, i.e. for some
$s>1,Q<\infty$,

$f_{\rho}\in W(s,Q,\\{\psi_{j}\\})$ (33)

  4. A4

(noise) One of the following two assumptions is satisfied by the noise
variable $\epsilon=Y-f_{\rho}(X)$:

     * •

$\epsilon$ is bounded by some $C_{\epsilon}$ almost surely.

     * •

$\epsilon$ is independent of the features, $X$, and has a finite second moment
$E[\epsilon^{2}]=C_{\epsilon}^{2}$.

Note: In assumption A3, we do not require $\psi_{j}$ to be orthonormal w.r.t.
$\rho_{X}$ (and it is in general not true), but only require them to be
orthonormal w.r.t. the known measure $\nu$. In many cases $\nu$ might be taken
to be Lesbesgue (or uniform) measure over a domain containing $\mathcal{X}$,
as this is the canonical measure under which function spaces such as Sobolev
spaces and Besov spaces are defined. As long as the density function $p_{X}$
satisfies A2, using the non-orthonormal (w.r.t. $\rho_{X}$) basis functions,
$\psi_{j}$, does not prevent Sieve-SGD from having rate-optimal MSE. Also, the
lower bound requirement of $p_{X}$ in A2 may be due to artifacts in our proof.
In reality, especially when the dimension of feature $X$ is higher, such an
requirement is hard to be satisfied. According to our simulation results,
Sieve-SGD still achieves the minimax rate even when $\rho_{X}$ has a strictly
smaller support than $\nu$. As compared with other work in nonparametric
online learning [16, 57, 65], our assumptions are more direct. We discuss this
in detail later in this section.

###  6.2 Rate optimality when $t_{j}=j^{-2s}$

In this section, we present the rate-optimality results of Sieve-SGD when we
choosing the component-specific learning rate to be $t_{j}=j^{-2s}$ (or
$\omega=s$ in (24)). In this setting, our theoretical analysis treats Sieve-
SGD as a truncated-version (in the basis expansion domain) of a “correct”
kernel SGD procedure (we will discuss the “incorrect” version very soon in
Section 6.3). Here is the main result in this setting:

######  Theorem 6.1.

Assume A1-A4. Set the component-specific learning rate as $t_{j}=j^{-2s}$.
Also set the overall learning rate to be
$\gamma_{n}=\gamma_{0}n^{-\frac{1}{2s+1}}$ with
$\gamma_{0}\leq(2M^{2}\zeta(2s))^{-1}$, where
$\zeta(k)~{}=~{}\sum_{i=1}^{\infty}i^{-k}$. Choose the number of basis
functions to be $J_{n}\geq n^{\alpha}\log^{2}n\vee 1$ for an arbitrary
$\alpha\geq\frac{1}{2s+1}$.

Then the MSE of Sieve-SGD (25) converges at the following rate

$E\|\bar{f}_{n}-f_{\rho}\|^{2}_{L^{2}_{\rho_{X}}}=O\left(n^{-\frac{2s}{2s+1}}\right).$
(34)

This implies that Sieve-SGD is a minimax rate-optimal estimator of $f_{\rho}$
over $W(s,Q,\\{\psi_{j}\\})$.

We now discuss our assumptions and results in more detail, and relate them to
what is currently in the literature.

Note 1: In the analysis of many reproducing kernel methods for nonparametric
estimation [16, 57, 66], the spectrum of the _covariance operator_ plays an
important role in controlling the statistical behavior of estimators. It is
conventional in the community to make assumptions associated with this
spectrum, which we find less natural than our related assumptions A2 and A3.
The covariance operator is an analog of the covariance matrix in infinite
dimensional spaces. For our problem setting, one of the natural covariance
operator $T_{X}$ is defined as:

$\displaystyle T_{X}:L_{\rho_{X}}^{2}$ $\displaystyle\rightarrow
L_{\rho_{X}}^{2}$ (35) $\displaystyle g$
$\displaystyle\mapsto\int_{\mathcal{X}}g(\tau)\left(\sum_{j=1}^{\infty}j^{-2s}\psi_{j}(\tau)\psi_{j}\right)d\rho_{X}(\tau).$

A direct analysis of the spectrum of $T_{X}$ is hard. However, there is a
simpler operator that we have in hand which we can relate $T_{X}$ to:

$\displaystyle T_{\nu}:L^{2}_{\nu}$ $\displaystyle\rightarrow L^{2}_{\nu}$
(36) $\displaystyle g$
$\displaystyle\mapsto\int_{\mathcal{X}}g(\tau)\left(\sum_{j=1}^{\infty}j^{-2s}\psi_{j}(\tau)\psi_{j}\right)d\nu(\tau).$

We know the eigensystem of $T_{\nu}$: It is exactly $(j^{-2s},\psi_{j})$
(eigenvalue, eigenfunction). It is direct to check because $\\{\psi_{j}\\}$’s
are orthonormal w.r.t. $\nu$, so
$\int\psi_{j}(\tau)\sum_{j=1}^{\infty}j^{-2s}\psi_{j}(\tau)\psi_{j}d\nu(\tau)=j^{-2s}\psi_{j}$.
As an additional contribution, our work shows that with the simple assumption
A2 & A3, we can get knowledge about $T_{X}$’s eigenvalues from those of
$T_{\nu}$.

######  Lemma 6.2.

Given assumptions A2, A3, the $j$-th eigenvalue, $\lambda_{j}$, (sorted in a
decreasing order) of $T_{X}$ satisfies $\lambda_{j}=\Theta(j^{-2s})$.

Moreover, the upper bound of the density in A2 ensures the upper bound in
Lemma 6.2 ($\lambda_{j}=O(j^{-2s})$), and the lower bound of the density
ensures the other half of the result. The proof of the above Lemma uses the
underlying connection between the eigenvalues of an operator and its metric
entropy. For rigorous definitions and proof of Lemma 6.2, see Appendix C.

Although the exact result of Lemma 6.2 is not used in the proof of Theorem 6.1
(or Theorem 6.3). We still present it here since it may be of interest itself
and the stated results is less technical and easier to comprehend. The proof
of more technical version (Lemma C.14) follows very closely to that of Lemma
6.2. In such a more general version, we investigate the spectrum of covariance
operators of form:$T_{X,J_{n}}^{\omega}(f)=\int
f(\tau)\left(\sum_{j=1}^{J_{n}}j^{-2\omega}\psi_{j}(x)\psi_{j}\right)d\rho_{X}(\tau).$

To prove Theorem 6.1, we need to engage with a series of RKHSs with kernels
given by

$\displaystyle K_{J_{n}}:\mathcal{X}\times\mathcal{X}$
$\displaystyle\rightarrow\mathbb{R}$ (37) $\displaystyle(s,t)$
$\displaystyle\mapsto\sum_{j=1}^{J_{n}}j^{-2s}\psi_{j}(s)\psi_{j}(t):=K_{J_{n}}(s,t).$

While we discuss our work in the context of Sobolev ellipsoids, there is an
equivalent formulation directly in RKHS. See Appendix C for more discussion.
Although an explicit form for $K_{J_{n}}$ is not in general necessary or
accessible for Sieve-SGD, the existence (i.e. the absolute convergence of the
infinite sum) of $K_{J_{n}}$ is a direct consequence of A3. This is enough for
theoretical analysis. For kernel SGD methods, a fixed kernel (with
$J_{n}=\infty$) is used and there is only one relevant RKHS. This means, on
average, kernel SGD is applying the same procedure each iteration; but for
Sieve-SGD, we need to engage with a series of increasing RKHSs (on average,
Sieve-SGD may not be doing the same thing between iterations). As a side
contribution, we present how to handle such a more technically involved case.

Note 2: In contrast to our assumption A3, the hypothesis spaces in [16, 57,
65, 40] are described in terms of “$T_{X}$” and its eigen-decomposition. This
unfortunately obfuscates difficulties related to verifying those conditions:
In particular because $\rho_{X}$ is involved in the definition of $T_{X}$
(35), we need knowledge of (generally unknown) $\rho_{X}$ to characterize
$T_{X}$, and understand its eigenvalues and eigenfunctions.

More specifically, in the literature we mentioned above, it is often assumed
that for some $r\in[1/2,1]$ (Definition C.6):

$\|T_{X}^{-r}(f_{\rho})\|_{L^{2}_{\rho_{X}}}^{2}<\infty$ (38)

This can be related to a Sobolev ellipsoid-type condition

$\|T_{X}^{-r}(f_{\rho})\|_{L^{2}_{\rho_{X}}}^{2}=\sum_{j=1}^{\infty}\lambda_{j}^{-2r}\theta_{j}^{2}<\infty\quad\text{
where }f_{\rho}=\sum_{j=1}^{\infty}\theta_{j}\phi_{j}$ (39)

where $(\lambda_{j},\phi_{j})_{j=1}^{\infty}$ are the eigenvalue and
eigenfunctions of operator $T_{X}$, and $\phi_{j}$’s are orthonormal w.r.t.
$L^{2}_{\rho_{X}}$. Unfortunately, we cannot directly engage with
$(\lambda_{j},\phi_{j})_{j=1}^{\infty}$, since calculating them requires
knowledge of $\rho_{X}$. Thus, assumptions formulated in the language of
$T_{X}^{-r}$ are difficult to directly understand. In contrast, our
assumptions translate to analyzing the spectrum of $T_{\nu}$, which has no
dependence on $\rho_{X}$, and its spectrum can been directly calculated (as
noted above).

Note 3: For parametric SGD methods, usually a bound on the second moment of
the gradient vector is required to guarantee rate-optimal performance (both
theoretically and in practice). Formally, for optimization problem (9), it is
usually required that $E[\|\nabla\ell(\beta)\|^{2}]\leq~{}R^{2}<~{}\infty$ for
all $\beta\in\mathbb{R}^{d}$ [6, 17].

For nonparametric stochastic approximation, there is a similar regularity
requirement for the “gradient”. The assumptions A2-A3 are enough to ensure
this for Sieve-SGD. In our proof, we show that there exists a number
$R<\infty$ such that for all $x\in\mathcal{X}$ and any $J_{n}$, we have
$\|K_{x,J_{n}}\|_{K}^{2}\leq R^{2}$. This result is listed in Lemma D.1 where
$R^{2}=M^{2}\zeta(2s)$ and $\zeta(k)=\sum_{i=1}^{\infty}i^{-k}$. In Theorem
6.1, we required $\gamma_{0}$ to be smaller than $(2M^{2}\zeta(2s))^{-1}$ to
ensure our theoretical guarantees.

Note 4: For completeness, here we state the minimax-rate of our nonparametric
regression problem over a Sobolev ellipsoid:

$\liminf_{n\rightarrow\infty}\inf_{\hat{f}}\sup_{f_{\rho}\in
W(s,Q,\\{\psi_{j}\\})}E\left[n^{\frac{2s}{2s+1}}\|\hat{f}-f_{\rho}\|_{L^{2}_{\rho_{X}}}^{2}\right]\geq
C$ (40)

where the infimum ranges over all possible functions $\hat{f}$ that are
sufficiently measurable. For a derivation of this lower bound, see [62,
Chapter 15].

###  6.3 Adaptivity to $t_{j}$ for Properly Chosen $J_{n}$

In section 6.2 we presented the optimality guarantee of Sieve-SGD: when the
component-specific learning rate is properly chosen, i.e. $t_{j}=j^{-2s}$,
Sieve-SGD is statistically optimal so long as the number of basis functions
does not increase too slow, that is, $J_{n}\geq
n^{\frac{1}{2s+1}}\log^{2}(n)$. Specifically, when $J_{n}=\infty$, the Sieve-
SGD updating rule reduces to the kernel SGD updating rule (19) with kernel
$K(X_{n},\cdot)=\sum_{j=1}^{\infty}j^{-2s}\psi_{j}(X_{n})\psi_{j}(\cdot)$. So
long as we have access to the closed-form of $K(X_{n},\cdot)$, the
corresponding kernel SGD is also optimal under the same conditions. In such a
scenario, Sieve-SGD can be seen as a truncated-version of a “correct” kernel
SGD method with much better computational properties.

However, if we choose $t_{j}=j^{-2\omega},\omega\neq s$, the corresponding
kernel SGD, using kernel
$K(X_{n},\cdot)=\sum_{j=1}^{\infty}j^{-2\omega}\psi_{j}(X_{n})\psi_{j}(\cdot)$,
is no longer optimal without modifying the learning rate $\gamma_{n}$
accordingly [16]. But for Sieve-SGD, so long as the truncation level $J_{n}$
is properly selected, the statistical performance of Sieve-SGD is near optimal
for a quite wide range of choice of $\omega$.

######  Theorem 6.3.

Assume A1-A4. Set the component-specific learning rate to be
$t_{j}=j^{-2\omega}$ with $\frac{1}{2}<\omega<s$. Choose the learning rate to
be $\gamma_{n}=\gamma_{0}n^{-\frac{1}{2s+1}}$, with $\gamma_{0}\leq
M^{2}\zeta(2\omega)/2$. Choose the number of basis functions to be
$J_{n}~{}=~{}n^{\frac{1}{2s+1}}\log^{2}n\vee 1$.

Then the MSE of Sieve-SGD (25) converges at the following near optimal rate

$E\|\bar{f}_{n}-f_{\rho}\|^{2}_{L^{2}_{\rho_{X}}}=O\left(n^{-\frac{2s}{2s+1}}\log^{2}n\right)$
(41)

Note 1: The requirement of $t_{j}<j^{-1}$ is to guarantee a finite “second
moment” of the gradient, recall the Note 3 under Theorem 6.1. Once such a
minimal requirement is satisfied, the decay rate of $t_{j}$ does not influence
neither the rate of statistical guarantee, nor the computational expense of
the estimators. As we will discuss very soon in section 6.4, the choice of
$J_{n}=n^{\frac{1}{2s+1}}\log^{2}n$ in Theorem 6.3 and Theorem 6.1 would
result in algorithms that are both statistically and computationally near-
optimal, which is very rare in the literature of online nonparametric learning
and could be of interest in practice.

Note 2: The most direct form of projection estimator determines the basis
functions’ coefficients by solving a (unpenalized) least square problem (7) in
which there are no learning rates involved. It is the truncation level $J_{n}$
that determines the bias-variance trade-off and statistical performance. In
Theorem 6.3 we present a stochastic approximation version of such a result.
From a reproducing-kernel methodology perspective, Theorem 6.1 investigates
the cases when the capacity of the kernel ($\omega$) matches the source
smoothness ($s$); in Theorem 6.3 we discussed under what conditions the
mismatch between these two quantities does not affect the statistical (and
computational) properties of Sieve-SGD. We also note that the overall proof
structures of Theorem 6.1 and Theorem 6.3 are similar; the difference is, in
the proof of Theorem 6.1 we need Lemma D.4 and related technical results, but
for Theorem 6.3 we use Lemma E.1 instead.

###  6.4 Near optimal space expense

In this section we will show that Sieve-SGD is asymptotically (nearly) space-
optimal for estimating $f_{\rho}$ in a Sobolev ellipsoid under the conditions
listed in Section 6.1. We will show that, even with computer round-off error,
Sieve-SGD only needs $\Theta(n^{\frac{1}{2s+1}}\log^{3}n)$ _bits_ to achieve
the minimax rate for MSE (or off by a $\log^{2}(n)$ term when $\omega\neq s$
as stated in Theorem 6.3), and further, that there is no estimator with
$o(n^{\frac{1}{2s+1}})$ bits of space expense that can achieve the minimax-
rate for estimating $f_{\rho}\in W(s,Q)$. In our analysis we note that
computers cannot store decimals in infinite precision, and formally deal with
a modified version of our algorithm that stores coefficients in fixed
precision (that grows in $n$): This necessitates the extra $\log(n)$ term
(compared with the number of basis function needed in Theorem 6.1 and 6.3).
The modified algorithm with fixed, but growing precision still results in the
same MSE when round-off error is not considered.

We first give a more formal definition of the space expense of an estimator in
our analysis. A regression estimator can be seen as a mapping $M_{n}$ from the
data $Z_{1}^{n}=\\{(X_{i},Y_{i})\ |\ i=1,2,...,n\\}$ to a function
$\hat{f}_{n}\in\mathcal{F}$. For any such $M_{n}$ that can be engaged by a
computer, must be decomposable into an “encoder-decoder” pair $(E_{n},D_{n})$.
Here $E_{n}$ represents the “encoder” that compresses the information into
computer memory. Formally, we define $E_{n}$ to be a mapping from $Z_{1}^{n}$
to a binary sequence of length $b_{n}$. And the corresponding $D_{n}$ is the
“decoder” of the binary sequence that translates the information saved in
memory back to a mathematical object $\hat{f}_{n}$. Generally, the binary
sequence length $b_{n}$ will increase with $n$: As more information is
contained in the data, we need more memory to store an increasingly accurate
estimate of our regression function.

Given an estimator that can decomposed into a pair $(E_{n},D_{n})$, one can
see that the decomposition is not unique. There are, in fact, infinitely many
pairs that are trivially different from each other for any such estimator.
Moreover, $E_{n},D_{n}$’s can be random mappings if we allow random
algorithms: For example, random forests include additional randomness due to
bootstrapping/variable selection. In order to be more precise regarding memory
complexity constraints, we introduce the following formalization.

######  Definition 6.4 ($b_{n}$-sized estimator).

Given a sequence of integers $(b_{i})_{i\in\mathbb{N}}$, we say an estimator
$M_{n}:(\mathcal{X}\times\mathbb{R})^{n}\rightarrow\mathcal{F}$ is a
$b_{n}$-sized estimator if it satisfies the following conditions:

  1. 1.

For every $n$, there exists an encoder mapping
$E_{n}:(\mathcal{X}\times\mathbb{R})^{n}\rightarrow\\{0,1\\}^{b_{n}}$, and a
decoder mapping $D_{n}:\\{0,1\\}^{b_{n}}\rightarrow\mathcal{F}$ such that

$M_{n}=D_{n}\circ E_{n}$ (42)

  2. 2.

The decoder $D_{n}$ is a known, fixed mapping. $E_{n}$ can be either a random
or fixed mapping.

We use the sample mean as a toy example to illustrate the above definition. In
practice, the sample mean is usually a 64-sized estimator of the population
mean. Here 64 stands for the number of bits needed to represent a double-
precision floating point number. In this case the size $b_{n}=64$ does not
increase with sample size $n$. However not every real number can be
arbitrarily precisely specified by a fixed-length floating-point number, so a
careful asymptotic analysis of estimation of the mean suggests that perhaps we
should store a sample mean with growing levels of precision, i.e. $b_{n}$
would need to grow with $n$. A binary sequence of length $s$ can specify
$2^{s}$ real numbers, so to achieve the $O(n^{-1})$ statistically optimal
bound for mean estimation, a $\log(n)$-sized version of sample mean is
formally required. In practice, 64-bit precision is generally more than enough
for mean estimation. Nevertheless, in this manuscript we aim to give a more
formal and precise asymptotic analysis of our Sieve-SGD estimator.

Readers who are more familiar with computational complexity theory may find
our definition similar to a (probabilistic) Turing machine. However, in our
framework the machine does not use binary sequences on tapes as input and
output; nor do we need to identify the basic operations on the ”machine”. We
aimed to remove unnecessary complexity for readers with a more statistical
background. Discussion of Turing machines using finite length working tape can
be found in [1, Chapter 4].

To construct Sieve-SGD estimators that achieve (near) optimal MSE, we only
need to store the coefficients of the
$J_{n}=\Theta(n^{\frac{1}{2s+1}}\log^{2}n)$ basis functions. However, as in
our example with the sample mean, we need to be careful about the precision
with which we store those coefficients. We need to determine: i) how small we
require the round-off error to be in order to maintain the statistical
optimality of Sieve-SGD, and ii) how much space expense is required to achieve
such precision. In Appendix F.1 we identify how round-off error is introduced
into the system and how it decreases as more bits are used to store each
coefficient. In Corollary F.2 we show that by using $\Theta(\log n)$ bits per
coefficient, a $O(n^{\frac{1}{2s+1}}\log^{3}n)$-sized version of Sieve-SGD can
achieve the same optimal convergence rate as in the infinite precision setting
(or equivalently round-off-error-free setting).

Combining the above result with the following theorem, we can conclude that no
MSE rate-optimal estimator can require less memory by a polynomial factor than
Sieve-SGD.

######  Theorem 6.5.

Let $b_{n}$ be a sequence of integers, and
$b_{n}=o\left(n^{\frac{1}{2s+1}}\right)$. Let $\mathcal{M}(b_{n})$ be the
collection of all $b_{n}$-sized estimators, then we have

$\lim_{n\rightarrow\infty}\inf_{M_{n}\in\mathcal{M}(b_{n})}\sup_{f_{\rho}\in
W(s,Q,\\{\psi_{j}\\})}E\left[n^{\frac{2s}{2s+1}}\|M_{n}(Z_{1}^{n})-f_{\rho}\|_{L^{2}_{\rho_{X}}}^{2}\right]=\infty$
(43)

i.e. no such $b_{n}$-sized estimators can be rate-optimal.

This theorem tells us we cannot find any minimax rate-optimal
$o(n^{\frac{1}{2s+1}})$-sized estimator. Thus the best rate-optimal estimator
one can expect to find is a $\Theta(n^{\frac{1}{2s+1}})$-sized estimator:
Sieve-SGD’s space expense only misses this lower bound by a poly-logarithmic
factor.

We give the proof of the above theorem in Appendix F.2. Although here we focus
on regression in Sobolev spaces, the technique used can be applied to other
hypothesis spaces. The proof is based on the concept that metric-entropy is
the minimal number of bits needed to represent an arbitrary function from a
function space up to $\epsilon$-error, which can be traced back to [32]. Also,
following a very similar argument, one can prove that no constant-sized
estimator can be rate-optimal (or even consistent) for parametric regression
problems. We discuss this further in the Appendix F.2. We also include some
discussion of the time expense in Section 8.

##  7 Simulation study

###  7.1 Sieve-SGD for online regression

In this section, we illustrate both the statistical and computational
properties of Sieve-SGD with simulated examples. The two examples we use have
different $f_{\rho}$, $W(s,Q,\\{\psi_{j}\\})$ and $\rho_{X}$. The user-
specified measure $\nu$ is taken as the uniform distribution over $[0,1]$ in
both. We provide the details of our simulation settings in Table 1. These two
examples are designed for verifying our theoretical guarantees: The $f_{\rho}$
we use have known explicit series expansion or is constructed explicitly using
the basis function $\psi_{j}$ (to ensure the truth is hard enough to learn in
the assumed Sobolev ellipsoid). In Appendix B we provide more numerical
examples to better mimic the practical application: we engage with multi-
variate features and compare Sieve-SGD with many popular machine learning
methods.

Example 1 In this example, we examine the empirical performance of Sieve-SGD
and compare it with two other methods in batch or online nonparametric
regression: kernel ridge regression (KRR) [62] and kernel SGD [16]. We will
see that the relationship between generalization error
$E\|\bar{f}_{n}-f_{\rho}\|_{2}^{2}$ and sample size corresponds well with our
theoretical expectations presented in Theorem 6.1 (Fig 1).

The true regression function we chose for Example 1 is also used in the
analysis of kernel SGD [16]. In that paper, kernel SGD with Polyak averaging
is compared with other (kernel-based) nonparametric online estimators [57,
65], and has been shown to have similar or better performance, so we include
only kernel SGD with averaging as the reference online-estimator. We also note
that although KRR performs slightly better than online methods, its time
expense (which is of order $\Theta(n^{3})$ per update) is dramatically more
than online-estimators (kernel SGD $\Theta(n)$, Sieve-SGD $\Theta(J_{n})$, per
update).

We compare the empirical performance of Sieve-SGD under two different
distributions of $X$. In Fig 1 panel (A), $X$ has an uniform distribution over
$[0,1]$ and in panel (B) it has a distribution with a strictly smaller support
(uniform over $[0.25,0.75]$). The trigonometric basis functions we use are
orthonormal w.r.t. $\nu$, the Lebesgue measure over $[0,1]$ (panel (A)) but
not w.r.t. the one in panel (B). Although only the feature distribution in
panel (A) satisfies the distribution assumption in A2, in both cases Sieve-SGD
achieves the optimal-rate. This is a heuristic evidence indicating the lower
bound requirement in A2 may be due to some artifacts in the proof.

Table 1: Settings of simulation studies. $B_{4}(x)=x^{4}-2x^{3}+x^{2}-\frac{1}{30}$ is the 4-th Bernoulli polynomial. $\\{x\\}$ indicates the fractional part of $x$.  | Example 1 | Example 2
---|---|---
True $f_{\rho}$ | $B_{4}(x)$ | $4\sqrt{2}\sum_{j=1}^{50}(-1)^{j+1}j^{-4}\sin((2j-1)\pi x/2)$
ellipsoid para. $s$ | 2 | 3
$J_{n}$ | $n^{0.21}$ |  $n^{0.10}$ & $n^{0.15}$ & $n^{0.43}$
$t_{j}$ |  $j^{-1.02}$ & $j^{-4}$ | $j^{-6}$
$\psi_{j}(x)$ |  $\sin(2\pi\lceil j/2\rceil x)$, $j$ is even | $\sqrt{2}\sin(\frac{(2j-1)\pi x}{2})$
$\cos(2\pi\lceil j/2\rceil x)$, $j$ is odd
Kernel $K(s,t)$ | $-\frac{1}{24}B_{4}(\\{s-t\\})$ | $\min(s,t)$
Noise | Unif[$-0.02,0.02$] or Unif[$-0.2,0.2$] | Normal(0,1)
$\gamma_{0}$ | 3 | 1
Figure 1: Example 1, $\log_{10}\|\bar{f}_{n}-f_{\rho}\|_{2}^{2}$ against
$\log_{10}n$. The Black line has slope $=-4/5$, which represents the optimal-
rate. Each curve is calculated as the average of 100 repetitions. (A) $X$ is
uniformly distributed over $[0,1]$. In this setting, SNR $\sim 3$. (B) $X$ has
a distribution in which $\psi_{j}$ are not orthonormal. We present the results
with very large noise, SNR$\sim 0.03$. Due to different computational costs,
we chose different maximum $n$ for different methods.

Example 2 In this example, we consider the performance of Sieve-SGD under
different $J_{n}=\lfloor n^{\alpha}\rfloor$ (number of basis functions). The
$f_{\rho}$ we use is explicitly constructed with basis functions
$\psi_{j}(x)=\sqrt{2}\sin\left((2j-1)\pi x/2\right)$ and we tune the proposed
method based on the (correct) assumption that it belongs to Sobolev ellipsoid
$W(3,Q,\\{\psi_{j}\\})$ (see Theorem 4.1 of [30, Chapter 1] for completeness
and orthonormality of $\\{\psi_{j}\\}$).

According to Theorem 6.1, in order to guarantee statistical optimality, we
need $\alpha\geq\frac{1}{(2s+1)}\sim 0.14$. We consider several values of
$\alpha$, one below the this threshold, and two above it:

$\mathbf{0.10}<\frac{1}{2s+1}\sim 0.14<\mathbf{0.15}<\mathbf{0.43}$ (44)

As we can see from Fig 2 (A), when $\alpha=0.15\ \&\ 0.43$, Sieve-SGD is rate-
optimal as expected. When $\alpha=0.10$, we are using too few basis functions,
which results in the sub-optimal statistical performance. Such a suboptimality
is because of the bias term: there are too few basis functions used. In fact,
the parameter setting $\alpha=0.1$ is so small that there are only $3$ basis
functions used when $n=10^{5}$. To verify the above statement, we can briefly
calculate when the second and the third basis functions are added in:
$(10^{3})^{0.1}\sim 2$, this corresponds to the first acceleration of the
learning rate around $\log_{10}(n)=3$; similarly, $(10^{4.8})^{0.1}\sim 3$,
which explains the second one.

In Fig 2 (B), we show the CPU time for reference. For Sieve-SGD, the
accumulated CPU time should be on the order of $\Theta(n^{1+\alpha})$: The
larger $\alpha$, the more basis functions required, the slower the algorithm.
We also include the CPU time of kernel SGD with averaging as a benchmark,
which has a cumulative computational expense of order $\Theta(n^{2})$. The
code is written in R (4.0.4), and runs on (the CPU of) a machine with 1 Intel
Core m3 processor, 1.2 GHz, with 8 GB of RAM.

Figure 2: Example 2, effect of truncation exponents $\alpha=0.10,0.15,0.43$.
(A) Statistical performance, $\log_{10}\|\bar{f}_{n}-f_{\rho}\|_{2}^{2}$
against $\log_{10}n$. The black line has slope $=-6/7$, which represents the
optimal-rate. Each curve is calculated as the average of 100 repetitions. (B)
The accumulated CPU time to process $n$ observations. The black line is the
CPU time of kernel SGD, included for benchmark.

###  7.2 Sieve-SGD for Alternative Convex Losses

In this section, we provide the results of an experiment applying Sieve-SGD to
online nonparametric logistic regression. Although this manuscript gives no
theoretical guarantees in this setting, it is still of interest to see the
empirical performance of Sieve-SGD for general convex loss. Here, the
distribution of class labels $Y$ was generated by $Y\sim
2\operatorname{Ber}(g(X))-1$, where $(g(x))^{-1}=1+\exp(-5(1-2|x-0.5|))$; and
the distribution of $X$ was uniform over $[0,1]$. Thus, the minimizer $f^{*}$
of loss $E[\ell(Y,f(X))]=E[\log(1+\exp(-Yf(X)))]$ is $f^{*}=5(1-2|x-0.5|)$.

When we apply the Sieve-SGD estimator (28) to this problem, we assume

$f^{*}\in W\left(1,Q,\left\\{\sqrt{2}\sin\left((2j-1)\pi
x/2\right)\right\\}\right)$ (45)

We try several $\alpha=0.10,0.33,0.50$, all with $\gamma_{0}=6$. As we can see
from Fig 3, the regret $E[\ell(\bar{f}_{n})-\ell(f^{*})]$ converges to zero at
an apparent rate of $n^{-2/3}$ when $\alpha=0.33,0.50$ (which would agree with
our result for squared error loss). When the number of basis functions
increases too slowly (here is $\alpha=0.10$), the regret decreases slowly
after $\sim 10$ observations (for similar reason of overflowing bias term as
we noted in section 7.1).

Figure 3: Example 3, empirical performance of Sieve-SGD in nonparametric
logistic regression problem. Plot $\log_{10}(l(\bar{f}_{n})-l(f^{*}))$ against
$\log_{10}n$. The Black line has slope $=-2/3$. Each curve is calculated as
the average of 100 repetitions.

##  8 Discussion

In this paper, we considered online nonparametric regression in a Sobolev
ellipsoid. We proposed the _Sieve Stochastic Gradient Descent estimator
(Sieve-SGD)_ , an online estimator inspired by both a) the nonparametric
projection estimator, which is a special realization of general sieve
estimators; and b) estimators constructed using stochastic gradient descent
algorithms. By using an increasing number of basis functions, Sieve-SGD has
rate-optimal estimation error and is computationally very efficient.

For online learning problems with general convex losses, the optimal
estimation rate depends on both the hypothesis space and loss function (e.g.
whether it is Lipschitz or strongly convex). In this paper we did not
establish theoretical guarantees for Sieve-SGD when applied to general convex
loss, however, we gave some empirical evidence that it can perform well there.
We believe our proof techniques might be extended beyond squared-error loss,
perhaps using ideas in [3, 10, 40, 41].

We’ve seen a rich collection of work in the past decade targeting the
optimality of estimators under computational (especially time expense)
constraints. A lot of those results are established in the context of sparse
PCA and related sparse-low-rank matrix problems, e.g. [8, 23, 24, 38, 63, 69].
The main focus of these work is usually comparing the statistical performance
of the best polynomial-time algorithm with that of the ”optimal” algorithm
without any computational restrictions. By relating their statistical problem
with a known NP problem [1], they can usually show the sub-optimality of
polynomial-time algorithms under the famous conjecture $P\neq NP$. However,
for the nonparametric regression problem in this paper, there is a polynomial-
time estimator that can achieve the global minimax-rate. It is of theoretical
interest to know if there are any statistically rate-optimal online estimators
that require less than $\Theta(n^{1+\frac{1}{2s+1}})$ time-expense: We
hypothesize that there are not.

## Acknowledgments

N.S and T.Z. were both supported by NIH grant R01HL137808.

##  Appendix A Algorithm of Sieve-SGD, Numerical Version

In the main text, section 5.1, we presented a functional form of the proposed
Sieve-SGD algorithm. To facilitate the comprehension of it, we also attach an
equivalent numerical version of it.

Proposed Algorithm: Sieve Stochastic Gradient Descent (Sieve-SGD)

Set $\alpha,\omega>0$, step size $\\{\gamma_{i}\\}$ and basis functions
$\\{\psi_{j}\\}$.

Initialize $\bar{\beta}_{j},\hat{\beta}_{j}=0$ for all $j\in\mathbb{N}^{+}$.

For $i=1,2,...$ :

  1. 1.

Calculate $J_{i}=\lfloor i^{\alpha}\rfloor$, collect data pair
$(X_{i},Y_{i})$.

  2. 2.

Update $\hat{f}_{i}$:

$\text{res}_{i}\leftarrow
Y_{i}-\sum_{j=1}^{J_{i-1}}\hat{\beta}_{j}\psi_{j}(X_{i})$ (46)

For $j=1,...,J_{i}$:

$\hat{\beta}_{j}\leftarrow\hat{\beta}_{j}+\gamma_{i}\text{res}_{i}\left(j^{-2\omega}\psi_{j}(X_{i})\right)$
(47)

  3. 3.

Update $\bar{f}_{i}$

For $j=1,...,J_{i}$:

$\bar{\beta}_{j}\leftarrow\frac{i}{i+1}\bar{\beta}_{j}+\frac{1}{i+1}\hat{\beta}_{j}$
(48)

##  Appendix B Multivariate Regression Problems

In this section, we will give additional discussion of the technical details
for multivariate regression using the Sieve-SGD estimator.

###  B.1 Hyperbolic cross and Sieve-SGD

In main text Section 5.4, we discussed using a tensor product basis to
approach multivariate problems. Here we go into more details and more
technical discussion. Given an univariate orthonormal basis
$\\{\psi_{j},j\in\mathbb{N}^{+}\\}$, the set of tensor product functions
$\\{\psi_{\mathbf{j}}(\mathbf{x})=\prod_{k=1}^{p}\psi_{\mathbf{j}^{(k)}}(\mathbf{x}^{(k)}),\mathbf{j}\in\left(\mathbb{N}^{+}\right)^{p}\\}$
is also an orthonormal basis ($\mathbf{v}^{(k)}$ is the $k$-th component of
$\mathbf{v}\in\mathbb{R}^{p}$). However, there are more choices of the order
in which we include basis functions when estimating an unknown regression
function. We propose using the index product $\prod_{k=1}^{p}\mathbf{j}^{(k)}$
to determine such an ordering. That is, basis functions with smaller index
product will be used earlier when constructing Sieve-SGD, such a choice is
called hyperbolic cross in the literature [18, 50]. Before we discuss the
intuition of such an ordering (which has been established in the literature),
we present some numerical examples of applying Sieve-SGD in multi-variate
regression problems.

We consider two dimension settings of the feature variable $\mathbf{X}$,
$p=2\text{ and }10$. The feature vector is defined as:
$\mathbf{X}^{(1)}=U_{1}$, $\mathbf{X}^{(k)}=(U_{k}-U_{k-1}+1)/2$ for
$k=2,...,p$. Here $U_{k}$ are independent $\text{Unif}[0,1]$ variables. The
true regression function is defined as

$f_{\rho}(\mathbf{x})=\sum_{k=1}^{p}\sum_{l=k}^{p}(0.5-|\mathbf{x}^{(k)}-0.5|)(0.5-|\mathbf{x}^{(l)}-0.5|)$
(49)

And the outcome $Y=f_{\rho}(\mathbf{x})+\epsilon$ is contaminated by a normal
distributed noise, $\text{SNR}=3$. The main update rule of Sieve-SGD we
applied here is

$\hat{f}_{i}=\hat{f}_{i-1}+\gamma_{0}i^{-1/(2s+1)}\left(Y_{i}-\hat{f}_{i-1}\left(\mathbf{X}_{i}\right)\right)\sum_{\mathbf{j}\in\mathcal{J}_{prod}(cpi^{1/(2s+1)})}\left(\prod_{k=1}^{p}\mathbf{j}^{(k)}\right)^{-2\omega}\psi_{\mathbf{j}}\left(\mathbf{X}_{i}\right)\psi_{\mathbf{j}},$
(50)

where $\omega=0.51$ and $s=2$. We use $\gamma_{0}\in\\{0.1,0.5\\}$ and
$c\in\\{4,8\\}$: the latter two parameters may be different in each
replication. In Figure B.1, we present average performance of each method
based on 100 replications. The index set
$\mathcal{J}_{prod}(J)\subset\left(\mathbb{N}^{+}\right)^{p}$ contains the
$p$-dimension index vectors of smallest product. For example, when $p=2$,
$\mathcal{J}_{prod}(5)=\\{(1,1),(1,2),(2,1),(1,3),(3,1)\\}$. Arbitrary choice
is used when there is a tie. The univariate basis functions we use are

$\psi_{1}(x)=1,\quad\psi_{j}=\sqrt{2}\cos((j-1)\pi x),\text{ for }j\geq 2.$
(51)

In Figure B.1, we compare the statistical performance of SieveSGD with other
popular methods in statistics and computer science communities. The Sobolev
tensor product kernel we use there is
$K(\mathbf{s},\mathbf{t})=\prod_{k=1}^{p}\left(1+\min\\{\mathbf{s}^{(k)},\mathbf{t}^{(k)}\\}\right)$.
The RKHS corresponding to this kernel is the tensor product of univariate
Sobolev space on $[0,1]$. Like many other learning methods trained with
stochastic gradient descent, it is possible to have several pass over the data
set to achieve better generalization ability. While processing the data
multiple times, we continue to increase the number of basis function of Sieve-
SGD. That is, after 5 passes we use $cp(5\times 10^{5})^{1/(2s+1)}$ basis
functions. This strategy is not feasible for kernel SGD methods so for include
relevant results for comparison.

Figure B.1: Multivariate numerical examples of applying SieveSGD. Other
benchmark methods: ExpKRR, kernel ridge regression with Gaussian kernel;
KernelSGD, [16] with tensor product Sobolev kernel; GBM, gradient boosting
machine; RF, random forest; SobolevKRR, kernel ridge regression with tensor
product Sobolev kernel. We present the result of each method under oracle
hyperparameters (best testing error). The shaded area corresponds to the
second to fifth pass of SieveSGD over the same training data ($10^{5}$ unique
observations). (A) $p=2$ (B) $p=10$.

Now we present some intuition behind our choice of ordering the multi-variate
tensor product basis functions $\psi_{\mathbf{j}}$. While there are many
(equivalent) ways to arrive at such a choice, we choose to use some basic
theory in RKHS for easier exposition. While RKHSs are usually defined with the
reproducing property (as we will see very soon in Theorem C.2), there is
another characterization of the same function space which we will introduce
now to help our exposition.

Let $\nu$ be a positive Borel measure on the feature domain $\mathcal{X}$ that
has support$(\nu)=\mathcal{X}$, i.e. $\nu(U)>0$ for any nonempty open subset
$U$ of $\mathcal{X}\subset\mathbb{R}$. Let $K$ be a Mercer kernel. We use
$(\lambda_{j},\psi_{j})$ to denote the eigensystem of the covariance operator
$T_{\nu}$ of $K$ ($\psi_{j}$’s are orthonormal w.r.t. $\nu$). It is known [15]
that the following Hilbert space is the same as the one we described in
Theorem C.2.

######  Theorem B.1.

Define a Hilbert space

$\mathcal{H}=\left\\{f\in L^{2}_{\nu}\mid
f=\sum_{k=1}^{\infty}\theta_{j}\psi_{j}\
\text{with}\sum_{j=1}^{\infty}\left(\frac{\theta_{j}}{\sqrt{\lambda_{j}}}\right)^{2}<\infty\right\\}$
(52)

equipped with inner product:

$\langle
f,g\rangle_{\mathcal{H}}=\sum_{j=1}^{\infty}\frac{\theta_{j}\eta_{j}}{\lambda_{j}}$
(53)

for $f=\sum_{j=1}^{\infty}\theta_{j}\psi_{j}$ and
$g=\sum_{j=1}^{\infty}\eta_{j}\psi_{j}$, where $\theta_{j}=\langle
f,\psi_{j}\rangle_{L^{2}_{\nu}}$, $\eta_{j}=\langle
g,\psi_{j}\rangle_{L^{2}_{\nu}}$. Then
$(\mathcal{H},\langle\cdot,\cdot\rangle_{\mathcal{H}})$ is the reproducing
Hilbert space of kernel $K$. And we have the following Mercer expansion of the
kernel $K$:

$K(x,z)=\sum_{j=1}^{\infty}\lambda_{j}\psi_{j}(x)\psi_{j}(z)$ (54)

where the convergence is uniform and absolute.

If $\lambda_{j}=j^{-2s}$, then a ball in the kernel’s RKHS is the same as the
ellipsoid

$W(s,Q,\\{\psi_{j}\\})=\left\\{f\in L^{2}_{\nu}\mid
f=\sum_{j=1}^{\infty}\theta_{j}\psi_{j},\sum_{j=1}^{\infty}(\theta_{j}j^{s})^{2}\leq
Q^{2}\right\\}.$ (55)

If we consider the two-dimensional tensor product kernel
$\tilde{K}:\mathbb{R}^{2}\times\mathbb{R}^{2}\rightarrow\mathbb{R}$
constructed from $K$, that is

$\tilde{K}(\mathbf{x},\mathbf{z})=K(\mathbf{x}^{(1)},\mathbf{z}^{(1)})K(\mathbf{x}^{(2)},\mathbf{z}^{(2)})$
(56)

where $\mathbf{x}^{(j)}$ is the $j$-th component of
$\mathbf{x}\in\mathbb{R}^{2}$. It is known ([62] Section 12.4.2) that

$\tilde{K}(\mathbf{x},\mathbf{z})=\sum_{j=1}^{\infty}\sum_{k=1}^{\infty}(jk)^{-2s}\psi_{j}(\mathbf{x}^{(1)})\psi_{j}(\mathbf{z}^{(1)})\psi_{k}(\mathbf{x}^{(2)})\psi_{k}(\mathbf{z}^{(2)})$
(57)

and a ball in the RKHS of $\tilde{K}$ takes the form

$\tilde{W}=\left\\{f\in L^{2}_{\nu}\otimes L^{2}_{\nu}\mid
f(x)=\sum_{j,k=1}^{\infty}\theta_{jk}\psi_{j}\left(\mathbf{x}^{(1)}\right)\psi_{k}\left(\mathbf{x}^{(2)}\right),\sum_{j,k=1}^{\infty}\left(\theta_{jk}(jk)^{s}\right)^{2}\leq
Q^{2}\right\\},$ (58)

or equivalently, the eigenvalues are
$\prod_{k=1}^{2}\left(\mathbf{j}^{(k)}\right)^{-2s}$,
$\mathbf{j}\in(\mathbb{N}^{+})^{2}$. Accroding to (58), we would intuitively
expect $\theta_{jk}$ to be smaller when the product of the index vector is
larger.

When the univariate RKHS is a Sobolev space, estimating with a tensor product
kernel is essentially assuming the true regression function is (or can be
well-approximated by a function) in the tensor product space of Sobolev space.
The latter is also characterized as a Sobolev space with (dominating) mixed
derivatives [49]. In statistics, such a model has been studied under the name
of nonparametric Tensor product ANOVA [36]. Although methods engaging with
hyperbolic cross have been actively studied in numerical analysis in the past
decade, there are few works adopting such an idea into statistics. Sobolev
spaces with mixed derivatives are not homogeneous spaces in the sense that
they contain functions of different smoothness in different directions.
Specifically, functions in such spaces can be less smooth along the directions
of coordinate axis than other directions. This can be useful in the case when
the features as a strong “main effect” on the outcome and a weaker
“interaction effect” (in the sense of [36]).

###  B.2 Additive model and Sieve-SGD

In the main text Section 5.4, we described the nonparametric additive model
and how to use Sieve-SGD to estimate it. We simplified things by omitting the
intercept term to streamline exposition. The additive model with intercept is
given by

$f_{\rho}(\mathbf{x})=\beta^{0}+\sum_{k=1}^{p}f_{\rho,k}\left(\mathbf{x}^{(k)}\right)$
(59)

for some $\beta^{0}\in\mathbb{R}$ and $f_{\rho,k}\in
W_{k}(s_{k},Q_{k},\\{\psi_{jk}\\})$ for some centered $\psi_{jk}$
$\left(\int\psi_{jk}(x)d\nu(x)=0\right)$ (for example the functions in (51)).
For the additive model with intercept (59), the updating rule (24) of Sieve-
SGD could be replaced by a two-step procedure:

$\displaystyle\hat{f}_{i}$
$\displaystyle=\hat{f}_{i-1}+\gamma_{i}\left(Y_{i}-\hat{\beta}_{i-1}^{0}-\sum_{k=1}^{p}\hat{f}_{i-1,k}\left(\mathbf{X}_{i}^{(k)}\right)\right)\sum_{k=1}^{p}\sum_{j=1}^{J_{ik}}j^{-2s_{k}}\psi_{jk}\left(\mathbf{X}_{i}^{(k)}\right)\psi_{jk}$
(60) $\displaystyle\hat{\beta}^{0}_{i}$
$\displaystyle=\hat{\beta}_{i-1}^{0}+\gamma_{i}\left(Y_{i}-\hat{\beta}_{i-1}^{0}-\sum_{k=1}^{p}\hat{f}_{i-1,k}\left(\mathbf{X}_{i}^{(k)}\right)\right)$

here $J_{ik}$ is the truncation level of the $k$-th covariate when the sample
size is equal to $i$; and $\hat{f}_{i-1,k}$ is the estimate of $f_{\rho,k}$.
After applying Polyak averaging (averaging $\hat{\beta}^{0}_{i}+\hat{f}_{i}$
with previous estimates), we will get the Sieve-SGD estimate of $f_{\rho}$.

##  Appendix C Proof of Lemma 6.2

In this section, we will prove Lemma 6.2, together with results regarding the
spectrum of some related operators. To this end, we need to prepare the reader
by reminding them about some established results and ideas in the literature.
In this section we will

  * •

Define a Reproducing Kernel Hilbert Space (RKHS) formally;

  * •

Define covariance operators characterized by a kernel and discuss related
geometric properties, and;

  * •

Define the entropy of a compact operator and relate it to the eigenvalues of
the operator.

After all these, we will be ready to give a proof of Lemma 6.2.
In this section, we need to distinguish functions (in the RKHS) and their
equivalence class (in $L^{2}$ spaces) for a more rigorous discussion. For a
given measure $\mu$ on $\mathcal{X}\subset\mathbb{R}^{p}$, we use
$\mathcal{L}^{2}_{\mu}$ to denote the Hilbert space of all $\mu$-square
integrable functions. The $L^{2}_{\mu}$ spaces should be understand as the
quotient spaces of $\mathcal{L}^{2}_{\mu}$ under the equivalence relation:

$f=g\Leftrightarrow\int_{\mathcal{X}}\left(f(\tau)-g(\tau)\right)^{2}d\mu(\tau)=0.$
(61)

For a function $g\in\mathcal{L}^{2}_{\mu}$, we use $[g]=[g]_{\mu}\in
L^{2}_{\mu}$ to denote its equivalence class. This mathematical framework [53]
we present here allows the discussion when the measure $\mu$ does not have a
full-support over $\mathcal{X}$ (which is weaker than our Assumption A2), or
when the RKHS is not dense in $\mathcal{L}^{2}_{\mu}$.

###  C.1 Mercer kernel and RKHS

We first introduce the definition of a Mercer kernel and its corresponding
RKHS.

######  Definition C.1 (Mercer kernel).

A symmetric bivariate function
$K:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ is positive semi-
definite (PSD) if for any $n\geq 1$ and $(x_{i})_{i=1}^{n}\subset\mathcal{X}$,
the $n\times n$ matrix $\mathbb{K}$ whose elements are
$\mathbb{K}_{ij}=K(x_{i},x_{j})$ is always a PSD matrix.

A continuous, bounded, PSD kernel function $K$ is called a _Mercer kernel_.

In Assumption A3 of the main manuscript, we assumed $\\{\psi_{j}\\}$ is a set
of bounded, continuous functions in $\mathcal{L}^{2}_{\nu}$. For each
$J\in\mathbb{N}^{+}\cup\\{\infty\\}$ and $\omega>0.5$, we can show that the
bivariate functions

$K_{J}^{\omega}(s,t):=\sum_{j=1}^{J}j^{-2\omega}\psi_{j}(s)\psi_{j}(t),\quad$
(62)

are Mercer kernels. We also use $K(s,t)=K_{\infty}^{s}(s,t)$ to denote the
canonical (untruncated) kernel in our analysis.

It is well-known that for any Mercer kernel, there is a unique associated
Hilbert space $(\mathcal{H}_{K},\langle\cdot,\cdot\rangle_{K})$ which has the
so-called reproducing property. The following theorem formally defines such a
Hilbert space and states its uniqueness.

######  Theorem C.2 ([15]).

For a Mercer Kernel $k:\mathbb{X}\times\mathbb{X}\rightarrow\mathbb{R}$, there
exists an unique Hilbert Space
$(\mathcal{H}_{k},\langle\cdot,\cdot\rangle_{k})$ of functions on $\mathbb{X}$
satisfying the following conditions. Let $k_{x}:z\mapsto k(x,z)$:

  1. 1.

For all $x\in\mathbb{X}$, $k_{x}\in\mathcal{H}_{k}$.

  2. 2.

The linear span of $\\{k_{x}\ |\ x\in\mathbb{X}\\}$ is dense (w.r.t
$\|\cdot\|_{k}$) in $\mathcal{H}_{k}$

  3. 3.

For all $f\in\mathcal{H}_{k},x\in\mathbb{X}$, $f(x)=\langle
f,k_{x}\rangle_{k}$

We call this Hilbert space the _Reproducing kernel Hilbert space (RKHS)_
associated with kernel $k$ .

Note that in the above definition, we did not mention any measures on
$\mathbb{X}$. The RKHS $\mathcal{H}_{k}$ of a kernel $k$ is defined
independently as a Hilbert space (a complete linear space equipped with an
inner product).

The inner product in Theorem C.2 is implicitly defined and appear to be quite
abstract. However, there is an equivalent definition of the RKHS corresponding
to $K$ in which our Sobolev ellipsoid assumption appears:

######  Theorem C.3 (p.37,Theorem 4 in [15]).

The Hilbert space defined in Theorem C.2 is identical (same function class
with the same inner product) to the following Hilbert space $\mathbb{H}_{K}$.

$\mathbb{H}_{K}=\left\\{f\in\mathcal{L}_{\nu}^{2}(X)\mid
f=\sum_{j=1}^{\infty}a_{j}\psi_{j}\quad\text{ with
}\quad\sum_{j=1}^{\infty}\left(j^{s}a_{j}\right)^{2}<\infty\right\\}$ (63)

Equipped with the inner product:

$\langle f,g\rangle_{K}=\sum_{j=1}^{\infty}j^{2s}a_{j}b_{j}$ (64)

for $f=\sum a_{j}\psi_{j},\text{ and }g=\sum b_{j}\psi_{j}$.

It is direct to check our assumption A3 is the same as assuming the
conditional mean $f_{\rho}$ belongs to a ball of radius $Q$ in the above
constructed RKHS.

In this section, we consider the RKHSs related to kernels $K_{J}$ with
$\mathbb{X}=\mathcal{X}$, where $\mathcal{X}$ is the support of $\nu$. Under
our assumption A2 and A3, the functions in $\mathcal{H}_{K}$ are all square-
integrable. We also define the identity mapping (w.r.t. measure $\mu$)
$\text{id}_{\mu}$ as

$\displaystyle\text{id}_{\mu}:\mathcal{H}_{K}$ $\displaystyle\rightarrow
L^{2}_{\mu}$ (65) $\displaystyle g$ $\displaystyle\mapsto[g]_{\mu}.$

###  C.2 Covariance operators

In Section C.2 and C.3, we engage with the RKHS of the canonical kernel $K$.
Once the properties of its related operators are clear, we can directly
generate our analysis techniques to other truncated kernels $K_{J}^{\omega}$.
Recall our definitions of $T_{X}$ and $T_{\nu}$ in the main text:

$\displaystyle T_{X}:L_{\rho_{X}}^{2}$ $\displaystyle\rightarrow
L_{\rho_{X}}^{2}$ (66) $\displaystyle g$
$\displaystyle\mapsto\int_{\mathcal{X}}g(\tau)K(\tau,\cdot),d\rho_{X}(\tau)$

and

$\displaystyle T_{\nu}:L^{2}_{\nu}$ $\displaystyle\rightarrow L^{2}_{\nu}$
(67) $\displaystyle g$ $\displaystyle\mapsto\int
g(\tau)K(\tau,\cdot)d\nu(\tau).$

Now we state several basic properties of $T_{X}$. Similar properties also hold
for $T_{\nu}$ and can be verified much easier without abstract analysis. For
proofs of Lemma C.4 and more properties, see [53, Section 2 & 3].

######  Lemma C.4.

Under Assumptions A2, A3 in the main text:

  * •

The operator $T_{X}$ is bounded, self-adjoint, positive.

  * •

There exists a at most countable set of functions
$\\{\phi_{j}\\}\subset\mathcal{H}_{K}$, $j\in\mathcal{J}$ and a at most
countable sequence of positive numbers $\lambda_{j}$ (decreasingly ordered)
such that

$T_{X}(g)=\sum_{j\in\mathcal{J}}\lambda_{j}\langle
g,\left[\phi_{j}\right]_{\rho_{X}}\rangle_{L^{2}_{\rho_{X}}}\left[\phi_{j}\right]_{\rho_{X}},\quad\text{for
}g\in L^{2}_{\rho_{X}}.$ (68)

  * •

The $\\{[\phi_{j}]_{\rho_{X}}\\}$ above is an orthonormal system in
$L^{2}_{\rho_{X}}$ and $\\{\sqrt{\lambda_{j}}\phi_{j}\\}$ is an orthonormal
system in $\mathcal{H}_{K}$. Therefore, $(\lambda_{j},[\phi]_{j})$ is an
eigensystem of $T_{X}$:

$T_{X}([\phi_{j}]_{\rho_{X}})=\lambda_{j}[\phi_{j}]_{\rho_{X}}$ (69)

  * •

$T_{X}$ is a trace class operator, i.e. $\sum_{j}\lambda_{j}<\infty$.

Under the assumptions A2, A3, we can actually say more about the properties of
$\phi_{j}$. But we list them as a separate lemma since they are not
necessarily true when we discuss truncated kernels later.

######  Lemma C.5 (Theorem 3.1, [53]).

Under the same assumptions as in Lemma C.4. Let $\phi_{j}$ and
$[\phi_{j}]_{\rho_{X}}$ denote the orthonormal systems in Lemma C.4. Then

  * •

The family $[\phi_{j}]_{\rho_{X}}$ is an orthonormal basis of
$L^{2}_{\rho_{X}}$.

  * •

The family $\sqrt{\lambda_{j}}\phi_{j}$ is an orthonormal basis of
$\mathcal{H}_{K}$.

Now we define several operators related to $T_{X}$ that will facilitate our
analysis of the spectrum of it.

######  Definition C.6.

Under the same assumptions as in Lemma C.4, with the same $\\{\lambda_{j}\\}$
and $\\{\phi_{j}\\}$:

  * •

We define the $r$-th power of $T_{X}$ as:

$\displaystyle T^{r}_{X}:L^{2}_{\rho_{X}}$ $\displaystyle\rightarrow
L^{2}_{\rho_{X}}$ (70) $\displaystyle g$
$\displaystyle\mapsto\sum_{j}\lambda_{j}^{r}\langle
g,\left[\phi_{j}\right]_{\rho_{X}}\rangle_{L^{2}_{\rho_{X}}}\left[\phi_{j}\right]_{\rho_{X}},\quad\text{for
}g\in L^{2}_{\rho_{X}},$

in this work, we are most interested in the square-root of $T_{X}$, i.e.
$T_{X}^{1/2}$.

  * •

Define the operator $S_{X}^{1/2}$ as:

$\displaystyle S^{1/2}_{X}:L^{2}_{\rho_{X}}$
$\displaystyle\rightarrow\mathcal{H}_{K}$ (71) $\displaystyle g$
$\displaystyle\mapsto\sum_{j}\lambda_{j}^{1/2}\langle
g,\left[\phi_{j}\right]_{\rho_{X}}\rangle_{L^{2}_{\rho_{X}}}\phi_{j},\quad\text{for
}g\in L^{2}_{\rho_{X}}$

The operator $S_{X}^{1/2}$ has a very importance geometric properties: it
preserves distance between two subspaces of $L^{2}_{\rho_{X}}$ and
$\mathcal{H}_{K}$, as stated in the following lemma.

######  Lemma C.7 ([53],Theorem 2.11).

Under the same assumptions as in Lemma C.4, let $S_{X}^{1/2}$ be the operator
defined in (71). Then

  * •

$S_{X}^{1/2}$ is bijective between
$\overline{\text{span}\\{[\phi_{j}]_{\rho_{X}},j\in\mathcal{J}\\}}=L^{2}_{\rho_{X}}$
and $\overline{\text{span}\\{\phi_{j},j\in\mathcal{J}\\}}=\mathcal{H}_{K}$
and,

  * •

$\|S_{X}^{1/2}(g)\|_{K}=\|g\|_{L^{2}_{\rho_{X}}}$, for
$g\in\overline{\text{span}\\{[\phi_{j}]_{\rho_{X}},j\in\mathcal{J}\\}}$.

That is, $S_{X}^{1/2}$ is an isometric isomorphism between $L^{2}_{\rho_{X}}$
and $\mathcal{H}_{K}$.

It is direct to check the following lemma by the equivalent definition of the
RKHS (Theorem C.3).

######  Lemma C.8.

Define $S_{\nu}^{1/2}$ similarly for the operator $T_{\nu}$:

$\displaystyle S^{1/2}_{\nu}:L^{2}_{\nu}$
$\displaystyle\rightarrow\mathcal{H}_{K}$ (72) $\displaystyle g$
$\displaystyle\mapsto\sum_{j}j^{-s}\langle
g,\left[\psi_{j}\right]_{\nu}\rangle_{L^{2}_{\nu}}\psi_{j},\quad\text{for
}g\in L^{2}_{\nu}$

Then

  * •

$S_{\nu}^{1/2}$ is bijective between $L^{2}_{\nu}$ and $\mathcal{H}_{K}$ and,

  * •

$\|S_{\nu}^{1/2}(g)\|_{K}=\|g\|_{L^{2}_{\nu}}$, for $g\in L^{2}_{\nu}$.

###  C.3 Entropy of an operator and its spectrum

The above Lemma C.7 and Lemma C.8 is one set of the elements we are going to
use in the proof of Lemma 6.2. Another important part of our proof is the
correspondence between the spectrum of an operator and its metric entropy. We
believe that these results might help show connections between proof methods
regarding nonparametric problems that using the spectrum of operators [16, 66]
and those using metric entropy [62, 25].

We first define the metric entropy of an operator. There is a correspondence
between our definition and the ”standard” definition of metric entropy for
compact sets. In the following we will use $B_{E}$ to denote the unit ball of
a function space $E$.

######  Definition C.9 (Entropy of an operator).

For $k\geq 1$ we define the $k$-th entropy number of a metric space $S$ to be

$e_{k}(S)=\inf\left\\{\epsilon>0|\exists\text{ closed balls
}D_{1},\ldots,D_{2^{k}-1}\text{ with radius }\epsilon\text{ covering
}S\right\\}$ (73)

If $T:E\rightarrow F$ is a linear map, then we define the $k$-th entropy
number of $T$ as

$e_{k}(T)=e_{k}\left(T\left(B_{E}\right)\right)$ (74)

The first result bounds the entropy of an operator by its eigenvalues.

######  Theorem C.10.

Let $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{j}\geq\ldots\geq 0$ be a
sequence of real numbers, $w_{j}$ be a sequence of elements in $l^{2}$ (space
of square-summable sequences). Consider the operator defined by

$\displaystyle T:l^{2}$ $\displaystyle\rightarrow l^{2}$ (75)
$\displaystyle(\xi_{1},\xi_{2},...,\xi_{j},...)$
$\displaystyle\mapsto(\lambda_{1}\xi_{1},\lambda_{2}\xi_{2},...,\lambda_{j}\xi_{j},...)$

If $\lambda_{j}\leq Cj^{-s}$ for some $C,s$ and all $j\geq 1$, then for all
$j\geq 2$

$e_{j}(T)\leq 12Cs^{s}j^{-s}$ (76)

For proof, see Proposition 9 in [15, Appendix A]. The proof there uses
Proposition 1.3.2 in [12], which engages with $l^{2}$ spaces as well. Theorem
C.10 is good enough for our purpose because we can diagonalize the covariance
operators of interest (Lemma C.4). For a general result of the same nature as
Proposition 1.3.2, one can use Proposition 3.4.2 combined with Proposition
4.4.1 in [12].

We also state a result in the other direction: Carl’s inequality, see [11,
Theorem 4].

######  Theorem C.11.

Let $E$ be a Banach space and let $T:E\rightarrow E$ be a linear compact
operator. Then

$\left|\lambda_{j}(T)\right|\leq\sqrt{2}e_{j}(T)$ (77)

where $\lambda_{j}(T)$ is the non-increasing sequence of eigenvalues of $T$.

Now we are ready to prove Lemma 6.2 in the main text.

###### Proof of Lemma 6.2.

We use $B_{X}$ to denote the unit ball in $L^{2}_{\rho_{X}}$ and $B_{\nu}$ for
the unit ball in $L^{2}_{\nu}$. The proof of this lemma can be divided into
three steps.
Step 1: bound entropy by spectrum. By the definition of the power of an
operator, the eigenvalue of $T^{1/2}_{\nu}$ is $j^{-s}$ (because the
eigenvalue of $T_{\nu}$ is $j^{-2s}$). And the corresponding eigenfunctions
are $\psi_{j}$. For any function $f\in L^{2}_{\nu}$, there is a unique
sequence $(f_{j})_{j=1}^{\infty}\in l^{2}$ such that
$f=\sum_{j=1}^{\infty}f_{j}\psi_{j}$. Apply Theorem C.10 to $T^{1/2}_{\nu}$ we
know:

$e_{j}(T^{1/2}_{\nu})\leq 12s^{s}j^{-s}$ (78)

Step 2: relate the entropy of operators. Now we are going to show that
$e_{j}(T^{1/2}_{X})\leq ue_{j}(T^{1/2}_{\nu})$.To investigate the entropy of
an operator, we only need to look at the entropy of $T^{1/2}_{X}(B_{X})$ and
$T^{1/2}_{\nu}(B_{\nu})$. By Lemma C.7 and Lemma C.8, we know
$S_{X}^{1/2}(B_{X})=S_{\nu}^{1/2}(B_{\nu})=B_{\mathcal{H}_{K}}$. When
measuring the entropy of $T^{1/2}_{X}(B_{X})$ we need to use $\text{id}_{X}$
to map it to $L^{2}_{\rho_{X}}$, but for $T^{1/2}_{\nu}(B_{\nu})$ we need to
use $\text{id}_{\nu}$ to map it to $L^{2}_{\nu}$.

By the assumption A2 of $\rho_{X}$, for any $f\in
L^{2}_{\rho_{X}}=L^{2}_{\nu}$, we have

$\int[f]_{\rho_{X}}^{2}(x)d\rho_{X}(x)\leq u^{2}\int[f]_{\nu}^{2}(x)d\nu(x)$
(79)

This means we can use the center of an $\epsilon$-cover of
$T^{1/2}_{\nu}(B_{\nu})$ to construct a $u\epsilon$-cover of
$T^{1/2}_{X}(B_{X})$ using the same center points. By the definition of
entropy,

$e_{j}(T^{1/2}_{X})\leq ue_{j}(T^{1/2}_{\nu})$ (80)

Step 3: bound spectrum by entropy. We use Theorem C.11 to translate our bound
from entropy back to the spectrum of our operator. For the eigenvalue
$\lambda_{j}$ of $T_{X}$, we have

$\lambda_{j}=\left(\lambda_{j}^{1/2}\right)^{2}\leq\left(\sqrt{2}e_{j}(T^{1/2}_{X})\right)^{2}\leq\left(\sqrt{2}ue_{j}(T^{1/2}_{\nu})\right)^{2}\leq\left(12\sqrt{2}us^{s}j^{-s}\right)^{2}$
(81)

This concludes our claim that the eigenvalues of $T_{X}$ satisfy
$\lambda_{j}\leq 288u^{2}s^{2s}j^{-2s}=O(j^{-2s})$.

Similarly, we can show $j^{-2s}\leq 288\ell^{2}s^{2s}\lambda_{j}$ by going in
the other direction. So we conclude that $\lambda_{j}=\Theta(j^{-2s})$. ∎

###  C.4 Technical Results

In this subsection we will present several results needed in the proof of
Theorem 6.1. Since they are related to the spectrum of covariance operators,
we put them here for rather than the technical section of Appendix D.

For a fixed $1\leq J\leq\infty$ and $\omega>\frac{1}{2}$, we consider the
kernel $K_{J}^{\omega}(s,t)=\sum_{j=1}^{J}j^{-2\omega}\psi_{j}(s)\psi_{j}(t)$.
By Theorem C.2, there is an unique Hilbert space with the reproducing
properties. According to Theorem C.3, we have the hierarchy relation that
$\mathcal{H}_{K_{J}^{\omega}}\subset\mathcal{H}_{K_{I}^{\omega}}$ for $I\geq
J$, equipped with the inner product $\langle
f,g\rangle_{K}=\sum_{j=1}^{I}j^{2\omega}a_{j}b_{j}$ for $f=\sum
a_{j}\psi_{j}\text{ and }g\sum b_{j}\psi_{j}$. Similarly, we can define the
covariance operators w.r.t. $\rho_{X}$ and $\nu$:

$\displaystyle T_{X,J}^{\omega}:L_{\rho_{X}}^{2}$ $\displaystyle\rightarrow
L_{\rho_{X}}^{2}$ (82) $\displaystyle g$
$\displaystyle\mapsto\int_{\mathcal{X}}g(\tau)K_{J}^{\omega}(\tau,\cdot),d\rho_{X}(\tau)$

and

$\displaystyle T_{\nu,J}^{\omega}:L^{2}_{\nu}$ $\displaystyle\rightarrow
L^{2}_{\nu}$ (83) $\displaystyle g$ $\displaystyle\mapsto\int
g(\tau)K_{J}^{\omega}(\tau,\cdot)d\nu(\tau).$

Similar to Lemma C.4, we can diagonalize $T_{X,J}^{\omega}$ with an
eigensystem $(\lambda_{J,j}^{\omega},\phi_{J,j}^{\omega})$. Different from
$\\{\phi_{j}\\}$, $\\{\phi_{J,j}^{\omega}\\}$ is a basis of
$\mathcal{H}_{K_{J}^{\omega}}$ but $\\{[\phi_{J,j}^{\omega}]_{\rho_{X}}\\}$ is
not a basis of $L^{2}_{\rho_{X}}$. We formally state it in the following
lemma.

######  Lemma C.12 (Theorem 3.1, [53]).

Under assumptions A2, A3. Let $(\lambda_{J,j}^{\omega},\phi_{J,j}^{\omega})$
denote the eigensystem of $T_{X,J}^{\omega}$. Then

  * •

The family $[\phi_{J,j}^{\omega}]_{\rho_{X}}$ is an orthonormal system of
$L^{2}_{\rho_{X}}$.

  * •

The family $\sqrt{\lambda_{J,j}^{\omega}}\phi_{J,j}^{\omega}$ is an
orthonormal basis of $\mathcal{H}_{K_{J}^{\omega}}$.

Related to Lemma C.12, the mapping $\left(S_{X,J}^{\omega}\right)^{1/2}$ is
not an isomorphism between $L^{2}_{\rho_{X}}$ and
$\mathcal{H}_{K_{J}^{\omega}}$ – it has a non-trivial kernel space.

######  Lemma C.13 ([53],Theorem 2.11).

Under the same assumptions as in Lemma C.12, define
$\left(S_{X,J}^{\omega}\right)^{1/2}$ as

$\displaystyle\left(S_{X,J}^{\omega}\right)^{1/2}:L^{2}_{\rho_{X}}$
$\displaystyle\rightarrow\mathcal{H}_{K_{J}^{\omega}}$ (84) $\displaystyle g$
$\displaystyle\mapsto\sum_{j}\left(\lambda_{J,j}^{\omega}\right)^{1/2}\langle
g,\left[\phi_{J,j}^{\omega}\right]_{\rho_{X}}\rangle_{L^{2}_{\rho_{X}}}\phi_{J,j}^{\omega},\quad\text{for
}g\in L^{2}_{\rho_{X}}$

Then

  * •

$\left(S_{X,J}^{\omega}\right)^{1/2}$ is bijective between
$\overline{\text{span}\\{[\phi_{J,j}^{\omega}]_{\rho_{X}},j\in\mathcal{J}\\}}\subset
L^{2}_{\rho_{X}}$ and
$\overline{\text{span}\\{\phi_{J,j}^{\omega},j\in\mathcal{J}\\}}=\mathcal{H}_{K_{J}^{\omega}}$
and,

  * •

$\left\|\left(S_{X,J}^{\omega}\right)^{1/2}(g)\right\|_{K}=\|g\|_{L^{2}_{\rho_{X}}}$,
for
$g\in\overline{\text{span}\\{[\phi_{J,j}^{\omega}]_{\rho_{X}},j\in\mathcal{J}\\}}$.

Now we state and proof the main result of this section:

######  Lemma C.14.

Let $1\leq J\leq\infty$, $\omega>\frac{1}{2}$, assume A2 & A3. We use
$(\lambda_{J,j}^{\omega},\phi_{J,j}^{\omega})$ to denote the eigensystem of
$T_{X,J}^{\omega}$ (similarly defined as in Lemma C.4). Then

$(288\ell^{2}\omega^{2\omega})^{-1}j^{-2\omega}\leq\lambda_{J,j}^{\omega}\leq
288u^{2}\omega^{2\omega}j^{-2\omega},\quad\text{for }1\leq j\leq J$ (85)

###### Proof.

The proof assembles that of Lemma 6.2. To investigate the eigenvalues of
$T_{X,J}^{\omega}$, we just need to compare the entropy of
$\left(S_{X,J}^{\omega}\right)^{1/2}(B_{X})$ and
$\left(S_{\nu,J}^{\omega}\right)^{1/2}(B_{\nu})$, where $B_{\mu}$ is the unit
ball in $L^{2}_{\mu}$. We know

$\left(S_{X,J}^{\omega}\right)^{1/2}(B_{X})\stackrel{{\scriptstyle(1)}}{{=}}\text{unit
ball in
}\mathcal{H}_{K_{J}^{\omega}}=\left(S_{\nu,J}^{\omega}\right)^{1/2}(B_{\nu})$
(86)

In (1) we used the distribution assumption A2, which ensures
$\phi_{J,j}^{\omega}$ is a basis of $\mathcal{H}_{K_{J}^{\omega}}$ (Lemma
C.13). After embedding the unit ball in $\mathcal{H}_{K_{J}^{\omega}}$ back to
the $L^{2}$ spaces, we know a $\epsilon$-covering of
$T_{\nu,J}^{\omega}(B_{\nu})$ can generate an $u\epsilon$-covering of
$T_{X,J}^{\omega}(B_{X})$, which gives the above upper bound by similar
argument as in the proof of Lemma 6.2. It is also similar to show the lower
bound in (85), noting that the feature distribution is assumed in A2 to have a
strictly positive density. ∎

Note: the constants show up in (85) does not depend on the truncation level
$J$. Therefore, for a given $\omega$, we can treat the result in Lemma C.14 as
an uniform bound which is applicable to all truncation levels.

##  Appendix D Proof of Theorem 6.1

In this section we are going to prove the main performance guarantees, results
Theorem 6.1. In our proof, we first split the error into two parts: one part
is noiseless and depends only on the the initial bias, the other is due to the
noise in our data. We will bound each term separately and choose the learning
rate $\gamma_{n}$ to balance the trade-off. The last part of this section will
give some technical lemmas that will be referred to in the proofs of Theorem
6.1. Some of the proof techniques are taken from [3] and [16].

###  D.1 Notation

In this section, the RKHS we are considering is the one associated with mercer
kernel $K(s,t)=\sum_{j=1}^{\infty}j^{-2s}\psi_{j}(s)\psi_{j}(t)$. We use
$\|\cdot\|_{K}$ and $\langle\cdot,\cdot\rangle_{K}$ to denote the RKHS-norm
and RKHS-inner product. For any elements $g,h\in\mathcal{H}_{K},$ we define
the operator $g\otimes h$ as a mapping from $\mathcal{H}_{K}$ to
$\mathcal{H}_{K}$ such that $(g\otimes h)f=\langle f,h\rangle_{K}\ g$. As we
will show in Lemma D.1, the quantity $\|K_{X_{n},J_{n}}\|_{K}^{2}$ is bounded
(and this bound only depends on $s$). We use $R^{2}$ to denote the smallest
bound for $\|K_{X_{n},J_{n}}\|_{K}^{2}$. And any $\gamma_{n}$ in this section
is assumed to satisfy: $\gamma_{0}R^{2}<1$ and $\gamma_{n}\asymp
n^{-\frac{1}{2s+1}}$.
We consider a filtration of $\sigma$-algebras $\\{\mathcal{F}_{n}\\}$, where
$\mathcal{F}_{n}$ is the $\sigma$-algebra generated by
$(X_{i},Y_{i})_{i=1}^{n}$.
The sign $\preccurlyeq$ denotes the order between self-adjoint operators over
the RKHS. That is, for self-adjoint operators
$A,B:E\rightarrow\mathcal{H}_{K}$, $A\preccurlyeq B$ means $\langle
f,(B-A)f\rangle_{K}\geq 0$ for any $f\in\mathcal{H}_{K}$. Intuitively we can
think of $B-A$ as a positive semi-definite matrix in a finite-dimensional
space. The expectation of random function/operator should be understood as the
Bochner integral, a generalization of the Lebesgue integral where the random
element takes value in a Banach space, see [43, 14] or Chapter 4 of [5]. The
$\zeta(s)$ function that shows up in this section is the Riemann-zeta function
$\zeta(s):=\sum_{k=1}^{\infty}k^{-s}$. We use it for simplifying the notation
— it’s neat that it shows up, but we do not need any of the exciting (and
difficult to show) properties that number theorists/combinatorists are
concerned with!

###  D.2 Separation of the error

In our theoretical analysis of the generalization error, rather than study how
$\bar{f}_{n}$ converge to $f_{\rho}$, we instead directly study how the
difference shrinks to zero. We define

$\displaystyle\Delta_{n}$ $\displaystyle=\hat{f}_{n}-f_{\rho}$ (87)
$\displaystyle\bar{\Delta}_{n}$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{n}\Delta_{j}=\bar{f}_{n}-f_{\rho}$

We can relate the $\|\cdot\|_{2}$-norm (the natural norm shows up in the
regression problem) and the $\|\cdot\|_{K}$-norm (which facilitates our
theoretical analysis) using $T_{X}^{1/2}$. We will be repeatedly using the
following equivalence in our proof:

$\left\|g\right\|_{2}^{2}=\left\|T^{1/2}_{X}g\right\|_{K}^{2}=\left\langle
g,T_{X}g\right\rangle_{K}\quad\text{ for }g\in\mathcal{H}_{K}$ (88)

Similarly,

$\left\|g\right\|_{2}^{2}=\left\|T^{1/2}_{X,J_{n}}g\right\|_{K}^{2}=\left\langle
g,T_{X,J_{n}}g\right\rangle_{K}\quad\text{ for
}g\in\mathcal{H}_{K}\cap\text{span}\\{\psi_{1},...,\psi_{J_{n}}\\}$ (89)

And we have a recursive relationship for $\eta_{n}$ based on the recursive
relationship for $\hat{f}_{n}$:

$\displaystyle\hat{f}_{n}$
$\displaystyle=\hat{f}_{n-1}-\gamma_{n}(Y_{n}-\hat{f}_{n-1}(X_{n}))K_{X_{n},J_{n}}$
(90) $\displaystyle\Rightarrow\hat{f}_{n}$
$\displaystyle=\left(I-\gamma_{n}T_{X_{n},J_{n}}\right)\hat{f}_{n-1}+\gamma_{n}Y_{n}K_{X_{n},J_{n}}$
$\displaystyle\Rightarrow\hat{f}_{n}-f_{\rho}$
$\displaystyle=\left(I-\gamma_{n}T_{X_{n},J_{n}}\right)\left(\hat{f}_{n-1}-f_{\rho}\right)+\gamma_{n}\Xi_{n}$

where

$\displaystyle T_{X_{n},J_{n}}(f)=f(X_{n})K_{X_{n},J_{n}}$ (91)
$\displaystyle\Xi_{n}=\left(Y_{n}-f_{\rho}\left(X_{n}\right)\right)K_{X_{n},J_{n}}$

Thus, we have a recursive formula for $\Delta_{n}$:

$\displaystyle\Delta_{0}$ $\displaystyle=-f_{\rho}$ (92)
$\displaystyle\Delta_{n}$
$\displaystyle=\left(I-\gamma_{n}T_{X_{n},J_{n}}\right)\Delta_{n-1}+\gamma_{n}\Xi_{n}$

We further decompose $\Delta_{n}$ into two parts
$\Delta_{n}=\eta_{n}+\vartheta_{n}$:

  1. 1.

$\left(\eta_{n}\right)\text{ is defined as }$:

$\displaystyle\eta_{0}$ $\displaystyle=-f_{\rho}$ (93) $\displaystyle\eta_{n}$
$\displaystyle=\left(I-\gamma_{n}T_{X_{n},J_{n}}\right)\eta_{n-1}$

It is the part of $\Delta_{n}$ due to an initial value not equal to
$f_{\rho}$. We note that it does not contain the noise term $\Xi_{n}$, so the
only randomness comes from features $X_{n}$.

  2. 2.

The pure noise component $\left(\vartheta_{n}\right)$ is defined as:

$\displaystyle\vartheta_{0}$ $\displaystyle=0$ (94)
$\displaystyle\vartheta_{n}$
$\displaystyle=\left(I-\gamma_{n}T_{X_{n},J_{n}}\right)\vartheta_{n-1}+\gamma_{n}\Xi_{n}$

We can directly show that
$\bar{\Delta}_{n}=\bar{\eta}_{n}+\bar{\vartheta}_{n}$. By Minkowski’s
inequality:

$\left(E\left[\left\|\bar{\Delta}_{n}\right\|_{2}^{2}\right]\right)^{1/2}\leq\left(E\left[\left\|\bar{\eta}_{n}\right\|_{2}^{2}\right]\right)^{1/2}+\left(E\left[\left\|\bar{\vartheta}_{n}\right\|_{2}^{2}\right]\right)^{1/2}$
(95)

Our job now is just to bound the two terms separately and then choose the
correct $\gamma_{n}$ to minimize the combined bound.

###  D.3 Bound on initial condition sub-process

In this section, we will engage with bounding $\eta_{n}$, which is the part of
error due to the imperfect initialization.

###### proof of bound on initial value.

By definition,

$\eta_{n}=\eta_{n-1}-\gamma_{n}\eta_{n-1}(X_{n})K_{X_{n},J_{n}}$ (96)

We square both sides w.r.t. the RKHS inner product

$\|\eta_{n}\|_{K}^{2}=\|\eta_{n-1}\|_{K}^{2}-2\gamma_{n}\eta_{n-1}(X_{n})\langle\eta_{n-1},K_{X_{n},J_{n}}\rangle_{K}+\gamma_{n}^{2}\eta_{n-1}^{2}(X_{n})\left\|K_{X_{n},J_{n}}\right\|_{K}^{2}$
(97)

We take the expectation on both sides: conditioned on $\mathcal{F}_{n-1}$
first, then unconditionally.

$\displaystyle E\|\eta_{n}\|_{K}^{2}=$ $\displaystyle
E\|\eta_{n-1}\|_{K}^{2}-2\gamma_{n}E\left[\eta_{n-1}(X_{n})\langle\eta_{n-1},K_{X_{n},J_{n}}\rangle_{K}\right]$
(98)
$\displaystyle+\gamma_{n}^{2}E\left[\eta_{n-1}^{2}(X_{n})\left\|K_{X_{n},J_{n}}\right\|_{K}^{2}\right]$
$\displaystyle\stackrel{{\scriptstyle(1)}}{{\leq}}$ $\displaystyle
E\|\eta_{n-1}\|_{K}^{2}-2\gamma_{n}E\left[\eta_{n-1}(X_{n})\langle\eta_{n-1},K_{X_{n},J_{n}}\rangle_{K}\right]+\gamma_{n}E\|\eta_{n-1}\|_{2}^{2}$

in $(1)$ we use the fact that for any $n$,
$\gamma_{n}\left\|K_{X_{n},J_{n}}\right\|_{K}^{2}\leq\gamma_{0}R^{2}<1$. This
is actually how we choose our $(\gamma_{n})$. If the $K_{X_{n},J_{n}}$ in the
middle term above is actually $K_{X_{n}}$, then the whole middle term will
become just $E\|\eta_{n-1}\|_{2}^{2}$, which makes the following algebra
easier. However, since we do not use exactly $K_{X_{n}}$ but truncated at
level $J_{n}$, we need some extra effort to deal with it.

$\displaystyle E\|\eta_{n}\|_{K}^{2}$ $\displaystyle\leq
E\|\eta_{n-1}\|_{K}^{2}-2\gamma_{n}E\left[\eta_{n-1}(X_{n})\langle\eta_{n-1},K_{X_{n}}+K_{X_{n},J_{n}}-K_{X_{n}}\rangle_{K}\right]+\gamma_{n}E\|\eta_{n-1}\|_{2}^{2}$
(99)
$\displaystyle=E\|\eta_{n-1}\|_{K}^{2}-\gamma_{n}E\|\eta_{n-1}\|_{2}^{2}+2\gamma_{n}E\left[\eta_{n-1}(X_{n})\langle\eta_{n-1},K_{X_{n}}-K_{X_{n},J_{n}}\rangle_{K}\right]$
$\displaystyle\stackrel{{\scriptstyle(1)}}{{\leq}}E\|\eta_{n-1}\|_{K}^{2}-\gamma_{n}E\|\eta_{n-1}\|_{2}^{2}+\frac{1}{2}\gamma_{n}E\|\eta_{n-1}\|_{2}^{2}+2\gamma_{n}E\left[\langle\eta_{n-1},K_{X_{n}}-K_{X_{n},J_{n}}\rangle_{K}^{2}\right]$

In (1) we use Young’s inequality. Now we bound the last term:

$\displaystyle
E\left[\langle\eta_{n-1},K_{X_{n}}-K_{X_{n},J_{n}}\rangle_{K}^{2}\right]$
$\displaystyle=E\left[\langle\eta_{n-1},\sum_{j=J_{n}+1}^{\infty}j^{-2s}\psi_{j}(X_{n})\psi_{j}\rangle_{K}^{2}\right]$
(100)
$\displaystyle=E\left[\left(\sum_{j=J_{n}+1}^{\infty}\eta_{n-1,j}\psi_{j}(X_{n})\right)^{2}\right]\text{\
where\ }\eta_{n-1,j}=\langle\eta_{n-1},\psi_{j}\rangle_{L^{2}_{\nu}}$
$\displaystyle\leq
uE\left[\int_{\mathcal{X}}\left(\sum_{j=J_{n}+1}^{\infty}\eta_{n-1,j}\psi_{j}(x)\right)^{2}d\nu(x)\right]$
$\displaystyle\leq
uJ_{n}^{-2s}E\left[\sum_{j=J_{n}+1}^{\infty}\left(j^{s}\eta_{n-1,j}\right)^{2}\right]\leq
uJ_{n}^{-2s}E\|\eta_{n-1}\|_{K}^{2}$

Continue (99):

$E\|\eta_{n}\|_{K}^{2}\leq
E\|\eta_{n-1}\|_{K}^{2}-\frac{1}{2}\gamma_{n}E\|\eta_{n-1}\|_{2}^{2}+2u\gamma_{n}J_{n}^{-2s}E\|\eta_{n-1}\|_{K}^{2}$
(101)

Now for each $i$ we have such a recursive relationship for
$\|\eta_{i}\|_{K}^{2}$, $\|\eta_{i-1}\|_{K}^{2}$ and $\|\eta_{i-1}\|_{2}^{2}$.
We can sum this from $i=1$ to $n$.

$\displaystyle E\|\eta_{n}\|_{K}^{2}$ $\displaystyle\leq
E\|\eta_{0}\|_{K}^{2}-\frac{1}{2}\sum_{i=1}^{n}\gamma_{i}E\|\eta_{i-1}\|_{2}^{2}+2u\sum_{i=1}^{n}\gamma_{i}J_{i}^{-2s}E\|\eta_{i-1}\|_{K}^{2}$
(102)
$\displaystyle\Rightarrow\frac{1}{2}\sum_{i=1}^{n}\gamma_{i}E\|\eta_{i-1}\|_{2}^{2}$
$\displaystyle\leq\|\eta_{0}\|_{K}^{2}+2u\sum_{i=1}^{n}\gamma_{i}J_{i}^{-2s}E\|\eta_{i-1}\|_{K}^{2}$
$\displaystyle\Rightarrow\frac{\gamma_{n}}{n}\sum_{i=1}^{n}E\|\eta_{i-1}\|_{2}^{2}$
$\displaystyle\leq
2\|\eta_{0}\|_{K}^{2}/n+4u\sum_{i=1}^{n}\gamma_{i}J_{i}^{-2s}E\|\eta_{i-1}\|_{K}^{2}/n$
$\displaystyle\stackrel{{\scriptstyle(1)}}{{\Rightarrow}}E\left\|\frac{1}{n}\sum_{i=1}^{n}\eta_{i-1}\right\|^{2}_{2}$
$\displaystyle\leq\frac{2\|f_{\rho}\|_{K}^{2}}{n\gamma_{n}}+\frac{4u\sum_{i=1}^{n}\gamma_{i}J_{i}^{-2s}E\|\eta_{i-1}\|_{K}^{2}}{n\gamma_{n}}$
$\displaystyle\Rightarrow\left(E[\|\bar{\eta}_{n}\|_{2}^{2}]\right)^{1/2}$
$\displaystyle\leq\left(\frac{2\|f_{\rho}\|_{K}^{2}+4u\sum_{i=1}^{n}\gamma_{i}J_{i}^{-2s}E\|\eta_{i-1}\|_{K}^{2}}{n\gamma_{n}}\right)^{1/2}$

Under assumptions A1-A3, we use Lemma D.2 to show that for $J_{i}\geq
i^{\alpha}\log^{2}i\vee 1$ for some $\alpha\geq\frac{1}{2s+1}$, then the
series $\sum_{i=1}^{n}\gamma_{i}J_{i}^{-2s}E\|\eta_{i-1}\|_{K}^{2}$ is
convergent and uniformly bounded, that is, it can be bounded by a constant
that does not depend on $n$. Recall that we chose
$\gamma_{n}=\gamma_{0}n^{-\frac{1}{2s+1}}$, so we conclude

$\left(E[\|\bar{\eta}_{n}\|_{2}^{2}]\right)^{1/2}=O\left(n^{-\frac{s}{2s+1}}\right)$
(103)

∎

###  D.4 Bound on noise sub-process

We remind the reader of the definition of our noise sub-process:

$\displaystyle\vartheta_{0}$ $\displaystyle=0$ (104)
$\displaystyle\vartheta_{n}$
$\displaystyle=\left(I-\gamma_{n}T_{X_{n},J_{n}}\right)\vartheta_{n-1}+\gamma_{n}\Xi_{n}$

where $\Xi_{n}=\left(Y_{n}-f_{\rho}\left(X_{n}\right)\right)K_{X_{n},J_{n}}$
(we also remind our reader $K_{X_{n},J_{n}}$ is the ”truncated kernel” at
level $J_{n}$). Also, $T_{X,J_{n}}$ and $T_{X_{n},J_{n}}$ are defined as:

$\displaystyle T_{X,J_{n}}(f)$ $\displaystyle=\int_{\mathcal{X}}\langle
f,K_{x,J_{n}}\rangle_{K}K_{x,J_{n}}d\rho_{X}(x)\quad\text{population
operator}$ (105) $\displaystyle T_{X_{n},J_{n}}(f)$ $\displaystyle=\langle
f,K_{X_{n},J_{n}}\rangle_{K}K_{X_{n},J_{n}}\quad\text{random operator}$

###### proof of bound on noise.

We need to define several sequences that are related to $\vartheta_{n}$ for
our technical analysis. The first sequence is:

$\displaystyle\eta_{0}^{0}$ $\displaystyle=0$ (106)
$\displaystyle\eta_{n}^{0}$
$\displaystyle=(I-\gamma_{n}T_{X,J_{n}})\eta_{n-1}^{0}+\gamma_{n}\Xi_{n}^{0}$

where
$\Xi_{n}^{0}=\Xi_{n}=\left(Y_{n}-f_{\rho}\left(X_{n}\right)\right)K_{X_{n},J_{n}}$.We
also define additional sequences, $\eta_{n}^{r}$, for each integer $r>0$, by

$\displaystyle\eta_{0}^{r}$ $\displaystyle=0$ (107)
$\displaystyle\eta_{n}^{r}$
$\displaystyle=(I-\gamma_{n}T_{X,J_{n}})\eta_{n-1}^{r}+\gamma_{n}\Xi_{n}^{r}$

where $\Xi_{n}^{r}=\left(T_{X,J_{n}}-T_{X_{n},J_{n}}\right)\eta_{n-1}^{r-1}$.
These sequences are easier to analyze than $\vartheta_{n}$ because the
operator in the recursive relationship $(I-\gamma_{n}T_{X,J_{n}})$ is a
deterministic (population) operator. In contrast the operator in the original
$\vartheta_{n}$ is random. We show in Lemma D.3 that as $r$ increases, the
amplitudes of ”noise” $\Xi_{n}^{r}$ get smaller. Additionally, using the fact
that all the sequences $\eta_{n}^{r}$ start with $\eta_{0}^{r}=0$, we can show
that $\eta_{n}^{r}$ becomes more concentrated about $0$ for larger $r$.

We split the noise process $\vartheta_{n}$ into two parts:

$\vartheta_{n}=\left(\vartheta_{n}-\sum_{k=0}^{r}\eta_{n}^{k}\right)+\sum_{k=0}^{r}\eta_{n}^{k}$
(108)

So its average satisfies

$\bar{\vartheta}_{n}=\left(\bar{\vartheta}_{n}-\sum_{k=0}^{r}\bar{\eta}_{n}^{k}\right)+\sum_{k=0}^{r}\bar{\eta}_{n}^{k}$
(109)

Here $\bar{\eta}_{n}^{k}$ is the averaged sequence of $\eta_{n}^{k}$ (i.e.
$\bar{\eta}_{n}^{k}=\frac{1}{n}\sum_{i=1}^{n}\eta_{i}^{k}$). Applying
Minkowski’s inequality gives us

$\left(E\|\bar{\vartheta}_{n}\|_{2}^{2}\right)^{1/2}\leq\sum_{k=0}^{r}\left(E\|\bar{\eta}_{n}^{k}\|_{2}^{2}\right)^{1/2}+\left(E\|\bar{\vartheta}_{n}-\sum_{k=0}^{r}\bar{\eta}_{n}^{k}\|_{2}^{2}\right)^{1/2}$
(110)

Now we define

$\alpha_{n}^{r}=\vartheta_{n}-\sum_{k=0}^{r}\eta_{n}^{k}$ (111)

We will show (Lemma D.7) that for $r\geq n$, we have
$\alpha_{n}^{r}=\eta_{n}^{r}=0$. We can see the second term in (110) is
exactly zero when we choose $r\geq n$, that is:

$\bar{\vartheta}_{n}-\sum_{k=0}^{r}\bar{\eta}_{n}^{k}=\frac{1}{n}\sum_{i=1}^{n}\left(\vartheta_{i}-\sum_{k=0}^{r}\eta_{i}^{k}\right)=\frac{1}{n}\sum_{i=1}^{n}\alpha_{i}^{r}=0$
(112)

Now our task is just bounding the first term
$\sum_{k=0}^{r}(E\|\bar{\eta}_{n}^{k}\|_{2}^{2})^{1/2}$, we analyze the
summation term by term. In Lemma D.4, we show that:

$E\left[\left\|\bar{\eta}_{n}^{k}\right\|_{2}^{2}\right]=O\left(\gamma_{0}^{k}R^{2k}C_{\epsilon}^{2}n^{-\frac{2s}{2s+1}}\right).$
(113)

To get this result, we first show in Lemma D.3 that the $\Xi_{j}^{k}$
variables are centered and satisfy some moment bounds. With these properties
in hand, we prove the above result in Lemma D.4, with the help of some
technical Lemma D.5 and D.6

Following (113), we have

$\displaystyle\sum_{k=0}^{r}\left(E\|\bar{\eta}_{n}^{k}\|_{2}^{2}\right)^{1/2}$
$\displaystyle\leq\sum_{k=0}^{r}C(\gamma_{0}R^{2})^{k/2}C_{\epsilon}n^{-\frac{s}{2s+1}}$
(114) $\displaystyle\leq
CC_{\epsilon}n^{-\frac{s}{2s+1}}\sum_{k=0}^{\infty}(\gamma_{0}R^{2})^{k/2}$
$\displaystyle=CC_{\epsilon}n^{-\frac{s}{2s+1}}\frac{1}{1-\sqrt{\gamma_{0}R^{2}}}=O\left(n^{-\frac{s}{2s+1}}\right)$

we note that $\gamma_{0}R^{2}<1$, which allows us to sum up the geometric
series.

Combining the results in (114) with $\alpha_{n}^{r}=0$ for $r\geq n$, we can
conclude from (110) that, for $r\geq n$:

$\left(E\left[\left\|\bar{\vartheta}_{n}\right\|_{2}^{2}\right]\right)^{1/2}\leq\sum_{k=0}^{r}\left(E\|\bar{\eta}_{n}^{k}\|_{2}^{2}\right)^{1/2}+0=O\left(n^{-\frac{s}{2s+1}}\right)$
(115)

∎

###  D.5 Combining the bounds

Plugging the final bounds (103) and (115) back into (95), we have the desired
result:

$E\left[\left\|\bar{f}_{n}-f_{\rho}\right\|_{2}^{2}\right]=O\left(n^{-\frac{2s}{2s+1}}\right)$
(116)

###  D.6 Technical Results

######  Lemma D.1.

There exists $R<\infty$, such that:

$\left\|K_{X_{n},J_{n}}\right\|_{K}^{2}\leq R^{2}$ (117)

for any $J_{n}\geq 1,X_{n}\in\mathcal{X}$.

###### Proof.

By the definition of $\|\cdot\|_{K}$ we have:

$\displaystyle\left\|K_{X_{n},J_{n}}\right\|_{K}^{2}$
$\displaystyle=\left\|\sum_{j=1}^{J_{n}}j^{-2s}\psi_{j}(X_{n})\psi_{j}\right\|_{K}^{2}$
(118)
$\displaystyle=\sum_{j=1}^{J_{n}}\frac{j^{-4s}\psi_{j}^{2}(X_{n})}{j^{-2s}}$
$\displaystyle\leq M^{2}\zeta(2s)=:R^{2}$

where $\zeta(\cdot)$ is the Riemann-zeta function. ∎

We use $R^{2}$ rather than $R$, in the bounds in this lemma because it
simplifies calculation where this lemma is applied

######  Lemma D.2.

Under A1-A3, and if we choose $J_{i}\geq i^{\alpha}\log^{2}i\vee 1$ for some
$\alpha\geq\frac{1}{2s+1}$, $\gamma_{i}=\gamma_{0}i^{-\frac{1}{2s+1}}$, then
there exists a number $C$ that does not depend on $n$ such that

$\sum_{i=1}^{n}\gamma_{i}J_{i}^{-2s}E\|\eta_{i-1}\|_{K}^{2}\leq C$ (119)

for all $n$.

###### Proof.

We first show the expectation of $\|\eta_{n}\|_{K}^{2}$ can be uniformly
bounded for $J_{n}$ that increases fast enough. Recall we have the following
recursive relationship (101):

$\displaystyle E\|\eta_{n}\|_{K}^{2}$ $\displaystyle\leq
E\|\eta_{n-1}\|_{K}^{2}-\frac{1}{2}\gamma_{n}E\|\eta_{n-1}\|_{2}^{2}+2u\gamma_{n}J_{n}^{-2s}E\|\eta_{n-1}\|_{K}^{2}$
(120) $\displaystyle\leq(1+2u\gamma_{n}J_{n}^{-2s})E\|\eta_{n-1}\|_{K}^{2}$
$\displaystyle\Rightarrow E\|\eta_{n}\|_{K}^{2}$
$\displaystyle\leq\prod_{i=1}^{n}(1+2u\gamma_{i}J_{i}^{-2s})\|\eta_{0}\|_{K}^{2}$

When we take $J_{i}\geq i^{\alpha}\log^{2}i\vee 1$, for some
$\alpha\geq\frac{1}{2s+1}$ we have
$\gamma_{i}J_{i}^{-2s}\leq\gamma_{0}\left((i\log^{2}i)^{-1}~{}\wedge~{}1\right)$.
Therefore $\prod_{i=1}^{n}(1+2u\gamma_{i}J_{i}^{-2s})$ converges as
$n\rightarrow\infty$:

$\displaystyle\prod_{i=1}^{n}(1+2u\gamma_{i}J_{i}^{-2s})$
$\displaystyle\leq\prod_{i=1}^{n}(1+C(i\log^{2}i)^{-1}\wedge 1)$ (121)
$\displaystyle=\exp\left(\log(\prod_{i=1}^{n}(1+C(i\log^{2}i)^{-1}\wedge
1)\right)$
$\displaystyle=\exp\left(\sum_{i=1}^{n}\log((1+C(i\log^{2}i)^{-1}\wedge
1)\right)$ $\displaystyle\leq\exp\left(\sum_{i=1}^{n}C(i\log^{2}i)^{-1}\wedge
1\right)$

We note that the series in the exponent of the last line converges. So at this
point we know, $E\|\eta_{n}\|_{K}^{2}$ can be uniformly bounded by
$K\|\eta_{0}\|_{K}=K\|f_{\rho}\|_{K}$, with a number $K$ that does not depend
on $n$. Now it direct to control the quantity of interest:

$\sum_{i=1}^{n}\gamma_{i}J_{i}^{-2s}E\|\eta_{i-1}\|_{K}^{2}\leq
K\|f_{\rho}\|\sum_{i=1}^{n}\gamma_{i}J_{i}^{-2s}<\infty,$ (122)

when $J_{i}\geq i^{\alpha}\log^{2}i\vee 1$ with $\alpha\geq\frac{1}{2s+1}$. ∎

######  Lemma D.3.

Assume A1-A3, for any integer $r,n\geq 0$ we have:

  * •

$\Xi_{n}^{r}$ is $\mathcal{F}_{n}$ measurable and
$\Xi_{n}^{r}\in\text{span}\\{\psi_{1},...,\psi_{J_{n}}\\}$,

  * •

$E[\Xi_{n}^{r}\ |\ \mathcal{F}_{n-1}]=0$, and

  * •

$E\left[\Xi_{n}^{r}\otimes\Xi_{n}^{r}\right]\preccurlyeq\gamma_{0}^{r}R^{2r}C_{\epsilon}^{2}T_{X,J_{n}}$.

Here $\mathcal{F}_{n}$ is the $\sigma$-algebra generated by
$(X_{i},Y_{i})_{i=1}^{n}$.

Note: because $\gamma_{0}R^{2}<1$ by our choice of $\gamma_{0}$, the upper
bound on $E\left[\Xi_{n}^{r}\otimes\Xi_{n}^{r}\right]$ is smaller for larger
$r$ (the third line in Lemma D.3).

###### Proof.

We remind our readers the definitions of $\Xi_{n}^{r}$ and $\eta_{n}^{r}$ are
given around (107). The first claim is direct. We first note that
$\eta^{0}_{n}$, $\Xi_{n}^{0}$ are both $\mathcal{F}_{n}$-measurable and belong
to $\text{span}\\{\psi_{1},...,\psi_{J_{n}}\\}$. Then we can show the
corresponding properties of $\eta_{n}^{r}$, $\Xi_{n}^{r}$ by induction.

For the second claim, we calculate directly:

$\displaystyle E[\Xi_{n}^{r}|\mathcal{F}_{n-1}]$
$\displaystyle=E[(T_{X,J_{n}}-T_{X_{n},J_{n}})\eta_{n-1}^{r-1}|\mathcal{F}_{n-1}]$
(123)
$\displaystyle=E[(T_{X,J_{n}}-T_{X_{n},J_{n}})|\mathcal{F}_{n-1}]\eta_{n-1}^{r-1}=0$

Now we show the last claim. We define

$D_{k}^{n}=(I-\gamma_{n}T_{X,J_{n}})(I-\gamma_{n-1}T_{X,J_{n-1}})\cdot\cdot\cdot(I-\gamma_{k}T_{X,J_{k}})=\prod_{i=k}^{n}\left(I-\gamma_{i}T_{X,J_{i}}\right).$
(124)

Note that each of the element in $D_{k}^{n}$ is self-adjoint, positive but
they in general do not commute. Because $T_{X,J_{i}}$ is a positive operator,
we have $I-\gamma_{i}T_{X,J_{i}}\preccurlyeq I$ for our choice of
$\gamma_{i}$. We are also going to use the following relationship in the rest
of this proof. Recall that we denote the adjoint of operator $A$ as $A^{*}$:

$\displaystyle\sum_{k=1}^{n}D_{k+1}^{n}\gamma_{k}^{2}T_{X,J_{k}}\left(D_{k+1}^{n}\right)^{*}$
$\displaystyle\preccurlyeq\gamma_{0}\sum_{k=1}^{n}D_{k+1}^{n}\gamma_{k}T_{X,J_{k}}\left(D_{k+1}^{n}\right)^{*}$
(125)
$\displaystyle\stackrel{{\scriptstyle(1)}}{{=}}\gamma_{0}\sum_{k=1}^{n}D_{k+1}^{n}\left(D_{k+1}^{n}\right)^{*}-D_{k}^{n}\left(D_{k+1}^{n}\right)^{*}$
$\displaystyle=\gamma_{0}\sum_{k=1}^{n}\left(I-\gamma_{n}T_{X,J_{n}}\right)\cdot\cdot\cdot\left(I-\gamma_{k+1}T_{X,J_{k+1}}\right)\left(I-\gamma_{k+1}T_{X,J_{k+1}}\right)\cdot\cdot\cdot\left(I-\gamma_{n}T_{X,J_{n}}\right)$
$\displaystyle\hskip 85.35826pt-D_{k}^{n}\left(D_{k+1}^{n}\right)^{*}$
$\displaystyle\stackrel{{\scriptstyle(2)}}{{\preccurlyeq}}\gamma_{0}\sum_{k=1}^{n}\left(I-\gamma_{n}T_{X,J_{n}}\right)\cdot\cdot\cdot\left(I-\gamma_{k+1}T_{X,J_{k+1}}\right)\cdot\cdot\cdot\left(I-\gamma_{n}T_{X,J_{n}}\right)-D_{k}^{n}\left(D_{k+1}^{n}\right)^{*}$
$\displaystyle=\gamma_{0}\sum_{k=1}^{n}D_{k+1}^{n}\left(D_{k+2}^{n}\right)^{*}-D_{k}^{n}\left(D_{k+1}^{n}\right)^{*}=\gamma_{0}\left(I-D_{1}^{n}\left(D_{2}^{n}\right)^{*}\right)\preccurlyeq\gamma_{0}I$

In (1) we used $D_{k+1}^{n}\gamma_{k}T_{X,J_{k}}=D_{k+1}^{n}-D_{k}^{n}$. In
(2) we used

$ABBA^{*}=AB^{1/2}BB^{1/2}A^{*}\preccurlyeq ABA$ (126)

for positive, self-adjoint $B\preccurlyeq I$.
Now we give an inductive argument (applying induction on $r$):
Initialisation: When $r=0$, recall that
$\Xi_{n}^{0}=(Y_{n}-f_{\rho}(X_{n}))K_{X_{n},J_{n}}$. Thus, we have that

$\displaystyle\langle
f,E\left[\Xi_{n}^{0}\otimes\Xi_{n}^{0}\right](f)\rangle_{K}$
$\displaystyle=\langle f,E\left[(Y_{n}-f_{\rho}(X_{n}))^{2}\langle
K_{X_{n},J_{n}},f\rangle_{K}K_{X_{n},J_{n}}\right]\rangle_{K}$ (127)
$\displaystyle\stackrel{{\scriptstyle(1)}}{{\leq}}\langle
f,C_{\epsilon}^{2}E\left[\langle
K_{X_{n},J_{n}},f\rangle_{K}K_{X_{n},J_{n}}\right]\rangle_{K}$
$\displaystyle=C_{\epsilon}^{2}\langle f,T_{X,J_{n}}f\rangle_{K}$

In (1) we used our noise assumption (A4). So we have

$E\left[\Xi_{n}^{0}\otimes\Xi_{n}^{0}\right]\preccurlyeq
C_{\epsilon}T_{X,J_{n}}$ (128)

To perform induction over $r$, we also need a bound on
$E\left[\eta_{n}^{0}\otimes\eta_{n}^{0}\right]$ as well. Recall that
$\eta_{n}^{0}=\sum_{k=1}^{n}\gamma_{k}D_{k+1}^{n}\Xi_{k}^{0}$ as defined in
(106). Thus we have:

$\displaystyle E\left[\eta_{n}^{0}\otimes\eta_{n}^{0}\right]$
$\displaystyle=\sum_{k=1}^{n}\sum_{j=1}^{n}\gamma_{j}\gamma_{k}D_{j+1}^{n}E\left[\Xi_{j}^{0}\otimes\Xi_{k}^{0}\right]\left(D_{k+1}^{n}\right)^{*}$
(129)
$\displaystyle\stackrel{{\scriptstyle(1)}}{{=}}\sum_{k=1}^{n}\gamma_{k}^{2}D_{k+1}^{n}E\left[\Xi_{k}^{0}\otimes\Xi_{k}^{0}\right]\left(D_{k+1}^{n}\right)^{*}$
$\displaystyle\preccurlyeq
C_{\epsilon}\sum_{k=1}^{n}D_{k+1}^{n}\gamma_{k}^{2}T_{X,J_{k}}\left(D_{k+1}^{n}\right)^{*}$
$\displaystyle\preccurlyeq C_{\epsilon}\gamma_{0}I$

In (1) the interaction terms vanish because the noise variables
$\Xi_{j}^{0},\Xi_{k}^{0}$ are mean-zero and independent when $j\neq k$.
Induction : If we assume for $r\geq 0$,

$E\left[\Xi_{n}^{r}\otimes\Xi_{n}^{r}\right]\preccurlyeq\gamma_{0}^{r}R^{2r}C_{\epsilon}T_{X,J_{n}}$
(130)

and

$E\left[\eta_{n}^{r}\otimes\eta_{n}^{r}\right]\preccurlyeq\gamma_{0}^{r+1}R^{2r}C_{\epsilon}I$
(131)

then for $r+1$:

$\displaystyle E\left[\Xi_{n}^{r+1}\otimes\Xi_{n}^{r+1}\right]$
$\displaystyle\stackrel{{\scriptstyle(1)}}{{=}}E\left[\left(T_{X,J_{n}}-T_{X_{n},J_{n}}\right)\eta_{n-1}^{r}\otimes\eta_{n-1}^{r}\left(T_{X,J_{n}}-T_{X_{n},J_{n}}\right)\right]$
(132)
$\displaystyle=E\left[\left(T_{X,J_{n}}-T_{X_{n},J_{n}}\right)E\left[\eta_{n-1}^{r}\otimes\eta_{n-1}^{r}\right]\left(T_{X,J_{n}}-T_{X_{n},J_{n}}\right)\right]$
$\displaystyle\preccurlyeq\gamma_{0}^{r+1}R^{2r}C_{\epsilon}E\left[\left(T_{X,J_{n}}-T_{X_{n},J_{n}}\right)^{2}\right]$
$\displaystyle=\gamma_{0}^{r+1}R^{2r}C_{\epsilon}\left(E\left[\left(T_{X_{n},J_{n}}\right)^{2}\right]-T_{X,J_{n}}^{2}\right)$
$\displaystyle\stackrel{{\scriptstyle(2)}}{{\preccurlyeq}}\gamma_{0}^{r+1}R^{2r+2}C_{\epsilon}T_{X,J_{n}}$

Here (1) is the definition of $\Xi_{n}^{r+1}$. For (2), it is sufficient to
show $E\left[\left(T_{X_{n},J_{n}}\right)^{2}\right]\preccurlyeq
R^{2}T_{X,J_{n}}$. This is true because:

$\langle
f,E\left[\left(T_{X_{n},J_{n}}\right)^{2}\right]f\rangle_{K}=E[\|T_{X_{n},J_{n}}(f)\|_{K}^{2}]=E[\langle
f,K_{X_{n},J_{n}}\rangle_{K}^{2}\|K_{X_{n},J_{n}}\|_{K}^{2}]\leq R^{2}\langle
f,T_{X,J_{n}}f\rangle_{K}$ (133)

Recall that $\eta_{n}^{r+1}=\sum_{k=1}^{n}D_{k+1}^{n}\gamma_{k}\Xi_{k}^{r+1}$,
then

$\displaystyle E\left[\eta_{n}^{r+1}\otimes\eta_{n}^{r+1}\right]$
$\displaystyle=E\left[\sum_{k=1}^{n}\gamma_{k}^{2}D_{k+1}^{n}\Xi_{k}^{r+1}\otimes\Xi_{k}^{r+1}\left(D_{k+1}^{n}\right)^{*}\right]$
(134)
$\displaystyle=\sum_{k=1}^{n}\gamma_{k}^{2}D_{k+1}^{n}E[\Xi_{k}^{r+1}\otimes\Xi_{k}^{r+1}]\left(D_{k+1}^{n}\right)^{*}$
$\displaystyle\preccurlyeq
C_{\epsilon}\gamma_{0}^{r+1}R^{2r}\sum_{k=1}^{n}D_{k+1}^{n}\gamma_{k}T_{X,J_{k}}\left(D_{k+1}^{n}\right)^{*}$
$\displaystyle\preccurlyeq C_{\epsilon}\gamma_{0}^{r+2}R^{2r+2}I$

∎

######  Lemma D.4.

Under assumptions A1-A3, we have

$E\left[\left\|\bar{\eta}_{n}^{r}\right\|_{2}^{2}\right]=O\left(\gamma_{0}^{r}R^{2r}C_{\epsilon}^{2}n^{-\frac{2s}{2s+1}}\right)$
(135)

###### Proof.

$\displaystyle n^{2}E\left[\left\|\bar{\eta}_{n}^{r}\right\|_{2}^{2}\right]$
$\displaystyle=E\left\|\sum_{j=1}^{n}\sum_{k=1}^{j}\left[\prod_{i=k+1}^{j}\left(I-\gamma_{i}T_{X,J_{i}}\right)\right]\gamma_{k}\Xi_{k}^{r}\right\|_{2}^{2}$
(136)
$\displaystyle=E\left\|\sum_{k=1}^{n}\sum_{j=k}^{n}\left[\prod_{i=k+1}^{j}\left(I-\gamma_{i}T_{X,J_{i}}\right)\right]\gamma_{k}\Xi_{k}^{r}\right\|_{2}^{2}$
$\displaystyle=\sum_{k=1}^{n}\gamma_{k}^{2}E\left\|\sum_{j=k}^{n}\left[\prod_{i=k+1}^{j}\left(I-\gamma_{i}T_{X,J_{i}}\right)\right]\Xi_{k}^{r}\right\|_{2}^{2}$
$\displaystyle=\sum_{k=1}^{n}\gamma_{k}^{2}E\left\langle\sum_{j=k}^{n}\left[\prod_{i=k+1}^{j}\left(I-\gamma_{i}T_{X,J_{i}}\right)\right]\Xi_{k}^{r},T_{X,J_{n}}\sum_{j=k}^{n}\left[\prod_{i=k+1}^{j}\left(I-\gamma_{i}T_{X,J_{i}}\right)\right]\Xi_{k}^{r}\right\rangle_{K}$
$\displaystyle=\sum_{k=1}^{n}\gamma_{k}^{2}E\operatorname{tr}\left(T_{X,J_{n}}M_{k}^{n}\Xi_{k}^{r}\otimes\Xi_{k}^{r}\left(M_{k}^{n}\right)^{*}\right)\quad\text{where
}M_{k}^{n}=\sum_{j=k}^{n}\left[\prod_{i=k+1}^{j}\left(I-\gamma_{i}T_{X,J_{i}}\right)\right]$