# Robust Functional Regression with Discretely Sampled Predictors

Ioannis Kalogridis1 and Stanislav Nagy2

(1 Department of Mathematics, KU Leuven
2 Department of Probability and Mathematical Statistics, Charles University
)

###### Abstract

The functional linear model is an important extension of the classical
regression model allowing for scalar responses to be modeled as functions of
stochastic processes. Yet, despite the usefulness and popularity of the
functional linear model in recent years, most treatments, theoretical and
practical alike, suffer either from (i) lack of resistance towards the many
types of anomalies one may encounter with functional data or (ii) biases
resulting from the use of discretely sampled functional data instead of
completely observed data. To address these deficiencies, this paper introduces
and studies the first class of robust functional regression estimators for
partially observed functional data. The proposed broad class of estimators is
based on thin-plate splines with a novel computationally efficient quadratic
penalty, is easily implementable and enjoys good theoretical properties under
weak assumptions. We show that, in the incomplete data setting, both the
sample size and discretization error of the processes determine the asymptotic
rate of convergence of functional regression estimators and the latter cannot
be ignored. These theoretical properties remain valid even with multi-
dimensional random fields acting as predictors and random smoothing
parameters. The effectiveness of the proposed class of estimators in practice
is demonstrated by means of a simulation study and a real-data example.

##  1 Introduction

In recent years, technological innovations and improved storage capabilities
have led practitioners to observe and record increasingly complex high-
dimensional data that are characterized by an underlying functional structure.
Such data are nowadays commonly referred to as functional data and relevant
research has been enjoying considerable popularity, following works such as
Ramsay, (1982), Ramsay and Dalzell, (1991) and Ramsay and Silverman, (2005).
While the field of functional data analysis (FDA) has become very broad with
many specialized subpaths, see, e.g., Ferraty and Vieu, (2006); Horváth and
Kokoszka, (2012); Kokoszka and Reimherr, (2017), the functional linear model
continues to occupy a prominent place within FDA. This model stipulates that
an $\mathcal{L}^{2}(\mathcal{I})$-process $\\{X(t),t\in\mathcal{I}\\}$, for
some bounded interval $\mathcal{I}\subset\mathbbm{R}$, influences the response
variable $Y$ through an unknown coefficient function
$\beta_{0}\in\mathcal{L}^{2}(\mathcal{I})$, viz,

$\displaystyle
Y=\alpha_{0}+\int_{\mathcal{I}}X(t)\beta_{0}(t)\,\mathrm{d}\,t+\epsilon,$ (1)

where $\alpha_{0}\in\mathbbm{R}$ is an unknown offset (intercept) term and
$\epsilon$ a random error, which is assumed to be independent of $X$.
Typically, $\epsilon$ is also assumed to possess finite second moments, but,
as we shall see, moment assumptions on the errors will not be needed for the
theoretical results of this paper.

Assuming that the researcher readily possesses a random sample
$(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$ following model (1), there is a wealth of
possible estimation methods that she could employ in order to estimate the
unknown parameters
$(\alpha_{0},\beta_{0})\in\mathbbm{R}\times\mathcal{L}^{2}(\mathcal{I})$.
Since the covariance operator of the process $\\{X(t),t\in\mathcal{I}\\}$ does
not possess a bounded inverse, such estimation methods are invariably based on
lower rank representations of either $X$ or $\beta_{0}$, regularization
through appropriate penalties or a combination thereof. We mention, in
particular, the popular approaches of functional principal component
regression and penalized basis expansions regression expounded by Ramsay and
Silverman, (2005) and Kokoszka and Reimherr, (2017) and theoretically
investigated by Hall and Horowitz, (2007) and Cardot et al., (2003); Li and
Hsing, (2007), respectively. More advanced approaches include the hybrid
method of Reiss and Ogden, (2007) which combines functional principal
components and penalized spline expansions, as well as the reproducing kernel
Hilbert space estimator of Yuan and Cai, (2010). To these methods we may add
the sparsity-inducing proposal of James et al., (2009), which ensures the
smoothness and interpretability of the estimates by taking advantage of the
sparsity of higher order derivatives of $\beta_{0}$ when those exist.

From the point of view of robustness, that is, resistance against atypical
observations and model misspecification, the aforementioned proposals are not
satisfactory, as they are all generalizations of the ordinary least-squares
estimator. The latter is known to be very susceptible even to mild deviations
from the ideal model assumptions (see, e.g., Huber and Ronchetti,, 2009). To
overcome this weakness of least-squares based estimators, several proposals
have been put forth in the literature over the years. These proposals
essentially amount to replacing the square loss with a more slowly increasing
loss function thereby ensuring that the influence of atypical observations and
model misspeficication on the estimates is better controlled. Examples of such
an approach include Qingguo, (2017) and Kalogridis and Van Aelst, (2019), who
generalized the functional principal component estimator of Hall and Horowitz,
(2007). Later, Boente et al., (2020); Kalogridis and Van Aelst, (2023)
proposed robust spline estimators based on the principle of MM-estimation, and
Shin and Lee, (2016) generalized the work of Yuan and Cai, (2010) by allowing
for a large number of loss functions.

Nevertheless, a severe drawback of all aforementioned theoretical
contributions, robust and non-robust alike, is their reliance on completely
observed functional predictors $X_{1},\ldots,X_{n}$ throughout $\mathcal{I}$.
In practice, this is an unrealistic assumption and in the vast majority of
cases either computational or resource constraints dictate that the curves are
only recorded at a finite grid of points within $\mathcal{I}$:
$t_{1},\ldots,t_{p}$, say. Functional regression with discretely sampled data
is a setting that has received much less attention despite the prevalence of
partially observed functional data in practice. To the best of our knowledge,
the only theoretical contributions for this setting are the works of Crambes
et al., (2009) and Kato, (2012), who proposed least-squares smoothing-spline
estimation with a slightly modified penalty and quantile regression on the
leading functional principal components, respectively. The estimator of
Crambes et al., (2009) is non-robust on account of the least-squares criterion
employed therein whereas the $L_{1}$-estimator of Kato, (2012), while robust
towards heavy-tailed error distributions, is inefficient in the case of clean
data, e.g., under light tailed sub-Gaussian errors. Moreover, both works rely
on the assumption that $p$, the number of discretization points, is large
relative to $n$ in order for the discretization error to be asymptotically
negligible and it is not clear how these estimators behave for sparsely
observed functional data.

To overcome the drawbacks of existing methods either with respect to lack of
robustness or reliance on densely observed functional data, this paper
introduces a flexible family of penalized thin-plate spline M-estimators that
can operate under both densely and sparsely observed functional data. This
flexibility is accomplished through the introduction of a novel penalty
functional that is inspired by norms on multi-dimensional Sobolev spaces,
which are widely used in the study of partial differential equations.
Robustness towards atypical observations and model misspecification can be
achieved through appropriate selection of the loss function. Possible loss
functions in that respect include not only the quantile loss, but also the
Huber loss, which can lead to estimators that are resistant to atypical
observations and model misspecification while remaining highly efficient under
light-tailed errors.

We study the proposed class of estimators not only under the classical one-
dimensional functional linear model with $\mathcal{I}\subset\mathbbm{R}$, but
also with general $\mathcal{I}\subset\mathbbm{R}^{d}$. This extension of the
functional linear model is particularly important, as it allows more general
random objects, such as images, to be used as explanatory variables for the
response variable $Y$. To the best of our knowledge, functional regression
with general, possibly multi-dimensional, random fields has not been
considered previously even for relatively simple least-squares estimators so
that the present contributions is the first of its kind. Moreover, unlike the
theoretical results of Crambes et al., (2009); Kato, (2012), our theoretical
results hold even with random smoothing parameters. Since smoothing parameters
are normally selected in a data-dependent way, treating them as random
variables leads to considerably more realistic and useful results.

The rest of the paper is structured as follows. Section 2 introduces the
proposed class of penalized M-estimators and establishes its existence under
general conditions. Section 3 is dedicated to the asymptotic study of our
estimators. We obtain a useful decomposition of the frequently used mean
squared error highlighting the delicate interplay between the sample size and
level of discretization, which jointly determine the asymptotic rate of
convergence of the proposed family of estimators. Section 4 outlines an
effective computational algorithm and a model selection procedure, which works
well in a variety of settings as we demonstrate by means of our numerical
experiments and real-data example in Section 5 and Section 6, respectively.
Finally, Section 7 briefly discusses two interesting directions for future
research.

##  2 Thin-Plate Splines for Functional Regression

###  2.1 The Proposed Family of Estimators

Throughout this section we suppose that $X_{1},\ldots,X_{n}$ are independent
and identically distributed (i.i.d.) copies of the random field
$\\{X(\mathbf{t}),\mathbf{t}\in\mathcal{I}\\}$ for some
$\mathcal{I}\subset\mathbbm{R}^{d}$. As commonly done in functional
regression, we also assume that $X$ is a second-order process over
$\mathcal{I}$, or, equivalently, $\mathbb{E}\\{\|X\|^{2}\\}<\infty$ with
$\|\cdot\|$ denoting the classical $\mathcal{L}^{2}(\mathcal{I})$-norm. The
response variables $Y_{1},\ldots,Y_{n}$, are assumed to be generated according
to the model

$\displaystyle
Y_{i}=\alpha_{0}+\int_{\mathcal{I}}X_{i}(\mathbf{t})\beta_{0}(\mathbf{t})\,\mathrm{d}\,\mathbf{t}+\epsilon_{i},\quad(i=1,\ldots,n),$
(2)

for unknown quantities
$(\alpha_{0},\beta_{0})\in\mathbbm{R}\times\mathcal{L}^{2}(\mathcal{I})$ that
are the objects of interest and i.i.d. errors $\\{\epsilon_{i}\\}_{i=1}^{n}$.

For completely observed random fields $\\{X_{i}\\}_{i=1}^{n}$ most current
estimation proposals would boil down to the minimization of

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\rho\left(Y_{i}-\alpha-\int_{\mathcal{I}}X_{i}(\mathbf{t})\beta(\mathbf{t})\,\mathrm{d}\,\mathbf{t}\right)+\lambda
J(\beta),$ (3)

over $(\alpha,\beta)\in\mathbbm{R}\times\Theta$, where $\Theta$ is a suitable
subspace of functions, $\rho:\mathbbm{R}\to\mathbbm{R}_{+}$ represents an
a-priori chosen loss function, e.g., the square loss $\rho(x)=x^{2}$ or the
$L_{1}$ loss $\rho(x)=|x|$, $J:\Theta\to\mathbbm{R}_{+}$ is a penalty
functional, usually a semi-norm, and $\lambda\geq 0$ a penalty parameter that
regulates the trade-off between smoothness and fidelity to the data. In the
one-dimensional case, that is, $d=1$, popular choices of $\Theta$ include the
subspace spanned by a small number of eigenfunctions of the covariance
operator of $X$ or a spline subspace, whereas $J(\beta)$ would be a function
of the derivatives of $\beta$ or it could even be omitted entirely whenever
$\Theta$ is a finite dimensional subspace and thus the risk of overfitting is
mitigated. The interested reader is referred to Ramsay and Silverman, (2005,
Chapter 15) for a discussion of the most popular estimation methods for
scalar-on-function regression in the one-dimensional setting.

In cases where the data is only discretely observed, at
$\mathbf{t}_{1},\ldots,\mathbf{t}_{p}\in\mathbbm{R}^{d}$, say, the above
recipe is not applicable, as the $\mathcal{L}^{2}(\mathcal{I})$-inner product
between the random fields and each candidate coefficient function $\beta$
cannot be computed. Instead, a popular method of estimation with discretely
sampled functional data consists of approximating the integral in (3) with a
Riemann sum and applying any of the aforementioned methods. In essence, this
strategy consists of ignoring the error associated with the discrete nature of
the data. However, while appealing due to its simplicity, the theoretical
effects of such a practice are not entirely understood. We take a different
approach in this paper by dealing with the Riemann sum directly and deriving
our estimators from this Riemann sum in combination with a suitable penalty
functional $J(\cdot)$ on the space of smooth functions on $\mathbbm{R}^{d}$.

To explain our method in detail, let $\\{A_{j}\\}_{j=1}^{p}$ denote a disjoint
partition of $\mathcal{I}$ such that $\mathbf{t}_{j}\in A_{j}$ and
$\mu(A_{j})>0$ for each $j=1,\dots,p$, where $\mu(A)$ denotes the Lebesgue-
measure (volume) of a measurable set $A$. With this notation we may rewrite
the functional linear model (2) as

$\displaystyle
Y_{i}=\alpha_{0}+\sum_{j=1}^{p}X_{i}(\mathbf{t}_{j})\beta_{0}(\mathbf{t}_{j})\mu(A_{j})+d_{i}+\epsilon_{i},\quad(i=1,\ldots,n),$

where $\\{d_{i}\\}_{i=1}^{n}$ represent the unobserved discretization errors,
i.e.,

$\displaystyle
d_{i}=\int_{\mathcal{I}}X_{i}(\mathbf{t})\beta_{0}(\mathbf{t})\,\mathrm{d}\,\mathbf{t}-\sum_{j=1}^{p}X_{i}(\mathbf{t}_{j})\beta_{0}(\mathbf{t}_{j})\mu(A_{j}),\quad(i=1,\ldots,n).$
(4)

In order to both deal with the discretization error theoretically and allow
for a great variety of coefficient functions $\beta_{0}$, we assume in this
paper that $\beta_{0}$ is a smooth function on $\mathbbm{R}^{d}$ in the sense
of possessing partial derivatives of a given order. A rich space of functions
fulfilling this property is the Sobolev-Hilbert space of functions of order
$m\geq 1$, $\mathcal{H}^{m}(\mathbbm{R}^{d})$, defined by

$\displaystyle\mathcal{H}^{m}(\mathbbm{R}^{d})=\left\\{f:\mathbbm{R}^{d}\to\mathbbm{R},\
\frac{\partial^{m}f(t_{1},\ldots,t_{d})}{\partial t_{1}^{m_{1}}\ldots\partial
t_{d}^{m_{d}}}\ \text{exists for all}\ m_{1}+\ldots+m_{d}=m\ \text{and}\
I_{m}^{2}(f)<\infty\right\\},$

where the semi-norm $I_{m}:\mathcal{H}^{m}(\mathbbm{R}^{d})\to\mathbbm{R}_{+}$
is given by

$\displaystyle
I_{m}^{2}(f)=\sum_{m_{1}+\ldots+m_{d}=m}\binom{m}{m_{1},\ldots,m_{d}}\int_{\mathbbm{R}}\ldots\int_{\mathbbm{R}}\left(\frac{\partial^{m}f(t_{1},\ldots,t_{d})}{\partial
t_{1}^{m_{1}}\ldots\partial
t_{d}^{m_{d}}}\right)^{2}\,\mathrm{d}\,t_{1}\ldots\,\mathrm{d}\,t_{d}.$ (5)

In other words, $\mathcal{H}^{m}(\mathbbm{R}^{d})$ consists of all functions
with square integrable partial derivatives of total order $m$. It can be shown
(see, e.g., Adams and Fournier,, 2003, Chapter 3) that
$\mathcal{H}^{m}(\mathbbm{R}^{d})$ is the completion of the space of functions
with continuous partial derivatives of total order $m$ under a suitable norm
and as such it affords us slightly greater generality. In what follows we
place no restrictions on the coefficient function $\beta_{0}$, except that it
is an element of the Sobolev-Hilbert space of functions, i.e.,
$\beta_{0}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$.

If we accept that $\beta_{0}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$, it is
intuitively appealing to minimize a suitably chosen objective function over
the whole space $\mathbbm{R}\times\mathcal{H}^{m}(\mathbbm{R}^{d})$ in order
to recover the unknown $(\alpha_{0},\beta_{0})$. To achieve both resistance
against atypical observations $(X_{i},Y_{i})$ and high efficiency in regular
data we propose to estimate $(\alpha_{0},\beta_{0})$ by
$(\widehat{\alpha}_{n},\widehat{\beta}_{n})$ solving

$\displaystyle\left(\widehat{\alpha}_{n},\widehat{\beta}_{n}\right)=\operatorname*{argmin}_{(\alpha,\beta)\in\mathbbm{R}\times\mathcal{H}^{m}\left(\mathbbm{R}^{d}\right)}\left[\frac{1}{n}\sum_{i=1}^{n}\rho\left(\frac{Y_{i}-\alpha-\sum_{j=1}^{p}X_{i}(\mathbf{t}_{j})\beta(\mathbf{t}_{j})\mu(A_{j})}{\widehat{\sigma}_{n}}\right)+\lambda
J_{m}^{2}(\beta)\right],$ (6)

where $\rho:\mathbbm{R}\to\mathbbm{R}_{+}$ is a convex loss function,
$\widehat{\sigma}_{n}$ is an auxiliary scale estimator for the scale of the
error $\epsilon$, $\lambda$ is a penalty parameter, and the penalty functional
$J_{m}:\mathcal{H}^{m}(\mathbbm{R}^{d})\to\mathbbm{R}_{+}$ is given by

$\displaystyle
J_{m}^{2}(\beta)=\sum_{j=1}^{p}\left|\beta(\mathbf{t}_{j})\right|^{2}\mu(A_{j})+I_{m}^{2}(\beta),\quad\beta\in\mathcal{H}^{m}(\mathbbm{R}^{d}).$
(7)

That is, our penalty functional $J_{m}^{2}$ is the sum of two semi-norms: the
empirical squared semi-norm at the discretization points
$\\{\mathbf{t}_{j}\\}_{j=1}^{p}$ and the Sobolev semi-norm $I_{m}^{2}$ from
(5).

It is worth commenting on two particular aspects of our estimation framework,
namely on the flexibility allowed with respect to the choice of the loss
function $\rho$ and the slightly unconventional form of the penalty functional
$J_{m}^{2}(\cdot)$. Indeed, our framework allows not only for the square loss
$\rho(x)=x^{2}$, which is known to lead to efficient but highly sensitive
estimates, and the $L_{1}$ loss $\rho(x)=|x|$, which leads to resistant but
inefficient estimates under light-tailed errors. It also allows for a
compromise between these two extremes through the use of many other loss
functions, such as the logistic loss $\rho(x)=2x+4\log(1+e^{-x})$ and the
celebrated Huber loss given by

$\displaystyle\rho_{k}(x)=\begin{cases}x^{2}&|x|<k\\\ 2k|x|-k^{2}&|x|\geq
k,\end{cases}$ (8)

for some $k>0$ regulating the mixing of the square and $L_{1}$ losses. For
loss functions that are not power functions, the scale estimate,
$\widehat{\sigma}_{n}$ in (6), is necessary to ensure approximate scale
equivariance of the estimates. That is, to ensure that our estimates for
$(\alpha_{0},\beta_{0})$ do not heavily depend on the measurement units of the
response variables $\\{Y_{i}\\}_{i=1}^{n}$. Standardization with auxiliary
scale estimates in order to achieve equivariance has been standard practice in
classical (unpenalized) robust regression (see, e.g. Maronna et al.,, 2019).
Such scale estimates can be defined without assuming the existence of any
moments of the error $\epsilon$, see Section 3 for an example.

The penalty functional (7) employed herein is a modification of the thin-plate
penalty function appearing in the classical monographs of Wahba, (1990) and
Green and Silverman, (1994) where it only involves $I_{m}^{2}(\beta)$. It will
be shown, however, that in our framework the addition of
$\sum_{j=1}^{p}|\beta(\mathbf{t}_{j})|^{2}\mu(A_{j})$ is essential in order to
ensure the existence of the proposed estimators without heavy assumptions on
the process $\\{X(\mathbf{t}),\mathbf{t}\in\mathcal{I}\\}$. We will also show
that despite the inclusion of this additional term the solution to the
minimization problem may still be found in the space of natural thin-plate
splines, as is the case when only employing the classical thin-plate penalty
functional. The practical interpretation of $J_{m}^{2}$ is similar to that of
$I_{m}^{2}$: a premium is placed on functions with “large” partial derivatives
of total order $m$ but, unlike $I_{m}^{2}$, $J_{m}^{2}$ also penalizes
functions assuming large values at the $\\{\mathbf{t}_{j}\\}_{j=1}^{p}$. Thus,
$J_{m}^{2}$ enforces the smoothness of $\widehat{\beta}_{n}$ not only by
penalizing roughness but also by shrinking it towards zero.

###  2.2 Existence of the Estimators

We now establish the existence of the penalized M-estimators defined in (6).
We shall require the following two general conditions involving the set
$\mathcal{I}$, its boundary $\partial\mathcal{I}$ and the discretization
points $\\{\mathbf{t}_{j}\\}_{j=1}^{p}$. We use
$\left\|\cdot\right\|_{\mathbbm{R}^{d}}$ to denote the Euclidean norm on
$\mathbbm{R}^{d}$.

  1. (A1)

The discretization points $\\{\mathbf{t}_{j}\\}_{j=1}^{p}$ are contained in a
bounded open set $\mathcal{I}\subset\mathbbm{R}^{d}$ whose boundary,
$\partial\mathcal{I}$, satisfies the uniform cone condition of Adams and
Fournier, (2003, p. 83).

  2. (A2)

Define the quantities

$\displaystyle h_{\max}$
$\displaystyle=\sup_{\mathbf{t}\in\mathcal{I}}\min_{1\leq j\leq
p}\|\mathbf{t}-\mathbf{t}_{j}\|_{\mathbbm{R}^{d}}$ $\displaystyle h_{\min}$
$\displaystyle=\min_{j\neq
k}\|\mathbf{t}_{j}-\mathbf{t}_{k}\|_{\mathbbm{R}^{d}}.$

Then, for all large $n$ and $p$, there exists a finite constant $B$ such that
$h_{\max}/h_{\min}\leq B$.

Condition (A1) is technical in nature and precludes very irregular boundaries.
It is satisfied quite generally. For example, it is valid for balls and
rectangles in $\mathbbm{R}^{d}$. Condition (A2) is essentially a density
condition for the discretization points; it is required that discretization
points are unique and cover $\mathcal{I}$ sufficiently well. Both of these
implications follow from the bound $h_{\max}/h_{\min}\leq B$. This condition
was introduced by Utreras, (1988) for the study of thin-plate splines in the
context of non-parametric regression.

Proposition 1 below establishes the existence of the estimators and provides a
useful characterization that will form the basis of our computational
algorithm in Section 4. Prior to the statement of Proposition 1, we remind the
reader that a thin-plate spline of order $m$ with knots at
$\\{\mathbf{t}_{j}\\}_{j=1}^{p}$ is any function
$g:\mathbbm{R}^{d}\to\mathbbm{R}$ of the form

$\displaystyle
g\left(\mathbf{t}\right)=\sum_{j=1}^{p}\gamma_{j}\eta_{m,d}\left(\left\|\mathbf{t}-\mathbf{t}_{j}\right\|_{\mathbbm{R}^{d}}\right)+\sum_{k=1}^{M}\delta_{k}\phi_{k}\left(\mathbf{t}\right),$
(9)

with $\eta_{m,d}:\mathbbm{R}_{+}\to\mathbbm{R}$ given by

$\displaystyle\eta_{m,d}(x)=\begin{cases}\frac{(-1)^{m+1+d/2}}{2^{2m-1}\pi^{d/2}(m-1)!(m-d/2)!}x^{2m-d}\log(x)&d\
\text{even}\\\ \frac{\Gamma(d/2-m)}{2^{2m}\pi^{d/2}(m-1)!}x^{2m-d}&d\
\text{odd},\end{cases}$

where $\Gamma(\cdot)$ denotes Euler’s gamma function. The functions
$\\{\phi_{k}\\}_{k=1}^{M}$ are (any) basis for the space of polynomials on
$\mathbbm{R}^{d}$ of total order less than $m$, which has dimension
$M=\binom{m+d-1}{d}$. Moreover, $g$ is called a natural thin-plate spline, if,
in addition to (9), $I_{m}(g)<\infty$ where $I_{m}(\cdot)$ is defined in (5).
As we explain in Section 4 below, this condition places an orthogonality
restriction on the coefficients
$\boldsymbol{\gamma}=(\gamma_{1},\ldots,\gamma_{p})$.

######  Proposition 1.

Suppose that $\rho$ is a convex loss function, $2m>d$ and (A1) and (A2) hold.
Then, there exists a solution to (6) in
$\mathbbm{R}\times\mathcal{H}^{m}(\mathbbm{R}^{d})$ denoted by
$(\widehat{\alpha}_{n},\widehat{\beta}_{n})$. Moreover, if the set
$\\{\mathbf{t}_{j}\\}_{j=1}^{p}$ contains a $\mathcal{P}_{m}$-unisolvent set,
then $\widehat{\beta}_{n}$ is necessarily a natural thin-plate spline of order
$m$ with knots at $\\{\mathbf{t}_{j}\\}_{j=1}^{p}$.

The condition $2m>d$ in the statement of Proposition 1 is both necessary and
sufficient for $\mathcal{H}^{m}(\mathbbm{R}^{d})$ to be a reproducing kernel
Hilbert space of functions wherein point evaluation is well-defined and
therefore this condition cannot be weakened. This condition is satisfied for
all integers $m\geq 1$ when $d=1$, but for larger $d$ it may require us to
increase $m$ accordingly, e.g., $m>2$ whenever $d=2$. The additional condition
that $\\{\mathbf{t}_{j}\\}_{j=1}^{p}$ contains a $\mathcal{P}_{m}$-unisolvent
set is a general condition requiring that there exists a subset of
$\\{\mathbf{t}_{j}\\}_{j=1}^{p}$ that allows for unique polynomial
interpolation. Equivalently, the $p\times M$ matrix $\boldsymbol{\Phi}$ with
elements $\phi_{j}(\mathbf{t}_{i})$, $i=1,\dots,p$, $j=1,\dots,M$, needs to
have full column rank. In the unidimensional case $d=1$ it is easy to show
that this condition is satisfied if the points $\\{t_{j}\\}_{j=1}^{p}$ are
distinct and $p\geq m$, which is a classical condition for the existence of
smoothing spline estimators (see, e.g. Green and Silverman,, 1994, Theorem
2.3). Lastly, we remark that although Proposition 1 does not establish the
uniqueness of the minimizers, uniqueness can be established rather easily
whenever $\rho$ is a strictly convex loss function, such as the square or
logistic losses.

######  Corollary 1.

Suppose that the condititions of Proposition 1 are satisfied and that $\rho$
is a strictly convex function. Then, the solution
$(\widehat{\alpha}_{n},\widehat{\beta}_{n})$ to (6) in
$\mathbbm{R}\times\mathcal{H}^{m}(\mathbbm{R}^{d})$ is unique.

It is interesting to relate the penalty functional $J_{m}$ employed herein
with the penalty functional earlier used by Crambes et al., (2009) in the
specific least-squares case, i.e., $\rho(x)=x^{2}$. In particular, for the
one-dimensional case, Crambes et al., (2009) proposed using

$\displaystyle\widetilde{J}_{m}^{2}(\beta)=\frac{1}{p}\sum_{j=1}^{p}(\mathcal{P}\beta)^{2}(t_{j})+\int_{0}^{1}|\beta^{(m)}(t)|^{2}\,\mathrm{d}\,t,$

where $\mathcal{P}$ is the projection operator onto the space of discretized
order $m$ polynomials, that is, the polynomials of order $m$ evaluated only at
$\\{t_{j}\\}_{j=1}^{p}$. Note that for $d=1$,
$I_{m}^{2}(\beta)=\int_{\mathbbm{R}}|\beta^{(m)}(t)|^{2}\,\mathrm{d}\,t$ but
for $\mathcal{I}=(0,1)$, as is assumed by Crambes et al., (2009), it can be
shown (see, e.g., Green and Silverman,, 1994, Theorem 2.3) that we can reduce
the range of integration from $\mathbbm{R}$ to $(0,1)$ without loss of
generality. This simplification is not valid for $d>1$. For
$\mathcal{I}=(0,1)$ it is easy to see that condition (A1) is satisfied while
assumption (A2) is satisfied, e.g., if $t_{j}=(j-1/2)/p$ in which case
$h_{\max}\lesssim p^{-1}$ and $h_{\min}=p^{-1}$. Although seemingly different,
$J_{m}(\beta)$ and $\widetilde{J}_{m}(\beta)$ turn out to be equivalent semi-
norms, as Proposition 2 shows.

######  Proposition 2.

Suppose that $\mathcal{I}=(0,1)$. Then, there exist $0<c_{1}\leq c_{2}<\infty$
with the property that $c_{1}\widetilde{J}_{m}(\beta)\leq J_{m}(\beta)\leq
c_{2}\widetilde{J}_{m}(\beta)$ for all $\beta\in\mathcal{H}^{m}(\mathbbm{R})$.

Hence, for $d=1$ it does not make a difference whether one employs
$J_{m}(\beta)$ or $\widetilde{J}_{m}(\beta)$, as these penalty functionals are
qualitatively similar and share the same null space, which, by the proof of
Proposition 1 in the Appendix, consists solely of the zero function. The
advantage of our penalty, however, is that it avoids the computation of the
projection onto the discretized polynomials so that the associated
computational burden is reduced. For higher dimensions the computational
simplicity of our penalty functional is even more appealing as in that case
the projection $\mathcal{P}$ would have to be on a space of polynomials whose
dimension grows essentially like $m^{d}/(d!)$.

##  3 Asymptotic Properties

###  3.1 Thin-Plate Splines without Scale Estimation

To simplify our notation we shall assume in this section that the response
variables are centered so that $\alpha_{0}=0$ and the object of interest is
the coefficient function $\beta_{0}$, as is typically the case. Furthermore,
we also assume for simplicity that $\mu(A_{j})=\mu(\mathcal{I})/p$, i.e.,
$\mathcal{I}$ is partitioned into $p$ sets of equal volume. We begin our study
of the asymptotic properties of our estimators by first examining the
asymptotic behavior of thin-plate spline estimators that do not require
standardization with a scale $\widehat{\sigma}_{n}$. In other words we can set
$\widehat{\sigma}_{n}=1$ in (6). As noted in Section 2, the $L_{1}$ estimator
with $\rho(x)=|x|$ belongs to this class of estimators and so does its
quantile generalization with $\rho_{\tau}(x)=2x(\tau-I(x<0))$ for
$\tau\in(0,1)$. Symmetry of the loss function is not required in our
treatment, hence this important class of estimators is also covered by our
theory.

Our aim is to establish the rate of convergence of our estimators as a
function of the sample size and the discretization error with respect to the
distance given by

$\displaystyle\left\|\widehat{\beta}_{n}-\beta_{0}\right\|_{n}^{2}=\frac{1}{n}\sum_{i=1}^{n}\left|\langle
X_{i},\widehat{\beta}_{n}-\beta_{0}\rangle\right|^{2},$ (10)

where $\langle\cdot,\cdot\rangle$ denotes the standard
$\mathcal{L}^{2}(\mathcal{I})$ inner product. This distance is related to the
centered empirical covariance operator
$\Gamma_{n}(\beta)=\sum_{i=1}^{n}X_{i}\langle X_{i},\beta\rangle/n$ of a
process $X$ through
$\|\widehat{\beta}_{n}-\beta_{0}\|_{n}^{2}=\langle\Gamma_{n}(\widehat{\beta}_{n}-\beta_{0}),\widehat{\beta}_{n}-\beta_{0}\rangle$.
It bears the intuitive interpretation of the mean squared error resulting from
using $\langle X_{i},\widehat{\beta}_{n}\rangle$ to predict
$\mathbbm{E}\\{Y_{i}|X_{i}\\}=\langle X_{i},\beta_{0}\rangle$ for
$i=1,\ldots,n$. As a first step in our analysis, we establish a rate of
convergence with respect to the discretized variant of (10) given by

$\displaystyle\left\|\widehat{\beta}_{n}-\beta_{0}\right\|_{n,p}^{2}=\frac{\mu(\mathcal{I})}{np^{2}}\sum_{i=1}^{n}\left|\sum_{j=1}^{p}X_{i}(\mathbf{t}_{j})\left(\widehat{\beta}_{n}(\mathbf{t}_{j})-\beta_{0}(\mathbf{t}_{j})\right)\right|^{2}.$
(11)

Clearly, (11) may be viewed as a Riemann approximation to (10) based on the
discretization points $\\{\mathbf{t}_{j}\\}_{j=1}^{p}$ and partition sets
$\\{A_{j}\\}_{j=1}^{p}$ with $\mu(A_{j})=\mu(\mathcal{I})/p$. The assumptions
that we need for our theoretical development are as follows.

  1. (A3)

There exist a $t\in(0,1)$ such that
$\mathbb{E}\\{e^{t\|X\|_{\infty}^{2}}\\}<\infty$ and a $C_{1}>0$ and
$\kappa\in(0,1]$ such that

$\displaystyle\mathbb{E}\left\\{\left|X\left(\mathbf{t}\right)-X\left(\mathbf{s}\right)\right|^{2}\right\\}\leq
C_{1}\left\|\mathbf{t}-\mathbf{s}\right\|^{2\kappa}_{\mathbbm{R}^{d}},$

for all $\mathbf{t},\mathbf{s}\in\mathcal{I}$.

  2. (A4)

The loss function $\rho:\mathbbm{R}\to\mathbbm{R}_{+}$ is convex and satisfies
a Lipschitz condition, i.e., there exists a $C_{2}>0$ such that

$\displaystyle|\rho(x)-\rho(y)|\leq
C_{2}|x-y|,\quad\forall(x,y)\in\mathbbm{R}^{2}.$

  3. (A5)

The function $g(t):=\mathbb{E}\\{\rho(\epsilon_{1}+t)\\}$ is uniquely
minimized at $t=0$ and is twice differentiable with a uniformly bounded second
derivative.

  4. (A6)

There exists an $\varepsilon>0$ such that for all $|t|\leq\varepsilon$,
$g(t)-g(0)\geq\varepsilon\,t^{2}$.

Assumption (A3) concerns the distribution and the sample paths of the
functional covariate $X$. We require that $X$ is a light-tailed but not
necessarily smooth process. The requirement of the existence of a squared
exponential moment of $\|X\|_{\infty}$ is clearly satisfied by bounded
processes but it is also satisfied by the vast majority of Gaussian processes
(see, e.g. Adler and Taylor,, 2007, Theorem 2.1.2). Thus, (A3) is a
considerable generalization of the boundedness requirement of (Boente et al.,,
2020; Kalogridis and Van Aelst,, 2023) for their respective robust functional
regression estimators. The required Hölder-continuity in mean-square implies
that $X$ has a modification with almost surely continuous sample paths. Beyond
continuity, however, the criterion does not imply any higher order smoothness
of the sample paths. For example, if $X$ is the Wiener process on
$\mathcal{I}=(0,1)$, it is well-known that $\mathbb{E}\\{|X(t)|^{2}\\}=t$ and
$\mathbb{E}\\{X(t)X(s)\\}=\min(t,s)$ for all $t,s\in(0,1)$. Hence, (A3) is
satisfied with $\kappa=1/2$, although the Wiener process is well-known to
possess sample paths that are nowhere differentiable.

Assumption (A4) is standard in robust regression (see, e.g., van de Geer,,
2000) and is satisfied by many loss functions including the quantile loss. It
should be noted that (A4) is the only assumption that is placed directly on
the loss function, as (A5) and (A6) are in effect placed on the distribution
of the errors. Specifically, in (A5) we require that the expectation
$\mathbb{E}\\{\rho(\epsilon_{1}+t)\\}$ is uniquely minimized at zero and is
smooth as a function of $t$. The former is required for identifiability while
the latter for convenience in the proofs. It is also worth noting that by
requiring $\mathbb{E}\\{\rho(\epsilon_{1}+t)\\}$ to be smooth instead of
$\rho$ itself to be smooth, we are effectively allowing smoothness to be
traded between $\rho$ and the distribution of the error, thereby permitting
even non-smooth $\rho$-functions within our theoretical framework. It is also
possible, but more notationally cumbersome, to relax the i.i.d. assumption and
instead consider merely independent errors in (2). In this case, (A5) and (A6)
would need to hold uniformly for all $\epsilon_{i}$.

Both parts of (A5) are satisfied quite generally. For instance, for
$\rho(x)=|x|$, (A5) is satisfied if the errors have unique median at zero.
Furthermore, as for this $\rho$,

$\displaystyle\mathbb{E}\left\\{\rho\left(\epsilon_{1}+t\right)-\rho\left(\epsilon_{1}\right)\right\\}=2\int_{0}^{-t}F(x)\,\mathrm{d}\,x+t,$
(12)

it is easy to see that the second part of (A5) is satisfied if the
distribution function of the error $F(x)=\mathbbm{P}(\epsilon_{1}\leq x)$ is
differentiable with a bounded derivative (density) $f=F^{\prime}$. This
assumption is less restrictive than the corresponding assumption in Kato,
(2012) (see assumption (A5) there), as we do not require the density to be
differentiable. Similarly to (A5), assumption (A6) is satisfied quite
generally. Continuing with the $L_{1}$ loss, it is easy to see that for all
small $|t|$ expression (12) becomes

$\displaystyle
g(t)-g(0)=\mathbb{E}\left\\{\rho\left(\epsilon_{1}+t\right)-\rho\left(\epsilon_{1}\right)\right\\}=f(-t)t^{2}+o(t^{2}),$

so that (A6) is ensured if $f$ is strictly positive in a neighborhood about
zero. These observations generalize straightforwardly to quantile estimation
with $\tau\in(0,1)$.

With assumptions (A1)–(A6) in place we may establish a rate of convergence for
the minimizer $\widehat{\beta}_{n}$ in (6), whose existence is ensured by
Proposition 1. In the statement of Theorem 1 below we use
$\operatorname{diam}(A)$ to denote the diameter of a set
$A\subset\mathbbm{R}^{d}$, that is,
$\operatorname{diam}(A)=\sup_{\mathbf{x},\mathbf{y}\in
A}\|\mathbf{x}-\mathbf{y}\|_{\mathbbm{R}^{d}}$.

######  Theorem 1.

Assume that (A1)–(A6) hold, $2m>d+1$ and $\lambda$ satisfies

  1. $(i)$

$\lambda=\log^{\frac{3}{2}}(n)O_{\mathbb{P}}(n^{-\frac{2m}{2m+d}}+\max_{1\leq
j\leq p}\operatorname{diam}^{2\kappa}(A_{j}))$.

  2. $(ii)$

$\lambda^{-1}\log^{\frac{3}{2}}(n)\left(n^{-\frac{2m}{2m+d}}+\max_{1\leq j\leq
p}\operatorname{diam}^{2\kappa}(A_{j})\right)=O_{\mathbb{P}}(1)$.

Then, for $\widehat{\beta}_{n}$ the solution to (6) with
$\widehat{\sigma}_{n}=1$, we have

$\displaystyle\left\|\widehat{\beta}_{n}-\beta_{0}\right\|_{n,p}=O_{\mathbb{P}}\left(\log^{2}(n)\left\\{n^{-\frac{m}{2m+d}}+\max_{1\leq
j\leq
p}\operatorname{diam}^{\kappa}(A_{j})\right\\}\right)\quad\text{and}\quad
J_{m}(\widehat{\beta}_{n})=O_{\mathbb{P}}(1).$

The side condition of the theorem $2m>d+1$ is a strengthening of our previous
condition $2m>d$, but likewise, it is not too restrictive. For instance, for
both $d=1$ and $d=2$ it is fulfilled by the popular choice $m=2$ corresponding
to cubic thin-plate splines. The limit conditions on $\lambda$ are more
complex than those for the smoothing parameters in Crambes et al., (2009) and
Kato, (2012), but these conditions ensure that our results remain valid even
if $\lambda$ is random. For deterministic $\lambda$ conditions $(i)$ and
$(ii)$ simplify to

$\displaystyle\lambda\asymp\log^{\frac{3}{2}}(n)\left(n^{-\frac{2m}{2m+d}}+\max_{1\leq
j\leq p}\operatorname{diam}^{2\kappa}(A_{j})\right).$

As noted in the introduction, the results of Crambes et al., (2009) and Kato,
(2012) do not hold if the smoothing parameters of their respective methods are
random. Hence, our theoretical results are novel in this respect. Theorem 1
shows that the rate of convergence for functional regression problems with
discretely sampled predictors depends on both the sample size and the grid
resolution, i.e., the discretization error. In particular, except for a
$\log^{2}(n)$-term, the rate of convergence is the optimal non-parametric rate
of convergence for $d$-dimensional data (Stone,, 1982) plus the Riemann
approximation error, which is governed by the regularity of the sample paths
of $X$. We see, in particular, that processes with $\kappa$ closer to $1$,
i.e., processes with smoother sample paths, will produce better rates of
convergence. It can be shown that the $\log^{2}(n)$-term in the result of
Theorem 2 does not appear if $X$ is assumed bounded, as other authors have
done (e.g., by Cardot et al.,, 2007; Boente et al.,, 2020). As noted earlier,
however, our assumptions are considerably more general, because they permit
both bounded and unbounded processes.

Another interesting observation that emerges from the rate of convergence
presented in Theorem 1 is the existence of a threshold value between the
sample size, $n$, and the discretization error, $\max_{j\leq
p}\operatorname{diam}^{\kappa}(A_{j})$, that determines the leading term in
the asymptotic error of the estimator, namely,
$T_{n,p}=n^{m/(2m+d)}\max_{j\leq p}\operatorname{diam}^{\kappa}(A_{j})$. For
$T_{n,p}<1$, that is, for relatively fast decay of the discretization error,
the asymptotic error, $\|\widehat{\beta}_{n}-\beta_{0}\|_{n,p}$, is solely
determined by the sample size through $n^{-m/(2m+d)}$ whereas in the case
$T_{n,p}\geq 1$ it is the discretization error, $\max_{j\leq
p}\operatorname{diam}^{\kappa}(A_{j})$, that determines up to a logarithmic
factor the rate of convergence. This finding provides support for the
intuitive idea that as long as the discretization grid is sufficiently dense,
then the discretization error becomes negligible in the limit. However,
Theorem 1 also demonstrates that the discretization error cannot be ignored
otherwise and indeed dominates the asymptotic error whenever the
discretization grid is relatively sparse.

The dimension $d$ plays an important role in the foregoing discussion, as, for
larger $d$, one needs a larger sample size relative to the discretization
error in order to obtain the optimal non-parametric rate of convergence
$n^{-m/(2m+d)}$. This may be viewed as a manifestation of the curse of
dimensionality that is often encountered in non-parametric estimation. For
$d=1$ and equispaced points $\\{t_{j}\\}_{j=1}^{p}$ within
$\mathcal{I}=(0,1)$, the result of Theorem 1 simplifies to

$\displaystyle\left\|\widehat{\beta}_{n}-\beta_{0}\right\|_{n,p}=O_{\mathbb{P}}\left(\log^{2}(n)\left\\{n^{-\frac{m}{2m+1}}+p^{-\kappa}\right\\}\right)\quad\text{and}\quad
J_{m}(\widehat{\beta}_{n})=O_{\mathbb{P}}(1),$

which, except for the $\log^{2}(n)$-term, is the same rate of convergence
obtained by Crambes et al., (2009) for their least-squares estimator without
any tail conditions on the eigenvalues of the covariance operator, that is,
with $q=0$ in the notation of those authors. It should be noted, nevertheless,
that our rate of convergence holds much more broadly, that is, for a much
greater collection of loss functions as well as for higher dimensions.

Theorem 1 establishes not only a rate of convergence with respect to the semi-
norm $\|\cdot\|_{n,p}$ but also the boundedness in probability of the penalty
functional $J_{m}(\widehat{\beta}_{n})$ evaluated at $\widehat{\beta}_{n}$. As
this penalty functional consists of two terms involving both
$\sum_{j=1}^{p}|\widehat{\beta}_{n}(\mathbf{t}_{j})|^{2}\mu(A_{j})$ and the
integrals of the squared partial derivatives,
$I_{m}^{2}(\widehat{\beta}_{n})$, the boundedness of
$J_{m}(\widehat{\beta}_{n})$ is essential in extending this rate of
convergence to a rate of convergence in the popular $\|\cdot\|_{n}$-norm given
in (10). Under the same assumptions as in Theorem 1, we obtain the following
useful corollary.

######  Corollary 2.

Suppose that the conditions of Theorem 1 hold. Then,
$\sup_{\mathbf{t}\in\mathcal{I}}|\widehat{\beta}_{n}(\mathbf{t})|=O_{\mathbb{P}}(1)$
and

$\displaystyle\left\|\widehat{\beta}_{n}-\beta_{0}\right\|_{n}=O_{\mathbb{P}}\left(\log^{2}(n)\left\\{n^{-\frac{m}{2m+d}}+\max_{1\leq
j\leq p}\operatorname{diam}^{\kappa}(A_{j})\right\\}\right).$

Much like Theorem 1, Corollary 2 reveals that unless
$n^{-m/(2m+d)}>\max_{1\leq j\leq p}\operatorname{diam}^{\kappa}(A_{j})$ or,
equivalently, $T_{n,p}<1$, the average squared error of prediction depends on
both the sample size and the discretization error and the latter cannot be
ignored in one’s analysis.

###  3.2 Thin-Plate Splines with Scale Estimation

We now turn to the theoretical investigation of thin-plate estimators
requiring an auxiliary scale estimate, such as the estimators based on the
logistic or Huber losses. We require the following assumptions in order to
establish the equivalent result of Theorem 1 in this setting.

  1. (B4)

The loss function $\rho:\mathbbm{R}\to\mathbbm{R}_{+}$ is convex with bounded
derivative $\psi$ that satisfies the following tail condition: for every
$\varepsilon_{1}>0$ there exists $D_{\varepsilon_{1}}>0$ such that

$\displaystyle\sup_{x\in\mathbbm{R}}\left|\psi(tx)-\psi(sx)\right|\leq
D_{\varepsilon_{1}}|t-s|,$

for any $t\geq\varepsilon_{1}$, $s\geq\varepsilon_{1}$.

  2. (B5)

There exists a $\sigma_{0}\in(0,\infty)$ such that
$\widehat{\sigma}_{n}\xrightarrow{\mathbb{P}}\sigma_{0}$.

  3. (B6)

There exists an $\varepsilon_{2}>0$ such that for all
$\sigma\in(\sigma_{0}-\varepsilon_{2},\sigma_{0}+\varepsilon_{2})$,
$\mathbb{E}\\{\psi(\epsilon_{1}/\sigma)\\}=0$ and the function
$g_{\sigma}(t):=\mathbb{E}\\{\psi(\epsilon_{1}/\sigma+t)\\}$ is continuously
differentiable with a uniformly bounded derivative $g_{\sigma}^{\prime}$ and
$g_{\sigma}^{\prime}(0)>0$.

Assumption (B4) is slightly more restrictive than (A4), as a bounded
derivative of $\rho$ ensures that the Lipschitz condition in (A4) is
satisfied. The tail condition in (B4) is satisfied whenever
$\psi(x)=\rho^{\prime}(x)$ changes slowly in the tail. For differentiable
$\psi$-functions it is easy to show that this tail condition is satisfied
provided that $\sup_{x\in\mathbbm{R}}|x\psi^{\prime}(x)|<\infty$, but
differentiability is not required. For instance, (B4) is also satisfied for
the Huber $\psi$-function derived from (8) although $\psi$ in that case is
only piecewise differentiable. It is important to note that the limiting value
$\sigma_{0}$ in (B5) needs not be the standard deviation of the error, which
may not exist. Lastly, (B6) parallels (A5) and (A6) and ensures the Fisher
consistency of the estimators while also allowing for smoothness to be traded
between $\psi$ and $F$, the distribution function of the error. As a specific
example, consider again the Huber $\rho$-function of (8). Then, if $F$ is
symmetric about $0$ with density $f$, direct calculation shows that

$\displaystyle g_{\sigma}^{\prime}(0)=2F\left(\frac{k}{\sigma}\right)-1,$

for every $\sigma>0$. Thus, $g_{\sigma}^{\prime}(0)$ is strictly positive
under symmetric $F$ and $k>0$ so that (B6) is satisfied. We are now ready to
state the main result of this section.

######  Theorem 2.

Assume that (A1)–(A3) and (B4)–(B6) hold, $2m>d+1$ and
$\lambda=o_{\mathbb{P}}(1)$ in such a way that

  1. $(i)$

$\lambda=\log^{\frac{5}{2}}(n)O_{\mathbb{P}}(n^{-\frac{2m}{2m+d}}+\max_{1\leq
j\leq p}\operatorname{diam}^{2\kappa}(A_{j}))$.

  2. $(ii)$

$\lambda^{-1}\log^{\frac{5}{2}}(n)\left(n^{-\frac{2m}{2m+d}}+\max_{1\leq j\leq
p}\operatorname{diam}^{2\kappa}(A_{j})\right)=O_{\mathbb{P}}(1)$.

Then, for $\widehat{\beta}_{n}$ the solution to (6), we have

$\displaystyle\left\|\widehat{\beta}_{n}-\beta_{0}\right\|_{n,p}=O_{\mathbb{P}}\left(\log^{3}(n)\left\\{n^{-\frac{m}{2m+d}}+\max_{1\leq
j\leq
p}\operatorname{diam}^{\kappa}(A_{j})\right\\}\right)\quad\text{and}\quad
J_{m}(\widehat{\beta}_{n})=O_{\mathbb{P}}(1).$

Comparing Theorem 2 to Theorem 1 reveals that after the addition of the random
$\widehat{\sigma}_{n}$ to (6), the the power on the $\log(n)$-term has been
raised from two to three. This small change is a consequence of the possibly
slow rate of convergence of $\widehat{\sigma}_{n}$ to $\sigma_{0}$ and it can
be shown that it does not occur provided that
$n^{1/2}(\widehat{\sigma}_{n}-\sigma_{0})=O_{\mathbb{P}}(1)$. It is worth
noting, nevertheless, that such high a rate is extremely difficult to achieve
in inverse statistical problems, such as functional regression. By contrast,
the generality of (B5) allows for explicit checking for many interesting scale
estimates, as we now demonstrate by means of an example.

Consider the popular M-scale estimator based on the residuals of an initial
estimator not requiring standardization with a scale estimate, such as the
$L_{1}$ estimator with $\rho(x)=|x|$. In particular, let
$\\{r_{i}\\}_{i=1}^{n}$ denote the residuals of the initial estimator, i.e.,
$r_{i}=Y_{i}-\widehat{\alpha}_{n}-\sum_{j=1}^{p}X_{i}(\mathbf{t}_{j})\widehat{\beta}_{n}(\mathbf{t}_{j})\mu(A_{j})$,
and let $\chi:\mathbbm{R}\to(0,1)$ denote an even loss function. Then, the
M-scale estimate $\widehat{\sigma}_{n}$ is the solution of

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\chi\left(\frac{r_{i}}{\widehat{\sigma}_{n}}\right)=\frac{1}{2}.$

Observe that for $\chi(x)=x^{2}$, $\widehat{\sigma}_{n}$ is $\sqrt{2}$ times
the standard deviation of the $\\{r_{i}\\}_{i=1}^{n}$. We refer to Maronna et
al., (2019, Chapters 2 and 10) for an extensive discussion of such estimates.
Define the population value $\sigma_{0}$ as the solution to

$\displaystyle\mathbb{E}\left\\{\chi\left(\frac{\epsilon_{1}}{\sigma_{0}}\right)\right\\}=\frac{1}{2}.$

It can be shown that $\sigma_{0}$ is unique and positive provided that $\chi$
is continuous and strictly increasing in all $x>0$ for which $\chi(x)<1$
(recall that $\chi$ is even). Thus, for bounded $\chi$ no moments of the error
$\epsilon$ are required for $\sigma_{0}$ to be well-defined. Now, the
residuals $\\{r_{i}\\}_{i=1}^{n}$ approximate $\\{\epsilon_{i}\\}_{i=1}^{n}$
and it may be expected that
$\widehat{\sigma}_{n}\xrightarrow{\mathbb{P}}\sigma_{0}$ under the right
conditions. This is indeed true, as Theorem 3 below shows.

######  Theorem 3.

Suppose that the conditions of Theorem 1 hold for the initial thin-plate
estimates $(\widehat{\alpha}_{n},\widehat{\beta}_{n})$. Suppose also that
$\chi:\mathbbm{R}\to[0,1]$ is even, strictly increasing for all $x>0$ for
whuch $\chi(x)<1$ and has a bounded derivative $\chi^{\prime}$. Then,
$\widehat{\sigma}_{n}\xrightarrow{\mathbb{P}}\sigma_{0}$, i.e. condition (B5)
is satisfied.

A popular loss function fulfilling the requirements of Theorem 3 is the Tukey
bisquare given by $\chi(x)=\min(1-(1-x^{2})^{3},1)$, which we will also use in
our numerical and real data examples. We close this section by remarking that
the rate of convergence of thin-plate spline estimators with auxiliary scale
in the $\|\cdot\|_{n}$-semi-norm given in (10) is exactly the same as the rate
obtained in Theorem 2 for the $\|\cdot\|_{n,p}$-norm. The proof of this
follows from completely analogous arguments as those employed in the proof of
Corollary 2 and is thus omitted.

##  4 Practical Implementation

By Proposition 1 the solution to (6) may be found in the space of natural
thin-plate splines with knots at $\\{\mathbf{t}_{j}\\}_{j=1}^{p}$ and each $g$
in that space admits the representation (9). As a basis
$\left\\{\phi_{j}\right\\}_{j=1}^{M}$ for the space of polynomials of total
order less than $m$, we use the monomials $t_{1}^{m_{1}}\ldots t_{d}^{m_{d}}$
for all nonnegative integers $m_{1},\ldots,m_{d}$ satisfying $0\leq
m_{1}+\ldots+m_{d}\leq m-1$. That space is of dimension $M=\binom{m+d-1}{d}$.
The requirement of $g$ to satisfy $I_{m}(g)<\infty$ is equivalent to
$\mathbf{\Phi}^{\top}\boldsymbol{\gamma}=\mathbf{0}$ where
$\mathbf{\Phi}_{i,j}=\phi_{j}(\mathbf{t}_{i})$, $i=1,\dots,p$, $j=1,\dots,M$
(see, e.g. Green and Silverman,, 1994, Chapter 7). Moreover, for
$\boldsymbol{\gamma}=(\gamma_{1},\ldots,\gamma_{M})$ satisfying this
constraint, Green and Silverman, (1994) show that
$I_{m}^{2}(g)=\boldsymbol{\gamma}^{\top}\mathbf{\Omega}\boldsymbol{\gamma}$
with $\mathbf{\Omega}_{i,j}=\eta_{m,d}(\|\mathbf{t}_{i}-\mathbf{t}_{j}\|)$,
$i,j=1,\dots,p$. In order to automatically incorporate the constraint into the
problem, we write $\boldsymbol{\gamma}=\mathbf{Q}\boldsymbol{\xi}$ for a
$p\times(p-M)$ matrix $\mathbf{Q}$ whose columns span the null space of
$\boldsymbol{\Phi}^{\top}$ and $\boldsymbol{\xi}\in\mathbbm{R}^{p-M}$.

Denoting the matrix of the discretized processes by $\mathbf{X}$, i.e.,
$\mathbf{X}_{i,j}=\mathbf{X}_{i}(\mathbf{t}_{j})$, $i=1,\dots,n$,
$j=1,\dots,p$, and combining the above facts, we may deduce that to obtain
$(\widehat{\alpha}_{n},\widehat{\beta}_{n})$ in (6), it suffices to minimize

$\displaystyle
L_{n}\left(\boldsymbol{\theta}\right)=\frac{1}{n}\sum_{i=1}^{n}\rho\left(\frac{Y_{i}-\mathbf{Z}_{i}^{\top}\boldsymbol{\theta}}{\widehat{\sigma}_{n}}\right)+\lambda\boldsymbol{\theta}^{\top}\mathbf{H}\boldsymbol{\theta},$
(13)

with
$\mathbf{Z}=[\mathbf{1},\mathbf{X}\mathbf{W}\mathbf{\Omega}\mathbf{Q},\mathbf{X}\mathbf{W}\boldsymbol{\Phi}]$,
$\mathbf{Z}_{i}$ its $i$th row,
$\boldsymbol{\theta}=(\alpha,\boldsymbol{\xi}^{\top},\boldsymbol{\delta}^{\top})^{\top}$,
and the $(p+1)\times(p+1)$ penalty matrix $\mathbf{H}$ given by

$\displaystyle\mathbf{H}=\begin{bmatrix}0&0&0\\\
0&\boldsymbol{Q}^{\top}\boldsymbol{\Omega}\mathbf{W}\boldsymbol{\Omega}\boldsymbol{Q}+\mathbf{Q}^{\top}\boldsymbol{\Omega}\mathbf{Q}&\mathbf{Q}^{\top}\boldsymbol{\Omega}\mathbf{W}\boldsymbol{\Phi}\\\
0&\boldsymbol{\Phi}^{\top}\mathbf{W}\boldsymbol{\Omega}\mathbf{Q}&\boldsymbol{\Phi}^{\top}\mathbf{W}\boldsymbol{\Phi}\end{bmatrix}.$

Here, $\mathbf{W}$ is a $p\times p$ diagonal matrix with
$\mu(A_{1}),\ldots,\mu(A_{p})$ in its diagonal. Thus, after suitable
simplification, the objective function in (6) may be reduced to the objective
function of a finite-dimensional penalized M-estimator, as given in (13).

Except for the least-squares case $\rho(x)=x^{2}$, the minimizer of (13)
cannot be obtained in closed form. Nevertheless, the minimizer can be
identified through the penalized variant of the well-known and efficient
iteratively reweighted least-squares (IRLS) algorithm, see Maronna et al.,
(2019, Chapter 5). The algorithm amounts to a penalized weighted least-squares
regression with weights at the $(k+1)$th step given by
$w_{i}=\rho^{\prime}(r_{i}^{k}/\widehat{\sigma}_{n})/(2\widehat{\sigma}_{n}r_{i}^{k})$
where $\rho^{\prime}$ is the derivative of $\rho$ and
$\left(r_{1}^{k},\dots,r_{n}^{k}\right)^{\top}=\mathbf{r}^{k}$ denote the
residuals from the $k$th step of the algorithm, i.e.,
$\mathbf{r}^{k}=\mathbf{Y}-\mathbf{Z}\boldsymbol{\theta}^{k}$, for
$\boldsymbol{\theta}^{k}$ the value of $\boldsymbol{\theta}$ in the $k$th step
of IRLS. The IRLS algorithm has the remarkable property of monotonically
decreasing the objective function until convergence to the minimizer is
reached, which for convex loss functions is guaranteed irrespective of the
starting point (Maronna et al.,, 2019, Chapter 9). The IRLS algorithm is not
directly applicable for the quantile loss function
$\rho_{\tau}(x)=2x(\tau-I(x<0))$, as $\rho_{\tau}$ is not differentiable at
zero. In this case, we replace $\rho_{\tau}$ with its smooth approximation

$\displaystyle\widetilde{\rho}_{\tau}(x)=\begin{cases}\rho_{\tau}(x)&|x|\geq\varepsilon\\\
2\tau x^{2}/\varepsilon&0\leq x<\varepsilon\\\
2(1-\tau)x^{2}/\varepsilon&-\varepsilon<x\leq 0,\end{cases}$

for some small $\varepsilon<0$. The IRLS algorithm is applicable for
$\widetilde{\rho}_{\tau}$ and, in our experience, converges fast to the
minimizer of (13).

In order to select the penalty parameter $\lambda$ from the data we rely on
the minimization of a robust scale of an approximation of the leave one out
residuals. Specifically, let
$\mathbf{r}_{-}(\lambda)=(r_{-1}(\lambda),\ldots,r_{-n}(\lambda))$ denote the
leave-one-out residuals as obtained with $\lambda$ as the smoothing parameter.
For the least-squares thin-plate estimator, $\mathbf{r}_{-}(\lambda)$ may be
obtained explicitly using the “leaving-out-one” lemma of (Wahba,, 1990,
Chapter 4):

$\displaystyle
r_{-i}(\lambda)=\frac{r_{i}(\lambda)}{1-\operatorname{diag}\left(\mathbf{H}(\lambda)\right)_{i}},\quad(i=1,\ldots,n),$

where $\\{r_{i}(\lambda)\\}_{i=1}^{n}$ denote the residuals and
$\mathbf{H}(\lambda)$ the “hat”-matrix, i.e.,
$\mathbf{H}(\lambda)=\mathbf{Z}\left(\mathbf{Z}^{\top}\mathbf{Z}+\lambda\mathbf{H}\right)^{-1}\mathbf{Z}^{\top}$
in our notation. For general thin-plate spline estimators a good approximation
of $\mathbf{r}_{-}(\lambda)$ may be obtained from the last step of the IRLS
algorithm, viz,

$\displaystyle
r_{-i}(\lambda)\approx\frac{r_{i}(\lambda)}{1-\operatorname{diag}(\mathbf{H}_{\mathbf{W}}(\lambda))_{i}},\quad(i=1,\ldots,n),$

where now $\mathbf{H}_{\mathbf{W}}(\lambda)$ is the weighted hat-matrix with
weighting matrix $\mathbf{W}=\operatorname{diag}(w_{1},\ldots,w_{n})$, i.e.,
$\mathbf{H}_{\mathbf{W}}(\lambda)=\mathbf{Z}\left(\mathbf{Z}^{\top}\mathbf{W}\mathbf{Z}+\lambda\mathbf{H}\right)^{-1}\mathbf{Z}^{\top}\mathbf{W}$.

With either the exact or approximate leave-one-out residuals, we propose to
select the value of $\lambda$ minimizing

$\displaystyle\operatorname{RCV}(\lambda)=\tau(\mathbf{r}_{-}(\lambda))^{2},$

where $\tau$ is the robust and efficient $\tau$-scale introduced by Yohai and
Zamar, (1988), which may be loosely interpreted as a weighted standard
deviation. This criterion has also been used by Maronna and Yohai, (2013) and
is a robustification of the classical ordinary cross validation criterion,
which is based on the mean of the squared leave-one-out residuals.
Implementations and illustrative examples of the least-squares, $L_{1}$, Huber
and logistic thin-plate spline estimators are provided in the R-package
RobustSpline (Nagy,, 2023) accompanying this paper.

##  5 Finite-Sample Performance

###  5.1 Numerical Experiments

In our numerical experiments we are interested in the sensitivity of the
proposed family of estimators to the grid size $p$, the noise-to-signal ratio
as well as the distributions of the functional covariate $X$ and the error
$\epsilon$. The thin-plate spline estimators that we consider for this
simulation study are:

  * •

The least-squares estimator with $\rho(x)=x^{2}$ abbreviated as S.

  * •

The $L_{1}$ estimator with $\rho(x)=|x|$ abbreviated as A.

  * •

The Huber estimator with $\rho(x)$ given in (8) and abbreviated as H.

  * •

The logistic estimator with $\rho(x)=2x+4\log(1+e^{-x})$ abbreviated as L.

For the Huber estimator in (8) we use $k=1.345$, which ensures $95\%$
efficiency in the Gaussian location model (see, e.g., Maronna et al.,, 2019,
Chapter 2). In order to cut down on computing times, in all the simulation
settings we calculate the scale $\widehat{\sigma}_{n}$ for the Huber and
logistic estimators from the residuals of initial undersmoothed $L_{1}$
estimator instead of a properly smoothed $L_{1}$ estimator that would be
computationally more intensive. As we shall see, this simplification
drastically reduces the computation times while not affecting the performance
of the Huber and logistic estimators. We consider one-dimensional functional
data in this study and present an example of two-dimensional functional data
in the following section. For all estimators considered we use $m=2$
corresponding to the popular cubic thin-plate splines.

(a)

(b)

Figure 1: $100$ curves generated with standard Gaussian and
$t_{2}$-distributed $\\{W_{ij}\\}_{i,j=1}^{100,50}$ on the left and right
panels, respectively.

We generate the functional covariates $\\{X_{i}\\}_{i=1}^{n}$ defined on
$\mathcal{I}=(0,1)$ according to the truncated Karhunen-Loève decomposition

$\displaystyle
X_{i}(t)=\sqrt{2}\sum_{j=1}^{50}W_{ij}\frac{\sin\left(\left(j-1/2\right)\pi
t\right)}{\left(j-1/2\right)\pi t},\quad(i=1,\ldots,n),$

where $\\{W_{ij}\\}_{i,j=1}^{n,50}$ are independent mean-zero random
variables. The response variables $\\{Y_{i}\\}_{i=1}^{n}$ are subsequently
generated from

$\displaystyle
Y_{i}=\frac{1}{100}\sum_{k=1}^{100}X_{i}(t_{k})\beta_{0}(t_{k})+\sigma\epsilon_{i},\quad(i=1,\ldots,n),$

where $\\{t_{k}\\}_{k=1}^{100}$ are equidistant within $(0,1)$,
$\beta_{0}(t)=-\sin(5\,t/1.2)/0.5-1$ is the coefficient function,
$\\{\epsilon_{i}\\}_{i=1}^{n}$ are random errors and $\sigma>0$ a constant
that is chosen so as to achieve a given noise-to-signal-ratio (NSR). The NSR
that we consider herein are $\\{0.1,0.2\\}$ corresponding to moderate and high
levels of noise in the data. We have also considered other coefficient
functions but found similar results so that we only report the results for
this $\beta_{0}$ here.

We test the estimators in the following three scenarios involving different
distributions of $\\{W_{ij}\\}_{i,j=1}^{n,50}$ and
$\\{\epsilon_{i}\\}_{i=1}^{n}$. The quantities $W_{ij}$ and $\epsilon_{i}$ are
always independent.

  1. Model 1

$\\{W_{ij}\\}_{i,j=1}^{n,50}$ and $\\{\epsilon_{i}\\}_{i=1}^{n}$ follow
standard Gaussian distributions.

  2. Model 2

$\\{W_{ij}\\}_{i,j=1}^{n,50}$ follow a standard Gaussian distribution and
$\\{\epsilon_{i}\\}_{i=1}^{n}$ follow a t-distribution with two degrees of
freedom.

  3. Model 3

$\\{W_{ij}\\}_{i,j=1}^{n,50}$ follow a $t$-distribution with two degrees of
freedom and $\\{\epsilon_{i}\\}_{i=1}^{n}$ follow the standard Gaussian
distribution.

The above scenarios reflect the settings of regular data (Model 1) and data
contaminated with vertical (Model 2) and leverage (Model 3) outliers,
respectively. For a better appreciation of the effect of the distribution of
the scores $\\{W_{ij}\\}_{i,j=1}^{n,50}$ on the shapes of the curves, Figure 1
plots 100 representative curves under Model 1 and Model 3 on the left and
right panels, respectively. It may be seen that switching from the Gaussian to
a t-distribution leads to several types of functional outliers among the
curves, e.g., shape or amplitude outliers, which in turn complicate the
estimation of $\beta_{0}$.

Table 1: Means and standard errors (in brackets) of the prediction/estimation
errors in Models 1–3 based on 1000 replications with $n=300$.

|  |  | $\operatorname{NSR}=0.1$ | $\operatorname{NSR}=0.2$
---|---|---|---|---
|  |  | $p=25$ | $p=50$ | $p=100$ | $p=25$ | $p=50$ | $p=100$
Model 1 |  $\operatorname{SPE}$ | S | 0.202 (0.004) | 0.045 (0.001) | 0.018 (0.000) | 0.258 (0.005) | 0.091 (0.001) | 0.068 (0.002)
A | 0.211 (0.004) | 0.052 (0.001) | 0.026 (0.001) | 0.281 (0.005) | 0.113 (0.002) | 0.090 (0.002)
H | 0.203 (0.004) | 0.046 (0.001) | 0.019 (0.000) | 0.260 (0.005) | 0.093 (0.002) | 0.070 (0.002)
L | 0.202 (0.004) | 0.045 (0.001) | 0.019 (0.000) | 0.260 (0.005) | 0.092 (0.002) | 0.070 (0.002)
3-9 |  $\operatorname{SEE}$ | S | 4.652 (0.247) | 3.414 (0.217) | 3.107 (0.200) | 11.451 (0.948) | 10.194 (0.886) | 11.312 (0.834)
| A | 4.910 (0.275) | 3.361 (0.238) | 3.403 (0.238) | 11.724 (0.863) | 10.016 (0.972) | 10.332 (0.798)
| H | 4.525 (0.226) | 3.535 (0.235) | 3.280 (0.217) | 11.659 (0.962) | 10.265 (0.818) | 11.430 (0.843)
| L | 4.533 (0.237) | 3.553 (0.233) | 3.313 (0.220) | 11.992 (0.960) | 10.358 (0.886) | 12.146 (0.923)
Model 2 |  $\operatorname{SPE}$ | S | 0.542 (0.118) | 0.591 (0.404) | 0.210 (0.035) | 0.824 (0.036) | 0.709 (0.045) | 0.610 (0.036)
A | 0.232 (0.004) | 0.069 (0.001) | 0.043 (0.001) | 0.362 (0.006) | 0.193 (0.004) | 0.162 (0.004)
H | 0.226 (0.004) | 0.067 (0.001) | 0.041 (0.001) | 0.347 (0.006) | 0.183 (0.004) | 0.152 (0.003)
L | 0.229 (0.004) | 0.070 (0.001) | 0.044 (0.001) | 0.359 (0.006) | 0.195 (0.004) | 0.162 (0.003)
3-9 |  $\operatorname{SEE}$ | S | 26.266 (2.830) | 20.125 (2.071) | 23.092 (2.132) | 88.298 (13.230) | 88.913 (10.978) | 78.883 (12.039)
| A | 13.615 (0.829) | 10.327 (0.540) | 9.963 (0.496) | 29.432 (1.601) | 35.711 (2.221) | 35.942 (2.173)
| H | 10.587 (0.753) | 8.601 (0.448) | 8.327 (0.477) | 25.757 (1.742) | 29.599 (1.924) | 28.171 (1.771)
| L | 10.097 (0.538) | 8.603 (0.457) | 8.128 (0.474) | 25.784 (1.603) | 30.830 (2.021) | 27.636 (1.705)
Model 3 |  $\operatorname{SPE}$ | S | 1.727 (0.047) | 0.520 (0.068) | 0.191 (0.014) | 2.638 (0.104) | 1.911 (0.726) | 0.746 (0.046)
A | 2.699 (0.239) | 0.679 (0.082) | 0.238 (0.013) | 3.352 (0.190) | 2.154 (0.670) | 0.933 (0.056)
H | 2.418 (0.201) | 0.564 (0.070) | 0.195 (0.014) | 2.816 (0.111) | 1.931 (0.733) | 0.769 (0.049)
L | 2.395 (0.200) | 0.556 (0.069) | 0.193 (0.014) | 2.777 (0.107) | 1.935 (0.733) | 0.754 (0.047)
3-9 |  $\operatorname{SEE}$ | S | 32.467 (1.898) | 8.766 (0.685) | 4.168 (0.473) | 44.604 (3.095) | 16.257 (1.326) | 16.348 (1.967)
| A | 10.076 (0.780) | 5.884 (0.627) | 3.617 (0.274) | 16.881 (1.386) | 10.832 (0.933) | 13.711 (1.629)
| H | 12.434 (0.806) | 5.774 (0.488) | 4.349 (0.569) | 24.784 (1.759) | 12.262 (0.962) | 16.448 (2.006)
| L | 12.631 (0.933) | 5.942 (0.509) | 4.227 (0.516) | 23.292 (1.600) | 12.099 (0.921) | 15.733 (1.671)

In order to test our estimators in the incomplete data setting, in each
simulation run we randomly choose $p$ points $0<t_{k_{1}}<\ldots<t_{k_{p}}<1$
from $\\{t_{k}\\}_{k=1}^{100}$ and we compute the thin-plate spline estimators
of $\beta_{0}$ from the incomplete data
$(X_{i}(t_{k_{1}}),\ldots,X_{i}(t_{k_{p}}),Y_{i})$, $i=1,\ldots,n$. We
consider values of $p$ in $\\{25,50,100\\}$ with $p=100$ corresponding to
fully observed data. To assess the performance of the competing estimators we
rely on the squared prediction and estimation error, SPE and SEE,
respectively, given by

$\displaystyle\operatorname{SPE}=\frac{1}{n}\sum_{i=1}^{n}\left|\frac{1}{100}\sum_{j=1}^{100}X_{i}(t_{k})\beta_{0}(t_{k})-\widehat{\alpha}_{n}-\sum_{j=2}^{p}X_{i}(t_{k_{j}})\widehat{\beta}_{n}(t_{k_{j}})(t_{k_{j}}-t_{k_{j-1}})\right|^{2},$

and

$\displaystyle\operatorname{SEE}=\sum_{j=2}^{p}\left|\beta_{0}(t_{k_{j}})-\widehat{\beta}_{n}(t_{k_{j}})\right|^{2}\left(t_{k_{j}}-t_{k_{j-1}}\right).$

Table 1 presents the average SPEs and SEEs and their standard errors for the
four competing estimators based on $1000$ replications with $n=300$.

There are several interesting observations emerging from Table 1. The most
notable is the rapid deterioration of the performance of the least-squares
estimator upon deviation from the ideal model conditions. Indeed, while under
Model 1 the least-squares estimator marginally outperforms its competitors,
under the heavy-tailed $t_{2}$-distributed errors in Model 2 its lead quickly
evaporates and the estimator ends up widely outperformed by all other
estimators. It is quite remarkable that these robust estimators maintain a
relatively stable performance under Model 1 and Model 2 with respect to
prediction. Their performance with respect to estimation deteriorates under
Model 2 but clearly not to the same extent as the performance of the least-
squares estimator. The Huber and logistic estimators perform quite comparably
in both situations and prove to be more efficient than the $L_{1}$ estimator
under light-tailed errors. To illustrate the differences in performance
between the least-squares and robust estimators, Figure 2 presents the least-
squares and Huber estimates for $\beta_{0}$ under Model 1 and Model 2 on the
top and bottom rows respectively. These plots show that while the Huber
estimates remain stable under contamination, the least-squares estimates can
become erratic leading to a noticeable deterioration in performance.

It is interesting to note that all estimators appear vulnerable to
contamination in the predictor space as given in Model 3, although not to the
same extent. The fact that estimators based on convex $\rho$-functions are
vulnerable to this type of contamination is well-known in the robustness
literature, see Maronna et al., (2019, Chapter 4). No estimator convincingly
outperforms the others with respect to prediction in this model, but the
robust estimators perform significantly better than the least-squares
estimator with respect to estimation. Among the robust estimators, the $L_{1}$
estimator offers the most protection against outlying observations in the
predictor space albeit not by a large margin. However, since the
$L_{1}$-estimator can perform considerably worse than the Huber and logistic
estimators under regular data, its use is, in our opinion, warranted only in
cases when one suspects heavy contamination within the data.

It is also worth noting that higher level of noise and larger $p$ affect all
estimators in the same way. In particular, higher NSR makes both estimation
and prediction harder for all estimators. The effect of a larger $p$, that is,
more completely observed functional data, is overall positive and particularly
noticeable for moderate NSR. For high NSR, the situation is less clear-cut and
a larger $p$ in this setting may even lead to deterioration in the performance
of the estimators. We attribute this difference to the fact that for a larger
$p$ more coefficients in (9) need to be estimated from the data. If the level
of noise is high, these coefficients cannot be accurately estimated and as a
result the performance of the estimators need not necessarily improve.

(a)  (b)

(c)  (d)

Figure 2: 1000 Least-squares (left) and Huber (right) estimates for Model 1
(top panels) and Model 2 (bottom panels) under $\operatorname{NSR}=0.1$ and
$p=100$. The lines ( , ) depict the true coefficient function $\beta_{0}$ and
the first $5$ estimated functions.

In summary, the results of our simulation study point towards the advantages
of thin-plate spline estimators based on more slowly increasing
$\rho$-functions relative to $\rho(x)=x^{2}$. With an appropriate choice of
$\rho$, these estimators are very versatile in the sense that they are
efficient in clean data and considerably safer in the presence of atypical
observations than the least-squares thin-plate spline estimator.

###  5.2 Computing Times

Table 2 provides the average computing times and standard errors for each
estimator under all the settings of our numerical study. From these numbers it
may be seen that, while the least-squares estimator is the fastest to compute,
the Huber and logistic estimators are not too far behind. In fact, for Model 1
and Model 2, for the Huber and logistic estimators in the vast majority of
cases, the IRLS algorithm converges in at most four to five iterations. Thus,
to obtain each one of these estimators it suffices to solve up to five least-
square problems. By contrast, owing to its less smooth $\rho$-function, the
$L_{1}$ estimator is much more computationally intensive.

Table 2: Means and standard errors (in brackets) of the computation times in
seconds.

|  | $\operatorname{NSR}=0.1$ | $\operatorname{NSR}=0.2$
---|---|---|---
|  | $p=25$ | $p=50$ | $p=100$ | $p=25$ | $p=50$ | $p=100$
Model 1 | S | 0.156 (0.015) | 0.337 (0.031) | 0.982 (0.101) | 0.156 (0.014) | 0.340 (0.036) | 0.986 (0.106)
A | 2.404 (1.115) | 8.083 (3.470) | 35.795 (14.486) | 2.632 (1.164) | 8.473 (3.543) | 36.316 (13.402)
H | 0.489 (0.712) | 0.772 (0.633) | 2.675 (1.520) | 0.752 (1.014) | 1.113 (1.673) | 2.588 (1.124)
L | 0.489 (0.757) | 0.732 (0.772) | 2.357 (0.766) | 0.788 (1.107) | 1.048 (1.520) | 2.328 (0.949)
1-8  Model 2 | S | 0.156 (0.014) | 0.348 (0.046) | 0.983 (0.106) | 0.157 (0.018) | 0.355 (0.047) | 0.963 (0.110)
A | 2.451 (1.192) | 8.333 (3.550) | 36.243 (13.968) | 2.733 (1.214) | 8.914 (3.621) | 35.693 (13.150)
H | 0.625 (0.877) | 0.949 (0.999) | 2.876 (1.174) | 1.023 (1.173) | 1.553 (2.228) | 2.719 (0.863)
L | 0.694 (1.017) | 0.928 (1.077) | 2.636 (0.701) | 1.188 (1.367) | 1.573 (2.319) | 2.534 (0.664)
1-8  Model 3 | S | 0.158 (0.015) | 0.352 (0.040) | 0.979 (0.111) | 0.154 (0.018) | 0.354 (0.049) | 0.937 (0.104)
A | 4.565 (1.612) | 13.268 (6.006) | 49.798 (22.925) | 4.614 (1.267) | 14.792 (5.286) | 48.045 (19.948)
H | 3.154 (1.680) | 8.140 (6.355) | 3.898 (3.896) | 3.215 (1.279) | 10.592 (5.431) | 5.475 (6.696)
L | 3.737 (1.868) | 9.104 (6.584) | 4.217 (6.063) | 3.785 (1.388) | 11.699 (5.262) | 9.147 (12.140)

For Model 3, i.e., contamination in the predictor space, the Huber and
logistic estimators require more iterations of the IRLS algorithm, but these
estimators remain still considerably easier to compute than the $L_{1}$
estimator. Overall, Table 2 indicates that the Huber and logistic estimators
are computationally feasible alternatives to the least-squares estimator
particularly in light of their increased resistance to atypical observations.

(a)

(b)

Figure 3: Contours of the least-squares and Huber thin-plate spline
coefficient function estimates on the convex hull of the data on the left and
right panels, respectively. Lighter colors correspond to larger values of the
estimates. The symbols ( , , ) respectively indicate the positions of the
monitoring stations, the position of our reference station and the outlying
observations detected by the Huber estimates.

##  6 Application: Ozone Concentration in Midwestern United States

It is well-known that stratospheric ozone is beneficial to living organisms,
as it protects them from ultraviolet sun radiation. High concentration of
ground ozone, on the other hand, has been linked with a variety of respiratory
problems, primarily for children and the elderly. Therefore, the ability to
predict high concentrations of ground ozone may be helpful. As an illustration
of the practical usefulness of the proposed family of estimators we analyze
ground level ozone concentration in south Chicago as a function of ozone
concentration in a broad area of the midwestern US. In particular, the data
for this analysis consists of 7-hour average measurements of ground level
ozone from 9am to 4pm in parts per billion (PPB) from $153$ stations during 89
days and are freely available as a part of the fields R-package (Nychka et
al.,, 2021) on CRAN.

To explain our approach in detail, let $\\{Y_{i}\\}_{i=1}^{89}$ denote the
ground level ozone concentration in south Chicago on each one of the 89 days
of measurement and let $\\{X_{i}(t,s)\\}_{i=1}^{89}$ denote the ground level
ozone concentration on each one of these days as a function of longitude and
latitude, denoted by $t$ and $s$, respectively. We consider the functional
linear model

$\displaystyle
Y_{i}=\alpha_{0}+\int_{37}^{45}\int_{-93}^{-83}X_{i}(t,s)\beta_{0}(t,s)\,\mathrm{d}\,t\,\mathrm{d}\,s+\epsilon_{i},\quad(i=1,\ldots,89),$

for some unknown
$(\alpha_{0},\beta_{0})\in\mathbbm{R}\times\mathcal{H}^{m}(\mathbbm{R}^{2})$.
Here, the boundaries of integration result from the geographical positions of
the measuring stations within the US. Of course, we are immediately faced with
an incomplete data problem, as each $X_{i}$ is only observed at $152$ stations
rather than in its entirety. Nevertheless, our methodology is applicable and
in order to estimate $(\alpha_{0},\beta_{0})$ we compute both the least-
squares and Huber thin-plate spline estimators, which yield the contour plots
in Figure 3. Note that, since the ozone concentration reported by the
monitoring station in south Chicago is used as the response variable, the
geographic location of this station is excluded from the dataset.

(a)

(b)

Figure 4: Contours of the least-squares and Huber thin-plate spline
coefficient function estimates on the convex hull of the data on the left and
right panels, respectively, after the outliers have been removed. Lighter
colors correspond to larger values of the estimates. The symbols ( , )
respectively indicate the positions of the monitoring stations and the
position of our reference station.

It is natural to hypothesize that surface ozone concentration in one place
would be positively correlated with surface ozone concentration in nearby
places and that this correlation would wane as the distance increases. In
other words, it is primarily ozone concentrations immediately in the vicinity
of south Chicago that should be the most relevant. However, an examination of
the contour plots in Figure 4 shows that while this is mostly true for the
Huber estimates, the least-squares estimates indicate that the ozone
concentrations in a large area in the north-west and south-east of St. Louis
are very important predictors for the ozone concentration in south Chicago.
This may seem puzzling, as the distance between these two places is more than
450 kilometers and so we may wonder whether the least-squares estimates are
being distorted by atypical observations.

In order to detect atypical and influential observations we may use the
residuals of the more resistant Huber estimates. These estimates are much less
drawn to atypical observations and therefore atypical observations show as
large residuals. A classical detection rule in that respect involves examining
the standardized residuals $r_{i}/\widehat{\sigma}_{n}$ and classifying as
outliers observations whose absolute standardized residuals exceed $2.6$
(Rousseeuw and Leroy,, 1987, Chapter 6). In our case, this results in $14$
observations being classified as outliers; these observations are depicted
with red triangles in the right panel of Figure 3.

It is interesting to observe that many of the outlying observations are
located in or around the areas on which the least-squares and Huber estimates
differ. As a sensitivity check, we remove the outlying observations from the
data and recompute the estimators. Doing so now yields the contours in Figure
4. These contours reveal that while the Huber estimates have barely changed
after the removal of the outlying observations, the least-squares estimates
have undergone significant adjustment and are now very similar to the Huber
estimates. Moreover, in the revised contour plot, the area around St. Louis
seems to be much less influential in the determination of ground level ozone
in Chicago and thus the estimates are now in line with prior expectations.
This example demonstrates the sensitivity of least-squares estimators as well
as the versatility of robust estimators, such as Huber estimators, which may
be used to good effect both in the presence and absence of outliers and can
even help detect outlying observations.

##  7 Concluding Remarks

The theoretical and practical results of this paper provide justification for
the use of a broad class of functional linear regression estimators in the
ubiquitous setting of discretely sampled, possibly high dimensional,
functional data. There are several research directions worth exploring from
here, the most notable of which would be the generalization of our methodology
to the setting of generalized linear models (GLM) allowing for different types
of response variables, such as binary or count. The extension of our thin-
plate spline estimator to the GLM context as well as the theoretical
investigation of its properties in that setting is a promising topic for
future research.

Another interesting and practically useful direction of future research would
be the investigation of the theoretical properties of robust model selection
procedures, such as the $\tau$-scale of the leave-one-out residuals employed
herein. To the best of our knowledge, there do not currently exist any
theoretical results with respect to robust model selection in the functional
data setting despite the fact that the vast majority of available estimation
procedures make use of robust model selection criteria. Thus, research in this
direction would fill a gap in the literature and provide practical guarantees
with respect to the use of robust model selection procedures.

## Acknowledgments

Most of the present research was carried out while I. Kalogridis was visiting
Charles University. I. Kalogridis gratefully acknowledges support from the
Research Foundation-Flanders (project 1221122N). The research of S. Nagy was
supported by Czech Science Foundation (project n. 23-05737S).

##  Appendix A Appendix: Proofs of our Theoretical Results

In our proofs, we occasionally use inequalities of the form $A\leq B$ for
$A,B$ random variables defined on the same probability space $\Omega$. These
inequalities are always to be interpreted in the almost sure sense, i.e.
$A(\omega)\leq B(\omega)$ for $\mathbb{P}$-almost all $\omega\in\Omega$.
Moreover, we use $\langle\cdot,\cdot\rangle_{p}$ denote the discrete inner
product, i.e., $\langle
f,g\rangle_{p}=\sum_{j=1}^{p}f(\mathbf{t}_{j})g(\mathbf{t}_{j})\mu(A_{j})$ and
$\|\cdot\|_{p}$ to denote the associated norm, i.e., $\|f\|_{p}=\langle
f,f\rangle_{p}^{1/2}$. For the proofs of Theorem 1, Corollary 2, Theorem 2 and
Theorem 3 we assume that $\mu(A_{j})=\mu(\mathcal{I})/p$, but this assumption
is not needed for the proofs of Proposition 1, Corollary 1 and Proposition 2.

###  A.1 Proof of Proposition 1

We may restrict attention to $\beta\in\mathcal{H}^{m}(\mathbbm{R}^{d})$, as
$\alpha$ is a finite-dimensional parameter and therefore the existence of a
minimizer $\widehat{\alpha}_{n}$ follows from the convexity of $\rho$. Our
proof consists of two steps related to the existence and characterization of
the minimizer $\widehat{\beta}_{n}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$.

_Step 1: Existence._ We aim to apply the existence theorem given by Gu, (2013,
Theorem 2.9) and for this we need to characterize the null space of the
penalty functional $J_{m}$. We claim that $J_{m}(g)=0\iff g\equiv 0$ so that
the null space of $J_{m}$ consists only of the zero element. To prove this we
shall show that $J_{m}(g)=0\iff g|_{\mathcal{I}}\equiv 0$ and $g$ is a
polynomial on $\mathbbm{R}^{d}$ of total order $m$. Necessity is obvious,
hence we focus on the sufficiency part, that is, we prove that (i)
$J_{m}(g)=0\implies g|_{\mathcal{I}}=0$, and (ii) $g$ is a polynomial on
$\mathbbm{R}^{d}$ of total order $m$. Implication (ii) is clear, as
$J_{m}(g)=0$ implies that $I_{m}(g)=0$, which in turn implies that $g$ is a
polynomial of total order $m$. Therefore, we need to show that, for every
$g\in\mathcal{H}^{m}\left(\mathbbm{R}^{d}\right)$, $J_{m}(g)=0$ implies that
$g|_{\mathcal{I}}=0$. To see this, recall that by (A1), $\mathcal{I}$ is a
bounded open set satisfying the uniform cone condition. One can thus consider
also the Sobolev space $\mathcal{H}^{m}\left(\mathcal{I}\right)$ of functions
$f\colon\mathcal{I}\to\mathbbm{R}$ defined analogously as
$\mathcal{H}^{m}\left(\mathbbm{R}^{d}\right)$, equipped with the norm
$\left\|\cdot\right\|_{\mathcal{H}^{m}\left(\mathcal{I}\right)}$ defined by

$\displaystyle\left\|g\right\|_{\mathcal{H}^{m}(\mathcal{I})}^{2}=\int_{\mathcal{I}}|g(\mathbf{t})|^{2}\,\mathrm{d}\,\mathbf{t}+\sum_{m_{1}+\ldots+m_{d}=m}\binom{m}{m_{1},\ldots,m_{d}}\int_{\mathcal{I}}\left(\frac{\partial^{m}g(\mathbf{t})}{\partial
t_{1}^{m_{1}}\ldots\partial
t_{d}^{m_{d}}}\right)^{2}\,\mathrm{d}\,\mathbf{t},$ (14)

where naturally $\mathbf{t}=(t_{1},\dots,t_{d})$. By Theorem 3.4 of Utreras,
(1988), there exist $C_{0}$ and $B_{1}$ depending only on $m,d,\mathcal{I}$
such that

$\displaystyle\int_{\mathcal{I}}|g(\mathbf{t})|^{2}d\mathbf{t}\leq\frac{C_{0}}{p}\sum_{j=1}^{p}|g(\mathbf{t}_{j})|^{2}+C_{0}B_{1}^{2m}\sum_{m_{1}+\ldots+m_{d}=m}\binom{m}{m_{1},\ldots,m_{d}}\int_{\mathcal{I}}\left(\frac{\partial^{m}g(\mathbf{t})}{\partial
t_{1}^{m_{1}}\ldots\partial
t_{d}^{m_{d}}}\right)^{2}\,\mathrm{d}\,\mathbf{t}.$

By definition of $J_{m}(g)$ and since $\mathcal{I}\subset\mathbbm{R}^{d}$, we
now find that

$\displaystyle\int_{\mathcal{I}}|g(\mathbf{t})|^{2}\,\mathrm{d}\,\mathbf{t}$
$\displaystyle\leq\frac{C_{0}}{p}\sum_{j=1}^{p}|g(\mathbf{t}_{j})|^{2}+C_{0}B_{1}^{2m}\sum_{m_{1}+\ldots+m_{d}=m}\binom{m}{m_{1},\ldots,m_{d}}\int_{\mathbbm{R}^{d}}\left(\frac{\partial^{m}g(\mathbf{t})}{\partial
t_{1}^{m_{1}}\ldots\partial t_{d}^{m_{d}}}\right)^{2}\,\mathrm{d}\,\mathbf{t}$
$\displaystyle=0,$

as $J_{m}(g)=0$ implies that $g(\mathbf{t}_{j})=0$ for $j=1,\ldots,p$, due to
the fact that, by construction, $\mu(A_{j})>0$. By hypothesis, $\mathcal{I}$
satisfies the uniform cone condition, so for $2m>d$ the Sobolev embedding
theorem (Adams and Fournier,, 2003, Theorem 6.3) implies the (compact)
embedding

$\displaystyle\mathcal{H}^{m}\left(\mathcal{I}\right)\to\mathcal{C}(\mathcal{I}),$

with the consequence that there exists a $c_{0}>0$ with the property that

$\displaystyle\sup_{\mathbf{t}\in\mathcal{I}}\left|g(\mathbf{t})\right|\leq
c_{0}\left\|g\right\|_{\mathcal{H}^{m}\left(\mathcal{I}\right)},$

for every $g\in\mathcal{H}^{m}(\mathcal{I})$. Clearly, for every
$g\in\mathcal{H}^{m}(\mathbbm{R}^{d})$,
$g|_{\mathcal{I}}\in\mathcal{H}^{m}(\mathcal{I})$, hence

$\displaystyle\sup_{\mathbf{t}\in\mathcal{I}}\left|g|_{\mathcal{I}}(\mathbf{t})\right|$
$\displaystyle\leq c_{0}\left\|g\right\|_{\mathcal{H}^{m}(\mathcal{I})}=0,$

proving the assertion. Thus, the null space of $J_{m}$ consists only of the
zero function so that the uniqueness condition of Theorem 2.9 in Gu, (2013) is
trivially satisfied. That theorem therefore applies and yields the existence
of a minimizer $\widehat{\beta}_{n}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$
concluding the first step in our proof.

_Step 2: Characterization._ We claim that under the conditions of the
proposition $\widehat{\beta}_{n}$ must necessarily be a natural thin-plate
spline with knots at the $\\{\mathbf{t}_{j}\\}_{j=1}^{p}$. To show this, we
apply Theorem 4 of Duchon, (1977) according to which there exists exactly one
natural thin-plate spline, $s$, taking prescribed values at the
$\\{\mathbf{t}_{j}\\}_{j=1}^{p}$ while minimizing $I_{m}^{2}(g)$. That is, for
every set of points $y_{1},\ldots,y_{p}\in\mathbbm{R}$,
$s:\mathbbm{R}^{d}\to\mathbbm{R}$ it solves the problem

$\displaystyle\min_{g\in\mathcal{H}^{m}\left(\mathbbm{R}^{d}\right)}I_{m}^{2}(g),\quad\text{subject
to}\quad g(\mathbf{t}_{j})=y_{j},\ (j=1,\ldots,p).$

By Step 1 of this proof, there exists a minimizer
$\widehat{\beta}_{n}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$. Denote
$\widehat{\boldsymbol{\beta}}_{n}=(\widehat{\beta}_{n}(\mathbf{t}_{1}),\ldots,\widehat{\beta}_{n}(\mathbf{t}_{p}))$.
Then, by Theorem 4 of Duchon, (1977), there exists a natural thin-plate
spline, $\widehat{s}_{n}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$, interpolating
$\widehat{\beta}_{n}$ at the $\\{\mathbf{t}_{j}\\}_{j=1}^{p}$ while
$I_{m}(\widehat{s}_{n})<I_{m}(\widehat{\beta}_{n})$. Because of these
properties, it follows that $\widehat{\beta}_{n}$ must itself be a natural
thin-plate spline. Indeed, let
$L_{n}(\widehat{\alpha}_{n},\widehat{\beta}_{n},\widehat{\sigma}_{n})$ denote
the minimal value of the objective function, that is,

$\displaystyle
L_{n}(\widehat{\alpha}_{n},\widehat{\beta}_{n},\widehat{\sigma}_{n})=\frac{1}{n}\sum_{i=1}^{n}\rho\left(\frac{Y_{i}-\widehat{\alpha}_{n}-\sum_{j=1}^{p}X_{i}(\mathbf{t}_{j})\widehat{\beta}_{n}(\mathbf{t}_{j})\mu(A_{j})}{\widehat{\sigma}_{n}}\right)+\lambda\left[\sum_{j=1}^{p}|\widehat{\beta}_{n}(\mathbf{t}_{j})|^{2}\mu(A_{j})+I_{m}^{2}\left(\widehat{\beta}_{n}\right)\right].$

If $\widehat{\beta}_{n}$ were not a natural thin-plate spline, then
$\widehat{s}_{n}$ would leave the first two terms of
$L_{n}(\widehat{\alpha}_{n},\widehat{\beta}_{n},\widehat{\sigma}_{n})$
unchanged, but, due to $I_{m}(\widehat{s}_{n})<I_{m}(\widehat{\beta}_{n})$, we
would have
$L_{n}(\widehat{\alpha}_{n},\widehat{s}_{n},\widehat{\sigma}_{n})<L_{n}(\widehat{\alpha}_{n},\widehat{\beta}_{n},\widehat{\sigma}_{n})$,
contradicting the fact that $\widehat{\beta}_{n}$ is a minimizer. The proof is
complete.

###  A.2 Proof of Corollary 1

Let $L_{n,\inf}$ denote the infimum of the objective function, that is,

$\displaystyle
L_{n,\inf}=\inf_{(\alpha,\beta)\in\mathbbm{R}\times\mathcal{H}^{m}(\mathbbm{R}^{d})}\left[\frac{1}{n}\sum_{i=1}^{n}\rho\left(\frac{Y_{i}-\alpha-\sum_{j=1}^{p}X_{i}(\mathbf{t}_{j})\beta(\mathbf{t}_{j})\mu(A_{j})}{\widehat{\sigma}_{n}}\right)+\lambda
J_{m}^{2}\left(\beta\right)\right].$

Note that $L_{n,\inf}$ is finite, as the objective function is bounded from
below by $0$. Define
$\mathcal{D}=\\{(\alpha,\beta)\in\mathbbm{R}\times\mathcal{H}^{m}(\mathbbm{R}^{d}):L_{n}(\alpha,\beta,\widehat{\sigma}_{n})=L_{n,\inf}\\}$.
Clearly, $\mathcal{D}=\bigcap_{k=1}^{\infty}\mathcal{D}_{k}$ with
$\mathcal{D}_{k}=\\{(\alpha,\beta)\in\mathbbm{R}\times\mathcal{H}^{m}(\mathbbm{R}^{d}):L_{n}(\alpha,\beta,\widehat{\sigma}_{n})\leq
L_{n,\inf}+k^{-1}\\}$. By the convexity of
$L_{n}(\alpha,\beta,\widehat{\sigma}_{n})$, each $\mathcal{D}_{k}$ is a convex
set. Hence $\mathcal{D}$, as the intersection of convex sets, is also convex.
By Proposition 1, it is also non-empty, as
$(\widehat{\alpha}_{n},\widehat{\beta}_{n})\in\mathcal{D}$. Suppose now that
there exist distinct $(\alpha_{1},\beta_{1})$ and $(\alpha_{2},\beta_{2})$ in
$\mathcal{D}$. By the convexity of $\mathcal{D}$ and the strict convexity of
$\rho(x)$, we would have

$\displaystyle
L_{n,\inf}=L_{n}(\alpha_{1}/2+\alpha_{2}/2,\beta_{1}/2+\beta_{2}/2,\widehat{\sigma}_{n})<L_{n}(\alpha_{1},\beta_{1},\widehat{\sigma}_{n})/2+L_{n}(\alpha_{2},\beta_{2},\widehat{\sigma}_{n})/2=L_{n,\inf},$

which is clearly impossible. It follows that $\mathcal{D}$ must consist of a
single element, as was to be shown.

###  A.3 Proof of Proposition 2

Define for simplicity
$\|\beta\|_{p}^{2}=p^{-1}\sum_{j=1}^{p}|\beta(t_{j})|^{2}$ so that

$\displaystyle\widetilde{J}_{m}^{2}(\beta)=\left\|\mathcal{P}\beta\right\|^{2}_{p}+\left\|\beta^{(m)}\right\|^{2}.$

By the properties of Hilbert projections (see, e.g., Conway,, 1990, Theorem
2.7), we have $\|\mathcal{P}\beta\|^{2}_{p}\leq\|\beta\|_{p}^{2}$ and from
this we can immediately see that $\widetilde{J}_{m}^{2}(\beta)\leq
J_{m}^{2}(\beta)$. Thus, Proposition 2 holds with $c_{1}=1$.

To prove the other inequality, use the Hilbert projection theorem (Conway,,
1990, Theorem 2.5) to write

$\displaystyle\beta=\mathcal{P}\beta+(\mathcal{I}-\mathcal{P})\beta.$

with $\mathcal{I}$ denoting the identity operator on $\mathcal{H}^{m}(0,1)$.
Clearly,

$\displaystyle J_{m}^{2}(\beta)$
$\displaystyle=\left\|\mathcal{P}\beta\right\|_{p}^{2}+\left\|\beta-\mathcal{P}\beta\right\|_{p}^{2}+2\langle\mathcal{P}\beta,\beta-\mathcal{P}\beta\rangle_{p}+\left\|\beta^{(m)}\right\|^{2}$
$\displaystyle\leq
2\left(\left\|\mathcal{P}\beta\right\|_{p}^{2}+\left\|\beta-\mathcal{P}\beta\right\|_{p}^{2}\right)+\left\|\beta^{(m)}\right\|^{2}.$

Now, again by properties of the projections (Conway,, 1990, Theorem 2.5), we
have $\|\beta-\mathcal{P}\beta\|_{p}^{2}\leq\|\beta-v\|_{p}^{2}$ for every
other polynomial $v$ of order $m$. Choose $v$ to be the Taylor polynomial
about zero of order $m$, viz,

$\displaystyle v(t)=\sum_{i=0}^{m-1}\beta^{(i)}(0)\frac{t^{i}}{i!},$

so that, using the integral form of the remainder for Taylor polynomials and
the Cauchy-Schwarz inequality,

$\displaystyle\left\|\beta-\mathcal{P}\beta\right\|_{p}^{2}\leq\left\|\beta-v\right\|_{p}^{2}=\frac{1}{p}\sum_{j=1}^{p}\left|\int_{0}^{1}\frac{(t_{j}-u)_{+}^{m-1}}{(m-1)!}\beta^{(m)}(u)\,\mathrm{d}\,u\right|^{2}\leq
C\left\|\beta^{(m)}\right\|^{2},$

for some $C>0$ depending only on $m$. It follows that

$\displaystyle
J_{m}^{2}(\beta)\leq(2+2C)\widetilde{J}_{m}^{2}(\beta)=c_{2}^{2}\widetilde{J}_{m}^{2}(\beta),$

for $c_{2}=\sqrt{2(1+C)}$, as was to be shown.

###  A.4 Auxiliary Results on Empirical Processes

Our next task is to derive the modulus of continuity of the empirical process
associated with our estimator, which we do in general terms. Therefore, the
notations used in this section are mostly independent of that used in the rest
of our paper.

Let $\mathcal{B}$ denote a subset of a (semi-)metric space
$(\Lambda,{\mathfrak{D}})$ endowed with a semi-metric ${\mathfrak{D}}$ and let
$N(\delta,\mathcal{B},{\mathfrak{D}})$ denote the number of balls required to
cover $\mathcal{B}$, that is,

$\displaystyle N(\delta,\mathcal{B},{\mathfrak{D}})=\min(N:\ \text{there
exist}\ \\{\beta_{1},\ldots,\beta_{N}\\}\subset\mathcal{B}\ \text{such that}\
\max_{\beta\in\mathcal{B}}\min_{j=1,\ldots,N}{\mathfrak{D}}(\beta,\beta_{j})\leq\delta).$

Further, let $H(\delta,\mathcal{B},{\mathfrak{D}})=\log
N(\delta,\mathcal{B},{\mathfrak{D}})$ denote the $\delta$-entropy of
$\mathcal{B}$. We assume that the semi-metric ${\mathfrak{D}}$ is of the form

$\displaystyle{\mathfrak{D}}^{2}=\frac{1}{n}\sum_{i=1}^{n}{\mathfrak{D}}_{i}^{2},$
(15)

where ${\mathfrak{D}}_{1},\ldots,{\mathfrak{D}}_{n}$ are also semi-metrics on
$\mathcal{B}$. Consider real-valued random variables
$U_{i,\beta},i=1,\ldots,n,\ \beta\in\mathcal{B}$. We assume that the processes
$\\{U_{i,\beta}:\beta\in\mathcal{B}\\}$, $i=1,\ldots,n$, are independent and
centered, and that

$\displaystyle\left|U_{i,\beta}-U_{i,\widetilde{\beta}}\right|\leq
V_{i}\,{\mathfrak{D}}_{i}(\beta,\widetilde{\beta}),\quad(i=1,\ldots,n),\
\beta,\widetilde{\beta}\in\mathcal{B}$ (16)

for uniformly sub-Gaussian $V_{i}$. That is,

$\displaystyle\max_{1\leq i\leq
n}K^{2}\mathbb{E}\left\\{e^{V_{i}^{2}/K^{2}}-1\right\\}\leq\sigma_{0}^{2},$
(17)

for positive $K$ and $\sigma_{0}$. Lemma 8.5 in van de Geer, (2000) presents
an exponential inequality for the supremum of the difference of these
processes.

######  Lemma 1 (Lemma 8.5 in van de Geer,, 2000).

Assume (16) and (17) and that
$\sup_{\beta\in\mathcal{B}}{\mathfrak{D}}(\beta,\beta_{0})\leq R$ for some
$\beta_{0}\in\Lambda$. Then for some constant $\widetilde{C}$ depending only
on $K$ and $\sigma_{0}$ and for all $\delta>0$ and $\sigma>0$ satisfying

$\displaystyle\delta\geq\widetilde{C}\max\left(\int_{\delta/(8\sigma)}^{R}H^{1/2}(u,\mathcal{B},{\mathfrak{D}})\,\mathrm{d}\,u,R\right),$
(18)

we have

$\displaystyle\mathbb{P}\left(\sup_{\beta\in\mathcal{B}}\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|\geq\delta,\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)\leq\widetilde{C}\exp\left[-\frac{\delta^{2}}{\widetilde{C}^{2}R^{2}}\right].$
(19)

This lemma is not directly applicable in our proofs, as in our setting we aim
to identify $\Lambda$ with the Sobolev space $\mathcal{H}^{m}(\mathcal{I})$
and $\mathcal{B}$ with a large subset of $\mathcal{H}^{m}(\mathcal{I})$ that
does not have finite entropy. That is, the entropy integral in (18) diverges.
But we will decompose $\mathcal{B}$ as

$\displaystyle\mathcal{B}=\bigcup_{M\geq 1}\mathcal{B}_{M},$

with $\mathcal{B}_{M}=\\{f\in\mathcal{B},1+J_{m}(f)\leq M\\}$, where in
general, $J_{m}$ can be any function
$J_{m}\colon\mathcal{B}\to\mathbbm{R}_{+}$, such that $\mathcal{B}_{M}$ has
finite entropy. In fact, under our assumptions it will be shown that

$\displaystyle H(\delta,\mathcal{B}_{M},{\mathfrak{D}})\leq
A_{0}\left(\frac{M}{\delta}\right)^{\frac{d}{m}},\quad\delta>0,$

for $A_{0}>0$ a constant, with the consequence that the corresponding entropy
integral behaves like

$\displaystyle\int_{0}^{\delta}H^{1/2}(u,\mathcal{B}_{M},{\mathfrak{D}})\,\mathrm{d}\,u\leq
A_{0}^{\star}M^{\frac{d}{2m}}\delta^{1-\frac{d}{2m}},\quad\delta>0,$ (20)

for some global $A_{0}^{\star}>0$ depending only on $m$. In (20), we have also
taken $\sigma\to\infty$ in the RHS of (18). This is permissible, because the
LHS of (20) converges for $2m>d$. Lemma 2 below greatly extends Lemma 1 by
establishing the modulus of continuity of

$\displaystyle\sup_{\beta\in\mathcal{B}}\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|,$

both in terms of ${\mathfrak{D}}$ and the function $J_{m}$. The result is not
only useful for our purposes but it can also be of general interest.

######  Lemma 2.

Assume (16), (17), (20) and that
$\sup_{\beta\in\mathcal{B}}{\mathfrak{D}}(\beta,\beta_{0})\leq R$ for some
$\beta_{0}\in\Lambda$. Then, there exists $c>0$ such that for all $T\geq c$ we
have

$\displaystyle\mathbb{P}\left(\sup_{\beta\in\mathcal{B}}\frac{\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|}{{\mathfrak{D}}^{1-\frac{d}{2m}}(\beta,\beta_{0})\\{1+J_{m}(\beta)\\}^{\frac{d}{2m}}}\geq
T,\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)\leq
c\exp\left[-\frac{T^{2}}{c^{2}}\right].$ (21)

###### Proof.

Under our assumptions, the conditions of Lemma 1 are satisfied for each
$\mathcal{B}_{M}$, hence the lemma is applicable. Our proof consists of
iterative application of (19) on each $\mathcal{B}_{M}$ and the peeling
technique (see, e.g., van de Geer,, 2000, p. 70). We break down the proof in
two steps where in the first step we deal with $\mathcal{B}_{M}$ and in the
second step pass to $\mathcal{B}$.

_Step 1._ We will prove that there exists a $c_{1}>0$ such that for all $M\geq
1$,

$\displaystyle\mathbb{P}\left(\sup_{\beta\in\mathcal{B}_{M}}\frac{\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|}{{\mathfrak{D}}^{1-\frac{d}{2m}}(\beta,\beta_{0})}\geq
T\left(\frac{M}{2}\right)^{\frac{d}{2m}},\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)\leq
c_{1}\exp\left[-\frac{T^{2}M^{\frac{d}{m}}}{c_{1}^{2}}\right].$ (22)

To show (22) let $T=2\widetilde{C}A_{0}^{\star}$ and use Boole’s inequality to
obtain

$\displaystyle\mathbb{P}$
$\displaystyle\left(\sup_{\beta\in\mathcal{B}_{M}}\frac{\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|}{{\mathfrak{D}}^{1-\frac{d}{2m}}(\beta,\beta_{0})}\geq
T\left(\frac{M}{2}\right)^{\frac{d}{2m}},\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)$
$\displaystyle\leq\sum_{s=1}^{\infty}\mathbb{P}\left(\sup_{\beta\in\mathcal{B}_{M},2^{-s}R\leq{\mathfrak{D}}(\beta,\beta_{0})\leq
2^{-s+1}R}\frac{\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|}{{\mathfrak{D}}^{1-\frac{d}{2m}}(\beta,\beta_{0})}\geq
T\left(\frac{M}{2}\right)^{\frac{d}{2m}},\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)$
$\displaystyle\leq\sum_{s=1}^{\infty}\mathbb{P}\left(\sup_{\beta\in\mathcal{B}_{M},{\mathfrak{D}}(\beta,\beta_{0})\leq
2^{-s+1}R}\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|\geq
T\left(2^{-s}R\right)^{1-\frac{d}{2m}}\left(\frac{M}{2}\right)^{\frac{d}{2m}},\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)$
$\displaystyle=\sum_{s=1}^{\infty}\mathbb{P}\left(\sup_{\beta\in\mathcal{B}_{M},{\mathfrak{D}}(\beta,\beta_{0})\leq
2^{-s+1}R}\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|\geq\widetilde{C}A_{0}^{\star}\left(2^{-s+1}R\right)^{1-\frac{d}{2m}}M^{\frac{d}{2m}},\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right).$

Recalling that $M\geq 1$, we may apply Lemma 1 on each of the summands to see
that

$\displaystyle\mathbb{P}$
$\displaystyle\left(\sup_{\beta\in\mathcal{B}_{M}}\frac{\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|}{{\mathfrak{D}}^{1-\frac{d}{2m}}(\beta,\beta_{0})}\geq
T\left(\frac{M}{2}\right)^{\frac{d}{2m}},\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)$
$\displaystyle\quad\leq\sum_{s=1}^{\infty}\widetilde{C}\exp\left[-\frac{|A_{0}^{\star}|^{2}(2^{-s+1}R)^{-\frac{d}{m}}M^{\frac{d}{m}}}{\widetilde{C}^{2}}\right]\leq
c_{1}\exp\left[-\frac{T^{2}M^{\frac{d}{m}}}{c_{1}^{2}}\right],$

for some $c_{1}>0$ depending only on $\widetilde{C}$, $m$ and $d$, as the
series converges. The proof of (22) is complete.

_Step 2._ In the second step we use again the peeling technique along with
(22) in order to complete the proof of the lemma. In particular, we have

$\displaystyle\mathbb{P}$
$\displaystyle\left(\sup_{\beta\in\mathcal{B}}\frac{\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|}{{\mathfrak{D}}^{1-\frac{d}{2m}}(\beta,\beta_{0})\\{1+J_{m}(\beta)\\}^{\frac{d}{2m}}}\geq
T,\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)$
$\displaystyle\leq\sum_{s=0}^{\infty}\mathbb{P}\left(\sup_{\beta\in\mathcal{B},2^{s}\leq
1+J_{m}(\beta)\leq
2^{s+1}}\frac{\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|}{{\mathfrak{D}}^{1-\frac{d}{2m}}(\beta,\beta_{0})\\{1+J_{m}(\beta)\\}^{\frac{d}{2m}}}\geq
T,\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)$
$\displaystyle\leq\sum_{s=0}^{\infty}\mathbb{P}\left(\sup_{\beta\in\mathcal{B},1+J_{m}(\beta)\leq
2^{s+1}}\frac{\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|}{{\mathfrak{D}}^{1-\frac{d}{2m}}(\beta,\beta_{0})}\geq
T2^{\frac{sd}{2m}},\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)$
$\displaystyle=\sum_{s=0}^{\infty}\mathbb{P}\left(\sup_{\beta\in\mathcal{B}_{2^{s+1}}}\frac{\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|}{{\mathfrak{D}}^{1-\frac{d}{2m}}(\beta,\beta_{0})}\geq
T\left(\frac{2^{s+1}}{2}\right)^{\frac{d}{2m}},\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\leq\sigma^{2}\right)$
$\displaystyle\leq\sum_{s=0}^{\infty}c_{1}\exp\left[-\frac{T^{2}2^{\frac{(s+1)d}{m}}}{c_{1}^{2}}\right]\leq
c\exp\left[-\frac{T^{2}}{c^{2}}\right],$

for some $c>0$, as the series converges. This completes the proof of the
lemma.

∎

###  A.5 Proof of Theorem 1

Before presenting the proof of Theorem 1, it will be helpful to state and
prove an auxiliary lemma involving ratios of sequences.

######  Lemma 3.

Let $\\{x_{n}\\}_{n=1}^{\infty}$ denote a sequence of non-negative real
numbers and for $a>b$ denote the sequences $y_{n}=x_{n}^{a}/(1+x_{n})^{b}$ and
$z_{n}=x_{n}/(1+x_{n})$. Then,

  1. A.

For any $M>0$ such that $\sup_{n}y_{n}\leq M$ there exists a finite $B=B(M)>0$
such that $\sup_{n}x_{n}\leq B$.

  2. B.

If there exists an $\varepsilon\in(0,1)$ such that $z_{n}\leq 1-\varepsilon$
for all large $n$, then we can find finite $B=B(\varepsilon)>0$ such that
$\sup_{n}x_{n}\leq B$.

###### Proof of Lemma 3.

We start by proving part A. Assume for contradiction that
$\\{x_{n}\\}_{n=1}^{\infty}$ is unbounded. Then $\\{x_{n}\\}_{n=1}^{\infty}$
contains a subsequence $\\{x_{n_{k}}\\}_{k=1}^{\infty}$ that diverges to
$\infty$. But then $y_{n_{k}}\sim x_{n_{k}}^{a-b}$ where $a_{n}\sim b_{n}$
means that $\lim a_{n}/b_{n}=1$. But then, since, by assumption, $a>b$,
$y_{n_{k}}$ also diverges contradicting its boundedness. It follows that
$\\{x_{n}\\}_{n=1}^{\infty}$ must be bounded, as was to be shown.

To prove part B, it suffices to note that if $\\{x_{n}\\}_{n=1}^{\infty}$
contained a divergent subsequence, $\\{x_{n_{k}}\\}_{k=1}^{\infty}$, then
$\lim_{k\to\infty}z_{n_{k}}=1$ contradicting the fact that $z_{n}\leq
1-\varepsilon$ eventually. ∎

The proof of Theorem 1 consists of the convexity step of van de Geer, (2002)
and the derivation of the modulus of continuity of the empirical process with
respect to the semi-metric $\left\|\cdot\right\|_{n,p}$ that is induced by our
penalized M-estimator. Let $L_{n}(\beta)$ denote the objective function, that
is,

$\displaystyle
L_{n}(\beta)=\frac{1}{n}\sum_{i=1}^{n}\rho\left(Y_{i}-\sum_{j=1}^{p}X_{i}\left(\mathbf{t}_{j}\right)\beta\left(\mathbf{t}_{j}\right)\mu\left(A_{j}\right)\right)+\lambda
J_{m}^{2}\left(\beta\right).$

By definition of the minimizer, we have $L_{n}(\widehat{\beta}_{n})\leq
L_{n}(\beta)$ for all $\beta\in\mathcal{H}^{m}(\mathbbm{R}^{d})$. Moreover,
$L_{n}$ is convex, as it is the sum of two convex functions. Define the convex
combination

$\widetilde{\beta}_{n}=\gamma_{n}\widehat{\beta}_{n}+(1-\gamma_{n})\beta_{0}\qquad\mbox{with}\qquad\gamma_{n}=1/(1+\|\widehat{\beta}_{n}-\beta_{0}\|_{n,p}).$
(23)

Clearly $\gamma_{n}\in(0,1)$ and, by convexity,

$\displaystyle
L_{n}(\widetilde{\beta}_{n})\leq\gamma_{n}L_{n}(\widehat{\beta}_{n})+(1-\gamma_{n})L_{n}(\beta_{0})\leq
L_{n}(\beta_{0}),$ (24)

as, by assumption, $\beta_{0}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$. Define for
simplicity

$\displaystyle\;M_{n}(\beta):=\frac{1}{n}\sum_{i=1}^{n}\rho\left(\epsilon_{i}+d_{i}+\sum_{j=1}^{p}X_{i}\left(\mathbf{t}_{j}\right)\left(\beta_{0}\left(\mathbf{t}_{j}\right)-\beta\left(\mathbf{t}_{j}\right)\right)\mu\left(A_{j}\right)\right),$

so that $L_{n}(\beta)=M_{n}(\beta)+\lambda J^{2}_{m}(\beta)$. Furthermore,
define

$\displaystyle
M(\beta):=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{\epsilon_{i}}\left\\{\rho\left(\epsilon_{i}+d_{i}+\sum_{j=1}^{p}X_{i}\left(t_{j}\right)\left(\beta_{0}\left(\mathbf{t}_{j}\right)-\beta\left(\mathbf{t}_{j}\right)\right)\mu\left(A_{j}\right)\right)\right\\},$

where $\mathbb{E}_{\epsilon_{i}}\\{\cdot\\}$ denotes expectation with respect
to the $\epsilon_{i}$ (recall that the $\epsilon_{i}$ and $X_{i}$ are assumed
to be independent and the $\epsilon_{i}$ are independent). With this notation,
rearranging (24) yields

$\displaystyle\left(M(\widetilde{\beta}_{n})-M(\beta_{0})\right)+\lambda
J_{m}^{2}(\widetilde{\beta}_{n})\leq\left(M_{n}(\beta_{0})-M(\beta_{0})\right)-\left(M_{n}(\widetilde{\beta}_{n})-M(\widetilde{\beta}_{n})\right)+\lambda
J_{m}^{2}(\beta_{0}).$ (25)

Our approach consists of deriving a lower bound on the LHS of (25) in terms of
$\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}$ and an upper bound on the RHS,
also in terms of $\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}$ on a set
$A_{\varepsilon}\subset\Omega$ with $\mathbb{P}(A_{\varepsilon})\geq
1-\varepsilon$ for any given $\varepsilon>0$. Combining these two bounds
appropriately will yield

$\displaystyle\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p}=O_{\mathbb{P}}\left(\log^{2}(n)\left\\{n^{-\frac{m}{2m+d}}+\max_{1\leq
j\leq
p}\operatorname{diam}^{\kappa}(A_{j})\right\\}\right)\quad\text{and}\quad
J_{m}\left(\widetilde{\beta}_{n}\right)=O_{\mathbb{P}}(1),$ (26)

under $(i)$ and $(ii)$ of the theorem. The definition of the convex
combination $\widetilde{\beta}_{n}$ will then allow us to derive the same rate
of convergence for $\|\widehat{\beta}_{n}-\beta_{0}\|_{n,p}$ as well as the
boundedness of $J_{m}(\widehat{\beta}_{n})$.

_Step 1._ We begin by deriving a lower bound on the LHS of (25). By (A4) and
Taylor’s theorem, for every $\beta\in\mathcal{H}^{m}(\mathbbm{R}^{d})$, we
find

$\displaystyle\mathbb{E}_{\epsilon_{i}}\left\\{\rho\left(\epsilon_{i}+d_{i}+\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)\right\\}$
$\displaystyle=g\left(d_{i}+\langle X_{i},\beta_{0}-\beta\rangle_{p}\right)$
$\displaystyle=g\left(\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)+d_{i}g^{\prime}\left(\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)+O(d_{i}^{2})$
$\displaystyle=g\left(\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)+d_{i}g^{\prime}\left(0\right)+O\left(\left|d_{i}\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right|\right)+O(d_{i}^{2})$
$\displaystyle=g\left(\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)+O\left(\left|d_{i}\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right|\right)+O(d_{i}^{2})$

where we have used the fact that, by (A4), $g^{\prime}(0)=0$ and
$\sup_{t\in\mathbbm{R}}|g^{\prime\prime}(t)|\leq B$ for some finite $B>0$.
Expanding $\mathbb{E}_{\epsilon_{i}}\\{\rho(\epsilon_{i}+d_{i})\\}$ in a
similar manner, we see that

$\displaystyle\mathbb{E}_{\epsilon_{i}}\\{\rho(\epsilon_{i}+d_{i})\\}=g(d_{i})=g(0)+d_{i}g^{\prime}(0)+O(d_{i}^{2})=g(0)+O(d_{i}^{2}).$

Combining these two expansions and averaging, we may deduce the existence of a
large enough $C>0$ such that

$\displaystyle M(\beta)-M(\beta_{0})$ (27)
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{\epsilon_{i}}\left\\{\rho\left(\epsilon_{i}+d_{i}+\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)-\rho\left(\epsilon_{i}+d_{i}\right)\right\\}$
$\displaystyle\geq\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{\epsilon_{i}}\left\\{\rho\left(\epsilon_{i}+\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)-\rho\left(\epsilon_{i}\right)\right\\}-\frac{C}{n}\sum_{i=1}^{n}d_{i}^{2}-\frac{C}{n}\sum_{i=1}^{n}\left|d_{i}\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right|$
$\displaystyle\geq\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{\epsilon_{i}}\left\\{\rho\left(\epsilon_{i}+\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)-\rho\left(\epsilon_{i}\right)\right\\}-\frac{C}{n}\sum_{i=1}^{n}d_{i}^{2}-C\left(\frac{1}{n}\sum_{i=1}^{n}d_{i}^{2}\right)^{\frac{1}{2}}\left\|\beta-\beta_{0}\right\|_{n,p},$

where to obtain the second inequality we have used the Cauchy-Schwarz
inequality. Now, again by the Cauchy-Schwarz inequality,

$\displaystyle\max_{1\leq i\leq n}\left|\langle
X_{i},\beta-\beta_{0}\rangle_{p}\right|\leq\max_{1\leq i\leq
n}\left\|X_{i}\right\|_{\infty}\left\|\beta-\beta_{0}\right\|_{p},$

for every $\beta\in\mathcal{H}^{m}(\mathcal{I})$. Now, by Boole’s and Markov’s
inequalities, for all $u>0$, we have

$\displaystyle\mathbb{P}\left(\max_{1\leq i\leq
n}\left\|X_{i}\right\|_{\infty}>u\right)\leq
n\mathbb{P}\left(\left\|X_{1}\right\|_{\infty}>u\right)\leq
n\mathbb{E}\left\\{\exp\left[t\left\|X\right\|_{\infty}^{2}\right]\right\\}e^{-tu^{2}}\leq
Dne^{-tu^{2}},$ (28)

for some $D>0$. Setting $u=\sqrt{C\log(n)}$ for sufficiently large $C>0$ we
now see that the RHS of (28) is summable whence, by the Borel-Cantelli lemma,
we may conclude that

$\displaystyle\max_{1\leq i\leq
n}\left\|X_{i}\right\|_{\infty}\leq\sqrt{C\log(n)},$

for all large $n$, almost surely. But, by definition of $J_{m}(\cdot)$, we now
see that

$\displaystyle\max_{1\leq i\leq n}\left|\langle
X_{i},\beta-\beta_{0}\rangle_{p}\right|\leq\sqrt{C\log(n)}J_{m}(\beta-\beta_{0}),$
(29)

almost surely. By our assumptions,
$\beta_{0}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$, hence $J_{m}(\beta_{0})$ is
bounded. Without loss of generality, we thus take $J_{m}(\beta_{0})\leq 1$.
From (29), it follows that we can find a large enough constant
$D_{\varepsilon}\geq 1$ such that

$\displaystyle\max_{1\leq i\leq n}\frac{\left|\langle
X_{i},\beta-\beta_{0}\rangle_{p}\right|}{D_{\varepsilon}\sqrt{\log(n)}\left(1+J_{m}(\beta)\right)}\leq\varepsilon$

with high probability, where $\varepsilon$ is the constant in (A6). It follows
that for all large $n$, $n\geq n_{0}$, say,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{\epsilon_{i}}\left\\{\rho\left(\epsilon_{i}+\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)-\rho\left(\epsilon_{i}\right)\right\\}$
$\displaystyle\geq\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{\epsilon_{i}}\left\\{\rho\left(\epsilon_{i}+\frac{\langle
X_{i},\beta_{0}-\beta\rangle_{p}}{D_{\varepsilon}\sqrt{\log(n)}\\{1+J_{m}(\beta)\\}}\right)-\rho\left(\epsilon_{i}\right)\right\\}$
(30)
$\displaystyle\geq\varepsilon\frac{\left\|\beta-\beta_{0}\right\|^{2}_{n,p}}{D_{\varepsilon}^{2}\log(n)\\{1+J_{m}(\beta)\\}^{2}}.$

The first inequality in (30) follows from the fact that the function
$g(t)=\mathbb{E}_{\epsilon_{i}}\left\\{\rho\left(\epsilon_{i}+t\right)\right\\}$
is non-decreasing in $|t|$ due to (A5) and the convexity of $\rho$ from (A4),
and because of our assumption that $D_{\varepsilon}\geq 1$, meaning that also
$D_{\varepsilon}\sqrt{\log(n)}(1+J_{m}(\beta))\geq 1$ for all large $n$ and
$\beta\in\mathcal{H}^{m}(\mathbbm{R}^{d})$. In the second inequality in (30),
we have used (A6), and the definition of the semi-norm $\|\cdot\|_{n,p}$ from
(10).

The inequality (30) holds for every fixed
$\beta\in\mathcal{H}^{m}\left(\mathbbm{R}^{d}\right)$, hence what we have
shown in (27) and (30) is

$\displaystyle\inf_{\beta\in\mathcal{H}^{m}\left(\mathbbm{R}^{d}\right)}\left[\frac{M(\beta)-M(\beta_{0})+\frac{C}{n}\sum_{i=1}^{n}d_{i}^{2}+C\left(\frac{1}{n}\sum_{i=1}^{n}d_{i}^{2}\right)^{\frac{1}{2}}\left\|\beta-\beta_{0}\right\|_{n,p}}{\varepsilon\frac{\left\|\beta-\beta_{0}\right\|^{2}_{n,p}}{D_{\varepsilon}^{2}\log(n)\left(1+J_{m}(\beta)\right)^{2}}}\right]\geq
1.$

To establish a similar inequality for the random
$\widetilde{\beta}_{n}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$ it suffices to note
that, by the last inequality, we have

$\displaystyle M(\widetilde{\beta}_{n})-M(\beta_{0})$
$\displaystyle=M(\widetilde{\beta}_{n})-M(\beta_{0})+\frac{C}{n}\sum_{i=1}^{n}d_{i}^{2}+C\left(\frac{1}{n}\sum_{i=1}^{n}d_{i}^{2}\right)^{1/2}\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p}$
$\displaystyle\quad-\frac{C}{n}\sum_{i=1}^{n}d_{i}^{2}-C\left(\frac{1}{n}\sum_{i=1}^{n}d_{i}^{2}\right)^{1/2}\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p}$
$\displaystyle\geq\varepsilon\frac{\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|^{2}_{n,p}}{D_{\varepsilon}^{2}\log(n)\left(1+J_{m}(\widetilde{\beta}_{n})\right)^{2}}-\frac{C}{n}\sum_{i=1}^{n}d_{i}^{2}-C\left(\frac{1}{n}\sum_{i=1}^{n}d_{i}^{2}\right)^{1/2}\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p},$
(31)

which provides the desired lower bound for (25).

_Step 2._ The next step in our proof is the derivation of an upper bound for
the RHS of (25) in terms of $\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}$. To
that end, for each $\beta\in\mathcal{H}^{m}(\mathbbm{R}^{d})$ define the
independent centered processes

$\displaystyle U_{i,\beta}=\rho\left(\epsilon_{i}+d_{i}+\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)-\mathbb{E}_{\epsilon_{i}}\left\\{\rho\left(\epsilon_{i}+d_{i}+\langle
X_{i},\beta_{0}-\beta\rangle_{p}\right)\right\\},\quad(i=1,\ldots,n).$

With this definition we find

$\displaystyle\left(M_{n}(\beta_{0})-M(\beta_{0})\right)-\left(M_{n}(\widetilde{\beta}_{n})-M(\widetilde{\beta}_{n})\right)=\frac{1}{n}\sum_{i=1}^{n}\left(U_{i,\beta_{0}}-U_{i,\widetilde{\beta}_{n}}\right).$

Notice that, by definition (23) of $\widetilde{\beta}_{n}$, we have

$\displaystyle\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p}=\gamma_{n}\left\|\widehat{\beta}_{n}-\beta_{0}\right\|_{n,p}=\frac{\left\|\widehat{\beta}_{n}-\beta_{0}\right\|_{n,p}}{1+\left\|\widehat{\beta}_{n}-\beta_{0}\right\|_{n,p}}\leq
1.$

Moreover, by (A1), $\\{\mathbf{t}_{j}\\}_{j=1}^{p}\subset\mathcal{I}$. These
two facts imply that we may restrict attention to the subset
$\mathcal{B}\subset\mathcal{H}^{m}(\mathcal{I})$ given by
$\mathcal{B}=\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):\|\beta-\beta_{0}\|_{n,p}\leq
1\\}$. We will apply Lemma 2 to these processes after identifying
${\mathfrak{D}}(\beta,\widetilde{\beta})$ in (15) with
$\|\beta-\widetilde{\beta}\|_{n,p}$, and
${\mathfrak{D}}_{i}(\beta,\widetilde{\beta})$ with $|\langle
X_{i},\beta-\widetilde{\beta}\rangle_{p}|$. By the Lipschitz continuity of
$\rho$ assumed in (A4), we have

$\displaystyle\left|U_{i,\beta}-U_{i,\widetilde{\beta}}\right|\leq
2C_{2}\left|\langle
X_{i},\beta-\widetilde{\beta}\rangle_{p}\right|,\quad(i=1,\ldots,n),$

and the constant random variables $\\{2\,C_{2}\\}_{i=1}^{n}$ are trivially
uniformly sub-Gaussian. In order to apply Lemma 2, we need to verify (20) by
showing that the entropy integral for the class of functions
$\mathcal{B}_{M}:=\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):\|\beta-\beta_{0}\|_{n,p}\leq
1,1+J_{m}(\beta)\leq M\\}$ for $M\geq 1$ behaves like

$\displaystyle\int_{0}^{\delta}H^{1/2}(u,\mathcal{B}_{M},d)\,\mathrm{d}\,u\lesssim
M^{d/(2m)}\delta^{1-d/(2m)},\quad\delta>0.$ (32)

To verify the bound in (32), notice that, by our assumptions, for every
$\beta\in\mathcal{H}^{m}(\mathcal{I})$, $\|\beta\|_{n,p}$ is up to a constant
dominated by $\|\beta\|_{p}$, with high probability. This is because, by the
Cauchy-Schwarz inequality,

$\displaystyle\left\|\beta\right\|_{n,p}^{2}=\frac{1}{n}\sum_{i=1}^{n}\left|\sum_{j=1}^{p}X_{i}(\mathbf{t}_{j})\beta(\mathbf{t}_{j})\mu(A_{j})\right|^{2}$
$\displaystyle\leq\frac{1}{n}\sum_{i=1}^{n}\left[\left\\{\sum_{j=1}^{p}\left|X_{i}(\mathbf{t}_{j})\right|^{2}\mu(A_{j})\right\\}\left\\{\sum_{j=1}^{p}\left|\beta(\mathbf{t}_{j})\right|^{2}\mu(A_{j})\right\\}\right].$

But, by (A3),

$\displaystyle\mathbb{E}\left\\{\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{p}\left|X_{i}(\mathbf{t}_{j})\right|^{2}\mu(A_{j})\right\\}\leq\sup_{\mathbf{t}\in\mathcal{I}}\mathbb{E}\left\\{|X_{1}(\mathbf{t})|^{2}\right\\}\mu(\mathcal{I}),$

where we have used the fact that $A_{1},\ldots,A_{p}$ form a partition of
$\mathcal{I}$. The finiteness of the RHS implies the existence of a $M_{2}$
such that $\|\beta\|_{n,p}\leq M_{2}\|\beta\|_{p}$ with high probability. It
follows that a cover of $\mathcal{B}_{M}$ in the $\|\cdot\|_{p}$-semimetric
also provides a cover in the $\|\cdot\|_{n,p}$-semimetric. In other words,

$\displaystyle H(u,\mathcal{B}_{M},\left\|\cdot\right\|_{n,p})\leq
H(u/M_{2},\mathcal{B}_{M},\left\|\cdot\right\|_{p}),\quad u>0,$

with high probability. Now, $1+J_{m}(\beta)\leq M$ implies both that
$1+\|\beta\|_{p}\leq M$ and that $1+I_{m}(\beta)\leq M$. Hence, we obtain the
following inclusions:

$\displaystyle\mathcal{B}_{M}$
$\displaystyle\subset\left\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):\left\|\beta-\beta_{0}\right\|_{n,p}\leq
1,1+\left\|\beta\right\|_{p}\leq M,1+I_{m}(\beta)\leq M\right\\}$
$\displaystyle\subset\left\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):1+\left\|\beta\right\|_{p}\leq
M,1+I_{m}(\beta)\leq M\right\\}.$

Therefore,

$\displaystyle H(u,\mathcal{B}_{M},\left\|\cdot\right\|_{p})\leq
H\left(u,\left\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):1+\left\|\beta\right\|_{p}\leq
M,1+I_{m}(\beta)\leq M\right\\},\left\|\cdot\right\|_{p}\right),\quad u>0.$

We now bound the latter entropy. Since, for any
$\beta\in\mathcal{H}^{m}(\mathcal{I})$, we have

$\displaystyle\sum_{j=1}^{p}|\beta(\mathbf{t}_{j})|^{2}\mu(A_{j})\leq\sup_{\mathbf{t}\in\mathcal{I}}|\beta(\mathbf{t})|^{2}\mu(\mathcal{I}),$

it follows that

$\displaystyle H$
$\displaystyle\left(u,\left\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):1+\left\|\beta\right\|_{p}\leq
M,1+I_{m}(\beta)\leq M\right\\},\left\|\cdot\right\|_{p}\right)$
$\displaystyle\quad\leq
H_{\infty}\left(u/\sqrt{\mu(\mathcal{I})},\left\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):1+\left\|\beta\right\|_{p}\leq
M,1+I_{m}(\beta)\leq M\right\\}\right),\quad u>0,$

with $H_{\infty}(u,\mathcal{F})$ denoting the entropy in the supremum norm of
a class of functions $\mathcal{F}$. We will now show that there exists a
constant $c_{0}>0$, independent of $M$, such that

$\displaystyle\left\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):1+\left\|\beta\right\|_{p}\leq
M,1+I_{m}(\beta)\leq
M\right\\}\subset\left\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):\left\|\beta\right\|_{\mathcal{H}^{m}(\mathcal{I})}\leq
c_{0}M\right\\}.$

Indeed, take a $\beta\in\mathcal{H}^{m}(\mathcal{I})$ such that
$1+\|\beta\|_{p}\leq M$ and $1+I_{m}(\beta)\leq M$. By Theorem 3.4 of Utreras,
(1988), there exists a global $c_{0}>0$ such that

$\displaystyle\int_{\mathcal{I}}|\beta(\mathbf{t})|^{2}\,\mathrm{d}\,\mathbf{t}$
$\displaystyle\leq\frac{c_{0}}{p}\sum_{j=1}^{p}|\beta(\mathbf{t}_{j})|^{2}+c_{0}I_{m}^{2}(\beta)\leq
c_{0}(M-1)^{2}+c_{0}(M-1)^{2}=2c_{0}(M-1)^{2}.$

Hence, for every $\beta\in\mathcal{H}^{m}(\mathcal{I})$, by definition of
$\|\cdot\|_{\mathcal{H}^{m}(\mathcal{I})}$,

$\displaystyle\left\|\beta\right\|^{2}_{\mathcal{H}^{m}(\mathcal{I})}\leq\int_{\mathcal{I}}|\beta(\mathbf{t})|^{2}\,\mathrm{d}\,\mathbf{t}+I_{m}^{2}(\beta)\leq
2c_{0}(M-1)^{2}+(M-1)^{2}\leq\max\\{2c_{0},1\\}M^{2},$

as $M\geq 1$. Combining all the above, we have shown that

$\displaystyle H\left(u,\mathcal{B}_{M},\left\|\cdot\right\|_{n,p}\right)\leq
H_{\infty}\left(u/(\sqrt{\mu(\mathcal{I})}M_{2}),\left\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):\left\|\beta\right\|_{\mathcal{H}^{m}(\mathcal{I})}\leq
c_{0}M\right\\}\right),\quad u>0,$ (33)

with high probability. But the RHS of (33) is the supremum entropy of the
closed ball in $\mathcal{H}^{m}(\mathcal{I})$ with radius $c_{0}M$. Since, by
our assumptions $2m>d$, Proposition 6 in Cucker and Smale, (2002) now yields

$\displaystyle
H_{\infty}\left(u,\left\\{\beta\in\mathcal{H}^{m}(\mathcal{I}):\left\|\beta\right\|_{\mathcal{H}^{m}(\mathcal{I})}\leq
c_{0}M\right\\}\right)\lesssim M^{\frac{d}{m}}u^{-\frac{d}{m}},\quad u>0,$

which, upon taking square roots and integrating, yields (32).

Since (32) holds, applying Lemma 2 next leads to

$\displaystyle\sup_{\beta\in\mathcal{H}^{m}\left(\mathbbm{R}^{d}\right):\left\|\beta-\beta_{0}\right\|_{n,p}\leq
1}\frac{\left|\frac{1}{n}\sum_{i=1}^{n}\left(U_{i,\beta}-U_{i,\beta_{0}}\right)\right|}{\left\|\beta-\beta_{0}\right\|_{n,p}^{1-\frac{d}{2m}}\\{1+J_{m}(\beta)\\}^{\frac{d}{m}}}=O_{\mathbb{P}}(n^{-\frac{1}{2}}).$

This inequality provides the desired upper bound for the RHS in (25), namely,

$\displaystyle\left(M_{n}(\widetilde{\beta}_{n})-M(\widetilde{\beta}_{n})\right)-\left(M_{n}(\beta_{0})-M(\beta_{0})\right)=O_{\mathbb{P}}(n^{-\frac{1}{2}})\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}^{1-\frac{d}{2m}}\\{1+J_{m}(\widetilde{\beta}_{n})\\}^{\frac{d}{2m}}.$
(34)

_Step 3._ The third step of the proof involves combining the lower and upper
bounds derived in the previous steps in order to obtain a rate of convergence
for $\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}$. Indeed, plugging (A.5) and
(34) into (25) yields

$\displaystyle\varepsilon\frac{\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|^{2}_{n,p}}{D_{\varepsilon}^{2}\log(n)\left(1+J_{m}(\widetilde{\beta}_{n})\right)^{2}}+\lambda
J_{m}^{2}\left(\widetilde{\beta}_{n}\right)$ $\displaystyle\leq
O_{\mathbb{P}}(n^{-\frac{1}{2}})\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}^{1-\frac{d}{2m}}\\{1+J_{m}(\widetilde{\beta}_{n})\\}^{\frac{d}{2m}}+\lambda
J_{m}^{2}\left(\beta_{0}\right)$
$\displaystyle\quad+\frac{C}{n}\sum_{i=1}^{n}d_{i}^{2}+C\left(\frac{1}{n}\sum_{i=1}^{n}d_{i}^{2}\right)^{\frac{1}{2}}\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p},$
(35)

with high probability.

We determine next the order of $n^{-1}\sum_{i=1}^{n}d_{i}^{2}$, for
$\\{d_{i}\\}_{i=1}^{n}$ defined in (4). It is easy to see that

$\displaystyle\left|d_{i}\right|$
$\displaystyle\leq\left|\sum_{j=1}^{p}\int_{A_{j}}\left(\beta_{0}(\mathbf{t})-\beta_{0}(\mathbf{t}_{j})\right)X_{i}(\mathbf{t})\,\mathrm{d}\,\mathbf{t}\right|+\left|\sum_{j=1}^{p}\int_{A_{j}}\beta_{0}(\mathbf{t}_{j})(X_{i}(\mathbf{t})-X_{i}(\mathbf{t}_{j}))\,\mathrm{d}\,\mathbf{t}\right|$
$\displaystyle=I+II,$

say. To bound $I$, let us first observe that for $2m>d+1$, by the Sobolev
embedding theorem (Adams and Fournier,, 2003, Theorem 4.12),
$\mathcal{H}^{m}(\mathcal{I})$ embeds on the space of Lipschitz continuous
functions on $\mathcal{I}$ denoted by $\mathcal{C}^{0,1}(\mathcal{I})$. In
other words, for every $\beta\in\mathcal{H}^{m}(\mathcal{I})$ we have
$\beta\in\mathcal{C}^{0,1}(\mathcal{I})$ and there exists a global $C>0$ such
that

$\displaystyle\sup_{\mathbf{t},\mathbf{s}\in\mathcal{I}}\frac{\left|\beta\left(\mathbf{t}\right)-\beta\left(\mathbf{s}\right)\right|}{\left\|\mathbf{t}-\mathbf{s}\right\|_{\mathbbm{R}^{d}}}$
$\displaystyle\leq
C\left\\{\int_{\mathcal{I}}|\beta\left(\mathbf{t}\right)|^{2}\,\mathrm{d}\,\mathbf{t}+\sum_{m_{1}+\ldots+m_{d}=m}\binom{m}{m_{1},\ldots,m_{d}}\int_{\mathcal{I}}\left(\frac{\partial^{m}\beta(\mathbf{t})}{\partial
t_{1}^{m_{1}}\ldots\partial
t_{d}^{m_{d}}}\right)^{2}\,\mathrm{d}\,\mathbf{t}\right\\}^{\frac{1}{2}}$
$\displaystyle=C\left\|\beta\right\|_{\mathcal{H}^{m}(\mathcal{I})}.$

Recalling that, by assumption, $\beta_{0}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$,
hence also $\beta_{0}\in\mathcal{H}^{m}(\mathcal{I})$, hence

$\displaystyle
I\leq\sum_{j=1}^{p}\int_{A_{j}}\left|\left(\beta_{0}(\mathbf{t})-\beta_{0}(\mathbf{t}_{j})\right)X_{i}(\mathbf{t})\right|\,\mathrm{d}\,\mathbf{t}$
$\displaystyle\leq
C\left\|\beta_{0}\right\|_{\mathcal{H}^{m}(\mathcal{I})}\sup_{\mathbf{t}\in\mathcal{I}}\left|X_{i}(\mathbf{t})\right|\sum_{j=1}^{p}\int_{A_{j}}\left\|\mathbf{t}-\mathbf{t}_{j}\right\|_{\mathbbm{R}^{d}}\,\mathrm{d}\,\mathbf{t}$
$\displaystyle\leq
C\left\|\beta_{0}\right\|_{\mathcal{H}^{m}(\mathcal{I})}\sup_{\mathbf{t}\in\mathcal{I}}\left|X_{i}(\mathbf{t})\right|\mu(\mathcal{I})\max_{1\leq
j\leq p}\operatorname{diam}(A_{j}),$

for some $C>0$ where the last inequality follows from the fact that the sets
$A_{j}$ form a partition of $\mathcal{I}$. Now note that, by (A3),
$\mathbb{E}\\{\|X_{i}\|_{\infty}\\}<\infty$, as exponential moments imply the
existence of polynomial moments so that

$\displaystyle I=O_{\mathbb{P}}\left(\max_{1\leq j\leq
p}\operatorname{diam}(A_{j})\right).$ (36)

To bound $II$, observe that

$\displaystyle
II\leq\sup_{\mathbf{t}\in\mathcal{I}}\left|\beta_{0}(\mathbf{t})\right|\sum_{j=1}^{p}\int_{A_{j}}\left|X_{i}\left(\mathbf{t}\right)-X_{i}\left(\mathbf{t}_{j}\right)\right|\,\mathrm{d}\,\mathbf{t}.$

Note that $\sup_{\mathbf{t}\in\mathcal{I}}\left|\beta_{0}(\mathbf{t})\right|$
is finite, as, by (A1), $\mathcal{I}$ is bounded and $\beta_{0}$ is continuous
for all $m\geq 1$. By the Cauchy-Schwarz inequality and (A3), we have

$\displaystyle\mathbb{E}\left\\{\left|X_{1}\left(\mathbf{t}\right)-X_{1}\left(\mathbf{t}_{j}\right)\right|\right\\}\leq\left[\mathbb{E}\left\\{\left|X_{1}\left(\mathbf{t}\right)-X_{1}\left(\mathbf{t}_{j}\right)\right|^{2}\right\\}\right]^{1/2}\leq\sqrt{C_{1}}\left\|\mathbf{t}-\mathbf{t}_{j}\right\|^{\kappa}_{\mathbbm{R}^{d}},$

which, by Markov’s inequality, leads to

$\displaystyle II=O_{\mathbb{P}}\left(\max_{1\leq j\leq
p}\operatorname{diam}^{\kappa}(A_{j})\right).$ (37)

As $\kappa\leq 1$, the bounds in (36) and (37) imply that

$\frac{1}{n}\sum_{i=1}^{n}d_{i}^{2}=O_{\mathbb{P}}(\max_{1\leq j\leq
p}\operatorname{diam}^{2\kappa}(A_{j}))$ (38)

so that

$\displaystyle\underbrace{\varepsilon\frac{\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|^{2}_{n,p}}{\log(n)D_{\varepsilon}^{2}\left\\{1+J_{m}\left(\widetilde{\beta}_{n}\right)\right\\}^{2}}+\lambda
J_{m}^{2}\left(\widetilde{\beta}_{n}\right)}_{a}$
$\displaystyle\leq\underbrace{O_{\mathbb{P}}(n^{-\frac{1}{2}})\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}^{1-\frac{d}{2m}}\\{1+J_{m}(\widetilde{\beta}_{n})\\}^{\frac{d}{2m}}}_{b}+\underbrace{\lambda
J_{m}^{2}\left(\beta_{0}\right)}_{c}$
$\displaystyle\quad+\underbrace{O_{\mathbb{P}}\left(\max_{1\leq j\leq
p}\operatorname{diam}^{2\kappa}(A_{j})\right)}_{d}$
$\displaystyle\quad+\underbrace{O_{\mathbb{P}}\left(\max_{1\leq j\leq
p}\operatorname{diam}^{\kappa}(A_{j})\right)\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p}}_{e}.$
(39)

We next show that the bound in (A.5) implies (26). The former is an inequality
of the form $a\leq b+c+d+e$ for real numbers $a,b,c,d,e$ and therefore we must
have either $a\leq 4b$ or $a\leq 4c$ or $a\leq 4d$ or $a\leq 4e$ (for if that
were not true we would have $a>b+c+d+e$). We consider each one of these four
possible cases and show that in all cases (26) holds.

_Case 1._ Starting from the first possible case, we have $a\leq 4\,b$, that
is,

$\displaystyle\varepsilon\frac{\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|^{2}_{n,p}}{\log(n)D_{\varepsilon}^{2}\left\\{1+J_{m}\left(\widetilde{\beta}_{n}\right)\right\\}^{2}}+\lambda
J_{m}^{2}\left(\widetilde{\beta}_{n}\right)$ $\displaystyle\leq
O_{\mathbb{P}}(n^{-\frac{1}{2}})\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p}^{1-\frac{d}{2m}}\left\\{1+J_{m}(\widetilde{\beta}_{n})\right\\}^{\frac{d}{2m}}.$

Since both terms on the LHS are non-negative, this inequality implies both

$\displaystyle\varepsilon\frac{\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|^{2}_{n,p}}{\log(n)D_{\varepsilon}^{2}\left\\{1+J_{m}\left(\widetilde{\beta}_{n}\right)\right\\}^{2}}\leq
O_{\mathbb{P}}(n^{-\frac{1}{2}})\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}^{1-\frac{d}{2m}}\left\\{1+J_{m}(\widetilde{\beta}_{n})\right\\}^{\frac{d}{2m}},$
(40)

and

$\displaystyle\lambda J_{m}^{2}\left(\widetilde{\beta}_{n}\right)\leq
O_{\mathbb{P}}(n^{-\frac{1}{2}})\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}^{1-\frac{d}{2m}}\left\\{1+J_{m}(\widetilde{\beta}_{n})\right\\}^{\frac{d}{2m}}.$
(41)

We will use such a separation into two inequalities also in the analysis of
the remaining cases that follow. From (40) we get

$\displaystyle\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p}\leq
O_{\mathbb{P}}\left(\log^{\frac{2m}{2m+d}}(n)n^{-\frac{m}{2m+d}}\right)\left\\{1+J_{m}(\widetilde{\beta}_{n})\right\\}^{\frac{4m+d}{2m+d}}$
(42)

which upon plugging into (41) yields

$\displaystyle\lambda J_{m}^{2}\left(\widetilde{\beta}_{n}\right)\leq
O_{\mathbb{P}}\left(n^{-\frac{2m}{2m+d}}\log^{\frac{2m-d}{2m+d}}(n)\right)\left\\{1+J_{m}(\widetilde{\beta}_{n})\right\\}^{\frac{4m}{2m+d}},$

or equivalently,

$\displaystyle\frac{J_{m}^{2}\left(\widetilde{\beta}_{n}\right)}{\left\\{1+J_{m}(\widetilde{\beta}_{n})\right\\}^{\frac{4m}{2m+d}}}\leq
O_{\mathbb{P}}\left(n^{-\frac{2m}{2m+d}}\log^{\frac{2m-d}{2m+d}}(n)\right)\lambda^{-1},$
(43)

whose RHS is $O_{\mathbb{P}}(1)$ provided that
$n^{-2m/(2m+d)}\log^{(2m-d)/(2m+d)}(n)\lambda^{-1}=O_{\mathbb{P}}(1)$. This is
ensured by $(ii)$, as, for all large $n$,

$\displaystyle n^{-\frac{2m}{2m+d}}\log^{\frac{2m-d}{2m+d}}(n)\lambda^{-1}\leq
n^{-\frac{2m}{2m+d}}\log^{\frac{3}{2}}(n)\lambda^{-1}=O_{\mathbb{P}}(1).$

Now, since the RHS of (43) is $O_{\mathbb{P}}(1)$ and $2>4m/(2m+d)$, an
application of part A of Lemma 3 reveals that
$J_{m}(\widetilde{\beta}_{n})=O_{\mathbb{P}}(1)$ also. From (42) we now see
that
$\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}=O_{\mathbb{P}}(\log^{2}(n)n^{-m/(2m+d)})$
and consequently (26) holds.

_Case 2._ Focusing now on the second possible case, i.e., on

$\displaystyle\varepsilon\frac{\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|^{2}_{n,p}}{\log(n)D_{\varepsilon}^{2}\left\\{1+J_{m}\left(\widetilde{\beta}_{n}\right)\right\\}^{2}}+\lambda
J_{m}^{2}\left(\widetilde{\beta}_{n}\right)\leq 4\lambda
J_{m}^{2}\left(\beta_{0}\right),$

we immediately see that $\lambda
J_{m}^{2}\left(\widetilde{\beta}_{n}\right)\leq 4\lambda
J_{m}^{2}\left(\beta_{0}\right)=O_{\mathbb{P}}(\lambda)$ entails
$J_{m}(\widetilde{\beta}_{n})=O_{\mathbb{P}}(1)$, which along with the
boundedness of $J_{m}^{2}\left(\beta_{0}\right)$ leads to

$\displaystyle\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}=O_{\mathbb{P}}(\lambda^{1/2}\log^{1/2}(n))=\log^{2}(n)O_{\mathbb{P}}\left(n^{-m/(2m+d)}+\max_{1\leq
j\leq p}\operatorname{diam}^{\kappa}(A_{j})\right),$

by our assumption $(i)$ on $\lambda$. Thus, (26) is again obtained.

_Case 3._ The third possible case is

$\displaystyle\varepsilon\frac{\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|^{2}_{n,p}}{D_{\varepsilon}^{2}\log(n)\left\\{1+J_{m}\left(\widetilde{\beta}_{n}\right)\right\\}^{2}}+\lambda
J_{m}^{2}\left(\widetilde{\beta}_{n}\right)\leq
O_{\mathbb{P}}\left(\max_{1\leq j\leq
p}\operatorname{diam}^{2\kappa}(A_{j})\right).$

It is clear from this inequality that
$J_{m}(\widetilde{\beta}_{n})=O_{\mathbb{P}}(1)$ provided that $\max_{1\leq
j\leq p}\operatorname{diam}^{2\kappa}(A_{j})\lambda^{-1}=O_{\mathbb{P}}(1)$,
which is ensured by $(ii)$, as, for all large $n$,

$\displaystyle\max_{1\leq j\leq
p}\operatorname{diam}^{2\kappa}(A_{j})\lambda^{-1}\leq\max_{1\leq j\leq
p}\operatorname{diam}^{2\kappa}(A_{j})\log^{\frac{3}{2}}(n)\lambda^{-1}=O_{\mathbb{P}}(1).$

Therefore,
$\|\widetilde{\beta}_{n}-\beta_{0}\|_{n,p}=O_{\mathbb{P}}(\log^{1/2}(n)\max_{1\leq
j\leq
p}\operatorname{diam}^{\kappa}(A_{j}))=O_{\mathbb{P}}(\log^{2}(n)\max_{1\leq
j\leq p}\operatorname{diam}^{\kappa}(A_{j}))$ and (26) is again established.

_Case 4._ The fourth and final possible case is

$\displaystyle\varepsilon\frac{\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|^{2}_{n,p}}{D_{\varepsilon}^{2}\log(n)\left\\{1+J_{m}\left(\widetilde{\beta}_{n}\right)\right\\}^{2}}+\lambda
J_{m}^{2}\left(\widetilde{\beta}_{n}\right)$ $\displaystyle\leq
O_{\mathbb{P}}\left(\max_{1\leq j\leq
p}\operatorname{diam}^{\kappa}(A_{j})\right)\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p}.$
(44)

Inequality (44) immediately gives

$\displaystyle\left\|\widetilde{\beta}_{n}-\beta_{0}\right\|_{n,p}\leq
O_{\mathbb{P}}\left(\max_{1\leq j\leq
p}\operatorname{diam}^{\kappa}(A_{j})\log(n)\right)\left\\{1+J_{m}\left(\widetilde{\beta}_{n}\right)\right\\}^{2},$
(45)

which, upon plugging this back into the RHS of (44), leads to