# Robust thin-plate splines for multivariate spatial smoothing

Ioannis Kalogridis<EMAIL_ADDRESS>Department of Mathematics,
KU Leuven

###### Abstract

A novel family of multivariate robust smoothers based on the thin-plate
(Sobolev) penalty that is particularly suitable for the analysis of spatial
data is proposed. The proposed family of estimators can be expediently
computed even in high dimensions, is invariant with respect to rigid
transformations of the coordinate axes and can be shown to possess optimal
theoretical properties under mild assumptions. The competitive performance of
the proposed thin-plate spline estimators relative to their non-robust
counterpart is illustrated in a simulation study and a real data example
involving two-dimensional geographical data on ozone concentration.

###### keywords:

Robustness, spatial data, thin-plate splines, asymptotics.
MSC 2020: 62G08, 62G35, G2H11, 62G20.

††journal: Econometrics and Statistics

##  1 Introduction

Consider the problem of estimating the regression function
$f_{0}:\mathbbm{R}^{d}\to\mathbbm{R}$ from $n$ observations
$(\mathbf{x}_{i},Y_{i})\in\mathbbm{R}^{d}\times\mathbbm{R},\ i=1,\ldots,n,$
following the model

$\displaystyle Y_{i}=f_{0}(\mathbf{x}_{i})+\epsilon_{i},\quad(i=1,\ldots,n),$
(1)

where the $\epsilon_{i}$ are random errors that are commonly assumed to be
independent and identically distributed (i.i.d.) with mean zero and finite
variance, but we will be able to considerably relax these assumptions over the
course of this paper. Models of this general type arise naturally throughout
the sciences, as very often empirical data cast doubt on parametric regression
assumptions (Wood, 2017, Chapter 5).

A popular method of estimation of $f_{0}$ that is expounded by Wahba (1990)
and Green & Silverman (1994) involves restricting $f_{0}$ to the multivariate
Sobolev space of functions of order $m$, $\mathcal{H}^{m}(\mathbb{R}^{d})$,
i.e., the space of all functions whose partial derivatives of total order $m$
for some $m\in\mathbbm{N}_{+}$ are in $\mathcal{L}^{2}(\mathbbm{R}^{d})$.
Mathematically, the space $\mathcal{H}^{m}(\mathbb{R}^{d})$ is defined as

$\displaystyle\mathcal{H}^{m}(\mathbb{R}^{d})=\left\\{f:\mathbbm{R}^{d}\to\mathbbm{R},\
\frac{\partial^{m}f(x_{1},\ldots,x_{d})}{\partial x_{1}^{m_{1}}\ldots\partial
x_{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}(\mathbb{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(x_{1},\ldots,x_{d})}{\partial
x_{1}^{m_{1}}\ldots\partial x_{d}^{m_{d}}}\right)^{2}dx_{1}\ldots dx_{d}.$ (2)

Classical least squares thin-plate spline estimators are defined as minimizers
of

$\displaystyle\frac{1}{2n}\sum_{i=1}^{n}\left(Y_{i}-f(\mathbf{x}_{i})\right)^{2}+\lambda
I_{m}^{2}(f),$

over $\mathcal{H}^{m}(\mathbb{R}^{d})$, for some $\lambda>0$ that governs the
trade-off between smoothness and goodness of fit to the data. It is easy to
see that for $d=1$ the penalty $I_{m}^{2}(f)$ reduces to
$\int_{\mathbbm{R}}|f^{(m)}(x)|^{2}dx$ and with standard arguments this may be
simplified further to $\int_{\min{x}_{i}}^{\max{x}_{i}}|f^{(m)}(x)|^{2}dx$,
which gives rise to classical smoothing spline estimators, see, e.g., Wahba
(1990). Such simplifications are not valid for $d>1$, but the problem still
admits an elegant solution provided only that $2m>d$. In fact, for $2m>d$, the
solution to this variational problem may be written in terms of $n+M$ radial
and polynomial basis functions with $M=\binom{m+d-1}{d}$, see Wood (2017, p.
216) for more details. The resulting least squares thin-plate spline can be
shown to converge at a fast rate under the aforementioned assumptions on the
error term (Cox, 1984; Györfi et al., 2010) with the result that least squares
thin-plate estimators combine computational efficiency with good theoretical
properties.

An important drawback of least squares based estimators is their lack of
resistance towards atypical observations. That is, these estimators tend to be
overly attracted towards observations that do not follow the bulk of the data,
which often leads to poor explanatory or predictive performance. In order to
remedy this deficiency, this paper introduces and studies a family of
generalized (M-type) thin-plate spline estimators defined as minimizers of

$\displaystyle
L_{n}(f)=\frac{1}{n}\sum_{i=1}^{n}\rho(Y_{i}-f(\mathbf{x}_{i}))+\lambda
I_{m}^{2}(f),$ (3)

over $\mathcal{H}^{m}(\mathbb{R}^{d})$, for some convex loss function
$\rho:\mathbbm{R}\to\mathbbm{R}_{+}$. The square loss $\rho(x)=x^{2}/2$ is a
typical example of such a loss function, but the generality of the above
formulation also permits loss functions that increase less sharply as
$|x|\to\infty$ thereby providing better protection against outliers. Notable
examples include Huber, quantile and $L^{q}$ loss functions. As we discuss
below, a minimizer of (3) in $\mathcal{H}^{m}(\mathbbm{R}^{d})$ exists under
fairly general conditions and, as a result, to identify this minimizer it
suffices to restrict attention to an $(n+M)$-dimensional subspace spanned by
radial and polynomial functions. Furthermore, we show that, for well-chosen
loss functions, optimal rates of convergence may be attained without placing
any moment conditions on the errors, thereby allowing for very heavy tailed
error distributions within our theoretical analysis.

Thin-plate splines possess a number of notable advantages over their more
popular tensor product counterparts. First, the objective function in (3) only
involves one smoothing parameter irrespective of the dimension of the
predictor space. By contrast, tensor product penalties require as many
smoothing parameters as the number of predictors. Since smoothing parameters
are usually selected from the data, tensor product smoothers entail a
considerable computational burden which becomes prohibitive for higher
dimensions. Another attractive property of thin-plate splines is the
invariance of the penalty $I_{m}^{2}(f)$ to rotations or reflections of the
coordinate axes. This fact implies that thin-plate splines are ideal for
spatial/geographic smoothing as in these cases the amount of smoothing does
not depend on which axes represent the relative positions, e.g., latitude and
longitude. To these advantages we may add that tensor product smooths often
rely on the subjective choice of the number of knots and their position,
whereas with thin-plate splines the user is absolved from this responsibility,
as both of these aspects emerge naturally from the mathematical problem that
thin-plate splines solve.

In view of these advantages it comes as a surprise that thin-plate splines,
with the exception of univariate splines with $d=1$, have not received enough
attention in the literature. In fact, both available treatments of thin-plate
spline estimators of arbitrary dimensions (Cox, 1984; Györfi et al., 2010)
concern the narrow class of least squares thin-plate estimators and as a
result the larger and more versatile class of M-type thin-plate estimators has
been overlooked. To address this gap in the literature, the present paper
establishes optimal rates of convergence for a wide variety of thin-plate
estimators under a mild set of conditions considerably extending our
theoretical understanding. Furthermore, while the main emphasis of this work
is on outlier-resistant loss functions, we also demonstrate how our
methodology can be used to improve upon the results of Györfi et al. (2010) by
sharpening the rate of convergence obtained by these authors for the least
squares thin-plate estimator.

##  2 Main results: existence of solutions and rates of convergence

We begin our study of generalized M-type thin-plate estimators by providing
sufficient conditions for the existence of a minimizer in
$\mathcal{H}^{m}(\mathbbm{R}^{d})$ for the objective function $L_{n}(f)$, as
given in (3). For the special least squares case with $\rho(x)=x^{2}/2$ such
conditions have already been presented in the literature (see, e.g., p. 31,
Wahba, 1990) and involve the uniqueness of the least squares estimator on the
null-space of the penalty functional $I_{m}$. This null-space is
$M$-dimensional with $M=\binom{m+d-1}{d}$ and consists of all polynomials of
total degree at most $m$. Let $\phi_{1},\ldots,\phi_{M}$ denote any basis for
this space of polynomials. Proposition 1 below shows that the least squares
requirement for existence generalizes nicely to arbitrary convex losses.

######  Proposition 1.

Assume that $\rho$ is a convex loss function, $2m>d$ and the covariates
$\mathbf{x}_{i}\in\mathbb{R}^{d}$ are such that the corresponding unpenalized
M-estimator on $\phi_{1},\ldots,\phi_{M}$ restricted to the $\mathbf{x}_{i}$
is unique. Then, $L_{n}(f)$ has a minimizer in
$\mathcal{H}^{m}(\mathbbm{R}^{d})$.

The uniqueness requirement of Proposition 1 is very mild. For example, for
$d=2$ and $m=2$ we may take $\phi_{1}(x_{1},x_{2})=1$,
$\phi_{2}(x_{1},x_{2})=x_{1}$ and $\phi_{3}(x_{1},x_{2})=x_{2}$, so that if
the estimated plane is unique then the existence of a minimizer is guaranteed.
The condition $2m>d$ is satisfied for all $m\geq 1$ when $d=1$, i.e.,
univariate data, but precludes small values of $m$ in higher dimensions.
Unfortunately, this condition cannot be weakened as it is both necessary and
sufficient for $\mathcal{H}^{m}(\mathbbm{R}^{d})$ to be a reproducing kernel
Hilbert space and consequently for point evaluation maps
$\mathcal{H}^{m}(\mathbbm{R}^{d})\to\mathbbm{R}:f\mapsto f(\mathbf{x}),\
\mathbf{x}\in\mathbbm{R}^{d}$, to be well defined. Proposition 1 is of primary
importance not only as a first step towards the theoretical analysis of M-type
thin-plate estimators but also for their efficient computation; Section 3
provides the details.

Having established that the problem is well-defined, we now investigate the
rate of convergence of the minimizer of $L_{n}(f)$ in
$\mathcal{H}^{m}(\mathcal{R}^{d})$, which we denote by $\widehat{f}_{n}$, to
the true function $f_{0}$. The distance metric to be used is the
$\mathcal{L}^{2}(Q_{n})$-distance given by

$\displaystyle\|f-g\|_{\mathcal{L}^{2}(Q_{n})}=\left\\{\int_{\mathbbm{R}^{d}}|f(\mathbf{x})-g(\mathbf{x})|^{2}dQ_{n}(\mathbf{x})\right\\}^{1/2}$

where $Q_{n}$ is the empirical measure placing mass $n^{-1}$ at each
$\mathbf{x}_{i}\in\mathbbm{R}^{d}$. As we shall see, under very mild
conditions on the covariates, rates of convergence in terms of the empirical
norm $\|\cdot\|_{\mathcal{L}^{2}(Q_{n})}$ also translate to rates of
convergence in the classical $\mathcal{L}^{2}$-distance. It is worth noting
that, while here and in the sequel we treat the covariates as non-random, all
subsequent results also hold for random covariates provided that we condition
on them. Our theoretical development relies on the following regularity
conditions.

  1. A1.

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

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

  2. A2.

The errors $\epsilon_{i},i=1,\ldots,n$ are independent random variables and
there exists a constant $\kappa>0$ such that for all $|t|\leq\kappa$,

$\displaystyle\inf_{n}\min_{1\leq i\leq
n}\mathbb{E}\\{\rho(\epsilon_{i}+t)-\rho(\epsilon_{i})\\}\geq\kappa t^{2}$

  3. A3.

The covariates $\mathbf{x}_{i},\ i=1,\ldots,n$ are contained in a bounded open
set $\mathcal{O}\subset\mathbb{R}^{d}$ whose boundary, $\partial\mathcal{O}$,
satisfies the uniform cone condition of Adams & Fournier (2003, p. 83). That
is, there exists a locally finite open cover $\\{U_{j}\\}$ of
$\partial\mathcal{O}$ and a corresponding sequence of finite cones
$\\{C_{j}\\}$, each congruent to some fixed cone $C$, such that

     * (i)

There exists an $M<\infty$ such that every $U_{j}$ has diameter less than $M$.

     * (ii)

$\\{\mathbf{x}\in\mathcal{O}:\inf_{\mathbf{y}\in\partial\mathcal{O}}\|\mathbf{x}-\mathbf{y}\|<\delta\\}\subset\bigcup_{j=1}^{\infty}U_{j}$
for some $\delta>0$.

     * (iii)

$Q_{j}\equiv\bigcup_{\mathbf{x}\in\mathcal{O}\cap
U_{j}}(\mathbf{x}+C_{j})\subset\mathcal{O}$ for every $j$.

     * (iv)

For some finite $R$, every collection of $R+1$ of the sets $Q_{j}$ has an
empty intersection.

  4. A4.

Define the quantities

$\displaystyle h_{\max,n}$
$\displaystyle=\sup_{\mathbf{x}\in\mathcal{O}}\min_{1\leq i\leq
n}\|\mathbf{x}-\mathbf{x}_{i}\|$ $\displaystyle h_{\min,n}$
$\displaystyle=\min_{i\neq j}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|.$

For all large $n$, there exist finite positive constants $B_{1},B_{2}$ such
that $h_{max,n}\leq B_{1}$ and $h_{\max,n}/h_{\min,n}\leq B_{2}$.

Assumption A1 is very general and a large number of loss functions fulfils
these conditions. The convexity requirement implies that $\rho$ is minimally
continuous, but, importantly, differentiability is not needed for Theorem 1
below to hold. Thus, this assumption also covers non-smooth loss functions,
such as the quantile loss $\rho_{\alpha}(x)=x(\alpha-\mathcal{I}(x<0)),\
\alpha\in(0,1)$. Assumption A2 requires the local quadratic behaviour of
$m_{i}(t):=\mathbb{E}\\{\rho(\epsilon_{i}+t)\\}$, $i=1,\ldots,n$ about zero.
Assumptions of this type have been extensively used in the asymptotics of
M-estimators (see, e.g., van der Vaart & Wellner, 1996, Theorem 3.2.5, p.
289). It is similarly a weak condition that is satisfied quite generally. To
see this, observe that by Fubini’s theorem

$\displaystyle
m_{i}(t)-m_{i}(0)=\int_{0}^{t}\mathbb{E}\\{\rho^{\prime}(\epsilon_{i}+x)\\}dx,$

where $\rho^{\prime}$ is any subgradient of $\rho$, so that all examples of
loss functions given by Kalogridis (2021) also satisfy our A2 under the
conditions given by that author. Our assumptions do not entail identically
distributed errors, as in our experience this is too strong of an assumption
for many practical settings.

Conditions A3 and A4 concerning the design points (knots) and their positions
have been previously used in the asymptotics of least squares thin-plate
spline estimators, see Utreras (1988). The uniform cone condition precludes
sets with very irregular boundaries, but is satisfied quite generally
otherwise; it is satisfied, e.g., by balls and rectangles. Condition A4 is
also very modest, as it essentially requires the design points to be distinct
and dense within the domain of interest, at least for large $n$. Both of these
requirements follow from the fact that the ratio $h_{\max,n}/h_{\min,n}$
remains bounded for all large $n$, as if either the observations are not
unique or dense enough as $n\to\infty$, $h_{\max,n}/h_{\min,n}$ could become
unbounded.

With these assumptions we can now state our first asymptotic result. It is
worth noting that in Theorem 1 below we treat the smoothing parameter
$\lambda$ as a random variable possibly depending on our sample
$(\mathbf{x}_{1},Y_{1}),\ldots,(\mathbf{x}_{n},Y_{n})$. We place no
restrictions on the way that $\lambda$ may depend on the sample and instead
merely require that $\lambda$ is a measurable function of our sample so that
$\lambda$ is a properly defined random variable. Our treatment of $\lambda$
needs to be contrasted with the treatment of the smoothing parameter by other
authors, e.g., Cox (1984); Utreras (1988), who regard it as a deterministic
sequence. In our view, our treatment constitutes an important extension of
existing results as $\lambda$ is most often selected from the data and is thus
random rather than fixed, see, e.g., the data-driven method of selecting
$\lambda$ presented in Section 3.

######  Theorem 1.

Assume that the conditions of Proposition 1 are met, A1–A4 hold and further
that $\lambda=O_{P}(n^{-2m/(2m+d)})$ as well as
$\lambda^{-1}=O_{P}(n^{2m/(2m+d)})$. Then, there exists a sequence,
$\widehat{f}_{n}$, of M-type thin-plate splines minimizing (3) such that

$\displaystyle\|\widehat{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}^{2}=O_{P}(n^{-2m/(2m+d)})\quad\text{and}\quad
I_{m}(\widehat{f}_{n})=O_{P}(1),$

as $n\to\infty$.

The limit conditions involving the smoothing parameter require that $\lambda$
tends to zero in probability as $n\to\infty$, but not too fast. The rate of
decay $n^{-2m/(2m+d)}$ for $\lambda$ ensures that the asymptotic variance and
bias of the estimator are balanced and this leads to the rate of convergence
$n^{-2m/(2m+d)}$ for thin-plate estimators. It is worth noting that
$n^{-2m/(2m+d)}$ is the optimal (squared) rate of convergence for functions in
$\mathcal{H}^{m}(\mathcal{O})$ (Stone, 1982) and therefore it cannot be
improved, except in trivial cases. For $d=1$ our rate of convergence is in
agreement with the rate obtained by Kalogridis (2021) for univariate smoothing
spline estimators, but for $d>1$ the result in Theorem 1 constitutes an
important generalization. The obtained rate of convergence suggests that thin-
plate spline estimators, like most nonparametric estimators, suffer from the
curse of dimensionality and consequently, for given sample size $n$,
estimation becomes less and less precise for larger predictor dimension $d$.

Theorem 1 establishes not only a rate of convergence in the empirical norm
$\|\cdot\|_{\mathcal{L}^{2}(Q_{n})}$, but also the boundedness of the semi-
norm $I_{m}(\widehat{f}_{n})$. The latter is crucial in extending this rate of
convergence to the $\mathcal{L}^{2}$-norm as well as in establishing optimal
rates of convergence of certain useful derivatives. These extensions are given
in Corollary 1 below.

######  Corollary 1.

Suppose that the assumptions of Theorem 1 hold and that
$h_{\max,n}=O(n^{-1/(2m+d)})$. Then, the sequence of minimizers of (3),
$\widehat{f}_{n}$, satisfies

$\displaystyle\int_{\mathcal{O}}|\widehat{f}_{n}(\mathbf{x})-f_{0}(\mathbf{x})|^{2}d\mathbf{x}=O_{P}(n^{-2m/(2m+d)}),$

and for every $(j_{1},\ldots,j_{d})\in\mathbbm{R}^{d}$ such that
$j_{1}+\ldots+j_{d}=j\leq m$

$\displaystyle\int_{\mathcal{O}}\left|\frac{\partial\widehat{f}_{n}^{j}(\mathbf{x})}{\partial
x_{1}^{j_{1}}\ldots\partial x_{d}^{j_{d}}}-\frac{\partial
f^{j}_{0}(\mathbf{x})}{\partial x_{1}^{j_{1}}\ldots\partial
x_{d}^{j_{d}}}\right|^{2}d\mathbf{x}=O_{P}(n^{-2(m-j)/(2m+d)}).$

The above rates of convergence are again optimal according to the results of
Stone (1982). To the best of our knowledge, these are the first convergence
results for derivatives of thin-plate estimators even in the relatively simple
least squares case.

It is interesting to compare the rate of convergence $n^{-2m/(2m+d)}$ obtained
herein with the rate $\log(n)n^{-2m/(2m+d)}$ for the least squares thin-plate
spline estimator emerging from the results of Györfi et al. (2010, Chapter
21). The square loss $\rho(x)=x^{2}/2$ is not covered by our previous set of
assumptions (it is not Lipschitz), but least squares estimators may be easily
treated on a separate basis. In particular, letting $\widehat{f}_{n}$ now
denote the least squares estimator, one may easily verify the inequality

$\displaystyle\|\widehat{f}_{n}-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}+\lambda
I_{m}^{2}(\widehat{f}_{n})\leq
2\langle\epsilon,\widehat{f}_{n}-f_{0}\rangle_{\mathcal{L}^{2}(Q_{n})}+\lambda
I_{m}^{2}(f_{0}),$ (4)

where $\langle\epsilon,\widehat{f}_{n}-f_{0}\rangle_{\mathcal{L}^{2}(Q_{n})}$
stands for
$n^{-1}\sum_{i=1}^{n}\epsilon_{i}(\widehat{f}_{n}(\mathbf{x}_{i})-f_{0}(\mathbf{x}_{i}))$.
This inequality combined with our sharper estimate for the local entropy
(compare our Lemma 1 below with Lemma 20.6 in Györfi et al. (2010)) and the
modulus of continuity of the empirical process derived in van de Geer (2000,
Chapter 10) now leads to Theorem 2.

######  Theorem 2.

Assume that the conditions of Proposition 1 are met, A3 holds and that the
errors $\epsilon_{i}$ are uniformly sub-Gaussian, i.e., there exist finite
constants $K_{1}$ and $K_{2}>0$ such that

$\displaystyle\sup_{n}\max_{1\leq i\leq
n}K_{1}^{2}\mathbb{E}\\{e^{|\epsilon_{i}|^{2}/K_{1}^{2}}-1\\}\leq K_{2}.$

If $\lambda=O_{P}(n^{-2m/(2m+d)})$ and $\lambda^{-1}=O_{P}(n^{2m/(2m+d)})$,
then there exists a sequence of least squares thin-plate spline estimators,
$\widehat{f}_{n}$, satisfying

$\displaystyle\|\widehat{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}^{2}=O_{P}(n^{-2m/(2m+d)})\quad\text{and}\quad
I_{m}(\widehat{f}_{n})=O_{P}(1),$

as $n\to\infty$.

It is important to emphasize that the sub-Gaussian requirement in Theorem 2 is
used exclusively for the treatment of the least-squares thin-plate estimator
and not for the robust estimators for the treatment of which we rely solely on
assumptions A1–A4. The sub-Gaussian requirement is met in practice, e.g.,
whenever the errors follow a Gaussian distribution or, more generally,
whenever they possess a squared exponential moment. As noted previously,
Theorem 2 improves upon the corresponding result of Györfi et al. (2010).
Rates of convergence for the least squares estimator in the
$\mathcal{L}^{2}(\mathcal{O})$-norm as well as rates of convergence for the
derivatives may now be established exactly as in the proof of Corollary 1; we
omit the details.

##  3 Computation and smoothing parameter selection

As hinted previously, Proposition 1 is crucial not only from a theoretical,
but also from a practical standpoint. In fact, with the help of Proposition 1
and reasoning along the same lines as in the proof of Green & Silverman (1994,
Theorem 7.3), it can be shown that in order to identify the minimizer of
$L_{n}(f)$ in (3) it suffices to restrict attention to functions of the form

$\displaystyle
f(\mathbf{x})=\sum_{i=1}^{n}\gamma_{i}\eta_{m,d}(\|\mathbf{x}-\mathbf{x}_{i}\|)+\sum_{j=1}^{M}\delta_{j}\phi_{j}(\mathbf{x}),$
(5)

where $\eta_{m,d}:\mathbb{R}_{+}\to\mathbb{R}$ is 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}$

with $\Gamma(\cdot)$ denoting Euler’s gamma function. The coefficient vector
$\boldsymbol{\gamma}$ is subject to the set of linear constraints
$\boldsymbol{\Phi}^{\top}\boldsymbol{\gamma}=\mathbf{0}$, where the $n\times
M$ matrix $\boldsymbol{\Phi}$ has $(i,j)$th entry $\phi_{j}(\mathbf{x}_{i})$.
In other words, the coefficient vector $\boldsymbol{\gamma}\in\mathbbm{R}^{n}$
needs to be perpendicular to the $M$-dimensional space spanned by the
restriction of the polynomial functions $\phi_{1},\ldots,\phi_{M}$ to the
design points.

The representation in (5) as well as the set of linear constraints are derived
from reproducing kernel Hilbert space arguments that split the space
$\mathcal{H}^{m}(\mathbb{R}^{d})$ into the closed finite-dimensional null
space of $I_{m}(f)$ and its orthogonal complement. It follows from these
arguments that plugging (5) into (3) yields the quadratic form
$I_{m}^{2}(f)=\boldsymbol{\gamma}^{\top}\mathbf{\Omega}\boldsymbol{\gamma}$
where the $n\times n$ matrix $\mathbf{\Omega}$ has $(i,j)$th entry
$\eta_{m,d}(\|\mathbf{x}_{i}-\mathbf{x}_{j}\|)$ and the minimization problem
becomes

$\displaystyle\min_{(\boldsymbol{\gamma},\boldsymbol{\delta})\in\mathbbm{R}^{n}\times\mathbbm{R}^{M}}\left[\frac{1}{n}\sum_{i=1}^{n}\rho(Y_{i}-\boldsymbol{\omega}_{i}^{\top}\boldsymbol{\gamma}-\boldsymbol{\phi}_{i}^{\top}\boldsymbol{\delta})+\lambda\boldsymbol{\gamma}^{\top}\mathbf{\Omega}\boldsymbol{\gamma}\right]$
(6) $\displaystyle\quad\quad\text{s.t.}\
\boldsymbol{\Phi}^{\top}\boldsymbol{\gamma}=\mathbf{0},$

with $\boldsymbol{\omega}_{i}\in\mathbbm{R}^{n}$ and
$\boldsymbol{\phi}_{i}\in\mathbbm{R}^{M}$ denoting row vectors of
$\mathbf{\Omega}$ and $\boldsymbol{\Phi}$, respectively. It can be shown that
for all $\boldsymbol{\gamma}\in\mathbbm{R}^{n}$ satisfying
$\boldsymbol{\Phi}^{\top}\boldsymbol{\gamma}=\mathbf{0}$ we have
$\boldsymbol{\gamma}^{\top}\mathbf{\Omega}\boldsymbol{\gamma}>0$, see Wahba
(1990, pp. 32–33) and the references therein. We may automatically incorporate
the linear constraints and further simplify the problem by putting
$\boldsymbol{\gamma}=\mathbf{Q}\boldsymbol{\beta}$ for a matrix $\mathbf{Q}$
whose columns span the null space of $\boldsymbol{\Phi}^{\top}$ and
$\boldsymbol{\beta}\in\mathbbm{R}^{n-M}$. The matrix $\mathbf{Q}$ can be
easily obtained through the QR or singular value decompositions of
$\boldsymbol{\Phi}$. With this simplification, the solution to (LABEL:eq:6)
may be identified using, for example, the variant of the iteratively
reweighted least squares (IRLS) algorithm employed by Kalogridis (2021) with
the radial basis functions and polynomials in (5) taking the place of B-spline
basis functions there.

To determine the penalty parameter $\lambda$ in a data-driven way, we propose
a suitable adaptation of the strategy of Maronna (2011). Let
$\mathbf{r}_{-}=(r_{-1},\ldots,r_{-n})^{\top}$ denote an approximation to the
leave-one out residuals, as obtained, for example, from the last step of the
IRLS algorithm. We propose to select the value of $\lambda$ that minimizes the
robust cross-validation (RCV) criterion

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

where $\tau$ denotes the robust and efficient $\tau$-scale introduced by Yohai
and Zamar (1988) with tuning constants equal to $c_{1}=3$ and $c_{2}=5$,
corresponding to the biweighting of the mean and standard deviation
respectively. This criterion may be viewed as a robustification of the
celebrated leave-one-out criterion (see, e.g., Wahba, 1990, pp. 47–52) in
which the $\tau$-scale is replaced by the mean of squares of the $r_{-i}$.
Thus, while in the classical leave-one-out criterion all the $|r_{-i}|^{2}$
contribute equally (with weight $n^{-1}$ each), the use of the robust
$\tau$-scale employed herein reduces the effect of large $|r_{-i}|^{2}$ on the
selection of $\lambda$ thereby leading to a robust automatic selection
procedure.

##  4 A Monte Carlo study

We now examine the practical performance of several thin-plate spline
estimators by means of a simulation study. We are interested in the
performance of the competing estimators not only in completely regular data,
but also in data that may contain a number of outlying observations. The
estimators to be considered are

  * 1.

The classical least squares thin-plate spline estimator, abbreviated as LS.

  * 2.

The least absolute deviations type thin-plate spline estimator with loss
function $\rho(x)=|x|$ abbreviated as LAD.

  * 3.

The Huber type thin-plate spline estimator with loss function

$\displaystyle\rho(x)=\begin{cases}x^{2}/2&|x|<1.345\\\
1.345|x|-1.345^{2}/2&|x|\geq 1.345\end{cases},$

abbreviated as Huber.

  * 4.

The logistic type thin-plate spline estimator with loss function

$\displaystyle\rho(x)=2x+4\log(1+e^{-x}),$

abbreviated as Logistic.

To compare the above estimators we generate observations from the model

$\displaystyle Y_{i}=f_{0}(\mathbf{x}_{i})+\epsilon_{i},\quad(i=1,\ldots,n),$

where the regression function $f_{0}(\mathbf{x})$ is either given by
$f_{1}(x_{1},x_{2})=\exp[-8|x_{1}-0.5|^{2}-8|x_{2}-0.5|^{2}]$,
$f_{2}(x_{1},x_{2})=\sin(2\pi x_{1})\cos(2\pi x_{2})$ or
$f_{3}(x_{1},x_{2},x_{3})=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}$,
$(x_{i1},x_{i2},x_{i3})$ are uniformly distributed random variables in the
unit square $(0,1)^{3}$ and the errors $\epsilon_{i}$ are generated according
to the following distributions: (i) the standard Gaussian distribution, (ii)
the t-distribution with $3$ degrees of freedom, (iii) the skewed
t-distribution with $3$ degrees of freedom and non-centrality parameter equal
to $1$ leading to right-skewed data, (iv) a mixture Gaussian distribution with
weights means equal to $1$ and $10$, variances equal to $1$ and $0.01$ and
weights equal to $0.85$ and $0.15$ respectively and (v) Tukey’s Slash
distribution defined as the distribution of the quotient of independent
standard normal and uniform random variables. We assess the performance of the
competing estimators through the mean-squared error (MSE) given by

$\displaystyle\operatorname{MSE}=\frac{1}{n}\sum_{i=1}^{n}|\widehat{f}_{n}(\mathbf{x}_{i})-f_{0}(\mathbf{x}_{i})|^{2},$

which is an approximation to the squared $\mathcal{L}^{2}$-error. Table 1
below reports average MSEs and their standard errors obtained from 1000
replications with datasets of size $n=100$.

|  | LS | Huber | Logistic | LAD
---|---|---|---|---|---
$f_{0}$ | Dist. | Mean | SE | Mean | SE | Mean | SE | Mean | SE
$f_{1}$ | Gaussian | 7.17 | 0.11 | 7.34 | 0.10 | 7.39 | 0.10 | 11.76 | 0.14
$t_{3}$ | 15.11 | 0.36 | 9.52 | 0.13 | 9.89 | 0.15 | 13.01 | 0.19
$st_{3,1}$ | 74.04 | 0.87 | 52.26 | 0.49 | 56.92 | 0.55 | 48.76 | 0.52
| M. Gaussian | 292.2 | 4.22 | 21.49 | 0.40 | 36.41 | 0.68 | 21.06 | 0.39
| Slash | 9e+07 | 3e+06 | 25.09 | 0.48 | 27.16 | 0.56 | 27.80 | 0.56
$f_{2}$ | Gaussian | 14.20 | 0.16 | 14.17 | 0.17 | 13.88 | 0.17 | 21.04 | 0.23
$t_{3}$ | 28.25 | 0.45 | 21.98 | 0.37 | 22.02 | 0.30 | 24.44 | 0.27
$st_{3,1}$ | 84.83 | 0.90 | 62.13 | 0.55 | 66.81 | 0.58 | 59.13 | 0.57
| M. Gaussian | 311.1 | 4.18 | 39.12 | 1.23 | 57.14 | 1.50 | 34.94 | 0.89
| Slash | 4e+07 | 2e+06 | 72.19 | 2.63 | 56.70 | 1.38 | 57.29 | 0.95
$f_{3}$ | Gaussian | 7.72 | 0.27 | 7.77 | 0.26 | 7.73 | 0.25 | 15.07 | 0.28
$t_{3}$ | 18.19 | 0.75 | 11.52 | 0.43 | 12.01 | 0.43 | 18.43 | 0.33
$st_{3,1}$ | 74.37 | 1.05 | 53.10 | 0.57 | 57.45 | 0.60 | 52.47 | 0.56
| M. Gaussian | 321.90 | 6.35 | 25.89 | 0.72 | 43.02 | 1.47 | 25.09 | 0.50
| Slash | 4e+07 | 5e+06 | 45.83 | 1.40 | 48.63 | 1.73 | 43.66 | 0.87
Table 1: Means and standard errors of 1000 MSEs $(\times 100)$ with $n=100$ of
the least squares, Huber, logistic and least absolute deviations type thin-
plate spline estimators.

The results in Table 1 indicate the sensitivity of least squares estimator to
departures from the (sub-)Gaussian assumption. In particular, even with
$t_{3}$ errors that lead to a few mild outliers, the least squares estimator
significantly loses ground compared to the Huber, Logistic and LAD estimators.
The former two estimators almost match the performance of the least squares
estimator in regular data but also exhibit a high degree of resistance towards
gross errors. The LAD estimator is inefficient relative to its competitors in
the situation of light-tailed Gaussian errors, but exhibits superior
performance to the Huber and Logistic thin-plate spline estimators in
situations of heavy-tailed asymmetric contamination, such as contamination
incurred by $st_{3,1}$ and mixture Gaussian errors.

To illustrate the key practical differences between the least squares thin-
plate spline and its robust counterparts, Figures 1 and 2 present two examples
of estimated surfaces by the least squares and least absolute deviations
estimators under Gaussian and mixture Gaussian errors, respectively. The
regression function $f_{1}$ has a unit-sized bump at $(0.5,0.5)$, which both
estimators get right in the absence of outliers, as evidenced in Figure 1.
Figure 2, however, indicates that least squares estimates can be heavily
distorted by the presence of outlying observations to the extent that the
estimated surface bears no resemblance to the true surface. Despite the heavy
contamination, the bump is still visible in the right panel of Figure 2, which
depicts the LAD estimated surface.

(a)  (b)

Figure 1: Typical estimated surfaces for
$f_{1}(x_{1},x_{2})=\exp[-8|x_{1}-0.5|^{2}-8|x_{2}-0.5|^{2}]$ under Gaussian
errors by the least squares and least absolute deviations estimators on the
left and right panels, respectively.

##  5 Application: Ozone Levels in Midwestern USA

While stratospheric ozone protects living organisms from ultraviolet radiation
from the sun, high ground-level concentrations of ozone can trigger a variety
of health problems, particularly for people with breathing difficulties,
children and the elderly. In this example we examine the concentration of
ground-level in the midwestern US as a function of geographical longitude and
latitude. That is, we consider the model

$\displaystyle\text{Ozone}_{i}=f_{0}(\text{Longitude}_{i},\text{Latitude}_{i})+\epsilon_{i},\quad(i=1,\ldots,n).$

The data for this analysis consists of $8$-hour average surface ozone from 9am
to 4pm in parts per billion (PPB) from 147 sites in midwestern US on July 3,
1987. The geographical coordinates of these sites constitute the predictor
variables. The present dataset is part of a much larger dataset which is
freely available as part of the fields package (Nychka et al., 2017) in CRAN
(R Core Team, 2022).

(a)  (b)

Figure 2: Typical estimated surfaces for
$f_{1}(x_{1},x_{2})=\exp[-8|x_{1}-0.5|^{2}-8|x_{2}-0.5|^{2}]$ under mixture
Gaussian errors by the least squares and least absolute deviations estimators
on the left and right panels, respectively.

Since ozone concentration tends to be a skewed random variable with several
potentially outlying values, we have estimated the regression function $f_{0}$
with both the least squares (LS) and least absolute deviations (LAD) thin-
plate spline estimators. The contours of the estimated surfaces on the convex
hull of the data are depicted in the left and right panels of Figure 3,
respectively. The panels of the figure suggest that while there is broad
agreement between these estimated surfaces around Indianapolis, the surfaces
are noticeably different on the central and western part of the data. In
particular, in the areas west and south of Milwaukee and St. Louis, LS
estimates tend to underestimate ozone concentrations relative to LAD
estimates. These differences are probably more consequential when it comes to
concentrations in the excess of 60 PPB as these are more detrimental to one’s
health.

In view of the lack of resistance of LS estimators, it may be conjectured that
the large differences between the estimates depicted in Figure 3 are
attributable to the presence of atypical observations within the data. Since
robust regression estimators are less attracted to outlying observations, such
observations result in large absolute residuals which we may then use for
their detection. A popular rule of thumb in that respect involves flagging the
$i$th observation as an outlier if its standardized absolute LAD residual
$r_{i}/\operatorname{MAD}(\mathbf{r})$ is larger than $2.5$. This rule results
in the detection of 14 outlying observations for the LAD estimator, but only 7
for the non-robust LS estimator.

(a)  (b)

Figure 3: Contours of the estimated surfaces. Left: contours of the least
squares thin-plate spline estimator. Right: contours of the least absolute
deviations thin-plate spline estimator. Darker colors indicate higher ozone
concentrations. The observation sites are depicted with solid black dots and
green rhombuses depending on whether the observation is classified as an
outlier or not.

(a)  (b)

Figure 4: Contours of the estimated surfaces after removing outlying
observations from the data. Left: contours of the least squares thin-plate
spline estimator. Right: contours of the least absolute deviations thin-plate
spline estimator. Darker colors indicate higher ozone concentrations. The
observation sites are depicted with solid black dots.

To demonstrate precisely how sensitive the LS estimator is towards outlying
observations, we have removed the LAD-detected outliers from the dataset and
recomputed the estimates. This results in the left and right panels in Figure
4 for the LS and LAD estimates, respectively. It is interesting to observe
that without these outliers the LAD and LS estimates are much closer to one
another. Moreover, comparing the right panels of Figure 3 and Figure 4 reveals
that, contrary to LS estimates, LAD estimates have undergone minimal change
after removing the outliers from the data. LAD estimates are thus better able
to describe the bulk of the data than their sensitive LS counterparts.

##  6 Conclusion

The present paper introduces robust estimators for multivariate nonparametric
regression models based on the highly advantageous thin-plate penalty. The
proposed class of estimators enjoys optimal theoretical properties under mild
assumptions and can be expediently computed even in high dimensions. There are
several research directions worth exploring from here. Our theoretical
treatment relies on a specific rate of decay of $\lambda$ and it is not known
whether $\lambda$ attains this rate of decay when determined by our robust
cross validation procedure. In order to ascertain whether this is the case, a
close investigation of robust model selection procedures can be worthwhile.
Additionally, our treatment rests upon the convexity of the loss function and
hence the important class of redescending thin-plate estimators is left out.
Future effort can therefore be directed towards a dedicated treatment of these
estimators.

An important generalization of our ideas would involve robust estimation of
multivariate nonparametric generalized linear models, where the distribution
of the response variable can be any member of the exponential family, thus
significantly extending the range of applications. Such an estimator may be
based on, for example, the thin-plate penalty proposed herein and the density
power divergence, as used by Kalogridis et al. (2023) in the univariate
setting.

Another important area of research where thin-plate splines are likely to be
successful is functional data analysis and, in particular, location and
dispersion estimation from discretely sampled functional data. For the
location case, thin-plate splines may be used to extend the robust estimator
of Kalogridis and Van Aelst (2022) for discretely sampled functional data on a
bounded interval to discretely sampled functional data on much more
complicated multivariate domains, such as a sphere. For dispersion estimation,
thin-plate splines can be combined with resistant loss functions and provide a
potent alternative to estimators based on tensor product penalties, see, e.g.,
(Hsing & Eubank, 2015, Chapter 8). We aim to study these important
generalizations as part of our future work.

## Appendix: Proofs of the theoretical results

###### Proof of Proposition 1.

The result would follow by direct application of Theorem 3.2 of Cox &
O’Sullivan (1985) provided that we can check the conditions of that theorem.
Since $I_{m}^{2}$ is a squared semi-norm on $\mathcal{H}^{m}(\mathbbm{R}^{d})$
and $I_{m}$ has a finite $M$-dimensional null space, these conditions entail
the convexity and weak lower semicontinuity of the map

$\displaystyle\mathcal{H}^{m}(\mathbbm{R}^{d})\to\mathbb{R}_{+}:f\mapsto
n^{-1}\sum_{i=1}^{n}\rho(Y_{i}-f(\mathbf{x}_{i})).$

We only need to check weak lower continuity, because convexity follows easily
from the convexity of $\rho$. Convexity also implies that we only need to
establish the lower semicontinuity of the map, as by the Hahn-Banach theorem
(Rynne & Youngston, 2008) convexity and lower semicontinuity imply weak lower
semicontinuity. To that end, let $\\{f_{k}\\}_{k}$ denote a sequence in
$\mathcal{H}^{m}(\mathbbm{R}^{d})$ converging to some $f^{\star}$. By the
Sobolev embedding theorem (Adams & Fournier, 2003, Theorem 4.12), we have

$\displaystyle\max_{1\leq i\leq
n}|f_{k}(\mathbf{x}_{i})-f^{\star}(\mathbf{x}_{i})|\leq
c_{0}\|f_{k}-f^{\star}\|_{\mathcal{H}^{m}(B_{r_{n}}^{d}(0))}\leq
c_{0}\|f_{k}-f^{\star}\|_{\mathcal{H}^{m}(\mathbb{R}^{d})},$

for some $c_{0}>0$ not depending on $k$, where $B_{r_{n}}^{d}(0)$ denotes the
smallest ball in $\mathbbm{R}^{d}$ containing all the $\mathbf{x}_{i}$ and
$\|\cdot\|_{\mathcal{H}^{m}(B_{r_{n}}^{d}(0))},\|\cdot\|_{\mathcal{H}^{m}(\mathbb{R}^{d})}$
denote the standard Sobolev norms on $\mathcal{B}_{r_{n}}^{d}(0)$ and
$\mathbbm{R}^{d}$, respectively. Letting $k\to\infty$, we find that
$\max_{1\leq i\leq n}|f_{k}(\mathbf{x}_{i})-f^{\star}(\mathbf{x}_{i})|\to 0$.
The continuity of $\rho$ concludes the proof. ∎

We now introduce some useful notation that will be used in the proof of
Theorem 1 and Theorem 2. Let $Q_{n}$ denote the probability measure on
$(\mathbb{R}^{d},\mathcal{B}(\mathbbm{R}^{d}))$ given by

$\displaystyle Q_{n}(A)=\frac{1}{n}\sum_{i=1}^{n}\mathcal{I}(\mathbf{x}_{i}\in
A),\ A\in\mathcal{B}(\mathbbm{R}^{d}),$

where $\mathcal{I}(\cdot)$ denotes the indicator function. Let $\mathcal{F}$
denote a class of real-valued functions on some domain
$\mathcal{X}\subset\mathbb{R}^{d}$. The $\delta$-entropy for $\mathcal{F}$ in
the $\mathcal{L}^{2}(Q_{n})$-norm,
$H(\delta,\mathcal{F},\mathcal{L}^{2}(Q_{n}))$, is defined as the logarithm of
the smallest number $N$ for which there exists a collection of functions
$f_{1},\ldots,f_{N}$ such that for every $f\in\mathcal{F}$ there exists a
$j=j(f)\in\\{1,\ldots,N\\}$ with the property

$\displaystyle\left\\{\int_{\mathcal{X}}|f(\mathbf{x})-f_{j}(\mathbf{x})|^{2}dQ_{n}(\mathbf{x})\right\\}^{1/2}\leq\delta.$

Similarly, we define the $\delta$-entropy with respect to the supremum norm,
$H_{\infty}(\delta,\mathcal{F})$, as the logarithm of the smallest $N$ for
which there exists $f_{1},\ldots,f_{N}$ such that for every $f\in\mathcal{F}$
there is a $j=j(f)$ with the property

$\displaystyle\sup_{\mathbf{x}\in\mathcal{X}}|f(\mathbf{x})-f_{j}(\mathbf{x})|\leq\delta.$

It is clear that $H(\delta,\mathcal{F},\mathcal{L}^{2}(Q_{n}))\leq
H_{\infty}(\delta,\mathcal{F})$ for all $\delta>0$. To lighten the notation in
all our proofs below we will use $c_{0}$ to denote generic constants. Thus,
the value of $c_{0}$ may change from appearance to appearance.

######  Lemma 1.

Assume A3 and A4 and that $2m>d$. Then, there exists a universal constant
$c_{0}$, independent of $n$, such that

$\displaystyle
H(\delta,\\{f\in\mathcal{H}^{m}(\mathbb{R}^{d}),\|f\|_{\mathcal{L}^{2}(Q_{n})}\leq
1,I_{m}(f)\leq M\\},\mathcal{L}^{2}(Q_{n})\\}\leq
c_{0}\left(\frac{M}{\delta}\right)^{d/m},\quad\delta>0,\ M\geq 1.$

###### Proof.

Observe that by A3, $Q_{n}$ is restricted to the open set $\mathcal{O}$ and
for any $A\in\mathcal{B}(\mathbbm{R}^{d})$ we have
$Q_{n}(A)=Q_{n}(A\cap\mathcal{O})$. Therefore, for any
$(f,g)\in\mathcal{H}^{m}(\mathbb{R}^{d})\times\mathcal{H}^{m}(\mathbb{R}^{d})$
we find

$\displaystyle\int_{\mathbb{R}^{d}}|f(\mathbf{x})-g(\mathbf{x})|^{2}dQ_{n}(\mathbf{x})$
$\displaystyle=\int_{\mathcal{O}}|f(\mathbf{x})-g(\mathbf{x})|^{2}dQ_{n}(\mathbf{x})$
$\displaystyle\leq\sup_{\mathbf{x}\in\mathcal{O}}|f(\mathbf{x})-g(\mathbf{x})|^{2}.$

Hence, to bound $H(\delta,\mathcal{F},\mathcal{L}^{2}(Q_{n})\\}$ it suffices
to obtain a bound on

$\displaystyle
H_{\infty}(\delta,\\{f\in\mathcal{H}^{m}(\mathcal{O}),\|f\|_{\mathcal{L}^{2}(Q_{n})}\leq
1,I_{m}(f)\leq M\\}),\quad\delta>0,\ M\geq 1.$ (7)

To bound (7) we start by showing the following inclusion

$\displaystyle\\{f\in\mathcal{H}^{m}(\mathcal{O}),\|f\|_{\mathcal{L}^{2}(Q_{n})}\leq
1,I_{m}(f)\leq
M\\}\subset\\{f\in\mathcal{H}^{m}(\mathcal{O}):\|f\|_{\mathcal{H}^{m}(\mathcal{O})}\leq
c_{0}M\\}$

for some $c_{0}>0$ and all $M\geq 1$, where
$\|\cdot\|_{\mathcal{H}^{m}(\mathcal{O})}$ is the Sobolev norm given by

$\displaystyle\|f\|_{\mathcal{H}^{m}(\mathcal{O})}=\left\\{\int_{\mathcal{O}}|f(\mathbf{x})|^{2}d\mathbf{x}+\sum_{m_{1}+\ldots+m_{d}=m}\binom{m}{m_{1},\ldots,m_{d}}\int_{\mathcal{O}}\left|\frac{\partial^{m}f(\mathbf{x})}{\partial
x_{1}^{m_{1}}\ldots\partial
x_{d}^{m_{d}}}\right|^{2}d\mathbf{x}\right\\}^{1/2}.$

To establish this inclusion, take an $f\in\mathcal{H}^{m}(\mathcal{O})$ such
that $\|f\|_{\mathcal{L}^{2}(Q_{n})}\leq 1$ and $I_{m}(f)\leq M$ for some
$M\geq 1$. Our assumptions imply those of Theorem 3.4 in Utreras (1988), hence
an application of that theorem reveals the existence of a constant $c_{0}$
such that

$\displaystyle\int_{\mathcal{O}}|f(\mathbf{x})|^{2}d\mathbf{x}$
$\displaystyle\leq
c_{0}\int_{\mathcal{O}}|f(\mathbf{x})|^{2}dQ_{n}(\mathbf{x})$
$\displaystyle\quad+c_{0}\sum_{m_{1}+\ldots+m_{d}=m}\binom{m}{m_{1},\ldots,m_{d}}\int_{\mathcal{O}}\left|\frac{\partial^{m}f(\mathbf{x})}{\partial
x_{1}^{m_{1}}\ldots\partial x_{d}^{m_{d}}}\right|^{2}d\mathbf{x}$
$\displaystyle\leq c_{0}\left\\{1+I_{m}^{2}(f)\right\\}$ $\displaystyle\leq
c_{0}M^{2},$

where the last inequality follows from our assumption that $M\geq 1$. With
this bound we now obtain

$\displaystyle\|f\|_{\mathcal{H}^{m}(\mathcal{O})}\leq\left\\{c_{0}M^{2}+I_{m}^{2}(f)\right\\}^{1/2}\leq
c_{0}M,$

as claimed.

The final step of the proof is a bound on

$\displaystyle
H_{\infty}\left(\delta,\\{f\in\mathcal{H}^{m}(\mathcal{O}):\|f\|_{\mathcal{H}^{m}(\mathcal{O})}\leq
c_{0}M\\}\right),\ \delta>0,M\geq 1.$

But this is the entropy of the closed $c_{0}M$-ball and Proposition 6 of
Cucker & Smale (2001) implies the existence of a universal $c_{0}$ such that

$\displaystyle
H_{\infty}\left(\delta,\\{f\in\mathcal{H}^{m}(\mathcal{O}):\|f\|_{\mathcal{H}^{m}(\mathcal{O})}\leq
c_{0}M\\}\right)\leq c_{0}\left(\frac{M}{\delta}\right)^{d/m},\ \delta>0,M\geq
1.$

The proof is complete.

∎

###### Proof of Theorem 1.

The proof of the theorem employs the convexity argument of van de Geer (2002)
combined with the Sobolev embedding theorem in order to localize the behaviour
of the objective function around $f_{0}$. A tight bound on the asymptotic
variance is established with the help of an exponential inequality based on
the improved entropy estimates obtained in Lemma 1.

Write $L_{n}(f)=M_{n}(f)+\lambda I_{m}^{2}(f)$ where
$M_{n}(f)=n^{-1}\sum_{i=1}^{n}\rho(Y_{i}-f(\mathbf{x}_{i}))$. Observe that
$L_{n}(f)$ is the sum of two convex functions and as such it is itself convex.
By Proposition 1 there exists a minimizer of $L_{n}(f)$ in
$\mathcal{H}^{m}(\mathbb{R}^{d})$, which we denote with $\widehat{f}_{n}$. Put

$\displaystyle\widetilde{f}_{n}=\alpha\widehat{f}_{n}+(1-\alpha)f_{0},$

for some $\alpha\in(0,1)$ to be chosen. As
$f_{0}\in\mathcal{H}^{m}(\mathbb{R}^{d})$ we have

$\displaystyle L_{n}(\widetilde{f}_{n})\leq\alpha
L_{n}(\widehat{f}_{n})+(1-\alpha)L_{n}(f_{0})\leq L_{n}(f_{0}),$

from where, after adding
$\mathbb{E}\\{M_{n}(\widetilde{f}_{n})-M_{n}(f_{0})\\}$ on both sides, we get

$\displaystyle\mathbb{E}\\{M_{n}(\widetilde{f}_{n})-M_{n}(f_{0})\\}+I_{m}^{2}(\widetilde{f}_{n})\leq\left[M_{n}(f_{0})-\mathbb{E}\\{M_{n}(f_{0})\\}-M_{n}(\widetilde{f}_{n})+\mathbb{E}\\{M_{n}(\widetilde{f}_{n})\\}\right]+I_{m}^{2}(f_{0}).$
(8)

We choose $\alpha=1/(1+\|\widehat{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})})$.
Clearly, $\alpha\in(0,1)$ and

$\displaystyle\|\widetilde{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}=\alpha\|\widehat{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}\leq
1.$

Our proof consists of deriving a lower bound on the left-hand side of (8) and
an upper bound on the right-hand side of (8), both in terms of
$\|\widetilde{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}$.

We begin with the lower bound. By assumption $2m>d$ and $\mathcal{O}$ is a
bounded set satisfying the uniform cone condition. Therefore, by Sobolev’s
theorem (Adams & Fournier, 2003, Theorem 4.12) we have the (compact) embedding

$\displaystyle\mathcal{H}^{m}(\mathcal{O})\to\mathcal{C}(\mathcal{O}),$

which implies the existence of a universal $c_{0}>0$ such that for any
$f\in\mathcal{H}^{m}(\mathcal{O})\subset\mathcal{C}(\mathcal{O})$,

$\displaystyle\sup_{\mathbf{x}\in\mathcal{O}}|f(\mathbf{x})|\leq
c_{0}\left\\{\int_{\mathcal{O}}|f(\mathbf{x})|^{2}d\mathbf{x}+I_{m}^{2}(f)\right\\}^{1/2}.$

Approximating $\int_{\mathcal{O}}|f(\mathbf{x})|^{2}d\mathbf{x}$ with
$\int_{\mathcal{O}}|f(\mathbf{x})|^{2}dQ_{n}(\mathbf{x})$, as in the proof of
Lemma 1, yields

$\displaystyle\sup_{\mathbf{x}\in\mathcal{O}}|f(\mathbf{x})|\leq
c_{0}\left\\{\int_{\mathcal{O}}|f(\mathbf{x})|^{2}dQ_{n}(\mathbf{x})+I_{m}^{2}(f)\right\\}^{1/2}.$

Now, since $I_{m}^{2}$ is a squared semi-norm and $I_{m}^{2}(f_{0})$ is
bounded, $I_{m}^{2}(f_{0})\leq 1$, say, this inequality reveals that for all
$f\in\mathcal{H}^{m}(\mathbbm{R}^{d})$ satisfying
$\|f-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}\leq 1$ we have

$\displaystyle\sup_{\mathbf{x}\in\mathcal{O}}|f(\mathbf{x})-f_{0}(\mathbf{x})|\leq
c_{0}\left\\{1+I_{m}^{2}(f-f_{0})\right\\}^{1/2}\leq
c_{0}\\{1+I_{m}^{2}(f)\\}^{1/2}.$

It follows that we can choose a large enough $D_{\kappa}>1$, not depending on
$f$, such that

$\displaystyle\sup_{\mathbf{x}\in\mathcal{O}}\frac{|f(\mathbf{x})-f_{0}(\mathbf{x})|}{D_{\kappa}\\{1+I_{m}^{2}(f)\\}^{1/2}}\leq\kappa,$

where $\kappa$ is the constant in assumption A2. Therefore, for all
$f\in\mathcal{H}^{m}(\mathbbm{R}^{d})$ satisfying
$\|f-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}\leq 1$, by A2, we find

$\displaystyle\mathbb{E}\\{M_{n}(f)-M_{n}(f_{0})\\}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left\\{\rho\left(\epsilon_{i}+f_{0}(\mathbf{x}_{i})-f(\mathbf{x}_{i})\right)-\rho\left(\epsilon_{i}\right)\right\\}$
$\displaystyle\geq\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left\\{\rho\left(\epsilon_{i}+\frac{f_{0}(\mathbf{x}_{i})-f(\mathbf{x}_{i})}{D_{k}\left\\{1+I_{m}^{2}(f)\right\\}^{1/2}}\right)-\rho\left(\epsilon_{i}\right)\right\\}$
$\displaystyle\geq\kappa\frac{\left\|f-f_{0}\right\|^{2}_{\mathcal{L}^{2}(Q_{n})}}{D_{k}^{2}\left\\{1+I_{m}^{2}(f)\right\\}}$

It follows that

$\displaystyle\inf_{f\in\mathcal{H}^{m}(\mathbbm{R}^{d}):\|f-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}\leq
1}\left[\frac{\mathbb{E}\\{M_{n}(f)-M_{n}(f_{0})\\}}{\kappa\frac{\|f-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}}{D_{\kappa}^{2}\\{1+I_{m}^{2}(f)\\}}}\right]\geq
1$

and from this, since, by construction,
$\widetilde{f}_{n}\in\mathcal{H}^{m}(\mathbbm{R}^{d})$ and
$\|\widetilde{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}\leq 1$, we see that

$\displaystyle\mathbb{E}\\{M_{n}(\widetilde{f}_{n})-M_{n}(f_{0})\\}$
$\displaystyle=\frac{\mathbb{E}\\{M_{n}(\widetilde{f}_{n})-M_{n}(f_{0})\\}}{\kappa\frac{\|\widetilde{f}_{n}-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}}{D_{\kappa}^{2}\\{1+I_{m}^{2}(\widetilde{f}_{n})\\}}}\kappa\frac{\|\widetilde{f}_{n}-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}}{D_{\kappa}^{2}\\{1+I_{m}^{2}(\widetilde{f}_{n})\\}}$
$\displaystyle\geq\kappa\frac{\|\widetilde{f}_{n}-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}}{D_{\kappa}^{2}\\{1+I_{m}^{2}(\widetilde{f}_{n})\\}}\inf_{f\in\mathcal{H}^{m}(\mathbbm{R}^{d}):\|f-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}\leq
1}\left[\frac{\mathbb{E}\\{M_{n}(f)-M_{n}(f_{0})\\}}{\kappa\frac{\|f-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}}{D_{\kappa}^{2}\\{1+I_{m}^{2}(f)\\}}}\right]$
$\displaystyle\geq\kappa\frac{\|\widetilde{f}_{n}-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}}{D_{\kappa}^{2}\\{1+I_{m}^{2}(\widetilde{f}_{n})\\}}$
(9)

This provides a lower bound for the left-hand side of (8) and completes the
first part of our derivation.

Next, we derive an upper bound for the right-hand side of (8). To accomplish
this, we need to derive the modulus of continuity of the mean-centered process
$M_{n}(f)-\mathbb{E}\\{M_{n}(f)\\}$. We will apply Lemma 8.5 of van de Geer
(2000) to this process. First, observe that by A1 for all
$(f,g)\in\mathcal{H}^{m}(\mathbbm{R}^{d})\times\mathcal{H}^{m}(\mathbbm{R}^{d})$
we have

$\displaystyle|\rho(Y_{i}-f(\mathbf{x}_{i}))-\rho(Y_{i}-g(\mathbf{x}_{i}))|\leq
c_{0}|f(\mathbf{x}_{i})-g(\mathbf{x}_{i})|,$

so that the lemma is applicable with
$d_{i}(f,g)=|f(\mathbf{x}_{i})-g(\mathbf{x}_{i})|$ in the notation of van de
Geer (2000). In combination with Lemma 1 we thus have

$\displaystyle\sup_{f\in\mathcal{H}^{m}(\mathbbm{R}^{d}):\|f-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}\leq
1}\left|\frac{M_{n}(f_{0})-M_{n}(f)-\mathbb{E}\\{M_{n}(f_{0})-M_{n}(f)\\}}{n^{-1/2}\|f-f_{0}\|^{1-d/2m}_{\mathcal{L}^{2}(Q_{n})}\\{1+I_{m}(f)\\}^{d/2m}}\right|=O_{P}(1),$

so that

$\displaystyle
M_{n}(f_{0})-M_{n}(\widetilde{f}_{n})-\mathbb{E}\\{M_{n}(f_{0})-M_{n}(\widetilde{f}_{n})\\}=O_{P}(n^{-1/2})\|\widetilde{f}_{n}-f_{0}\|^{1-d/2m}_{\mathcal{L}^{2}(Q_{n})}\\{1+I_{m}(\widetilde{f}_{n})\\}^{d/2m},$
(10)

which provides the desired upper bound for the right-hand side of (8).

Plugging the lower bound in (Proof of Theorem 1.) and the upper bound in (10)
into (8), we finally obtain

$\displaystyle
c_{0}\frac{\|\widetilde{f}_{n}-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}}{1+I_{m}^{2}(\widetilde{f}_{n})}+\lambda
I_{m}^{2}(\widetilde{f}_{n})\leq
O_{P}(n^{-1/2})\|\widetilde{f}_{n}-f_{0}\|^{1-d/2m}_{\mathcal{L}^{2}(Q_{n})}\\{1+I_{m}(\widetilde{f}_{n})\\}^{d/2m}+\lambda
I_{m}^{2}(f_{0}).$ (11)

This inequality implies that

$\displaystyle\|\widetilde{f}_{n}-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}=O_{P}(n^{-2m/(2m+d)})\quad\text{and}\quad
I_{m}(\widetilde{f}_{n})=O_{P}(1).$ (12)

To see this implication, note that if for real numbers $a,b,c$ we have $a\leq
b+c$ then either $a\leq 2b$ or $a\leq 2c$. Applying this on (11) leads to
either

$\displaystyle
c_{0}\frac{\|\widetilde{f}_{n}-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}}{1+I_{m}^{2}(\widetilde{f}_{n})}+\lambda
I_{m}^{2}(\widetilde{f}_{n})\leq
O_{P}(n^{-1/2})\|\widetilde{f}_{n}-f_{0}\|^{1-d/2m}_{\mathcal{L}^{2}(Q_{n})}\\{1+I_{m}(\widetilde{f}_{n})\\}^{d/2m},$
(13)

or

$\displaystyle
c_{0}\frac{\|\widetilde{f}_{n}-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}}{1+I_{m}^{2}(\widetilde{f}_{n})}+\lambda
I_{m}^{2}(\widetilde{f}_{n})\leq 2\lambda I_{m}^{2}(f_{0}).$ (14)

If (14) holds, (12) is easily verified. On the other hand, if (13) holds,
solving it we get

$\displaystyle\|\widetilde{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}=O_{P}(n^{-m/(2m+d)})\\{1+I_{m}(\widetilde{f}_{n})\\}^{\frac{d}{2m+d}}\\{1+I_{m}^{2}(\widetilde{f}_{n})\\}^{\frac{2m}{2m+d}},$

as well as

$\displaystyle\frac{I_{m}^{2}(\widetilde{f}_{n})}{\left\\{1+I_{m}(\widetilde{f}_{n})\right\\}^{\frac{2d}{2m+d}}\left\\{1+I_{m}^{2}(\widetilde{f}_{n})\right\\}^{\frac{2m-d}{2m+d}}}=O_{P}(n^{-2m/(2m+d)})\lambda^{-1}.$

By our assumptions on $\lambda$, the latter implies that
$I_{m}(\widetilde{f}_{n})=O_{P}(1)$. Hence,
$\|\widetilde{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}=O_{P}(n^{-m/(2m+d)})$
verifying (12) again.

The last step in our proof involves passage from $\widetilde{f}_{n}$ to
$\widehat{f}_{n}$. For this, first note that by definition of the convex
combination $\widetilde{f}_{n}$ and (12),

$\displaystyle\|\widetilde{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}=\frac{\|\widehat{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}}{1+\|\widehat{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}}=O_{P}(n^{-m/(2m+d)}),$

whence also
$\|\widehat{f}_{n}-f_{0}\|_{\mathcal{L}^{2}(Q_{n})}=O_{P}(n^{-m/(2m+d)})$.
Finally, since again by (12) $I_{m}(\widetilde{f}_{n})=O_{P}(1)$, by the
triangle inequality we get

$\displaystyle I_{m}(\alpha(\widehat{f}_{n}-f_{0}))\leq
I_{m}(\widetilde{f}_{n})+I_{m}(f_{0})=O_{P}(1).$

But then also $I_{m}(\widehat{f}_{n}-f_{0})=O_{P}(1)$, which implies the
result. The proof is complete.

∎

###### Proof of Corollary 1.

The first part of the Corollary follows from Theorem 3.4 of Utreras (1988)
which in our case reads

$\displaystyle\int_{\mathcal{O}}|\widehat{f}_{n}(\mathbf{x})-f_{0}(\mathbf{x})|^{2}d\mathbf{x}$
$\displaystyle\leq
c_{0}\int_{\mathcal{O}}|\widehat{f}_{n}(\mathbf{x})-f_{0}(\mathbf{x})|^{2}dQ_{n}(\mathbf{x})$
$\displaystyle\quad+c_{0}h_{\max,n}^{2m}\sum_{m_{1}+\ldots+m_{d}=m}\binom{m}{m_{1},\ldots,m_{d}}\int_{\mathcal{O}}\left|\frac{\partial^{m}(\widehat{f}_{n}(\mathbf{x})-f_{0}(\mathbf{x}))}{\partial
x_{1}^{m_{1}}\ldots\partial x_{d}^{m_{d}}}\right|^{2}d\mathbf{x},$

for some constant $c_{0}$ that does not depend on either $\widehat{f}_{n}$ or
$f_{0}$. Applying now Theorem 1 and our assumption
$h_{\max,n}=O(n^{-1/(2m+d)})$ yields

$\displaystyle\int_{\mathcal{O}}|\widehat{f}_{n}(\mathbf{x})-f_{0}(\mathbf{x})|^{2}d\mathbf{x}$
$\displaystyle=O_{P}(n^{-2m/(2m+d)})+c_{0}h_{\max,n}^{2m}O_{P}(1)$
$\displaystyle=O_{P}(n^{-2m/(2m+d)})+O_{P}(n^{-2m/(2m+d)})$
$\displaystyle=O_{P}(n^{-2m/(2m+d)}),$

as asserted.

For the second part of the Corollary we apply an interpolation inequality due
to Nirenberg (1959), a simplified version of which states

$\displaystyle\int_{\mathcal{O}}\left|\frac{\partial
f^{j}(\mathbf{x})}{\partial x_{1}^{j_{1}}\ldots\partial
x_{d}^{j_{d}}}\right|^{2}d\mathbf{x}$ $\displaystyle\leq
c_{0}\left\\{\int_{\mathcal{O}}\left|\frac{\partial
f^{m}(\mathbf{x})}{\partial x_{1}^{m_{1}}\ldots\partial
x_{d}^{m_{d}}}\right|^{2}d\mathbf{x}\right\\}^{j/m}\left\\{\int_{\mathcal{O}}|f(\mathbf{x})|^{2}d\mathbf{x}\right\\}^{1-j/m}$
$\displaystyle\quad+c_{0}\int_{\mathcal{O}}|f(\mathbf{x})|^{2}d\mathbf{x},$

for every $f\in\mathcal{H}^{m}(\mathcal{O})$, where $c_{0}$ is a universal
constant and the inequality holds for all tuples $(j_{1},\ldots,j_{d})$ and
$(m_{1},\ldots,m_{d})$ such that $j_{1}+\ldots+j_{d}=j$ and
$m_{1}+\ldots+m_{d}=m$, respectively. Now apply this inequality with $f$
replaced by $\widehat{f}_{n}-f_{0}$ and use the first part of the corollary
and Theorem 1 to get

$\displaystyle\int_{\mathcal{O}}\left|\frac{\partial\widehat{f}_{n}^{j}(\mathbf{x})}{\partial
x_{1}^{j_{1}}\ldots\partial x_{d}^{j_{d}}}-\frac{\partial
f^{j}_{0}(\mathbf{x})}{\partial x_{1}^{j_{1}}\ldots\partial
x_{d}^{j_{d}}}\right|^{2}d\mathbf{x}$
$\displaystyle=O_{P}(1)O_{P}(n^{-2(m-j)/(2m+d)})+O_{P}(n^{-2m/(2m+d)})$
$\displaystyle=O_{P}(n^{-2(m-j)/(2m+d)}),$

which completes the proof. ∎

###### Proof of Theorem 2.

For uniformly sub-Gaussian errors $\epsilon_{i}$ and under Lemma 1, the
derivation on van de Geer (2000, p. 168) shows that

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}(\widehat{f}_{n}(\mathbf{x}_{i})-f_{0}(\mathbf{x}_{i}))=O_{P}(n^{-1/2})\|\widehat{f}_{n}-f_{0}\|^{1-d/2m}_{\mathcal{L}^{2}(Q_{n})}\\{1+I_{m}(\widehat{f}_{n})\\}^{d/2m}.$

Plug this into (4) to obtain

$\displaystyle\|\widehat{f}_{n}-f_{0}\|^{2}_{\mathcal{L}^{2}(Q_{n})}+\lambda
I_{m}^{2}(\widehat{f}_{n})\leq
O_{P}(n^{-1/2})\|\widehat{f}_{n}-f_{0}\|^{1-d/2m}_{\mathcal{L}^{2}(Q_{n})}\\{1+I_{m}(\widehat{f}_{n})\\}^{d/2m}+\lambda
I_{m}^{2}(f_{0}).$

Now argue as in the proof of Theorem 1 to complete the proof.

∎

## Acknowledgements

The author is grateful to two anonymous referees, the associate editor and the
editor for suggestions that lead to a much improved manuscript. He is also
grateful to Stefan Van Aelst for helpful discussions. This research has been
supported by the Research Foundation-Flanders (project 1221122N).

## References

  * Adams & Fournier (2003) authorAdams, R.A. & authorFournier, J.J.F. (year2003) titleSobolev Spaces, 2nd ed., publisherElsevier/Academic Press, Amsterdam.
  * Cox (1984) authorCox, D.D. (year1984). titleMultivariate Smoothing Spline Functions. journalSIAM J. Numer. Anal. volume21, pages789–813.
  * Cox & O’Sullivan (1985) authorCox, D.D. & authorO’Sullivan, F. (year1985). titleAnalysis of penalized likelihood type estimators with applications to generalized smoothing in Sobolev Spaces. journalTech. Rep. No. 51, Dept. of Statistics, University of California, Berkeley.
  * Cucker & Smale (2001) authorCucker, F. & authorSmale, S. (year2001). titleOn the Mathematical Foundations of Learning. journalBul. Amer. Math. Soc. volume39 pages1–49.
  * Green & Silverman (1994) authorGreen, P.J. & authorSilverman, B.W. (year1994). titleNonparametric Regression and Generalized Linear Models. publisherChapman & Hall, London.
  * Györfi et al. (2010) authorGyörfi, L., authorKohler, M., authorKrzyżak, A. & authorWalk, H. (year2010). titleA Distribution-Free Theory of Nonparametric Regression. publisherSpringer, New York.
  * Hsing & Eubank (2015) authorHsing, T. & authorEubank, R. (year2015). titleTheoretical foundations of functional data analysis, with an introduction to linear operators. publisherWiley, New York.
  * Kalogridis (2021) authorKalogridis, I. (year2021). titleAsymptotics for M-type smoothing splines with non-smooth objective functions. journalTEST volume31, pages373–389.
  * Kalogridis and Van Aelst (2022) authorKalogridis, I. & authorVan Aelst, S. (year2022). titleRobust optimal estimation of location from discretely sampled functional data. journalScand. J. Stat., volumeappeared online.
  * Kalogridis et al. (2023) authorKalogridis, I., authorClaeskens, G. & authorVan Aelst, S. (year2023). titleRobust and efficient estimation of nonparametric generalized linear models. journalTEST, appeared online.
  * Maronna (2011) authorMaronna, R.A. (year2011). titleRobust ridge regression for high-dimensional data. journalTechnometrics volume53, pages44–53.
  * Nirenberg (1959) authorNirenberg, L. (year1959). titleOn elliptic partial differential equations. journalAnnali della Scuola Normale Superiore di Pisa - Classe di Scienze volumeSerie 3, Tome 13, pages115–162.
  * Nychka et al. (2017) authorNychka, D., authorFurrer, R., authorPaige, J. & authorSain, S. (year2017). titlefields: Tools for spatial data. journalR-package version 14.0. https://www.r-project.org/.
  * R Core Team (2022) authorR Core Team (year2022). titleR: A language and environment for statistical computing. organizationR Foundation for Statistical Computing, addressVienna, Austria. https://www.R-project.org/.
  * Rynne & Youngston (2008) authorRynne, B. & authorYoungston, M.A. (year2008). titleLinear functional analysis. publisherSpringer, London.
  * Stone (1982) authorStone, C.J. (year1982). titleOptimal Global Rates of Convergence for Nonparametric Regression. journalAnn. Statist. volume10, pages1040–1053.
  * Utreras (1988) authorUtreras, F. (year1988). titleConvergence Rates for Multivariate Smoothing Spline Functions. journalJ. Approx. Theory volume52, pages1–27.
  * van de Geer (2000) authorvan de Geer, S. (year2000). titleEmpirical processes in M-estimation. publisherCambridge University Press, New York, New York.
  * van de Geer (2002) authorvan de Geer, S. (year2002). titleM-estimation using penalties or sieves. journalJ. Statist. Plann. Inference volume108, pages55–69.
  * van der Vaart & Wellner (1996) authorvan der Vaart, A. W. & authorWellner, J. A. (year1996). titleWeak convergence and empirical processes. publisherSpringer, New York.
  * Wahba (1990) authorWahba, G. (year1990). titleSpline models for observational data. publisherSiam, Philadelphia, Pen.
  * Wood (2017) authorWood, S.N. (year2017). titleGeneralized Additive Models, 2nd ed.. publisherCRC Press, Boca Raton, FL.
  * Yohai and Zamar (1988) authorYohai, V.J. & authorZamar, R.H. (year1988) titleHigh Breakdown-Point Estimates of Regression by Means of the Minimization of an Efficient Scale. journalJ. Amer. Statist. Assoc. volume83 pages406–413.


# The galactic tooth-fairy and a cosmic bullet: Amateur discoveries and a call
for further research

Muhammad Shaheer Niazi,1
1Institute of space technology, Islamabad, Pakistan
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

There are countless digital sky surveys and automated scans of the night sky
which use computer algorithms to detect and categorize objects. With the
advent of Artificial Intelligence such surveys will become even more efficient
in the near future. Despite this some objects are missed by surveys or pose no
initial interest. At times such missed objects are unique in nature and of
decent angular sizes, demanding research, unlike the billions of tiny specs of
galaxies that would be too tedious to name and study. In this scenario the
amateur astronomer and their spirit for old school astronomical discovery
steps in, to manually comb the sky and catalogue unique objects as was done in
the early days of astronomy. In this paper two unique, previously uncatalogued
galaxy candidates, namely Shaheer I and Shaheer II are identified and studied.
Both galaxies lay at a distance of 6.67 arc-minutes from each other in the
constellation of Camelopardalis. One boasts an unusual morphological profile,
akin to a molar tooth, while the other seems to be shooting through space at
tremendous velocities. The objects were discovered during visual inspection of
digital surveys and then imaged from amateur telescopes at Taqwa observatory,
Pakistan’s first and only dark sky observatory (bortle 1). We perform
photometry using PetroFit to discuss the potential nature of the galaxies and
implore further collaborative research to fully uncover their characteristics.

###### keywords:

methods: observational –catalogues – Galaxy: general

††pubyear: 2023††pagerange: The galactic tooth-fairy and a cosmic bullet:
Amateur discoveries and a call for further research –A

##  1 Introduction

Amateur astronomy and the digitization of sky surveys (e.g. pan-STARRS,
Chambers et al., 2016) has brought in a surge of new discoveries over the last
decade. In the past, the astronomer would manually survey photographic plates
and come upon a new discovery such as with Bothun et al. (1987). Digitized
surveys now make it possible for anyone to explore the cosmos and engage in
citizen science. Today the dedicated amateur astronomer can build a decent
observatory with large aperture SCT’s (schmidt-Cassegrain telescope), a dark
sky and a high resolution, high Quantum efficiency imaging sensor enabling
them to contribute to science with their observations. The rapid development
and accessibility of CMOS imaging technologies has a significant role to play
here. Examples of contributions include DGSAT: Dwarf Galaxy Survey with
Amateur Telescopes (Javanmardi et al., 2016) . Their team has made many low
surface brightness dwarf galaxy discoveries (Martínez-Delgado et al., 2021;
Collins et al., 2022) in the local group using preliminary amateur discoveries
and then larger professional setups to probe deeper into the galaxies.
Similarly KPS-1b was the first transiting exoplanet discovered through an
amateur setup, Burdanov et al. (2018). Due to the trend of deep space astro-
photography and the associated long integration times (sometimes reaching over
100 hours), we have begun to uncover hidden objects which were previously not
visible due to the magnitude limitations of sky surveys. Planetary Nebulae are
one such class of objects being discovered primarily by amateur astronomers,
e.g. Dû et al. (2022). Being almost invisible in sky surveys, they are brought
out by deep imaging of sectors with potential candidates using narrow-band
filters and advanced processing techniques. The aim of this paper, alongside
describing two previously uncatalogued galaxies, is to promote scientific
discovery and research for the amateur community and to introduce the latest
accessible computational methods which can be easily used for preliminary
analysis. The related intermediate processing steps are also detailed to make
it easy to understand.

##  2 Discovery

The galaxy candidate Shaheer I was found during visual inspection of the
digital pan-STARRS survey near the MCG+12-06-003 group of galaxies in the
camelopardalis constellation. Due to its unusual shape and dim profile it
caught interest leading to its coordinates being noted down. On searching
through various catalogues and online databases, it was established that
automated surveys had not detected the object as a potential galaxy. After the
field of view around Shaheer I was imaged at Taqwa observatory (see3.1) and
inspected, Shaheer II was found just 6.67 arc minutes to the south. Despite
its brighter surface profile and obvious galactic nucleus, it had also not
been detected in galactic surveys however some point source astrometric and
photometric data (see table 2) is available through the extended Gaia-PS1-SDSS
(GPS1+) proper motion catalog (Tian et al., 2020) for Shaheer II and the pan-
STARRS release (PS1) survey for Shaheer I. See table 1 for coordinates.

Table 1: Coordinates of both galaxy candidates. | Shaheer I | Shaheer II
---|---|---
RA(J2000) | 05 25 15.26 | 05 24 54.28
Dec(J2000) | +72 56 31.06 | +72 49 58.52

##  3 Imaging

Imaging was done at Taqwa Space Observatory (26°27’35.0"N 66°18’28.0"E) (see
appendix A) using two telescopes. Luminance data was taken using a Meade lx200
acf 16 inch SCT with an FLI ML29052 CCD operating at -20°C. RBG data was taken
through a Celestron edge HD 8 inch SCT and an ASI294mc CMOS colour camera
operating at -15°C. The data was stacked in APP (Astro Pixel Processor) and no
further editing was done.

###  3.1 Luminance

Luminance here will refer to filter-less, monochrome imaging (Fig. 1, 2) of
the objects to try and resolve surface features. We imaged a large field of
view containing Shaheer I and II for a total integration time for 37,080
second (10.3 hours) over multiple nights in February 2023. Single exposures of
180 second were used taking in consideration tracking errors. With this setup
we had an effective angular resolution of 0.28"/px.

Figure 1: cropped frame showcasing Shaheer I. The tooth like structure is
evident. 20 percent dynamic stretch. FOV 2.4’ x 2.4’  Figure 2: cropped frame
showcasing Shaheer II. Obvious nucleus and cluster regions on rim visible. 15
percent dynamic stretch as this object is brighter. FOV 2.4’ x 2.4’

###  3.2 RGB colour

Colour data (Fig. 3) was imaged for a total integration time of 29,760 seconds
(8.2 hours) with single exposures of 120 seconds. The CMOS camera used and the
smaller aperture telescope limited the resolution available for the colour
images in comparison to the high resolution luminance images.

Figure 3: Shaheer I (left), bluish colour is visible. Stretched to 20 percent
of the dynamic range as it is a dim object. Saturation applied to enhance
colour. Shaheer II (right), whitish with hints of faint blue. Stretched to 15
percent of the dynamic range.

##  4 Photometry

###  4.1 Source selection

Photometry was performed using the newly developed PetroFit (Geda et al.,
2022) which helps create a petrosian profile (Petrosian, 1976) and fitting
model. It was used to calculate various structural parameters of the galaxies,
namely; Ellipticity, elongation, Sersic index and Half-light radius
($r_{eff}$). PetroFit also calculated the noise levels of the luminance images
which serves to estimate detection thresholds and background subtraction. The
sources in each image were then identified and segmented (fig 4) using a
detection threshold of $2\sigma$ ( where $\sigma$ is the data standard
deviation.) and a minimum source size of 5 pixels. A gaussian Kernel size of 3
was chosen with FWHM values of 5 and 3 for Shaheer II and Shaheer I
respectively. This was due to the differing surface brightness of both
objects. Segment de-blending was used for Shaheer II to isolate the star near
its rim while no de-blending was necessary for Shaheer I. A point spread
function (PSF) of the sensor used was estimated by selecting one bright star
from the FOV and normalizing it, see fig 5.

Figure 4: FITS data and isolated sources for Shaheer II (left) and Shaheer I
(right). Note no deblending of sources in the right image. Figure 5:
Normalized star PSF cutout from the FOV in odd dimensions of 53x53 pixels.

###  4.2 Curve of Growth

The radial profiles (fig 7, 8) of both candidates in the x and y axis along
with the enclosed flux curve of growth (fig 6) was modelled using PetroFit’s
source photometry. Note that the center of the object is automatically
determined, in the case of Shaheer I, its asymmetric shape causes some degree
of error. The background is subtracted to ensure the light profile of the
galaxy is measured only.

Figure 6: Curve of growths and PetroFit’s apertures visualized on the targets.
Shaheer I bottom, Shaheer II on top. Note each pixel is 0.28 arcsec. Shaheer I
is stretched to a greater degree to make it more apparent. Figure 7: Radial
profile along x and y axis for Shaheer I. Figure 8: Radial profile along x and
y axis for Shaheer II.

###  4.3 Petrosian profile

Using the fluxes found through the appertures, the petrosian profile (fig 9)
was created from which the petrosian radius, effective half light radius,
total flux radius and concentration indexes for both were found. These values
are listed in table 2. Values for $\epsilon$ were changed from the standard
value of 2 used by the Sloan Digital Sky Survey (SDDS, Strauss et al. (2002)).
This was done to keep the total flux radius within where the petrosian values
lie $1\geq\eta\geq 0$. For $\epsilon=2$, the total flux radius would lie where
$\eta<0$. PetroFit plots the curves of growth (in green) along with the
petrosian profile and marks other important values as well. Fig 10 visualizes
the various radii on the images of each galaxy.

Figure 9: Petroskian profiles for Shaheer I (top, $\epsilon=1.1$) and Shaheer
II (bottom, $\epsilon=1.2$). Green dashed curves are the curve of growth.
Figure 10: Visualized radii for both galaxies.

###  4.4 Sérsic model

Finally Sérsic models were created (fig 11), which is done using the PSF
(shown earlier in fig 6). No oversampling or weights were used. The sérsic
index, $n$, was estimated using the fitted models, see table 2. Residuals were
produced by subtracting the model from the data.

Figure 11: Sérsic fitted models and residuals. Table 2: Important values, note no error map was available to estimate uncertainties. Radial motion and r-Mag data for Shaheer II was obtained from GPS1+. r-Mag data for Shaheer 1 from pan-STARRS DR1 (PS1).  | Shaheer I | Shaheer II
---|---|---
$\epsilon$ | 1.1 | 1.2
$r_{petro}(arcmin)$ | 15.97 | 14.25
$r_{total}(arcmin)$ | 17.57 | 17.08
$r_{eff}(arcmin)$ | 7.91 | 6.89
$a/b$ | 0.382 | 0.177
$elong$ | 1.619 | 1.216
$n$ | 0.857 | 0.925
$r_{90}(arcmin)$ | 12.24 | 14.19
$c_{5090}$ | 1.25 | 1.27
$\mu_{\xi}(RA,mas)$ | - | -40.739
$\mu_{\eta}(Dec,mas)$ | - | -27.696
$M_{r}(mag)$ | 20.538 | 19.105

##  5 Discussion

Shaheer I is definitely the more elusive of the galaxy (or dwarf galaxy)
candidates. Its shape is asymmetric, as is evident in the radial profile given
in fig 7. On first look it resembles a molar tooth (thus the title galactic
tooth-fairy). With no apparent nucleus its low surface brightness could be the
reason it has not been investigated prior. Its blue color profile would
indicate that it is mostly comprised of young hot stars. Among the two,
Shaheer II is one which contains an obvious galactic nucleus and possesses a
tear drop like shape. There are also star cluster regions towards what can be
considered the top rim of the galaxy. These regions are more evident in the
residual after the sérsic model is subtracted in fig 11. Young stars might
populate the regions owing to the blue-white colour profile. The presence of
these brighter regions are only visible through current images in the top rim
of the galaxy, while the southern region is a diffuse halo. Despite its
apparent symmetry, the diffuse halo might extend a bit further outwards. The
large proper motion of the source found in GPS1+ (see table 2) suggests that
it is a high propermotion/velocity galaxy (thus the title cosmic bullet)
creating a cometary phenomena similar to the "Comet galaxy" introduced in
Cortese et al. (2007). In terms of gravitational influences, there are two
potential attractor’s which might be pulling Shaheer II creating its teardrop
shape. The first is a galaxy, 2MASX J05251648+7253484, (2MASS; Skrutskie et
al., 2006), which is closest in terms of angular distances (see fig 12). The
direction of the top rim coincides with the direction of 2MASX
J05251648+7253484. The second attractor candidate would be MCG+12-06-003 (fig
13), which is part of a group of galaxies 11.18 arc-minutes away, which also
appear to be interacting with each other. Directions to both galaxies are
shown in fig 14. Spectral analysis will be needed to calculate further
parameter. Deeper, higher resolution imaging through a larger ground or space
based telescope would be required to resolve the regions and complete
structure of each galaxy.

Figure 12: The galaxy 2MASX J05251648+7253484 and Shaheer II, line drawn
between them, the angular distance between them is 4.164 arcsec. Image cropped
from the total FOV taken from Taqwa observatory. Figure 13: MCG+12-06-003 and
associated galaxy group. MCG+12-06-003 seems to be interacting with the galaxy
to its right with visible tidal streams, albeit faint. The black arrow points
towards the direction of Shaheer II. Image cropped from the total FOV taken
from Taqwa observatory. Figure 14: Directions from Shaheer II towards
MCG+12-06-003 (green) and 2MASX J05251648+7253484 (red).

##  6 Conclusions

This paper presents the discovery and analysis of two new galaxies, Shaheer I
and Shaheer II with images taken through amateur telescopes at the Taqwa space
observatory. The galaxies were found through visual inspection of the pan-
STARRS survey. The images which have resolutions and limiting magnitudes
comparable to many optical surveys of the night sky, help prove that amateur
contributions towards astronomical research can go a long way. Using the data
taken through our observatory, photometric analysis was performed using
PetroFit (while detailing each step for reproducability). I extracted as much
quantitative data along with qualitative data, through visual and photometric
analysis, as was possible. Deeper imaging of both galaxies will be necessary
to estimate distances and to probe further into the nature of both galaxies,
especially for Shaheer I. I propose further collaborative research
incorporating the use of larger telescopes and photometric filters to obtain
more precise data. I also invite the amateur community to take deeper, long
integration, images of the galaxies and perform any kinds of analyses
possible. Data to be found in future studies would include detailed
morphological classification and spectra along with investigating and proving
the high proper motion of Shaheer II.

## Acknowledgements

All images and data was capture at Taqwa space observatory in Bela,
Balochistan, Pakistan, see Appendix A. Thanks to the Taqwa team for their
great setup and for the imaging time.

This research made use of Photutils, an Astropy package for detection and
photometry of astronomical sources (Bradley et al., 2022). This research made
use of PetroFit (Geda et al., 2022), a package based on Photutils, for
calculating Petrosian properties and fitting galaxy light profiles.

## Data Availability

PetroFit can be found at github.com/PetroFit/petrofit/tree/main. Extracted
astrometry and photometry (Magnitudes and proper motions) are mentioned in the
article and can be found from VizieR.

PS1 (pan-STARRS data release 1), obj id:195530813054621280 and GSP1+ obj id:
-9050757507386836921.

Astro pixel processor (APP) can be purchased from
www.astropixelprocessor.com/.

## References

  * Bothun et al. (1987) Bothun G. D., Impey C. D., Malin D. F., Mould J. R., 1987, AJ, 94, 23
  * Bradley et al. (2022) Bradley L., et al., 2022, astropy/photutils:1.5.0, doi:10.5281/zenodo.6825092, https://doi.org/10.5281/zenodo.6825092
  * Burdanov et al. (2018) Burdanov A., et al., 2018, PASP, 130
  * Chambers et al. (2016) Chambers K. C., et al., 2016, arXiv e-prints, p. arXiv:1612.05560
  * Collins et al. (2022) Collins M. L., Charles E. J., Martínez-Delgado D., Monelli M., Karim N., Donatiello G., Tollerud E. J., Boschin W., 2022, MNRAS: Letters, 515, L72
  * Cortese et al. (2007) Cortese L., et al., 2007, MNRAS, 376, 157
  * Dû et al. (2022) Dû P. L., et al., 2022, A&A, 666, A152
  * Geda et al. (2022) Geda R., Crawford S. M., Hunt L., Bershady M., Tollerud E., Randriamampandry S., 2022, AJ, 163, 202
  * Javanmardi et al. (2016) Javanmardi B., et al., 2016, A&A, 588
  * Martínez-Delgado et al. (2021) Martínez-Delgado D., et al., 2021, A&A, 652
  * Petrosian (1976) Petrosian V., 1976, ApJ, 210, L53
  * Skrutskie et al. (2006) Skrutskie M. F., et al., 2006, AJ, 131, 1163
  * Strauss et al. (2002) Strauss M. A., et al., 2002, AJ, 124, 1810
  * Tian et al. (2020) Tian H.-J., Xu Y., Liu C., Rix H.-W., Sesar B., Goldman B., 2020, AJ Supplement Series, 248, 28

##  Appendix A Taqwa Space observatory

Taqwa Space observatory located under the bortle 1 skies of Bela, Balochistan,
is the first and only dark sky observatory in Pakistan. Pakistan’s
astronomical community is mostly comprised of amateurs with their private
telescopes. We have no research grade observatory in the country to produce
data but our dark skies have a lot of potential. Taqwa Space observatory is a
privately funded observatory that is gearing itself with the best equipment
available to utilize the dark and clean skies of Balochistan, Pakistan and
begin in house observational research through the acquisition of larger
telescopes in the future. See fig 15 and 16 for images. Webpage:
https://www.taqwaspaceobservatory.com/

Figure 15: Equipment at Taqwa Space observatory. Figure 16: Aerial view of the
observatory and surroundings.


# Topological Data Analysis of Neural Network Layer Representations

Archie Shahidullah
Computing + Mathematical Sciences
California Institute of Technology
Pasadena, CA 91125
<EMAIL_ADDRESS>
###### Abstract

This paper is a cursory study on how topological features are preserved within
the internal representations of neural network layers. Using techniques from
topological data analysis, namely persistent homology, the topological
features of a simple feedforward neural network’s layer representations of a
modified torus with a Klein bottle-like twist were computed. The network
appeared to approximate homeomorphisms in early layers, before significantly
changing the topology of the data in deeper layers. The resulting noise
hampered the ability of persistent homology to compute these features, however
similar topological features seemed to persist longer in a network with a
bijective activation function.

##  1 Introduction

In recent years, deep neural networks have revolutionized many computing
problems once thought to be difficult. However, they are often referred to as
black boxes and the mechanisms through which they learn have remained elusive.
Insight into how a neural network internally represents the dataset it is
trained on will provide understanding into what makes for effective training
data, and how a neural network extracts relevant features from a dataset.

###  1.1 Layer Representations

A feedforward neural network, $N$, that acts on input data $x$ can be viewed
as a composition of $n$ functions, where each function represents a layer,
$L_{i}$,

$N(x)=L_{n}(L_{n-1}(\dots L_{1}(x)))$

Each layer, acting on input $z$, consists of a weights matrix $W$, bias vector
$b$, and (usually) nonlinear activation function $\sigma$,

$L_{i}(z)=\sigma(Wz+b)$

The output of $L_{i}$ is what we will refer to as a layer representation. It
is important to note that more complicated architectures are not restricted to
a linear chain structure and thus can have skip and recurrent connections that
complicate their graph structure. However, skip connections simply take in
multiple inputs and recurrent (and recursive) connections can be unrolled and
layer representations can be recovered.

###  1.2 Overview of Paper

Section 2 gives relevant background about topology and its relation to neural
networks. Section 3 gives an overview of persistent homology, a technique to
compute the topological features of a point-cloud dataset. Section 4 outlines
the experiment performed and its results. Finally, Section 5 has the
discussion of results and directions for future work.

##  2 Background and Previous Work

We will define topology, give motivation for the homology groups of a
topological space, its applicability to neural networks, and review previous
work on the subject.

###  2.1 Definition of Topology

Topology is the study of geometric objects with a particular structure that
allows for a rigorous treatment of the concepts of "bending" and "twisting" a
space and how properties of the space are preserved under these
transformations. Formally, a topological space is the tuple $(X,\tau)$, where
$X$ is some set and $\tau$ is a multiset consisting of subsets of $X$.
Armstrong (1983) gives the following definition of a topological space,

######  Definition 2.1 (Topological space).

Given a set $X$ and a multiset $\tau$ that consists of subsets of $X$, a
topological space is the tuple $(X,\tau)$ that fulfills the following axioms,

  1. 1.

$\emptyset\in\tau$ and $X\in\tau$

  2. 2.

$\forall S_{i}\in\tau,\bigcup_{i}S_{i}\in\tau$ for finite or infinite unions

  3. 3.

$\forall S_{i}\in\tau,\bigcap_{i}S_{i}\in\tau$ for only finite intersections

Any member of $\tau$ is termed an open set (and its complement is a closed
set), and $\tau$ is the topology on $X$.

Another central concept in topology is continuous deformation between
topological spaces, and this is formalized in the idea of a homeomorphism,
defined as,

######  Definition 2.2 (Homeomorphism).

A homeomorphism $f:X\to Y$ is an isomorphism between two topological spaces
$X$ and $Y$ and therefore fulfills the following critera,

  1. 1.

$f$ is bijective

  2. 2.

$f$ is continuous

  3. 3.

$f^{-1}$ is continuous

###  2.2 Homology

The homology of a topological space informally characterizes the number of
"holes" in the space. If we take a cycle to be a generalization of a closed
loop on some space, such as on the surface of the sphere $S^{2}$, we can
classify cycles by whether or not they can be continuously deformed into each
other. If two cycles cannot be deformed into each other, it is said that there
exists a hole on the topological space.

A 0-dimensional hole is a connected component, a 1-dimensional hole is a loop,
a 2-dimensional hole is a shell, and so on. The study of these holes requires
a formal description of the boundaries on topological spaces. Boundaries in
general are linear combinations of more basic geometric objects, which
motivates us to introduce some sort of structure to allow for this. The
homology of a topological space $X$ is formally represented by its homology
groups,

$H_{0}(X),H_{1}(X),H_{2}(X),\dots$

Each homology group $H_{k}(X)$ has an abelian group structure, and as such we
naturally use $\mathbb{Z}$. We refer to the rank of a homology group as the
Betti number $b_{k}$. In general,

$H_{k}(X)=\mathbb{Z}\times\cdots\times\mathbb{Z}=Z_{b_{k}}$

A Betti number counts the number of holes in a topological space and is
topologically invariant, which makes it an ideal candidate to judge the
topological features of a space. Under homeomorphism, the Betti numbers of our
domain and codomain remain unchanged. It is important to note that the Betti
numbers do not account for all topological invariants, such as torsion.
Torsion refers to features such as the twist of a Möbius strip.

The Betti numbers of $S^{1}$, the circle, are $b_{0}=1,b_{1}=1$, and 0
otherwise. This is because there is one connected component, and one loop
(1-dimensional boundary). The Betti numbers of $S^{2}$ are $b_{0}=1,b_{2}=1$,
and 0 otherwise. This is because there is 1 connected component and one
2-dimensional boundary (around the interior of the sphere). Notably, there are
no loops that can be drawn on the surface that cannot be deformed to a single
point. The torus $T^{2}$, a doughnut-shaped object, has Betti numbers
$b_{0}=1,b_{1}=2,b_{2}=1$, and 0 otherwise. There is one connected component,
and two loops (one around the ring and another around the "main hole").
Lastly, there is a 2-dimensional boundary around the interior. It is
interesting to note $T^{2}=S^{1}\times S^{1}$.

###  2.3 Application to Neural Networks

We are interested in whether the layer representation $L_{i}$ acts similar to
a homeomorphism in that it allows a neural network to represent the
topological features of a dataset. Unfortunately, neural network layers are
not in general homeomorphisms,.

######  Theorem 2.1.

A neural network layer, $L(x)=\sigma(Wx+b)$, need not be a homeomorphism.

###### Proof.

If $W$ is not a member of the general linear group
$\mathrm{GL}_{n}(\mathbb{R})$, it is not invertible and therefore no
homeomorphism exists by 2.2. Additionally, if $\sigma$ is not bijective, no
homeomorphism can exist also by 2.2. ∎

However, the topological approach is still useful. If layer representations
are examined and are found to resemble the topological features of the
dataset, this suggests a reason neural networks are effective is that they
learn a robust representation of the topological features of a space.

###  2.4 Previous Work

Studying neural network learning with topology has been a popular approach in
theoretical machine learning. This is best exemplified with the concept of the
manifold hypothesis. Fefferman, et al. (2013) gives the following definition
of the manifold hypothesis,

######  Definition 2.3 (Manifold Hypothesis).

High-dimensional data tends to lie near low-dimensional manifolds.

This concept is best explained with an example. Imagine there exists a model
to classify $m\times n$ images between cats and dogs. The data space is
$\mathbb{R}^{mn}$. However, this space clearly contains data not relevant to
the task at hand (such as images of flowers or random noise), and the manifold
hypothesis conjectures that there exists a much lower-dimensional submanifold
of the data space that approximates the relevant data.

Given that manifolds are topological spaces, it is reasonable to assume
methods from topology can be used to gain insight into the learning process,
should the manifold hypothesis be correct. There exists a field called
topological data analysis (TDA) precisely focused on extracting topological
features from data. The most common technique used in TDA is persistent
homology. Persistent homology will be described in the next section, but it
essentially computes whether the features corresponding to each homology group
exist in a point-cloud, e.g. if one can draw loops via linear combinations of
datapoints. Montúfar, et al. (2020) showed that a neural network can be
trained to approximate topological features similar to persistent homology.
While this paper will use persistent homology to compute topological features,
it is encouraging that a neural network can be trained to compute the
topological features of a dataset.

##  3 Persistent Homology

Persistent homology is a technique to compute the topological features of a
point-cloud dataset and returns these results in a persistence diagram.

###  3.1 Method

Given a point-cloud, persistent homology grows $n$-dimensional balls around
each point and marks when these balls intersect with the balls of another
point. By connecting these balls together, one can represent the data as
simplicial complexes by chaining together simplices (generalizations of
points, lines, and triangles in higher-dimensional space). The validity of
this is a result of the abelian group structure imposed on the homology of our
dataset. These simplicial complexes can then represent topological features
such as loops. At some point, all points will intersect as the balls will grow
to encompass the whole dataset. We record the birth and death times of these
features and plot them in a persistence diagram.

###  3.2 Persistence Diagrams

We generated a variation of a torus with a Klein bottle-like twist, $X$,
according to the following parametrization,

$x=(R+2P\cos\theta)\cos\phi\;,$ $y=(R+2P\cos\theta)\sin\phi\;,$
$z=2P\sin\theta\cos\frac{\phi}{2}\;,$ $w=2P\sin\theta\sin\frac{\phi}{2}$

We display the point-cloud (only the first three dimensions) and its
persistence diagram in figure 1. For $H_{0}(X)$, we can see that all features
but one die soon after they are birthed, which is indicated by points residing
close to the diagonal of the diagram. The single point that persists at
infinity corresponds to the single connected component obtained when all the
balls intersect in persistent homology. For $H_{1}(X)$, we similarly see most
features die close to birth. However, there are two features that appear to
persist for a significant time before ultimately dying during their merger.
These correspond to loops found by chaining together simplices generated from
datapoints to form a simplicial complex by the algorithm. In general, the
Betti numbers are difficult to extract for complicated persistence diagrams,
however we can make a reasonable guess that were this manifold to be taken in
the continuum limit, we would get Betti numbers $b_{0}=1$ and $b_{1}=2$.

((a)) The three-dimensional projection of the point cloud.

((b)) The persistence diagram of the point-cloud.

Figure 1: A torus with a Klein bottle-like twist along with its persistence
diagram.

##  4 Experiment and Results

Table 1: Architecture of Neural Network Name | Hidden Units | Activation Function
---|---|---
Input | 4 | ReLU
Layer 1 | 10 | ReLU
Layer 2 | 30 | ReLU
Layer 3 | 10 | ReLU
Output | 1 | Sigmoid

This section outlines the experiment conducted, the layer representations
obtained, their persistence diagrams, and a PCA projection of them.

###  4.1 Experiment Design

The experiment consisted of the modified torus parametrized in the previous
section along with noise sampled from the uniform distribution. The total
dataset length was 9800. The network, whose architecture is given in table 1,
was trained on the binary classification task of distinguishing whether or not
a point lay on the modified torus. The model was built using TensorFlow.

Using the Adam optimizer, the network was trained for 300 epochs when an
accuracy of approximately 98% was obtained. Using the HDBSCAN clustering
algorithm, the output of each layer was clustered and projected back onto the
data space. Persistence diagrams on the clusters were calculated using Ripser.
Lastly, principal component analysis (PCA) was applied on the layer
representations corresponding to the torus to project them onto three
dimensions.

The experiment was then repeated using the Tanh activation function on the
network architecture described in table 1.

###  4.2 Results

The results from the experiments are summarized in figures 2 to 13. We notice
that the clustering becomes progressively less defined throughout the layers
as the network’s representation becomes more consolidated. Notably, despite
the visual resemblance to the data, persistent homology appeared unable to
calculate topological features with the same degree of accuracy as on the raw
data. The PCA projections obtained appear to resemble the original data less
as the layers become deeper. It is noted that we were unable to compute the
persistence diagrams of the clusterings of the third layer of the ReLU network
due to the these clusterings being extremely noisy.

Figure 2: First two clusters of the layer 1 representations of the ReLU
network.

((a)) Persistence diagrams of clusters.

((b)) PCA projection of torus.

Figure 3: Persistence diagrams of layer 1 clusters for the ReLU network and
PCA projection.

Figure 4: First two clusters of the layer 2 representations of the ReLU
network.

((a)) Persistence diagrams of clusters.

((b)) PCA projection of torus.

Figure 5: Persistence diagrams of layer 2 clusters for the ReLU network and
PCA projection.

Figure 6: First two clusters of the layer 3 representations of the ReLU
network. Figure 7: PCA projection of the layer 3 representation of the ReLU
network.

Figure 8: First two clusters of the layer 1 representations of the Tanh
network.

((a)) Persistence diagrams of clusters.

((b)) PCA projection of torus.

Figure 9: Persistence diagrams of layer 1 clusters for the Tanh network and
PCA projection.

Figure 10: First two clusters of the layer 2 representations of the Tanh
network.

((a)) Persistence diagrams of clusters.

((b)) PCA projection of torus.

Figure 11: Persistence diagrams of layer 2 clusters for the Tanh network and
PCA projection.

Figure 12: First two clusters of the layer 3 representations of the Tanh
network.

((a)) Persistence diagrams of clusters.

((b)) PCA projection of torus.

Figure 13: Persistence diagrams of layer 3 clusters for the Tanh network and
PCA projection.

##  5 Discussion and Further Work

The visual results from the experiments seem to indicate that the network is
approximating homeomorphisms in early layers, before deeper layer
representations obtain a topology markedly different than the original data.
This could result from repeated application of non-homeomorphic layers making
recovery of the underlying topology difficult. This appears to be correlated
with inability of HDBSCAN to properly cluster the layer representations of
deeper layers. However, we notice that while not as defined as the raw data,
the persistence diagrams indicate the presence of topological features similar
to the original data for the first two layers of both networks, albeit with a
much shorter lifetime due to interference from the noise of the clusterings.
This brings up the possibility that the network is attempting to isolate the
relevant topological features for classification, and excising superfluous
features. This would appear to explain the progressively sparser PCA
projections obtained as the layers get deeper. The deep PCA projections seem
to cluster points in more localized areas, leaving more empty portions in the
space unlike the spread-out representations of earlier layers.

Furthermore, the PCA projections seem to indicate that the network
approximates a continuous deformation for the first two layers of the Tanh
network and first layer of the ReLU network, before the topology becomes
unrecognizable. This could a result of ReLU being a surjective function,
rather than bijective function like Tanh, making it harder for the layer to
approximate a homeomorphism.

It is noted that nothing in this study is conclusive, and these are
observations gathered from a preliminary experiment. Further work must be done
in testing neural networks on more varied topologies and perhaps using
alternative clustering algorithms as the noise included in these clusters
appear to disrupt persistent homology.

## References

[1] Abadi, M. et al. (2015). TensorFlow: Large-Scale Machine Learning on
Heterogenous Systems. https://doi.org/10.5281/zenodo.4724125.

[2] Armstrong, M. A. (1983). Basic Topology. Springer,
doi:10.1007/978-1-4757-1793-8.

[3] Fefferman, C., Mitter, S., & Narayanan, H. (2013). Testing the Manifold
Hypothesis. arXiv, https://doi.org/10.48550/arXiv.1310.0425.

[4] McInnes, L., Healy, J., & Astels, S. (2017). hdbscan: Hierarchical density
based clustering. The Journal of Open Source Software, 2.
https://doi.org/10.21105/joss.00205.

[5] Montúfar, G., Otter, N., & Wang, Y. (2020). Can neural networks learn
persistent homology features?. arxiv,
https://doi.org/10.48550/arXiv.2011.14688.

[6] Tralie et al. (2018). Ripser.py: A Lean Persistent Homology Libary for
Python. Journal of Open Software, 3(29), 925,
https://doi.org/10.21105/joss.00925.