Smooth Lasso Estimator for the Function-on-Function Linear Regression Model

A new estimator, named S-LASSO, is proposed for the coefficient function of the Function-on-Function linear regression model. The S-LASSO estimator is shown to be able to increase the interpretability of the model, by better locating regions where the coefficient function is zero, and to smoothly estimate non-zero values of the coefficient function. The sparsity of the estimator is ensured by a \textit{functional LASSO penalty}, which pointwise shrinks toward zero the coefficient function, while the smoothness is provided by two roughness penalties that penalize the curvature of the final estimator. The resulting estimator is proved to be estimation and pointwise sign consistent. Via an extensive Monte Carlo simulation study, the estimation and predictive performance of the S-LASSO estimator are shown to be better than (or at worst comparable with) competing estimators already presented in the literature before. Practical advantages of the S-LASSO estimator are illustrated through the analysis of the \textit{Canadian weather}, \textit{Swedish mortality} and \textit{ship CO\textsubscript{2} emission data}. The S-LASSO method is implemented in the \textsf{R} package \textsf{slasso}, openly available online on CRAN.


Introduction
Functional linear regression (FLR) is the generalization of the classical multivariate regression to the context of the functional data analysis (FDA) (e.g. Ramsay and Silverman (2005); Horváth and Kokoszka (2012); Hsing and Eubank (2015); Kokoszka and Reimherr (2017)), where either the predictor or the response or both have a functional form. In particular, we study the Functionon-Function (FoF) linear regression model, where both the predictor and the response variable are functions and each value of the latter, for any domain point, depends on the full trajectory of the former. The model is as follows for i = 1, . . . , n. The pairs (X i , Y i ) are independent realizations of the predictor X and the response Y , which are assumed to be smooth random process with realizations in L 2 (S) and L 2 (T ), i.e., the Hilbert spaces of square integrable functions defined on the compact sets S and T , respectively. Without loss of generality, the latter are also assumed with functional mean equal to zero.
The functions ε i are zero-mean random errors, independent of X i , and have autocovariance function K (t 1 , t 2 ), t 1 and t 2 ∈ T . The function β is smooth in L 2 (S ×T ), i.e., the Hilbert space of bivariate square integrable functions defined on the compact set S × T , and is hereinafter referred to as coefficient function.
FLR analysis is a hot topic in the FDA literature. A comprehensive review of the main results is provided by Morris (2015) as well as by Ramsay and Silverman (2005); Horváth and Kokoszka (2012) and Cuevas (2014) who give worthwhile modern perspectives. Although the research efforts have been focused mainly on the case where either the predictor or the response have functional form (Cardot et al., 2003;Li et al., 2007;Hall et al., 2007;Abramowicz et al., 2018), the interest in the FoF linear regression has increased in the very last years. In particular, Besse and Cardot (1996) developed a spline based approach to estimate the coefficient function β, while Ramsay and Silverman (2005) proposed an estimator assumed to be in a finite dimension tensor space spanned by two basis sets and where regularization is achieved by either truncation or roughness penalties. Yao et al. (2005b) built up an estimation method based on the principal component decomposition of the autovariance function of both the predictor and the response based on the principal analysis by conditional expectation (PACE) method (Yao et al., 2005a). This estimator was extended by Chiou et al. (2014) to the case of multivariate functional responses and predictors. A general framework for the estimation of the coefficient function was proposed by Ivanescu et al. (2015) by means of the mixed model representation of the penalized regression. An extension of the ridge regression (Hastie et al., 2009) to the FoF linear regression with an application to the Italian gas market was presented in Canale and Vantini (2016). To take into account the case when the errors ε i are correlated, in Scheipl et al. (2015) the authors developed a general framework for additive mixed models by extending the work of Ivanescu et al. (2015).
Analogously to the classical multivariate setting, in Equation (1) the functional predictor X contributes linearly to the response Y through the coefficient function β, which works as a continuous weight function. Trivially note that, in the domain regions over which β is equal to zero (if any), changes in the functional predictor X do not affect the conditional value of Y . Because of the infinite dimensional nature of the FLR problem, coefficient functions that are sparse, i.e., zero valued over large parts of domain, arise very often in real applications. When the aim of the analysis is descriptive, that is the interest relies on understanding the relationship between X and Y , rather than predictive only, methods that are able to capture the sparse nature of the coefficient function may be of great practical interest. These methods are referred to as sparse or interpretable because they allow better interpreting the effects of the predictor on the response and reveal the sparse nature of β. On the contrary, the interpretation of the relationship between X and Y is often cumbersome for methods that do not focus on the sparsity of the coefficient function. In particular, here the interpretability of the model in Equation (1) is ultimately related to the knowledge of the parts of the domain S × T where β is equal or different to zero, which are hereinafter referred to as null and non-null regions, respectively.
Few works address the issue of the interpretability in FLR. In the scalar-onfunction setting,  proposed the FLiRTI (Functional Linear Regression That's Interpretable) estimator that is able to recover the sparseness of the coefficient function, by imposing L 1 -penalties on the coefficient function itself and its first two derivatives. Zhou et al. (2013) introduced an estimator obtained in two stages where the initial estimate is obtained by means of a Dantzig selector (Candes et al., 2007)

refined via the group Smoothly Clipped
Absolute Deviation (SCAD) penalty (Fan and Li, 2001). The most recent work that addresses the issue of interpretability is that of Lin et al. (2017), who pro- However, to the best of the author knowledge, no effort has been made in the literature to obtain an interpretable estimator for the FoF linear regression model. In this work, we try to fill this gap by developing an interpretable estimator of the coefficient function β, named S-LASSO (Smooth plus LASSO) that is locally sparse (i.e., is zero on the null region) and, at the same time, smooth on the non-null region. The property of sparseness of the S-LASSO estimator is provided by a functional LASSO penalty (FLP), which is the functional generalization of the classical Least Absolute Shrinkage and Selection Operator (LASSO) (Tibshirani, 1996). Whereas, two roughness penalties, introduced in the objective function, ensure smoothness of the estimator on the non-null region. From a computational point of view, the S-LASSO estimator is obtained as the solution of a single optimization problem by means of a new version of the orthant-wise limited-memory quasi-Newton (OWL-QN) algorithm (Andrew and Gao, 2007), which is specifically designed to solve optimization problems involving L 1 penalties. The method presented in this article is implemented in the R package slasso, openly available online on CRAN.
To give an idea of the properties of the proposed estimator, in Figure 1 the S-LASSO estimator is applied to four different scenarios, whose data generation is detailed in Section 4. In particular, in each plot the S-LASSO estimate, the true coefficient function, and the smoothing spline estimate proposed by Ramsay and Silverman (2005), referred to as SMOOTH, are shown for t = 0.5.
For Scenario I, the true coefficient function is zero all over the domain, which means that the predictor X is independent of the response. It is clear from while simultaneously estimating the relationship between X and Y over the nonnull region. The performance of the S-LASSO estimator will be deeply analysed in Section 4.
The paper is organized as follows. In Section 2, the S-LASSO estimator is presented. In Section 3, asymptotic properties of the S-LASSO estimator are discussed in terms of consistency and pointwise sign consistency. In Section 4, by means of a Monte Carlo simulation study, the S-LASSO estimator is compared, in terms of estimation error and prediction accuracy, with competing estimators already proposed in the literature. In Section 5, the potential of the S-LASSO estimator are demonstrated with respect of three datasets: the Canadian weather, Swedish mortality and ship CO 2 emission data. Proofs, data generation procedures in the simulation study as well as additional simulation results, and, algorithm description are given in the Supplementary Material.

Methodology
In Section 2.1, we briefly describe the smoothing spline estimator. Readers who are already familiar with this approach may skip to the next subsection.
In Section 2.2, the S-LASSO estimator definition is given along with details on both computational issues and model selection.

The smoothing spline estimator
The smoothing spline estimator of the FoF linear regression model (Ramsay and Silverman, 2005) is the first key component of the S-LASSO estimator.
It is based on the assumption that the coefficient function β may be well approximated by an element in the tensor product space generated by two spline function spaces, where a spline is a function defined piecewise by polynomials.
Well-known basis functions for the spline space are the B-splines. A B-spline basis is a set of spline functions uniquely defined by an order k and a non-decreasing sequence of M + 2 knots, that we hereby assume to be equally spaced in a general domain D. Cubic B-splines are B-splines of order k = 4. Each B-spline function is a positive polynomial of degree k − 1 over each subinterval defined by the knot sequence and is non-zero over no more than k of these subintervals (i.e., the compact support property). In our setting, besides the computational advantage (Hastie et al., 2009), the compact support property is fundamental because it allows one to link the values of β over a given region to the B-splines with support in the same region and to discard all the B-splines that are outside that region. Thorough descriptions of splines and B-splines are in De Boor (2001) and Schumaker (2007).
The smoothing spline estimator (Ramsay and Silverman, 2005) is defined aŝ where S k1,k2,M1,M2 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 S and T , respectively. L ms s and L mt t , with m s ≤ k 1 − 1 and m t ≤ k 2 − 1, are the m s -th and m t -th order linear differential operators applied to α with respect to the variables s and t, respectively. The symbol || · || denotes the L 2 -norm corresponding to the inner product < f, g >= f g.
The parameters λ s , λ t ≥ 0 are generally referred to as roughness parameters.
The aim of the second and third terms on the right-hand side of Equation (2) is that of penalizing features along s and t directions. A common practice, when dealing with cubic splines, is to choose m s = 2 and m t = 2, which results in the penalization of the curvature of the final estimator. When λ s = λ t = 0, the wiggliness of the estimator is not penalized and the resulting estimator is the one that minimizes the sum of squared errors. On the contrary, as λ s → ∞ and λ t → ∞,β S converges to a bivariate polynomial with degree equal to | max (m s , m t )−1|. However, there is no guarantee thatβ S is a sparse estimator, i.e., it is exactly equal to zero in some part of the domain S × T .

The S-LASSO Estimator
Based on the smoothing spline estimator of Equation (2), the S-LASSO estimator is defined as followŝ The last term in the right-hand side of Equation (3) is the extension of the LASSO penalty (Tibshirani, 1996) to the FoF linear regression setting, which is referred to as functional LASSO penalty (FLP). In particular, by starting from the multivariate LASSO penalty, the FLP is built by substituting summation with integration and by taking into account the continuity of the domain. The FLP is able to pointwise shrink the value of the coefficient function for each s, t ∈ S × T . Due to the property of the absolute value function of being singular at zero, some of these values are shrunken exactly to zero. Thus, the FLP allowsβ SL to be exactly zero over the null region. The constant λ L ≥ 0 is usually called regularization parameter and controls the degree of shrinkage towards zero of the FLP. The larger this value, the larger the shrinkage effect and the domain portion where the resulting estimator is zero. Moreover, the FLP is expected to be able to improve the prediction accuracy (in terms of expected mean square error) by introducing a bias in the final estimator.
The second and third terms on the right-hand side of Equation (3) represent the two roughness penalties introduced in Equation (2) to control the smoothness of the coefficient function estimator.
It is worth noting that, in general, the estimator smoothness may be also controlled by opportunely choosing the dimension of the space S k1,k2,M1,M2 , that is, by fixing k 1 and k 2 , and choosing M 1 and M 2 (Ramsay and Silverman, 2005).
However, this strategy is not suitable in this case. To obtain a sparse estimator, the dimension of the space S k1,k2,M1,M2 must be in fact as large as possible. In this way, the value of β in a given region is strictly related to the coefficients of the B-spline functions defined on the same part of the domain and, thus, they tend to be zero in the null region. On the contrary, when the dimension of S k1,k2,M1,M2 is small, there is a larger probability that some B-spline functions have support both in the null and non-null regions and, thus the corresponding B-spline coefficients result different from zero. Therefore, we find suitable the use of the two roughness penalty terms in Equation (3).
To compute the S-LASSO estimator, let us consider the space S k1,k2,M1,M2  (3) over α ∈ S k1,k2,M1,M2 , we implicitly assume that β can be suitably approximated by an element in S k1,k2,M1,M2 , that is where B = {b ij } ∈ R M1+k1×M2+k2 and b ij are scalar coefficients. So stated, the problem of estimating β has been reduced to the estimation of the unknown co- Then, the first term of the right-hand side of Equation (3) may be rewritten as whereas, the roughness penalties on the left side of Equation (3) become

The term Tr [A] denotes the trace of a square matrix A.
Standard optimization algorithms for L 1 -regularized objective functions are designed for the case where the absolute value enters the problem as a linear function of the parameters. However, in the optimization problem in Equation (3), the FLP particularizes as follows which is not a linear function of the absolute value of the coefficient matrix |B α |, because the absolute value is instead applied to a linear combination of the parameters. Therefore, in order to be able to use optimization algorithms for L 1 -regularized objective functions by the following theorem, we provide a practical way to approximate the FLP as a linear function of |B α |, and thus extremely simplify the computation.
, with {ψ ij } the set of B-splines of orders k j and non-decreasing evenly spaced knots sequences and ∆ j,e the extended partitions corresponding to ∆ j defined as ∆ j,e = {y j,l } Mj +2kj l=1 where y j,1 , . . . , y j,kj = x j,0 , y j,1+kj , . . . , y j,Mj +kj = x j,1 , . . . , x j,Mj and y j,Mj +1+kj , . . . , y j,Mj +2kj = where and The interpretation of the above theorem is quite simple. It basically says that for large values of M 1 and M 2 , ||f || L 1 is well approximated from the top by ||f || 1 ,∆2,e,∆1,e and the approximation error tends to zero as M 1 , M 2 → ∞.
By using this result, the FLP can be approximated as follows where {s e i } and {t e i } are the extended partitions associated with ∆ s and ∆ t , re- Therefore, upon using the approximation in Equation (11), Equation (5) and Equation (6), the optimization problem in Equation (3) becomeŝ Then, the coefficient β is estimated byβ SL (s, t) = ψ s (s) TB SL ψ t (t) for s ∈ S and t ∈ T . Note that, in Equation (12), the FLP is approximated through a weighted linear combination of the absolute values of the coefficients which strictly resembles the multivariate LASSO penalty applied to the basis expan- M2+k2 j=1 |b α,ij |. However, the presence of w s and w t in the FLP approximation in Equation (12) is crucial because it allows the penalty to differently shrink coefficients among B-splines. That is, it avoids that the absolute values of coefficients corresponding to B-splines strictly localized are weighted as the absolute value of coefficients of spreader basis in the computation of the penalty. This is the direct consequence of the fact that the proposed approximation is a better approximation of the FLP than the multivariate LASSO penalty applied to the coefficients.
Remark 1. Theorem 1 states that the L 1 norm of a bivariate function in S k1,k2,∆1,e∆2,e defined in Equation (10) (44)) that provides guidelines to choose an appropriate knot sequence to construct a set of B-splines for any particular polynomial spline space defined on a univariate closed interval. Specifically, the Curry-Schoenberger theorem requires the knot sequence to be an extended knot partition (Schumaker (2007), Definition 4.8). That is, being k the order of the set of B-splines, the first k and last k knots should be respectively chosen smaller than or equal to lower bound and larger than or equal to the upper bound of the definition interval of the polynomial spline space. A convenient choice (De Boor, 2001;Ramsay and Silverman, 2005) is to set them exactly equal to the lower and upper bounds, respectively. In this way, we avoid to impose continuity condition (De Boor, 2001) at the boundaries where, thus, differentiability of the coefficient function is lost.
This choice is consistent with the fact that the set of B-splines provides a valid representation of the corresponding spline space only in the definition interval, as there are generally no information about the behaviour of the function to be estimated outside the definition interval, where the coefficient function is not constrained to be continuous. This intuition is directly translated to Theorem 1, where no continuity condition is imposed at the definition interval boundaries to the bivariate functions generated by the two sets of B-splines.
The optimization problem with L 1 -regularized loss in Equation (12) is (i) convex, being sum or integral of convex function; and (ii) has a unique solution given some general conditions on the matrix W t ⊗ X T X (with ⊗ the Kronecker product). See Section 3 for further details. Unfortunately, the objective function is not differentiable in zero, and thus it has not a closed-form solution.
In view of this, general purpose gradient-based optimization algorithms -as for instance the L-BFGS quasi-Newton method (Nocedal and Wright, 2006)and classical optimization algorithms for solving LASSO problems -such as coordinate descent (Friedman et al., 2010) and least-angle regression (LARS) (Efron et al., 2004)-are not suitable. In contrast, we found very promising a modified version of the orthant-wise limited-memory quasi-Newton (OWL-QN) algorithm proposed by Andrew and Gao (2007). The OWL-QN algorithm is based on the fact that the L 1 norm is differentiable for the set of points named orthant in which each coordinate never changes sign, being a linear function of its argument. In each orthant, the second-order behaviour of an objective function of the form f (x) = l (x) + C||x|| 1 , to be minimized, is determined by given positive constant, and || · || 1 is the usual L 1 norm. Therefore, Andrew and Gao (2007) propose to derive a quadratic approximation of the function l that is valid for some orthant containing the current point and then to search for the minimum of the approximation, by constraining the solution in the orthant where the approximation is valid. There may be several orthants containing or adjacent to a given point. The choice of the orthant to explore is based on the pseudo-gradiant f (x) of f at x, whose components are defined as where sign (·) denotes the usual sign function. However, the objective function of the optimization problem in Equation (12) is in the form f * (x) = l (x) + C||Dx|| 1 , with D = {d i } ∈ R r×r a diagonal matrix of positive weights. To take into account these weights, the OWL-QN algorithm must be implemented with a different pseudo-gradiant f * (x) whose components are defined as A more detailed description of the OWL-QN algorithm for objective functions in the form l (x) + C||Dx|| 1 is given in the Supplementary Material A. Note that, the optimization problem in Equation (12) can be rewritten by vectorization aŝ is the diagonal matrix whose diagonal elements are w T s ⊗ w T t . Moreover, for generic a matrix A ∈ R j×k , vec(A) indicates the vector of length jk obtained by writing the matrix A as a vector column-wise. Therefore, the OWL-QN with pseudo-gradient as in Equation (14) can be straightforwardly applied.
In the following, we summarize all the parameters that need to be set to obtain the S-LASSO estimator. The orders k 1 and k 2 should be chosen with respect to the degree of smoothness we want to achieve, and the computational efforts. The larger the values of k 1 and k 2 , the smoother the resulting estimator will be. Following Remark 1, a standard choice are cubic B-splines, with equally spaced knot sequences, i.e., k 1 = k 2 = 4. Moreover, M 1 and M 2 should be as large as possible to ensure that the null region is correctly captured and the approximation in Equation (11) is valid, with respect to the maximum computational efforts. Finally, at given k 1 , k 2 , M 1 , and M 2 , the optimal values of λ s , λ t and λ L can be selected as those that minimize the the estimated prediction error function CV (λ s , λ t , λ L ), i.e., CV (λ s , λ t , λ L ), over a grid of candidate values (Hastie et al., 2009). However, although this choice could be optimal for the prediction performance, it may affect the interpretability of the model. Much more interpretable models, with a slight decrease in predictive performance, may in fact exist. The k-standard error rule, which is a generalization of the one-standard error rule (Hastie et al., 2009), may be a more reasonable choice.
That is, to choose the most sparse model whose error is no more than k standard errors above the error of the best model. In practice, as spareness is controlled by the parameter λ L , we first find the best model in terms of estimated prediction error at given λ L and then, among the selected models, we apply the k-standard error rule. This rule may be particularly useful when CV

Theoretical Properties of the S-LASSO Estimator
In this section, we provide some theoretical results on the S-LASSO estimator, under some regularity assumptions, i.e., the estimation consistency (Theorem 2) and the pointwise sign consistency (Theorem 3) ofβ SL . All proofs are in the Supplementary Material B.
C 2. The null space of K X , i.e., the covariance operator of X i , is empty where K X f (s ) = S E(X(s )X(s))f (s)ds, s ∈ S, and E||X|| 2 2 < ∞. = p + ν such that 3/2 < p ≤ k 1 − 1 and 3/2 < p ≤ k 2 − 1. , where a n = ω (b n ) means an bn → ∞ for n → ∞, C 5. There exist two positive constants b and B such that where Λ min (M ) and Λ max (M ) denote the minimum and maximum eigenvalues of the matrix M .  (1) (Cardot et al., 2003;Prchal and Sarda, 2007;Scheipl and Greven, 2016). When this assumption is not verified, the results of this section are analogously valid by considering the unique β that satisfies Equation (1), and, for each t ∈ T , belongs to the closure of Im(K X ) = {K X f : f ∈ L 2 (S)}, i.e., the image of the covariance operator of X i (Cardot et al., 2003). C.3 ensures that β is sufficiently smooth. C.4 provides information on the growth rate of the number of knots M 1 and M 2 , which are strictly related to the sample size n. C.5 is the anolugus of condition (F) of Fan et al. (2004) and assumes that the matrix W t ⊗ n −1 X T X has reasonably good behaviour, whereas, C.6 provides guidance on the choice of the parameters λ s and λ t .
Theorem 2 shows that with probability tending to one there exists a solution of the optimization problem in Equation (3) that converges toβ, chosen such To prove Theorem 2, in addition to C.1-C.6, the following condition is considered.
The first result is about the convergence rate ofβ SL to β in terms of L ∞ -norm.
Theorem 2. Under assumptions C.1-C.7, there exists a unique solutionβ SL of the optimization problem in Equation (3), such that According to the above theorem, there exists an estimatorβ SL of β that is root-n/M 1 M 2 consistent.
Before stating Theorem 3, let us define with b (1) the vector whose entries are the q non-zero elements of b that are and with b (2) the vector whose entries are the (M 1 + k 1 )(M 2 + k 2 ) − q elements of b that are equal to zero. In what follows, we assume, without loss of generality, that To prove Theorem 3, in addition to C.1-C.6, the following conditions are considered.
C 8. (S-LASSO irrepresentable condition (SL-IC)) There exists λ s , λ t , λ L , and a constant η > 0 such that, element-wise, W t ⊗ n −1 X T X 21 + n −1 λ s L wr,21 + n −1 λ t L rw,21 W t ⊗ n −1 X T X 11 + n −1 λ s L wr,11 + n −1 λ t L rw,11 C 9. The functions ε i (t) in Equation (1) are zero mean Gaussian random processes with autocovariance function K (t 1 , t 2 ), t 1 and t 2 ∈ T , independent of The SL-IC in C.8 is the straightforward generalization to the problem in Equation (3) of the elastic irrepresentable condition described in Jia and Yu (2010).

It is a consequence of the standard Karush−Kuhn−Tucker (KKT) conditions
applied to the optimization problem in Equation (12). C.9 gives some conditions on the relationship of λ s , λ t , and λ L with M 1 , M 2 and n. In the classical setting, an estimator is sign selection consistent if it has the same sign of the true parameter with probability tending to one. Analogously, we say that an estimator of the coefficient function β is pointwise sign consistent if, in each point of the domain, it has the same sign of β with probability tending to one.
The following theorem states that, under opportune assumptions, the S-LASSO estimator is pointwise sign consistent.
Theorem 3. Under assumptions C.1-C.6 and C.8-C.10,β SL is pointwise sign consistent, that is, for all s ∈ S and t ∈ T , as n → ∞.
This theorem is the functional extension of the sign consistency result for the multivariate LASSO estimator (Zou and Zhang, 2009).
Remark 3. Assumption C1 is in fact not in contradiction with the assumption that the null space of K X , i.e., the covariance operator of X i , is empty. While the assumption that the null space of K X is related to the joint variability of X over the domain S, the assumption C1 refers to the magnitude of X, and, thus, broadly speaking, to its expected value. Following Bosq (2000), K X can be decomposed as where (λ j , v j ), j ≥ 1 is a complete sequence of eigenelements of K X such that ∞ j=1 λ j = E||X|| 2 2 < ∞. From assumption C1 we know that E||X|| 2 2 < c 2 almost surely. Thus, assumption C1 provides information ust on the eigenvalues sum and not on their single values. The null space of K X is empty when all the eigenvalues are strictly larger than zero, while assumption C1 is not related to single eigenvalues. For these reasons, these assumptions are simultaneously used in several papers (Cardot et al., 2003;Prchal and Sarda, 2007;Zhou and Chen, 2012).

Simulation Study
In this section, we conduct a Monte Carlo simulation study to explore the performance of the S-LASSO estimator. We consider four scenarios whose corresponding coefficient functions are depicted in Figure 2. Note that the coefficient function for Scenario I is not shown because it is zero all over the domain. In Scenario II and III, β is sparse, indeed, it is zero on the edge and in the central part of the domain, respectively. Scenario IV corresponds to a non-sparse setting, which is not expected to be favourable to the S-LASSO estimator. The independent observations of the covariates X i are generated as X i = For each scenario, we generate 100 datasets composed of a training set with sample size n and a test set T with size N equal to 4000 that are used to estimate the coefficient function and to test its predictive performance. This is repeated for three different sample sizes n = 150, 500, 1000. As in Lin et al. (2017), we consider the integrated squared error (ISE) to asses the quality of the estimator β of the coefficient function β. In particular, the ISE over the null region (ISE 0 ) and the non-null region (ISE 1 ) are defined as where A 0 and A 1 are the measures of the null (N (β)) and non-null (N N (β)) regions, respectively. The ISE 0 and the ISE 1 are indicators of the estimation error ofβ over both the null and the non-null regions. Moreover, predictive performance is measured through the prediction mean squared error (PMSE), defined as whereβ is obtained through the observations in the training set. The observations in the test set are centred by means of the sample mean functions estimated through the training set observations.
As a remark, the coefficient function β is not identifiable in L 2 (S × T ), because the X i is generated as a finite linear combination of basis functions, i.e., not empty null space of K X . Thus, as stated in Section 3, β is identifiable in the closure of Im(K x ). To obtain a meaningful measure of the estimation accuracy, The performance of the estimators in terms of ISE 0 is displayed in Figure   3. It is not surprising that the estimation error of β over N (β) of the S-LASSO estimator is significantly smaller than those of the other estimators, being the capability of recovering sparseness of β its main feature. In Scenario I, the RIDGE estimator is the only one that performs comparably to the S-LASSO estimator. This is in accordance with the multivariate setting where it is well known that, when the response is independent of the covariates, the ridge estimator is able to shrink all the coefficients towards zero. The TRU, SMOOTH, PFFR, and PCA estimators have difficulties to correctly identify N (β) for all sample sizes. Nevertheless, their performance is very poor at n = 150. In Scenario II, the S-LASSO estimator is still the best one to estimate β over N (β). However, in this case, the RIDGE estimator performance is unsatisfactory and is mainly caused by the lack of smoothness control that makes the estimator over-rough, especially at small n. Among the competitor estimators, the SMOOTH one has the best performance, readily followed by the PFFR estimator. In Scenario III, results are similar to those of Scenario II, even if the TRU estimator appears as the best alternative. Both PCA and RIDGE estimators are not able to successfully recover sparseness of β for n = 150. For the former, the cause is the number of observations not sufficient to capture the covariance structure of the data, whereas for the latter, it is due to the excessive roughness of the estimator. For n = 500 and m = 1000, the PFFR estimator performs very poorly. This is due to the incapacity of the REML approach to appropriately select the smoothing parameters.
Results in terms of ISE 1 are summarized in Figure 4. It is worth noting that, in this case, as expected the performance of the S-LASSO estimator is generally worse than that of the SMOOTH estimator. In some cases, it is worse   overshrink the estimators of the non-null coefficients (Fan et al., 2004;James and Radchenko, 2009). By looking at the result for Scenario II and III, we surmise that this phenomenon arises in the FoF linear regression model as well.
Finally, in Scenario IV, where β is always different from zero, the S-LASSO estimator,performs comparably to the SMOOTH (i.e., the S-LASSO estimator with λ L = 0). In this case β is not sparse and, thus, the FLP does not help. Figure 5 shows PMSE averages and corresponding standard errors for all the considered estimators. Since PMSE is strictly related to the ISE 0 and the ISE 1 , results are totally consistent with those of Figure 3 and Figure 4.
In particular, the S-LASSO estimator outperforms all the competitor ones in favorable scenarios (viz., Scenario I, II, and III), being the corresponding PMSE lower than that achieved by the other competing estimators. In these scenarios, although the performance of the S-LASSO estimator in terms of ISE 1 is not excellent, the clear superiority in terms of ISE 0 compensates and gives rise to smaller PMSE. Otherwise, for Scenario IV, where the coefficient function is not sparse, the performance of the S-LASSO estimator is very similar to that of the SMOOTH estimator, which is the best one in this case. This is encouraging, because, it proves that the performance of the S-LASSO estimator does not dramatically decline in less favourable scenarios.
In summary, the S-LASSO estimator outperforms the competitors both in terms of estimation error on the null region and prediction accuracy on a new dataset, as well as that it is able to estimate competitively the coefficient function on the non-null region. On the other hand, in order to achieve sparseness, the S-LASSO tends to overshrink the estimator of the coefficient function on the non-null region. This means that, as in the classical setting (James and Radchenko, 2009), there is a trade-off between the ability of recovering sparseness and the estimation accuracy on the non-null region of the final estimator.
Moreover, even when the coefficient function is not sparse (Scenario IV), the proposed estimator demonstrates to have both good prediction and estimation performance. This is another key property of the proposed estimator that, encourages practitioners to use the S-LASSO estimator even when there is not prior knowledge about the shape of the coefficient function. Finally, it should be noticed that, in scenarios similar to those analysed, the PCA and RIDGE estimators should not be preferred with respect to the TRU, SMOOTH and S-LASSO ones, and that the PFFR performance is strictly related to the REML approach to select the smoothing parameters.

Real-Data Examples
In this section, we analyse two real-data examples. We aim to confirm that the S-LASSO estimator has advantages in terms of both prediction accuracy and interpretability, over the SMOOTH estimator, which has been demonstrated in Section 4 to be the best alternative among the competitors. The datasets used in the examples are the Canadian weather and Swedish mortality. Both are classical benchmark functional data sets thoroughly studied in the literature.

Canadian Weather Data
The Canadian weather data have been studied by Ramsay and Silverman (2005) and Sun et al. (2018). The data set contains the daily mean temperature curves, measured in Celsius degree, and the log-scale of the daily rainfall profiles, measured in millimeter, recorded at 35 cities in Canada. Both temperature and rainfall profiles are obtained by averaging over the years 1960 through 1994. Figure 6 shows the profiles. The aim is to predict the log-daily rainfall based on the daily temperature using the model reported in Equation (1). Figure 7 shows the S-LASSO and SMOOTH estimates of the coefficient function β. The SMOOTH estimate is obtained using a Fourier basis-to take into account the periodicity of the data-and roughness penalties were chosen by using 10-fold cross-validation over an opportune grid of values. 10-fold cross-validation is used to set the parameters λ s , λ t and λ L as well.
The S-LASSO estimates is roughly zero over large domain portions. In particular, except for values from July through August, it is always zero in summer It is worth noting that the S-LASSO estimate is more interpretable than the SMOOTH estimates, which does not allow for a straightforward interpretation.
Moreover, the S-LASSO estimate appears to have, even if slightly, better prediction performance than the SMOOTH one. Indeed, 10-fold cross-validation mean squared errors are 22.314 and 22.365, respectively.
Finally, we perform two permutation tests to asses the statistical significance of the S-LASSO estimator. The first test is based on the global functional coefficient of determination defined as R 2 dt (Horváth and Kokoszka, 2012), with T = [0, 365]. In Figure 8(a) the solid black line indicates the observed R 2 g that is equal to 0.55. The bold points represent 500 R 2 g values obtained by means of random permutations of the response variable. Whereas, the grey line correspond to the 95th sample percentile. All 100 values of R 2 as well as the value of the 95th sample percentile is far below 0.55, which gives a strong evidence of a significant relationship between rainfalls and temperature, globally.
By a second test, we aim to analyse the pointwise statistical significance, i.e., for each t ∈ T . It is based on the pointwise functional coefficient of determination for t ∈ T (Horváth and Kokoszka, 2012). 95th sample percentile curve. The latter has been obtained by means of 500 R 2 values produced by randomly permuting the response variable. The observed R 2 is far above the 95th sample percentile curve, except for some summer months (viz., July and August). As global conclusion, we can state that the temperature has a large influence on the rainfalls in autumn, winter and spring.

Swedish Mortality Data
The Swedish mortality data, available from the Human Mortality Database (http://mortality.org), are regarded as a very reliable dataset on long-term longitudinal mortalities. In particular, we focus on 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-1894 with ages 0-80. The value of a log-hazard rate function at a given age is the natural logarithm of the ratio of the number of females who died at that age and the number of females alive with the same age. Note that, those data have been analysed also by Chiou and Müller (2009) and Ramsay et al. (2009). Figure 9 shows the 144 log-hazard functions.
The aim of the analysis is to explore the relationship of the log-hazard rate  In particular, at given t, the S-LASSO estimate is different from zero in an interval located right after that age. This can likely support the conjecture that if an event influences the mortality of the Swedish female at a given age, it impacts on the the death rate below that age born in the following years. Nevertheless, this expected dependence is poorly pointed out by the SMOOTH estimator, where this behaviour is confounded by less meaningful periodic components. It is interesting to note that the S-LASSO estimate at high values of t is slightly different from zero for ages ranging from 40 to 60. This shows that if an event affecting the death rate occurs in that range, the log-hazard functions of the following cohorts will be influenced at high ages (i.e., corresponding to high values of t). On the contrary, the wiggle of the SMOOTH estimate does not allow drawing such conclusions.
Finally, we perform the two permutation tests already described in the Canadian weather data example. Figure 11 shows the results. Both the observed R 2 g and R 2 are far above the 95th sample percentile (Figure 11(a)) and the pointwise 95th sample percentile curve (Figure 11(b)) respectively. This significantly evidences a relation between two consecutive log-hazard rate functions for all ages.

Ship CO 2 Emission Data
The ship CO 2 emission data have been previously studied in (Lepore et al., 2018;Reis et al., 2019;Capezza et al., 2020;Centofanti et al., 2020b,a). The dataset, provided by the shipping company Grimaldi Group, regards the issue of monitoring fuel consumptions or CO 2 emissions for a Ro-Pax ship that sails along a route in the Mediterranean Sea. The aim of the analysis is to study the relation between the fuel consumption per hour (FCPH), which is a proxy of the CO 2 emissions, and the speed over ground (SOG), assumed as predictor.
The observations considered were recorded from 2015 to 2017. Figure 12 shows the 44 available observations of SOG and FCPH (Centofanti et al., 2020b). Then, Figure 13 displays the S-LASSO and SMOOTH estimates of coefficient function β estimated as described in the Swedish mortality example. The S-LASSO coefficient function estimate also in this example is different from zero over a small portion of domain. Specifically, the FCPH during the navigation phase (i.e., t ∈ [0.1, 0.9]) is positively influenced by the SOG, in three specific voyage instants, viz., s ≈ 0.15, 0.55, 0.85, where the coefficient function estimate is positive. Thus, the FCPH during the navigation phase positively depends on the SOG at the start, in the middle and at the end of the navigation phase. Differently, during acceleration (t ∈ [0, 0.1]) and deceleration (t ∈ [0.9, 1]) phases, the relationship between the FCPH and the SOG is mainly concurrent. That is, the FCPH observed at a given time instant is strictly related to the SOG observed at the same time. On the contrary, such interpretation of the FCPH-SOG relationship is not easily obtained through the analysis of the SMOOTH coefficient function estimate, which is overall different from zero. Moreover, the S-LASSO and the SMOOTH estimates achieve 10-fold cross-validation mean squared errors of 0.077 and 0.093, respectively. Thus, the proposed estimator achieves slightly better prediction performance than the competitor. Figure 14 shows the results for the two permutation tests already described in the Canadian weather and Swedish mortality data examples. The relation between the FCPH and the SOG can be then considered as significant, as both the observed R 2 g and R 2 are above the 95th sample percentile (Figure 14(a)) and the pointwise 95th sample percentile curve (Figure 14(b)), respectively.

Conclusion
The LASSO is one of the most used and popular method to estimate coefficients in classical linear regression models as it ensures both prediction accuracy and interpretability of the phenomenon under study (by simultaneously performing variable selection). In this paper, we propose the S-LASSO estimator for the coefficient function of the Function-on-Function (FoF) linear regression model, which is an extension to the functional setting of the multivariate LASSO estimator. As the latter, the S-LASSO estimator is able to increase both the prediction accuracy of the estimated model, via continuous shrinking, and the interpretability, by identifying the null region of the regression coefficient, i.e., the region where the coefficient function is exactly zero.
The S-LASSO estimator is obtained by combining several elements: the functional LASSO penalty (FLP), which has the task of shrinking towards zero the estimator on the null region; the B-splines, which are essential to ensure sparsity of the estimator because of the compact support property; and two roughness penalties, which are needed to ensure smoothness of the estimator on the non-null region, also when the number of B-splines escalates. We proved that the S-LASSO estimator is both estimation and point-wise sign consistent, i.e., the estimation error in terms of L 2 -norm goes to zero in probability and the S-LASSO estimator has the same sign of the true coefficient function with probability tending to one. Moreover, we showed via an extensive Monte Carlo simulation study that, with respect to other methods that have already appeared in the literature, the S-LASSO estimator is much more interpretable, on the one hand, and has still good estimation and appealing predictive performance, on the other. However, consistently with the behaviour of the classical LASSO estimator (Fan et al., 2004), the S-LASSO estimator is found sometimes to over-shrink the coefficient function over the non-null region.
Note that the S-LASSO method is in the spirit of the fused LASSO of Tibshirani et al. (2005) in the multivariate context. Both methods rely on penalties that encourage sparsity in the coefficients and impose smoothness in the coef-ficient profile. The smoothness penalties of the proposed S-LASSO estimator and the fused LASSO approach have a different nature. In the fused LASSO, the differences of two consecutive coefficients are in fact penalized in their absolute value thus shrinking consecutive coefficients to the same value. On the contrary, the S-LASSO estimator relies on quadratic smoothness penalties that do not enjoy the sparseness property. Moreover, while the fused LASSO considers penalization of the first differences, the S-LASSO estimator allows the penalization of several types of smoothness, through the ms-th and mt-th order linear differential operators applied to the coefficient function.
To the best of the authors knowledge, this is the first work that addresses the issue of interpretability, intended as sparseness of the coefficient function, in the FoF linear regression setting. However, although the FLP produces an estimator with good properties, other penalties, e.g. the SCAD (Fan and Li, 2001) and adaptive LASSO (Zou, 2006), properly adapted to the functional setting, may guarantee even better performance both in terms of interpretabilty and prediction accuracy, and are, indeed, subjects of ongoing research.