# Adaptive Smoothing Spline Estimator for the Function-on-Function Linear
Regression Model

Fabio Centofanti<EMAIL_ADDRESS>[ Department of Industrial
Engineering, University of Naples Federico II, Piazzale Tecchio 80, 80125,
Naples, Italy
Antonio Lepore<EMAIL_ADDRESS>[ Department of Industrial
Engineering, University of Naples Federico II, Piazzale Tecchio 80, 80125,
Naples, Italy
Alessandra Menafoglio<EMAIL_ADDRESS>[ MOX -
Modelling and Scientific Computing, Department of Mathematics, Politecnico di
Milano, Piazza Leonardo da Vinci 32, 20133, Milan, Italy
Biagio Palumbo<EMAIL_ADDRESS>[ Department of Industrial
Engineering, University of Naples Federico II, Piazzale Tecchio 80, 80125,
Naples, Italy
Simone Vantini<EMAIL_ADDRESS>[ MOX - Modelling and
Scientific Computing, Department of Mathematics, Politecnico di Milano, Piazza
Leonardo da Vinci 32, 20133, Milan, Italy
Some University and Another University

(2020)

###### Abstract

In this paper, we propose an adaptive smoothing spline (AdaSS) estimator for
the function-on-function linear regression model where each value of the
response, at any domain point, depends on the full trajectory of the
predictor. The AdaSS estimator is obtained by the optimization of an objective
function with two spatially adaptive penalties, based on initial estimates of
the partial derivatives of the regression coefficient function. This allows
the proposed estimator to adapt more easily to the true coefficient function
over regions of large curvature and not to be undersmoothed over the remaining
part of the domain. A novel evolutionary algorithm is developed ad hoc to
obtain the optimization tuning parameters. Extensive Monte Carlo simulations
have been carried out to compare the AdaSS estimator with competitors that
have already appeared in the literature before. The results show that our
proposal mostly outperforms the competitor in terms of estimation and
prediction accuracy. Lastly, those advantages are illustrated also on two
real-data benchmark examples.

62R10,

62J05 ,

62G08,

Functional data analysis,

Function-on-function linear regression,

Adaptive smoothing,

Functional regression,

###### keywords:

[class=MSC]

###### keywords:

††volume: 0††issue: 0

T1Corresponding Author

###### Contents

  1. 1 Introduction
  2. 2 The Adaptive Smoothing Spline Estimator
    1. 2.1 The Estimator
    2. 2.2 The Derivation of the AdaSS Estimator
    3. 2.3 The Algorithm for the Parameter Selection
  3. 3 Simulation Study
    1. 3.1 Mexican Hat Function
    2. 3.2 Dampened Harmonic Motion Function
    3. 3.3 Rapid Change Function
  4. 4 Real-data Examples
    1. 4.1 Swedish Mortality Dataset
    2. 4.2 Ship CO2 Emission Dataset
  5. 5 Conclusion
  6. A Approximation of the Two Penalty Terms for the AdaSS Estimator Derivation

##  1 Introduction

Complex datasets are increasingly available due to advancements in technology
and computational power and have stimulated significant methodological
developments. In this regard, functional data analysis (FDA) addresses the
issue of dealing with data that can be modeled as functions defined on a
compact domain. FDA is a thriving area of statistics and, for a comprehensive
overview, the reader could refer to [30, 18, 17, 21, 12]. In particular, the
generalization of the classical multivariate regression analysis to the case
where the predictor and/or the response have a functional form is referred to
as functional regression and is illustrated e.g., in [25] and [30]. Most of
the functional regression methods have been developed for models with scalar
response and functional predictors (scalar-on-function regression) or
functional response and scalar predictors (function-on-scalar regression).
Some results may be found in [8, 20, 39, 26]. Models where both the response
and the predictor are functions, namely function-on-function (FoF) regression,
have been far less studied until now. In this work, we study FoF linear
regression models, where the response variable function, at any domain point,
depends linearly on the full trajectory of the predictor. That is,

$Y_{i}\left(t\right)=\int_{\mathcal{S}}X_{i}\left(s\right)\beta\left(s,t\right)ds+\varepsilon_{i}\left(t\right)\quad
t\in\mathcal{T},$ (1.1)

for $i=1,\dots,n$. The pairs $\left(X_{i},Y_{i}\right)$ are independent
realizations of the predictor $X$ and the response $Y$, which are assumed to
be smooth random process with realizations in $L^{2}(\mathcal{S})$ and
$L^{2}(\mathcal{T})$, i.e., the Hilbert spaces of square integrable functions
defined on the compact sets $\mathcal{S}$ and $\mathcal{T}$, respectively.
Without loss of generality, the latter are also assumed with functional mean
equal to zero. The functions $\varepsilon_{i}$ are zero-mean random errors,
independent of $X_{i}$. The function $\beta$ is smooth in
$L^{2}(\mathcal{S}\times\mathcal{T})$, i.e., the Hilbert space of bivariate
square integrable functions defined on the closed intervals
$\mathcal{S}\times\mathcal{T}$, and is hereinafter referred to as coefficient
function. For each $t\in\mathcal{T}$, the contribution of
$X_{i}\left(\cdot\right)$ to the conditional value of $Y_{i}\left(t\right)$ is
generated by $\beta\left(\cdot,t\right)$, which works as continuous set of
weights of the predictor evaluations. Different methods to estimate $\beta$ in
(1.1) have been proposed in the literature. Ramsay and Silverman [30] assume
the estimator of $\beta$ to be in a finite dimension tensor space spanned by
two basis sets and where regularization is achieved by either truncation or
roughness penalties. (The latter is the foundation of the method proposed in
this article as we will see below.) Yao et al. [41] assume the estimator of
$\beta$ to be in a tensor product space generated by the eigenfunctions of the
covariance functions of the predictor $X$ and the response $Y$, estimated by
using the principal analysis by conditional expectation (PACE) method [40].
More recently, Luo and Qi [23] propose an estimation method of the FoF linear
model with multiple functional predictors based on a finite-dimensional
approximation of the mean response obtained by solving a penalized generalized
functional eigenvalue problem. Qi and Luo [28] generalize the method in [23]
to the high dimensional case, where the number of covariates is much larger
than the sample size (i.e., $p>>n$). In order to improve model flexibility and
prediction accuracy, Luo and Qi [24] consider a FoF regression model with
interaction and quadratic effects. A nonlinear FoF additive regression model
with multiple functional predictors is proposed by Qi and Luo [29].

One of the most used estimation method is the smoothing spline estimator
$\hat{\beta}_{SS}$ introduced by Ramsay and Silverman [30]. It is obtained as
the solution of the following optimization problem

$\hat{\beta}_{SS}=\operatorname*{argmin}_{\alpha\in\mathbb{S}_{k_{1},k_{2},M_{1},M_{2}}}\Big{\\{}\sum_{i=1}^{n}\int_{\mathcal{T}}\left[Y_{i}\left(t\right)-\int_{\mathcal{S}}X_{i}\left(s\right)\alpha\left(s,t\right)ds\right]^{2}dt\\\
\hskip
113.81102pt+\lambda_{s}\int_{\mathcal{S}}\int_{\mathcal{T}}\left(\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)\right)^{2}dsdt+\lambda_{t}\int_{\mathcal{S}}\int_{\mathcal{T}}\left(\mathcal{L}_{t}^{m_{t}}\alpha\left(s,t\right)\right)^{2}dsdt\Big{\\}},$
(1.2)

where $\mathbb{S}_{k_{1},k_{2},M_{1},M_{2}}$ is the tensor product space
generated by the sets of B-splines of orders $k_{1}$ and $k_{2}$ associated
with the non-decreasing sequences of $M_{1}+2$ and $M_{2}+2$ knots defined on
$\mathcal{S}$ and $\mathcal{T}$, respectively. The operators
$\mathcal{L}_{s}^{m_{s}}$ and $\mathcal{L}_{t}^{m_{t}}$, with $m_{s}\leq
k_{1}-1$ and $m_{t}\leq k_{2}-1$, are the $m_{s}$th and $m_{t}$th order linear
differential operators applied to $\alpha$ with respect to the variables $s$
and $t$, respectively. The two penalty terms on the right-hand side of (1.2)
measure the roughness of the function $\alpha$. The positive constants
$\lambda_{s}$ and $\lambda_{t}$ are generally referred to as roughness
parameters and trade off smoothness and goodness of fit of the estimator. The
higher their values, the smoother the estimator of the coefficient function.

Note that the two penalty terms on the right-side hand of (1.2) do not depend
on $s$ and $t$. Therefore, the estimator $\hat{\beta}_{SS}$ may suffer from
over and under smoothing when, for instance, the true coefficient function
$\beta$ is wiggly or peaked only in some parts of the domain. To solve this
problem, we consider two adaptive roughness parameters that are allowed to
vary on the domain $\mathcal{S}\times\mathcal{T}$. In this way, more flexible
estimators can be obtained to improve the estimation of the coefficient
function.

Methods that use adaptive roughness parameters are very popular and well
established in the field of nonparametric regression, and are referred to as
adaptive methods. In particular, the smoothing spline estimator for
nonparametric regression [36, 13, 11, 14] has been extended by different
authors to take into account the non-uniform smoothness along the domain of
the function to be estimated [33, 27, 35, 37, 38].

In this paper, a spatially adaptive estimator is proposed as the solutions of
the following minimization problem

$\operatorname*{argmin}_{\alpha\in\mathbb{S}_{k_{1},k_{2},M_{1},M_{2}}}\Big{\\{}\sum_{i=1}^{n}\int_{\mathcal{T}}\left[Y_{i}\left(t\right)-\int_{\mathcal{S}}X_{i}\left(s\right)\alpha\left(s,t\right)ds\right]^{2}dt\\\
\hskip
85.35826pt+\int_{\mathcal{S}}\int_{\mathcal{T}}\lambda_{s}\left(s,t\right)\left(\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)\right)^{2}dsdt+\int_{\mathcal{S}}\int_{\mathcal{T}}\lambda_{t}\left(s,t\right)\left(\mathcal{L}_{t}^{m_{t}}\alpha\left(s,t\right)\right)^{2}dsdt\Big{\\}},$
(1.3)

where the two roughness parameters $\lambda_{s}\left(s,t\right)$ and
$\lambda_{t}\left(s,t\right)$ are functions that produce different amount of
penalty, and, thus, allow the estimator to spatially adapt, i.e., to
accommodate varying degrees of roughness over the domain
$\mathcal{S}\times\mathcal{T}$. Therefore, the model may accommodate the local
behavior of $\beta$ by imposing a heavier penalty in regions of lower
smoothness. Because $\lambda_{s}\left(s,t\right)$ and
$\lambda_{t}\left(s,t\right)$ are intrinsically infinite dimensional, their
specification could be rather complicated without further assumptions.

The proposed estimator is applied to FoF linear regression model reported in
(1.1), and is referred to as adaptive smoothing spline (AdaSS) estimator. It
is obtained as the solution of the optimization problem in (1.3), with
$\lambda_{s}\left(s,t\right)$ and $\lambda_{t}\left(s,t\right)$ chosen based
on an initial estimate of the partial derivatives
$\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)$ and
$\mathcal{L}_{t}^{m_{t}}\alpha\left(s,t\right)$. The rationale behind this
choice is to allow the contribution of $\lambda_{s}\left(s,t\right)$ and
$\lambda_{t}\left(s,t\right)$, to the penalties in (1.3), to be small over
regions where the initial estimate has large $m_{s}$th and $m_{t}$th
curvatures (i.e., partial derivatives), respectively. This can be regarded as
an extension to the FoF linear regression model of the idea of Storlie et al.
[35] and Abramovich and Steinberg [1]. Moreover, to overcome some limitations
of the most famous grid-search method [4], a new evolutionary algorithm is
proposed for the choice of the unknown parameters, needed to compute the AdaSS
estimator.

The rest of the paper is organized as follows. In Section 2.1, the proposed
estimator is presented. Computational issues involved in the AdaSS estimator
calculation are discussed in Section 2.2 and Section 2.3. In Section 3, by
means of a Monte Carlo simulation study, the performance of the proposed
estimator are compared with those achieved by competing estimators already
appeared in the literature. Lastly, two real-data examples are presented in
Section 4 to illustrate the practical applicability of the proposed estimator.
The conclusion is in Section 5.

##  2 The Adaptive Smoothing Spline Estimator

###  2.1 The Estimator

The AdaSS estimator $\hat{\beta}_{AdaSS}$ is defined as the solution of the
optimization problem in (1.3) where the two roughness parameters
$\lambda_{s}\left(s,t\right)$ and $\lambda_{t}\left(s,t\right)$ are as follows

$\displaystyle\lambda_{s}\left(s,t\right)=\lambda^{AdaSS}_{s}\frac{1}{\left(|\widehat{\beta_{s}^{m_{s}}}\left(s,t\right)|+\delta_{s}\right)^{\gamma_{s}}}$
$\displaystyle\lambda_{t}\left(s,t\right)=\lambda^{AdaSS}_{t}\frac{1}{\left(|\widehat{\beta_{t}^{m_{t}}}\left(s,t\right)|+\delta_{t}\right)^{\gamma_{t}}}$

that is,

$\hat{\beta}_{AdaSS}=\operatorname*{argmin}_{\alpha\in\mathbb{S}_{k_{1},k_{2},M_{1},M_{2}}}\Big{\\{}\sum_{i=1}^{n}\int_{\mathcal{T}}\left[Y_{i}\left(t\right)-\int_{\mathcal{S}}X_{i}\left(s\right)\alpha\left(s,t\right)ds\right]^{2}dt\\\
+\lambda^{AdaSS}_{s}\int_{\mathcal{S}}\int_{\mathcal{T}}\frac{1}{\left(|\widehat{\beta_{s}^{m_{s}}}\left(s,t\right)|+\delta_{s}\right)^{\gamma_{s}}}\left(\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)\right)^{2}dsdt\\\
+\lambda^{AdaSS}_{t}\int_{\mathcal{S}}\int_{\mathcal{T}}\frac{1}{\left(|\widehat{\beta_{t}^{m_{t}}}\left(s,t\right)|+\delta_{t}\right)^{\gamma_{t}}}\left(\mathcal{L}_{t}^{m_{t}}\alpha\left(s,t\right)\right)^{2}dsdt\Big{\\}},$
(2.1)

for some tuning parameters
$\lambda^{AdaSS}_{s},\delta_{s},\gamma_{s},\lambda^{AdaSS}_{t},\delta_{t},\gamma_{t}\geq
0$ and $\widehat{\beta_{s}^{m_{s}}}$ and $\widehat{\beta_{t}^{m_{t}}}$ initial
estimates of $\mathcal{L}_{s}^{m_{s}}\beta$ and
$\mathcal{L}_{t}^{m_{t}}\beta$, respectively. Note that the two roughness
parameters $\lambda_{s}$ and $\lambda_{t}$ assume large values over domain
regions where $\widehat{\beta_{s}^{m_{s}}}$ and $\widehat{\beta_{t}^{m_{t}}}$
are small. Therefore, in the right-hand side of (2.1),
$(\mathcal{L}_{s}^{m_{s}}\alpha)^{2}$ and
$(\mathcal{L}_{t}^{m_{t}}\alpha)^{2}$ are weighted through the inverse of
$|\widehat{\beta_{s}^{m_{s}}}|$ and $|\widehat{\beta_{t}^{m_{t}}}|$. That is,
over domain regions where $\widehat{\beta_{s}^{m_{s}}}$ and
$\widehat{\beta_{t}^{m_{t}}}$ are small, $(\mathcal{L}_{s}^{m_{s}}\alpha)^{2}$
and $(\mathcal{L}_{t}^{m_{t}}\alpha)^{2}$ have larger weights than over those
regions where $\widehat{\beta_{s}^{m_{s}}}$ and $\widehat{\beta_{t}^{m_{t}}}$
are large. For this reasons, the final estimator is able to adapt to the
coefficient function over regions of large curvature without over smoothing it
over regions where the $m_{s}$th and $m_{t}$th curvatures are small.

The constants $\delta_{s}$ and $\delta_{t}$ allow $\hat{\beta}_{AdaSS}$ not to
have $m_{s}$th and $m_{t}$th-order inflection points at the same location of
$\widehat{\beta_{s}^{m_{s}}}$ and $\widehat{\beta_{t}^{m_{t}}}$, respectively.
Indeed, when $\delta_{s}$ and $\delta_{t}$ are set to zero, where
$\widehat{\beta_{s}^{m_{s}}}=0$ and $\widehat{\beta_{t}^{m_{t}}}=0$ ($m_{s}$th
and $m_{t}$th-order inflection points), the corresponding penalties go to
infinite, and, thus, $\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)$ and
$\mathcal{L}_{t}^{m_{t}}\alpha\left(s,t\right)$ become zero in accordance with
the minimization problem. Therefore, the presence of $\delta_{s}$ and
$\delta_{t}$ makes $\hat{\beta}_{AdaSS}$ more robust against the choice of the
initial estimate of the linear differential operators applied to $\beta$ with
respect to $s$ and $t$. Finally, $\gamma_{s}$ and $\gamma_{t}$ control the
amount of weight placed in $\widehat{\beta_{s}^{m_{s}}}$ and
$\widehat{\beta_{t}^{m_{t}}}$, whereas $\lambda^{AdaSS}_{s}$ and
$\lambda^{AdaSS}_{t}$ are smoothing parameters. The solution of the
optimization problem in (2.1) can be obtained in closed form if the penalty
terms are approximated as described in Section 2.2. There are several choices
for the initial estimates $\widehat{\beta_{s}^{m_{s}}}$ and
$\widehat{\beta_{t}^{m_{t}}}$. As in [1], we suggest to apply the $m_{s}$th
and $m_{t}$th order linear differential operator to the smoothing spline
estimator $\hat{\beta}_{SS}$ in (1.2).

###  2.2 The Derivation of the AdaSS Estimator

The minimization in (2.1) is carried out over
$\alpha\in\mathbb{S}_{k_{1},k_{2},M_{1},M_{2}}$. This implicitly means that we
are approximating $\beta$ as follows

$\beta\left(s,t\right)\approx\tilde{\beta}\left(s,t\right)=\sum_{i=1}^{M_{1}+k_{1}}\sum_{j=1}^{M_{2}+k_{2}}b_{ij}\psi^{s}_{i}\left(s\right)\psi^{t}_{j}\left(t\right)=\bm{\psi}^{s}\left(s\right)^{T}\bm{B}\bm{\psi}^{t}\left(t\right)\quad
s\in\mathcal{S},t\in\mathcal{T},$ (2.2)

where $\bm{B}=\\{b_{ij}\\}\in\mathbb{R}^{M_{1}+k_{1}\times M_{2}+k_{2}}$. The
two sets
$\bm{\psi}^{s}=\left(\psi^{s}_{1},\dots,\psi^{s}_{M_{1}+k_{1}}\right)^{T}$ and
$\bm{\psi}^{t}=\left(\psi^{t}_{1},\dots,\psi^{t}_{M_{2}+k_{2}}\right)^{T}$ are
B-spline functions of order $k_{1}$ and $k_{2}$ and non-decreasing knots
sequences $\Delta^{s}=\\{s_{0},s_{1},\dots,s_{M_{1}},s_{M_{1}+1}\\}$ and
$\Delta^{t}=\\{t_{0},t_{1},\dots,t_{M_{2}},t_{M_{2}+1}\\}$, defined on
$\mathcal{S}=\left[s_{0},s_{M_{1}+1}\right]$ and
$\mathcal{T}=\left[t_{0},t_{M_{2}+1}\right]$, respectively, that generate
$\mathbb{S}_{k_{1},k_{2},M_{1},M_{2}}$. Thus, estimating $\beta$ in (2.1)
means estimating $\bm{B}$. Let
$\alpha\left(s,t\right)=\bm{\psi}^{s}\left(s\right)^{T}\bm{B}_{\alpha}\bm{\psi}^{t}\left(t\right)$,
$s\in\mathcal{S},t\in\mathcal{T}$, in $\mathbb{S}_{k_{1},k_{2},M_{1},M_{2}}$,
where $\bm{B}_{\alpha}=\\{b_{\alpha,ij}\\}\in\mathbb{R}^{M_{1}+k_{1}\times
M_{2}+k_{2}}$. Then, the first term of the right-hand side of (2.1) may be
rewritten as (see [30], pag 291-293, for the derivation)

$\sum_{i=1}^{n}\int_{\mathcal{T}}\left[Y_{i}\left(t\right)-\int_{\mathcal{S}}X_{i}\left(s\right)\alpha\left(s,t\right)ds\right]^{2}dt=\sum_{i=1}^{n}\int_{\mathcal{T}}Y_{i}\left(t\right)^{2}dt-2\operatorname{Tr}\left[\bm{X}\bm{B}_{\alpha}\bm{Y}^{T}\right]+\operatorname{Tr}\left[\bm{X}^{T}\bm{X}\bm{B}_{\alpha}\bm{W}_{t}\bm{B}_{\alpha}^{T}\right],$
(2.3)

where $\bm{X}=\left(\bm{X}_{1},\dots,\bm{X}_{n}\right)^{T}$, with
$\bm{X}_{i}=\int_{\mathcal{S}}X_{i}\left(s\right)\bm{\psi}^{s}\left(s\right)ds$,
$\bm{Y}=\left(\bm{Y}_{1},\dots,\bm{Y}_{n}\right)^{T}$ with
$\bm{Y}_{i}=\int_{\mathcal{T}}Y_{i}\left(t\right)\bm{\psi}^{t}\left(t\right)dt$,
and
$\bm{W}_{t}=\int_{\mathcal{T}}\bm{\psi}^{t}\left(t\right)\bm{\psi}^{t}\left(t\right)^{T}dt$.
The term $\operatorname{Tr}\left[\bm{A}\right]$ denotes the trace of a square
matrix $\bm{A}$.

In order to simplify the integrals in the two penalty terms on the right-hand
side of (2.1), and thus obtain a linear form in $\bm{B}_{\alpha}$, we
consider, for $s\in\mathcal{S}$ and $t\in\mathcal{T}$, the following
approximations of $\widehat{\beta_{s}^{m_{s}}}$ and
$\widehat{\beta_{t}^{m_{t}}}$

$\widehat{\beta_{s}^{m_{s}}}\left(s,t\right)\approx\sum_{i=0}^{L_{s}}\sum_{j=0}^{L_{t}}\widehat{\beta_{s}^{m_{s}}}\left(\tau_{s,i+1},\tau_{t,j+1}\right)I_{\left[\left(\tau_{s,i},\tau_{s,i+1}\right)\times\left(\tau_{t,j},\tau_{t,j+1}\right)\right]}\left(s,t\right),$
(2.4)

and

$\widehat{\beta_{t}^{m_{t}}}\left(s,t\right)\approx\sum_{i=0}^{L_{s}}\sum_{j=0}^{L_{t}}\widehat{\beta_{t}^{m_{t}}}\left(\tau_{s,i+1},\tau_{t,j+1}\right)I_{\left[\left(\tau_{s,i},\tau_{s,i+1}\right)\times\left(\tau_{t,j},\tau_{t,j+1}\right)\right]}\left(s,t\right),$
(2.5)

where
$\Theta^{s}=\\{\tau_{s,0},\tau_{s,1},\dots\tau_{s,L_{s}},\tau_{s,L_{s}+1}\\}$
and
$\Theta^{t}=\\{\tau_{t,0},\tau_{t,1},\dots\tau_{t,L_{t}},\tau_{t,L_{t}+1}\\}$
are non increasing knot sequences with $\tau_{s,0}=s_{0}$,
$\tau_{s,L_{s}+1}=s_{M_{1}+1}$, $\tau_{t,0}=t_{0}$,
$\tau_{t,L_{t}+1}=t_{M_{2}+1}$, and $I_{\left[a\times
b\right]}\left(z_{1},z_{2}\right)=1$ for
$\left(z_{1},z_{2}\right)\in\left[a\times b\right]$ and zero elsewhere. In
(2.4) and (2.5), we are assuming that $\widehat{\beta_{s}^{m_{s}}}$ and
$\widehat{\beta_{t}^{m_{t}}}$ are well approximated by a piecewise constant
function, whose values are constant on rectangles defined by the two knot
sequences $\Theta^{s}$ and $\Theta^{t}$. It can be easily proved, by following
Schumaker [34] (pag. 491, Theorem 12.7), that the approximation error in both
cases goes to zero as the mesh widths
$\overline{\delta}^{s}=\max_{i}\left(\tau_{s,i+1}-\tau_{s,i}\right)$ and
$\overline{\delta}^{t}=\max_{j}\left(\tau_{t,j+1}-\tau_{t,g}\right)$ go to
zero. Therefore, $\widehat{\beta_{s}^{m_{s}}}$ and
$\widehat{\beta_{t}^{m_{t}}}$ can be exactly recovered by uniformly increasing
the number of knots $L_{s}$ and $L_{t}$. In this way, the two penalties on the
right-hand side of (2.1) can be rewritten as (A)

$\lambda^{AdaSS}_{s}\int_{\mathcal{S}}\int_{\mathcal{T}}\frac{1}{\left(|\widehat{\beta_{s}^{m_{s}}}\left(s,t\right)|+\delta_{s}\right)^{\gamma_{s}}}\left(\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)\right)^{2}dsdt\\\
\hskip
56.9055pt\approx\lambda^{AdaSS}_{s}\sum_{i=1}^{L_{s}+1}\sum_{j=1}^{L_{t}+1}d^{s}_{ij}\operatorname{Tr}\left[\bm{B}_{\alpha}^{T}\bm{R}_{s,i}\bm{B}_{\alpha}\bm{W}_{t,j}\right]$
(2.6)

and

$\lambda^{AdaSS}_{t}\int_{\mathcal{S}}\int_{\mathcal{T}}\frac{1}{\left(|\widehat{\beta_{s}^{m_{s}}}\left(s,t\right)|+\delta_{t}\right)^{\gamma_{t}}}\left(\mathcal{L}_{t}^{m_{t}}\alpha\left(s,t\right)\right)^{2}dsdt\\\
\hskip
56.9055pt\approx\lambda^{AdaSS}_{t}\sum_{i=1}^{L_{s}+1}\sum_{j=1}^{L_{t}+1}d^{t}_{ij}\operatorname{Tr}\left[\bm{B}_{\alpha}^{T}\bm{W}_{s,i}\bm{B}_{\alpha}\bm{R}_{t,j}\right],$
(2.7)

where
$\bm{W}_{s,i}=\int_{\left[\tau_{s,i-1},\tau_{s,i}\right]}\bm{\psi}^{s}\left(s\right)\bm{\psi}^{s}\left(s\right)^{T}ds$,
$\bm{W}_{t,j}=\int_{\left[\tau_{t,j-1},\tau_{t,j}\right]}\bm{\psi}^{t}\left(t\right)\bm{\psi}^{t}\left(t\right)^{T}dt$,
$\bm{R}_{s,i}=\int_{\left[\tau_{s,i-1},\tau_{s,i}\right]}\mathcal{L}_{s}^{m_{s}}\left[\bm{\psi}^{s}\left(s\right)\right]\mathcal{L}_{s}^{m_{s}}\left[\bm{\psi}^{s}\left(s\right)\right]^{T}ds$
and
$\bm{R}_{t,j}=\int_{\left[\tau_{t,j-1},\tau_{t,j}\right]}\mathcal{L}_{t}^{m_{t}}\left[\bm{\psi}^{t}\left(t\right)\right]\mathcal{L}_{t}^{m_{t}}\left[\bm{\psi}^{t}\left(t\right)\right]^{T}dt$,
and
$d^{s}_{ij}=\Big{\\{}\frac{1}{\left(|\widehat{\beta_{s}^{m_{s}}}\left(\tau_{s,i},\tau_{t,j}\right)|+\delta_{s}\right)^{\gamma_{s}}}\Big{\\}}$
and
$d^{t}_{ij}=\Big{\\{}\frac{1}{\left(|\widehat{\beta_{t}^{m_{t}}}\left(\tau_{s,i},\tau_{t,j}\right)|+\delta_{t}\right)^{\gamma_{t}}}\Big{\\}}$,
for $i=1,\dots,L_{s}+1$ and $j=1,\dots,L_{t}+1$.

The optimization problem in (2.1) can be then approximated with the following

$\displaystyle\hat{\bm{B}}_{AS}\approx\operatorname*{argmin}_{\bm{B}_{\alpha}\in\mathbb{R}^{\left(M_{1}+k_{1}\right)\times\left(M_{2}+k_{2}\right)}}\Big{\\{}$
$\displaystyle\sum_{i=1}^{n}\int_{\mathcal{T}}Y_{i}\left(t\right)^{2}dt-2\operatorname{Tr}\left[\bm{X}\bm{B}_{\alpha}\bm{Y}^{T}\right]+\operatorname{Tr}\left[\bm{X}^{T}\bm{X}\bm{B}_{\alpha}\bm{W}_{t}\bm{B}_{\alpha}^{T}\right]$
$\displaystyle+\sum_{i=1}^{L_{s}+1}\sum_{j=1}^{L_{t}+1}\left(\lambda^{AdaSS}_{s}d^{s}_{ij}\operatorname{Tr}\left[\bm{B}_{\alpha}^{T}\bm{R}_{s,i}\bm{B}_{\alpha}\bm{W}_{t,j}\right]+\lambda^{AdaSS}_{t}d^{t}_{ij}\operatorname{Tr}\left[\bm{B}_{\alpha}^{T}\bm{W}_{s,i}\bm{B}_{\alpha}\bm{R}_{t,j}\right]\right)\Big{\\}},$
(2.8)

or by vectorization as

$\displaystyle\hat{\bm{b}}_{AS}\approx\operatorname*{argmin}_{\bm{b}_{\alpha}\in\mathbb{R}^{\left(M_{1}+k_{1}\right)\left(M_{2}+k_{2}\right)}}\Big{\\{}$
$\displaystyle-2\operatorname*{vec}\left(\bm{X}^{T}\bm{Y}\right)^{T}\bm{b}_{\alpha}+\bm{b}_{\alpha}^{T}\left(\bm{W}_{t}\otimes\bm{X}^{T}\bm{X}\right)\bm{b}_{\alpha}$
$\displaystyle\sum_{i=1}^{L_{s}+1}\sum_{j=1}^{L_{t}+1}\left(\lambda^{AdaSS}_{s}d^{s}_{ij}\bm{b}_{\alpha}^{T}\bm{L}_{wr,ij}\bm{b}_{\alpha}+\lambda^{AdaSS}_{t}d^{t}_{ij}\bm{b}_{\alpha}^{T}\bm{L}_{rw,ij}\bm{b}_{\alpha}\right)\bigg{\\}},$
(2.9)

where $\hat{\bm{b}}_{AS}=\operatorname*{vec}\left(\hat{\bm{B}}_{AS}\right)$,
$\bm{L}_{rw,ij}=\left(\bm{R}_{t,j}\otimes\bm{W}_{s,i}\right)$ and
$\bm{L}_{wr,ij}=\left(\bm{W}_{t,j}\otimes\bm{R}_{s,i}\right)$, for
$i=1,\dots,L_{s}+1$ and $j=1,\dots,L_{t}+1$. For a matrix
$\bm{A}\in\mathbb{R}^{j\times k}$, $\operatorname*{vec}(\bm{A})$ indicates the
vector of length $jk$ obtained by writing the matrix $\bm{A}$ as a vector
column-wise, and $\otimes$ is the Kronecker product. Because the matrices
$\bm{W}_{t}$, $\bm{L}_{wr,ij}$ and $\bm{L}_{rw,ij}$ for $i=1,\dots,L_{s}+1$
and $j=1,\dots,L_{t}+1$ are positive definite and by assuming that
$\bm{X}^{T}\bm{X}$ is positive definite, then the minimizer of the
optimization problem in (2.2) exists, is unique and has the following
expression [5]

$\displaystyle\hat{\bm{b}}_{AdaSS}\approx\left[\left(\bm{W}_{t}\otimes\bm{X}^{T}\bm{X}\right)+\sum_{i=1}^{L_{s}+1}\sum_{j=1}^{L_{t}+1}\left(\lambda^{AdaSS}_{s}d^{s}_{ij}\bm{L}_{wr,ij}+\lambda^{AdaSS}_{t}d^{t}_{ij}\bm{L}_{rw,ij}\right)\right]^{-1}\operatorname*{vec}\left(\bm{X}^{T}\bm{Y}\right).$
(2.10)

To obtain $\hat{\bm{b}}_{AdaSS}$ in (2.10) the tuning parameters
$\lambda^{AdaSS}_{s},\delta_{s},\gamma_{s},\lambda^{AdaSS}_{t},\delta_{t},\gamma_{t}$
must be opportunely chosen. This issue is discussed in Section 2.3.

###  2.3 The Algorithm for the Parameter Selection

There are some tuning parameters in the optimization problem (2.2) that must
be chosen to obtain the AdaSS estimator. Usually, the tensor product space
$\mathbb{S}_{k_{1},k_{2},M_{1},M_{2}}$ is chosen with $k_{1}=k_{2}=4$, i.e.,
cubic B-splines, and equally spaced knot sequences. Although the choice of
$M_{1}$ and $M_{2}$ is not crucial [8], it should allow the final estimator to
capture the local behaviour of the coefficient function $\beta$, that is,
$M_{1}$ and $M_{2}$ should be sufficiently large. The smoothness of the final
estimator is controlled by the two penalty terms on the right-hand side of
(2.2).

The tuning parameters
$\lambda^{AdaSS}_{s},\delta_{s},\gamma_{s},\lambda^{AdaSS}_{t},\delta_{t},\gamma_{t}$
could be fixed by using the conventional $K$-fold cross validation (CV) [15],
where the combination of parameters to be explored is chosen by means of the
classic grid search method [15]. That is an exhaustive searching through a
manually specified subset of the tuning parameter space [3]. Although, in our
setting, grid search is embarrassingly parallel [16], it is not scalable
because it suffers from the curse of dimensionality. However, even if this is
beyond the scope of the present work, note that the number of combinations to
explore grows exponentially with the number of tuning parameters and makes
unsuitable the application of the proposed method to the FoF linear model in
the case of multiple predictors. Then, to facilitate the use of the proposed
method by practitioners, in what follows, we proposed a novel evolutionary
algorithm for tuning parameter selection, referred to as evolutionary
algorithm for adaptive smoothing estimator (EAASS) inspired by the population
based training (PBT) introduced by Jaderberg et al. [19]. The PBT algorithm
was introduced to address the issue of hyperparameter optimization for neural
networks. It bridges and extends parallel search method (e.g., grid search and
random search) with sequential optimization method (e.g., hand tuning and
Bayesian optimization). The former runs many parallel optimization processes,
for different combinations of hyperparameter values, and, then chooses the
combination that shows the best performance. The latter performs several steps
of few parallel optimizations, where, at each step, information coming from
the previous step is used to identify the combinations of hyperparameter
values to explore. For further details on the PBT algorithm the readers should
refer to [19], where the authors demonstrated its effectiveness and wide
applicability. The pseudo code of the EAASS algorithm is given in Algorithm 1.

Algorithm 1 EAASS algorithm

1:Choose the initial population $\mathcal{P}=\\{p_{i}\\}$ of combinations of
tuning parameter values

2:Obtain the set $\mathcal{V}=\\{v_{i}\\}$ of estimated prediction errors
corresponding to $\mathcal{P}$

3:repeat

4: Identify the set $\mathcal{Q}\subseteq\mathcal{P}$ and the corresponding
$\mathcal{Z}\subseteq\mathcal{V}$$\blacktriangleright$ exploitation

5: for $p_{i}\in\mathcal{Q}$ do$\blacktriangleright$ exploration

6: Obtain the new combination of tuning parameter values, $p^{\prime}_{i}$

7: Obtain the new estimated prediction error $v^{\prime}_{i}$ corresponding to
$p^{\prime}_{i}$

8: end for

9: Define $\mathcal{Q}^{\prime}=\\{p^{\prime}_{i}\\}$ and
$\mathcal{Z}^{\prime}=\\{v^{\prime}_{i}\\}$

10: Set $\mathcal{P}=\mathcal{P}\setminus\mathcal{Q}\cup\mathcal{Q}^{\prime}$
and $\mathcal{V}=\mathcal{V}\setminus\mathcal{Z}\cup\mathcal{Z}^{\prime}$

11:until The stopping condition is met

12:Return $p_{i}\in\mathcal{P}$ with the highest $v_{i}\in\mathcal{V}$

The first step is the identification of an initial population $\mathcal{P}$ of
tuning parameter combinations $p_{i}$s. This can be done, for each combination
and each tuning parameter, by randomly selecting a value in a pre-specified
range. Then, the set $\mathcal{V}$ of estimated prediction errors $v_{i}$s
corresponding to $\mathcal{P}$ is obtained by means of $K$-fold CV. We choose
a subset $\mathcal{Q}$ of $\mathcal{P}$, by following a given exploitation
strategy and, thus, the corresponding subset $\mathcal{Z}$ of $\mathcal{V}$. A
typical exploitation strategy is the truncation selection, where the worse
$r\%$, for $0\leq r\leq 100$, of $\mathcal{P}$, in terms of estimated
prediction error, is substituted by elements randomly sampled from the
remaining $(100-r)\%$ part of the current population [19]. Then the following
step consists of an exploration strategy where the tuning parameter
combinations in $\mathcal{Q}$ are substituted by new ones. The simulation
study in Section 3 and the real-data Examples in Section 4 are based on a
perturbation where each tuning parameter value of the given combination is
randomly perturbed by a factor of 1.2 or 0.8. The exploitation and exploration
phases are repeated until a stopping condition is met, e.g, maximum number of
iterations. Other exploration and exploitation strategies can be found in [2].
At last, the selected tuning parameter combination is obtained as an element
of $\mathcal{P}$ that achieves the lowest estimated prediction error. As a
remark, in our trials the AdaSS estimator works quite well with
$\delta_{s}=\delta_{s}^{*}\max|\widehat{\beta_{s}^{m_{s}}}\left(s,t\right)|$
and
$\delta_{t}=\delta_{t}^{*}\max|\widehat{\beta_{t}^{m_{t}}}\left(s,t\right)|$,
for $0\leq\delta_{s}^{*},\delta_{t}^{*}\leq 0.1$.

##  3 Simulation Study

In this section, the performance of the AdaSS estimator is assessed on several
simulated datasets. In particular, we compare the AdaSS estimator with cubic
B-splines and $m_{s}=m_{t}=2$ with five competing methods that represent the
state of the art in the FoF liner regression model estimation. The first two
are those proposed by Ramsay and Silverman [30]. The first one, hereinafter
referred to as SMOOTH estimator, is the smoothing spline estimator described
in (1.2), whereas, the second one, hereinafter referred to as TRU estimator,
assumes that the coefficient function is in a finite dimensional tensor
product space generate by two sets of B-splines with regularization achieved
by choosing the space dimension. Then, we consider also the estimator proposed
by Yao et al. [41] and Canale and Vantini [6]. The former is based on the
functional principal component decomposition, and is hereinafter referred to
as PCA estimator, while the latter relies on a ridge type penalization,
hereinafter referred to as RIDGE estimator. Lastly, as the fifth alternative,
we explore the estimator proposed by Luo and Qi [23], hereinafter referred to
as SIGCOMP. Moreover, the AdaSS estimator with cubic B-splines and
$m_{s}=m_{t}=2$ is considered. For illustrative purposes, we also consider a
version of the AdaSS estimator, referred to AdaSStrue, whose roughness
parameters are calculated by assuming that the true coefficient function is
known. Obviously, the AdaSStrue has not a practical meaning because the true
coefficient function is never known. However, it allows one to better
understand the influence of the initial estimates of the partial derivatives
on the AdaSS performance. All the unknown parameters of the competing methods
considered are chosen by means of $10$-fold CV. The tuning parameters of the
AdaSS and AdaSStrue estimators are chosen through the EAASS algorithm. The set
$\mathcal{P}$ is obtained by using $10$-fold CV, the exploitation and
exploration phases are as described in Section 2.3 and a maximum number of
iterations equal to 15 is set as stopping condition. For each simulation, a
training sample of $n$ observations is generated along with a test set $T$ of
size $N=4000$. They are used to estimate $\beta$ and to test the predictive
performance of the estimated model, respectively. Three different sample sizes
are considered, viz., $n=100,500,1000$. The estimation accuracy of the
estimators are assessed by using the integrated squared error  (ISE) defined
as

$\text{ISE}=\frac{1}{A}\int_{\mathcal{S}}\int_{\mathcal{T}}\left(\hat{\beta}\left(s,t\right)-\beta\left(s,t\right)\right)^{2}dsdt,$
(3.1)

where $A$ is the measure of $\mathcal{S}\times\mathcal{T}$. The ISE aims to
measure the estimation error of $\hat{\beta}$ with respect to $\beta$.
Whereas, the predictive accuracy is measured through the prediction mean
squared error (PMSE) defined as

$\text{PMSE}=\frac{1}{N}\sum_{\left(X,Y\right)\in
T}\int_{\mathcal{T}}\left(Y\left(t\right)-\int_{\mathcal{S}}X\left(s\right)\hat{\beta}\left(s,t\right)ds\right)^{2}dt.$
(3.2)

The observations in the test set are centred by subtracting to each
observation the corresponding sample mean function estimated in the training
set. The observations in the training and test sets are obtained as follows.
The covariate $X_{i}$ and the errors $\varepsilon_{i}$ are generated as linear
combination of cubic B-splines, $\Psi_{i}^{x}$ and $\Psi_{i}^{\varepsilon}$,
with evenly spaced knots, i.e., $X_{i}=\sum_{j=1}^{32}x_{ij}\Psi_{i}^{x}$ and
$\varepsilon_{i}=k\sum_{j=1}^{20}e_{ij}\Psi_{i}^{\varepsilon}$. The
coefficients $x_{ij}$ and $e_{ij}$, for $i=1,\dots,n$, $j=1,\dots,32$ and
$j=1,\dots,20$, are independent realizations of standard normal random
variable and the numbers of basis have been randomly chosen between 10 and 50.
The constant $k$ is chosen such that the signal-to-noise ratio
$SN\doteq\int_{\mathcal{T}}\operatorname{Var}_{X}[\operatorname{E}\left(Y_{i}|X_{i}\right)]/\int_{\mathcal{T}}\operatorname{Var}\left(\varepsilon_{i}\right)$
is equal to 4, where $\operatorname{Var}_{X}$ is the variance with respect to
the random covariate $X$. Then, given the coefficient function $\beta$, the
response $Y_{i}$ is obtained.

###  3.1 Mexican Hat Function

The Mexican hat function is a linear function with a sharp smoothness
variation in central part of the domain. In this case, the coefficient
function $\beta$ is defined as

$\beta\left(s,t\right)=-1+1.5s+1.5t+0.05\phi\left(s,t\right),\quad
s,t\in\left[0,1\right]\times\left[0,1\right]$

where $\phi$ is a multivariate normal distribution with mean
$\bm{\mu}=\left(0.6,0.6\right)^{T}$ and diagonal covariance matrix
$\bm{\Sigma}=\operatorname{diag}\left(0.001,0.001\right)$. Figure 1 displays
the AdaSS and the SMOOTH estimates along with the true coefficient function
for a randomly selected simulation run.

Figure 1:  AdaSS (solid line) and SMOOTH (dashed line) estimates of the
coefficient functions and the TRUE coefficient function $\beta$ (dotted line)
for different values of $t$ in the case of the Mexican hat function.

The proposed estimator tends to be smoother on the flat region and is able to
better capture the peak in the coefficient function (at $t\approx 0.6$) than
the SMOOTH estimate. The latter, to perform reasonably well along the whole
domain, selects tuning parameters that are not sufficiently small (large) on
the peaky (flat) region. This is also confirmed by the graphical appeal of the
AdaSS estimate with respect to the competitor ones. In Figure 2 and top of
Table 1, the values of ISE and PMSE achieved by the AdaSS, AdaSStrue, and
competitor estimators are shown as functions of the sample size $n$. Without
considering the AdaSStrue estimator, the AdaSS estimator yields the lowest ISE
for all sample sizes, and thus has the lowest estimation error. In terms of
PMSE, it is the best one for $n=150$, whereas for $n=500,1000$ it performs
comparably with SIGCOMP and PCA estimators. The performance of the AdaSStrue
and AdaSS estimators is very similar in terms of ISE, whereas the AdaSStrue
shows a lower PMSE. However, as expected, the effect of the knowledge of the
true coefficient function tends to disappear as $n$ increases, because the
partial derivatives estimates become more accurate.

(a)

(b)

Figure 2:  4(a) The integrated squared error (ISE) and 2(b) the prediction
mean squared error (PMSE) $\pm standard\hskip 3.41418pterror$ for the TRU,
SMOOTH, PCA, RIDGE, SIGCOMP, AdaSS and AdaSStrue estimators in the case of the
Mexican hat function.  Table 1: The integrated squared error (ISE) and the
prediction mean squared error (PMSE) for the TRU, SMOOTH, PCA, RIDGE, SIGCOMP,
AdaSS and AdaSStrue estimators. The numbers outside the parentheses are the
averages over $100$ Monte Carlo replications, and the numbers inside
parentheses are the corresponding standard errors. The values corresponding to
the AdaSStrue estimator are emphasized to underline the fact that they rely on
the knowledge of the true coefficient function, which is unlikely in real
applications. In bold are marked the lowest values among the AdaSS and the
competitors.

| $n=100$ | $n=500$ | $n=1000$
---|---|---|---
| ISE ($\times 10^{-1}$) | PMSE ($\times 10^{-2}$) | ISE ($\times 10^{-1}$) | PMSE ($\times 10^{-2}$) | ISE ($\times 10^{-1}$) | PMSE ($\times 10^{-2}$)
| Mexican hat
TRU | 0.4063(0.0059) | 0.3575(0.0011) | 0.1384(0.0020) | 0.3143(0.0007) | 0.0660(0.0011) | 0.3031 (0.0005)
SMOOTH | 0.2191(0.0020) | 0.3382(0.0007) | 0.0917(0.0008) | 0.3088(0.0005) | 0.0564(0.0006) | 0.3027 (0.0005)
PCA | 0.2519(0.0068) | 0.3234(0.0007) | 0.0681(0.0013) | 0.3030(0.0005) | 0.0368(0.0008) | 0.2995 (0.0005)
RIDGE | 0.8813(0.0083) | 0.3629(0.0008) | 0.3542(0.0041) | 0.3157(0.0006) | 0.1847(0.0022) | 0.3056 (0.0005)
SIGCOMP | 0.1465(0.0026) | 0.3192(0.0006) | 0.0532(0.0006) |  0.3026(0.0005) | 0.0358(0.0004) | 0.2999 (0.0005)
AdaSS |  0.0856(0.0023) |  0.3171(0.0007) |  0.0359(0.0010) | 0.3027(0.0005) |  0.0217(0.0007) |  0.2994 (0.0005)
AdaSStrue | 0.0726(0.0176) | 0.3080(0.0007) | 0.0399(0.0153) | 0.2994(0.0005) | 0.0188(0.0048) | 0.2977 (0.0005)
| Dampened harmonic
TRU | 0.2851 (0.0050) | 0.5403 (0.0014) | 0.0983 (0.0010) | 0.5051 (0.0010) | 0.0651 (0.0009) | 0.4960 (0.0010)
SMOOTH | 0.2288 (0.0042) | 0.5391 (0.0013) | 0.0836 (0.0007) | 0.5032 (0.0010) | 0.0555 (0.0005) | 0.4936 (0.0010)
PCA | 0.3710 (0.0093) | 0.5259 (0.0012) | 0.1100 (0.0020) |  0.4994 (0.0010) | 0.0594 (0.0011) |  0.4915 (0.0010)
RIDGE | 1.4221 (0.0135) | 0.5925 (0.0016) | 0.6082 (0.0076) | 0.5203 (0.0011) | 0.3271 (0.0038) | 0.5014 (0.0010)
SIGCOMP | 0.2541 (0.0045) |  0.5221 (0.0012) | 0.1235 (0.0013) | 0.5018 (0.0010) | 0.0942 (0.0009) | 0.4950 (0.0010)
AdaSS |  0.1749 (0.0038) | 0.5241 (0.0012) |  0.0695 (0.0012) | 0.4997 (0.0010) |  0.0461 (0.0008) | 0.4918 (0.0010)
AdaSStrue | 0.1504 (0.0030) | 0.5179 (0.0012) | 0.0744 (0.0018) | 0.4985 (0.0010) | 0.0582 (0.0022) | 0.4912 (0.0010)
| Rapid change
TRU | 1.9910(0.0278) | 4.0461(0.0001) | 0.9178(0.0100) | 3.7583(0.0001) | 0.6020(0.0074) | 3.6989 (0.0001)
SMOOTH | 1.2961(0.0133) | 3.9427(0.0001) | 0.5738(0.0046) | 3.7205(0.0001) | 0.3590(0.0027) | 3.6787 (0.0001)
PCA | 5.1052(0.0971) | 4.3070(0.0001) | 1.5870(0.0271) | 3.7978(0.0001) | 0.8383(0.0125) | 3.7141 (0.0001)
RIDGE | 10.4781(0.1059) | 4.4295(0.0001) | 4.1991(0.0537) | 3.8459(0.0001) | 2.2250(0.0278) | 3.7356 (0.0001)
SIGCOMP | 1.7129(0.0209) | 4.0352(0.0001) | 0.8615(0.0234) | 3.7702(0.0001) | 0.8552(0.0167) | 3.7428 (0.0001)
AdaSS |  1.0482(0.0166) |  3.8737(0.0001) |  0.4526(0.0077) |  3.6928(0.0001) |  0.2916(0.0044) |  3.6662 (0.0001)
AdaSStrue | 0.8181(0.0191) | 3.8274(0.0001) | 0.3434(0.0080) | 3.6759(0.0001) | 0.2114(0.0050) | 3.6541 (0.0001)

###  3.2 Dampened Harmonic Motion Function

This simulation scenario considers as coefficient function $\beta$ the
dampened harmonic motion function, also known as the spring function in the
engineering literature. It is characterized by a sinusoidal behaviour with
exponentially decreasing amplitude, that is

$\beta\left(s,t\right)=1+5\exp\left[-5\left(s+t\right)\right]\left[\cos\left(10\pi
s\right)+\cos\left(10\pi t\right)\right],\quad
s,t\in\left[0,1\right]\times\left[0,1\right].$

Figure 3 displays the AdaSS and the SMOOTH estimates along with the true
coefficient function. Also in this scenario, the AdaSS estimates is smoother
than the SMOOTH estimates in regions of small curvature. But, it is more
flexible where the coefficient function is more wiggly. Note that intuitively,
the SMOOTH estimator trades off its smoothness over the whole domain. Indeed,
it over-smooths at small values of $s$ and $t$ and under-smooths elsewhere.

Figure 3:  AdaSS (solid line) and SMOOTH (dashed line) estimates of the
coefficient functions and the TRUE coefficient function $\beta$ (dotted line)
for different values of $t$ in the case of the dampened harmonic motion
function.

In Figure 4 and in the second tier of Table 1, values of the ISE and PMSE for
the AdaSS, AdaSStrue, and competitor estimators are shown as function of the
sample size $n$, in the case of the dampened harmonic motion function. Even in
this case, the AdaSS estimator achives the lowest ISE for all sample sizes,
and thus, the lowest estimation error, without taking into account the
AdaSStrue estimator. Strictly speaking, in terms of PMSE, note that the
proposed estimator is not always the best choice, but it shows only a small
difference with best methods, viz., PCA and SIGCOMP estimators. In this case,
the AdaSS and AdaSStrue performance is very similar for $n=500,1000$, whereas,
for $n=150$, the AdaSStrue performs slightly better especially in terms of
PMSE.

(a)

(b)

Figure 4:  4(a) The integrated squared error (ISE) and 2(b) the prediction
mean squared error (PMSE) $\pm standard\hskip 3.41418pterror$ for the TRU,
SMOOTH, PCA, SIGCOMP, AdaSS and AdaSStrue estimators in the case of the
dampened harmonic motion function. The Ridge estimator is not considered due
to its too different performance.

###  3.3 Rapid Change Function

In this scenario the true coefficient function $\beta$ is obtained by the
rapid change function, that is

$\beta\left(s,t\right)=1-\frac{5}{1+\exp\left[10\left(s+t-0.2\right)\right]}+\frac{5}{1+\exp\left[75\left(s+t-0.8\right)\right]},\quad
s,t\in\left[0,1\right]\times\left[0,1\right].$

Figure 5 shows the AdaSS and SMOOTH estimate when $\beta$ is the rapid change
function. The SMOOTH estimate is rougher than the AdaSS one in regions that
are far from the rapid change point. On the contrary, the AdaSS estimate is
able to be smoother in the flat region and to be as accurate as the SMOOTH
estimate near the rapid change point.

Figure 5:  AdaSS (solid line) and SMOOTH (dashed line) estimates of the
coefficient functions and the TRUE coefficient function $\beta$ (dotted line)
for different values of $t$ in the case of the rapid change function.

In Figure 6 and the third tier of Table 1, values of the ISE and PMSE for the
AdaSS, AdaSStrue, and competitor estimators are shown for sample sizes
$n=150,500,1000$. In this case, the AdaSS estimator outperforms the
competitors, both in terms of ISE and PMSE. Also in this case, the performance
of the AdaSStrue estimator is slightly better than that of the AdaSS one and
this difference in performance reduces as $n$ increases.

(a)

(b)

Figure 6:  4(a) The integrated squared error (ISE) and 2(b) the prediction
mean squared error (PMSE) $\pm standard\hskip 3.41418pterror$ for the TRU,
SMOOTH, PCA, SIGCOMP, AdaSS and AdaSStrue estimators in the case of the rapid
change function function. The Ridge estimator is not considered due to its too
different performance.

##  4 Real-data Examples

In this section, two real datasets, namely Swedish mortality and ship CO2
emission datasets, are considered in order to asses the performance of the
AdaSS estimator in real applications.

###  4.1 Swedish Mortality Dataset

The Swedish mortality dataset (available from the Human Mortality Database
—http://mortality.org—) is very well known in the functional literature as
benchmark dataset. It has been analysed by Chiou and Müller [10] and Ramsay et
al. [31], among others. In this analysis, we consider the log-hazard rate
functions of the Swedish females mortality data for year-of-birth cohorts that
refer to females born in the years 1751-1935 with ages 0-80. The value of a
log-hazard rate function at a given age is the natural logarithm of the ratio
of females died at that age and the number of females alive with the same age.
The 184 considered log-hazard rate functions [10] are shown in Figure 7.
Without loss of generality they have been normalized to the domain
$\left[0,1\right]$.

Figure 7: Log-hazard rate functions for Swedish female cohorts from 1751 to
1935.

The functions from 1751 (1752) to 1934 (1935) are considered as observations
$X_{i}$ ($Y_{i}$) of the predictor (response) in (1.1), $i=1,\dots,184$. In
this way, the relationship between two consecutive log-hazard rate functions
becomes the focus of the analysis. To asses the predictive performance of the
methods considered in the simulation study (Section 3), for 100 times, 166
observations out of 184 are randomly chosen, as training set, to fit the
model. The 18 remaining ones are used as test set to calculate the PMSE. The
averages and standard deviations of PMSEs are shown in the first line of Table
2. The AdaSS estimator outperforms all the competitors. Only the RIDGE
estimator has comparable predictive performance.

Table 2: The prediction mean squared error (PMSE) for the TRU, SMOOTH, PCA,
RIDGE, SIGCOMP, and AdaSS estimators. The numbers outside the parentheses are
the averages of the PMSE over $100$ replications, and the numbers inside
parentheses are the corresponding standard errors.

| TRU | SMOOTH | PCA | RIDGE | SIGCOMP | AdaSS
---|---|---|---|---|---|---
Swedish mortality ($\times 10^{-2}$) | 0.7373 (0.0000) | 0.5938 (0.0000) | 0.6131 (0.0000) | 0.5749 (0.0000) | 1.0173 (0.0000) |  0.5706 (0.0000)
Ship CO2 emission | 0.1019 (0.0008) | 0.0814 (0.0007) | 0.0689 (0.0008) |  0.0625 (0.0007) | 0.1033 (0.0013) | 0.0771 (0.0007)

Figure 8 shows the AdaSS estimates along with the RIDGE estimates that
represents the best competitor methods in terms of PMSE. The proposed
estimator has slightly better performance than the competitor, but, at the
same time, it is much more interpretable. In fact, it is much smoother where
the coefficient function seem to be mostly flat and successfully captures the
pattern of $\beta$ in the peak region. On the contrary, the RIDGE estimates is
particularly rough over region of low curvature.

Figure 8:  AdaSS (solid line) and RIDGE (dashed line) estimates of the
coefficient functions for different values of $t$ in the Swedish Mortality
dataset.

###  4.2 Ship CO2 Emission Dataset

The ship CO2 emission dataset has been thoroughly studied in the very last
years [22, 32, 7, 9]. It was provided by the shipping company Grimaldi Group
to address some aspects that are related to the issue of monitoring fuel
consumptions or CO2 emissions for Ro-Pax ship that sail along a route in the
Mediterranean Sea. In particular, we focus on the study of the relation
between the fuel consumption per hour (FCPH), assumed as the response, and the
speed over ground (SOG), assumed as predictor. The observations considered
were recorded from 2015 to 2017. Figure 9 shows the 44 available observations
of SOG and FCPH [9].

Figure 9: SOG and FCPH observations from a Ro-Pax ship.

Similarly to the Swedish mortality dataset, The prediction performance of the
methods are assessed by randomly chosen 40 out of 44 observations to fit the
model and by using the 4 remaining observations to compute the PMSE. This is
repeated 100 times. The averages and standard deviations of the PMSEs are
listed in the second line of Table 2. The AdaSS estimator is in this case
outperformed by the RIDGE estimator, which achieves the lowest PMSE. However,
as shown in Figure 10, it is able both to well estimate the coefficient
function over peaky regions, as the RIDGE estimator, and to smoothly adapt
over the remaining part of the domain. In this case, also the PCA estimator
achieves smaller PMSE than that of the proposed estimator. However, the PCA
estimator is even rougher than the RIDGE estimator and, thus, it is not shown
in Figure 10.

Figure 10:  AdaSS (solid line) and RIDGE (dashed line) estimates of the
coefficient functions for different values of $t$ in the ship CO2 emission
dataset.

##  5 Conclusion

In this article, the AdaSS estimator is proposed for the function-on-function
linear regression model where each value of the response, for any domain
point, depends linearly on the full trajectory of the predictor. The
introduction of two adaptive smoothing penalties, based on initial estimate of
its partial derivatives, allows the proposed estimator to better adapt to the
coefficient function. By means of a simulation study, the proposed estimator
has proven favourable performance with respect to those achieved by the five
competitors already appeared in the literature before, both in terms of
estimation and prediction error. The adaptive feature of the AdaSS estimator
is advantageous for the interpretability of the results with respect to the
competitors. Moreover, its performance has shown to be competitive also with
respect to the case where the true coefficient function is known. Finally, the
proposed estimator has been successfully applied to real-data examples, viz.,
the Swedish mortality and ship CO2 emission datasets. However, some challenges
are still open. Even though the proposed evolutionary algorithm has shown to
perform particularly well both in the simulation study and the real-data
examples, the choice of the tuning parameters still remains in fact a critical
issue, because of the curse of dimensionality. This could be even more
problematic in the perspective to extend the AdaSS estimator to the FoF
regression model with multiple predictors.

##  Appendix A Approximation of the Two Penalty Terms for the AdaSS Estimator
Derivation

In this section the approximations of (2.6) and (2.7) are obtained. For the
first penalty, by using (2.4), we have

$\lambda^{AdaSS}_{s}\int_{\mathcal{S}}\int_{\mathcal{T}}\frac{1}{\left(|\widehat{\beta_{s}^{m_{s}}}\left(s,t\right)|+\delta_{s}\right)^{\gamma_{s}}}\left(\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)\right)^{2}dsdt\\\
\approx\lambda^{AdaSS}_{s}\int_{\mathcal{S}}\int_{\mathcal{T}}\frac{1}{\left(|\sum_{i=1}^{L_{s}+1}\sum_{j=1}^{L_{t}+1}\widehat{\beta_{s}^{m_{s}}}\left(\tau_{s,i},\tau_{t,j}\right)I_{\left[\left(\tau_{s,i-1},\tau_{s,i}\right)\times\left(\tau_{t,j-1},\tau_{t,j}\right)\right]}\left(s,t\right)|+\delta_{s}\right)^{\gamma_{s}}}\left(\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)\right)^{2}dsdt\\\
=\lambda^{AdaSS}_{s}\sum_{i=1}^{L_{s}+1}\sum_{j=1}^{L_{t}+1}\int_{\left[\tau_{s,i-1},\tau_{s,i}\right]}\int_{\left[\tau_{t,j-1},\tau_{t,j}\right]}\frac{1}{\left(|\widehat{\beta_{s}^{m_{s}}}\left(\tau_{s,i},\tau_{t,j}\right)|+\delta_{s}\right)^{\gamma_{s}}}\left(\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)\right)^{2}dsdt\\\
=\lambda^{AdaSS}_{s}\sum_{i=1}^{L_{s}+1}\sum_{j=1}^{L_{t}+1}d^{s}_{ij}\int_{\left[\tau_{s,i-1},\tau_{s,i}\right]}\int_{\left[\tau_{t,j-1},\tau_{t,j}\right]}\left(\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)\right)^{2}dsdt,$
(A.1)

where
$d^{s}_{ij}=\Big{\\{}\frac{1}{\left(|\widehat{\beta_{s}^{m_{s}}}\left(\tau_{s,i},\tau_{t,j}\right)|+\delta_{s}\right)^{\gamma_{s}}}\Big{\\}}$.
Then, for (2.2), and following Ramsay and Silverman [30], pag. 292,

$\lambda^{AdaSS}_{s}\sum_{i=1}^{L_{s}+1}\sum_{j=1}^{L_{t}+1}d^{s}_{ij}\int_{\left[\tau_{s,i-1},\tau_{s,i}\right]}\int_{\left[\tau_{t,j-1},\tau_{t,j}\right]}\left(\mathcal{L}_{s}^{m_{s}}\alpha\left(s,t\right)\right)^{2}dsdt\\\
=\lambda^{AdaSS}_{t}\sum_{i=1}^{L_{s}+1}\sum_{j=1}^{L_{t}+1}d^{t}_{ij}\operatorname{Tr}\left[\bm{B}_{\alpha}^{T}\bm{W}_{s,i}\bm{B}_{\alpha}\bm{R}_{t,j}\right]$
(A.2)

where
$\bm{W}_{s,i}=\int_{\left[\tau_{s,i-1},\tau_{s,i}\right]}\bm{\psi}^{s}\left(s\right)\bm{\psi}^{s}\left(s\right)^{T}ds$,
$\bm{R}_{t,j}=\int_{\left[\tau_{t,j-1},\tau_{t,j}\right]}\mathcal{L}_{t}^{m_{t}}\left[\bm{\psi}^{t}\left(t\right)\right]\mathcal{L}_{t}^{m_{t}}\left[\bm{\psi}^{t}\left(t\right)\right]^{T}dt$,
and
$d^{s}_{ij}=\Big{\\{}\frac{1}{\left(|\widehat{\beta_{s}^{m_{s}}}\left(\tau_{s,i},\tau_{t,j}\right)|+\delta_{s}\right)^{\gamma_{s}}}\Big{\\}}$,
for $i=1,\dots,L_{s}+1$ and $j=1,\dots,L_{t}+1$. Thus, (2.6) is demonstrated,
the arguments are analogous for (2.7).

## References

  * [1] Abramovich, F. and Steinberg, D. M. (1996).  Improved inference in nonparametric regression using lk-smoothing splines.  _Journal of Statistical Planning and Inference_ 49, 3, 327–341.
  * [2] Bäck, T., Fogel, D. B., and Michalewicz, Z. (1997).  _Handbook of evolutionary computation_.  CRC Press.
  * [3] Bergstra, J. and Bengio, Y. (2012).  Random search for hyper-parameter optimization.  _Journal of Machine Learning Research_ 13, Feb, 281–305.
  * [4] Bergstra, J. S., Bardenet, R., Bengio, Y., and Kégl, B. (2011).  Algorithms for hyper-parameter optimization.  In _Advances in neural information processing systems_. 2546–2554.
  * [5] Boyd, S., Boyd, S. P., and Vandenberghe, L. (2004).  _Convex optimization_.  Cambridge university press.
  * [6] Canale, A. and Vantini, S. (2016).  Constrained functional time series: Applications to the italian gas market.  _International Journal of Forecasting_ 32, 4, 1340–1351.
  * [7] Capezza, C., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2020).  Control charts for monitoring ship operating conditions and co2 emissions based on scalar-on-function regression.  _Applied Stochastic Models in Business and Industry_.
  * [8] Cardot, H., Ferraty, F., and Sarda, P. (2003).  Spline estimators for the functional linear model.  _Statistica Sinica_ , 571–591.
  * [9] Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2020).  Functional regression control chart.  _Technometrics_ , 1–14.
  * [10] Chiou, J.-M. and Müller, H.-G. (2009).  Modeling hazard rates as functional data for the analysis of cohort lifetables and mortality forecasting.  _Journal of the American Statistical Association_ 104, 486, 572–585.
  * [11] Eubank, R. L. (1999).  _Nonparametric regression and spline smoothing_.  CRC press.
  * [12] Ferraty, F. and Vieu, P. (2006).  _Nonparametric functional data analysis: theory and practice_.  Springer Science & Business Media.
  * [13] Green, P. J. and Silverman, B. W. (1993).  _Nonparametric regression and generalized linear models: a roughness penalty approach_.  Chapman and Hall/CRC.
  * [14] Gu, C. (2013).  _Smoothing spline ANOVA models_. Vol. 297.  Springer Science & Business Media.
  * [15] Hastie, T., Tibshirani, R., and Friedman, J. (2009).  _The elements of statistical learning: data mining, inference, and prediction_.  Springer series in statistics New York, NY, USA:.
  * [16] Herlihy, M. and Shavit, N. (2011).  _The art of multiprocessor programming_.  Morgan Kaufmann.
  * [17] Horváth, L. and Kokoszka, P. (2012).  _Inference for functional data with applications_. Vol. 200.  Springer Science & Business Media.
  * [18] Hsing, T. and Eubank, R. (2015).  _Theoretical foundations of functional data analysis, with an introduction to linear operators_.  John Wiley & Sons.
  * [19] Jaderberg, M., Dalibard, V., Osindero, S., Czarnecki, W. M., Donahue, J., Razavi, A., Vinyals, O., Green, T., Dunning, I., Simonyan, K., and others. (2017).  Population based training of neural networks.  _arXiv preprint arXiv:1711.09846_.
  * [20] James, G. M. (2002).  Generalized linear models with functional predictors.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 64, 3, 411–432.
  * [21] Kokoszka, P. and Reimherr, M. (2017).  _Introduction to functional data analysis_.  CRC Press.
  * [22] Lepore, A., Palumbo, B., and Capezza, C. (2018).  Analysis of profiles for monitoring of modern ship performance via partial least squares methods.  _Quality and Reliability Engineering International_ 34, 7, 1424–1436.
  * [23] Luo, R. and Qi, X. (2017).  Function-on-function linear regression by signal compression.  _Journal of the American Statistical Association_ 112, 518, 690–705.
  * [24] Luo, R. and Qi, X. (2019).  Interaction model and model selection for function-on-function regression.  _Journal of Computational and Graphical Statistics_ 28, 2, 1–14.
  * [25] Morris, J. S. (2015).  Functional regression.  _Annual Review of Statistics and Its Application_  _2_ , 321–359.
  * [26] Müller, H.-G., Stadtmüller, U., and others. (2005).  Generalized functional linear models.  _the Annals of Statistics_ 33, 2, 774–805.
  * [27] Pintore, A., Speckman, P., and Holmes, C. C. (2006).  Spatially adaptive smoothing splines.  _Biometrika_ 93, 1, 113–125.
  * [28] Qi, X. and Luo, R. (2018).  Function-on-function regression with thousands of predictive curves.  _Journal of Multivariate Analysis_  _163_ , 51–66.
  * [29] Qi, X. and Luo, R. (2019).  Nonlinear function on function additive model with multiple predictor curves.  _Statistica Sinica_  _29_ , 719–739.
  * [30] Ramsay, J. and Silverman, B. (2005).  _Functional Data Analysis_.  Springer Series in Statistics. Springer.
  * [31] Ramsay, J. O., Hooker, G., and Graves, S. (2009).  _Functional data analysis with R and MATLAB_.  Springer Science & Business Media.
  * [32] Reis, M. S., Rendall, R., Palumbo, B., Lepore, A., and Capezza, C. (2019).  Predicting ships’ co2 emissions using feature-oriented methods.  _Applied Stochastic Models in Business and Industry_.
  * [33] Ruppert, D. and Carroll, R. J. (2000).  Theory & methods: Spatially-adaptive penalties for spline fitting.  _Australian & New Zealand Journal of Statistics_ 42, 2, 205–223.
  * [34] Schumaker, L. (2007).  _Spline functions: basic theory_.  Cambridge University Press.
  * [35] Storlie, C. B., Bondell, H. D., and Reich, B. J. (2010).  A locally adaptive penalty for estimation of functions with varying roughness.  _Journal of Computational and Graphical Statistics_ 19, 3, 569–589.
  * [36] Wahba, G. (1990).  _Spline models for observational data_. Vol. 59.  Siam.
  * [37] Wang, X., Du, P., and Shen, J. (2013).  Smoothing splines with varying smoothing parameter.  _Biometrika_ 100, 4, 955–970.
  * [38] Yang, L. and Hong, Y. (2017).  Adaptive penalized splines for data smoothing.  _Computational Statistics & Data Analysis_  _108_ , 70–83.
  * [39] Yao, F. and Müller, H.-G. (2010).  Functional quadratic regression.  _Biometrika_ 97, 1, 49–64.
  * [40] Yao, F., Müller, H.-G., and Wang, J.-L. (2005a).  Functional data analysis for sparse longitudinal data.  _Journal of the American Statistical Association_ 100, 470, 577–590.
  * [41] Yao, F., Müller, H.-G., and Wang, J.-L. (2005b).  Functional linear regression analysis for longitudinal data.  _The Annals of Statistics_ , 2873–2903.