# Random Smoothing Regularization in Kernel Gradient Descent Learning

Liang Ding
Fudan University The authors’ names are sorted alphabetically. Corresponding
author: Wenjia Wang Tianyang Hu
Purdue University  Jiahang Jiang
The Hong Kong University of Science and Technology  Donghao Li
The Hong Kong University of Science and Technology  Wenjia Wang
The Hong Kong University of Science and Technology (Guangzhou)
and The Hong Kong University of Science and Technology  Yuan Yao
The Hong Kong University of Science and Technology

###### Abstract

Random smoothing data augmentation is a unique form of regularization that can
prevent overfitting by introducing noise to the input data, encouraging the
model to learn more generalized features. Despite its success in various
applications, there has been a lack of systematic study on the regularization
ability of random smoothing. In this paper, we aim to bridge this gap by
presenting a framework for random smoothing regularization that can adaptively
and effectively learn a wide range of ground truth functions belonging to the
classical Sobolev spaces. Specifically, we investigate two underlying function
spaces: the Sobolev space of low intrinsic dimension, which includes the
Sobolev space in $D$-dimensional Euclidean space or low-dimensional sub-
manifolds as special cases, and the mixed smooth Sobolev space with a tensor
structure. By using random smoothing regularization as novel convolution-based
smoothing kernels, we can attain optimal convergence rates in these cases
using a kernel gradient descent algorithm, either with early stopping or
weight decay. It is noteworthy that our estimator can adapt to the structural
assumptions of the underlying data and avoid the curse of dimensionality. This
is achieved through various choices of injected noise distributions such as
Gaussian, Laplace, or general polynomial noises, allowing for broad adaptation
to the aforementioned structural assumptions of the underlying data. The
convergence rate depends only on the effective dimension, which may be
significantly smaller than the actual data dimension. We conduct numerical
experiments on simulated data to validate our theoretical results.

##  1 Introduction

Random smoothing data augmentation is a technique used to improve the
generalization and robustness of machine learning models, particularly in the
context of deep learning. This method involves adding random noise, such as
Gaussian or Laplace noise, to the input data during the training process. The
idea behind random smoothing is to make the model more robust to small
perturbations in the input data, as the added noise simulates variations that
may occur naturally in real-world data. This augmentation approach has proven
to be an effective regularization technique, contributing to the empirical
success of deep learning models across various applications. For instance,
random flip, random crop, and color jitter can significantly improve the
classification accuracy in natural images (Goodfellow et al.,, 2016; Shorten
and Khoshgoftaar,, 2019). Random smoothing has been proven effective for
improving model robustness and generalization (Blum et al.,, 2020; Rosenfeld
et al.,, 2020; Mehra et al.,, 2021; Wang et al.,, 2020; Gao et al.,, 2020).
For example, random smoothing with Gaussian noise injection is introduced to
address the adversarial vulnerability (Cohen et al.,, 2019; Salman et al.,,
2019), and by encouraging the feature map to be invariant under data
augmentations, self-supervised contrastive learning methods (He et al.,, 2020;
Chen et al.,, 2020; Grill et al.,, 2020; Chen and He,, 2021; He et al.,, 2021)
can achieve state-of-the-art performance for various downstream tasks.

Random smoothing can be viewed as a form of regularization (Grandvalet et
al.,, 1997). Regularization techniques generally aim to reduce the complexity
of a model, making it less prone to fitting the noise in the training data
and, consequently, improving its performance on unseen data. Random smoothing
can be considered an implicit form of regularization, as it does not directly
modify the model’s parameters or loss function, unlike explicit regularization
techniques such as $\ell_{1}$ or $\ell_{2}$ regularization. Instead, it
indirectly influences the model’s behavior by altering the input data during
training. By adding random noise to the input data, random smoothing forces
the model to focus on the underlying structure of the data rather than
memorizing specific instances. This leads to more robust and generalizable
models that can better handle variations in real-world data. As a result,
random smoothing acts as a regularizer, improving the model’s ability to
generalize from the training set to unseen data. Such a regularization
perspective at least starts with Grandvalet et al., (1997). However, in spite
of the empirical success of random smoothing in various applications, there is
a lack of systematic research on the regularization effect of random smoothing
in the literature.

In this paper, we address this gap by examining the classic nonparametric
regression problem from the perspective of random smoothing regularization. In
nonparametric regression, the primary objective is to uncover the functional
relationship between input and output variables. By making appropriate
assumptions about the underlying truth function and selecting the appropriate
estimator, we focus on understanding the efficiency of the estimation,
specifically, the rate at which the estimation error converges to zero as the
sample size $n$ increases. The optimal convergence rate is typically dictated
by the problem’s inherent complexity. The actual achievable convergence rates
depend on the specific estimation methods employed. Among various techniques,
we consider kernel methods that have been extensively investigated in the
research literature (Wahba,, 1990; Hastie et al.,, 2001).

In this study, we present a unified framework that can learn a wide range of
$D$-dimensional ground truth functions belonging to the classical Sobolev
spaces ($\mathcal{W}^{m_{f}}$) in an effective and adaptive manner. The
framework incorporates random smoothing as a central component. Our hypothesis
space is a reproducing kernel Hilbert space that is associated with a kernel
function of smoothness denoted by $m_{0}$. Random smoothing regularization
leads to a novel convolution between the kernel function and a probability
density function for the injected input noise. This injected noise is governed
by either short or long-tail distributions, namely Gaussian and polynomial
(including Laplace) noises, respectively. The resulting convolution-based
random smoothing kernel enables us to adapt to the smoothness of the target
functions more efficiently. Notably, we establish that for any $m_{0}$ and
$m_{f}$ greater than $D/2$, optimal convergence rates can be achieved by
utilizing random smoothing regularization and appropriate early stopping
and/or weight decay techniques.

To be specific, we investigate two possible function spaces that may contain
the target function. In Section 4.2, we analyze the Sobolev space with a low
intrinsic dimension, which is denoted by $d$. This space covers both
$D$-dimensional Euclidean spaces (when $d=D$) and low-dimensional sub-
manifolds as specific examples. In Section 4.3, we explore the mixed smooth
Sobolev spaces, which possess a tensor structure. Our principal findings are
summarized below.

  * •

In case of Sobolev space of low intrinsic dimensionality $d\leq D$:

When using Gaussian random smoothing, an upper bound of the convergence rate
is achieved at $n^{-m_{f}/(2m_{f}+d)}(\log n)^{D+1}$, which recovers the
results presented in Hamm and Steinwart, 2021a  and is hypothetically optimal
up to a logarithmic factor. However, in contrast to Hamm and Steinwart, 2021a
, we present a different approach that allows us to analyze polynomial
smoothing;

When using polynomial random smoothing with data size adaptive smoothing
degree, a convergence rate of $n^{-m_{f}/(2m_{f}+d)}(\log n)^{2m_{f}+1}$ is
achieved, which is again, hypothetically optimal up to a logarithmic factor.

  * •

In case of mixed smooth Sobolev spaces, using polynomial random smoothing of
degree $m_{\varepsilon}$, a fast convergence rate of
$n^{-2m_{f}/(2m_{f}+1)}(\log
n)^{\frac{2m_{f}}{2m_{f}+1}\left(D-1+\frac{1}{2(m_{0}+m_{\varepsilon})}\right)}$
is achieved, which is optimal up to a logarithmic factor.

To the best of our knowledge, such results have not been studied in the
literature so far. They have various implications below.

First of all, these results enhance the convergence rates in the context of
kernel ridge regression by incorporating random smoothing data augmentation
with two other popular techniques, early stopping and weight decay. In kernel
ridge regression, it is crucial to balance the smoothness of the kernel
function ($m_{0}$) with that of the ground truth ($m_{f}$). In practice, it is
common for $m_{0}$ to be unequal to $m_{f}$. In cases of mismatch,
regularization becomes essential. Specifically, if $m_{0}\in[m_{f}/2,\infty)$,
the optimal convergence rate $n^{-m_{f}/(2m_{f}+D)}$ can be achieved by
employing an appropriate ridge penalty strength. This result can be
generalized to low intrinsic dimensionality $d\leq D$, where the
hypothetically optimal convergence rate is $n^{-m_{f}/(2m_{f}+d)}$ (Hamm and
Steinwart, 2021a, ). However, when the chosen kernel has a smoothness $m_{0}$
less than $m_{f}/2$, the optimal adaptation is not well studied in kernel
ridge regression. In contrast, our findings demonstrate optimal adaptation for
arbitrary $m_{0}$ and $m_{f}\geq D/2$ without such a constraint. This
highlights the broad adaptation ability of random smoothing regularization.

Moreover, the optimal adaptation of polynomial random smoothing has an
implication for neural networks via the (generalized) Laplace random
smoothing. It is known that the training of neural networks, with enough
overparametrization, can be characterized by kernel methods with a special
family of kernels called the “neural tangent kernel” (NTK). Due to the low
smoothness of the ReLU activation function, the corresponding NTK also has a
low smoothness that is the same as a Laplace kernel (Chen and Xu,, 2020;
Geifman et al.,, 2020). To the best of our knowledge, the estimation error is
at the rate $n^{-\frac{D}{2D-1}}$ (Hu et al.,, 2021). Our results, using the
polynomial random smoothing with (generalized) Laplace distributions, show
that the convergence rate can be improved, which sheds light on understanding
non-smooth augmentations such as random crop and mask. Based on this
understanding, numerical experiments with neural networks are conducted on
simulated data to corroborate our theoretical results.

Finally, it is worth mentioning that with random smoothing, the convergence
rates mentioned above can be obtained by early stopping. However, if one
applies weight decay, the number of iterations can be reduced from
polynomial$(n)$ to polynomial$(\log n)$. Additionally, our estimator can adapt
to the low-dimensional assumptions mentioned earlier, as the convergence rates
depend on $D$ at most logarithmically, alleviating the curse of
dimensionality. It is also important to note that we do not employ the
spectrum of integral operator technique (Yao et al.,, 2007; Lin et al.,, 2016;
Lin and Rosasco,, 2017), but instead use Fourier analysis, which provides a
universal basis for kernels of different smoothness, and avoids imposing
conditions on the eigenvalues and eigenfunctions of the kernel function. This
is because there is no clear relationship between the low intrinsic dimension
and the eigenvalues of the integral operator. Furthermore, our theoretical
analysis can be applied to the widely used Matérn kernel functions.

The remainder of this paper is structured as follows. In Section 2, we provide
a review of related works. Section 3 introduces the settings considered in
this work, which include early stopping with a random smoothing kernel, as
well as the conditions and assumptions utilized in this work. The main
theoretical results are presented in Section 4, and numerical studies are
conducted in Section 5. Conclusions and a discussion are provided in Section
6. Technical proofs are included in the Appendix.

##  2 Related Works

Various means of regularization have been proposed for kernel methods to
better recover the underlying function, among which, ridge penalty and early
stopping are the most popular. Kernel ridge regression has been extensively
studied in the literature, see Blanchard and Mücke, (2018); Dicker et al.,
(2017); Guo et al., (2017); Lin et al., (2017); Steinwart et al., (2009); Tuo
et al., (2020); Wu et al., (2006) for example. Early stopping treats the
number of training iterations as a hyperparameter in the optimization process,
which has been extensively studied by the applied mathematics community
(Dieuleveut and Bach,, 2016; Yao et al.,, 2007; Pillaud-Vivien et al.,, 2018;
Raskutti et al.,, 2014). Various forms of early stopping also have been
studied including boosting (Zhang and Yu,, 2005; Bartlett and Traskin,, 2007),
conjugate gradient algorithm (Blanchard and Krämer,, 2016) and kernel gradient
descent (Bühlmann and Yu,, 2002; Caponnetto and Yao,, 2006; Yao et al.,, 2007;
Wei et al.,, 2017; Lin et al.,, 2016). Some works (e.g. Lin et al., (2016);
Lin and Rosasco, (2017); Pillaud-Vivien et al., (2018)) have explored early
stopping by employing the integral operator induced by the kernel, imposing
conditions on the eigenvalues and eigenfunctions of the kernel function.
Smoothness or regularity of functions thus implicitly depends on the measure
that defines the spectrum of the integral operator, whereas classical
smoothness like Sobolev spaces is not explicitly handled.

In kernel regression with gradient descent, Raskutti et al., (2014) showed
that early stopping and ridge penalty both can achieve the optimal convergence
rate if the smoothness is well-specified. Yet, kernel ridge regression might
suffer the “saturation issues” while early stopping does not (Engl et al.,,
1996; Yao et al.,, 2007). In regression problems, it is usually assumed that
the domain of interest has a positive Lebesgue measure, while in practice, the
data generating distribution is supported on some low-dimensional smooth sub-
manifold (Scott and Nowak,, 2006; Yang and Dunson,, 2016; Ye and Zhou,, 2008,
2009; Hamm and Steinwart, 2021b, ; Hamm and Steinwart, 2021a, ). Kernel
methods can circumvent the curse of dimensionality and adapt to various low-
dimensional assumptions of the underlying function. In particular, Hamm and
Steinwart, 2021b ; Hamm and Steinwart, 2021a  generalized the manifold
assumption by applying the box-counting dimension of the support of the data
distribution, and derived upper bounds on the convergence rate of the
prediction error. Another simplifying assumption is tensor product kernels
(Gretton,, 2015; Szabó and Sriperumbudur,, 2017), whose product forms allow
efficient computation of Gaussian process regression (Saatçi,, 2012; Wilson
and Nickisch,, 2015; Ding and Zhang,, 2022; Chen et al.,, 2022) and analysis
of independent component (Bach and Jordan,, 2002; Gretton et al.,, 2005,
2007). The RKHS induced by a tensor product kernel is simply tensored RKHS
(Paulsen and Raghupathi,, 2016). Tensor product kernels we consider induce the
tensored Sobolev spaces (Dung et al.,, 2018).

For complicated high-dimensional data, deep learning models seem to perform
extremely well, which has sparked numerous investigations into their
generalization ability. As it turns out, the training of neural networks has
deep connections to kernel methods with neural tangent kernels (NTK). Under
proper initialization, training sufficiently wide DNN with gradient descent
equates to kernel regression using NTK. First introduced by Jacot et al.,
(2018), the correspondence has been significantly extended (Du et al.,, 2018;
Li and Liang,, 2018; Arora et al., 2019a, ; Cao and Gu,, 2020; Arora et al.,
2019b, ; Li et al.,, 2019; Huang et al.,, 2020; Kanoh and Sugiyama,, 2021; Hu
et al.,, 2022). From the NTK point of view, ridge penalty and early stopping
are also vital in training neural networks. The former is equivalent to weight
decay (Hu et al.,, 2021), which is applied by default in training deep
learning models for better generalization, so is early stopping (Prechelt,,
1998). Zhang et al., (2021); Hardt et al., (2016) revealed that longer
training can harm the generalization performance of deep models. Li et al.,
(2020); Bai et al., (2021) utilized early stopping to improve robustness to
label noises.

Besides NTK, various data augmentation techniques in deep learning that are
proven effective in improving model generalization can also provide
inspiration for kernel methods. Grandvalet et al., (1997) studied from a
regularization perspective how noise injection can improve generalization.
Data augmentation is particularly important for handling natural images
(Shorten and Khoshgoftaar,, 2019), where horizontal flip, random crop, color
jitter can significantly improve the classification accuracy. By applying the
above augmentations, self-supervised contrastive learning methods (He et al.,,
2020; Chen et al.,, 2020; Grill et al.,, 2020; Chen and He,, 2021; He et al.,,
2021) can achieve state-of-the-art performance for various downstream tasks.
Randomized smoothing (Cohen et al.,, 2019; Salman et al.,, 2019) is a special
data augmentation, first proposed to address the adversarial vulnerability
(Goodfellow et al.,, 2014; Carlini and Wagner,, 2017) of deep learning models.
The key idea is to perturb the input with random noise injection and make
predictions by aggregating the outputs from all augmented inputs. Random
smoothing has been proven effective for improving model robustness and
generalization (Rosenfeld et al.,, 2020; Mehra et al.,, 2021; Wang et al.,,
2020; Gao et al.,, 2020). Our proposed framework incorporates random
smoothing, together with weight decay and early stopping, to provide a unified
solution for the smoothness mismatch problem in kernel regression. It is worth
clarifying the difference between our method and the “errors in variables”
literature (Zhou et al.,, 2019; Wang et al.,, 2022; Cressie and Kornak,, 2003;
Cervone and Pillai,, 2015). Though the formulations seem similar, i.e., the
inputs in both cases are corrupted with noises, the two are fundamentally
different. In our setting, both the input $\bm{x}$ and added noise
$\bm{\varepsilon}$ are known (we control the noises in our estimator) while in
the other setting, the input is noisy and only $\bm{x}+\bm{\varepsilon}$ is
observed.

##  3 Random Smoothing Kernel Regression

In this section, we introduce the problem of interest, our methodology, and
the necessary conditions used in this work.

###  3.1 Problem Setting

Suppose we have observed data $(\bm{x}_{j},y_{j})$ for $j=1,...,n$, which
follows the relationship given by

$\displaystyle y_{j}=f^{*}(\bm{x}_{j})+\epsilon_{j}.$ (1)

Here, $\bm{x}_{j}$’s are independent and identically distributed (i.i.d.)
following a marginal distribution $P_{\mathbf{X}}$ with support
supp$(P_{\mathbf{X}})=\Omega\subset\mathbb{R}^{D}$. The function
$f^{*}\in\mathcal{H}(\Omega)$, where $\mathcal{H}(\Omega)$ denotes a function
space, and $\epsilon_{j}$’s are i.i.d. noise variables with mean zero and
finite variance. Our objective is to recover the function $f^{*}$ based on the
noisy observations.

In this work, we consider two cases. In the first case (Section 4.2), the
function space $\mathcal{H}(\Omega)$ is a Sobolev space with smoothness $m$,
denoted by $\mathcal{W}^{m}(\Omega)$, and the data is of low intrinsic
dimension. In the second case (Section 4.3), the function space
$\mathcal{H}(\Omega)$ is a tensor Sobolev space. Throughout this work, we
assume without loss of generality that $P_{\mathbf{X}}$ follows a uniform
distribution. Note that our theoretical analysis can be easily extended to the
case where $P_{\mathbf{X}}$ is upper and lower bounded by positive constants.

In order to recover the function $f^{*}$, we use reproducing kernel Hilbert
spaces (RKHSs). We briefly introduce the RKHSs and their relationship with
Sobolev spaces in the following, and refer to Wendland, (2004) and Adams and
Fournier, (2003) for details. Let $K:\Omega\times\Omega\rightarrow\mathbb{R}$
be a symmetric positive definite kernel function. Define the linear space

$\displaystyle
F_{K}(\Omega)=\left\\{\sum_{k=1}^{n}\beta_{k}K(\cdot,\bm{x}_{k}):\beta_{k}\in\mathbb{R},\bm{x}_{k}\in\Omega,n\in\mathbb{N}\right\\},$
(2)

and equip this space with the bilinear form

$\displaystyle\left\langle\sum_{k=1}^{n}\beta_{k}K(\cdot,\bm{x}_{k}),\sum_{j=1}^{m}\gamma_{j}K(\cdot,\bm{x}^{\prime}_{j})\right\rangle_{K}:=\sum_{k=1}^{n}\sum_{j=1}^{m}\beta_{k}\gamma_{j}K(\bm{x}_{k},\bm{x}^{\prime}_{j}).$

Then the reproducing kernel Hilbert space $\mathcal{H}_{K}(\Omega)$ generated
by the kernel function $K$ is defined as the closure of $F_{K}(\Omega)$ under
the inner product $\langle\cdot,\cdot\rangle_{K}$, and the norm of
$\mathcal{H}_{K}(\Omega)$ is $\|f\|_{\mathcal{H}_{K}(\Omega)}=\sqrt{\langle
f,f\rangle_{\mathcal{H}_{K}(\Omega)}}$, where
$\langle\cdot,\cdot\rangle_{\mathcal{H}_{K}(\Omega)}$ is induced by
$\langle\cdot,\cdot\rangle_{K}$. The following theorem gives another
characterization of the reproducing kernel Hilbert space when $K$ is
stationary, via the Fourier transform. Our notion of the Fourier transform is

$\displaystyle\mathcal{F}(g)(\bm{\omega})=(2\pi)^{-D/2}\int_{\mathbb{R}^{D}}g(\bm{x})e^{-i\bm{\omega}^{T}\bm{x}}{\rm
d}\bm{x},$

for a function $g\in L_{1}(\mathbb{R}^{D})$. Note that a kernel function $K$
is said to be stationary if the value $K(\bm{x},\bm{x}^{\prime})$ only depends
on the difference $\bm{x}-\bm{x}^{\prime}$. Thus, we can write
$K(\bm{x}-\bm{x}^{\prime}):=K(\bm{x},\bm{x}^{\prime})$.

######  Theorem 3.1 (Theorem 10.12 of Wendland, (2004)).

Let $K$ be a positive definite kernel function that is stationary, continuous,
and integrable in $\mathbb{R}^{D}$. Define

$\mathcal{G}:=\\{f\in L_{2}(\mathbb{R}^{D})\cap
C(\mathbb{R}^{D}):\mathcal{F}(f)/\sqrt{\mathcal{F}(K)}\in
L_{2}(\mathbb{R}^{D})\\},$

with the inner product

$\langle
f,g\rangle_{\mathcal{H}_{K}(\mathbb{R}^{D})}=(2\pi)^{-d/2}\int_{\mathbb{R}^{d}}\frac{\mathcal{F}(f)(\bm{\omega})\overline{\mathcal{F}(g)(\bm{\omega})}}{\mathcal{F}(K)(\bm{\omega})}{\rm
d}\bm{\omega}.$

Then $\mathcal{G}=\mathcal{H}_{K}(\mathbb{R}^{D})$, and both inner products
coincide.

For $m>D/2$, the (fractional) Sobolev norm for function $g$ on
$\mathbb{R}^{D}$ is defined by

$\displaystyle\|g\|^{2}_{\mathcal{W}^{m}(\mathbb{R}^{D})}=\int_{\mathbb{R}^{d}}|\mathcal{F}(g)(\bm{\omega})|^{2}(1+\|\bm{\omega}\|_{2}^{2})^{m}{\rm
d}\bm{\omega},$ (3)

and the inner product of a Sobolev space $\mathcal{W}^{m}(\mathbb{R}^{D})$ is
defined by

$\langle
f,g\rangle_{\mathcal{W}^{m}(\mathbb{R}^{D})}=\int_{\mathbb{R}^{D}}\mathcal{F}(f)(\bm{\omega})\overline{\mathcal{F}(g)(\bm{\omega})}(1+\|\omega\|_{2}^{2})^{m}{\rm
d}\bm{\omega}.$

######  Remark 3.1.

In this work, we are only interested in Sobolev spaces with $m>D/2$ because
these spaces contain only continuous functions according to the Sobolev
embedding theorem.

It can be shown that if $m$ is an integer, the norm defined in (3) is
equivalent to that of the usual Sobolev space (Adams and Fournier,, 2003). If
$m$ is not an integer, then the corresponding Sobolev space is called a Bessel
potential space (Almeida and Samko,, 2006; Gurka et al.,, 2007). The Sobolev
space on a region $\tilde{\Omega}$ with a positive Lebesgue measure can be
defined via restrictions as

$\displaystyle\|f\|_{\mathcal{W}^{m}(\tilde{\Omega})}=\inf\\{\|f_{E}\|_{\mathcal{W}^{m}(\mathbb{R}^{D})}:f_{E}\in\mathcal{W}^{m}(\mathbb{R}^{D}),f_{E}|_{\tilde{\Omega}}=f\\},$

where $f_{E}|_{\tilde{\Omega}}$ denotes the restriction of $f_{E}$ to
$\tilde{\Omega}$.

Comparing Theorem 3.1 and (3), it can be seen that if

$\displaystyle
c_{1}(1+\|\bm{\omega}\|_{2}^{2})^{-m}\leq\mathcal{F}(K)(\bm{\omega})\leq
c_{2}(1+\|\bm{\omega}\|_{2}^{2})^{-m},\forall\bm{\omega}\in\mathbb{R}^{D},$

for some two constants $c_{1},c_{2}>0$, then $\mathcal{W}^{m}(\mathbb{R}^{D})$
coincides with the reproducing kernel Hilbert space
$\mathcal{H}_{K}(\mathbb{R}^{D})$ with equivalent norms (also see Wendland,
(2004), Corollary 10.13). By the extension theorem (DeVore and Sharpley,,
1993), $\mathcal{H}_{K}(\Omega)$ also coincides with
$\mathcal{W}^{m}(\Omega)$, and two norms are equivalent.

###  3.2 Random Smoothing Kernel Regression with Early Stopping

In this study, we systematically investigate the efficiency of random
smoothing data augmentation, which is a widely used technique in deep
learning, in improving the estimation efficiency (i.e., convergence rate) for
$f^{*}\in\mathcal{H}(\Omega)$ without assuming any relationship between
$\mathcal{H}(\Omega)$ and $\mathcal{H}_{K}(\Omega)$ and considering a wide
context of $\Omega$ that may have Lebesgue measure zero. To overcome the lack
of smoothness in $\mathcal{H}_{K}(\Omega)$, we construct $N$ augmentations for
each observed input point $\bm{x}_{j}$ by adding i.i.d. noise
$\bm{\varepsilon}_{jk}$ with a continuous probability density function
$p_{\varepsilon}$. We can generate $\bm{\varepsilon}_{jk}$ independently for
each $j$, or we can generate $\bm{\varepsilon}_{k}$ for $k=1,...,N$, and apply
them to all $\bm{x}_{j}$, $j=1,...,n$ simultaneously. While the latter is
easier to implement, the former is easier to theoretically justify. Due to its
lower computational complexity, we only consider the latter method in this
work.

######  Remark 3.2 (Adding non-smooth noise and practical data augmentation
techniques).

It should be noted that we do not assume $p_{\varepsilon}$ to be Gaussian, and
can be non-smooth. While applying Gaussian noise is a common practice, not all
data augmentation techniques involve smooth noise, such as random crop, random
mask, and random flip. In this work, we investigate various types of noise,
including non-smooth Laplace noise and smooth Gaussian noise. Although adding
non-smooth noise still cannot capture the effects of complex data augmentation
techniques such as random mask or random crop, we aim to use it as a tool to
gain insights into the success of these more complicated data augmentations.

With augmented data, we proceed to the estimation of the function $f^{*}$. For
any point $\bm{x}\in\Omega$, we obtain the estimator by computing the average
of the function values evaluated at the $N$ augmented inputs. Specifically,
the estimator is constructed as

$\displaystyle
f(\bm{x})=\frac{1}{N}\sum_{k=1}^{N}h(\bm{x}+\bm{\varepsilon}_{k})$ (4)

for $h\in\mathcal{H}_{K}(\Omega)$. By properties of the RKHS, $f$ as in (4) is
also inside $\mathcal{H}_{K}(\Omega)$. We consider the following $l_{2}$ loss
function defined as

$\displaystyle
L_{n}(f)=\frac{1}{2n}\sum_{j=1}^{n}\left(f(\bm{x}_{j})-y_{j}\right)^{2},$ (5)

or equivalently,

$\displaystyle
L_{n}(h)=\frac{1}{2n}\sum_{j=1}^{n}\left(\frac{1}{N}\sum_{k=1}^{N}h(\bm{x}_{j}+\bm{\varepsilon}_{k})-y_{j}\right)^{2}.$

######  Remark 3.3.

The loss function $L_{n}(h)$ is slightly different from the loss function used
in practice, i.e.,

$\displaystyle
L_{n}^{\prime}(h)=\frac{1}{2n}\sum_{j=1}^{n}\frac{1}{N}\sum_{k=1}^{N}\left(h(\bm{x}_{j}+\bm{\varepsilon}_{k})-y_{j}\right)^{2}.$

However, it can be shown that $L_{n}(h)$ is close to $L_{n}^{\prime}(h)$. To
see this, note that

$\displaystyle L_{n}^{\prime}(h)-L_{n}(h)=$
$\displaystyle\frac{1}{2n}\sum_{j=1}^{n}\frac{1}{2N^{2}}\sum_{k=1}^{N}\sum_{l=1}^{N}\left(h(\bm{x}_{j}+\bm{\varepsilon}_{k})-h(\bm{x}_{j}+\bm{\varepsilon}_{l})\right)^{2}.$
(6)

As we will see later in Section 4, we require that the variance of
$\bm{\varepsilon}_{k}$ to converge to zero, which implies that the right-hand
side in (6) is close to zero.

In order to minimize (5), we apply the gradient descent method. Since we
impose a restriction that the estimator $f$ is in the RKHS
$\mathcal{H}_{K}(\Omega)$, by the representer theorem, it suffices to consider
the function space

$\displaystyle\mathcal{F}_{0}=\left\\{f:f(\cdot)=\sum_{j=1}^{n}\sum_{k=1}^{N}w_{jk}K(\cdot-(\bm{x}_{j}+\bm{\varepsilon}_{k})),w_{jk}\in\mathbb{R}\right\\}.$

Because the number of parameters in $\mathcal{F}_{0}$ scales as $n\times N$,
which can be prohibitively large if there are too many augmentations, it is
often necessary to reduce the flexibility of $\mathcal{F}_{0}$ in order to
minimize the loss function (5). To achieve this, we consider a subspace of
$\mathcal{F}_{0}$, denoted by

$\displaystyle\mathcal{F}=\left\\{f:f(\cdot)=\sum_{j=1}^{n}\sum_{k=1}^{N}w_{j}K(\cdot-(\bm{x}_{j}+\bm{\varepsilon}_{k})),w_{j}\in\mathbb{R}\right\\},$

i.e., all the weights for the different augmented data from the same input
$\bm{x}_{j}$ are the same. Define an empirical random smoothing kernel
function by

$\displaystyle
K_{S}(\bm{x}_{l}-\bm{x}_{j}):=\frac{1}{N^{2}}\sum_{k_{1}=1}^{N}\sum_{k_{2}=1}^{N}K(\bm{x}_{l}+\bm{\varepsilon}_{k_{1}}-(\bm{x}_{j}+\bm{\varepsilon}_{k_{2}})),$
(7)

whose expectation leads to the following random smoothing kernel function,
which plays an important role in the convergence analysis.

######  Definition 3.1 (Random smoothing kernel function).

The kernel function $K_{S}$ defined in (7) is the empirical random smoothing
kernel function corresponding to the original kernel $K$. The expectation of
$K_{S}$ with respect to the noise $\bm{\varepsilon}_{k}$ is the convoluted
kernel function $K*p_{\varepsilon}$, where $*$ is a convolution operator
defined by

$\displaystyle(g_{1}*g_{2})(\bm{s})=\int g_{1}(\bm{t})g_{2}(\bm{s}-\bm{t}){\rm
d}\bm{t},$

for two functions $g_{1}$ and $g_{2}$. We call the convoluted kernel function
$K*p_{\varepsilon}$ as the random smoothing kernel function.

Now we can rewrite the loss function $L_{n}(f)$ in (5) as

$\displaystyle
L_{n}(\bm{w})=\frac{1}{2}\left\|\bm{y}-\mathbf{K}\bm{w}\right\|_{2}^{2},$ (8)

where $\mathbf{K}=(K_{S}(\bm{x}_{j}-\bm{x}_{k}))_{jk}$,
$\bm{w}=(w_{1},...,w_{n})^{T}$, and $\bm{y}=(y_{1},...,y_{n})^{T}$. As stated
in Raskutti et al., (2014), it is more natural to perform gradient descent on
the transformed vector $\bm{\theta}=\sqrt{\mathbf{K}}\bm{w}$, where the square
root can be taken because $\mathbf{K}$ is positive (semi-)definite. Then, we
apply gradient descent on the square loss (8) with the transformed vector
$\bm{\theta}$. Initialize $\bm{\theta}_{0}=\bm{w}_{0}=0$. Taking gradient with
respect to $\bm{\theta}$, direct computation shows that the gradient update is

$\displaystyle\bm{\theta}_{t+1}=\bm{\theta}_{t}-\beta_{t}\left(\mathbf{K}\bm{\theta}_{t}-\sqrt{\mathbf{K}}\bm{y}\right),$
(9)

where $\beta_{t}>0$, $t=0,1,2,\ldots$ is the learning rate (step size). With
parameter $\bm{w}_{t}$ obtained at the $t$-th iteration, the corresponding
estimator of $f^{*}(\bm{x})$ for any point $\bm{x}\in\Omega$ is defined by

$\displaystyle f_{t}(\bm{x})=\bm{w}_{t}^{T}\mathbf{k}(\bm{x}),$ (10)

where
$\mathbf{k}(\bm{x})=(K_{S}(\bm{x}-\bm{x}_{1}),\ldots,K_{S}(\bm{x}-\bm{x}_{n}))^{T}$.

In practice, gradient descent is often paired with weight decay (Krogh and
Hertz,, 1992) to prevent overfitting and improve generalization (Hu et al.,,
2021). Therefore, we also consider the gradient descent with weight decay,
where the parameter $\bm{\theta}$ is updated by

$\displaystyle\bm{\theta}_{t+1}=\bm{\theta}_{t}-\beta_{t}\left(\mathbf{K}\bm{\theta}_{t}-\sqrt{\mathbf{K}}\bm{y}\right)-\alpha_{t}\bm{\theta}_{t},$
(11)

with $\alpha_{t}>0$, $t=0,1,2,\ldots$ being the strength of weight decay. The
learning rate $\beta_{t}$ and weight decay parameter $\alpha_{t}$ can be
varied with $t$, but for mathematical convenience, we assume that the step
sizes $\beta_{t}$ and the weights decay parameter $\alpha_{t}$ are not related
to the iteration number $t$, i.e., $\beta_{t}=\beta$ and $\alpha_{t}=\alpha$
for all $t=0,1,2,\ldots$.

In this work, we are interested in the prediction error

$\displaystyle\|f^{*}-f_{t}\|_{L_{2}(P_{\mathbf{X}})}.$ (12)

In the rest of this paper, the following definitions are used. For two
positive sequences $a_{n}$ and $b_{n}$, we write $a_{n}\asymp b_{n}$ if, for
some $C,C^{\prime}>0$, $C\leq a_{n}/b_{n}\leq C^{\prime}$. Similarly, we write
$a_{n}\gtrsim b_{n}$ if $a_{n}\geq Cb_{n}$ for some constant $C>0$, and
$a_{n}\lesssim b_{n}$ if $a_{n}\leq C^{\prime}b_{n}$ for some constant
$C^{\prime}>0$. Also, $C,C^{\prime},c_{j},C_{j},j\geq 0$ are generic positive
constants, of which value can change from line to line.

##  4 Main Results

In this section, we present our main theoretical results. We begin by
collecting all the assumptions that will be used throughout the paper in
Section 4.1. Then, in Section 4.2, we consider the case where $\Omega$ has a
finite intrinsic dimension. Finally, in Section 4.3, we consider the case
where $\mathcal{H}(\Omega)$ is a tensor RKHS.

###  4.1 Assumptions

In this work, we will use the following assumptions.

######  Assumption 4.1.

The error $\epsilon_{j}$’s in (1) are i.i.d. sub-Gaussian (van de Geer,,
2000), i.e., satisfying

$\displaystyle C^{2}(\mathbb{E}e^{|\epsilon_{j}|^{2}/C^{2}}-1)\leq
C^{\prime},\quad j=1,...,n.$

######  Assumption 4.2.

There exists $m_{0}>D/2$ such that

$\displaystyle
c_{1}(1+\|\bm{\omega}\|_{2}^{2})^{-m_{0}}\leq\mathcal{F}(K)(\bm{\omega})\leq
c_{2}(1+\|\bm{\omega}\|_{2}^{2})^{-m_{0}},\forall\bm{\omega}\in\mathbb{R}^{D}.$
(13)

######  Assumption 4.3 (Tensor kernel function).

The kernel function $K$ can be expressed as $K=\prod_{j=1}^{D}K_{j}$, where
$K_{j}$’s are one-dimensional kernel functions. There exists $m_{0}>1/2$ such
that for $j=1,\ldots,D$,

$\displaystyle
c_{1}(1+\omega_{j}^{2})^{-m_{0}}\leq\mathcal{F}(K_{j})(\omega)\leq
c_{2}(1+\omega_{j}^{2})^{-m_{0}},\forall\omega_{j}\in\mathbb{R}.$ (14)

######  Example 4.1.

A class of kernel functions satisfying Assumption 4.2 is the isotropic Matérn
kernel functions (Williams and Rasmussen,, 2006). With reparameterization, the
Matérn kernel function is given by

$\displaystyle
K(\bm{x})=\frac{(2\phi\sqrt{m_{0}-D/2}\|\bm{x}\|_{2})^{m_{0}-D/2}}{\Gamma(m_{0}-D/2)2^{m_{0}-D/2-1}}B_{m_{0}-D/2}(2\phi\sqrt{m_{0}-D/2}\|\bm{x}\|_{2}),$
(15)

with the Fourier transform (Tuo and Wu,, 2016)

$\displaystyle\mathcal{F}(K)(\bm{\omega})=\pi^{-D/2}\frac{\Gamma(m_{0})}{\Gamma(m_{0}-D/2)}(4\phi^{2}(m_{0}-D/2))^{m_{0}-D/2}(4\phi^{2}(m_{0}-D/2)+\|{\bm{\omega}}\|^{2})^{-m_{0}},$
(16)

where $\phi>0$, and $B_{m_{0}-D/2}$ is the modified Bessel function of the
second kind. It can be seen that (16) is bounded above and below by
$(1+\|\bm{\omega}\|_{2}^{2})^{-m_{0}}$, up to a constant multiplier.

Another example satisfying Assumption 4.2 is the generalized Wendland kernel
function (Wendland,, 2004; Gneiting,, 2002; Chernih and Hubbert,, 2014;
Bevilacqua et al.,, 2019; Fasshauer and McCourt,, 2015), defined as

$\displaystyle K_{GW}(\bm{x})=\left\\{\begin{array}[]{lc}\frac{1}{{\rm
Beta}(2\kappa,\mu+1)}\int_{\|\phi\bm{x}\|_{2}}^{1}u(u^{2}-\|\phi\bm{x}\|_{2}^{2})^{\kappa-1}(1-u)^{\mu}{\rm
d}u,&0\leq\|\bm{x}\|<\frac{1}{\phi},\\\
0,&\|\bm{x}\|_{2}\geq\frac{1}{\phi},\end{array}\right.$ (19)

where $\phi,\kappa>0$ and $\mu\geq(D+1)/2+\kappa$, and ${\rm Beta}$ denotes
the beta function. Theorem 1 of Bevilacqua et al., (2019) shows that (19)
satisfies Assumption 4.2 with $m_{0}=(D+1)/2+\kappa$.

If the kernel function $K=\prod_{j=1}^{D}K_{j}$, and each $K_{j}$ is a one-
dimensional Matérn kernel function or generalized Wendland kernel function,
then Assumption 4.3 is satisfied.

######  Assumption 4.4 (Random smoothing noise).

The elements of $\bm{\varepsilon}_{k}$ are i.i.d. mean zero sub-Gaussian
random variables. Furthermore, we consider three cases of
$\bm{\varepsilon}_{k}$ as follows, where $\sigma_{n}^{2}$’s are positive
parameters to be specified later in Section 4.

  * (C1)

(Polynomial noise) There exists $m_{\varepsilon}>D/2$ such that the
characteristic function of $\bm{\varepsilon}_{k}$ satisfies

$\displaystyle
c_{1}(1+\sigma_{n}^{2}\|\bm{\omega}\|_{2}^{2})^{-m_{\varepsilon}}\leq\mathbb{E}(e^{i\bm{\omega}^{T}\bm{\varepsilon}_{k}})\leq
c_{2}(1+\sigma_{n}^{2}\|\bm{\omega}\|_{2}^{2})^{-m_{\varepsilon}},\forall\bm{\omega}\in\mathbb{R}^{D}.$

  * (C2)

(Tensor Polynomial noise) There exists $m_{\varepsilon}>1/2$ such that the
characteristic function of $\bm{\varepsilon}_{k}$ satisfies

$\displaystyle
c_{1}\prod_{j=1}^{D}(1+\sigma_{n}^{2}\omega_{j}^{2})^{-m_{\varepsilon}}\leq\mathbb{E}(e^{i\bm{\omega}^{T}\bm{\varepsilon}_{k}})\leq
c_{2}\prod_{j=1}^{D}(1+\sigma_{n}^{2}\omega_{j}^{2})^{-m_{\varepsilon}},\forall\bm{\omega}=(\omega_{1},\ldots,\omega_{D})\in\mathbb{R}^{D}.$

  * (C3)

(Gaussian noise) The elements of $\bm{\varepsilon}_{k}$ are normally
distributed with variance $\sigma_{n}^{2}$.

Here the constants $c_{1}$ and $c_{2}$ do not depend on $\sigma_{n}$ and
$m_{\varepsilon}$. We call $\sigma_{n}$ the smoothing scale in this work.

######  Example 4.2.

It is easy to construct distributions satisfying (C1) or (C2). For example,
the generalized Laplace distribution with parameter $s$ has a density function
(Kozubowski et al.,, 2013; Kotz et al.,, 2001)

$\displaystyle
p_{\varepsilon}(\bm{x})=\frac{2^{1-s}}{(2\pi)^{D/2}\Gamma(s)}(\sqrt{2}\|\bm{x}\|_{2})^{s+D/2}B_{s-D/2}\left(\sqrt{2}\|\bm{x}\|_{2}\right),$
(20)

where $\Gamma$ is the Gamma function, and $B_{s-D/2}$ is the modified Bessel
function of the second kind. It can be shown that the generalized Laplace
distribution has the characteristic function

$\displaystyle\mathbb{E}_{\bm{X}}(e^{i\bm{\omega}^{T}\bm{X}})=\left(1+\frac{1}{2}\bm{\omega}^{T}\bm{\omega}\right)^{-s}.$

Then $\bm{\varepsilon}_{k}=\sigma_{n}\bm{X}$ satisfies Assumption 4.4 (C1).

If each component of $\bm{\varepsilon}_{k}/\sigma_{n}$ has a univariate
generalized Laplace distribution and all components are independent, then
Assumption 4.4 (C2) is satisfied.

Assumption 4.1 assumes that the observation error is sub-Gaussian, which is a
standard assumption in nonparametric literature. See van de Geer, (2000) for
example. Assumption 4.2 assumes that the Fourier transform of the kernel
function $K(\cdot-\cdot)$ has an algebraic decay. Under this assumption,
Corollary 10.13 of Wendland, (2004) shows that the reproducing kernel Hilbert
space $\mathcal{H}_{K}(\mathbb{R}^{D})$ coincides with the Sobolev space
$\mathcal{W}^{m_{0}}(\mathbb{R}^{D})$, with equivalent norms. More details on
this can be found in Section 3.1. Assumption 4.3 states that the kernel
function $K$ has a tensor structure, and the Fourier transform of each
component $K_{j}$ has an algebraic decay. Assumptions 4.2 and 4.3 will be used
in Sections 4.2 and 4.3, respectively. Assumption 4.4 imposes conditions on
the noise $\bm{\varepsilon}_{k}$’s and considers three types of augmentations:
polynomial noise, tensor polynomial noise, and Gaussian noise. The
corresponding smoothing techniques are referred to as polynomial smoothing,
tensor polynomial smoothing, and Gaussian smoothing, respectively.

###  4.2 Low Intrinsic Dimension Space

We first consider $\Omega$ with finite intrinsic dimension. The intrinsic
dimension provides a “measure of the complexity” for the region of interest
$\Omega$. The definition of the intrinsic dimension depends on the covering
number; see Definition 2.1 of van de Geer, (2000) for example.

######  Definition 4.1 (Covering number).

Consider a subset $\mathcal{A}\subset\mathcal{G}$ where $\mathcal{G}$ is a
normed space. For a given $\delta>0$, the covering number of $\mathcal{A}$,
denoted by $\mathcal{N}_{\mathcal{G}}(\delta,\mathcal{A})$, is defined by the
smallest integer $M$ such that $\mathcal{A}$ can be covered by $M$ balls with
radius $\delta$ and centers $\bm{x}_{1},...,\bm{x}_{M}\in\mathcal{G}$.

######  Assumption 4.5 (Low intrinsic dimension).

There exist positive constants $c_{1}$ and $d\leq D$ such that for all
$\delta\in(0,1)$, we have

$\displaystyle\mathcal{N}_{\ell_{\infty}^{D}}(\delta,\Omega)\leq
c_{1}\delta^{-d},$

where $\ell_{\infty}^{D}$ is the $\mathbb{R}^{D}$ space equipped with
$\ell_{\infty}$ norm.

For discussion and examples of regions that satisfy Assumption 4.5, we refer
to Hamm and Steinwart, 2021a . In particular, if $\Omega\subset\mathbb{R}^{D}$
is a bounded region with positive Lebesgue measure or a bounded
$D^{\prime}$-dimensional differentiable manifold, then Assumption 4.5 holds
with $d=D$ and $d=D^{\prime}$, respectively.

Besides the low intrinsic dimension, our theoretical results depend on the
smoothness of the underlying function. Because we are considering function
space on a finite intrinsic dimensional space, which may have Lebesgue measure
zero, the usual definition of (fractional) Sobolev space via Fourier transform
stated in Section 3.1 cannot be directly applied in our case. Thus, we need to
introduce our notion of the smoothness of functions on finite intrinsic
dimension space. Specifically, we impose the following assumption on the
underlying true function $f^{*}$.

######  Assumption 4.6.

There exists a region $\Omega_{1}$ with positive Lebesgue measure and a
Lipschitz boundary such that $\Omega\subset\Omega_{1}$. The underlying true
function $f^{*}$ is well-defined on $\Omega_{1}$ and
$m_{f}=\mathop{\mathrm{arginf}}_{m>D/2}\\{m:f^{*}\in\mathcal{W}^{m}(\Omega_{1})\\}$
with $f^{*}\in\mathcal{W}^{m_{f}}(\Omega_{1})$, and $m_{f}>D/2$.

In Assumption 4.6, we assume that the boundary of $\Omega_{1}$ is
“sufficiently regular” (see Leoni, (2017) for the definition of Lipschitz
boundary) and $\Omega$ can be contained by $\Omega_{1}$. Thus, the extension
theorem (DeVore and Sharpley,, 1993) ensures that there exists an extension
operator from $L_{2}(\Omega_{1})$ to $L_{2}(\mathbb{R}^{D})$ and the
smoothness of each function is maintained. With Assumption 4.6, we use $m_{f}$
to denote the smoothness of $f^{*}$. By some well-known extension theorems
(see, for example, DeVore and Sharpley, (1993); Evans, (2009); Stein, (1970)),
if $D=d$, then our notion of smoothness coincides with the smoothness of
functions on the whole space $\mathbb{R}^{D}$.

Now we are ready to present the main theorems in this subsection. Theorems 4.1
and 4.2 state the convergence rates when applying polynomial smoothing and
Gaussian smoothing, respectively.

######  Theorem 4.1 (Polynomial smoothing).

Suppose Assumptions 4.1, 4.2, 4.4 (C1), 4.5 and 4.6 are satisfied. Let
$f_{t}(\bm{x})$ be as in (10) and $\beta=n^{-1}C_{1}$ with
$C_{1}\leq(2\sup_{\bm{x}\in\mathbb{R}^{D}}K_{S}(\bm{x}))^{-1}$. Suppose the
smoothing scale $\sigma_{n}\asymp n^{\nu}$ with $\nu\leq 0$. Suppose one of
the following holds:

  1. 1.

There is no weight decay in the gradient descent, and the iteration number $t$
satisfies $t\asymp
n^{\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}\sigma_{n}^{2m_{\varepsilon}}$

  2. 2.

There is weight decay in the gradient descent with $\alpha\asymp
n^{-1-\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}\sigma_{n}^{-2m_{\varepsilon}}$,
and the iteration number satisfies $t\geq
C_{2}(\frac{m_{f}}{2m_{f}+d}+1/2)\log n/(\log(1-\alpha))$.

Then by setting $m_{\varepsilon}=2d^{-1}(2D\max(m_{0},m_{f})+m_{0}d)\log
n-m_{0}$ and

$\displaystyle\nu=\left\\{\begin{array}[]{ll}-\frac{2(2m_{0}+2m_{\varepsilon})D-(2m_{0}+2m_{\varepsilon}-D)d}{(2m_{f}+d)(4m_{\varepsilon}D-(2m_{0}+2(1-(\log
n)^{-1})m_{\varepsilon}-D)d)}<0,&D>d,\\\ 0,&D=d,\end{array}\right.$

we have

$\displaystyle\|f_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}^{2}=$ $\displaystyle
O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}}(\log n)^{2m_{f}+1}\right).$

for $N>N_{0}$, where $N$ is the number of augmentations, and $N_{0}$ depends
on $n$ (specified in (45)).

######  Theorem 4.2 (Gaussian smoothing).

Suppose Assumptions 4.1, 4.2, 4.4 (C3), 4.5, and 4.6 are satisfied. Let
$f_{t}(\bm{x})$ be as in (10), $\beta=n^{-1}C_{1}$ with
$C_{1}\leq(2\sup_{\bm{x}\in\mathbb{R}^{D}}K_{S}(\bm{x}))^{-1}$, and
$\sigma_{n}\asymp n^{-\frac{1}{2m_{f}+d}}$. Suppose one of the following
holds:

  1. 1.

There is no weight decay in the gradient descent, and the iteration number $t$
satisfies $t\asymp n^{\frac{2m_{0}+2m_{f}}{2m_{f}+d}}$

  2. 2.

There is weight decay in the gradient descent with $\alpha\asymp
n^{-1-\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}$, and the iteration number
satisfies $t\geq C_{2}(\frac{m_{f}}{2m_{f}+d}+1/2)\log n/(\log(1-\alpha))$.

Then we have

$\displaystyle\|f^{*}-\hat{f}_{t}\|_{L_{2}(P_{\mathbf{X}})}^{2}=O_{\mathbb{P}}(n^{-\frac{2m_{f}}{2m_{f}+d}}(\log
n)^{D+1}),$ (21)

when $N>N_{0}$, where $N$ is the number of augmentations, and $N_{0}$ depends
on $n$ (specified in (77)).

######  Remark 4.1.

We require $\beta=n^{-1}C_{1}$ with
$C_{1}\leq(2\sup_{\bm{x}\in\mathbb{R}^{D}}K_{S}(\bm{x}))^{-1}$ in both
Theorems 4.1 and 4.2 is because by Gershgorin’s theorem (Varga,, 2010), we
have for sufficiently large $n$,

$\displaystyle\beta\eta_{1}(\mathbf{K})+\alpha\leq\beta
n\max_{j,k}|K_{S}(\bm{x}_{j},\bm{x}_{k})|+\alpha<1,$

which ensures that the gradient descent algorithm can converge.

If the region $\Omega$ has a positive Lebesgue measure, then it has been shown
that the optimal convergence rate is $n^{-m_{f}/(2m_{f}+D)}$ (Stone,, 1982).
By random smoothing, the gradient descent with early stopping can achieve the
optimal convergence rate in this case, up to a logarithm term. Furthermore, it
can adapt to the low intrinsic dimension case, where $\Omega$ can have
Lebesgue measure zero. In Hamm and Steinwart, 2021a , it is strongly
hypothesized that the convergence rate $n^{-m_{f}/(2m_{f}+d)}$ is optimal.
Although our definition of the smoothness is different, we have the same
hypothesis and leave its exploration as a future work.

It is worth noting that our approach differs from that in Hamm and Steinwart,
2021a , and therefore, we can investigate the effects of polynomial smoothing,
which may have its own interest. Such non-smooth noise can shed light on non-
smooth augmentations commonly used in practice. Furthermore, we obtain an
identical result as in Hamm and Steinwart, 2021a  if we use Gaussian
smoothing. Comparing the convergence rates in Theorems 4.1 and 4.2, we find
that the convergence rate by polynomial smoothing is slightly worse than that
of Gaussian smoothing, since $m_{f}>D/2$ (Assumption 4.6). In comparison,
Eberts and Steinwart, (2013) achieved convergence rate of the similar form
$n^{-2m_{f}/(2m_{f}+d)+\xi}$ by applying kernel ridge regression with Gaussian
kernel functions, where $\xi$ can be any value strictly larger than zero.
Clearly, this rate is slower than those in Hamm and Steinwart, 2021a  and
ours. Under additional assumptions such as a compact Riemannian manifold input
space and the underlying function having Lipschitz continuity $m_{f}\in(0,1]$,
Ye and Zhou, (2008) derived convergence rates of the form
$\left(\log^{2}(n)/n\right)^{m_{f}/(8m_{f}+4d)}$. Instead of kernel ridge
regression, Yang and Dunson, (2016) focused on Bayesian regression with
Gaussian process and proved the convergence rate $n^{-2m_{f}/(2m_{f}+d)}(\log
n)^{d+1}$. However, their theorem is limited by a compact low dimensional
differentiable manifold input space, and the condition $m_{f}\leq 2$. As a
comparison, we do not require such restrictive assumptions.

From a different perspective of early stopping, we consider both cases with
and without weight decay, while existing studies only consider the case
without weight decay. With weight decay, one can achieve the same convergence
rate but with a much smaller iteration number. Specifically, the iteration
number should be polynomial in $n$ without weight decay, which can be reduced
to polynomial in $\log n$ if one applies weight decay. This also justifies the
use of weight decay in practice. Besides, the random smoothing kernel enables
us to establish connections with data augmentation and we further explain the
effectiveness of using augmentation, which may lead to a new interpretation of
using augmentations in deep learning.

Our approach to studying early stopping is distinct from previous studies in
the literature (see, e.g., Dieuleveut and Bach, (2016); Yao et al., (2007);
Pillaud-Vivien et al., (2018); Raskutti et al., (2014)), which typically use
integral operator techniques and impose assumptions on the eigenvalues of the
kernel function (which always exists by Mercer’s theorem). However, such
assumptions cannot be easily applied to the low intrinsic dimension case, as
it is unclear how eigenvalues behave in this regime. Additionally, previous
studies often impose a “source condition” that requires the kernel function to
have finite smoothness, which is not satisfied when using Gaussian smoothing
to construct the random smoothing kernel. Therefore, even for the special case
where the intrinsic dimension is equal to the ambient dimension, Theorems 4.1
and 4.2 improve upon previous results in the early stopping literature.

######  Remark 4.2.

In general, the Bessel potential space used in our work is different from the
Besov space used in Hamm and Steinwart, 2021a . Specifically, the Bessel
potential space is obtained via complex interpolation, while the Besov space
is constructed by real interpolation. For a more thorough explanation, readers
may refer to Edmunds and Triebel, (2008). We chose to use the Bessel potential
space because of its natural connection to the Fourier transform and the
characteristic function of a random variable, which allowed us to study the
impact of the augmentations considered in our work.

######  Remark 4.3.

There are some other notions of smoothness in the literature. For example,
Hamm and Steinwart, 2021a  define the smoothness induced by the Besov spaces,
and Yang and Dunson, (2016) assume $f^{*}$ has $k$-th continuous derivatives.
Another alternative definition of the Sobolev space on $\Omega$ is via
Sobolev–Slobodeckij spaces. For simplicity, let
$\Omega\subset\mathbb{R}^{D-1}$. For a function $f$, $\theta\in(0,1)$, and
$s>0$, define the Slobodeckij seminorm

$\displaystyle|f|_{\theta,\Omega}=\left(\int_{\Omega\times\Omega}\frac{|f(\bm{x})-f(\bm{x}^{\prime})|^{2}}{\|\bm{x}-\bm{x}^{\prime}\|_{2}^{2\theta+D-1}}{\rm
d}\bm{x}{\rm d}\bm{x}^{\prime}\right)^{1/2}.$

Then the Sobolev–Slobodeckij space on $\Omega$, denoted by $W^{s}(\Omega)$, is
defined by

$\displaystyle W^{s}(\Omega)=\left\\{f:f\in W^{\lfloor
s\rfloor}(\Omega):\sup_{\alpha=\lfloor
s\rfloor}|D^{\alpha}f|_{\theta,\Omega}<\infty\right\\},$

with norm

$\displaystyle\|f\|_{W^{s}(\Omega)}=\|f\|_{W^{\lfloor
s\rfloor}(\Omega)}+\sup_{\alpha=\lfloor
s\rfloor}|D^{\alpha}f|_{\theta,\Omega},$

and $D^{\alpha}f:=\frac{\partial^{|\alpha|}}{\partial x_{1}^{\alpha_{1}}\ldots
x_{d}^{\alpha_{d}}}f$ denotes the $\alpha$-th (weak) derivative of a function
$f$ with $|\alpha|=\alpha_{1}+\ldots+\alpha_{d}$ for a multi-index
$\alpha=(\alpha_{1},\ldots,\alpha_{d})\in\mathbb{N}_{0}^{d}$. By the trace
extension theorem (Triebel,, 2010), there exists an extension operator such
that the extended function $f_{E}\in W^{s+1/2}(\mathbb{R}^{D})$ and
$f_{E}|_{\mathbb{R}^{D-1}}=f$, which implies $m_{f}=s+1/2$ if $\Omega$ has a
positive Lebesgue measure in $\mathbb{R}^{D-1}$.

###  4.3 Tensor Reproducing Kernel Hilbert Space

In this section, we consider a low-dimensional structure for the function
class, specifically a tensor reproducing kernel Hilbert space. Let
$K=\prod_{j=1}^{D}K_{j}$ be kernel functions that satisfy Assumption 4.3,
while $\Omega$ can have a low intrinsic dimensional structure, as discussed in
Section 4.2, or have a positive Lebesgue measure in $\mathbb{R}^{D}$.

Our theoretical results in this section are based on mixed smooth Sobolev
spaces, denoted by $\mathcal{MW}^{m}(\mathbb{R}^{D})$, where $m>1/2$. For a
function $f$ defined on $\mathbb{R}^{D}$, the mixed smooth Sobolev norm is
defined as

$\displaystyle\|f\|_{\mathcal{MW}^{m}(\mathbb{R}^{D})}=\left(\int_{\mathbb{R}^{D}}|\mathcal{F}(f)(\bm{\omega})|^{2}\prod_{j=1}^{D}(1+|\omega_{j}|^{2})^{m}{\rm
d}\bm{\omega}\right)^{1/2},$ (22)

and the mixed smooth Sobolev spaces on $\Omega$ can be defined via restriction
similar to the Sobolev spaces. In fact, the mixed smooth Sobolev space is a
tensor product of one-dimensional Sobolev spaces, and it can be shown that
$\mathcal{MW}^{m_{0}}(\mathbb{R}^{D})$ is equivalent to the tensor reproducing
kernel Hilbert space generated by kernel function $K=\prod_{j=1}^{D}K_{j}$
satisfying Assumption 4.3. Because of such a tensor structure, it is often
considered as a reasonable model reducing the complexity in high-dimensional
spaces (Kühn et al.,, 2015; Dũng,, 2021). For instance, the mixed smooth
Sobolev spaces are utilized in high-dimensional approximation and numerical
methods of PDE (Bungartz and Griebel,, 1999), data mining (Garcke et al.,,
2001), and deep neural networks (Dũng,, 2021).

If the underlying function belongs to some mixed smooth Sobolev space, then it
can be shown that by applying appropriate augmentations, we can achieve a fast
convergence rate, which nearly coincides with the minimax rate in the one-
dimensional case, up to a logarithmic term. Similar to Assumption 4.6, we
assume that $f^{*}$ can be extended to some “regular space” with positive
Lebesgue measure, as follows.

######  Assumption 4.7.

There exists a region $\Omega_{1}$ with positive Lebesgue measure and a
Lipschitz boundary such that $\Omega\subset\Omega_{1}$, and the underlying
true function $f^{*}$ is well-defined on $\Omega_{1}$ and
$f^{*}\in\mathcal{MW}^{m_{f}}(\Omega_{1})$.

The following theorem states the convergence rate when applying tensor
polynomial smoothing in the tensor RKHS case.

######  Theorem 4.3 (Tensor polynomial smoothing).

Suppose Assumptions 4.1, 4.3, 4.4 (C2), 4.5, and 4.7 are satisfied. Let
$f_{t}(\bm{x})$ be as in (10) and $\beta=n^{-1}C_{1}$ with
$C_{1}\leq(2\sup_{\bm{x}\in\mathbb{R}^{D}}K_{S}(\bm{x}))^{-1}$. Let
$m_{\varepsilon}+m_{0}\geq m_{f}$, and the smoothing scale $\sigma_{n}\asymp
1$.

Then the following statements are true with $N>N_{0}$, where $N$ is the number
of augmentations, and $N_{0}$ depends on $n$ (specified in (87)). Suppose one
of the following holds:

  1. 1.

There is no weight decay in the gradient descent, and the iteration number $t$
satisfies $t\asymp n^{\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+1}}(\log
n)^{\frac{2(D-1)(m_{0}+m_{\varepsilon})+1}{2m_{f}+1}}$

  2. 2.

There is weight decay in the gradient descent with $\alpha\asymp
n^{-1-\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}(\log
n)^{\frac{2(D-1)(m_{0}+m_{\varepsilon})+1}{2m_{f}+1}}$, and the iteration
number satisfies $t\geq C_{2}(\frac{m_{f}}{2m_{f}+1}+1/2)\log
n/(\log(1-\alpha))$.

Then we have

$\displaystyle\|f_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}^{2}=$ $\displaystyle
O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+1}}(\log
n)^{\frac{2m_{f}}{2m_{f}+1}\left(D-1+\frac{1}{2(m_{0}+m_{\varepsilon})}\right)}\right).$
(23)

Based on Theorem 4.3, tensor polynomial smoothing leads to a convergence rate
of tensor RKHS, which is $O_{\mathbb{P}}(n^{-\frac{2m_{f}}{2m_{f}+1}}(\log
n)^{\frac{2m_{f}}{2m_{f}+1}(D-1+\frac{1}{2(m_{0}+m_{\varepsilon})})})$. This
convergence rate is almost the same as the optimal convergence rate in the
one-dimensional case $O_{\mathbb{P}}(n^{-\frac{2m_{f}}{2m_{f}+1}})$, differing
only by a logarithmic term.

Moreover, compared to Theorem 4.1, Theorem 4.3 has less stringent requirements
for tensor polynomial smoothing when Assumption 4.7 holds. Specifically,
Theorem 4.3 allows for $m_{\varepsilon}$ to be a constant as long as
$m_{\varepsilon}+m_{0}\geq m_{f}$, whereas Theorem 4.1 requires
$m_{\varepsilon}$ to be comparable to $\log n$. Additionally, while the
smoothing scale $\sigma_{n}$ in Theorem 4.1 demands careful selection, Theorem
4.3 permits a constant smoothing scale $\sigma_{n}$. These differences suggest
that the tensor RKHS has a simpler structure than the RKHS even in a low
intrinsic dimension space. The convergence rate in Theorem 4.3 does not depend
on the low intrinsic dimension of $\Omega$, and is almost dimension-free.
Moreover, because the power of the logarithmic term in (23) decreases as
$m_{\varepsilon}$ increases, the convergence rate in Theorem 4.3 decreases as
$m_{\varepsilon}$ increases, encouraging the use of a smoother tensor
polynomial smoothing for faster convergence. This aligns with the results in
Theorem 4.1 and Theorem 4.2, as Gaussian smoothing may yield faster
convergence rates than polynomial smoothing. Few studies have explored tensor
RKHSs with early stopping, and our findings can provide valuable insights into
this area.

######  Remark 4.4.

For any $\mathcal{W}^{m_{f}}(\mathbb{R}^{D})$ with $m_{f}>D/2$, there exist
$m^{*}>1/2$ such that
$\mathcal{W}^{m_{f}}(\mathbb{R}^{D})\hookrightarrow\mathcal{MW}^{m^{*}}(\mathbb{R}^{D})$
and
$\mathcal{MW}^{m^{*}}(\mathbb{R}^{D})\hookrightarrow\mathcal{C}(\mathbb{R}^{D})$.
Thus, the capacity of $\mathcal{MW}^{m^{*}}(\mathbb{R}^{D})$ is high-enough
for any approximation problem which can be solved by assuming that the
underlying true function lies in some Sobolev space.

##  5 Numerical Studies

In this section, we enhance our theoretical findings by experimentally
validating the effectiveness of the random smoothing kernel with data
augmentation and early stopping on synthetic datasets. We focus on three data
spaces with dimensions $D=1,2,3$, as illustrated in Figure 1, where
$\bm{x}_{j}$ samples are uniformly drawn.

In our experiments, the underlying function $f^{*}$ is obtained by drawing
random sample paths from the Gaussian process with the Matérn covariance
function. This covariance function is widely used in Gaussian process
modeling. We adopt the Matérn covariance function with the following form:

$\displaystyle
K_{\nu}(\bm{x})=\sigma^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{\|\bm{x}\|_{2}}{\rho}\right)^{\nu}B_{\nu}\left(\sqrt{2\nu}\frac{\|\bm{x}\|_{2}}{\rho}\right),$
(24)

where $\sigma,\phi,\nu>0$, $\Gamma$ is the Gamma function, and $B_{\nu}$ is
the modified Bessel function of the second kind. In order to make $f^{*}$
smoother, we set the smoothness parameter $\nu=5.0$ for Matérn kernel (24).
The error $\epsilon_{j}$’s are i.i.d. Gaussian with mean zero and variance
0.01.

We utilize two-hidden-layer neural networks with ReLU activation (Nair and
Hinton,, 2010) as our predictor. Each hidden layer of the neural network
comprises 100 nodes, and all weights are initialized using Kaiming
Initialization (He et al.,, 2015). For random smoothing, we experiment with
both non-smooth Laplace noise and smooth Gaussian noise. To be precise, each
element of $\bm{\varepsilon}_{k}$ is randomly sampled from either
$\mathcal{N}(0,\sigma^{2})$ or $Laplace(0,b)$. For more experiment details and
additional results, we refer to Appendix N.

Figure 1: Simulated data spaces in the forms of: line ($D=1$), circle ($D=2$)
and sphere ($D=3$).

Figure 2 presents a visualization of the underlying truth (blue curve),
training data (blue dots), and neural network predictions (orange dots) when
the training size is 50. The underlying truth is smooth since we use a smooth
kernel. However, the neural network predictions without random smoothing are
not smooth due to the low smoothness of the ReLU activation function and tend
to overfit the noise. Upon applying random smoothing, the neural network
predictions become smoother and approach the underlying truth.

Figure 2: Visualization of the underlying truth (blue curve), training data
(blue dots), and neural network predictions (orange dots) when training size
is 50, where the first and second rows represent cases with weight decay and
early stopping, respectively. It is obvious to see that the optimization
without random smoothing will be more vulnerable to noise.

Figure 3 and Figure 4 further show the underlying truth (blue curve), training
data (blue dots), and neural network predictions (orange dots) when the
training size is 100 and 200, respectively. Although increasing the training
size improves smoothness in cases like size 200 with weight decay, the fitted
curve still experiences a perturbation from overfitted noise compared to
examples where random smoothing is applied.

Figure 3: Underlying truth (blue curve), training data (blue dots), and neural
network predictions (orange dots) when training size is 100. Figure 4:
Underlying truth (blue curve), training data (blue dots), and neural network
predictions (orange dots) when training size is 200.

Table 1 presents a summary of the test $l_{2}$ loss under different settings.
Both Gaussian smoothing and polynomial smoothing (random smoothing with
Laplacian noise) improve the $l_{2}$ loss in all settings, demonstrating the
effectiveness of random smoothing. Figure 5 further investigates how the
$l_{2}$ loss changes concerning the smoothing scale $\sigma_{n}$ when $D=1$.
The plot shows a U-shaped curve, indicating that an optimal smoothing can
minimize the $l_{2}$ loss, while either smaller or larger values will result
in a larger $l_{2}$ loss. It is worth noting that when the training size is
small, such as size 50, the U-shape curve in Figure 5 may be less distinct due
to noise introduced by early stopping based on a small validation set. Another
observation from Figure 5 is that the optimal smoothing scales exhibit a
decreasing trend as the sample size increases, as indicated by Theorem 4.1 and
Theorem 4.2. Additionally, Figure 6 and Figure 7 depict the U-shaped curves of
$l_{2}$ loss changes concerning smoothing scale when $D=2$ and $D=3$,
respectively. While it is possible that some red points may not be accurately
placed due to a small validation set, the optimal smoothing scales exhibit a
decreasing trend with respect to training size, which is consistent with the
trend observed in $D=1$ as depicted in Figure 5.

Table 1: Test $l_{2}$ loss of SGD with early stopping. “G”, “L”, and “N” correspond to random smoothing with Gaussian noise, random smoothing with Laplacian noise, and no random smoothing. The smallest losses are underlined. Dim | Type | With weight decay | Early stopping
---|---|---|---
Training size | Training size
50 | 100 | 200 | 50 | 100 | 200
D=1 | G | 1.7466e-03 | 9.8343e-04 | 9.1924e-04 | 1.3468e-03 | 7.5579e-04 | 5.8775e-04
L | 1.6765e-03 | 9.3367e-04 | 8.2806e-04 | 2.0638e-03 | 9.2128e-04 | 6.5118e-04
N | 1.9381e-03 | 1.3045e-03 | 1.1135e-03 | 2.2168e-03 | 1.2985e-03 | 8.4292e-04
D=2 | G | 6.4208e-03 | 3.1423e-03 | 2.1842e-03 | 6.7205e-03 | 3.5027e-03 | 1.7132e-03
L | 6.4676e-03 | 2.9491e-03 | 2.2136e-03 | 8.2725e-03 | 3.9418e-03 | 1.7674e-03
N | 9.2474e-03 | 4.5782e-03 | 2.5810e-03 | 1.2628e-02 | 6.2301e-03 | 3.1396e-03
D=3 | G | 1.6498e-02 | 7.2578e-03 | 3.9938e-03 | 1.4852e-02 | 7.1306e-03 | 3.7147e-03
L | 1.6599e-02 | 6.9336e-03 | 4.4334e-03 | 1.5167e-02 | 6.6471e-03 | 3.8615e-03
N | 2.0987e-02 | 8.1158e-03 | 4.5752e-03 | 2.0178e-02 | 8.4932e-03 | 4.9460e-03
With Weight Decay | Early Stopping
---|---
|
Figure 5: Loss changes according to smoothing scale with training size increase from 50 to 200 in 1d data space. The red points represent the optimal smoothing scales selected based on the validation set. With Weight Decay | Early Stopping
---|---
|
Figure 6: Loss changes according to smoothing scale with training size increase from 50 to 200 in 2d data space. The red points represent the optimal smoothing scales selected based on the validation set. With Weight Decay | Early Stopping
---|---
|
Figure 7: Loss changes according to smoothing scale with training size
increase from 50 to 200 in 3d data space. The red points represent the optimal
smoothing scales selected based on the validation set.

##  6 Conclusions and Discussion

This work studies random smoothing kernel and random smoothing regularization,
which have a natural relationship with data augmentations. We consider two
cases: when the region $\Omega$ has a low intrinsic dimension, or when the
kernel function can be presented as a product of one-dimensional kernel
functions. In both cases, we show that by applying random smoothing, with
appropriate early stopping and/or weight decay techniques, the resulting
estimator can achieve fast convergence rates, regardless of the kernel
function used in the construction of the random smoothing kernel estimator.

There are several directions that could be pursued in future research. First,
while we consider noise injection to construct augmentations and use non-
smooth noise to interpret practical non-smooth augmentation techniques, such
as random crop, random mask, and random flip, this interpretation may not be
perfect. For example, the behavior of adding noise may differ from that of
random crop. Furthermore, these practical techniques may also introduce some
prior knowledge on the geometry of the low intrinsic dimension. A sharper
characterization of practical augmentation techniques is needed and will be
pursued in future work.

Second, while we consider gradient descent, we believe that our results can be
generalized to the stochastic gradient descent method. However, the discussion
of the latter is beyond the scope of the current work.

Third, we mainly consider regression in this work, where the square loss is a
natural choice. An interesting extension is to study whether the results
remain true when considering classification, which requires the study of other
loss functions, such as cross-entropy loss and hinge loss.

## References

  * Adams and Fournier, (2003) Adams, R. A. and Fournier, J. J. (2003).  Sobolev Spaces, volume 140.  Academic press.
  * Almeida and Samko, (2006) Almeida, A. and Samko, S. (2006).  Characterization of Riesz and Bessel potentials on variable lebesgue spaces.  Journal of Function Spaces and Applications, 4(2):113–144.
  * (3) Arora, S., Du, S. S., Hu, W., Li, Z., and Wang, R. (2019a).  Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks.  arXiv preprint arXiv:1901.08584.
  * (4) Arora, S., Du, S. S., Li, Z., Salakhutdinov, R., Wang, R., and Yu, D. (2019b).  Harnessing the power of infinitely wide deep nets on small-data tasks.  arXiv preprint arXiv:1910.01663.
  * Bach and Jordan, (2002) Bach, F. R. and Jordan, M. I. (2002).  Kernel independent component analysis.  Journal of machine learning research, 3(Jul):1–48.
  * Bai et al., (2021) Bai, Y., Yang, E., Han, B., Yang, Y., Li, J., Mao, Y., Niu, G., and Liu, T. (2021).  Understanding and improving early stopping for learning with noisy labels.  Advances in Neural Information Processing Systems, 34:24392–24403.
  * Bartlett and Traskin, (2007) Bartlett, P. L. and Traskin, M. (2007).  Adaboost is consistent.  The Journal of Machine Learning Research, 8:2347–2368.
  * Bevilacqua et al., (2019) Bevilacqua, M., Faouzi, T., Furrer, R., Porcu, E., et al. (2019).  Estimation and prediction using generalized Wendland covariance functions under fixed domain asymptotics.  The Annals of Statistics, 47(2):828–856.
  * Blanchard and Krämer, (2016) Blanchard, G. and Krämer, N. (2016).  Convergence rates of kernel conjugate gradient for random design regression.  Analysis and Applications, 14(06):763–794.
  * Blanchard and Mücke, (2018) Blanchard, G. and Mücke, N. (2018).  Optimal rates for regularization of statistical inverse learning problems.  Foundations of Computational Mathematics, 18(4):971–1013.
  * Blum et al., (2020) Blum, A., Dick, T., Manoj, N., and Zhang, H. (2020).  Random smoothing might be unable to certify $\ell_{\infty}$ robustness for high-dimensional images.  Journal of Machine Learning Research, 21(211):1–21.
  * Bühlmann and Yu, (2002) Bühlmann, P. and Yu, B. (2002).  Boosting with the $l_{2}$-loss: Regression and classification.  Journal of American Statistical Association, 98:324–340.
  * Bungartz and Griebel, (1999) Bungartz, H.-J. and Griebel, M. (1999).  A note on the complexity of solving Poisson’s equation for spaces of bounded mixed derivatives.  Journal of Complexity, 15(2):167–199.
  * Cao and Gu, (2020) Cao, Y. and Gu, Q. (2020).  Generalization error bounds of gradient descent for learning over-parameterized deep relu networks.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 3349–3356.
  * Caponnetto and Yao, (2006) Caponnetto, A. and Yao, Y. (2006).  Adaptation for regularization operators in learning theory.
  * Carlini and Wagner, (2017) Carlini, N. and Wagner, D. (2017).  Towards evaluating the robustness of neural networks.  In 2017 ieee symposium on security and privacy (sp), pages 39–57. Ieee.
  * Cervone and Pillai, (2015) Cervone, D. and Pillai, N. S. (2015).  Gaussian process regression with location errors.  arXiv preprint arXiv:1506.08256.
  * Chen et al., (2022) Chen, H., Ding, L., and Tuo, R. (2022).  Kernel packet: An exact and scalable algorithm for gaussian process regression with matérn correlations.  Journal of Machine Learning Research, 23(127):1–32.
  * Chen and Xu, (2020) Chen, L. and Xu, S. (2020).  Deep neural tangent kernel and laplace kernel have the same rkhs.  arXiv preprint arXiv:2009.10683.
  * Chen et al., (2020) Chen, T., Kornblith, S., Norouzi, M., and Hinton, G. (2020).  A simple framework for contrastive learning of visual representations.  arXiv preprint arXiv:2002.05709.
  * Chen and He, (2021) Chen, X. and He, K. (2021).  Exploring simple siamese representation learning.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 15750–15758.
  * Chernih and Hubbert, (2014) Chernih, A. and Hubbert, S. (2014).  Closed form representations and properties of the generalised Wendland functions.  Journal of Approximation Theory, 177:17–33.
  * Cohen et al., (2019) Cohen, J., Rosenfeld, E., and Kolter, Z. (2019).  Certified adversarial robustness via randomized smoothing.  In International Conference on Machine Learning, pages 1310–1320. PMLR.
  * Cressie and Kornak, (2003) Cressie, N. and Kornak, J. (2003).  Spatial statistics in the presence of location error with an application to remote sensing of the environment.  Statistical Science, 18(4):436–456.
  * Csörgő, (1985) Csörgő, S. (1985).  Rates of uniform convergence for the empirical characteristic function.  Acta Sci. Math.(Szeged), 48(1–4):97–102.
  * DeVore and Sharpley, (1993) DeVore, R. A. and Sharpley, R. C. (1993).  Besov spaces on domains in ${R}^{d}$.  Transactions of the American Mathematical Society, 335(2):843–864.
  * Dicker et al., (2017) Dicker, L. H., Foster, D. P., Hsu, D., et al. (2017).  Kernel ridge vs. principal component regression: Minimax bounds and the qualification of regularization operators.  Electronic Journal of Statistics, 11(1):1022–1047.
  * Dieuleveut and Bach, (2016) Dieuleveut, A. and Bach, F. (2016).  Nonparametric stochastic approximation with large step-sizes.  The Annals of Statistics, 44(4):1363–1399.
  * Ding and Zhang, (2022) Ding, L. and Zhang, X. (2022).  Sample and computationally efficient stochastic kriging in high dimensions.  Operations Research.
  * Du et al., (2018) Du, S. S., Zhai, X., Poczos, B., and Singh, A. (2018).  Gradient descent provably optimizes over-parameterized neural networks.  arXiv preprint arXiv:1810.02054.
  * Dũng, (2021) Dũng, D. (2021).  Deep ReLU neural networks in high-dimensional approximation.  Neural Networks, 142:619–635.
  * Dung et al., (2018) Dung, D., Temlyakov, V., and Ullrich, T. (2018).  Hyperbolic cross approximation.  Springer.
  * Eberts and Steinwart, (2013) Eberts, M. and Steinwart, I. (2013).  Optimal regression rates for svms using gaussian kernels.  Electronic Journal of Statistics, 7:1–42.
  * Edmunds and Triebel, (2008) Edmunds, D. E. and Triebel, H. (2008).  Function Spaces, Entropy Numbers, Differential Operators, volume 120.  Cambridge University Press.
  * Engl et al., (1996) Engl, H. W., Hanke, M., and Neubauer, A. (1996).  Regularization of Inverse Problems.  Kluwer Academic Publishers.
  * Evans, (2009) Evans, L. C. (2009).  Partial differential equations (graduate studies in mathematics, vol. 19).  Instructor, 67.
  * Fasshauer and McCourt, (2015) Fasshauer, G. E. and McCourt, M. J. (2015).  Kernel-based approximation methods using MATLAB, volume 19.  World Scientific Publishing Company.
  * Gao et al., (2020) Gao, Z., Hu, R., and Gong, Y. (2020).  Certified robustness of graph classification against topology attack with randomized smoothing.  In GLOBECOM 2020-2020 IEEE Global Communications Conference, pages 1–6. IEEE.
  * Garcke et al., (2001) Garcke, J., Griebel, M., and Thess, M. (2001).  Data mining with sparse grids.  Computing, 67(3):225–253.
  * Geifman et al., (2020) Geifman, A., Yadav, A., Kasten, Y., Galun, M., Jacobs, D., and Ronen, B. (2020).  On the similarity between the laplace and neural tangent kernels.  Advances in Neural Information Processing Systems, 33:1451–1461.
  * Gneiting, (2002) Gneiting, T. (2002).  Stationary covariance functions for space-time data.  Journal of the American Statistical Association, 97:590–600.
  * Goodfellow et al., (2016) Goodfellow, I., Bengio, Y., Courville, A., and Bengio, Y. (2016).  Deep learning, volume 1.  MIT Press.
  * Goodfellow et al., (2014) Goodfellow, I. J., Shlens, J., and Szegedy, C. (2014).  Explaining and harnessing adversarial examples.  arXiv preprint arXiv:1412.6572.
  * Grandvalet et al., (1997) Grandvalet, Y., Canu, S., and Boucheron, S. (1997).  Noise injection: Theoretical prospects.  Neural Computation, 9(5):1093–1108.
  * Gretton, (2015) Gretton, A. (2015).  A simpler condition for consistency of a kernel independence test.  Technical report, University College London.
  * Gretton et al., (2005) Gretton, A., Bousquet, O., Smola, A., and Schölkopf, B. (2005).  Measuring statistical dependence with hilbert-schmidt norms.  In International conference on algorithmic learning theory, pages 63–77. Springer.
  * Gretton et al., (2007) Gretton, A., Fukumizu, K., Teo, C., Song, L., Schölkopf, B., and Smola, A. (2007).  A kernel statistical test of independence.  Advances in neural information processing systems, 20.
  * Grill et al., (2020) Grill, J.-B., Strub, F., Altché, F., Tallec, C., Richemond, P. H., Buchatskaya, E., Doersch, C., Pires, B. A., Guo, Z. D., Azar, M. G., et al. (2020).  Bootstrap your own latent: A new approach to self-supervised learning.  arXiv preprint arXiv:2006.07733.
  * Guo et al., (2017) Guo, Z.-C., Lin, S.-B., and Zhou, D.-X. (2017).  Learning theory of distributed spectral algorithms.  Inverse Problems, 33(7):074009.
  * Gurka et al., (2007) Gurka, P., Harjulehto, P., and Nekvinda, A. (2007).  Bessel potential spaces with variable exponent.  Mathematical Inequalities and Applications, 10(3):661.
  * (51) Hamm, T. and Steinwart, I. (2021a).  Adaptive learning rates for support vector machines working on data with low intrinsic dimension.  The Annals of Statistics, 49(6):3153–3180.
  * (52) Hamm, T. and Steinwart, I. (2021b).  Intrinsic dimension adaptive partitioning for kernel methods.  arXiv preprint arXiv:2107.07750.
  * Hardt et al., (2016) Hardt, M., Recht, B., and Singer, Y. (2016).  Train faster, generalize better: Stability of stochastic gradient descent.  In International conference on machine learning, pages 1225–1234. PMLR.
  * Hastie et al., (2001) Hastie, T., Tibshirani, R., and Friedman, J. (2001).  The Elements of Statistical Learning.  Springer.
  * He et al., (2021) He, K., Chen, X., Xie, S., Li, Y., Dollár, P., and Girshick, R. (2021).  Masked autoencoders are scalable vision learners.  arXiv preprint arXiv:2111.06377.
  * He et al., (2020) He, K., Fan, H., Wu, Y., Xie, S., and Girshick, R. (2020).  Momentum contrast for unsupervised visual representation learning.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 9729–9738.
  * He et al., (2015) He, K., Zhang, X., Ren, S., and Sun, J. (2015).  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In Proceedings of the IEEE international conference on computer vision, pages 1026–1034.
  * Hu et al., (2022) Hu, T., Wang, J., Wang, W., and Li, Z. (2022).  Understanding square loss in training overparametrized neural network classifiers.  Advances in Neural Information Processing Systems, 35:16495–16508.
  * Hu et al., (2021) Hu, T., Wang, W., Lin, C., and Cheng, G. (2021).  Regularization matters: A nonparametric perspective on overparametrized neural network.  In International Conference on Artificial Intelligence and Statistics, pages 829–837. PMLR.
  * Huang et al., (2020) Huang, K., Wang, Y., Tao, M., and Zhao, T. (2020).  Why do deep residual networks generalize better than deep feedforward networks?—a neural tangent kernel perspective.  Advances in neural information processing systems, 33:2698–2709.
  * Jacot et al., (2018) Jacot, A., Gabriel, F., and Hongler, C. (2018).  Neural tangent kernel: Convergence and generalization in neural networks.  In Advances in Neural Information Processing Systems, pages 8571–8580.
  * Kanoh and Sugiyama, (2021) Kanoh, R. and Sugiyama, M. (2021).  A neural tangent kernel perspective of infinite tree ensembles.  arXiv preprint arXiv:2109.04983.
  * Kotz et al., (2001) Kotz, S., Kozubowski, T., and Podgórski, K. (2001).  The Laplace distribution and generalizations: a revisit with applications to communications, economics, engineering, and finance.  Number 183. Springer Science & Business Media.
  * Kozubowski et al., (2013) Kozubowski, T. J., Podgórski, K., and Rychlik, I. (2013).  Multivariate generalized laplace distribution and related random fields.  Journal of Multivariate Analysis, 113:59–72.
  * Krogh and Hertz, (1992) Krogh, A. and Hertz, J. A. (1992).  A simple weight decay can improve generalization.  In Advances in Neural Information Processing Systems, pages 950–957.
  * Kühn, (2011) Kühn, T. (2011).  Covering numbers of gaussian reproducing kernel hilbert spaces.  Journal of Complexity, 27(5):489–499.
  * Kühn et al., (2015) Kühn, T., Sickel, W., and Ullrich, T. (2015).  Approximation of mixed order Sobolev functions on the d-torus: Asymptotics, preasymptotics, and d-dependence.  Constructive Approximation, 42(3):353–398.
  * L Coifman and Weiss, (1977) L Coifman, R. R. and Weiss, G. L. (1977).  Transference methods in analysis, volume 31.  American Mathematical Soc.
  * Leoni, (2017) Leoni, G. (2017).  A first course in Sobolev spaces.  American Mathematical Soc.
  * Li et al., (2020) Li, M., Soltanolkotabi, M., and Oymak, S. (2020).  Gradient descent with early stopping is provably robust to label noise for overparameterized neural networks.  In International conference on artificial intelligence and statistics, pages 4313–4324. PMLR.
  * Li and Liang, (2018) Li, Y. and Liang, Y. (2018).  Learning overparameterized neural networks via stochastic gradient descent on structured data.  In Advances in Neural Information Processing Systems, pages 8157–8166.
  * Li et al., (2019) Li, Z., Wang, R., Yu, D., Du, S. S., Hu, W., Salakhutdinov, R., and Arora, S. (2019).  Enhanced convolutional neural tangent kernels.  arXiv preprint arXiv:1911.00809.
  * Lin and Rosasco, (2017) Lin, J. and Rosasco, L. (2017).  Optimal rates for multi-pass stochastic gradient methods.  The Journal of Machine Learning Research, 18(1):3375–3421.
  * Lin et al., (2016) Lin, J., Rosasco, L., and Zhou, D.-X. (2016).  Iterative regularization for learning with convex loss functions.  The Journal of Machine Learning Research, 17(1):2718–2755.
  * Lin et al., (2017) Lin, S.-B., Guo, X., and Zhou, D.-X. (2017).  Distributed learning with regularized least squares.  The Journal of Machine Learning Research, 18(1):3202–3232.
  * Massart, (2007) Massart, P. (2007).  Concentration Inequalities and Model Selection, volume 6.  Springer.
  * Mehra et al., (2021) Mehra, A., Kailkhura, B., Chen, P.-Y., and Hamm, J. (2021).  How robust are randomized smoothing based defenses to data poisoning?  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 13244–13253.
  * Nair and Hinton, (2010) Nair, V. and Hinton, G. E. (2010).  Rectified linear units improve restricted boltzmann machines.  In Icml.
  * Paulsen and Raghupathi, (2016) Paulsen, V. I. and Raghupathi, M. (2016).  An introduction to the theory of reproducing kernel Hilbert spaces, volume 152.  Cambridge university press.
  * Pillaud-Vivien et al., (2018) Pillaud-Vivien, L., Rudi, A., and Bach, F. (2018).  Statistical optimality of stochastic gradient descent on hard learning problems through multiple passes.  Advances in Neural Information Processing Systems, 31.
  * Prechelt, (1998) Prechelt, L. (1998).  Early stopping-but when?  In Neural Networks: Tricks of the trade, pages 55–69. Springer.
  * Raskutti et al., (2014) Raskutti, G., Wainwright, M. J., and Yu, B. (2014).  Early stopping and non-parametric regression: an optimal data-dependent stopping rule.  The Journal of Machine Learning Research, 15(1):335–366.
  * Rosenfeld et al., (2020) Rosenfeld, E., Winston, E., Ravikumar, P., and Kolter, Z. (2020).  Certified robustness to label-flipping attacks via randomized smoothing.  In International Conference on Machine Learning, pages 8230–8241. PMLR.
  * Saatçi, (2012) Saatçi, Y. (2012).  Scalable inference for structured Gaussian process models.  PhD thesis, University of Cambridge.
  * Salman et al., (2019) Salman, H., Li, J., Razenshteyn, I., Zhang, P., Zhang, H., Bubeck, S., and Yang, G. (2019).  Provably robust deep learning via adversarially trained smoothed classifiers.  Advances in Neural Information Processing Systems, 32.
  * Scott and Nowak, (2006) Scott, C. and Nowak, R. D. (2006).  Minimax-optimal classification with dyadic decision trees.  IEEE transactions on information theory, 52(4):1335–1353.
  * Shorten and Khoshgoftaar, (2019) Shorten, C. and Khoshgoftaar, T. M. (2019).  A survey on image data augmentation for deep learning.  Journal of big data, 6(1):1–48.
  * Stein, (1970) Stein, E. M. (1970).  Singular integrals and differentiability properties of functions, volume 2.  Princeton university press.
  * Steinwart et al., (2009) Steinwart, I., Hush, D. R., and Scovel, C. (2009).  Optimal rates for regularized least squares regression.  In COLT, pages 79–93.
  * Stone, (1982) Stone, C. J. (1982).  Optimal global rates of convergence for nonparametric regression.  The Annals of Statistics, pages 1040–1053.
  * Szabó and Sriperumbudur, (2017) Szabó, Z. and Sriperumbudur, B. K. (2017).  Characteristic and universal tensor product kernels.  J. Mach. Learn. Res., 18:233–1.
  * Triebel, (2010) Triebel, H. (2010).  Theory of Function Spaces II.  Springer Science & Business Media.
  * Tuo et al., (2020) Tuo, R., Wang, Y., and Wu, C. (2020).  On the improved rates of convergence for Matérn-type kernel ridge regression, with application to calibration of computer models.  arXiv preprint arXiv:2001.00152.
  * Tuo and Wu, (2016) Tuo, R. and Wu, C. F. J. (2016).  A theoretical framework for calibration in computer models: Parametrization, estimation and convergence properties.  SIAM/ASA Journal on Uncertainty Quantification, 4(1):767–795.
  * van de Geer, (2000) van de Geer, S. (2000).  Empirical Processes in M-estimation.  Cambridge University Press.
  * Varga, (2010) Varga, R. S. (2010).  Gershgorin and his circles, volume 36.  Springer Science & Business Media.
  * Wahba, (1990) Wahba, G. (1990).  Spline Models for Observational Data, volume 59.  SIAM.
  * Wang et al., (2020) Wang, B., Cao, X., Gong, N. Z., et al. (2020).  On certifying robustness against backdoor attacks via randomized smoothing.  arXiv preprint arXiv:2002.11750.
  * Wang, (2021) Wang, W. (2021).  On the inference of applying Gaussian process modeling to a deterministic function.  Electronic Journal of Statistics, 15(2):5014–5066.
  * Wang et al., (2022) Wang, W., Yue, X., Haaland, B., and Jeff Wu, C. (2022).  Gaussian processes with input location error and applications to the composite parts assembly process.  SIAM/ASA Journal on Uncertainty Quantification, 10(2):619–650.
  * Wei et al., (2017) Wei, Y., Yang, F., and Wainwright, M. J. (2017).  Early stopping for kernel boosting algorithms: A general analysis with localized complexities.  Advances in Neural Information Processing Systems, 30.
  * Wendland, (2004) Wendland, H. (2004).  Scattered Data Approximation, volume 17.  Cambridge University Press.
  * Williams and Rasmussen, (2006) Williams, C. K. and Rasmussen, C. E. (2006).  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA.
  * Wilson and Nickisch, (2015) Wilson, A. and Nickisch, H. (2015).  Kernel interpolation for scalable structured gaussian processes (kiss-gp).  In International conference on machine learning, pages 1775–1784. PMLR.
  * Wu et al., (2006) Wu, Q., Ying, Y., and Zhou, D.-X. (2006).  Learning rates of least-square regularized regression.  Foundations of Computational Mathematics, 6(2):171–192.
  * Yang and Dunson, (2016) Yang, Y. and Dunson, D. B. (2016).  Bayesian manifold regression.  The Annals of Statistics, 44(2):876–905.
  * Yao et al., (2007) Yao, Y., Rosasco, L., and Caponnetto, A. (2007).  On early stopping in gradient descent learning.  Constructive Approximation, 26(2):289–315.
  * Ye and Zhou, (2008) Ye, G.-B. and Zhou, D.-X. (2008).  Learning and approximation by gaussians on riemannian manifolds.  Advances in Computational Mathematics, 29(3):291–310.
  * Ye and Zhou, (2009) Ye, G.-B. and Zhou, D.-X. (2009).  Svm learning and lp approximation by gaussians on riemannian manifolds.  Analysis and Applications, 7(03):309–339.
  * Zhang et al., (2021) Zhang, C., Bengio, S., Hardt, M., Recht, B., and Vinyals, O. (2021).  Understanding deep learning (still) requires rethinking generalization.  Communications of the ACM, 64(3):107–115.
  * Zhang and Yu, (2005) Zhang, T. and Yu, B. (2005).  Boosting with early stopping: Convergence and consistency.  Annals of Statistics, 33(4):1538–1579.
  * Zhou et al., (2019) Zhou, S., Pati, D., Wang, T., Yang, Y., and Carroll, R. J. (2019).  Gaussian processes with errors in variables: Theory and computation.  arXiv preprint arXiv:1910.06235.

##  Appendix A Analysis of Gradient Update and Error Decomposition

Let $\mathbf{X}=(\bm{x}_{1},...,\bm{x}_{n})$, $\alpha>0$ if there is weight
decay, and $\alpha=0$ if there is no weight decay. By the gradient update
rule, we have

$\displaystyle f_{t+1}(\mathbf{X})=$
$\displaystyle\mathbf{K}\bm{w}_{t+1}=\sqrt{\mathbf{K}}\bm{\theta}_{t+1}$
$\displaystyle=$
$\displaystyle\sqrt{\mathbf{K}}\bm{\theta}_{t}-\beta\sqrt{\mathbf{K}}\left(\mathbf{K}\bm{\theta}_{t}-\sqrt{\mathbf{K}}\bm{y}\right)-\alpha\sqrt{\mathbf{K}}\bm{\theta}_{t}$
$\displaystyle=$
$\displaystyle((1-\alpha)\mathbf{I}-\beta\mathbf{K})f_{t}(\mathbf{X})+\beta\mathbf{K}\bm{y},$

which implies

$\displaystyle
f_{t+1}(\mathbf{X})-\beta(\alpha\mathbf{I}+\beta\mathbf{K})^{-1}\mathbf{K}\bm{y}=((1-\alpha)\mathbf{I}-\beta\mathbf{K})(f_{t}(\mathbf{X})-\beta(\alpha\mathbf{I}+\beta\mathbf{K})^{-1}\mathbf{K}\bm{y})$
$\displaystyle=$
$\displaystyle\ldots=-((1-\alpha)\mathbf{I}-\beta\mathbf{K})^{t+1}\beta(\alpha\mathbf{I}+\beta\mathbf{K})^{-1}\mathbf{K}\bm{y},$
(25)

where we recall $f_{0}(\mathbf{X})=\mathbf{0}$. If there is weight decay
(i.e., $\alpha>0$), then it can be seen that

$\displaystyle
f_{t+1}(\mathbf{X})-\mathbf{K}(\alpha/\beta\mathbf{I}+\mathbf{K})^{-1}\bm{y}=-((1-\alpha)\mathbf{I}-\beta\mathbf{K})^{t+1}\beta(\alpha\mathbf{I}+\beta\mathbf{K})^{-1}\mathbf{K}\bm{y}.$
(26)

If there is no weight decay (i.e., $\alpha=0$), then by rearrangement of (A),
we obtain

$\displaystyle f_{t+1}(\mathbf{X})=$
$\displaystyle\left(\mathbf{I}-(\mathbf{I}-\beta\mathbf{K})^{t+1}\right)\bm{y}.$
(27)

The estimator after $t$-th iteration can be obtained by

$\displaystyle f_{t}(\bm{x})=$
$\displaystyle\bm{w}_{t}^{T}\mathbf{k}(\bm{x})=\mathbf{k}(\bm{x})^{T}\mathbf{K}^{-1}f_{t}(\mathbf{X}).$
(28)

Note that the kernel matrix $\mathbf{K}$ is generated by the empirical kernel
$K_{S}$ defined in (7). By taking the expectation with respect to
$\bm{\varepsilon}_{k_{1}}$ and $\bm{\varepsilon}_{k_{2}}$, we define the
expected smoothing kernel $\tilde{K}_{S}$ as

$\displaystyle\tilde{K}_{S}(\bm{x},\bm{x}^{\prime})=\int_{\mathbb{R}^{D}}\int_{\mathbb{R}^{D}}K(\bm{x}+\bm{\epsilon}-(\bm{x}^{\prime}+\bm{\epsilon}^{\prime}))p_{\varepsilon}(\bm{\epsilon})p_{\varepsilon}(\bm{\epsilon}^{\prime}){\rm
d}\bm{\epsilon}{\rm d}\bm{\epsilon}^{\prime}.$ (29)

Since $\tilde{K}_{S}$ is close to the empirical version of the smoothing
kernel $K_{S}$, we can consider the gradient flow with respect to the kernel
function $\tilde{K}_{S}$. The error analysis between $\tilde{K}_{S}$ and
$K_{S}$ is provided in Appendix B.

Let $g_{t}$ be the function obtained at $t$-th iteration by the gradient
update rule with respect to the kernel function $\tilde{K}_{S}$. Analogous to
(26) and (27), we have

$\displaystyle
g_{t}(\mathbf{X})=\tilde{\mathbf{K}}(\alpha/\beta\mathbf{I}+\tilde{\mathbf{K}})^{-1}\bm{y}-((1-\alpha)\mathbf{I}-\beta\tilde{\mathbf{K}})^{t}\beta(\alpha\mathbf{I}+\beta\tilde{\mathbf{K}})^{-1}\mathbf{K}\bm{y},$
(30)

if there is weight decay, and

$\displaystyle g_{t}(\mathbf{X})=$
$\displaystyle\left(\mathbf{I}-(\mathbf{I}-\beta\tilde{\mathbf{K}})^{t}\right)\bm{y},$
(31)

if there is no weight decay, where
$\tilde{\mathbf{K}}=(\tilde{K}_{S}(\bm{x}_{j}-\bm{x}_{k}))_{jk}$. Similarly,
the predictor of $f^{*}(\bm{x})$ using the kernel function $\tilde{K}$ can be
obtained by

$\displaystyle
g_{t}(\bm{x})=\tilde{\mathbf{k}}(\bm{x})^{T}\tilde{\mathbf{K}}^{-1}g_{t}(\mathbf{X}).$
(32)

Thus, the empirical error $\|f_{t}(\mathbf{X})-f^{*}(\mathbf{X})\|_{2}$ can be
decomposed by

$\displaystyle\|f_{t}(\mathbf{X})-f^{*}(\mathbf{X})\|_{2}\leq\|f_{t}(\mathbf{X})-g_{t}(\mathbf{X})\|_{2}+\|g_{t}(\mathbf{X})-f^{*}(\mathbf{X})\|_{2}.$
(33)

##  Appendix B Error of Data Augmentation

We first consider bounding the difference between the empirical smoothing
kernel function

$\displaystyle
K_{S}(\bm{x}-\bm{x}^{\prime})=\frac{1}{N^{2}}\sum_{k=1}^{N}\sum_{j=1}^{N}K(\bm{x}+\bm{\varepsilon}_{j}-(\bm{x}^{\prime}+\bm{\varepsilon}_{k}))$

and the expected smoothing kernel function

$\displaystyle\tilde{K}_{S}(\bm{x}-\bm{x}^{\prime})=\mathbb{E}_{\bm{\varepsilon},\bm{\varepsilon}^{\prime}}\big{(}K(\bm{x}+\bm{\varepsilon}-(\bm{x}^{\prime}+\bm{\varepsilon}^{\prime}))\big{)}=\int_{\mathbb{R}^{D}}\int_{\mathbb{R}^{D}}K(\bm{x}+\bm{\varepsilon}-(\bm{x}^{\prime}+\bm{\varepsilon}^{\prime}))p_{\varepsilon}(\bm{\varepsilon})p_{\varepsilon}(\bm{\varepsilon}^{\prime}){\rm
d}\bm{\varepsilon}{\rm d}\bm{\varepsilon}^{\prime}.$

Specifically, we have the following lemma.

######  Lemma B.1.

If Assumption 4.2 or 4.3, and Assumption 4.4 are satisfied , then

$\displaystyle\sup_{\bm{x},\bm{x}^{\prime}\in\Omega}\left|\mathbb{E}_{\bm{\varepsilon},\bm{\varepsilon}^{\prime}}\big{(}K(\bm{x}+\bm{\varepsilon}-(\bm{x}^{\prime}+\bm{\varepsilon}^{\prime}))\big{)}-\frac{1}{N^{2}}\sum_{k=1}^{N}\sum_{j=1}^{N}K(\bm{x}+\bm{\varepsilon}_{j}-(\bm{x}^{\prime}+\bm{\varepsilon}_{k}))\right|=O_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{N}}\right).$

Based on Lemma B.1, we can obtain an upper bound of
$\|f_{t}-g_{t}\|_{L_{\infty}(\Omega)}$ as follows. Recall that
$\mathbf{K}=(K_{S}(\bm{x}_{j}-\bm{x}_{k}))_{j,k=1}^{n}$,
$\tilde{\mathbf{K}}=(\tilde{K}_{S}(\bm{x}_{j}-\bm{x}_{k}))_{j,k=1}^{n}$. Let
$\eta_{1}(\mathbf{K})$ and $\eta_{n}(\mathbf{K})$ be the largest and smallest
eigenvalues of $\mathbf{K}$, respectively. Let $\eta_{n}(\tilde{\mathbf{K}})$
be the smallest eigenvalue of $\tilde{\mathbf{K}}$.

######  Lemma B.2.

Suppose Assumption 4.2 or 4.3, and Assumption 4.4 are satisfied. Furthermore,
assume that

$\displaystyle\frac{1}{2}\eta_{n}(\tilde{\mathbf{K}})\geq n\sqrt{\frac{\log
N}{N}},$ (34)

and the learning rate $\beta$ satisfies $\beta\eta_{1}(\mathbf{K})+\alpha<1$,
where $\alpha=0$ if there is no weight decay, and $\alpha>0$ if there is
weight decay. Then we have

$\displaystyle\sup_{t\geq
1}\|f_{t}-g_{t}\|_{L_{\infty}(\Omega)}=O_{\mathbb{P}}\left(\frac{n^{2}\sqrt{\log
N/N}}{\eta_{n}(\tilde{\mathbf{K}})^{2}}\right),$

where the probability is with respect to the augmentation $\bm{\varepsilon}$.

Since $\tilde{\mathbf{K}}$ and $\eta_{n}(\tilde{\mathbf{K}})$ are determined
by the data $(\bm{x}_{j},y_{j})$, $j=1,...,n$, the left-hand side of (34) is
not depending on $N$. Therefore, the condition (34) can be fulfilled if we add
sufficient augmentations. In the next lemma, we provide a more explicit lower
bound of $\eta_{n}(\tilde{\mathbf{K}})$ in (34) in terms of $\bm{x}_{j}$’s.

######  Lemma B.3.

Let $q_{\mathbf{X}}$ be the separation distance defined as

$q_{\mathbf{X}}=\frac{1}{2}\min_{j\neq k}\|\bm{x}_{j}-\bm{x}_{k}\|_{2}.$

The minimum eigenvalue of $\tilde{\mathbf{K}}$, denoted by
$\eta_{n}(\tilde{\mathbf{K}})$, is lower bounded as follows.

  1. 1.

if Assumption 4.2 and Assumption 4.4 (C1) are satisfied, then

$\eta_{n}(\tilde{\mathbf{K}})\geq
C_{1}\big{(}1+4M^{2}\big{)}^{-m_{0}}\big{(}1+4\sigma_{n}^{2}M^{2}\big{)}^{-m_{\varepsilon}}M^{D};$

  2. 2.

if Assumption 4.3 and Assumption 4.4 (C2) are satisfied, then

$\eta_{n}(\tilde{\mathbf{K}})\geq
C_{2}\big{(}1+4M^{2}\big{)}^{-m_{0}D}\big{(}1+4\sigma_{n}^{2}M^{2}\big{)}^{-m_{\varepsilon}D}M^{D};$

  3. 3.

if Assumption 4.2 and Assumption 4.4 (C3) are satisfied, then

$\eta_{n}(\tilde{\mathbf{K}})\geq
C_{3}\big{(}1+4M^{2}\big{)}^{-m_{0}}e^{-8\sigma_{n}^{2}M^{2}}M^{D}$

where $C_{i}$’s are constants only depending on $D$,
$M=\frac{12}{q_{\mathbf{X}}}\big{(}\frac{\pi\Gamma^{2}(\frac{D}{2}+1)}{9}\big{)}^{\frac{1}{D+1}}$,
and $\Gamma(\cdot)$ denotes the Gamma function.

The proofs of the above three lemmas are put in Appendix H.

##  Appendix C A Comparison Theorem

In this section, we provide a byproduct, which is a generic comparison theorem
between the early-stopping without weight decay and the kernel ridge
regression estimator. Let $K_{1}$ be a positive definite kernel function. The
kernel ridge regression is defined by

$\displaystyle\tilde{g}=\mathop{\mathrm{argmin}}_{f\in\mathcal{H}_{K_{1}}(\Omega)}\|f-\bm{y}\|_{n}^{2}+\lambda\|f\|_{\mathcal{H}_{K_{1}}(\Omega)}^{2},$
(35)

where $\bm{y}=(y_{1},...,y_{n})^{T}$, $y_{j}$’s are as in (1), and $\lambda>0$
is a regularization parameter. The main theorem in this subsection is as
follows.

######  Theorem C.1.

Let $(\beta t)^{-1}=n\lambda$. Suppose $\epsilon_{j}$’s are i.i.d. random
noise with mean zero and finite variance $\sigma_{\epsilon}^{2}$. Let
$\tilde{g}_{t}(\bm{x})=\bm{w}_{t}^{T}\mathbf{k}(\bm{x})$, which is similar to
$\hat{f}_{t}(\bm{x})$ in (10) but with $K_{1}$ instead of $K_{S}$ and with
update rule (9). Then there exists a constant $C>0$ such that

$\displaystyle\mathbb{E}\|\tilde{g}_{t}-f^{*}\|_{n}^{2}\leq
C\mathbb{E}\|\tilde{g}-f^{*}\|_{n}^{2},$ (36)

and

$\displaystyle\mathbb{E}\|\tilde{g}_{t}\|_{\mathcal{H}_{K_{1}}(\Omega)}^{2}\leq
2\mathbb{E}\|\tilde{g}\|_{\mathcal{H}_{K_{1}}(\Omega)}^{2},$ (37)

where the expectation is taken with respect to the noises $\epsilon_{j}$,
$j=1,...,n$.

Theorem C.1 states that the mean squared prediction error of the early-
stopping without weight decay is smaller than (at most the same as) that of
the kernel ridge regression estimator, up to a multiplicative constant. This
explains why the upper bounds on the early-stopping without weight decay and
the kernel ridge regression estimator derived in Raskutti et al., (2014) are
identical, in a more explicit way. Note that the conditions of Theorem C.1 are
quite mild. We do not assume any relationship between $f^{*}$ and
$\mathcal{H}_{K_{1}}(\Omega)$, and do not require any particular structure of
the RKHS $\mathcal{H}_{K_{1}}(\Omega)$. Furthermore, we do not impose any
conditions on $\lambda$, and we only require that $\epsilon_{j}$’s are i.i.d.
with finite variance (not necessarily sub-Gaussian and can be even heavy-
tailed).

It is worth noting that the complexity (i.e., the RKHS norm) of the early-
stopping without weight decay is also bounded by the complexity of the kernel
ridge regression estimator, up to a constant multiplier. Since the difference
between the empirical norm $\|\cdot\|_{n}$ and the $L_{2}$ norm depends on the
complexity of the estimator, it can be expected that (36) still holds if we
replace the empirical norm by the $L_{2}$ norm.

##  Appendix D Proof of Theorem 4.1

In this section, we show the proof of the following theorem. Note that the
second statement in Theorem D.1 is Theorem 4.1.

######  Theorem D.1 (Polynomial smoothing).

Suppose Assumptions 4.1, 4.2, 4.4 (C1), and 4.5 are satisfied. Suppose there
exists $\Omega_{1}$ with positive Lebesgue measure and a Lipschitz boundary
such that $\Omega\subset\Omega_{1}$ and
$f^{*}\in\mathcal{W}^{m_{f}}(\Omega_{1})$. Let $f_{t}(\bm{x})$ be as in (10)
and $\beta=n^{-1}C_{1}$ with $C_{1}\leq
2^{-1}\sup_{\bm{x}\in\mathbb{R}^{D}}K_{S}(\bm{x})$. Suppose the smoothing
scale $\sigma_{n}\asymp n^{\nu}$ with $\nu\leq 0$. Suppose one of the
following holds:

  1. 1.

There is no weight decay in the gradient descent, and the iteration number $t$
satisfies $t\asymp
n^{\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}\sigma_{n}^{2m_{\varepsilon}}$

  2. 2.

There is weight decay in the gradient descent with $\alpha\asymp
n^{-1-\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}\sigma_{n}^{-2m_{\varepsilon}}$,
and the iteration number satisfies $t\geq
C_{2}(\frac{m_{f}}{2m_{f}+d}+1/2)\log n/(\log(1-\alpha))$.

Then the following statements are true with $N>N_{0}$, where $N$ is the number
of augmentations, and $N_{0}$ depends on $n$ and the iteration number $t$.

  1. 1.

For any $a>0$, there exists an $m_{\varepsilon}$ such that when

$\displaystyle\nu=\left\\{\begin{array}[]{ll}-\frac{2(2m_{0}+2m_{\varepsilon})D-(2m_{0}+2m_{\varepsilon}-D)d}{(2m_{f}+d)(4m_{\varepsilon}D-(2m_{0}+2(1-d^{-1}(2m_{f}+d)a)m_{\varepsilon}-D)d)},&D>d,\\\
0,&D=d,\end{array}\right.$

we have

$\displaystyle\|f_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}^{2}=$ $\displaystyle
O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}+a}\right).$

  2. 2.

Set $m_{\varepsilon}=2d^{-1}(2D\max(m_{0},m_{f})+m_{0}d)\log n-m_{0}$. Then by
choosing

$\displaystyle\nu=\left\\{\begin{array}[]{ll}-\frac{2(2m_{0}+2m_{\varepsilon})D-(2m_{0}+2m_{\varepsilon}-D)d}{(2m_{f}+d)(4m_{\varepsilon}D-(2m_{0}+2(1-(\log
n)^{-1})m_{\varepsilon}-D)d)}<0,&D>d,\\\ 0,&D=d,\end{array}\right.$

we have

$\displaystyle\|f_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}^{2}=$ $\displaystyle
O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}}(\log n)^{2m_{f}+1}\right).$

We first present several lemmas used in this proof. The proof of these lemmas
can be found in Appendix I.

######  Lemma D.1.

Suppose the conditions of Theorem 4.1 are fulfilled. Let $f_{n}^{*}$ be the
solution to the optimization problem

$\displaystyle\min_{g\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|f^{*}-g\|_{L_{2}(P_{\mathbf{X}})}^{2}+\lambda_{n}\|g\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}.$
(38)

Then if $m_{0}\leq m_{f}$, we have

$\displaystyle\|f^{*}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}^{2}+\lambda_{n}\|f_{n}^{*}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}\leq
C_{1}\max\left((\lambda_{n}(m_{\varepsilon}+1)^{m_{\varepsilon}}\sigma_{n}^{2m_{\varepsilon}})^{\frac{m_{f}}{m_{0}+m_{\varepsilon}}},\lambda_{n}\right).$
(39)

and if $m_{0}>m_{f}$, we have

$\displaystyle\|f^{*}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}^{2}+\lambda_{n}\|f_{n}^{*}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}\leq
C_{2}\max\left((\lambda_{n}(m_{\varepsilon}+1)^{m_{\varepsilon}}\sigma_{n}^{2m_{\varepsilon}})^{\frac{m_{f}}{m_{0}+m_{\varepsilon}}},\lambda_{n}^{\frac{m_{f}}{m_{0}}}\right).$
(40)

Here the constants $C_{1}$ and $C_{2}$ are independent with $m_{\varepsilon}$.

######  Lemma D.2.

Suppose the conditions of Theorem 4.1 are fulfilled. Let $f_{n}^{*}$ be as in
Lemma D.1. Suppose there exists $T>0$ (depending on $n$) such that

$\displaystyle\|f^{*}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}^{2}+\lambda_{n}\|f_{n}^{*}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}\leq
T.$

Let $\hat{f}_{n}$ be the solution to the optimization problem

$\displaystyle\min_{g\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|\bm{y}-g\|_{n}^{2}+\lambda_{n}\|g\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2},$
(41)

where $\bm{y}=(y_{1},...,y_{n})^{T}$. Suppose

$\displaystyle\sigma_{n}^{-d/2}m^{\frac{mD}{2m-D}+\frac{1}{2}}\log p$

converges to zero as $n$ goes to infinity, where $p=\frac{4D}{2m-D}$, and
$m=m_{0}+m_{\varepsilon}$. Then we have

$\displaystyle M_{1}=$
$\displaystyle\max\bigg{(}(T+n^{-1/2}T^{1/2})^{1/2},\lambda_{n}^{-\frac{p}{2(4-p)}}\left(\sigma_{n}^{-d/2}n^{-1/2}m^{\frac{mD}{2m-D}+\frac{1}{2}}(T+n^{-1/2}T^{1/2})^{\frac{1}{2}-\frac{p}{4}}\right)^{\frac{2}{4-p}},$
$\displaystyle\sigma_{n}^{-d/2}n^{-1/2}m^{\frac{mD}{2m-D}+\frac{1}{2}}\lambda_{n}^{-\frac{p}{4}},\left(\sigma_{n}^{-d/2}n^{-1/2}m^{\frac{mD}{2m-D}+\frac{1}{2}}(\lambda_{n}^{-1}T)^{\frac{p}{2}}(T+n^{-1/2}T^{1/2})^{1-\frac{p}{2}}\right)^{1/2},$
$\displaystyle(\sigma_{n}^{-d/2}n^{-1/2}m^{\frac{mD}{2m-D}+\frac{1}{2}})^{\frac{2}{2+p}}(\lambda_{n}^{-1}T)^{\frac{p}{2(2+p)}}\bigg{)},$
$\displaystyle M_{2}=$
$\displaystyle\max\bigg{(}(\lambda_{n}^{-1}(T+n^{-1/2}T^{1/2}))^{1/2},\left(\lambda_{n}^{-1}\sigma_{n}^{-d/2}n^{-1/2}m^{\frac{mD}{2m-D}+\frac{1}{2}}(T+n^{-1/2}T^{1/2})^{\frac{1}{2}-\frac{p}{4}}\right)^{\frac{2}{4-p}},$
$\displaystyle\sigma_{n}^{-d/2}n^{-1/2}m^{\frac{mD}{2m-D}+\frac{1}{2}}\lambda_{n}^{-\frac{2+p}{4}},\left(\lambda_{n}^{-1}\sigma_{n}^{-d/2}n^{-1/2}m^{\frac{mD}{2m-D}+\frac{1}{2}}(\lambda_{n}^{-1}T)^{\frac{p}{2}}(T+n^{-1/2}T^{1/2})^{1-\frac{p}{2}}\right)^{1/2},$
$\displaystyle\lambda_{n}^{-1/2}(\sigma_{n}^{-d/2}n^{-1/2}m^{\frac{mD}{2m-D}+\frac{1}{2}})^{\frac{2}{2+p}}(\lambda_{n}^{-1}T)^{\frac{p}{2(2+p)}}\bigg{)},$

Then we have

$\displaystyle\|f^{*}-\hat{f}_{n}\|_{n}=O_{\mathbb{P}}(M_{1}),\|\hat{f}_{n}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}=O_{\mathbb{P}}(M_{2}).$

Furthermore, if $\tilde{f}_{n}$ be the solution to the optimization problem

$\displaystyle\min_{f\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|f^{*}-f\|_{n}^{2}+\lambda_{n}\|f\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)},$
(42)

then

$\displaystyle\|f^{*}-\tilde{f}_{n}\|_{n}=O_{\mathbb{P}}((T+n^{-1/2}T^{1/2})^{1/2}),\|\tilde{f}_{n}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}=O_{\mathbb{P}}((\lambda_{n}^{-1}(T+n^{-1/2}T^{1/2}))^{1/2}).$

######  Lemma D.3 (Lemma F.5 of Wang, (2021)).

Assume for class $\mathcal{G}$,
$\sup_{g\in\mathcal{G}}\|g\|_{L_{\infty}(\Omega)}\leq c<1$, and the bracket
entropy
$H_{B}(\delta_{n},\mathcal{G},\|\cdot\|_{L_{2}(P_{\mathbf{X}})})\leq\frac{n\delta_{n}^{2}}{1200c^{2}}$,
and $n\delta_{n}^{2}\rightarrow\infty$, where $0<\delta_{n}<1$. Then we have

$\displaystyle P\bigg{(}\inf_{\|g\|_{L_{2}(P_{\mathbf{X}})}\geq
2\delta_{n},g\in\mathcal{G}}\frac{\|g\|^{2}_{n}}{\|g\|_{L_{2}(P_{\mathbf{X}})}^{2}}<C_{3}\bigg{)}\leq
C_{5}\exp(-C_{6}n\delta_{n}^{2}/c^{2}),$

and

$\displaystyle P\bigg{(}\sup_{\|g\|_{L_{2}(P_{\mathbf{X}})}\geq
2\delta_{n},g\in\mathcal{G}}\frac{\|g\|^{2}_{n}}{\|g\|_{L_{2}(P_{\mathbf{X}})}^{2}}>C_{4}\bigg{)}\leq
C_{7}\exp(-C_{8}n\delta_{n}^{2}/c^{2}),$

for some constants $C_{3},C_{4}>0$ and $C_{i}$’s ($i=5,6,7,8$) are only
depending on $\Omega$.

######  Lemma D.4 (Interpolation inequality for Polynomial RKHS).

Let $g\in\mathcal{W}^{m}(\mathbb{R}^{D})$. When
$r=\frac{D}{2(m_{0}+m_{\varepsilon})}$ and $D>1$, we have

$\displaystyle\|g\|_{L_{\infty}(\mathbb{R}^{D})}\leq
C_{9}\|g\|_{L_{2}(\mathbb{R}^{D})}^{1-r}\|g\|_{\mathcal{W}^{m}(\mathbb{R}^{D})}^{r},$

where
$C_{9}=\left(\int_{\mathbb{R}^{D}}(1+\|\bm{\omega}\|_{2}^{2})^{-\frac{D}{2}}d\bm{\omega}\right)^{\frac{1}{2}}<\infty$.

###  D.1 Without weight decay

By the triangle inequality, it can be seen that

$\displaystyle\|f_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}\leq$
$\displaystyle\|f_{t}-g_{t}\|_{L_{2}(P_{\mathbf{X}})}+\|g_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})},$
(43)

where $g_{t}$ is as in (31).

By Lemma B.2, the first term $\|f_{t}-g_{t}\|_{L_{2}(P_{\mathbf{X}})}$ in (43)
can be bounded by

$\displaystyle\|f_{t}-g_{t}\|_{L_{2}(P_{\mathbf{X}})}\leq
C_{10}\|f_{t}-g_{t}\|_{L_{\infty}(\Omega)}=O_{\mathbb{P}}\left(\frac{n^{2}\sqrt{\log
N/N}}{\eta_{n}(\tilde{\mathbf{K}})^{2}}\right),$

as long as

$\displaystyle\frac{1}{2}\eta_{n}(\tilde{\mathbf{K}})\geq n\sqrt{\frac{\log
N}{N}}.$ (44)

Choose

$\displaystyle N_{0}=\frac{4n^{2}}{\eta_{n}(\tilde{\mathbf{K}})^{2}}.$ (45)

Then it holds that when $N\geq N_{0}$,

$\displaystyle\|f_{t}-g_{t}\|_{L_{2}(P_{\mathbf{X}})}=O_{\mathbb{P}}\left(n^{-1/2}\right).$
(46)

It remains to consider $\|g_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}$ in (43). In
order to do so, we consider the empirical version of
$\|g_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}$, and let

$\displaystyle
J_{2}=\|g_{t}-f^{*}\|_{n}^{2}=\frac{1}{n}\|g_{t}(\mathbf{X})-f^{*}(\mathbf{X})\|_{2}^{2}.$
(47)

Let $(\beta t)^{-1}=n\lambda_{n}$. Consider the kernel ridge regression

$\displaystyle\tilde{g}=\mathop{\mathrm{argmin}}_{f\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|f-\bm{y}\|_{n}^{2}+\lambda_{n}\|f\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}.$
(48)

By the representer theorem,
$\tilde{g}(\bm{x})=\tilde{\mathbf{k}}(\bm{x})^{T}(\tilde{\mathbf{K}}+n\lambda_{n}\mathbf{I})^{-1}\bm{y}$
for all $\bm{x}\in\Omega$, where
$\tilde{\mathbf{k}}(\bm{x})=(\tilde{K}_{S}(\bm{x}-\bm{x}_{1}),...,\tilde{K}_{S}(\bm{x}-\bm{x}_{n}))^{T}$.
Then it can be seen that

$\displaystyle\tilde{g}(\mathbf{X})-f^{*}(\mathbf{X})=n\lambda_{n}(\tilde{\mathbf{K}}+n\lambda_{n}\mathbf{I})^{-1}f^{*}(\mathbf{X})+\tilde{\mathbf{K}}(\tilde{\mathbf{K}}+n\lambda_{n}\mathbf{I})^{-1}\bm{\epsilon}=\bm{q}_{1}+\bm{q}_{2}$

Recall that (see (31))

$\displaystyle g_{t}(\mathbf{X})=$
$\displaystyle\left(\mathbf{I}-(\mathbf{I}-\beta\tilde{\mathbf{K}})^{t}\right)\bm{y},$

which implies

$\displaystyle
g_{t}(\mathbf{X})-f^{*}(\mathbf{X})=-(\mathbf{I}-\beta\tilde{\mathbf{K}})^{t}f^{*}(\mathbf{X})+\left(\mathbf{I}-(\mathbf{I}-\beta\tilde{\mathbf{K}})^{t}\right)\bm{\epsilon}.$
(49)

By the Cauchy-Schwarz inequality, (47), and (49), it can be seen that

$\displaystyle nJ_{2}\leq$ $\displaystyle
2(f^{*}(\mathbf{X}))^{T}(\mathbf{I}-\beta\tilde{\mathbf{K}})^{2t}f^{*}(\mathbf{X})+2\bm{\epsilon}^{T}(\mathbf{I}-(\mathbf{I}-\beta\tilde{\mathbf{K}})^{t})^{2}\bm{\epsilon}$
$\displaystyle=$ $\displaystyle 2nJ_{21}+2nJ_{22},$ (50)

and

$\displaystyle n\|\tilde{g}-f^{*}\|_{n}^{2}\leq$ $\displaystyle
2(n\lambda_{n})^{2}(f^{*}(\mathbf{X}))^{T}(\tilde{\mathbf{K}}+n\lambda_{n}\mathbf{I})^{-2}f^{*}(\mathbf{X})+2\bm{\epsilon}^{T}(\tilde{\mathbf{K}}+n\lambda_{n}\mathbf{I})^{-1}\tilde{\mathbf{K}}^{2}(\tilde{\mathbf{K}}+n\lambda_{n}\mathbf{I})^{-1}\bm{\epsilon}$
$\displaystyle=$ $\displaystyle
2\|\bm{q}_{1}\|_{2}^{2}+2\|\bm{q}_{2}\|_{2}^{2}.$ (51)

Similar to (G), it can be seen that

$\displaystyle 2nJ_{21}\leq C_{11}\|\bm{q}_{1}\|_{2}^{2},$ (52)

and similar to (G), the term $2nJ_{22}$ can be further bounded by

$\displaystyle 2nJ_{22}=$ $\displaystyle
2\sum_{j=1}^{n}(1-(1-\beta\eta_{j})^{t})^{2}(\bm{v}_{j}^{T}\bm{\epsilon})^{2}\leq
2\sum_{j=1}^{n}\frac{4(\beta t\eta_{j})^{2}}{(1+\beta
t\eta_{j})^{2}}(\bm{v}_{j}^{T}\bm{\epsilon})^{2}$ $\displaystyle=$
$\displaystyle 8\bm{\epsilon}^{T}(\tilde{\mathbf{K}}+(\beta
t)^{-1}\mathbf{I})^{-1}\tilde{\mathbf{K}}^{2}(\tilde{\mathbf{K}}+(\beta
t)^{-1}\mathbf{I})^{-1}\bm{\epsilon}=8\|\bm{q}_{2}\|_{2}^{2},$ (53)

where $\eta_{1}\geq\ldots\geq\eta_{n}>0$ and $\bm{v}_{j}$, $j=1,\ldots,n$ be
the eigenvalues and corresponding eigenvectors of $\tilde{\mathbf{K}}$,
respectively. In the last inequality of (D.1), we note $(\beta
t)^{-1}=n\lambda_{n}$.

Plugging (52) and (D.1) into (D.1), we obtain

$\displaystyle
J_{2}\leq\frac{2C_{12}}{n}\left(\|\bm{q}_{1}\|_{2}^{2}+\|\bm{q}_{2}\|_{2}^{2}\right).$
(54)

The term $\|\bm{q}_{1}\|_{2}^{2}$ and $\|\bm{q}_{2}\|_{2}^{2}$ can be directly
bounded by Lemma D.2. To see this, let $f_{0}(\bm{x})=0$ for all
$\bm{x}\in\Omega$. Then it can be checked that

$\displaystyle\frac{1}{n}\|\bm{q}_{1}\|_{2}^{2}=$
$\displaystyle\|\tilde{f}_{n}-f\|_{n}^{2}$

and

$\displaystyle\frac{1}{n}\|\bm{q}_{2}\|_{2}^{2}=$
$\displaystyle\|\hat{f}_{0,n}-f_{0}\|_{n}^{2},$

where $\tilde{f}_{n}$ is as in (42), and $\hat{f}_{0,n}$ is the solution to
the optimization problem

$\displaystyle\min_{g\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|\bm{\epsilon}-g\|_{n}^{2}+\lambda_{n}\|g\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}.$

Let $\delta_{0}\in(0,1)$ such that
$4m_{\varepsilon}D-(2m_{0}+2(1-\delta_{0})m_{\varepsilon}-D)d>0$. Take

$\displaystyle\lambda_{n}\asymp
n^{-\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}\sigma_{n}^{-2m_{\varepsilon}},\sigma_{n}\asymp
n^{-\frac{2(2m_{0}+2m_{\varepsilon})D-(2m_{0}+2m_{\varepsilon}-D)d}{(2m_{f}+d)(4m_{\varepsilon}D-(2m_{0}+2(1-\delta_{0})m_{\varepsilon}-D)d)}},n^{-1}(\beta
t)^{-1}\asymp\lambda_{n},\beta\asymp n^{-1}.$

Therefore, if $m_{\varepsilon}=O((\log n)^{C})$ for some constant $C$, and

$\displaystyle\lambda_{n}\leq
C_{13}(\lambda_{n}(m_{\varepsilon}+1)^{m_{\varepsilon}}\sigma_{n}^{2m_{\varepsilon}})^{\frac{m_{f}}{m_{0}+m_{\varepsilon}}}$
$\displaystyle\Leftrightarrow$ $\displaystyle
n^{-\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}n^{\frac{4m_{\varepsilon}(2m_{0}+2m_{\varepsilon})D-2m_{\varepsilon}(2m_{0}+2m_{\varepsilon}-D)d}{(2m_{f}+d)(4m_{\varepsilon}D-(2m_{0}+2(1-\delta_{0})m_{\varepsilon}-D)d)}}\leq
C_{14}n^{-\frac{2m_{f}}{2m_{f}+d}}(m_{\varepsilon}+1)^{\frac{m_{\varepsilon}m_{f}}{m_{0}+m_{\varepsilon}}}$
$\displaystyle\Leftarrow$ $\displaystyle
m_{\varepsilon}^{2}\delta_{0}d>m_{\varepsilon}(2m_{f}D+(m_{0}-m_{f})(1-\delta_{0})d)$
$\displaystyle\Leftarrow$ $\displaystyle
m_{\varepsilon}>\frac{2m_{f}D+m_{0}d}{\delta_{0}d},$ (55)

when $m_{0}\leq m_{f}$, or

$\displaystyle\lambda_{n}^{\frac{m_{f}}{m_{0}}}\leq
C_{15}(\lambda_{n}(m_{\varepsilon}+1)^{m_{\varepsilon}}\sigma_{n}^{2m_{\varepsilon}})^{\frac{m_{f}}{m_{0}+m_{\varepsilon}}}$
$\displaystyle\Leftrightarrow$ $\displaystyle
n^{-\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}n^{\frac{4m_{\varepsilon}(2m_{0}+2m_{\varepsilon})D-2m_{\varepsilon}(2m_{0}+2m_{\varepsilon}-D)d}{(2m_{f}+d)(4m_{\varepsilon}D-(2m_{0}+2(1-\delta_{0})m_{\varepsilon}-D)d)}}\leq
C_{16}n^{-\frac{2m_{0}}{2m_{f}+d}}(m_{\varepsilon}+1)^{\frac{m_{\varepsilon}m_{0}}{m_{0}+m_{\varepsilon}}}$
$\displaystyle\Leftarrow$ $\displaystyle
m_{\varepsilon}^{2}\delta_{0}d>2m_{0}m_{\varepsilon}D$
$\displaystyle\Leftarrow$ $\displaystyle
m_{\varepsilon}>\frac{2m_{0}D+m_{0}d}{\delta_{0}d},$ (56)

when $m_{0}>m_{f}$, we have

$\displaystyle T\leq
C_{17}n^{-\frac{2m_{f}}{2m_{f}+d}}(m_{\varepsilon}+1)^{\frac{m_{\varepsilon}m_{f}}{m_{0}+m_{\varepsilon}}}\leq
C_{17}n^{-\frac{2m_{f}}{2m_{f}+d}}(m_{\varepsilon}+1)^{m_{f}},$

where $T$ is as in Lemma D.2. Suppose $D>1$, long but tedious calculation
shows that

$\displaystyle M_{1}\leq
C_{18}(m_{\varepsilon}+m_{0})^{m_{f}+\frac{1}{2}}n^{-\frac{m_{f}}{2m_{f}+d}+\delta^{\prime}},$

where $M_{1}$ is as in Lemma D.2, and

$\displaystyle\delta^{\prime}=$
$\displaystyle\frac{((4m_{0}+4m_{\varepsilon})D-(2m_{0}+2m_{\varepsilon}-D)d)m_{\varepsilon}d}{(2m_{f}+d)(2m_{\varepsilon}+2m_{0}-D)(4m_{\varepsilon}D-(2m_{0}+2(1-\delta_{0})m_{\varepsilon}-D)d)}\delta_{0}\leq\frac{d}{2(2m_{f}+d)}\delta_{0},$

where the inequality is because of (D.1) (if $m_{0}\leq m_{f}$) or (D.1) (if
$m_{0}>m_{f}$). Therefore, by taking $\delta_{0}=d^{-1}(2m_{f}+d)a$ and
$m_{\varepsilon}=(\delta_{0}d)^{-1}(2D\max(m_{0},m_{f})+m_{0}d)+1$, we have

$\displaystyle\frac{1}{n}\|\bm{q}_{1}\|_{2}^{2}=$
$\displaystyle\|\tilde{f}_{n}-f\|_{n}^{2}=O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}+a}\right),$
$\displaystyle\frac{1}{n}\|\bm{q}_{2}\|_{2}^{2}=$
$\displaystyle\|\hat{f}_{0,n}-f_{0}\|_{n}^{2}=O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}+a}\right).$
(57)

Then by (54) and (D.1), we obtain

$\displaystyle
J_{2}=O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}+a}\right),$ (58)

which corresponds to the first statement of Theorem 4.1.

Taking $\delta_{0}=(\log n)^{-1}$, we obtain that

$\displaystyle M_{1}\leq
C_{18}n^{-\frac{m_{f}}{2m_{f}+d}}e^{\frac{d}{2(2m_{f}+d)}}(m_{\varepsilon}+m_{0})^{m_{f}+\frac{1}{2}}\leq
C_{19}n^{-\frac{m_{f}}{2m_{f}+d}}(m_{\varepsilon}+m_{0})^{m_{f}+\frac{1}{2}},$

where we require $m_{\varepsilon}>d^{-1}(2D\max(m_{0},m_{f})+m_{0}d)\log n$.
Thus, we can directly take
$m_{\varepsilon}=2d^{-1}(2D\max(m_{0},m_{f})+m_{0}d)\log n-m_{0}$ such that

$\displaystyle\frac{1}{n}\|\bm{q}_{1}\|_{2}^{2}=$
$\displaystyle\|\tilde{f}_{n}-f\|_{n}^{2}=O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}}(\log
n)^{2m_{f}+1}\right),$ $\displaystyle\frac{1}{n}\|\bm{q}_{2}\|_{2}^{2}=$
$\displaystyle\|\hat{f}_{0,n}-f_{0}\|_{n}^{2}=O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}}(\log
n)^{2m_{f}+1}\right).$ (59)

Thus, by (54) and (D.1), we have

$\displaystyle J_{2}=O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}}(\log
n)^{2m_{f}+1}\right),$ (60)

which corresponds to the second statement of Theorem 4.1.

It remains to bound $\|g_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}$. Note that

$\displaystyle\|g_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}\leq\|g_{t}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}+\|f_{n}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}\leq\|g_{t}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}+T^{1/2},$

and

$\displaystyle\|g_{t}-f_{n}^{*}\|_{n}\leq$
$\displaystyle\|g_{t}-f^{*}\|_{n}+\|f_{n}^{*}-f^{*}\|_{n}\leq\|g_{t}-f^{*}\|_{n}+O_{\mathbb{P}}\left(\left(T+n^{-1/2}T^{1/2}\right)^{1/2}\right)$
$\displaystyle\leq$ $\displaystyle
O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}}(\log n)^{2m_{f}+1}\right),$

where the second inequality is because of (139). Therefore, it suffices to
bound the difference between $\|g_{t}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}$ and
$\|g_{t}-f_{n}^{*}\|_{n}$.

By (G) and Lemma D.2, we have

$\displaystyle\|g_{t}\|_{\mathcal{N}_{\sigma}(\Omega)}^{2}\leq\sigma_{n}^{-2m_{0}}\|g_{t}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}\leq
C_{20}\sigma_{n}^{-2m_{0}}\|\tilde{g}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}=O_{\mathbb{P}}\left(n^{\nu_{1}}(\log
n)^{2m_{f}+1}\right),$ (61)

where

$\displaystyle\nu_{1}=$
$\displaystyle\frac{2(m_{0}+m_{\varepsilon}-m_{f})}{2m_{f}+d}+2(m_{\varepsilon}-m_{0})\nu,$
$\displaystyle\mbox{and }\nu=$
$\displaystyle-\frac{2(2m_{0}+2m_{\varepsilon})D-(2m_{0}+2m_{\varepsilon}-D)d}{(2m_{f}+d)(4m_{\varepsilon}D-(2m_{0}+2(1-\delta_{0})m_{\varepsilon}-D)d)}.$
(62)

Consider function class
$\mathcal{G}=\\{h:h=(g_{t}-f_{n}^{*})/(C_{21}n^{\nu_{1}/2}(\log
n)^{m_{f}+1/2})\\}$, where the constant $C_{21}$ is taken such that
$\|h_{1}\|_{\mathcal{N}_{\sigma}(\Omega)}<1$ for all $h_{1}\in\mathcal{G}$.
Then lemma D.4 leads to

$\displaystyle\|h_{1}\|_{L_{\infty}(\Omega)}\leq
C_{22}\|h_{1}\|_{L_{2}(P_{\mathbf{X}})}^{1-\frac{D}{2(m_{0}+m_{\varepsilon})}}\|h_{1}\|_{\mathcal{N}_{\sigma}(\Omega)}^{\frac{D}{2(m_{0}+m_{\varepsilon})}},$

for all $h_{1}\in\mathcal{G}$, which implies

$\displaystyle
c_{1}:=\sup_{h_{1}\in\mathcal{G}}\|h_{1}\|_{L_{\infty}(\Omega)}\leq
C_{22}R_{1}^{1-\frac{D}{2(m_{0}+m_{\varepsilon})}},$

where
$R_{1}=\sup_{h_{1}\in\mathcal{G}}\|h_{1}\|_{L_{2}(P_{\mathbf{X}})}\leq\sup_{h_{1}\in\mathcal{G}}\|h_{1}\|_{L_{\infty}(\Omega)}\leq\sup_{h_{1}\in\mathcal{G}}\|h_{1}\|_{\mathcal{N}_{\sigma}(\Omega)}<1$,
because of the reproducing property. Let $m=m_{0}+m_{\varepsilon}$. Taking
$c=C_{22}R_{1}^{1-\frac{D}{2m}}<1$, and
$\delta_{n}=C_{23}(\sigma_{n}^{-d}n^{-1}c^{2}m^{\frac{2mD}{2m-D}})^{\frac{2m-D}{4m}}$
in Lemma D.3, it can be checked that

$\displaystyle C_{24}n\delta_{n}^{2}c^{-2}\geq
H(\delta,\mathcal{B}_{\mathcal{H}_{\sigma}(\Omega)},\|\cdot\|_{L_{\infty}(\Omega)}),$

which implies the conditions of Lemma D.3 are fulfilled. Applying Lemma D.3 to
the case
$\|g_{t}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}^{2}\geq\delta_{n}^{2}n^{\nu_{1}}$,
together with (60), we have

$\displaystyle
R_{1}=O_{\mathbb{P}}\left(\max\\{n^{-\frac{m_{f}}{2m_{f}+d}-\nu_{1}/2}(\log
n)^{m_{f}+1/2},\delta_{n}\\}\right).$ (63)

If $\delta_{n}\geq n^{-\frac{m_{f}}{2m_{f}+d}-\nu_{1}/2}(\log n)^{m_{f}+1/2}$,
we have $R_{1}\leq C_{25}\delta_{n}$, which implies

$\displaystyle R_{1}\leq
C_{26}(\sigma_{n}^{-d}n^{-1}c^{2}m^{\frac{2mD}{2m-D}})^{\frac{2m-D}{4m}}.$

Therefore, we have

$\displaystyle\|g_{t}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}\leq
C_{21}n^{\nu_{1}/2}R_{1}\leq C_{27}n^{\nu_{2}}(\log n)^{D/2},$

where

$\displaystyle\nu_{2}=$
$\displaystyle\frac{(m_{0}+m_{\varepsilon}-m_{f})}{2m_{f}+d}+(m_{\varepsilon}-m_{0})\nu-\frac{2m-D}{4m}(d\nu+1)<-\frac{m_{f}}{2m_{f}+d}.$

If $\delta_{n}<n^{-\frac{m_{f}}{2m_{f}+d}-\nu_{1}/2}(\log n)^{m_{f}+1/2}$,
then $R_{1}=O_{\mathbb{P}}(n^{-\frac{m_{f}}{2m_{f}+d}-\nu_{1}/2}(\log
n)^{m_{f}+1/2})$, which implies
$\|g_{t}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}=O_{\mathbb{P}}(n^{-\frac{m_{f}}{2m_{f}+d}}(\log
n)^{m_{f}+1/2})$. Here we note that the proof is still valid if we replace
$g_{t}$ with $\tilde{g}$. Therefore, in both cases we have
$\|g_{t}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}=O_{\mathbb{P}}(n^{-\frac{m_{f}}{2m_{f}+d}}(\log
n)^{m_{f}+1/2})$, which, together with (D.1) and (46), finishes the proof.

###  D.2 With weight decay

If $\alpha>0$, we decompose the error by

$\displaystyle\|f_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}\leq$
$\displaystyle\|f_{t}-g_{t}\|_{L_{2}(P_{\mathbf{X}})}+\|\tilde{\mathbf{k}}(\cdot)^{T}(\alpha/\beta\mathbf{I}+\tilde{\mathbf{K}})^{-1}\bm{y}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}$
$\displaystyle+\|\beta\mathbf{k}(\cdot)^{T}((1-\alpha)\mathbf{I}-\beta\tilde{\mathbf{K}})^{t}(\alpha\mathbf{I}+\beta\tilde{\mathbf{K}})^{-1}\bm{y}\|_{L_{2}(P_{\mathbf{X}})}$
$\displaystyle=$ $\displaystyle I_{1}+I_{2}+I_{3}.$ (64)

As in (46), there exists an $N_{0}$ (depending on $n$) such that when $N\geq
N_{0}$,

$\displaystyle I_{1}=O_{\mathbb{P}}\left(n^{-1/2}\right).$ (65)

The second term is the error
$\|\tilde{f}_{n}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}$, where $\tilde{f}_{n}$ is as
in (42). Lemma D.2 gives us that

$\displaystyle\|\tilde{f}_{n}-f^{*}\|_{n}=O_{\mathbb{P}}(n^{-\frac{m_{f}}{2m_{f}+d}}).$

Following a similar approach in Appendix D.1, it can be further shown that

$\displaystyle I_{2}=O_{\mathbb{P}}(n^{-\frac{m_{f}}{2m_{f}+d}}),$ (66)

where we let $\alpha\asymp
n^{-1-\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}\sigma_{n}^{-2m_{\varepsilon}}$,
and $\beta$ and $\sigma_{n}$ are as in Theorem 4.1.

It remains to bound $I_{3}$ in (D.2). By Cauchy-Schwarz inequality,

$\displaystyle\|\beta\mathbf{k}(\cdot)^{T}((1-\alpha)\mathbf{I}-\beta\tilde{\mathbf{K}})^{t}(\alpha\mathbf{I}+\beta\tilde{\mathbf{K}})^{-1}\bm{y}\|_{L_{2}(P_{\mathbf{X}})}$
$\displaystyle\leq$ $\displaystyle\left\|\left({\rm
tr}\left(\left((\alpha/\beta\mathbf{I}+\tilde{\mathbf{K}})^{-1}\bm{y}\mathbf{k}(\cdot)^{T}\right)^{2}\right){\rm
tr}\left(((1-\alpha)\mathbf{I}-\beta\tilde{\mathbf{K}})^{2t}\right)\right)^{1/2}\right\|_{L_{2}(P_{\mathbf{X}})}$
$\displaystyle\leq$
$\displaystyle\left\|\mathbf{k}(\cdot)^{T}(\alpha/\beta\mathbf{I}+\tilde{\mathbf{K}})^{-1}\bm{y}\right\|_{L_{2}(P_{\mathbf{X}})}\left({\rm
tr}\left(((1-\alpha)\mathbf{I}-\beta\tilde{\mathbf{K}})^{2t})\right)\right)^{1/2}$
$\displaystyle\leq$
$\displaystyle\|(\mathbf{k}(\cdot)^{T}\mathbf{k}(\cdot))^{1/2}\|_{L_{2}(P_{\mathbf{X}})}\|\bm{y}\|_{2}\beta/\alpha$
$\displaystyle=$ $\displaystyle
O_{\mathbb{P}}\left(n^{1+\frac{2(m_{0}+m_{\varepsilon})}{2m_{f}+d}}\sigma_{n}^{2m_{\varepsilon}}(1-\alpha)^{t}\right).$
(67)

Thus, there exists $t_{0}>0$ such that as long as $t>t_{0}$, $I_{2}$ dominates
$I_{3}$. Combining (65), (66), and (D.2), we finish the proof.

##  Appendix E Proof of Theorem 4.2

We first present some lemmas, whose proofs can be found in Appendix J.

######  Lemma E.1.

Let $k_{\sigma}(\bm{x}-\bm{x}^{\prime})$ be a Gaussian kernel defined by

$\displaystyle
k_{\sigma}(\bm{x}-\bm{x}^{\prime})=\exp\left(-\frac{\|\bm{x}-\bm{x}^{\prime}\|_{2}^{2}}{4\sigma^{2}}\right),$
(68)

and $\mathcal{H}_{\sigma}(\mathbb{R}^{D})$ be the RKHS generated by
$k_{\sigma}(\bm{x}-\bm{x}^{\prime})$. Then we have

$\displaystyle\|h_{1}\|_{\mathcal{H}_{\sigma_{n}/\sqrt{2}}(\mathbb{R}^{D})}\leq
C_{1}\sigma_{n}^{-D/2}\|h_{1}\|_{\mathcal{H}_{\tilde{K}_{S}}(\mathbb{R}^{D})}$
(69)

and

$\displaystyle\|h_{2}\|_{\mathcal{H}_{\tilde{K}_{S}}(\mathbb{R}^{D})}\leq
C_{2}\sigma_{n}^{-m_{0}-D/2}\|h_{2}\|_{\mathcal{H}_{\sqrt{3}\sigma_{n}}(\mathbb{R}^{D})},$
(70)

for $h_{1}\in\mathcal{H}_{\tilde{K}_{S}}(\mathbb{R}^{D})$ and
$h_{2}\in\mathcal{H}_{\sqrt{3}\sigma_{n}}(\mathbb{R}^{D})$, where $C_{1}$ and
$C_{2}$ does not depend on $\sigma_{n}$.

######  Lemma E.2.

Let $f_{n}^{*}$ be the solution to the optimization problem

$\displaystyle\min_{g\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|f^{*}-g\|_{L_{2}(P_{\mathbf{X}})}^{2}+\lambda_{n}\|g\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}.$
(71)

Then

$\displaystyle\|f^{*}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}^{2}\leq
C_{3}\max(\lambda_{n}\sigma_{n}^{-2m_{0}},\sigma_{n}^{2m_{f}}),$

and

$\displaystyle\|f_{n}^{*}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}\leq
C_{3}\lambda_{n}^{-1}\max(\lambda_{n}\sigma_{n}^{-2m_{0}},\sigma_{n}^{2m_{f}}).$

######  Lemma E.3.

Let $f_{n}^{*}$ be the solution to the optimization problem

$\displaystyle\min_{g\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|f^{*}-g\|_{L_{2}(P_{\mathbf{X}})}^{2}+\lambda_{n}\|g\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}.$
(72)

Suppose there exists $T>0$ (depending on $n$) such that

$\displaystyle\|f^{*}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}^{2}+\lambda_{n}\|f_{n}^{*}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}\leq
T.$

Let $\hat{f}_{n}$ be the solution to the optimization problem

$\displaystyle\min_{g\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|\bm{y}-g\|_{n}^{2}+\lambda_{n}\|g\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}.$
(73)

Let $p=(\log n)^{-1}$,

$\displaystyle M_{1}=$
$\displaystyle\max\bigg{(}(T+n^{-1/2}T^{1/2})^{1/2},\sigma_{n}^{-d/2-\frac{pD}{4}}p^{-(D+1)/2}n^{-1/2}\lambda_{n}^{-\frac{p}{4}},$
$\displaystyle\lambda_{n}^{-\frac{p}{2(4-p)}}\left(\sigma_{n}^{-d/2-\frac{pD}{4}}p^{-(D+1)/2}n^{-1/2}(T+n^{-1/2}T^{1/2})^{\frac{1}{2}-\frac{p}{4}}\right)^{\frac{2}{4-p}},$
$\displaystyle\left(\sigma_{n}^{-d/2-\frac{pD}{4}}p^{-(D+1)/2}n^{-1/2}(\lambda_{n}^{-1}T)^{\frac{p}{2}}(T+n^{-1/2}T^{1/2})^{1-\frac{p}{2}}\right)^{1/2},$
$\displaystyle(\sigma_{n}^{-d/2-\frac{pD}{4}}p^{-(D+1)/2}n^{-1/2})^{\frac{2}{2+p}}(\lambda_{n}^{-1}T)^{\frac{p}{2+p}}\bigg{)},$
$\displaystyle M_{2}=$
$\displaystyle\max\bigg{(}(\lambda_{n}^{-1}(T+n^{-1/2}T^{1/2}))^{1/2},\sigma_{n}^{-d/2-\frac{pD}{4}}p^{-(D+1)/2}n^{-1/2}\lambda_{n}^{-\frac{2+p}{4}},$
$\displaystyle\left(\lambda_{n}^{-1}\sigma_{n}^{-d/2-\frac{pD}{4}}p^{-(D+1)/2}n^{-1/2}(T+n^{-1/2}T^{1/2})^{\frac{1}{2}-\frac{p}{4}}\right)^{\frac{2}{4-p}},$
$\displaystyle\left(\lambda_{n}^{-1}\sigma_{n}^{-d/2-\frac{pD}{4}}p^{-(D+1)/2}n^{-1/2}(\lambda_{n}^{-1}T)^{\frac{p}{2}}(T+n^{-1/2}T^{1/2})^{1-\frac{p}{2}}\right)^{1/2},$
$\displaystyle\lambda_{n}^{-1/2}(\sigma_{n}^{-d/2-\frac{pD}{4}}p^{-(D+1)/2}n^{-1/2})^{\frac{2}{2+p}}(\lambda_{n}^{-1}T)^{\frac{p}{2+p}}\bigg{)}.$

Then we have

$\displaystyle\|f^{*}-\hat{f}_{n}\|_{n}=O_{\mathbb{P}}(M_{1}),\|\hat{f}_{n}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}=O_{\mathbb{P}}(M_{2}).$

Furthermore, if $\tilde{f}_{n}$ be the solution to the optimization problem

$\displaystyle\min_{f\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|f^{*}-f\|_{n}^{2}+\lambda_{n}\|f\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)},$
(74)

then

$\displaystyle\|f^{*}-\tilde{f}_{n}\|_{n}=O_{\mathbb{P}}((T+n^{-1/2}T^{1/2})^{1/2}),\|\tilde{f}_{n}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}=O_{\mathbb{P}}((\lambda_{n}^{-1}(T+n^{-1/2}T^{1/2}))^{1/2}).$

######  Lemma E.4 (Interpolation inequality for Gaussian RKHS).

Let $g\in\mathcal{H}_{\sigma}(\mathbb{R}^{D})$. For any $1>r>0$, we have

$\displaystyle\|g\|_{L_{\infty}(\mathbb{R}^{D})}\leq
C_{4}r^{-\frac{D}{4}}\sigma^{\frac{D(r-1)}{2}}\|g\|_{L_{2}(\mathbb{R}^{D})}^{1-r}\|g\|_{\mathcal{H}_{\sigma}(\mathbb{R}^{D})}^{r},$

where $C_{4}$ is a constant not related to $r,\sigma$ and $g$.

###  E.1 Without weight decay

We first decompose the error as

$\displaystyle\|f_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}\leq$
$\displaystyle\|f_{t}-g_{t}\|_{L_{2}(P_{\mathbf{X}})}+\|g_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})},$
(75)

where $g_{t}$ is as in (31).

By Lemma B.2, the first term $\|f_{t}-g_{t}\|_{L_{2}(P_{\mathbf{X}})}$ in (75)
can be bounded by

$\displaystyle\|f_{t}-g_{t}\|_{L_{2}(P_{\mathbf{X}})}\leq
C_{5}\|f_{t}-g_{t}\|_{L_{\infty}(\Omega)}=O_{\mathbb{P}}\left(\frac{n^{2}\sqrt{\log
N/N}}{\eta_{n}(\tilde{\mathbf{K}})^{2}}\right),$

as long as

$\displaystyle\frac{1}{2}\eta_{n}(\tilde{\mathbf{K}})\geq n\sqrt{\frac{\log
N}{N}}.$ (76)

Choose

$\displaystyle N_{0}=\frac{4n^{2}}{\eta_{n}(\tilde{\mathbf{K}})^{2}}.$ (77)

Then it holds that when $N\geq N_{0}$,

$\displaystyle\|f_{t}-g_{t}\|_{L_{2}(P_{\mathbf{X}})}=O_{\mathbb{P}}\left(n^{-1/2}\right).$
(78)

It remains to consider $\|g_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}$. We consider
the empirical version of $\|g_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}$, and let

$\displaystyle
J_{2}=\|g_{t}-f^{*}\|_{n}^{2}=\frac{1}{n}\|g_{t}(\mathbf{X})-f^{*}(\mathbf{X})\|_{2}^{2}.$
(79)

Let $(\beta t)^{-1}=n\lambda_{n}$. Consider the kernel ridge regression

$\displaystyle\tilde{g}=\mathop{\mathrm{argmin}}_{f\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|f-\bm{y}\|_{n}^{2}+\lambda_{n}\|f\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}.$

By the representer theorem,
$\tilde{g}(\bm{x})=\tilde{\mathbf{k}}(\bm{x})^{T}(\tilde{\mathbf{K}}+n\lambda_{n}\mathbf{I})^{-1}\bm{y}$
for all $\bm{x}\in\Omega$. Then it can be seen that

$\displaystyle\tilde{g}(\mathbf{X})-f^{*}(\mathbf{X})=n\lambda_{n}(\tilde{\mathbf{K}}+n\lambda_{n}\mathbf{I})^{-1}f^{*}(\mathbf{X})+\tilde{\mathbf{K}}(\tilde{\mathbf{K}}+n\lambda_{n}\mathbf{I})^{-1}\bm{\epsilon}=\bm{q}_{1}+\bm{q}_{2}$

Following the arguments in Appendix D.1, the term $J_{2}$ can be bounded by

$\displaystyle
J_{2}\leq\frac{2}{n}\left(2C_{6}\|\bm{q}_{1}\|_{2}^{2}+8\|\bm{q}_{2}\|_{2}^{2}\right),$
(80)

and

$\displaystyle\frac{1}{n}\|\bm{q}_{1}\|_{2}^{2}=$
$\displaystyle\|\tilde{f}_{n}-f^{*}\|_{n}^{2}$

and

$\displaystyle\frac{1}{n}\|\bm{q}_{2}\|_{2}^{2}=$
$\displaystyle\|\hat{f}_{0,n}-f_{0}\|_{n}^{2},$

where $f_{0}(\bm{x})=0$ for all $\bm{x}\in\Omega$, $\tilde{f}_{n}$ is as in
(74), and $\hat{f}_{0,n}$ is the solution to the optimization problem

$\displaystyle\min_{g\in\mathcal{H}_{\tilde{K}_{S}}(\Omega)}\|\bm{\epsilon}-g\|_{n}^{2}+\lambda_{n}\|g\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}.$

By setting $\beta t\asymp n^{\frac{2m_{0}-d}{2m_{f}+d}}$ (which implies
$\lambda_{n}\asymp n^{-\frac{2m_{0}+2m_{f}}{2m_{f}+d}}$), $\sigma_{n}\asymp
n^{-\frac{1}{2m_{f}+d}}$, Lemma E.2 implies that $T\asymp
n^{-\frac{2m_{f}}{2m_{f}+d}}$, which, together with Lemma E.3, implies

$\displaystyle\frac{1}{n}\|\bm{q}_{1}\|_{2}^{2}=$
$\displaystyle\|\tilde{f}_{n}-f^{*}\|_{n}^{2}=O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}}(\log
n)^{D+1}\right),$ $\displaystyle\frac{1}{n}\|\bm{q}_{2}\|_{2}^{2}=$
$\displaystyle\|\hat{f}_{0,n}-f_{0}\|_{n}^{2}=O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}}(\log
n)^{D+1}\right).$ (81)

By (E.1) and (80), we obtain

$\displaystyle J_{2}=O_{\mathbb{P}}\left(n^{-\frac{2m_{f}}{2m_{f}+d}}(\log
n)^{D+1}\right).$ (82)

Next, we consider bounding $\|g_{t}-f^{*}\|_{L_{2}(P_{\mathbf{X}})}$. Similar
to the proof in Appendix D.1, it suffices to consider bounding the difference
between $\|g_{t}-f_{n}^{*}\|_{L_{2}(P_{\mathbf{X}})}$ and
$\|g_{t}-f_{n}^{*}\|_{n}$. Lemma E.1 implies that

$\displaystyle\|\tilde{g}\|_{\mathcal{H}_{\sigma_{n}/\sqrt{2}}(\Omega)}^{2}\leq
C_{7}\sigma_{n}^{-D}\|\tilde{g}\|_{\mathcal{H}_{\tilde{K}_{S}}(\Omega)}^{2}=O_{\mathbb{P}}\left(n^{\frac{2m_{0}+D}{2m_{f}+d}}(\log
n)^{D+1}\right).$ (83)

Consider function class