#  Robust Nonparametric Regression with
Deep Neural Networks

Guohao Shenlabel=e1 [<EMAIL_ADDRESS>Yuling Jiaolabel=e2 [
<EMAIL_ADDRESS>Yuanyuan Linlabel=e3 [
<EMAIL_ADDRESS>Jian Huanglabel=e4 [<EMAIL_ADDRESS>Department of Statistics, The Chinese University of Hong Kong, Hong Kong,
China  School of Mathematics and Statistics, Wuhan University, China
Department of Statistics, The Chinese University of Hong Kong, Hong Kong,
China  Department of Statistics and Actuarial Science, University of Iowa,
Iowa, USA

###### Abstract

In this paper, we study the properties of robust nonparametric estimation
using deep neural networks for regression models with heavy tailed error
distributions. We establish the non-asymptotic error bounds for a class of
robust nonparametric regression estimators using deep neural networks with
ReLU activation under suitable smoothness conditions on the regression
function and mild conditions on the error term. In particular, we only assume
that the error distribution has a finite $p$-th moment with $p$ greater than
one. We also show that the deep robust regression estimators are able to
circumvent the curse of dimensionality when the distribution of the predictor
is supported on an approximate lower-dimensional set. An important feature of
our error bound is that, for ReLU neural networks with network width and
network size (number of parameters) no more than the order of the square of
the dimensionality $d$ of the predictor, our excess risk bounds depend sub-
linearly on the $d$. Our assumption relaxes the exact manifold support
assumption, which could be restrictive and unrealistic in practice. We also
relax several crucial assumptions on the data distribution, the target
regression function and the neural networks required in the recent literature.
Our simulation studies demonstrate that the robust methods can significantly
outperform the least squares method when the errors have heavy-tailed
distributions and illustrate that the choice of loss function is important in
the context of deep nonparametric regression.

Keywords: Circumventing the curse of dimensionality; deep neural networks;
heavy-tailed error; non-asymptotic error bound; low-dimensional manifolds.

111Guohao Shen and Yuling Jiao contributed equally to this work.

∗ and

##  1 Introduction

Consider a nonparametric regression model

$Y=f_{0}(X)+\eta,$ (1.1)

where $Y\in\mathbb{R}$ is a response, $X\in\mathcal{X}\subseteq\mathbb{R}^{d}$
is a $d$-dimensional vector of predictors, $f_{0}:\mathcal{X}\to\mathbb{R}$ is
an unknown regression function, $\eta$ is an unobservable error independent of
$X$, whose $p$-th moment is assumed finite. A basic problem in statistics and
machine learning is to estimate the unknown target regression function $f_{0}$
based on observations that are independent and identically distributed
(i.i.d.) copies of $(X,Y)$, $(X_{i},Y_{i}),i=1,\ldots,n,$ where $n$ is the
sample size. In this paper, we establish the non-asymptotic error bounds for a
class of robust nonparametric regression estimators using deep neural networks
under mild conditions on the error term, allowing heavy-tailed error
distributions. Through simulation studies, we demonstrate that it is important
to use a robust method for heavy-tailed or contaminated data in nonparametric
regression using deep neural networks. We also show that the deep robust
regression estimators are able to circumvent the curse of dimensionality when
the predictor is supported on an approximate lower-dimensional set.

To illustrate the importance of using robust methods for nonparametric
regression using deep neural networks in the presence of heavy-tailed errors,
we look at the fitting of the regression functions “Blocks”, “Bumps”,
“Heavisine” and “Dopller” (Donoho and Johnstone, 1994), when error follows a
contaminated normal distribution (Huber, 1964), represented by a mixture of
two normals, $\eta\sim 0.8\mathcal{N}(0,1)+0.2\mathcal{N}(0,10^{4})$. The
functional form of the models are given in Section 6. Fig 1 shows the fitting
results using LS, LAD, Huber, Cauchy and Tukey loss functions. We see that the
least squares method breaks down in the presence of contaminated normal
errors, while the robust methods using LAD, Huber and Cauchy loss functions
yield reasonable fitting results. In particular, the Huber method works best
among the methods considered in this example.

Figure 1: The fitting of four univariate models. The data points are displayed
as grey dots, the target functions are depicted as black curves and the
colored curves represent the estimated curves under different loss functions.
Red line: LS, Orange line: LAD, Green line: Huber, Purple line: Cauchy; Blue
line: Tukey. The data are generated based on the univariate regression models
in Section 6 with contaminated normal error.

There is a vast literature on nonparametric regression, many important methods
have been developed. Examples include kernel regression (Nadaraya, 1964;
Watson, 1964), tree based regression (Breiman et al., 1984), local polynomial
regression (Cleveland, 1979; Fan, 1992, 1993), various spline-based regression
methods (Stone, 1982; Wahba, 1990; Fridedman, 1991; Green and Silverman,
1993), regression in reproducing kernel Hilbert spaces (Wahba, 1990; Scholkopf
and Smola, 2018), among others. Most of the existing methods have been
developed using the least squares loss as the objective function over
different function spaces. See also the books (Hastie, Tibshirani and
Friedman, 2001; Györfi et al., 2002; Tsybakov, 2008) for detailed discussions
on the general methodologies and theories on nonparametric regression and
related topics.

Recently, several inspiring works studied the properties of nonparametric
regression using deep neural networks (Bauer and Kohler, 2019; Schmidt-Hieber,
2019; Chen et al., 2019; Schmidt-Hieber, 2020; Farrell, Liang and Misra,
2021). In these works, convergence rates of the empirical risk minimizers
based on the least squares objective function or a convex loss are established
under suitable smoothness conditions on the underlying regression function and
certain regularity conditions on the model. In particular, the response
variable or the error term in the regression model is assumed to be bounded
(Györfi et al., 2002; Farrell, Liang and Misra, 2021), or sub-Gaussian (Bauer
and Kohler, 2019; Chen et al., 2019; Schmidt-Hieber, 2019, 2020).

The aforementioned results are not applicable to the problems with heavy-
tailed errors that commonly arise in practice, due to the sub-Gaussian
assumption on the error distribution or the boundedness assumption on the
response variable. Also, it is well known that nonparametric estimators based
on least squares loss function are sensitive to outliers. To address the non-
robustness problem of the methods based on least squares and to deal with
heavy- tailed errors, many robust regression methods had been developed in the
context of linear regression models. Important examples include the least
absolute deviation (LAD) (Bassett and Koenker, 1978), Huber regression (Huber,
1973), Tukey’s biweight regression (Beaton and Tukey, 1974), among others.
Indeed, there is a rich literature on robust statistics going back many
decades (Anscombe, 1960; Tukey, 1960; Huber, 1964; Tukey, 1975). For
systematic discussions on the basic problems in robust statistics, we refer to
Huber (2004) and Hampel et al. (2011). Recently, there is a resurgent interest
in developing polynomial time robust estimators achieving near optimal
recovery guarantees for many parametric estimation problems including mean
estimation, linear regression and generalized linear models, where an analyst
is given access to samples, in which a fraction of the samples may have been
adversarially corrupted, see for example, Cherapanamjeri et al. (2020),
Jambulapati et al. (2021) and the references therein.

As every aspects of scientific research and the society are becoming
digitized, statistics and machine learning methods are applied increasingly to
problems with more complex and noisier data. It is important that data
modeling methodologies and algorithms are robust to potentially nonnormal
noises. In particular, it has been recognized that robustness is an important
and desirable property in applications using deep learning methods
(Belagiannis et al., 2015; Wang et al., 2016; Jiang et al., 2018; Barron,
2019). However, there has been only limited studies on the statistical
guarantees of robust nonparametric regression using deep neural networks in
the literature. Lederer (2020) studied non-asymptotic error bounds for the
empirical risk minimizer with robust loss functions under the assumption that
the response has a finite second moment. Padilla, Tansey and Chen (2020)
derived the convergence rates for quantile regression with ReLU networks. In
both works, the authors assumed that the true target function belongs to the
function class of neural networks without considering the approximation error.
However, it is unrealistic to assume that the underlying target regression
function would be a member of the function class used in estimation. More
recently, Farrell, Liang and Misra (2021) considered general Lipschitz loss
functions under possibly mis-specified case (i.e., the true target function
may not belong to the function class used for estimation), but they assume
that the response is bounded, which is restrictive and not compatible with
heavy-tailed error models. Also, the existing error bound results for the
nonparametric regression estimators using deep neural networks contain
prefactors that depend exponentially on the ambient dimension $d$ of the
predictor (Schmidt-Hieber, 2020; Farrell, Liang and Misra, 2021; Padilla,
Tansey and Chen, 2020). This adversely affects the quality of the error
bounds, especially in the high-dimensional settings when $d$ is large
(Ghorbani et al., 2020). In particular, such type of error bounds lead to a
sample complexity that depends exponentially on $d$. Such a sample size
requirement is difficult to be met even for a moderately large $d$.

In this paper, we study the properties of robust nonparametric regression
using deep neural networks. We establish non-asymptotic upper bounds of the
excess risk of the empirical risk minimizer for the nonparametric neural
regression using feedforward neural networks with Rectified Linear Unit (ReLU)
activation. Specifically, our main contributions are as follows:

  1. (i)

For a heavy-tailed response $Y$ with a finite $p$-th moment for some $p>1$, we
establish non-asymptotic bounds on the excess risk of robust nonparametric
regression estimators using deep neural networks under continuous Lipschitz
loss functions. Moreover, for ReLU neural networks with network width and
network size (number of parameters) no more than $O(d^{2})$, our excess risk
bounds depend sub-linearly on the dimensionality $d$, instead of exponentially
in terms of a factor $a^{d}$ (for some constant $a\geq 2$) as in the existing
results.

  2. (ii)

We explicitly describe how the error bounds depend on the neural network
parameters, including the width, the depth and the size of the network. The
notion of network relative efficiency between two types of neural networks,
defined as the ratio of the logarithms of the network sizes needed to achieve
the optimal rate of convergence, is used to evaluate the relative merits of
network structures.

  3. (iii)

We show that the curse of dimensionality can be alleviated if the distribution
of $X$ is assumed to be supported on an approximate low-dimensional manifold.
This assumption relaxes the exact manifold support assumption, which could be
restrictive and unrealistic in practice. Under such an approximate low-
dimensional manifold support assumption, we show that the rate of convergence
$n^{-(1-1/p)\alpha/(\alpha+d)}$ can be improved to
$n^{-(1-1/p)\alpha/(\alpha+d_{0})}$ for $d_{0}=O(d_{\mathcal{M}}\log(d))$,
where $d_{\mathcal{M}}$ is the intrinsic dimension of the low-dimensional
manifold and $\alpha\in(0,1]$ is the assumed order of the Hölder-continuity of
the target $f^{*}$.

  4. (iv)

We relax several crucial assumptions on the data distribution, the target
regression function and the neural networks required in the recent literature.
First, we do not assume that the response $Y$ is bounded, instead only its
$p$-th moment is finite. Second, we do not require the network to be sparse or
have uniformly bounded weights and biases. Third, we largely relax the
regularity conditions on the true target function $f^{*}$. In particular, we
only require $f^{*}$ to be uniformly continuous, as our results are derived in
terms of its modulus of continuity, which is general enough to cover nearly
all cases that $f^{*}$ belongs to different smooth classes of functions.

  5. (v)

We carry our numerical experiments to evaluate the performance of robust
nonparametric regression estimators using LAD, Huber, Cauchy and Tukey loss
functions and compare with the performance of the least squares estimator. The
results demonstrate that, with heavy-tailed errors, the robust methods with
LAD, Huber and Cauchy loss functions tend to significantly outperform the
least squares method.

The remainder of the paper is organized as follows. In Section 2 we introduce
the setup of the problem, present some popular robust loss functions, and
describe the class of ReLU activated feedforward neural networks used in
estimating the regression function. In Section 3 we present a basic inequality
for the excess risk in terms of the stochastic and approximation errors, and
describe our approach in analyzing the two types of errors. In Section 4 we
provide sufficient conditions under which the robust neural regression
estimator is shown to possess the basic consistency property. We also
establish non-asymptotic error bounds for the neural regression estimator
using deep feedforward neural networks. Moreover, we show that the neural
regression estimator can circumvent the curse of dimensionality if the data
distribution is supported on an approximate low-dimensional manifold. In
Section 5 we present the results on how the error bounds depend on the network
structures and consider a notion of network relative efficiency between two
types of neural networks, defined as the ratio of the logarithms of the
network sizes needed to achieve the optimal convergence rate. This can be used
as a quantitative measure for evaluating the relative merits of different
network structures. Extensive simulation studies are presented in section 6.
Concluding remarks are given in section 7.

##  2 Preliminaries

In this section, we introduce the robust nonparametric regressions, including
the commonly-used robust loss functions, and the definition of ReLU
feedforward neural networks.

###  2.1 Robust regression

To accommodate heavy-tailed errors, we consider a general robust paradigm for
estimating $f_{0}$ by using a robust loss function
$L(\cdot,\cdot):\mathbb{R}^{2}\to\mathbb{R}$. For example, our results apply
to the following robust loss functions.

  * •

Least absolute deviation (LAD) loss: $L(a,y)=|a-y|,\ (a,y)\in\mathbb{R}^{2}.$

  * •

Quantile loss: $L(a,y)=\rho_{\tau}(a-y),\ (a,y)\in\mathbb{R}^{2},$ where
$\rho_{\tau}(x)=\tau x$ if $x\geq 0$ and $\rho_{\tau}(x)=(\tau-1)x$ if $x<0$
for some $\tau\in(0,1)$.

  * •

Huber loss: $L(a,y)=h_{\zeta}(a-y),\ \ (a,y)\in\mathbb{R}^{2},\text{ for some
}\zeta>0,$ where where $h_{\zeta}(x)=x^{2}/2$ if $|x|\leq\zeta$ and
$h_{\zeta}(x)=\zeta|x|-\zeta^{2}/2$ otherwise.

  * •

Cauchy loss: $L(a,y)=\log\\{1+\kappa^{2}(a-y)^{2}\\},\
(a,y)\in\mathbb{R}^{2},$ for some $\kappa>0.$

  * •

Tukey’s biweight loss: $L(a,y)=T_{t}(a-y),\ (a,y)\in\mathbb{R}^{2},\text{ for
some }t>0,$ where $T_{t}(x)=t^{2}[1-\\{1-(x/t)^{2}\\}^{3}]/6$ if $|a|\leq t$
and $T_{t}(x)=t^{2}/6$ otherwise.

These popular robust loss functions are all continuous and Lipschitz in both
of its arguments. Throughout this paper, we consider a general and broad class
of loss functions that are continuous and Lipschitz in both of its arguments,
as stated in the following assumption.

######  Assumption 1.

The loss function $L(\cdot,\cdot):\mathbb{R}^{2}\to\mathbb{R}$ is continuous,
and $L(a,y)=0$ if $a=y$ for $(a,y)\in\mathbb{R}^{2}$. Besides, $L$ is
$\lambda_{L}$Lipschitz in both of its arguments, i.e.

$\displaystyle|L(a_{1},\cdot)-L(a_{2},\cdot)|$
$\displaystyle\leq\lambda_{L}|a_{1}-a_{2}|,$
$\displaystyle|L(\cdot,y_{1})-L(\cdot,y_{2})|$
$\displaystyle\leq\lambda_{L}|y_{1}-y_{2}|,$

for any $a_{1},a_{2},y_{1},y_{2}\in\mathbb{R}$.

The aforementioned robust loss functions all satisfy Assumption 1. To be
specific, for the LAD loss, $\lambda_{L}=1$; for the quantile loss,
$\lambda_{L}=\max\\{\tau,1-\tau\\}$; for the Huber loss, $\lambda_{L}=\zeta$;
for the Cauchy loss, $\lambda_{L}=\kappa$; for the Tukey’s biweight loss,
$\lambda_{L}=16t/25\sqrt{5}$. In particular, the LAD, the quantile check loss
and the Huber loss functions can be recast into a more simplified form, i.e.
$L(a,y)=\psi(a-y)$ for for some convex or Lipschitz function $\psi$ with
$\psi(0)=0$. A detailed summary is given in Table 1, and plots of the above
loss functions (except quantile loss function) and their derivatives in
comparison with least squares are displayed in Figure 2.

Table 1: Summary of different robust loss functions

| LAD | Quantile | Huber | Cauchy | Tukey
---|---|---|---|---|---
Hyper parameter | N/A | $\tau\in(0,1)$ | $\zeta\in(0,\infty)$ | $\kappa\in(0,\infty)$ | $t\in(0,\infty)$
$\lambda_{L}$ | 1 | $\max\\{\tau,1-\tau\\}$ | $\zeta$ | $\kappa$ | $16t/25\sqrt{5}$
Continuous | TRUE | TRUE | TRUE | TRUE | TRUE
Convex | TRUE | TRUE | TRUE | FALSE | FALSE
Differentiable | FALSE | FALSE | TRUE | TRUE | TRUE

Note: ”N/A” stands for not applicable.

Figure 2: Different loss functions (solid lines in the top row) and their
corresponding derivatives (dashed lines in the bottom row). For depicted Huber
loss function, $\zeta=1$; Cauchy loss function, $\kappa=1$; Tukey loss
function, $t=4.685$.

For any (random) function $f$, let $Z\equiv(X,Y)$ be a random vector
independent of $f$. For any loss function $L$ satisfying Assumption 1, the $L$
risk is defined by $\mathcal{R}(f)=\mathbb{E}_{Z}L(f(X),Y)$. At the population
level, an estimation procedure is to find a measurable function
$f^{*}:\mathbb{R}^{d}\to\mathbb{R}$ satisfying

$f^{*}:=\arg\min_{f}\mathcal{R}(f)=\arg\min_{f}\mathbb{E}_{Z}L(f(X),Y).$ (2.1)

Under proper assumptions on $L$, $X$ and $\eta$, the true regression function
$f_{0}$ in model (1.1) is the optimal solution $f^{*}$ on $\mathcal{X}$. For
instance, when $L(f(X),Y)=|f(X)-Y|$ and under the assumption that the
conditional median of $\eta$ given $X$ is $0$, $f^{*}$ in (2.1) is precisely
the true regression function $f_{0}$ on $\mathcal{X}$. For the rest of the
paper, under a general robust nonparametric estimation framework without
specific assumptions on $L,X$ or $\eta$ to ensure $f^{*}=f_{0}$ on
$\mathcal{X}$, we would focus on estimating $f^{*}$ instead of $f_{0}$.

In real applications, the distribution of $(X,Y)$ is typically unknown and
only a random sample $S\equiv\\{(X_{i},Y_{i})\\}_{i=1}^{n}$ with size $n$ is
available. Define the empirical risk of $f$ on the sample $S$ as

$\mathcal{R}_{n}(f)=\frac{1}{n}\sum_{i=1}^{n}L(f(X_{i}),Y_{i}).$ (2.2)

Based on the data, a primary goal is to construct estimators of $f^{*}$ within
a certain class of functions $\mathcal{F}_{n}$ by minimizing the empirical
risk. Thus, the empirical risk minimizer (ERM) is defined by

$\hat{f}_{n}\in\arg\min_{f\in\mathcal{F}_{n}}\mathcal{R}_{n}(f).$ (2.3)

In this paper, we set $\mathcal{F}_{n}$ to be a function class consisting of
feedforward neural networks. For any estimator $\hat{f}_{n}$, we evaluate its
quality via its excess risk, which is defined as the difference between the
$L$ risks of $\hat{f}_{n}$ and $f^{*}$:

$\displaystyle\mathcal{R}(\hat{f}_{n})-\mathcal{R}(f^{*})=\mathbb{E}_{Z}L(\hat{f}_{n}(X),Y)-\mathbb{E}_{Z}L(f^{*}(X),Y).$

A good estimator $\hat{f}_{n}$ shall have a small excess risk
$\mathcal{R}(\hat{f}_{n})-\mathcal{R}(f^{*})$. Thereafter, we shall focus on
deriving the non-asymptotic upper bounds of the excess risk
$\mathcal{R}(\hat{f}_{n})-\mathcal{R}(f^{*})$. Note that
$\mathcal{R}(\hat{f}_{n})$ is random as $\hat{f}_{n}$ is produced based on the
sample $S$. Thus, we will derive bounds for the expected excess risk
$\mathbb{E}_{S}\big{\\{}\mathcal{R}(\hat{f}_{n})-\mathcal{R}(f^{*})\big{\\}}$
as well.

###  2.2 ReLU Feedforward neural networks

In recent years, deep neural network modeling has achieved impressive
successes in many applications. Also, neural network functions have proven to
be an effective tool to approximate high-dimensional functions. We consider
regression function estimators based on the feedforward neural networks with
rectified linear unit (ReLU) activation function. Specifically, we set the
function class $\mathcal{F}_{n}$ to be
$\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$, a
class of feedforward neural networks $f_{\phi}:\mathbb{R}^{d}\to\mathbb{R}$
with parameter $\phi$, depth $\mathcal{D}$, width $\mathcal{W}$, size
$\mathcal{S}$, number of neurons $\mathcal{U}$ and $f_{\phi}$ satisfying
$\|f_{\phi}\|_{\infty}\leq\mathcal{B}$ for some $0<B<\infty$, where
$\|f\|_{\infty}$ is the sup-norm of a function $f$. Note that the network
parameters may depend on the sample size $n$, but the dependence is omitted in
the notation for simplicity. A brief description of the feedforward neural
networks are given below.

We begin with the multi-layer perceptron (MLP), an important subclass of
feedforward neural networks. The architecture of a MLP can be expressed as a
composition of a series of functions

$f_{\phi}(x)=\mathcal{L}_{\mathcal{D}}\circ\sigma\circ\mathcal{L}_{\mathcal{D}-1}\circ\sigma\circ\cdots\circ\sigma\circ\mathcal{L}_{1}\circ\sigma\circ\mathcal{L}_{0}(x),\
x\in\mathbb{R}^{d},$

where $\sigma(x)=\max(0,x)$ is the rectified linear unit (ReLU) activation
function (defined for each component of $x$ if $x$ is a vector) and

$\mathcal{L}_{i}(x)=W_{i}x+b_{i},\quad i=0,1,\ldots,\mathcal{D},$

where $W_{i}\in\mathbb{R}^{d_{i+1}\times d_{i}}$ is a weight matrix, $d_{i}$
is the width (the number of neurons or computational units) of the $i$-th
layer, and $b_{i}\in\mathbb{R}^{d_{i+1}}$ is the bias vector in the $i$-th
linear transformation $\mathcal{L}_{i}$. The input data consisting of
predictor values $X$ is the first layer and the output is the last layer. Such
a network $f_{\phi}$ has $\mathcal{D}$ hidden layers and $(\mathcal{D}+1)$
layers in total. We use a $(\mathcal{D}+1)$-vector
$(d_{0},d_{1},\ldots,d_{\mathcal{D}})^{\top}$ to describe the width of each
layer; particularly, $d_{0}=d$ is the dimension of the input $X$ and
$d_{\mathcal{D}}=1$ is the dimension of the response $Y$ in model (1.1). The
width $\mathcal{W}$ is defined as the maximum width of hidden layers, i.e.,
$\mathcal{W}=\max\\{d_{1},...,d_{\mathcal{D}}\\}$; the size $\mathcal{S}$ is
defined as the total number of parameters in the network $f_{\phi}$, i.e
$\mathcal{S}=\sum_{i=0}^{\mathcal{D}}\\{d_{i+1}\times(d_{i}+1)\\}$; the number
of neurons $\mathcal{U}$ is defined as the number of computational units in
hidden layers, i.e., $\mathcal{U}=\sum_{i=1}^{\mathcal{D}}d_{i}$. Note that
the neurons in consecutive layers of a MLP are connected to each other via
linear transformation matrices $W_{i}\in\mathbb{R}^{d_{i+1}\times d_{i}}$,
$i=0,1,\ldots,\mathcal{D}$. In other words, an MLP is fully connected between
consecutive layers and has no other connections. For an MLP
$\mathcal{F}_{\mathcal{D},\mathcal{U},\mathcal{W},\mathcal{S},\mathcal{B}}$,
its parameters satisfy the simple relationship

$\max\\{\mathcal{W},\mathcal{D}\\}\leq\mathcal{S}\leq\mathcal{W}(d+1)+(\mathcal{W}^{2}+\mathcal{W})(\mathcal{D}-1)+\mathcal{W}+1=O(\mathcal{W}^{2}\mathcal{D}).$

Different from multilayer perceptron, a general feedforward neural network may
not be fully connected. For such a network, each neuron in layer $i$ may be
connected to only a small subset of neurons in layer $i+1$. The total number
of parameters
$\mathcal{S}\leq\sum_{i=0}^{\mathcal{D}}\\{d_{i+1}\times(d_{i}+1)\\}$ is
reduced and the computational cost required to evaluate the network will also
be reduced.

For notational simplicity, we write
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
for short to denote the class of feedforward neural networks with parameters
$\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S}$ and $\mathcal{B}$. In the
following, we will present our main results: nonasymptotic upper bounds on the
excess risk for general feedforward neural networks with piecewise linear
activation function.

##  3 Basic error analysis

For general robust estimation, we start by introducing a generic upper bound
on the excess risk of the empirical risk minimizer $\hat{f}_{n}$ defined in
(2.3) with a general function class $\mathcal{F}_{n}$. To this end, for a
general loss function $L$ and any estimator $f$ belonging to
$\mathcal{F}_{n}$, the excess risk of $f$ can be written as

$\mathcal{R}(f)-\mathcal{R}(f^{*})=\left\\{\mathcal{R}(f)-\inf_{f\in\mathcal{F}_{n}}\mathcal{R}(f)\right\\}+\left\\{\inf_{f\in\mathcal{F}_{n}}\mathcal{R}(f)-\mathcal{R}(f^{*})\right\\}.$

The first term on the right hand side is the so-called stochastic error, which
measures the difference between the error of $f$ and the best one in the class
$\mathcal{F}_{n}$; the second term is the commonly-known approximation error,
which measures how well the function class $\mathcal{F}_{n}$ with respect to
the loss function $L$ can be used to approximate the function $f^{*}$.

It follows from Lemma 3.1 in Jiao et al. (2021) that, for any random sample
$S=\\{(X_{i},Y_{i})_{i=1}^{n}\\}$, the excess risk of the empirical risk
minimizer $\hat{f}_{n}$ satisfies

$\displaystyle\mathcal{R}(\hat{f}_{n})-\mathcal{R}(f^{*})\leq$ $\displaystyle
2\sup_{f\in\mathcal{F}_{n}}|\mathcal{R}(f)-\mathcal{R}_{n}(f)|+\inf_{f\in\mathcal{F}_{n}}\mathcal{R}({f})-\mathcal{R}(f^{*}).$
(3.1)

The upper bound of the excess risk of $\hat{f}_{n}$ in (3.1) does not depend
on $\hat{f}_{n}$ any more, but the function class $\mathcal{F}_{n}$, the loss
function $L$ and the random sample $S$. The first term in the upper bound
$2\sup_{f\in\mathcal{F}_{n}}|\mathcal{R}(f)-\mathcal{R}_{n}(f)|$ can be
properly handled by the classical empirical process theory (van der Vaart and
Wellner, 1996; Anthony and Bartlett, 1999; Bartlett et al., 2019), whose upper
bound is determined by the complexity of the function class $\mathcal{F}_{n}$.
The second term
$\inf_{f\in\mathcal{F}_{n}}\mathcal{R}({f})-\mathcal{R}(f^{*})$ is the
approximation error of the function class $\mathcal{F}_{n}$ for $f^{*}$ under
loss $L$. The neural networks have been shown to have strong power to
approximate high-dimensional functions in both theory and practice. Important
recent results on the approximation theory of neural networks can be found in
Yarotsky (2017, 2018); Shen, Yang and Zhang (2019); Shen, Yang and Zhang
(2020), among others. In order to apply those novel approximation theories to
our problem, we need to write
$\inf_{f\in\mathcal{F}_{n}}\mathcal{R}({f})-\mathcal{R}(f^{*})$ into a form of
$\inf_{f\in\mathcal{F}_{n}}\|f-f^{*}\|$ for some functional norm $\|\cdot\|$.

###  3.1 Stochastic error

We now proceed to analyze the stochastic error of ERM $\hat{f}_{\phi}$
constructed by the feedforward neural networks and provide an upper bound on
the expected excess risk. Here, the ERM $\hat{f}_{\phi}$ is defined by

$\hat{f}_{\phi}\in\arg\min_{f\in\mathcal{F}_{\phi}}\mathcal{R}_{n}(f),$ (3.2)

where
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
is the class of feedforward neural networks with parameters
$\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S}$ and $\mathcal{B}$. The
stochastic error of $\hat{f}_{\phi}$ can be upper bounded by the uniform
results (on the class $\mathcal{F}_{\phi}$) in empirical process theory. In
light of this, we next introduce the conceptions of metric entropy of a class
of functions which are needed in stochastic error bounds.

Denote the pseudo dimension for a class of functions $\mathcal{F}$:
$\mathcal{X}\to\mathbb{R}$ by $\text{Pdim}(\mathcal{F}),$, defined as the
largest integer $m$ for which there exists
$(x_{1},\ldots,x_{m},y_{1},\ldots,y_{m})\in\mathcal{X}^{m}\times\mathbb{R}^{m}$
such that for any $(b_{1},\ldots,b_{m})\in\\{0,1\\}^{m}$ there exists
$f\in\mathcal{F}$ such that $\forall i:f(x_{i})>y_{i}\iff b_{i}=1$ (Anthony
and Bartlett, 1999; Bartlett et al., 2019). Pseudo dimension is a natural
measure of the complexity of a function class, e.g., a class of real-valued
functions generated by neural networks. Especially, if $\mathcal{F}$ is the
class of functions generated by a neural network with a fixed architecture and
fixed activation functions, it follows from Theorem 14.1 in Anthony and
Bartlett (1999) that $\text{Pdim}(\mathcal{F})=\text{VCdim}(\mathcal{F})$,
where $\text{VCdim}(\mathcal{F})$ is the VC dimension of $\mathcal{F}$. In our
results, it is required that the sample size $n$ is larger than the pseudo
dimension of the class of neural networks under consideration.

For a given sequence $x=(x_{1},\ldots,x_{n})\in\mathcal{X}^{n},$ we define
$\mathcal{F}_{\phi}|_{x}=\\{(f(x_{1}),\ldots,f(x_{n}):f\in\mathcal{F}_{\phi}\\}$,
a subset of $\mathbb{R}^{n}$. For a positive number $\delta$, let
$\mathcal{N}(\delta,\|\cdot\|_{\infty},\mathcal{F}_{\phi}|_{x})$ be the
covering number of $\mathcal{F}_{\phi}|_{x}$ under the norm
$\|\cdot\|_{\infty}$ with radius $\delta$. And the uniform covering number
$\mathcal{N}_{n}(\delta,\|\cdot\|_{\infty},\mathcal{F}_{\phi})$ is defined by

$\mathcal{N}_{n}(\delta,\|\cdot\|_{\infty},\mathcal{F}_{\phi})=\max\\{\mathcal{N}(\delta,\|\cdot\|_{\infty},\mathcal{F}_{\phi}|_{x}):x\in\mathcal{X}\\},$
(3.3)

the maximum of the covering number
$\mathcal{N}(\delta,\|\cdot\|_{\infty},\mathcal{F}_{\phi}|_{x})$ over all
$x\in\mathcal{X}$.

We are now ready to present a lemma concerning the upper bound of the
stochastic error under a heavy-tailed assumption of the response $Y$.

######  Assumption 2.

The response variable $Y$ has finite $p$-th moment with $p>1$, i.e., there
exists a constant $M>0$ such that $\mathbb{E}|Y|^{p}\leq M<\infty$.

######  Lemma 3.1.

Consider the $d$-variate nonparametric regression model in (1.1) with an
unknown regression function $f_{0}$ and target function $f^{*}$ defined in
(2.1). Assume that Assumptions 1-2 hold, and
$\|f^{*}\|_{\infty}\leq\mathcal{B}$ for $\mathcal{B}\geq 1$. Let
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
be the class of feedforward neural networks with a continuous piecewise-linear
activation function of finite pieces and $\hat{f}_{\phi}$ be the ERM defined
in (3.2). Then, for $\text{Pdim}(\mathcal{F}_{\phi})\leq 2n$,

$\sup_{f\in\mathcal{F}_{\phi}}|\mathcal{R}(f)-\mathcal{R}_{n}(f)|\leq
c_{0}\frac{\lambda_{L}\mathcal{B}}{n^{1-1/p}}\log\mathcal{N}_{2n}(n^{-1},\|\cdot\|_{\infty},\mathcal{F}_{\phi}),$
(3.4)

where $c_{0}>0$ is a constant independent of
$n,d,\lambda_{L},\mathcal{B},\mathcal{S},\mathcal{W}$ and $\mathcal{D}$. In
addition,

$\mathbb{E}\big{\\{}\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\big{\\}}\leq
C_{0}\frac{\lambda_{L}\mathcal{B}\mathcal{S}\mathcal{D}\log(\mathcal{S})\log(n)}{n^{1-1/p}}+2\inf_{f\in\mathcal{F}_{\phi}}\big{\\{}\mathcal{R}(f)-\mathcal{R}(f^{*})\big{\\}},$
(3.5)

where $C_{0}>0$ is a constant independent of
$n,d,\lambda_{L},\mathcal{B},\mathcal{S},\mathcal{W}$ and $\mathcal{D}$.

Lemma 3.1 indicates that, the stochastic error is bounded by a term determined
by the metric entropy of $\mathcal{F}_{\phi}$ in (3.4). The metric entropy is
measured by the covering number of $\mathcal{F}_{\phi}$. To obtain the bound
for excess risk in (3.5), we next bound the covering number of
$\mathcal{F}_{\phi}$ by its pseudo dimension (VC dimension), and apply the
result in Bartlett et al. (2019) that ${\rm
Pdim}(\mathcal{F}_{\phi})=O(\mathcal{S}\mathcal{D}\log(\mathcal{S}))$, which
gives the upper bound for the expected excess risk as the sum of the
stochastic error and the approximation error of $\mathcal{F}_{\phi}$ to
$f^{*}$ in (3.5). A few important remarks are in order.

######  Remark 1.

If $Y$ is assumed to be sub-exponentially distributed, i.e., there exists a
constant $\sigma_{Y}>0$ such that $\mathbb{E}\exp(\sigma_{Y}|Y|)<\infty$, the
denominator $n^{1-1/p}$ in (3.4) and (3.5) can be improved to be $n$,
corresponding to the case that $p=+\infty$. Notably, the finite $p$-th moment
condition and sub-exponential tail condition have intrinsically different
impacts on the resulting error bounds in nonparametric regressions.

######  Remark 2.

The non-asymptotic stochastic error bounds in Lederer (2020) are derived under
finite second-moment assumption of both covariates and response. Our result
covers their result as a special case by letting $p=2$ in (3.4) in Lemma 3.1.
Assuming finite $p$-th moment ($p>2$) would give better rate. Furthermore, our
result is more general than that of Lederer (2020), since only the well-
specified case, i.e. $f^{*}\in\mathcal{F}_{n}$, is considered and no
approximation error is involved in Lederer (2020).

######  Remark 3.

Upper bounds for least-squares estimator with neural nonparametric regression
were studied by Györfi et al. (2002) and Farrell, Liang and Misra (2021) under
the assumption that the response $Y$ is bounded. Bauer and Kohler (2019) and
Schmidt-Hieber (2020) derived error bounds for the least squares nonparametric
regression estimator with deep neural networks under the assumption that the
response $Y$ is sub-exponential.

######  Remark 4.

Han and Wellner (2019) studied general nonparametric least squares regression
with heavy-tailed errors whose $p$-th moment ($p\geq 1$) is assumed finite and
independent of the covariates. Under the well-specified scenario that
$f^{*}\in\mathcal{F}_{n}$ they studied the convergence rate for nonparametric
least squares estimator. They showed that the rate of convergence with respect
to the sample size $n$ is $O(\max\\{n^{-1/(2+\gamma)},n^{-1/2+1/(2p)}\\})$
when the model satisfies a standard entropy condition with exponent
$\gamma\in(0,2)$, i.e
$\log\mathcal{N}(\delta,\|\cdot\|_{L^{2}(\nu)},\mathcal{F}_{n})\leq
C\delta^{-\gamma}$ for some constant $C>0$, where $\nu$ is the probability
measure of $X$. Their results shed light on the fact that, when the order of
moment $p\geq 1+2/\gamma$, the $L_{2}$ risk of the nonparametric least squares
estimator converges at the same rate as the case with Gaussian errors;
otherwise for $p<1+2/\gamma$, the convergence rate of the $L_{2}$ risk of
least squares estimator can be strictly slower than other robust estimators.
On the other hand, the entropy condition with exponent $\gamma\in(0,2)$ on the
function class may not be satisfied in some scenarios, especially when the
dimensionality $d$ is high. For instance, for the Hölder class
$\mathcal{F}=C^{\alpha}_{1}(\mathcal{X})$ ($f:\mathcal{X}\to\mathbb{R}$ with
$\|f\|_{\alpha}\leq 1$ for all $f\in\mathcal{F}$), its entropy satisfies
$\log\mathcal{N}(\delta,\|\cdot\|_{L^{2}(\nu)},\mathcal{F})\leq
K_{\alpha,d,\mathcal{X}}\delta^{-d/\alpha}$ for some constant
$K_{\alpha,d,\mathcal{X}}$; for Lipschitz functions with Lipschitz constant
$L$ on $[a,b]^{d}$, its entropy satisfies
$\log\mathcal{N}(\delta,\|\cdot\|_{L^{2}(\nu)},\mathcal{F})\leq
C\\{(b-a)L\\}^{-d}\delta^{-d}$. Besides, for VC-class with envelope function
$F$ and $r\geq 1$, and any probability measure $Q$, its entropy satisfies
$\mathcal{N}(\delta\|F\|_{L^{r}(Q)},\|\cdot\|_{L^{r}(Q)},\mathcal{F})\leq
CV(\mathcal{F})(16e)^{V(\mathcal{F})}\delta^{-r\\{V(\mathcal{F})-1\\}}$ for
any $0<\delta<1$. In some situations above, the entropy condition is not
satisfied.

###  3.2 Approximation error

The approximation error depends on the target function $f^{*}$ and the
function class
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$.
We assume that the distribution of the predictor $X$ is supported on a bounded
set, and for simplicity, we assume this bounded set to be $[0,1]^{d}$.
Thereafter, we consider the target functions $f^{*}$ defined on
$\mathcal{X}=[0,1]^{d}$.

To derive an upper bound of the approximation error
$\inf_{f\in\mathcal{F}_{\phi}}\\{R(f)-R^{\phi}(f^{*})\\}$ based on existing
approximation theories, we need to fill the gap between
$\inf_{f\in\mathcal{F}_{\phi}}\\{R(f)-R^{\phi}(f^{*})\\}$ and
$\inf_{f\in\mathcal{F}_{\phi}}\|f-f^{*}\|$. Indeed, this can be done under
Assumption 1 as follows.

######  Lemma 3.2.

Assume that Assumption 1 holds. Let $f^{*}$ be the target function defined in
(2.1) and $\mathcal{R}(f^{*})$ be its risk, then we have

$\inf_{f\in\mathcal{F}_{\phi}}\\{R(f)-R(f^{*})\\}\leq\lambda_{L}\inf_{f\in\mathcal{F}_{\phi}}\mathbb{E}|f(X)-f^{*}(X)|=:\lambda_{L}\inf_{f\in\mathcal{F}_{\phi}}\|f-f^{*}\|_{L^{1}(\nu)},$

where $\nu$ denotes the marginal probability measure of $X$ and
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
denotes the class of feedforward neural networks with parameters
$\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S}$ and $\mathcal{B}$.

As a consequence of Lemma 3.2, we only need to give upper bounds on the
approximation error $\inf_{f\in\mathcal{F}_{\phi}}\|f-f^{*}\|_{L^{1}(\nu)}$,
so as to bound the excess risk of the ERM $\hat{f}_{\phi}$ defined in (3.2).

Next, the approximation error bound depends on $f^{*}$ through its modulus of
continuity and is related to the function class
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
via its parameters. Recall that the modulus of continuity $\omega_{f}$ of a
function $f:[0,1]^{d}\to\mathbb{R}$ is defined as

$\omega_{f}(r):=\sup\\{|f(x)-f(y)|:x,y\in[0,1]^{d},\|x-y\|_{2}\leq r\\},{\rm
for\ any\ }r\geq 0.$ (3.6)

A uniformly continuous function $f$ has the property that $\lim_{r\to
0}\omega_{f}(r)=\omega_{f}(0)=0$. Based on the modulus of continuity, one can
define different equicontinuous families of functions. For instance, the
modulus $\omega_{f}(r)=\theta r$ describes the $\theta$-Lipschitz continuity;
the modulus $\omega_{f}(r)=\theta r^{\alpha}$ with $\theta,\alpha>0$ defines
the Hölder continuity.

In this paper, we only assume that $f^{*}$ is uniformly continuous (in
Assumption 3 given below), which is a mild assumption on the continuity of the
unknown target function $f^{*}$, as the existing works generally assume
stronger smoothness assumptions on $f^{*}$.

######  Assumption 3.

The target function $f^{*}:\mathcal{X}\to\mathbb{R}$ is uniformly continuous.

Assumption 3 implies that $\lim_{r\to
0}\omega_{f^{*}}(r)=\omega_{f^{*}}(0)=0$. It provides a necessary condition
for us to derive an explicit upper bounds of the approximation error in terms
of its modulus of continuity $\omega_{f^{*}}$, and give sufficient conditions
for the consistency of some estimators.

The DNN approximation theories are actively studied in recent years. We
provide a brief review on the most relevant literature in Appendix A for
reader’s convenience, see also Jiao et al. (2021) for some additional
references. We will use the recent approximation results obtained by Shen,
Yang and Zhang (2020), who established an explicit error bounds for the
approximation using deep neural networks. Their results clearly describe how
the error bounds depend on the network parameters and are non-asymptotic. Let
$\mathbb{N}^{+}$ be be the set of positive integers. For any
$N,M\in\mathbb{N}^{+}$ and $p\in[1,\infty]$, Shen, Yang and Zhang (2020)
showed that

$\inf_{f\in\mathcal{F}_{\phi}}\|f-f^{*}\|_{L^{p}([0,1]^{d})}\leq
19\sqrt{d}\,\omega_{f^{*}}(N^{-2/d}M^{-2/d})$ (3.7)

for a Hölder continuous function $f^{*}$ with the modulus of continuity
$\omega_{f^{*}}(\cdot)$ defined in (3.6) and the class of functions
$\mathcal{F}_{\phi}$ consisting of ReLU activated feedforward networks with

$\text{width }\mathcal{W}=C_{1}\max\\{d\lfloor N^{1/d}\rfloor,N+1\\}\ \text{
and depth }\ \mathcal{D}=12M+C_{2},$

where $C_{1}=12,C_{2}=14$ if $p\in[1,\infty)$, and $C_{1}=3^{d+3},C_{2}=14+2d$
if $p=\infty$. Such an approximation rate in terms of the width and depth is
more informative and more useful than those characterized by just the size
$\mathcal{S}$, as the upper bound represented by width and depth can imply the
one in terms of the size $\mathcal{S}$. Moreover, the approximation error
bound is explicit in the sense that it does not involve unclearly defined
constant, in contrast to other existing approximation error bounds that
involve an unknown prefactor or require the network width $\mathcal{W}$ and
depth $\mathcal{D}$ greater than some unknown constants.

Note that the approximation error bound in Lemma 3.2 only requires the
Lipschitz property of the loss function; if the curvature information of the
risk around $f^{*}$ can be used, the approximation error bound can be further
improved.

######  Assumption 4 (Local quadratic bound of the excess risk).

There exist some constants $\Lambda_{L,f^{*}}>0$ and $\varepsilon_{L,f^{*}}>0$
which may depend on the loss function $L(\cdot)$, the risk minimizer $f^{*}$
and the distribution of $(X,Y)$ such that

$\mathcal{R}(f)-\mathcal{R}(f^{*})\leq\Lambda_{L,f^{*}}\|f-f_{0}\|^{2}_{L^{2}(\nu)},$

for any $f$ satisfying
$\|f-f_{0}\|_{L^{\infty}(\nu)}\leq\varepsilon_{L,f^{*}}$, where $\nu$ is the
probability measure of $X$.

######  Lemma 3.3.

Assume that Assumption 1 and 4 hold. Let $f^{*}$ be the target function
defined in (2.1) and $\mathcal{R}(f^{*})$ be its risk, then we have

$\inf_{f\in\mathcal{F}_{\phi}}\\{R(f)-R(f^{*})\\}\leq\lambda_{L,f^{*}}\inf_{f\in\mathcal{F}_{\phi}}\|f-f^{*}\|^{2}_{L^{2}(\nu)},$

where
$\lambda_{L,f^{*}}\geq\max\\{\Lambda_{L,f^{*}},\lambda_{L}/\varepsilon_{L,f^{*}}\\}>0$
is a constant, $\nu$ is the probability measure of $X$ and
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
denotes the class of feedforward neural networks with parameters
$\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S}$ and $\mathcal{B}$.

##  4 Non-asymptotic error bounds

Lemma 3.1 gives the theoretical basis for establishing the consistency and
non-asymptotic error bounds. To guarantee consistency, the two terms on the
right hand side of (3.5) should go to zero as $n\to\infty$. Regarding the non-
asymptotic error bound, the exact rate of convergence will be determined by a
trade-off between the stochastic error and the approximation error.

###  4.1 Consistency and non-asymptotic error bounds

In the following, we present the consistency and non-asymptotic error bound of
the general robust nonparametric regression estimator using neural networks.

######  Theorem 4.1 (Consistency).

For model (1.1), suppose that Assumptions 1-2 hold. Assume that the target
function $f^{*}$ is continuous function on $[0,1]^{d}$,
$\|f^{*}\|_{\infty}\leq\mathcal{B}$ for some $\mathcal{B}\geq 1$, and the
function class of feedforward neural networks
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
with continuous piecewise-linear activation function of finite pieces
satisfies

$\mathcal{S}\to\infty\quad{\rm
and}\quad\mathcal{B}\frac{\log(n)}{n^{1-1/p}}\mathcal{S}\mathcal{D}\log(\mathcal{S})\to
0,$

as $n\to\infty$. Then, the excess risk of the empirical risk minimizer
$\hat{f}_{\phi}$ defined in (3.2) is consistent in the sense that

$\mathbb{E}\\{\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\\}\to 0\ \ {\rm
as\ }n\to\infty.$

Theorem 4.1 is a direct result from Lemma 3.1 and Theorem 1 in Yarotsky
(2018). Those conditions in Theorem 4.1 are sufficient conditions to ensure
the consistency of the deep neural regression, and they are relatively mild in
terms of the assumptions on the target function $f^{*}$ and the distribution
of $Y$.

######  Theorem 4.2 (Non-asymptotic excess risk bound).

Suppose that Assumptions 1-3 hold, $\nu$ is absolutely continuous with respect
to the Lebesgue measure, and $\|f^{*}\|_{\infty}\leq\mathcal{B}$ for some
$\mathcal{B}\geq 1$. Then, for any $N,M\in\mathbb{N}^{+}$, the function class
of ReLU multi-layer perceptrons
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
with depth $\mathcal{D}=12M+14$ and width $\mathcal{W}=\max\\{4d\lfloor
N^{1/d}\rfloor+3d,12N+8\\}$, for $2n\geq\text{Pdim}(\mathcal{F}_{\phi})$, the
expected excess risk of the ERM $\hat{f}_{\phi}$ defined in (3.2) satisfies

$\mathbb{E}\\{\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\\}\leq
C\lambda_{L}\mathcal{B}\frac{\mathcal{S}\mathcal{D}\log(\mathcal{S})\log(n)}{n^{1-1/p}}+18\lambda_{L}\sqrt{d}\omega_{f^{*}}(N^{-2/d}M^{-2/d}),$

where $C>0$ is a constant which does not depend on
$n,d,\mathcal{B},\mathcal{S},\mathcal{D},\mathcal{W},\lambda_{L},N$ or $M$.
Here $\lambda_{L}$ is the Lipschitz constant of the robust loss function $L$
under Assumption 1, and $\omega_{f^{*}}$ is the modulus of continuity of the
uniformly continuous target $f^{*}$ under Assumption 3.

Moreover, aside from Assumptions 1-3, if Assumption 4 also holds, then we have

$\mathbb{E}\\{\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\\}\leq
C\lambda_{L}\mathcal{B}\frac{\mathcal{S}\mathcal{D}\log(\mathcal{S})\log(n)}{n^{1-1/p}}+384\lambda_{L,f^{*}}d\big{\\{}\omega_{f^{*}}(N^{-2/d}M^{-2/d})\big{\\}}^{2},$

where
$\lambda_{L,f^{*}}=\max\\{\Lambda_{L,f^{*}},\lambda_{L}/\varepsilon_{L,f^{*}}\\}>0$
is a constant defined in Lemma 3.3.

Theorem 4.2 implies that the upper bound of the excess risk is a sum of two
terms: the upper bound of the stochastic error
$C\lambda_{L}\mathcal{B}\mathcal{S}\mathcal{D}\log(\mathcal{S})\log(n)/{n^{1-1/p}}$
and that of the approximation error
$18\lambda_{L}\sqrt{d}\omega_{f^{*}}(N^{-2/d}M^{-2/d})$. A few important
comments are in order. First, our error bound is non-asymptotic and explicit
in the sense that no unclearly defined constant is involved. Importantly, the
prefactor $18\lambda_{L}\sqrt{d}$ in the upper bound of approximation error
depends on $d$ through $\sqrt{d}$ (the root of the dimension $d$), which is
quite different from the exponential dependence in many existing results.
Second, the approximation rate $\omega_{f^{*}}(N^{-2/d}M^{-2/d})$ is in terms
of the network width $\mathcal{W}=\max\\{4d\lfloor N^{1/d}\rfloor+3d,12N+8\\}$
and depth $\mathcal{D}=12M+14$, rather than only the network size
$\mathcal{S}$. This provides insights into the relative merits of different
network architecture and provides certain qualitative guidance on the network
design. Third, our result on the approximation rate accommodates all kinds of
uniformly continuous target functions $f^{*}$ and the result is expressed in
terms of the modulus of continuity of $f^{*}$.

To achieve the optimal rate, we need to balance the trade-off between the
stochastic error and the approximation error. On one hand, the upper bound on
the stochastic error
$C\lambda_{L}\mathcal{B}\mathcal{S}\mathcal{D}\log(\mathcal{S})\log(n)/{n^{1-1/p}}$
increases in the complexity and richness of the function class
$\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$;
larger $\mathcal{D}$, $\mathcal{S}$ and $\mathcal{B}$ lead to a larger upper
bound on the stochastic error. On the other hand, the upper bound for the
approximation error $18\lambda_{L}\sqrt{d}\omega_{f^{*}}(N^{-2/d}M^{-2/d})$
decreases as the size of
$\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
increases, larger $\mathcal{D}$ and $\mathcal{W}$ result in smaller upper
bound on the approximation error.

Apart form the excess risk bounds, if the loss function is self-calibrated,
the distance between the estimate and the target can also be upper bounded via
the self-calibration inequality.

######  Assumption 5 (Self-calibration).

There exists a constant $C_{L,f^{*}}>0$ which may depend on the loss function
$L(\cdot)$, the risk minimizer $f^{*}$ and the distribution of $(X,Y)$ such
that

$\Delta^{2}(f,f^{*})\leq
C_{L,f^{*}}\big{\\{}\mathcal{R}(f)-\mathcal{R}(f^{*})\big{\\}}$

for any $f:\mathcal{X}\to\mathbb{R}$, where $\Delta^{2}(\cdot,\cdot)$ is some
distance defined in the functional space $\\{f:\mathcal{X}\to\mathbb{R}\\}$
such that $\Delta^{2}(f_{1},f_{2})=0$ if
$\|f_{1}-f_{2}\|_{L^{\infty}(\nu)}=0$.

######  Remark 5.

For least squares estimation, Assumption 5 naturally holds with
$C_{L,f^{*}}=1$ and
$\Delta^{2}(f,f^{*})=\|f-f^{*}\|^{2}_{L^{2}(\nu)}=\mathcal{R}(f)-\mathcal{R}(f^{*})$.
For least absolute deviation or quantile check loss, under mild conditions,
the self-calibration assumption is satisfied with
$\Delta^{2}(f,f^{*})=\mathbb{E}\big{[}\min\\{|f(X)-f^{*}(X)|,|f(X)-f^{*}(X)|^{2}\\}\big{]}$
and some $C_{L,f^{*}}>0$. More details can be found in Christmann and
Steinwart (2007); Steinwart et al. (2011); Lv et al. (2018) and Padilla,
Tansey and Chen (2020).

With the self-calibration condition satisfied, the distance between the
estimate and the target is upper bounded as a direct consequence of Theorem
4.2.

######  Corollary 1.

Suppose that Assumptions 1-3 and 5 hold, $\nu$ is absolutely continuous with
respect to the Lebesgue measure, and $\|f^{*}\|_{\infty}\leq\mathcal{B}$ for
some $\mathcal{B}\geq 1$. Then, for any $N,M\in\mathbb{N}^{+}$, the function
class of ReLU multi-layer perceptrons
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
with depth $\mathcal{D}=12M+14$ and width $\mathcal{W}=\max\\{4d\lfloor
N^{1/d}\rfloor+3d,12N+8\\}$, for $2n\geq\text{Pdim}(\mathcal{F}_{\phi})$, the
expected distance between $f^{*}$ and the ERM $\hat{f}_{\phi}$ defined in
(3.2) satisfies

$\mathbb{E}\big{\\{}\Delta^{2}(\hat{f}_{\phi},f^{*})\big{\\}}\leq
C_{L,f^{*}}\Big{\\{}C\lambda_{L}\mathcal{B}\frac{\mathcal{S}\mathcal{D}\log(\mathcal{S})\log(n)}{n^{1-1/p}}+18\lambda_{L}\sqrt{d}\omega_{f^{*}}(N^{-2/d}M^{-2/d})\Big{\\}},$

where $C_{L,f^{*}}$ is defined in Assumption 5 and $C>0$ is a constant which
does not depend on
$n,d,\mathcal{B},\mathcal{S},\mathcal{D},\mathcal{W},\lambda_{L},N$ or $M$.

Additionally, other than Assumptions 1-3 and 5, if Assumption 4 also holds, we
have

$\mathbb{E}\big{\\{}\Delta^{2}(\hat{f}_{\phi},f^{*})\big{\\}}\leq
C_{L,f^{*}}\Big{[}C\lambda_{L}\mathcal{B}\frac{\mathcal{S}\mathcal{D}\log(\mathcal{S})\log(n)}{n^{1-1/p}}+384\lambda_{L,f^{*}}d\big{\\{}\omega_{f^{*}}(N^{-2/d}M^{-2/d})\big{\\}}^{2}\Big{]},$

where
$\lambda_{L,f^{*}}=\max\\{\Lambda_{L,f^{*}},\lambda_{L}/\varepsilon_{L,f^{*}}\\}>0$
is a constant defined in Lemma 3.3.

###  4.2 On the curse of dimensionality

In many contemporary statistical and machine learning applications, the
dimensionality $d$ of the input $X$ can be large, which leads to an extremely
slow rate of convergence, even with huge sample size. This problem is well
known as the curse of dimensionality. A promising and effective way to
mitigate the curse of dimensionality is to impose additional conditions on $X$
and the target function $f^{*}$. Though the domain of $f^{*}$ is high
dimensional, when the support of $X$ is concentrated on some neighborhood of a
low-dimensional manifold, the upper bound of the approximation error can be
substantially improved in terms of the exponent of the convergence rate (Shen,
Yang and Zhang, 2020).

######  Assumption 6.

The input $X$ is supported on $\mathcal{M}_{\rho}$, a $\rho$-neighborhood of
$\mathcal{M}\subset[0,1]^{d}$, where $\mathcal{M}$ is a compact
$d_{\mathcal{M}}$-dimensional Riemannian submanifold (Lee, 2006) and

$\mathcal{M}_{\rho}=\\{x\in[0,1]^{d}:\inf\\{\|x-y\|_{2}:y\in\mathcal{M}\\}\leq\rho\\}$

for $\rho\in(0,1)$.

Assumption 6 is more realistic than the exact manifold support assumption
assumed in Schmidt-Hieber (2019),Nakada and Imaizumi (2019) and Chen et al.
(2019). In practical application, one can rarely observe data that are located
on an exact manifold. It is more reasonable to assume that the data are
concentrated on a neighborhood of a low-dimensional manifold.

######  Theorem 4.3 (Non-asymptotic error bound).

Under model (1.1), suppose that Assumptions 1-3 and 6 hold, the probability
measure $\nu$ of $X$ is absolutely continuous with respect to the Lebesgue
measure and $\|f_{0}\|_{\infty}\leq\mathcal{B}$ for some $\mathcal{B}\geq 1$.
Then for any $N,M\in\mathbb{N}^{+}$, the function class of ReLU multi-layer
perceptrons
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
with depth $\mathcal{D}=12M+14$ and width
$\mathcal{W}=\max\\{4d_{\delta}\lfloor
N^{1/d_{\delta}}\rfloor+3d_{\delta},12N+8\\}$, the expected excess risk of the
ERM $\hat{f}_{\phi}$ defined in (3.2) satisfies

$\mathbb{E}\\{\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\\}\leq
C_{1}\lambda_{L}\mathcal{B}\frac{\mathcal{S}\mathcal{D}\log(\mathcal{S})\log(n)}{n^{1-1/p}}+\lambda_{L}(2+18\sqrt{d_{\delta}}\,)\,\omega_{f^{*}}\big{\\{}(C_{2}+1)(NM)^{-2/d_{\delta}}\big{\\}},$

for $2n\geq\text{Pdim}(\mathcal{F}_{\phi})$ and $\rho\leq
C_{2}(NM)^{-2/d_{\delta}}(1-\delta)/\\{2(\sqrt{d/d_{\delta}}+1-\delta)\\}$,
where $d_{\delta}=O(d_{\mathcal{M}}{\log(d/\delta)}/{\delta^{2}})$ is an
integer such that $d_{\mathcal{M}}$$\leq d_{\delta}<d$ for any
$\delta\in(0,1)$ and $C_{1},C_{2}>0$ are constants which do not depend on
$n,d,\mathcal{B},\mathcal{S},\mathcal{D},\mathcal{W},\lambda_{L},N$ or $M$.
Here $\lambda_{L}$ is the Lipschitz constant of the general robust loss
function $L$ under Assumption 1, and $\omega_{f^{*}}$ is the modulus of
continuity of the uniformly continuous target $f^{*}$ under Assumption 3.

Additionally, aside from Assumptions 1-3 and 6, if Assumption 4 also holds,
then

$\displaystyle\mathbb{E}\\{\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\\}\leq$
$\displaystyle
C_{1}\lambda_{L}\mathcal{B}\frac{\mathcal{S}\mathcal{D}\log(\mathcal{S})\log(n)}{n^{1-1/p}}$
$\displaystyle\qquad+\lambda_{L,f^{*}}\Big{[}(2+18\sqrt{d_{\delta}}\,)\,\omega_{f^{*}}\big{\\{}(C_{2}+1)(NM)^{-2/d_{\delta}}\big{\\}}\Big{]}^{2},$

where
$\lambda_{L,f^{*}}=\max\\{\Lambda_{L,f^{*}},\lambda_{L}/\varepsilon_{L,f^{*}}\\}>0$
is a constant defined in Lemma 3.3.

Furthermore, if Assumption 5 also holds, then

$\displaystyle\mathbb{E}\\{\Delta^{2}(\hat{f}_{\phi},f^{*})\\}\leq
C_{L,f^{*}}\Bigg{(}C_{1}\lambda_{L}\mathcal{B}$
$\displaystyle\frac{\mathcal{S}\mathcal{D}\log(\mathcal{S})\log(n)}{n^{1-1/p}}$
$\displaystyle+\lambda_{L,f^{*}}\Big{[}\big{(}2+18\sqrt{d_{\delta}}\,\big{)}\,\omega_{f^{*}}\big{\\{}(C_{2}+1)(NM)^{-2/d_{\delta}}\big{\\}}\Big{]}^{2}\Bigg{)},$

where $C_{L,f^{*}}>0$ is a constant defined in Assumption 5.

Theorem 4.3 is established under the assumption that the distribution of $X$
is supported on an approximate Riemannian manifold with an intrinsic dimension
lower than the dimension $d$ of the ambient space $\mathbb{R}^{d}$. This is
different from the hierarchical structure assumption on the target function in
Bauer and Kohler (2019) and Schmidt-Hieber (2020). These two assumptions are
two different lines of assumptions on the underlying target function, and
either one can mitigate the curse of dimensionality.

##  5 Error bounds for different network structures

Theorem 4.2 provides a general description of how the non-asymptotic error
bounds depend on the network parameters. In this section, we consider the
error bounds for some specific network structures. We also consider the
relative efficiency of networks with different shapes in terms of the size of
the network needed to achieve the optimal error bound. The calculations given
below demonstrate the advantages of deep networks over shallow ones in the
sense that deep networks can achieve the same error bound as the shallow
networks with a fewer total number of parameters in the network.

To give concrete examples on the convergence rate of the excess risk, we
assume the target function $f^{*}$ to be Hölder continuous, i.e. there exist
$\theta\geq 0$ and $\alpha>0$ such that

$|f^{*}(x)-f^{*}(y)|\leq\theta\|x-y\|_{2}^{\alpha},\quad
x,y\in\mathcal{X}=[0,1]^{d}.$ (5.1)

Let $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ be the sizes of two neural
networks $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ needed to achieve the same
non-asymptotic error bound as given in Theorem 4.2. We define the relative
efficiency between two networks $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ as

$\text{REN}(\mathcal{N}_{1},\mathcal{N}_{2})=\frac{\log\mathcal{S}_{2}}{\log\mathcal{S}_{1}}.$
(5.2)

Here we use the logarithm of the size because the size of the network for
achieving the optimal error rate usually has the form $\mathcal{S}=n^{s}$ for
some $s>0$ up to a factor only involving the power of $\log n$, as will be
seen below. Let $r=\text{REN}(\mathcal{N}_{1},\mathcal{N}_{2})$. In terms of
sample sizes, this definition of relative efficiency implies that, if it takes
a sample of size $n$ for network $\mathcal{N}_{1}$ to achieve the optimal
error rate, then it will take a sample of size $n^{r}$ for network
$\mathcal{N}_{2}$ to achieve the same error rate.

For any multilayer neural network in
$\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$,
its parameters naturally satisfy

$\max\\{\mathcal{W},\mathcal{D}\\}\leq\mathcal{S}\leq\mathcal{W}(d+1)+(\mathcal{W}^{2}+\mathcal{W})(\mathcal{D}-1)+\mathcal{W}+1=O(\mathcal{W}^{2}\mathcal{D}).$
(5.3)

Based on the relationship in (5.3) and Theorem 4.2, we next present Corollary
2 for the networks with some specific structures without assuming Assumption
4. Similar corollaries with Assumption 4 satisfied can be established
accordingly, and the relative efficiency of different networks can be
compared. Denote

$n_{*}=n_{d,p,\alpha}=n^{\left(1-\frac{1}{p}\right)\frac{d}{d+\alpha}}.$ (5.4)

Recall if the response $Y$ is sub-exponential, we take $p=+\infty.$ For the
function class of ReLU multi-layer perceptrons
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
with depth $\mathcal{D}$, width $\mathcal{W}$ and size $\mathcal{S}$ given by

  1. (a)

Deep with fixed width networks ($\mathcal{N}_{\text{DFW}}$):

$\displaystyle\mathcal{D}_{1}$ $\displaystyle=12\lfloor n_{*}^{1/2}(\log
n)^{-1}\rfloor+14,\mathcal{W}_{1}=\max\\{7d,20\\},\mathcal{S}_{1}=O(n_{*}^{1/2}(\log
n)^{-1}),$

  2. (b)

Wide and fixed depth networks ($\mathcal{N}_{\text{WFD}}$):

$\displaystyle\mathcal{D}_{2}$
$\displaystyle=26,\mathcal{W}_{2}=\max\\{4d\lfloor
n_{*}^{1/(2d)}\rfloor+3d,12\lfloor
n_{*}^{1/2}\rfloor+8\\},\mathcal{S}_{2}=O(n_{*}(\log n)^{-1}),$

  3. (c)

Deep and wide networks ($\mathcal{N}_{\text{DAW}}$):

$\displaystyle\mathcal{D}_{3}$
$\displaystyle=O(n_{*}^{1/4}),\mathcal{W}_{3}=O(n_{*}^{1/4}),\mathcal{S}_{3}=O(n_{*}^{3/4}(\log
n)^{-2}).$

######  Corollary 2 (Error bounds).

Let $f^{*}$ defined in (2.1) be Hölder continuous of order $\alpha>0$ with
Hölder constant $\theta\geq 0$. Suppose that Assumptions 1-3 hold, $\nu$ is
absolutely continuous with respect to the Lebesgue measure, and
$\|f^{*}\|_{\infty}\leq\mathcal{B}$ for some $\mathcal{B}\geq 1$. Then, for
any one of the three types of networks
$\mathcal{N}_{\text{DFW}},\mathcal{N}_{\text{WFD}}$ and
$\mathcal{N}_{\text{DAW}}$ specified above, the ERM $\hat{f}_{\phi}$ defined
in (3.2) satisfies

$\mathbb{E}\\{\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\\}\leq\lambda_{L}\left\\{C\mathcal{B}+18\theta\sqrt{d}\right\\}n^{-(1-\frac{1}{p})\frac{\alpha}{d+\alpha}},$

for $2n\geq\text{Pdim}(\mathcal{F}_{\phi})$, where $C>0$ is a constant which
does not depend on $n,\lambda,\alpha$ or $\mathcal{B}$, and $C$ does not
depend on $d$ and $N$ if $d\lfloor N^{1/d}\rfloor+3d/4\leq 3N+2$, otherwise
$C=O(d^{2}\lfloor N^{2/d}\rfloor)$. Additionally, if Assumption 5 holds, then

$\mathbb{E}\\{\Delta^{2}(\hat{f}_{\phi},f^{*})\\}\leq\lambda_{L}C_{L,f^{*}}\left\\{C\mathcal{B}+18\theta\sqrt{d}\right\\}n^{-(1-\frac{1}{p})\frac{\alpha}{d+\alpha}},$

where $C_{L,f^{*}}>0$ is a constant defined in Assumption 5.

Corollary 2 is a direct consequence of Theorem 4.2 and Corollary 1. By
Corollaries 2, the sizes of the three types of networks
$\mathcal{N}_{\text{DFW}},\mathcal{N}_{\text{WDF}}$ and
$\mathcal{N}_{\text{DAW}}$ to achieve the same error rate are

$\mathcal{S}_{\text{DFW}}=O(n_{*}^{1/2}(\log n)^{-1}),\
\mathcal{S}_{\text{WFD}}=O(n_{*}(\log n)^{-1}),\ \text{ and }\
\mathcal{S}_{\text{DAW}}=O(n_{*}^{3/4}(\log n)^{-1}),$

respectively. So we have the relationship

$\mathcal{S}_{\text{DFW}}\approx\mathcal{S}_{\text{WFD}}^{1/2}\approx\mathcal{S}_{\text{DAW}}^{2/3}.$

The relative efficiency of these two networks as defined in (5.2) are

$\text{REN}(\mathcal{N}_{\text{DAW}},\mathcal{N}_{\text{DFW}})=\frac{1/2}{3/4}=\frac{2}{3},\text{REN}(\mathcal{N}_{\text{DAW}},\mathcal{N}_{\text{WFD}})=\frac{1}{3/4}=\frac{4}{3}.$

Thus the deep networks are twice as efficient as wide networks. If the sample
size for a deep with fixed width network to achieve the optimal error rate is
$n$, then it is about $n^{2}$ for a wide with fixed depth network and
$n^{4/3}$ for a deep and wide network. An explanation is that deep neural
networks are generally of greater approximation power than shallow networks.
In Telgarsky (2016), it was shown that for any integer $k\geq 1$ and dimension
$d\geq 1$, there exists a function computed by a ReLU neural network with
$2k^{3}+8$ layers, $3k^{2}+12$ neurons and $4+d$ different parameters such
that it can not be approximated by networks activated by piecewise polynomial
functions with no more than $k$ layers and $O(2^{k})$ neurons.

###  5.1 Optimal design of rectangle networks

We now discuss the optimal design of rectangle networks, i.e., networks with
equal width for each hidden layer. For such networks with a regular shape, we
have an exact relationship between the size of the network, the depth and the
width:

$\mathcal{S}=\mathcal{W}(d+1)+(\mathcal{W}^{2}+\mathcal{W})(\mathcal{D}-1)+\mathcal{W}+1=O(\mathcal{W}^{2}\mathcal{D}).$
(5.5)

Based on this relationship and Theorem 4.2, we can determine the depth and the
width of the network to achieve the optimal rate with the minimal size. Below,
we consider the best selection of network depth and width under Assumption 4.

Under Assumption 4 holds, to achieve the optimal rate with respect to the
sample size $n$ with a minimal network size, we can set

$\displaystyle\mathcal{W}$ $\displaystyle=$ $\displaystyle\max(7d,20),$
$\displaystyle\mathcal{D}$ $\displaystyle=$ $\displaystyle 12\lfloor
n^{(1-1/p)d/{(2d+4\alpha)}}(\log n)^{-1}\rfloor+14,$
$\displaystyle\mathcal{S}$ $\displaystyle=$
$\displaystyle\big{\\{}\max(49d^{2},400)+\max(7d,20)\big{\\}}\times\big{\\{}12\lfloor
n^{(1-1/p)d/{(2d+4\alpha)}}(\log n)^{-1}\rfloor+14\big{\\}}$
$\displaystyle+(d+2)\max(7d,20)+1$ $\displaystyle=$ $\displaystyle
O(49d^{2}n^{(1-1/p)d/{(2d+4\alpha)}}(\log n)^{-1}).$

Furthermore, if Assumption 6 is satisfied, we can take $\mathcal{F}_{\phi}$ to
be consisting of fixed-width networks with

$\displaystyle\mathcal{W}$ $\displaystyle=$
$\displaystyle\max(7d_{\delta},20),$ $\displaystyle\mathcal{D}$
$\displaystyle=$ $\displaystyle 12\lfloor
n^{(1-1/p)d_{\delta}/{(2d_{\delta}+4\alpha)}}(\log n)^{-1}\rfloor(\log
n)^{-1}+14,$ $\displaystyle\mathcal{S}$ $\displaystyle=$
$\displaystyle\big{\\{}\max(49d_{\delta}^{2},400)+\max(7d_{\delta},20)\big{\\}}\times\big{\\{}12\lfloor
n^{(1-1/p)d_{\delta}/{(2d_{\delta}+4\alpha)}}\rfloor+14\big{\\}}$
$\displaystyle+(d_{\delta}+2)\max(7d_{\delta},20)+1$ $\displaystyle=$
$\displaystyle
O(49d_{\delta}^{2}n^{(1-1/p)d_{\delta}/{(2d_{\delta}+4\alpha)}}(\log
n)^{-1}).$

Then the excess risk of $\hat{f}_{\phi}$ in Theorem 4.3 satisfies

$\mathbb{E}\\{\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\\}\leq\\{C_{1}\lambda_{L}\mathcal{B}+\lambda_{L,f^{*}}(2+18\sqrt{d_{\delta}})(C_{2}+1)^{\alpha}\theta\\}n^{-\left(1-\frac{1}{p}\right)\frac{2\alpha}{d_{\delta}+2\alpha}},$

where $C_{2}>0$ is a constant independent of
$n,d_{\delta},\mathcal{B},\mathcal{S},\mathcal{D},\lambda$ and $\alpha$, and
$C_{1}>0$ is a constant that does not depend on
$n,\mathcal{B},\mathcal{S},\mathcal{D}$ or $\lambda$, and does not depend on
$d_{\delta}$ if $7d_{\delta}\leq 20$, otherwise $C_{1}=O(d_{\delta}^{2})$.

Additionally, if Assumption 5 holds, then

$\mathbb{E}\\{\Delta^{2}(\hat{f}_{\phi},f^{*})\\}\leq
C_{L,f^{*}}\\{C_{1}\lambda_{L}\mathcal{B}+\lambda_{L,f^{*}}(2+18\sqrt{d_{\delta}})(C_{2}+1)^{\alpha}\theta\\}n^{-\left(1-\frac{1}{p}\right)\frac{2\alpha}{d_{\delta}+2\alpha}},$

where $C_{L,f^{*}}>0$ is a constant defined in Assumption 5.

##  6 Experiments

We apply the deep nonparametric regression with different robust loss
functions to the simulated data. We compare the results with the nonparametric
least squares estimates under different settings of error distributions. All
the deep network computations are implemented on Python and mainly via
PyTorch.

###  6.1 Network configuration

In all settings, we implement the empirical risk minimization by ReLU
activated fixed width multilayer perceptrons. We consider a class of functions
(implemented by fixed width networks) consisting of ReLU activated multilayer
perceptrons with 5 hidden layers and network width being
$(d,256,256,256,256,256,1)$, denoted by Nets-256. All weights and bias in each
layer are initialized by uniformly samples on bounded intervals according to
the default initialization mechanism in PyTorch.

The optimization algorithm is Adam (Kingma and Ba, 2014) with default learning
rate $0.01$ and default $\beta=(0.9,0.99)$ (coefficients used for computing
running averages of gradient and its square). We set the batch size to be
$n/4$ in all cases where $n$ is the size of the training data, and train the
network for at least $400$ epochs until the training loss converges or becomes
sufficient steady.

###  6.2 Estimations and Evaluations

We consider $5$ different loss functions: (i) Least square loss function,
denoted by LS; (ii) Least absolute deviation loss function, denoted by LAD;
(iii) Huber loss function with parameter $\zeta=1.345$, denoted by Huber; (iv)
Cauchy loss function with parameter $\kappa=1$, denoted by Cauchy; (v) Tukey’s
biweight loss function with parameter $t=4.685$.

For each configuration, we use each loss function above as the training loss
$L_{train}$ to obtain the ERM $\hat{f}_{train}$ within Nets-256,

$\hat{f}_{train}=\arg\min_{f\in{\rm\textit{Nets-256}}}\frac{1}{n}\sum_{i=1}^{n}L_{train}(f(X^{train}_{i}),Y^{train}_{i}),$

where $(X^{train}_{i},Y^{train}_{i})_{i=1}^{n}$ is the training data. For each
training loss function, we calculate the testing risks under five different
choices of loss function $L_{test}$ on the testing data
$(X^{test}_{t},Y^{test}_{t})_{t=1}^{T}$, that is

$\mathcal{R}(\hat{f}_{train},L_{test})=\frac{1}{T}\sum_{t=1}^{T}L_{test}(\hat{f}_{train}(X^{test}_{t}),Y^{test}_{t}).$

We report the averaged testing risk and its standard deviation over $R=10$
replications under different scenarios. For each data generation model, the
testing data $(X^{test}_{t},Y^{test}_{t})_{t=1}^{T}$, $T=100,000$, are
generated to have the same distribution as the training data
$(X^{train}_{i},Y^{train}_{i})_{i=1}^{n}$.

###  6.3 Data generation: univariate models

In our simulation, we first consider four univariate target functions,
including “Blocks”, “Bumps”, “Heavisine” and “Dopller” (Donoho and Johnstone,
1994). The formulae for these functions are given below.

  1. 1.

Blocks:

$f_{0}(x)=\sum h_{i}I(x>t_{i}),$

where

$\displaystyle(h_{i})$
$\displaystyle=(4,-5,-2.5,4,-3,2.1,4.3,-1.1,-2.1,-4.2),$
$\displaystyle(t_{i})$
$\displaystyle=(0.1,0.15,0.23,0.28,0.40,0.44,0.65,0.76,0.78,0.81).$

  2. 2.

Bumps:

$f_{0}(x)=\sum h_{i}(1+|x-t_{i}|/w_{i})^{-4},$

where

$\displaystyle(h_{i})$ $\displaystyle=(4,5,2.5,4,3,2.1,4.3,1.1,2.1,4.2),$
$\displaystyle(t_{i})$
$\displaystyle=(0.1,0.15,0.23,0.28,0.40,0.44,0.65,0.76,0.78,0.81),$
$\displaystyle(w_{i})$
$\displaystyle=(0.005,0.005,0.006,0.01,0.01,0.03,0.01,0.01,0.005,0.008).$

  3. 3.

Heavisine:

$f_{0}(x)=4\sin(4\pi x)-{\rm sgn}(x-0.3)-{\rm sgn}(0.72-x).$

  4. 4.

Dopller:

$f_{0}(x)=\\{x(1-x)\\}^{1/2}\sin\\{2.2\pi/(x+0.15)\\}.$

The four target functions are depicted in Fig 3.

Figure 3: Four univariate target functions including “Block”, “Bump”,
“Heavisine” and “Dopller”.

We generate the training data $(X_{i}^{train},Y_{i}^{train})_{i=1}^{n}$,
$n=128$ or $512$ according to model (1.1). The predictor $X^{train}$ is
generated from Uniform$[0,1]$. For the error $\eta$, we have tried several
distributions: (i) Standard normal distribution, denoted by
$\eta\sim\mathcal{N}(0,1)$; (ii) Student’s $t$ distribution with $2$ degrees
of freedom, denoted by $\eta\sim t(2)$; (iii) Standard Cauchy
distribution(with location parameter 0 and scale parameter 1), denoted by
$\eta\sim{\rm Cauchy(0,1)}$; (iv) Normal mixture distribution, denoted by
“Mixture”, where $\eta\sim\xi\mathcal{N}(0,1)+(1-\xi)\mathcal{N}(0,10^{4})$
with $\xi=0.8$ This is the contaminated normal model often used to depict the
situation when $100\xi\%$ of the observations are corrupt (Huber, 1964). Note
that Assumption 2 is satisfied respectively in the way of (i)
$\eta\sim\mathcal{N}(0,1)$, $\mathbb{E}|\eta|^{p}<\infty$ for
$p\in[0,\infty)$; (ii) $\eta\sim t(2)$, $\mathbb{E}|\eta|^{p}<\infty$ for
$p\in[0,2)$; (iii) $\eta\sim{\rm Cauchy(0,1)}$, $\mathbb{E}|\eta|^{p}<\infty$
for $p\in[0,1)$; (iv) $\eta$ follows normal mixture distribution,
$\mathbb{E}|\eta|^{p}<\infty$ for $p\in[0,\infty)$. In addition, the error
distributions are symmetric; together with the symmetric property of the five
loss functions, the underlying true function $f_{0}$ in (1.1) coincides with
the target function $f^{*}$ in (2.1) according to the definition of the risk
minimizer.

The simulation results are shown in Table 2-5. For normal errors, except for
“Heavisine” model, the least squares estimation generally outperforms other
methods in the sense that least squares estimator has smaller mean testing
excess risk and standard deviation. However, with heavier-tailed errors or
contaminated error, the performance of the least squares estimator becomes
worse or even breaks down, while robust regression methods with LAD, Huber or
Cauchy loss give more accurate and stable results. Therefore, these robust
methods possesses the desired robustness property with respect to heavy-tailed
or contaminated error terms in nonparametric estimations. The fitted curves
are presented in Figure 4 and the training processes are displayed in Figures
5-8 in Appendix.

###  6.4 Data generation: multivariate models

To simulate a multivariate target function $f^{*}:[0,1]^{d}\to\mathbb{R}$, we
refer to the Kolmogorov-Arnold (KA) representation theorem (Hecht-Nielsen,
1987), which says that every continuous function can be represented by a
specific network with two hidden layers. We formulate the original version of
the KA representation theorem according to Schmidt-Hieber (2021): for any
continuous function $f^{*}:[0,1]^{d}\to\mathbb{R}$, there exist univariate
continuous functions $g_{k},\psi_{m,k}$ such that

$f^{*}(x_{1},\ldots,x_{d})=\sum_{k=0}^{2d}g_{k}\big{(}\sum_{m=1}^{d}\psi_{m,k}(x_{m})\big{)}.$

In light of this, these $(2d+1)(d+1)$ univariate functions suffices for an
exact representation of a $d$-variate function. We consider $f^{*}$ according
to the KA representation by choosing $g_{k},\psi_{m,k},k=0,\ldots,2d,\
m=1,\ldots,d$ from a pool of univariate functions
$\mathcal{H}=\\{h_{i}(x),i=1,\ldots,7\\}$, where
$h_{1}(x)=-2.2x+0.3,h_{2}(x)=0.7x^{3}-0.2x^{2}+0.3x-0.3,h_{3}(x)=0.3\times{\rm
sign}(x)\sqrt{|x|},h_{4}(x)=0.8\times\log(|x|+0.01),h_{5}(x)=\exp(\min\\{0.2x-0.1,4\\}),h_{6}(x)=\sin(6.28x),h_{7}(x)={2}/({|x|+0.1}).$
For a given dimensionality $d$, we first generate $(2d+1)(d+1)$ random
integers uniformly distributed on $\\{1,2,3,4,5,6,7\\}$ as indices of the
functions $g_{k},\psi_{m,k}$; then we specify the choices of
$g_{k},\psi_{m,k}$ based on these randomly generated indices. To make our
simulation reproducible, the seed is set as $2021$ in the random number
generation, so the choices of $g_{k},\psi_{m,k}$ can be reproducible.

We generate $X$ from the uniform distribution on $[0,1]^{d}$, i.e. $X\sim{\rm
Unif}[0,1]^{d}$. We set the dimensionality $d$ of the input $X$ to be $4$ and
generate the training data $(X^{train}_{i},Y^{train}_{i})_{i=1}^{n}$ with
sample size $n=128$ and $512$ for each of $R=10$ replication. The results are
summarized in Tables 6. It can be seen that except for the case of normally
distributed error, estimators based on robust loss functions always outperform
least squares estimator in terms of the mean excess risk and standard
deviation.

Table 2: The testing excess risks and standard deviation (in parentheses)
across different loss functions of the ERMs trained with respect to different
loss functions using Nets-256. Data is generated from ”Blocks” model with
training sample size $n=128$ or $512$ and the number of replications $R=10$.

Blocks, $n=128$ | Testing Loss
---|---
Error | Training Loss | LS | LAD | Huber | Cauchy | Tukey
$N(0,1)$ | LS | 13.534(4.078) | 1.408(0.348) | 1.363(0.340) | 0.766(0.180) | 0.784(0.182)
LAD | 18.192(6.686) | 1.708(0.285) | 1.655(0.281) | 0.890(0.128) | 0.891(0.153)
Huber | 18.178(6.367) | 1.732(0.449) | 1.681(0.445) | 0.890(0.196) | 0.871(0.195)
Cauchy | 42.140(13.466) | 2.571(0.655) | 2.531(0.653) | 0.967(0.182) | 0.755(0.112)
Tukey | 159.659(15.485) | 8.089(0.683) | 8.020(0.677) | 2.516(0.281) | 1.645(0.225)
$t(2)$ | LS | 19.057(11.877) | 1.770(0.487) | 1.713(0.479) | 0.932(0.198) | 0.935(0.202)
LAD | 28.137(15.167) | 2.160(0.695) | 2.110(0.694) | 0.974(0.180) | 0.873(0.095)
Huber | 21.393(7.411) | 1.725(0.338) | 1.679(0.336) | 0.807(0.132) | 0.741(0.134)
Cauchy | 35.695(9.973) | 2.093(0.450) | 2.057(0.449) | 0.782(0.126) | 0.612(0.081)
Tukey | 161.804(18.701) | 7.989(0.585) | 7.925(0.584) | 2.383(0.112) | 1.459(0.092)
$Cauchy(0,1)$ | LS | 1006121.625(2722768.500) | 187.531(486.105) | 187.457(486.096) | 2.914(3.326) | 1.481(0.516)
LAD | 26.255(16.659) | 1.778(0.692) | 1.735(0.687) | 0.782(0.213) | 0.687(0.149)
Huber | 15.564(17.321) | 1.717(0.475) | 1.672(0.473) | 0.774(0.142) | 0.692(0.108)
Cauchy | 31.284(13.132) | 2.057(0.399) | 2.018(0.400) | 0.785(0.101) | 0.611(0.066)
Tukey | 156.523(18.607) | 7.206(0.624) | 7.145(0.622) | 2.147(0.135) | 1.314(0.108)
Mixture | LS | 372.860(109.658) | 9.262(1.380) | 9.170(1.376) | 2.758(0.276) | 2.138(0.134)
LAD | 42.559(16.447) | 2.264(0.458) | 2.218(0.451) | 0.943(0.163) | 0.832(0.164)
Huber | 29.233(12.558) | 1.853(0.529) | 1.806(0.526) | 0.864(0.173) | 0.814(0.131)
Cauchy | 34.202(8.009) | 1.842(0.324) | 1.807(0.322) | 0.744(0.117) | 0.626(0.111)
Tukey | 148.123(19.258) | 6.439(0.683) | 6.386(0.679) | 1.974(0.227) | 1.280(0.171)
Blocks, $n=512$ | LS | LAD | Huber | Cauchy | Tukey
$N(0,1)$ | LS | 7.547(2.020) | 0.948(0.231) | 0.908(0.224) | 0.561(0.134) | 0.596(0.148)
LAD | 10.319(3.478) | 1.151(0.281) | 1.109(0.275) | 0.640(0.134) | 0.655(0.128)
Huber | 9.470(3.630) | 1.051(0.349) | 1.013(0.344) | 0.581(0.176) | 0.590(0.175)
Cauchy | 29.639(4.070) | 1.720(0.229) | 1.693(0.229) | 0.626(0.087) | 0.474(0.070)
Tukey | 166.059(10.031) | 7.957(0.281) | 7.903(0.282) | 2.276(0.016) | 1.322(0.027)
$t(2)$ | LS | 8.003(2.626) | 1.047(0.290) | 1.000(0.277) | 0.644(0.174) | 0.708(0.191)
LAD | 8.197(2.587) | 1.045(0.279) | 1.003(0.271) | 0.624(0.161) | 0.669(0.174)
Huber | 8.649(3.074) | 1.043(0.359) | 1.002(0.346) | 0.611(0.221) | 0.654(0.253)
Cauchy | 28.553(5.855) | 1.702(0.358) | 1.674(0.356) | 0.634(0.127) | 0.497(0.090)
Tukey | 171.790(5.422) | 8.442(0.169) | 8.379(0.166) | 2.499(0.093) | 1.546(0.101)
$Cauchy(0,1)$ | LS | 214.717(529.612) | 3.321(2.587) | 3.256(2.580) | 1.291(0.396) | 1.191(0.280)
LAD | 11.284(3.994) | 1.004(0.272) | 0.969(0.267) | 0.530(0.121) | 0.525(0.114)
Huber | 14.093(3.489) | 0.981(0.180) | 0.947(0.180) | 0.506(0.068) | 0.491(0.067)
Cauchy | 26.687(5.734) | 1.322(0.256) | 1.298(0.255) | 0.491(0.070) | 0.381(0.047)
Tukey | 183.687(23.299) | 7.052(0.809) | 7.007(0.804) | 1.897(0.259) | 1.058(0.164)
Mixture | LS | 102.667(24.698) | 5.587(0.766) | 5.498(0.763) | 2.313(0.227) | 2.019(0.160)
LAD | 13.765(4.882) | 1.096(0.307) | 1.060(0.301) | 0.592(0.145) | 0.606(0.153)
Huber | 10.115(2.084) | 0.852(0.182) | 0.820(0.175) | 0.486(0.112) | 0.510(0.125)
Cauchy | 24.557(4.408) | 1.254(0.199) | 1.233(0.198) | 0.469(0.068) | 0.375(0.052)
Tukey | 164.275(27.099) | 6.745(0.605) | 6.692(0.606) | 2.027(0.107) | 1.297(0.117)

Table 3: The testing excess risks and standard deviation (in parentheses)
under different loss functions of the ERMs trained with respect to different
loss functions using Nets-256. Data is generated from ”Bumps” model with
training sample size $n=128$ or $512$ and the number of replications $R=10$.

Bumps, $n=128$ | Testing Loss
---|---
Error | Training Loss | LS | LAD | Huber | Cauchy | Tukey
$N(0,1)$ | LS | 5.991(1.141) | 0.798(0.114) | 0.762(0.111) | 0.493(0.063) | 0.541(0.065)
LAD | 6.859(1.065) | 0.821(0.091) | 0.785(0.088) | 0.487(0.057) | 0.520(0.070)
Huber | 6.214(1.289) | 0.799(0.116) | 0.763(0.112) | 0.487(0.065) | 0.530(0.071)
Cauchy | 8.011(0.631) | 0.898(0.043) | 0.862(0.039) | 0.515(0.041) | 0.538(0.057)
Tukey | 8.480(0.519) | 0.909(0.057) | 0.874(0.055) | 0.508(0.036) | 0.525(0.042)
$t(2)$ | LS | 8.697(2.977) | 1.018(0.223) | 0.974(0.218) | 0.593(0.103) | 0.629(0.096)
LAD | 7.491(0.719) | 0.869(0.086) | 0.829(0.081) | 0.507(0.062) | 0.535(0.077)
Huber | 7.543(0.924) | 0.866(0.080) | 0.827(0.078) | 0.498(0.051) | 0.518(0.064)
Cauchy | 8.743(0.582) | 0.901(0.028) | 0.865(0.028) | 0.492(0.022) | 0.497(0.033)
Tukey | 9.477(0.951) | 0.967(0.099) | 0.930(0.096) | 0.524(0.059) | 0.526(0.068)
$Cauchy(0,1)$ | LS | 774.661(1487.473) | 3.991(3.920) | 3.931(3.918) | 1.031(0.228) | 0.974(0.187)
LAD | 11.500(2.969) | 0.810(0.076) | 0.773(0.073) | 0.451(0.048) | 0.458(0.053)
Huber | 11.087(2.643) | 0.852(0.142) | 0.814(0.136) | 0.478(0.081) | 0.488(0.089)
Cauchy | 9.486(2.771) | 0.855(0.126) | 0.820(0.122) | 0.463(0.077) | 0.463(0.092)
Tukey | 11.173(3.318) | 0.889(0.118) | 0.853(0.116) | 0.473(0.064) | 0.463(0.065)
Mixture | LS | 432.363(307.985) | 10.093(4.175) | 10.000(4.171) | 2.845(0.623) | 2.106(0.278)
LAD | 13.203(9.829) | 1.064(0.378) | 1.027(0.374) | 0.569(0.142) | 0.584(0.123)
Huber | 21.442(31.043) | 1.232(0.627) | 1.194(0.625) | 0.601(0.135) | 0.598(0.096)
Cauchy | 8.017(1.069) | 0.821(0.109) | 0.788(0.106) | 0.467(0.060) | 0.489(0.064)
Tukey | 6.976(0.435) | 0.730(0.032) | 0.699(0.030) | 0.424(0.025) | 0.448(0.034)
Bumps, $n=512$ | LS | LAD | Huber | Cauchy | Tukey
$N(0,1)$ | LS | 4.248(1.836) | 0.602(0.183) | 0.573(0.176) | 0.382(0.102) | 0.425(0.110)
LAD | 4.377(0.748) | 0.570(0.061) | 0.544(0.059) | 0.350(0.030) | 0.381(0.029)
Huber | 3.810(0.725) | 0.530(0.070) | 0.505(0.068) | 0.332(0.040) | 0.366(0.043)
Cauchy | 5.605(0.738) | 0.641(0.050) | 0.616(0.050) | 0.368(0.022) | 0.385(0.023)
Tukey | 6.499(0.872) | 0.687(0.073) | 0.663(0.071) | 0.375(0.037) | 0.378(0.037)
$t(2)$ | LS | 4.192(1.204) | 0.606(0.134) | 0.573(0.129) | 0.386(0.076) | 0.424(0.081)
LAD | 4.704(0.773) | 0.572(0.062) | 0.546(0.060) | 0.334(0.030) | 0.349(0.031)
Huber | 4.231(1.020) | 0.531(0.097) | 0.504(0.094) | 0.318(0.053) | 0.336(0.056)
Cauchy | 6.140(0.614) | 0.650(0.051) | 0.624(0.049) | 0.355(0.030) | 0.356(0.036)
Tukey | 6.503(0.750) | 0.653(0.049) | 0.628(0.049) | 0.345(0.020) | 0.336(0.017)
$Cauchy(0,1)$ | LS | 281.707(421.853) | 3.520(2.033) | 3.457(2.027) | 1.136(0.269) | 1.053(0.221)
LAD | 5.259(1.686) | 0.493(0.054) | 0.468(0.052) | 0.287(0.032) | 0.296(0.037)
Huber | 5.266(1.538) | 0.493(0.093) | 0.467(0.091) | 0.293(0.049) | 0.306(0.050)
Cauchy | 6.621(0.998) | 0.575(0.051) | 0.552(0.050) | 0.301(0.026) | 0.290(0.027)
Tukey | 7.514(0.964) | 0.642(0.041) | 0.616(0.040) | 0.328(0.017) | 0.309(0.017)
Mixture | LS | 142.978(89.599) | 5.354(1.038) | 5.270(1.037) | 2.072(0.207) | 1.818(0.141)
LAD | 8.067(3.051) | 0.599(0.135) | 0.574(0.133) | 0.342(0.054) | 0.360(0.048)
Huber | 5.834(0.673) | 0.502(0.049) | 0.478(0.047) | 0.311(0.033) | 0.342(0.039)
Cauchy | 7.254(0.561) | 0.541(0.037) | 0.519(0.036) | 0.310(0.022) | 0.325(0.026)
Tukey | 8.339(1.042) | 0.585(0.054) | 0.564(0.053) | 0.317(0.024) | 0.321(0.024)

Table 4: The testing excess risks and standard deviation (in parentheses)
under different loss functions of the ERMs trained with respect to different
loss functions using Nets-256. Data is generated from ”Heavisine” model with
training sample size $n=128$ or $512$ and replications $R=10$.

Heavisine, $n=128$ | Testing Loss
---|---
Error | Training Loss | LS | LAD | Huber | Cauchy | Tukey
$N(0,1)$ | LS | 2.623(2.612) | 0.582(0.455) | 0.542(0.436) | 0.432(0.305) | 0.518(0.365)
LAD | 1.742(0.354) | 0.428(0.070) | 0.393(0.065) | 0.330(0.054) | 0.392(0.065)
Huber | 1.730(0.609) | 0.441(0.166) | 0.405(0.153) | 0.343(0.133) | 0.414(0.170)
Cauchy | 1.823(0.564) | 0.394(0.096) | 0.366(0.091) | 0.291(0.070) | 0.340(0.085)
Tukey | 288.203(541.200) | 9.012(9.163) | 8.934(9.144) | 2.522(1.385) | 1.747(0.722)
$t(2)$ | LS | 3.899(3.483) | 0.641(0.227) | 0.597(0.220) | 0.448(0.127) | 0.521(0.147)
LAD | 2.173(1.005) | 0.442(0.152) | 0.408(0.145) | 0.319(0.099) | 0.367(0.113)
Huber | 1.839(0.341) | 0.412(0.066) | 0.378(0.061) | 0.307(0.049) | 0.358(0.059)
Cauchy | 2.422(1.596) | 0.442(0.192) | 0.410(0.185) | 0.308(0.111) | 0.347(0.117)
Tukey | 9102.458(26252.125) | 38.803(84.576) | 38.721(84.566) | 3.504(2.490) | 1.857(0.550)
$Cauchy(0,1)$ | LS | 1090.375(3103.706) | 4.247(5.722) | 4.184(5.715) | 1.156(0.562) | 1.049(0.340)
LAD | 3.975(5.956) | 0.214(0.051) | 0.193(0.047) | 0.153(0.036) | 0.170(0.040)
Huber | 4.150(9.158) | 0.332(0.102) | 0.302(0.096) | 0.244(0.074) | 0.277(0.089)
Cauchy | 0.650(4.657) | 0.224(0.038) | 0.203(0.035) | 0.160(0.027) | 0.177(0.031)
Tukey | 452.600(1232.173) | 8.023(19.907) | 7.989(19.885) | 1.139(2.016) | 0.576(0.682)
Mixture | LS | 349.663(180.878) | 8.992(2.705) | 8.905(2.701) | 2.652(0.427) | 2.015(0.199)
LAD | 3.748(1.782) | 0.609(0.239) | 0.571(0.230) | 0.436(0.153) | 0.514(0.179)
Huber | 18.776(25.589) | 0.888(0.527) | 0.851(0.525) | 0.461(0.129) | 0.496(0.109)
Cauchy | 3.719(3.786) | 0.493(0.283) | 0.464(0.277) | 0.332(0.132) | 0.377(0.123)
Tukey | 1659.819(2865.520) | 19.828(25.013) | 19.752(25.000) | 3.055(1.834) | 1.749(0.528)
Heavisine, $n=512$ | LS | LAD | Huber | Cauchy | Tukey
$N(0,1)$ | LS | 1.097(0.194) | 0.287(0.060) | 0.263(0.054) | 0.225(0.050) | 0.269(0.060)
LAD | 1.063(0.275) | 0.274(0.079) | 0.251(0.071) | 0.214(0.064) | 0.255(0.075)
Huber | 1.492(0.729) | 0.379(0.182) | 0.346(0.166) | 0.296(0.144) | 0.352(0.172)
Cauchy | 1.108(0.299) | 0.280(0.089) | 0.256(0.080) | 0.217(0.073) | 0.256(0.087)
Tukey | 24153.281(72343.156) | 47.157(135.209) | 47.110(135.182) | 1.862(3.286) | 0.909(0.918)
$t(2)$ | LS | 4.782(9.469) | 0.451(0.300) | 0.419(0.294) | 0.290(0.126) | 0.325(0.136)
LAD | 1.456(0.464) | 0.320(0.107) | 0.291(0.098) | 0.237(0.082) | 0.272(0.097)
Huber | 1.434(0.458) | 0.328(0.117) | 0.298(0.106) | 0.245(0.091) | 0.283(0.106)
Cauchy | 1.356(0.271) | 0.291(0.048) | 0.264(0.044) | 0.213(0.034) | 0.239(0.038)
Tukey | 4.210(7.350) | 0.434(0.470) | 0.411(0.461) | 0.248(0.205) | 0.254(0.193)
$Cauchy(0,1)$ | LS | 1090.375(3103.706) | 4.247(5.722) | 4.184(5.715) | 1.156(0.562) | 1.049(0.340)
LAD | 3.975(5.956) | 0.214(0.051) | 0.193(0.047) | 0.153(0.036) | 0.170(0.040)
Huber | 4.150(9.158) | 0.332(0.102) | 0.302(0.096) | 0.244(0.074) | 0.277(0.089)
Cauchy | 0.650(4.657) | 0.224(0.038) | 0.203(0.035) | 0.160(0.027) | 0.177(0.031)
Tukey | 452.600(1232.173) | 8.023(19.907) | 7.989(19.885) | 1.139(2.016) | 0.576(0.682)
Mixture | LS | 100.342(40.824) | 5.224(0.974) | 5.137(0.972) | 2.182(0.229) | 1.928(0.155))
LAD | 0.864(0.483) | 0.297(0.090) | 0.271(0.082) | 0.235(0.073) | 0.281(0.087)
Huber | 0.885(0.371) | 0.267(0.072) | 0.244(0.065) | 0.207(0.059) | 0.246(0.070)
Cauchy | 0.860(0.291) | 0.235(0.038) | 0.214(0.034) | 0.182(0.032) | 0.213(0.037)
Tukey | 51982.227(155652.234) | 70.234(200.608) | 70.193(200.588) | 1.984(3.050) | 0.878(0.689)

Table 5: The testing excess risks and standard deviation (in parentheses)
under different loss functions of the ERMs trained with respect to different
loss functions using Nets-256. Data is generated from ”Dopller” model with
training sample size $n=128$ or $512$ and the number of replications $R=10$.

Dopller, $n=128$ | Testing Loss
---|---
Error | Training Loss | LS | LAD | Huber | Cauchy | Tukey
$N(0,1)$ | LS | 14.256(6.635) | 1.385(0.359) | 1.338(0.358) | 0.731(0.130) | 0.724(0.118)
LAD | 30.378(8.560) | 2.192(0.583) | 2.139(0.575) | 0.983(0.242) | 0.895(0.202)
Huber | 18.495(6.667) | 1.630(0.378) | 1.580(0.373) | 0.821(0.149) | 0.795(0.133)
Cauchy | 67.333(33.892) | 3.578(1.406) | 3.519(1.398) | 1.328(0.411) | 1.088(0.270)
Tukey | 4923.002(13946.816) | 30.952(60.211) | 30.857(60.199) | 3.884(2.361) | 2.253(0.653)
$t(2)$ | LS | 18.765(9.483) | 1.711(0.563) | 1.657(0.557) | 0.868(0.226) | 0.843(0.212)
LAD | 21.285(4.814) | 1.762(0.239) | 1.709(0.237) | 0.857(0.109) | 0.815(0.117)
Huber | 31.808(11.372) | 2.304(0.588) | 2.247(0.581) | 1.031(0.223) | 0.937(0.204)
Cauchy | 78.870(25.555) | 3.930(0.650) | 3.875(0.649) | 1.357(0.140) | 1.034(0.110)
Tukey | 180.530(36.958) | 8.965(1.507) | 8.883(1.497) | 2.837(0.495) | 1.904(0.339)
$Cauchy(0,1)$ | LS | 146.669(177.462) | 4.229(2.136) | 4.161(2.129) | 1.463(0.405) | 1.238(0.248)
LAD | 22.448(8.524) | 2.279(0.438) | 2.225(0.431) | 0.978(0.164) | 0.859(0.155)
Huber | 20.960(15.823) | 2.005(0.780) | 1.956(0.778) | 0.851(0.192) | 0.744(0.122)
Cauchy | 78.082(54.890) | 4.079(1.420) | 4.023(1.417) | 1.344(0.289) | 0.981(0.139)
Tukey | 197.141(64.730) | 8.053(1.696) | 7.982(1.690) | 2.375(0.419) | 1.531(0.242)
Mixture | LS | 450.810(193.680) | 10.148(1.856) | 10.056(1.856) | 2.830(0.259) | 2.113(0.128)
LAD | 51.893(70.423) | 2.092(0.928) | 2.042(0.927) | 0.894(0.164) | 0.834(0.115)
Huber | 31.130(11.430) | 1.836(0.474) | 1.791(0.471) | 0.837(0.166) | 0.776(0.128)
Cauchy | 79.468(26.835) | 3.426(0.508) | 3.377(0.507) | 1.214(0.123) | 0.956(0.109)
Tukey | 7063.935(20676.158) | 25.706(55.641) | 25.634(55.634) | 2.770(1.459) | 1.684(0.315)
Dopller, $n=512$ | LS | LAD | Huber | Cauchy | Tukey
$N(0,1)$ | LS | 13.193(5.793) | 1.267(0.369) | 1.226(0.364) | 0.654(0.149) | 0.633(0.128)
LAD | 17.459(9.004) | 1.458(0.529) | 1.414(0.524) | 0.715(0.191) | 0.682(0.154)
Huber | 19.160(5.569) | 1.558(0.389) | 1.516(0.383) | 0.743(0.184) | 0.699(0.183)
Cauchy | 35.604(17.747) | 2.243(0.649) | 2.198(0.645) | 0.922(0.157) | 0.791(0.097)
Tukey | 137.879(102.632) | 5.814(2.568) | 5.756(2.558) | 1.787(0.492) | 1.267(0.265)
$t(2)$ | LS | 13.319(5.055) | 1.374(0.276) | 1.320(0.274) | 0.743(0.120) | 0.738(0.146)
LAD | 16.387(7.046) | 1.302(0.367) | 1.266(0.362) | 0.606(0.147) | 0.551(0.133)
Huber | 20.022(8.077) | 1.558(0.477) | 1.513(0.468) | 0.720(0.195) | 0.651(0.176)
Cauchy | 28.573(11.641) | 1.847(0.498) | 1.808(0.496) | 0.753(0.165) | 0.630(0.137)
Tukey | 537.837(1257.923) | 11.205(16.875) | 11.143(16.858) | 2.317(1.715) | 1.366(0.543)
$Cauchy(0,1)$ | LS | 463.967(861.745) | 5.217(4.085) | 5.150(4.076) | 1.484(0.552) | 1.238(0.345)
LAD | 12.482(6.446) | 1.015(0.327) | 0.985(0.326) | 0.473(0.101) | 0.424(0.069)
Huber | 16.607(5.793) | 1.214(0.307) | 1.180(0.305) | 0.540(0.094) | 0.471(0.070)
Cauchy | 29.285(16.817) | 1.800(0.620) | 1.767(0.620) | 0.671(0.142) | 0.525(0.081)
Tukey | 133.563(95.067) | 5.897(3.034) | 5.843(3.021) | 1.728(0.705) | 1.123(0.376)
Mixture | LS | 127.257(44.378) | 6.153(1.271) | 6.063(1.269) | 2.394(0.303) | 2.014(0.150)
LAD | 20.396(9.476) | 1.420(0.446) | 1.382(0.443) | 0.668(0.156) | 0.619(0.123)
Huber | 14.738(4.872) | 1.149(0.248) | 1.114(0.245) | 0.582(0.088) | 0.571(0.068)
Cauchy | 53.697(28.314) | 2.410(0.814) | 2.373(0.811) | 0.873(0.201) | 0.697(0.125)
Tukey | 137.383(83.918) | 5.270(1.839) | 5.215(1.831) | 1.679(0.452) | 1.221(0.338)

Figure 4: The fitted curves of univariate models. The training data is
depicted as grey dots, the target function is depicted as black curves. The
colored curves represent the estimated curves under different loss functions.
Red line: LS; Orange line: LAD; Green line: Huber; Purple line: Cauchy; Blue
line: Tukey. From the top to the bottom, each row corresponds a certain type
of error: $N(0,1),t(2)$, $Cauchy(0,1)$ and normal mixture. From the left to
right, each column corresponds a certain univariate model: “Blocks”, “Bumps”,
“Heavisine” and “Dopller”. Table 6: The testing excess risks and standard
deviation (in parentheses) under different loss functions of the ERMs trained
with respect to different loss functions using Nets-256. The dimension of the
input $d=4$, the training sample size $n=128,512$ and the number of
replications $R=10$.

$d=4,n=128$ | Testing Loss
---|---
Error | Training Loss | LS | LAD | Huber | Cauchy | Tukey
$N(0,1)$ | LS | 13.678(0.674) | 1.336(0.046) | 1.272(0.045) | 0.814(0.027) | 0.923(0.031)
LAD | 14.832(1.483) | 1.368(0.119) | 1.304(0.114) | 0.821(0.079) | 0.926(0.096)
Huber | 14.033(0.463) | 1.325(0.066) | 1.262(0.064) | 0.799(0.041) | 0.903(0.047)
Cauchy | 22.141(2.693) | 1.505(0.111) | 1.441(0.108) | 0.834(0.065) | 0.923(0.079)
Tukey | 59.162(6.536) | 2.340(0.185) | 2.275(0.183) | 0.980(0.062) | 0.977(0.058)
$t(2)$ | LS | 19.227(3.166) | 1.725(0.242) | 1.660(0.238) | 0.961(0.112) | 1.029(0.109)
LAD | 18.103(2.433) | 1.483(0.149) | 1.422(0.146) | 0.819(0.074) | 0.879(0.076)
Huber | 20.011(1.943) | 1.684(0.114) | 1.620(0.112) | 0.922(0.059) | 0.981(0.062)
Cauchy | 21.666(1.591) | 1.570(0.136) | 1.510(0.134) | 0.842(0.072) | 0.897(0.074)
Tukey | 65.676(14.472) | 2.509(0.381) | 2.446(0.379) | 0.998(0.103) | 0.955(0.077)
$Cauchy(0,1)$ | LS | 2837.600(5740.160) | 6.181(5.959) | 6.116(5.955) | 1.386(0.380) | 1.223(0.223)
LAD | 18.400(113.676) | 2.098(0.963) | 2.041(0.962) | 0.947(0.150) | 0.942(0.096)
Huber | 130.800(422.469) | 2.623(1.663) | 2.563(1.663) | 1.033(0.160) | 1.012(0.112)
Cauchy | 32.026(8.921) | 1.656(0.283) | 1.601(0.280) | 0.787(0.097) | 0.802(0.089)
Tukey | 44.800(43.074) | 2.541(0.236) | 2.485(0.236) | 0.971(0.052) | 0.908(0.033)
$d=4,n=512$ | LS | LAD | Huber | Cauchy | Tukey
$N(0,1)$ | LS | 5.096(0.385) | 0.678(0.040) | 0.634(0.037) | 0.459(0.027) | 0.527(0.033)
LAD | 5.821(0.561) | 0.729(0.056) | 0.683(0.052) | 0.484(0.042) | 0.552(0.052)
Huber | 5.172(0.449) | 0.678(0.025) | 0.634(0.025) | 0.457(0.015) | 0.525(0.018)
Cauchy | 8.733(0.820) | 0.834(0.037) | 0.787(0.035) | 0.517(0.035) | 0.582(0.045)
Tukey | 23.338(5.598) | 1.166(0.134) | 1.120(0.134) | 0.561(0.031) | 0.581(0.031)
$t(2)$ | LS | 11.257(8.160) | 1.028(0.149) | 0.976(0.147) | 0.624(0.050) | 0.696(0.054)
LAD | 7.113(1.641) | 0.753(0.104) | 0.710(0.101) | 0.464(0.057) | 0.511(0.064)
Huber | 5.934(0.493) | 0.764(0.077) | 0.718(0.073) | 0.496(0.055) | 0.557(0.068)
Cauchy | 9.490(0.809) | 0.820(0.026) | 0.777(0.025) | 0.477(0.018) | 0.515(0.024)
Tukey | 25.927(5.004) | 1.273(0.145) | 1.227(0.143) | 0.590(0.051) | 0.594(0.054)
$Cauchy(0,1)$ | LS | 3705.178(7424.642) | 7.270(5.723) | 7.201(5.717) | 1.567(0.767) | 1.292(0.355)
LAD | 9.489(9.208) | 0.917(0.128) | 0.872(0.127) | 0.520(0.060) | 0.558(0.069)
Huber | 10.836(5.023) | 0.966(0.101) | 0.919(0.100) | 0.552(0.047) | 0.590(0.047)
Cauchy | 11.392(6.910) | 0.873(0.113) | 0.830(0.112) | 0.476(0.049) | 0.501(0.058)
Tukey | 34.442(6.141) | 1.445(0.152) | 1.399(0.149) | 0.615(0.064) | 0.598(0.069)

##  7 Concluding remarks

In this work, we establish the non-asymptotic error bounds for a class of
robust nonparametric regression estimators using deep neural networks with
ReLU activation under suitable smoothness conditions on the regression
function and mild conditions on the error term. In particular, we only assume
that the error distribution has a finite $p$-th moment with $p>1$. We also
show that the deep robust regression is able to circumvent the curse of
dimensionality when the predictor is supported on an approximate lower-
dimensional set. An important feature of our error bound is that, for ReLU
neural networks with network width and network size (number of parameters) no
more than the order of the square of the dimensionality $d$ of the predictor,
our excess risk bounds depend sub-linearly on the $d$. Our assumption relaxes
the exact manifold support assumption, which could be restrictive and
unrealistic in practice. We also relax several crucial assumptions on the data
distribution, the target regression function and the neural networks required
in the recent literature. Our simulation studies demonstrate the advantages of
using robust loss functions over the least squares loss function in the
context of deep nonparametric regression in terms of testing excess risks.

We only addressed the issue of robustness with respect to the heavy tailed
errors or outliers in the response variable. There are some other important
sources of non- robustness in nonparametric regression. For example, it is
possible that outliers exist in the observed values of the predictors, or that
there are contaminations in both the responses and the predictors. How to
address these important practical issues in the context of nonparametric
regression using deep neural networks deserve further studies in the future.

## Acknowledgements

The work of Y. Jiao is supported in part by the National Science Foundation of
China grant 11871474 and by the research fund of KLATASDSMOE of China. The
work of Y. Lin is supported by the Hong Kong Research Grants Council (Grant
No. 14306219 and 14306620) and Direct Grants for Research, The Chinese
University of Hong Kong. The work of J. Huang is partially supported by the
U.S. NSF grant DMS-1916199.

## References

  * Anscombe (1960) [author] Anscombe, F.F. (1960). Rejection of outliers. Technometrics 2 123-146.
  * Anthony and Bartlett (1999) [author] Anthony, MartinM. and Bartlett, Peter L.P. L. (1999). Neural Network Learning: Theoretical Foundations. Cambridge University Press, Cambridge. 10.1017/CBO9780511624216
  * Baraniuk and Wakin (2009) [author] Baraniuk, Richard G.R. G. and Wakin, Michael B.M. B. (2009). Random projections of smooth manifolds. Found. Comput. Math. 9 51–77. 10.1007/s10208-007-9011-z
  * Barron (2019) [author] Barron, Jonathan TJ. T. (2019). A general and adaptive robust loss function. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition 4331–4339.
  * Bartlett et al. (2019) [author] Bartlett, Peter L.P. L., Harvey, NickN., Liaw, ChristopherC. and Mehrabian, AbbasA. (2019). Nearly-tight VC-dimension and pseudodimension bounds for piecewise linear neural networks. J. Mach. Learn. Res. 20 Paper No. 63, 17.
  * Bassett and Koenker (1978) [author] Bassett, GilbertG. and Koenker, RogerR. (1978). Asymptotic theory of least absolute error regression. Journal of the American Statistical Association 73 618-622. 10.1080/01621459.1978.10480065
  * Bauer and Kohler (2019) [author] Bauer, BenediktB. and Kohler, MichaelM. (2019). On deep learning as a remedy for the curse of dimensionality in nonparametric regression. Annals of Statistics 47 2261–2285.
  * Beaton and Tukey (1974) [author] Beaton, Albert EA. E. and Tukey, John WJ. W. (1974). The fitting of power series, meaning polynomials, illustrated on band-spectroscopic data. Technometrics 16 147–185.
  * Belagiannis et al. (2015) [author] Belagiannis, VasileiosV., Rupprecht, ChristianC., Carneiro, GustavoG. and Navab, NassirN. (2015). Robust optimization for deep regression. In Proceedings of the IEEE international conference on computer vision 2830–2838.
  * Breiman et al. (1984) [author] Breiman, L.L., Friedman, J. H.J. H., Olshen, R. A.R. A. and Stone, C. J.C. J. (1984). Classification and Regression Trees. Wadsworth and Brooks, Monterey, CA.
  * Chen, Jiang and Zhao (2019) [author] Chen, MinshuoM., Jiang, HaomingH. and Zhao, TuoT. (2019). Efficient approximation of deep relu networks for functions on low dimensional manifolds. Advances in Neural Information Processing Systems.
  * Chen et al. (2019) [author] Chen, XirongX., Li, DeguiD., Li, QiQ. and Li, ZhengZ. (2019). Nonparametric estimation of conditional quantile functions in the presence of irrelevant covariates. Journal of Econometrics 212 433–450.
  * Cherapanamjeri et al. (2020) [author] Cherapanamjeri, YeshwanthY., Aras, EfeE., Tripuraneni, NileshN., Jordan, Michael I.M. I., Flammarion, NicolasN. and Bartlett, Peter L.P. L. (2020). Optimal Robust Linear Regression in Nearly Linear Time.
  * Christmann and Steinwart (2007) [author] Christmann, AndreasA. and Steinwart, IngoI. (2007). How SVMs can estimate quantiles and the median. In Advances in neural information processing systems 305–312.
  * Cleveland (1979) [author] Cleveland, William S.W. S. (1979). Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association 74 829-836. 10.1080/01621459.1979.10481038
  * Devore and Ron (2010) [author] Devore, RonaldR. and Ron, AmosA. (2010). Approximation using scattered shifts of a multivariate function. Trans. Amer. Math. Soc. 362 6205–6229. 10.1090/S0002-9947-2010-05070-6
  * Donoho and Johnstone (1994) [author] Donoho, David LD. L. and Johnstone, Iain MI. M. (1994). Ideal spatial adaptation by wavelet shrinkage. Biometrika 81 425-455. 10.1093/biomet/81.3.425
  * Fan (1992) [author] Fan, JianqingJ. (1992). Design-adaptive nonparametric regression. Journal of the American Statistical Association 87 998-1004. 10.1080/01621459.1992.10476255
  * Fan (1993) [author] Fan, JianqingJ. (1993). Local linear regression smoothers and their minimax efficiencies. Annals of Statistics 21 196 – 216. 10.1214/aos/1176349022
  * Farrell, Liang and Misra (2021) [author] Farrell, Max H.M. H., Liang, TengyuanT. and Misra, SanjogS. (2021). Deep neural networks for estimation and inference. Econometrica 89 181–213. 10.3982/ecta16901
  * Fridedman (1991) [author] Fridedman, JHJ. (1991). Multivariate adaptive regression splines (with discussion). Annals of Statistics 19 79–141.
  * Ghorbani et al. (2020) [author] Ghorbani, BehroozB., Mei, SongS., Misiakiewicz, TheodorT. and Montanari, AndreaA. (2020). Discussion of:“Nonparametric regression using deep neural networks with ReLU activation function”. The Annals of Statistics 48 1898–1901.
  * Green and Silverman (1993) [author] Green, P. J.P. J. and Silverman, B. W.B. W. (1993). Nonparametric Regression and Generalized Linear Models: A roughness penalty approach. Chapman and Hall/CRC.
  * Györfi et al. (2002) [author] Györfi, LászlóL., Kohler, MichaelM., Krzy˙zak, AdamA. and Walk, HarroH. (2002). A Distribution-free Theory of Nonparametric Regression. Springer Series in Statistics. Springer-Verlag, New York. 10.1007/b97848
  * Hampel et al. (2011) [author] Hampel, Frank RF. R., Ronchetti, Elvezio ME. M., Rousseeuw, Peter JP. J. and Stahel, Werner AW. A. (2011). Robust statistics: the approach based on influence functions 196\. John Wiley & Sons.
  * Han and Wellner (2019) [author] Han, QiyangQ. and Wellner, Jon A.J. A. (2019). Convergence rates of least squares regression estimators with heavy-tailed errors. Annals of Statistics 47 2286–2319.
  * Hangelbroek and Ron (2010) [author] Hangelbroek, ThomasT. and Ron, AmosA. (2010). Nonlinear approximation using Gaussian kernels. J. Funct. Anal. 259 203–219. 10.1016/j.jfa.2010.02.001
  * Hastie, Tibshirani and Friedman (2001) [author] Hastie, TrevorT., Tibshirani, RobertR. and Friedman, JeromeJ. (2001). The Elements of Statistical Learning. Springer Series in Statistics. Springer New York Inc., New York, NY, USA.
  * Hecht-Nielsen (1987) [author] Hecht-Nielsen, RobertR. (1987). Kolmogorov’s mapping neural network existence theorem. In Proceedings of the international conference on Neural Networks 3 11–14. IEEE Press New York.
  * Huber (1964) [author] Huber, Peter J.P. J. (1964). Robust estimation of a location parameter. Annals of Mathematical Statistics 35 73–101. 10.1214/aoms/1177703732
  * Huber (1973) [author] Huber, Peter JP. J. (1973). Robust regression: asymptotics, conjectures and Monte Carlo. Annals of statistics 1 799–821.
  * Huber (2004) [author] Huber, Peter JP. J. (2004). Robust Statistics 523\. John Wiley & Sons.
  * Jambulapati et al. (2021) [author] Jambulapati, ArunA., Li, JerryJ., Schramm, TselilT. and Tian, KevinK. (2021). Robust Regression Revisited: Acceleration and Improved Estimation Rates.
  * Jiang et al. (2018) [author] Jiang, LuL., Zhou, ZhengyuanZ., Leung, ThomasT., Li, Li-JiaL.-J. and Fei-Fei, LiL. (2018). Mentornet: Learning data-driven curriculum for very deep neural networks on corrupted labels. In International Conference on Machine Learning 2304–2313. PMLR.
  * Jiao et al. (2021) [author] Jiao, YulingY., Shen, GuohaoG., Lin, YuanyuanY. and Huang, JianJ. (2021). Deep nonparametric regression on approximately low-dimensional manifolds. arXiv preprint arXiv:2104.06708.
  * Kingma and Ba (2014) [author] Kingma, Diederik PD. P. and Ba, JimmyJ. (2014). Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
  * LeCun, Cortes and Burges (2010) [author] LeCun, YannY., Cortes, CorinnaC. and Burges, CJC. (2010). MNIST handwritten digit database. AT&T Labs [Online]. Available: http://yann. lecun. com/exdb/mnist 2\.
  * Lederer (2020) [author] Lederer, JohannesJ. (2020). Risk bounds for robust deep learning. arXiv preprint arXiv:2009.06202.
  * Lee (2006) [author] Lee, John MJ. M. (2006). Riemannian Manifolds: An Introduction to Curvature 176\. Springer Science & Business Media.
  * Liang and Srikant (2016) [author] Liang, ShiyuS. and Srikant, RayadurgamR. (2016). Why deep neural networks for function approximation? arXiv preprint arXiv:1610.04161.
  * Lin et al. (2014) [author] Lin, ShaoboS., Liu, XiaX., Rong, YuanhuaY. and Xu, ZongbenZ. (2014). Almost optimal estimates for approximation and learning by radial basis function networks. Mach. Learn. 95 147–164. 10.1007/s10994-013-5406-z
  * Lu et al. (2017) [author] Lu, ZhouZ., Pu, HongmingH., Wang, FeichengF., Hu, ZhiqiangZ. and Wang, LiweiL. (2017). The expressive power of neural networks: A view from the width. arXiv preprint arXiv:1709.02540.
  * Lv et al. (2018) [author] Lv, ShaogaoS., Lin, HuazhenH., Lian, HengH., Huang, JianJ. et al. (2018). Oracle inequalities for sparse additive quantile regression in reproducing kernel Hilbert space. The Annals of Statistics 46 781–813.
  * Mohri, Rostamizadeh and Talwalkar (2018) [author] Mohri, MehryarM., Rostamizadeh, AfshinA. and Talwalkar, AmeetA. (2018). Foundations of Machine Learning. Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA.
  * Nadaraya (1964) [author] Nadaraya, Elizbar AE. A. (1964). On estimating regression. Theory of Probability & Its Applications 9 141–142.
  * Nakada and Imaizumi (2019) [author] Nakada, RyumeiR. and Imaizumi, MasaakiM. (2019). Adaptive approximation and estimation of deep neural network with intrinsic dimensionality. arXiv preprint arXiv:1907.02177.
  * Padilla, Tansey and Chen (2020) [author] Padilla, Oscar Hernan MadridO. H. M., Tansey, WesleyW. and Chen, YanzhenY. (2020). Quantile regression with deep ReLU networks: Estimators and minimax rates. arXiv preprint arXiv:2010.08236v5.
  * Raghu et al. (2017) [author] Raghu, MaithraM., Poole, BenB., Kleinberg, JonJ., Ganguli, SuryaS. and Sohl-Dickstein, JaschaJ. (2017). On the expressive power of deep neural networks. In international conference on machine learning 2847–2854. PMLR.
  * Schmidt-Hieber (2019) [author] Schmidt-Hieber, JohannesJ. (2019). Deep relu network approximation of functions on a manifold. arXiv preprint arXiv:1908.00695.
  * Schmidt-Hieber (2020) [author] Schmidt-Hieber, JohannesJ. (2020). “Nonparametric regression using deep neural networks with ReLU activation function”. Annals of Statistics 48 1916–1921. 10.1214/19-AOS1931
  * Schmidt-Hieber (2021) [author] Schmidt-Hieber, JohannesJ. (2021). The Kolmogorov–Arnold representation theorem revisited. Neural Networks 137 119–126.
  * Scholkopf and Smola (2018) [author] Scholkopf, B.B. and Smola, A. J.A. J. (2018). Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT Press.
  * Shen, Yang and Zhang (2019) [author] Shen, ZuoweiZ., Yang, HaizhaoH. and Zhang, ShijunS. (2019). Nonlinear approximation via compositions. Neural Networks 119 74–84.
  * Shen, Yang and Zhang (2020) [author] Shen, ZuoweiZ., Yang, HaizhaoH. and Zhang, ShijunS. (2020). Deep network approximation characterized by number of neurons. Commun. Comput. Phys. 28 1768–1811. 10.4208/cicp.oa-2020-0149
  * Steinwart et al. (2011) [author] Steinwart, IngoI., Christmann, AndreasA. et al. (2011). Estimating conditional quantiles with the help of the pinball loss. Bernoulli 17 211–225.
  * Stone (1982) [author] Stone, Charles J.C. J. (1982). Optimal global rates of convergence for nonparametric regression. Annals of Statistics 10 1040–1053.
  * Suzuki (2018) [author] Suzuki, TaijiT. (2018). Adaptivity of deep ReLU network for learning in Besov and mixed smooth Besov spaces: optimal rate and curse of dimensionality. arXiv preprint arXiv:1810.08033.
  * Telgarsky (2016) [author] Telgarsky, MatusM. (2016). Benefits of depth in neural networks. In Conference on learning theory 1517–1539. PMLR.
  * Tsybakov (2008) [author] Tsybakov, AlexandreA. (2008). Introduction to Nonparametric Estimation. Springer Science & Business Media.
  * Tukey (1960) [author] Tukey, J. W.J. W. (1960). A survey of sampling from contaminated distributions. Contributions to Probability and Statistics 448-485.
  * Tukey (1975) [author] Tukey, J. W.J. W. (1975). Mathematics and the picturing of data. Proceedings of the International Congress of Mathematicians, Vancouver, 1975 2 523-531.
  * van der Vaart and Wellner (1996) [author] van der Vaart, Aad W.A. W. and Wellner, Jon A.J. A. (1996). Weak Convergence and Empirical Processes: With Applications to Statistics. Springer Series in Statistics. Springer-Verlag, New York. 10.1007/978-1-4757-2545-2
  * Wahba (1990) [author] Wahba, GraceG. (1990). Spline Models for Observational Data. Society for Industrial and Applied Mathematics. 10.1137/1.9781611970128
  * Wang et al. (2016) [author] Wang, ZhangyangZ., Chang, ShiyuS., Yang, YingzhenY., Liu, DingD. and Huang, Thomas ST. S. (2016). Studying very low resolution recognition using deep networks. In Proceedings of the IEEE conference on computer vision and pattern recognition 4792–4800.
  * Watson (1964) [author] Watson, Geoffrey SG. S. (1964). Smooth regression analysis. Sankhyā: The Indian Journal of Statistics, Series A 359–372.
  * Yarotsky (2017) [author] Yarotsky, DmitryD. (2017). Error bounds for approximations with deep ReLU networks. Neural Networks 94 103–114.
  * Yarotsky (2018) [author] Yarotsky, DmitryD. (2018). Optimal approximation of continuous functions by very deep ReLU networks. In Conference on Learning Theory 639–649. PMLR.

##  Appendix A A brief review of some recent results on neural network
approximation

In this paper, we only assume that $f^{*}$ is uniformly continuous (Assumption
3), which is a mild assumption on the continuity of the unknown target
function $f^{*}$, as the existing works generally assume stronger smoothness
assumptions on $f^{*}$. For example, in nonparametric least squares
regression, Stone (1982) and Bauer and Kohler (2019) require that $f^{*}$ is
$\beta$-Hölder smooth with $\beta\geq 1$, i.e., all its partial derivatives up
to order $\lfloor\beta\rfloor$ exist and the partial derivatives of order
$\lfloor\beta\rfloor$ are $\lfloor\beta\rfloor-\beta$ Hölder continuous.
Farrell, Liang and Misra (2021) assume that $f^{*}$ lies in a Sobolev ball
with smoothness $\beta\in\mathbb{N}^{+}$, i.e.
$f^{*}(x)\in\mathcal{W}^{\beta,\infty}([-1,1]^{d})$. Recent approximation
theories on different functional spaces, such as Korobov spaces, Besev spaces
or function space with $f^{*}\in C^{\alpha}[0,1]^{d}$ with $\alpha\geq 1$ can
be found in Liang and Srikant (2016), Lu et al. (2017), Yarotsky (2017),
Mohri, Rostamizadeh and Talwalkar (2018), Suzuki (2018), among others.

The function class $\mathcal{F}_{\phi}$ in this paper is constructed by the
feedforward neural networks with the ReLU activation function. Some recent
results on the approximation theory of deep neural networks can be found in
Devore and Ron (2010); Hangelbroek and Ron (2010); Lin et al. (2014); Yarotsky
(2017, 2018); Lu et al. (2017); Raghu et al. (2017); Shen, Yang and Zhang
(2019); Shen, Yang and Zhang (2020); Nakada and Imaizumi (2019); Chen, Jiang
and Zhao (2019).

An important result proved by Yarotsky (2017) is the following, for any
$\varepsilon\in(0,1)$, any $d,\beta$, and any $f_{0}$ in the Sobolev ball
$\mathcal{W}^{\beta,\infty}([0,1]^{d})$ with $\beta>0$, there exists a ReLU
network $\hat{f}$ with depth $\mathcal{D}$ at most
$c\\{\log(1/\varepsilon)+1\\}$, size $\mathcal{S}$ and number of neurons
$\mathcal{U}$ at most $c\varepsilon^{-d/\beta}\\{\log(1/\varepsilon)+1\\}$
such that
$\|\hat{f}-f_{0}\|_{\infty}=\max_{x\in[0,1]^{d}}|\hat{f}(x)-f_{0}(x)|\leq\varepsilon$,
where $c$ is some constant depending on $d$ and $\beta$. In particular, it is
required that the constant $c=O(2^{d})$, an exponential rate of $d$, due to
the technicality in the proof. The main idea of Yarotsky (2017) is to show
that, small neural networks can approximate polynomials locally, and stacked
neural networks (by $2^{d}$ small sub-networks) can further approximate smooth
function by approximating its Taylor expansions. Later, Yarotsky (2018) gave
the optimal rate of approximation for general continuous functions by deep
ReLU networks, in terms of the network size $\mathcal{S}$ and the modulus of
continuity of $f^{*}$. It was shown that
$\inf_{f\in\mathcal{F}_{\phi}}\|f-f^{*}\|_{\infty}\leq
c_{1}\omega_{f^{*}}(c_{2}\mathcal{S}^{-p/d})$ for some $p\in[1,2]$ and some
constants $c_{1},c_{2}$ possibly depending on $d,p$ but not
$\mathcal{S},f^{*}$. The upper bound holds for any $p\in(1,2]$ if the network
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
satisfies $\mathcal{D}\geq c_{3}\mathcal{S}^{p-1}/\log(\mathcal{S})$ for some
constant $c_{3}$ possibly depending on $p$ and $d$.

Based on a method similar to that in Yarotsky (2017), approximation of a
target function on a low-dimensional manifold $\mathcal{M}$ of $[0,1]^{d}$
with deep neural nets was studied by Nakada and Imaizumi (2019); Chen, Jiang
and Zhao (2019). The approximation error of ReLU networks to $f^{*}$ defined
on an exact manifold $\mathcal{M}$ will no longer depend on $d$, instead it
depends on the intrinsic dimension $d_{0}<d$ of the manifold $\mathcal{M}$.
Specifically, the approximation rate is shown to be
$c\mathcal{S}^{-\beta/d_{0}}$, where $f^{*}$ is often assumed to be a
$\beta$-Hölder smooth function defined on $[0,1]^{d}$ and $c$ is some constant
depending on $d,d_{0}$ and $\beta$. This helps alleviate the curse of
dimensionality in approximating high-dimensional functions with deep ReLU
neural networks.

It is worth noting that other related works, including Chen et al. (2019),
Chen, Jiang and Zhao (2019), Nakada and Imaizumi (2019), Schmidt-Hieber (2019)
and Schmidt-Hieber (2020), rely on a similar approximation construction of
Yarotsky (2017). A common feature of the results from these works is that, the
prefactor of the approximation error is of the order $O(2^{d})$ unless the
size $\mathcal{S}$ of the network has to grow exponentially with respect to
the dimension $d$. Unfortunately, a prefactor of the order $O(2^{d})$ is
extremely big for a large $d$ in the high-dimensional settings, which can
destroy the approximation error bound even for networks with a large size. For
example, each handwritten digit picture in the MNIST dataset (LeCun, Cortes
and Burges, 2010) is of dimension $d=28\times 28$ and the sample size $n$ is
about $70,000$, but $2^{d}$ is approximately $10^{236}$.

Lastly, it follows from Proposition 1 of Yarotsky (2017) that, for a neural
network, in terms of its computational power and complexity, there is no
substantial difference in using the ReLU activation function and other piece-
wise linear activation function with finitely many breakpoints. To elaborate,
let $\xi:\mathbb{R}\to\mathbb{R}$ be any continuous piece-wise linear function
with $M$ breakpoints ($1\leq M<\infty$). If a network $f_{\xi}$ is activated
by $\xi$, of depth $\mathcal{D}$, size $\mathcal{S}$ and the number of neurons
$\mathcal{U}$, then there exists a ReLU activated network with depth
$\mathcal{D}$, size not more than $(M+1)^{2}\mathcal{S}$, the number of
neurons not more than $(M+1)\mathcal{U}$, that computes the same function as
$f_{\xi}$. Conversely, let $f_{\sigma}$ be a ReLU activated network of depth
$\mathcal{D}$, size $\mathcal{S}$ and the number of neurons $\mathcal{U}$,
then there exists a network with activation function $\xi$, of depth
$\mathcal{D}$, size $4\mathcal{S}$ and the number of neurons $2\mathcal{U}$
that computes the same function $f_{\sigma}$ on a bounded subset of
$\mathbb{R}^{d}$.

##  Appendix B Proofs

In this part, we provide the proofs of Lemmas 3.1, 3.2 and 3.3, Theorems 4.2
and 4.3 and  Corollary 2.

###  B.1 Proof of Lemma 3.1

###### Proof.

Let $S=\\{Z_{i}=(X_{i},Y_{i})\\}_{i=1}^{n}$ be a sample form the distribution
of $Z=(X,Y)$ and
$S^{\prime}=\\{Z_{i}^{\prime}=(X^{\prime}_{i},Y^{\prime}_{i})\\}_{i=1}^{n}$ be
another sample independent of $S$. Define
$g(f,Z_{i})=L(f(X_{i}),Y_{i})-L(f^{*}(X_{i}),Y_{i})$ for any $f$ and sample
$Z_{i}$. Note that the empirical risk minimizer $\hat{f}_{\phi}$ defined in
Lemma 3.1 depends on the sample $S$, and its excess risk is
$\mathbb{E}_{S^{\prime}}\\{\sum_{i=1}^{n}g(\hat{f}_{\phi},Z_{i}^{\prime})/n\\}$
and its expected excess risk is

$\mathbb{E}\big{\\{}\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\big{\\}}=\mathbb{E}_{S}[\mathbb{E}_{S^{\prime}}\\{\frac{1}{n}\sum_{i=1}^{n}g(\hat{f}_{\phi},Z_{i}^{\prime})\\}].$
(B.1)

Next we will take 3 steps to complete the proof of Lemma 3.1.

#### Step 1: Excess risk decomposition

Define the ‘best in class’ estimator $f^{*}_{\phi}$ as the estimator in the
function class
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
with minimal $L$ risk:

$f^{*}_{\phi}=\arg\min_{f\in\mathcal{F}_{\phi}}\mathcal{R}(f).$

The approximation error of $f^{*}$ is
$\mathcal{R}(f^{*}_{\phi})-\mathcal{R}(f^{*})$. Note that the approximation
error only depends on the function class
$\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$
and the distribution of data. By the definition of empirical risk minimizer,
we have

$\mathbb{E}_{S}\\{\frac{1}{n}\sum_{i=1}^{n}g(\hat{f}_{\phi},Z_{i})\\}\leq\mathbb{E}_{S}\\{\frac{1}{n}\sum_{i=1}^{n}g(f^{*}_{\phi},Z_{i})\\}.$
(B.2)

Multiply 2 by the both sides of (B.2) and add it up with (B.1), we have

$\displaystyle\mathbb{E}\big{\\{}\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\big{\\}}$
$\displaystyle\leq\mathbb{E}_{S}\Big{[}\frac{1}{n}\sum_{i=1}^{n}\big{\\{}-2g(\hat{f}_{\phi},Z_{i})+\mathbb{E}_{S^{\prime}}g(\hat{f}_{\phi},Z_{i}^{\prime})\big{\\}}\Big{]}+2\mathbb{E}_{S}\\{\frac{1}{n}\sum_{i=1}^{n}g(f^{*}_{\phi},Z_{i})\\}$
$\displaystyle\leq\mathbb{E}_{S}\Big{[}\frac{1}{n}\sum_{i=1}^{n}\big{\\{}-2g(\hat{f}_{\phi},Z_{i})+\mathbb{E}_{S^{\prime}}g(\hat{f}_{\phi},Z_{i}^{\prime})\big{\\}}\Big{]}+2\big{\\{}\mathcal{R}(f^{*}_{\phi})-\mathcal{R}(f^{*})\big{\\}}.$
(B.3)

It is seen that the excess risk is upper bounded by the sum of a expectation
of a stochastic term and approximation error.

#### Step 2: Bounding the stochastic term

Next, we will focus on giving an upper bound of the first term on the right-
hand side in (B.3), and handle it with truncation and classical chaining
technique of empirical processes. In the following, for ease of presentation,
we write
$G(f,Z_{i})=\mathbb{E}_{S^{\prime}}\\{g(f,Z_{i}^{\prime})\\}-2g(f,Z_{i})$ for
$f\in\mathcal{F}_{\phi}$.

Given a $\delta$-uniform covering of $\mathcal{F}_{\phi}$, we denote the
centers of the balls by $f_{j},j=1,2,\ldots,\mathcal{N}_{2n},$ where
$\mathcal{N}_{2n}=\mathcal{N}_{2n}(\delta,\|\cdot\|_{\infty},\mathcal{F}_{\phi})$
is the uniform covering number with radius $\delta$ ($\delta<\mathcal{B}$)
under the norm $\|\cdot\|_{\infty}$, where
$\mathcal{N}_{2n}(\delta,\|\cdot\|_{\infty},\mathcal{F}_{\phi})$ is defined in
(3.3). By the definition of covering, there exists a (random) $j^{*}$ such
that $\|\hat{f}_{\phi}(x)-f_{j^{*}}(x)\|_{\infty}\leq\delta$ on
$x=(X_{1},\ldots,X_{n},X_{1}^{\prime},\ldots,X_{n}^{\prime})\in\mathcal{X}^{2n}$.
By the assumptions that
$\|f^{*}\|_{\infty},\|f_{j}\|_{\infty}\leq\mathcal{B}$,
$\mathbb{E}|Y_{i}|^{p}<\infty$, and Assumption 1, we have

$\displaystyle|g(\hat{f}_{\phi},Z_{i})-g(f_{j^{*}},Z_{i})|$
$\displaystyle=|L(\hat{f}_{\phi}(X_{i}),Y_{i})-L(f_{j^{*}}(X_{i}),Y_{i})|$
$\displaystyle\leq\lambda_{L}|\hat{f}_{\phi}(X_{i})-f_{j^{*}}(X_{i})|$
$\displaystyle\leq\lambda_{L}\delta,$

and

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S}\big{\\{}g(\hat{f}_{\phi},Z_{i})\\}\leq\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S}\big{\\{}$
$\displaystyle
g(f_{j^{*}},Z_{i})\\}+\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S}\big{\\{}g(\hat{f}_{\phi},Z_{i})-g(f_{j^{*}},Z_{i})\\}$
$\displaystyle\leq\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S}\big{\\{}$
$\displaystyle
g(f_{j^{*}},Z_{i})\\}+\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S}|g(\hat{f}_{\phi},Z_{i})-g(f_{j^{*}},Z_{i})|$
$\displaystyle\leq\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{S}\big{\\{}$
$\displaystyle g(f_{j^{*}},Z_{i})\\}+\lambda_{L}\delta.$

Then we have

$\displaystyle|G(\hat{f}_{\phi},Z_{i})-G(f_{j^{*}},Z_{i})|$
$\displaystyle=|\mathbb{E}_{S^{\prime}}\\{g(\hat{f}_{\phi},Z_{i}^{\prime})\\}-2g(\hat{f}_{\phi},Z_{i})-\mathbb{E}_{S^{\prime}}\\{g(f_{j^{*}},Z_{i}^{\prime})\\}+2g(f_{j^{*}},Z_{i})|$
$\displaystyle\leq|\mathbb{E}_{S^{\prime}}\\{g(\hat{f}_{\phi},Z_{i}^{\prime})\\}-\mathbb{E}_{S^{\prime}}\\{g(f_{j^{*}},Z_{i}^{\prime})\\}|+|2g(f_{j^{*}},Z_{i})-2g(\hat{f}_{\phi},Z_{i})|$
$\displaystyle\leq 3\lambda_{L}\delta,$

and

$\mathbb{E}_{S}\Big{[}\frac{1}{n}\sum_{i=1}^{n}G(\hat{f}_{\phi},Z_{i})\Big{]}\leq\mathbb{E}_{S}\Big{[}\frac{1}{n}\sum_{i=1}^{n}G(f_{j^{*}},Z_{i})\Big{]}+3\lambda_{L}\delta.$
(B.4)

Let $\beta_{n}\geq\mathcal{B}\geq 1$ be a positive number who may depend on
the sample size $n$. Denote $T_{\beta_{n}}$ as the truncation operator at
level $\beta_{n}$, i.e., for any $Y\in\mathbb{R}$, $T_{\beta_{n}}Y=Y$ if
$|Y|\leq\beta_{n}$ and $T_{\beta_{n}}Y=\beta_{n}\cdot{\rm sign}(Y)$ otherwise.
Define the function $f^{*}_{\beta_{n}}:\mathcal{X}\to\mathbb{R}$ pointwisely
by

$f^{*}_{\beta_{n}}(x)=\arg\min_{f(x):\|f\|_{\infty}\leq\beta_{n}}\mathbb{E}\big{\\{}L(f(X),T_{\beta_{n}}Y)|X=x\big{\\}},$

for each $x\in\mathcal{X}$. Besides, recall that
$\|f^{*}\|_{\infty}\leq\mathcal{B}\leq\beta_{n}$ and

$f^{*}(x)=\arg\min_{f(x):\|f\|_{\infty}\leq\beta_{n}}\mathbb{E}\big{\\{}L(f(X),Y)|X=x\big{\\}}.$

Then for any $f$ satisfying $\|f\|\infty\leq\beta_{n}$, the definition above
implies that
$\mathbb{E}\\{L(f^{*}_{\beta_{n}}(X_{i}),T_{\beta_{n}}Y_{i})\\}\leq\mathbb{E}\\{L(f(X_{i}),T_{\beta_{n}}Y_{i})\\}$
and
$\mathbb{E}\\{L(f^{*}(X_{i}),Y_{i})\\}\leq\mathbb{E}\\{L(f(X_{i}),Y_{i})\\}$.

For any $f\in\mathcal{F}_{\phi}$, we let
$g_{\beta_{n}}(f,Z_{i})=L(f(X_{i}),T_{\beta_{n}}Y_{i})-L(f^{*}_{\beta_{n}}(X_{i}),T_{\beta_{n}}Y_{i})$.
Then we have

$\displaystyle\mathbb{E}\\{g(f,Z_{i})\\}$
$\displaystyle=\mathbb{E}\\{g_{\beta_{n}}(f,Z_{i})\\}+\mathbb{E}\\{L(f(X_{i}),Y_{i})-L(f(X_{i}),T_{\beta_{n}}Y_{i})\\}$
$\displaystyle\qquad\qquad\qquad\quad+\mathbb{E}\\{L(f^{*}_{\beta_{n}}(X_{i}),T_{\beta_{n}}Y_{i})-L(f^{*}(X_{i}),T_{\beta_{n}}Y_{i})\\}$
$\displaystyle\qquad\qquad\qquad\quad+\mathbb{E}\\{L(f^{*}(X_{i}),T_{\beta_{n}}Y_{i})-L(f^{*}(X_{i}),Y_{i})\\}$
$\displaystyle\leq\mathbb{E}\\{g_{\beta_{n}}(f,Z_{i})\\}+\mathbb{E}|L(f(X_{i}),Y_{i})-L(f(X_{i}),T_{\beta_{n}}Y_{i})|$
$\displaystyle\qquad\qquad\qquad\quad+\mathbb{E}|L(f^{*}(X_{i}),T_{\beta_{n}}Y_{i})-L(f^{*}(X_{i}),Y_{i})|$
$\displaystyle\leq\mathbb{E}\\{g_{\beta_{n}}(f,Z_{i})\\}+2\lambda_{L}\mathbb{E}\\{|T_{\beta_{n}}Y_{i}-Y_{i}|\\}$
$\displaystyle\leq\mathbb{E}\\{g_{\beta_{n}}(f,Z_{i})\\}+2\lambda_{L}\mathbb{E}\big{\\{}||Y_{i}|I(|Y_{i}|>\beta_{n})\big{\\}}$
$\displaystyle\leq\mathbb{E}\\{g_{\beta_{n}}(f,Z_{i})\\}+2\lambda_{L}\mathbb{E}\\{|Y_{i}||Y_{i}|^{p-1}/\beta_{n}^{p-1}\\}$
$\displaystyle\leq\mathbb{E}\\{g_{\beta_{n}}(f,Z_{i})\\}+2\lambda_{L}\mathbb{E}|Y_{i}|^{p}/\beta_{n}^{p-1}.$

By Assumption 2, the response $Y$ has finite $p$-moment and thus
$\mathbb{E}|Y_{i}|^{p}<\infty$. Similarly,

$\displaystyle\mathbb{E}\\{g_{\beta_{n}}(f,Z_{i})\\}$
$\displaystyle=\mathbb{E}\\{g(f,Z_{i})\\}+\mathbb{E}\\{L(f^{*}(X_{i}),Y_{i})-L(f^{*}_{\beta_{n}}(X_{i}),Y_{i})\\}$
$\displaystyle\qquad\qquad\qquad+\mathbb{E}\\{L(f(X_{i}),T_{\beta_{n}}Y_{i})-L(f(X_{i}),Y_{i})\\}$
$\displaystyle\qquad\qquad\qquad+\mathbb{E}\\{L(f^{*}_{\beta_{n}}(X_{i}),Y_{i})-L(f^{*}_{\beta_{n}}(X_{i}),T_{\beta_{n}}Y_{i})\\}$
$\displaystyle\leq\mathbb{E}\\{g(f,Z_{i})\\}+\mathbb{E}|L(f(X_{i}),T_{\beta_{n}}Y_{i})-L(f(X_{i}),Y_{i})|$
$\displaystyle\qquad\qquad\qquad+\mathbb{E}|L(f^{*}_{\beta_{n}}(X_{i}),Y_{i})-L(f^{*}_{\beta_{n}}(X_{i}),T_{\beta_{n}}Y_{i})|$
$\displaystyle\leq\mathbb{E}\\{g(f,Z_{i})\\}+2\lambda_{L}\mathbb{E}|Y_{i}|^{p}/\beta_{n}^{p-1}.$

Note that above inequalities also hold for $g(f,Z_{i}^{\prime})$ and
$g_{\beta_{n}}(f,Z_{i}^{\prime})$. Then for any $f\in\mathcal{F}_{\phi}$,
define
$G_{\beta_{n}}(f,Z_{i})=\mathbb{E}_{S^{\prime}}\\{g_{\beta_{n}}(f,Z_{i}^{\prime})\\}-2g_{\beta_{n}}(f,Z_{i})$
and we have

$\mathbb{E}_{S}\Big{[}\frac{1}{n}\sum_{i=1}^{n}G(f_{j^{*}},Z_{i})\Big{]}\leq\mathbb{E}_{S}\Big{[}\frac{1}{n}\sum_{i=1}^{n}G_{\beta_{n}}(f_{j^{*}},Z_{i})\Big{]}+6\lambda_{L}\mathbb{E}|Y_{i}|^{p}/\beta_{n}^{p-1}.$
(B.5)

Besides, by Assumption 1, for any $f\in\mathcal{F}_{\phi}$ we have
$|g_{\beta_{n}}(f,Z_{i})|\leq 4\lambda_{L}\beta_{n}$ and
$\sigma^{2}_{g}(f):={\rm
Var}(g_{\beta_{n}}(f,Z_{i}))\leq\mathbb{E}\\{g_{\beta_{n}}(f,Z_{i})^{2}\\}\leq
4\lambda_{L}\beta_{n}\mathbb{E}\\{g_{\beta_{n}}(f,Z_{i})\\}$. For each $f_{j}$
and any $t>0$, let $u=t/2+{\sigma_{g}^{2}(f_{j})}/(8\lambda_{L}\beta_{n})$, by
the Bernstein inequality,

$\displaystyle
P\Big{\\{}\frac{1}{n}\sum_{i=1}^{n}G_{\beta_{n}}(f_{j},Z_{i})>t\Big{\\}}$
$\displaystyle=$ $\displaystyle
P\Big{\\{}\mathbb{E}_{S^{\prime}}\\{g_{\beta_{n}}(f_{j},Z_{i}^{\prime})\\}-\frac{2}{n}\sum_{i=1}^{n}g_{\beta_{n}}(f_{j},Z_{i})>t\Big{\\}}$
$\displaystyle=$ $\displaystyle
P\Big{\\{}\mathbb{E}_{S^{\prime}}\\{g_{\beta_{n}}(f_{j},Z_{i}^{\prime})\\}-\frac{1}{n}\sum_{i=1}^{n}g_{\beta_{n}}(f_{j},Z_{i})>\frac{t}{2}+\frac{1}{2}\mathbb{E}_{S^{\prime}}\\{g_{\beta_{n}}(f_{j},Z_{i}^{\prime})\\}\Big{\\}}$
$\displaystyle\leq$ $\displaystyle
P\Big{\\{}\mathbb{E}_{S^{\prime}}\\{g_{\beta_{n}}(f_{j},Z_{i}^{\prime})\\}-\frac{1}{n}\sum_{i=1}^{n}g_{\beta_{n}}(f_{j},Z_{i})>\frac{t}{2}+\frac{1}{2}\frac{\sigma_{g}^{2}(f_{j})}{4\lambda_{L}\beta_{n}}\\}\Big{\\}}$
$\displaystyle\leq$
$\displaystyle\exp\Big{(}-\frac{nu^{2}}{2\sigma_{g}^{2}(f_{j})+16u\lambda_{L}\beta_{n}/3}\Big{)}$
$\displaystyle\leq$
$\displaystyle\exp\Big{(}-\frac{nu^{2}}{16u\lambda_{L}\beta_{n}+16u\beta_{n}/3}\Big{)}$
$\displaystyle\leq$
$\displaystyle\exp\Big{(}-\frac{1}{16+16/3}\cdot\frac{nu}{\lambda_{L}\beta_{n}}\Big{)}$
$\displaystyle\leq$
$\displaystyle\exp\Big{(}-\frac{1}{32+32/3}\cdot\frac{nt}{\lambda_{L}\beta_{n}}\Big{)}.$

This leads to a tail probability bound of
$\sum_{i=1}^{n}G_{\beta_{n}}(f_{j^{*}},Z_{i})/n$, that is

$P\Big{\\{}\frac{1}{n}\sum_{i=1}^{n}G_{\beta_{n}}(f_{j^{*}},Z_{i})>t\Big{\\}}\leq
2\mathcal{N}_{2n}\exp\Big{(}-\frac{1}{43}\cdot\frac{nt}{\lambda_{L}\beta_{n}}\Big{)}.$

Then for $a_{n}>0$,

$\displaystyle\mathbb{E}_{S}\Big{[}\frac{1}{n}\sum_{i=1}^{n}G_{\beta_{n}}(f_{j^{*}},Z_{i})\Big{]}\leq$
$\displaystyle
a_{n}+\int_{a_{n}}^{\infty}P\Big{\\{}\frac{1}{n}\sum_{i=1}^{n}G_{\beta_{n}}(f_{j^{*}},Z_{i})>t\Big{\\}}dt$
$\displaystyle\leq$ $\displaystyle
a_{n}+\int_{a_{n}}^{\infty}2\mathcal{N}_{2n}\exp\Big{(}-\frac{1}{43}\cdot\frac{nt}{\lambda_{L}\beta_{n}}\Big{)}dt$
$\displaystyle\leq$ $\displaystyle
a_{n}+2\mathcal{N}_{2n}\exp\Big{(}-a_{n}\cdot\frac{n}{43\lambda_{L}\beta_{n}}\Big{)}\frac{43\lambda_{L}\beta_{n}}{n}.$

Choosing $a_{n}=\log(2\mathcal{N}_{2n})\cdot{43\lambda_{L}\beta_{n}}/{n}$, we
have

$\mathbb{E}_{S}\Big{[}\frac{1}{n}\sum_{i=1}^{n}G_{\beta_{n}}(f_{j^{*}},Z_{i})\Big{]}\leq\frac{43\lambda_{L}\beta_{n}(\log(2\mathcal{N}_{2n})+1)}{n}.$
(B.6)

Setting $\delta=1/n$, $\beta_{n}=c_{1}\max\\{\mathcal{B},n^{1/p}\\}$ and
combining (B.3), (B.4), (B.5) and (B.6), we get

$\mathbb{E}\big{\\{}\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\big{\\}}\leq\frac{c_{2}\lambda_{L}\mathcal{B}\log\mathcal{N}_{2n}(\frac{1}{n},\|\cdot\|_{\infty},\mathcal{F}_{\phi})}{n^{1-1/p}}+2\big{\\{}\mathcal{R}(f^{*}_{\phi})-\mathcal{R}(f^{*})\big{\\}},$
(B.7)

where $c_{2}>0$ is a constant not depending on $n,d,\mathcal{B}$ and
$\lambda_{L}$. This shows (3.4).

#### Step 3: Bounding the covering number

Lastly, we will give an upper bound on the covering number by the VC dimension
of $\mathcal{F}_{\phi}$ through its parameters. Denote ${\rm
Pdim}(\mathcal{F}_{\phi})$ by the pseudo dimension of $\mathcal{F}_{\phi}$, by
Theorem 12.2 in Anthony and Bartlett (1999), for $2n\geq{\rm
Pdim}(\mathcal{F}_{\phi})$

$\mathcal{N}_{2n}(\frac{1}{n},\|\cdot\|_{\infty},\mathcal{F}_{\phi})\leq\Big{(}\frac{2e\mathcal{B}n^{2}}{{\rm
Pdim}(\mathcal{F}_{\phi})}\Big{)}^{{\rm Pdim}(\mathcal{F}_{\phi})}.$

Besides, according to Theorem 3 and Theorem 6 in Bartlett et al. (2019), there
exist universal constants $c$ and $C$ such that

$c\cdot\mathcal{S}\mathcal{D}\log(\mathcal{S}/\mathcal{D})\leq{\rm
Pdim}(\mathcal{F}_{\phi})\leq C\cdot\mathcal{S}\mathcal{D}\log(\mathcal{S}).$

Combining the upper bound of the covering number and pseudo dimension,
together with (B.7), we have

$\mathbb{E}\big{\\{}\mathcal{R}(\hat{f}_{\phi})-\mathcal{R}(f^{*})\big{\\}}\leq
c_{3}\lambda_{L}\mathcal{B}\frac{\log(n)\mathcal{S}\mathcal{D}\log(\mathcal{S})}{n^{1-1/p}}+2\big{\\{}\mathcal{R}(f^{*}_{\phi})-\mathcal{R}(f^{*})\big{\\}},$
(B.8)

for some constant $c_{3}>0$ not depending on
$n,d,\lambda_{L},\mathcal{B},\mathcal{S}$ and $\mathcal{D}$, which leads to
(3.5). This completes the proof of Lemma 3.1. ∎

###  B.2 Proof of Lemma 3.2

The proof is straightforward under Assumption 1. For the risk minimizer
$f^{*}$ and any $f\in\mathcal{F}_{\phi}$, we have

$\displaystyle R(f)-R^{\phi}(f^{*})$
$\displaystyle=\mathbb{E}\\{L(f(X),Y)-L(f^{*}(X),Y)\\}\leq\mathbb{E}\\{\lambda_{L}|f(X)-f^{*}(X)|\\},$

thus

$\inf_{f\in\mathcal{F}_{\phi}}\\{R(f)-R^{\phi}(f^{*})\\}\leq\lambda_{L}\inf_{f\in\mathcal{F}_{\phi}}\mathbb{E}|f(X)-f^{*}(X)|=:\lambda_{L}\inf_{f\in\mathcal{F}_{\phi}}\|f-f^{*}\|_{L^{1}(\nu)},$

where $\nu$ denotes the marginal probability measure of $X$ and
$\mathcal{F}_{\phi}=\mathcal{F}_{\mathcal{D},\mathcal{W},\mathcal{U},\mathcal{S},\mathcal{B}}$