Sparse and Smooth Functional Data Clustering

A new model-based procedure is developed for sparse clustering of functional data that aims to classify a sample of curves into homogeneous groups while jointly detecting the most informative portions of domain. The proposed method is referred to as sparse and smooth functional clustering (SaS-Funclust) and relies on a general functional Gaussian mixture model whose parameters are estimated by maximizing a log-likelihood function penalized with a functional adaptive pairwise penalty and a roughness penalty. The former allows identifying the noninformative portion of domain by shrinking the means of separated clusters to some common values, whereas the latter improves the interpretability by imposing some degree of smoothing to the estimated cluster means. The model is estimated via an expectation-conditional maximization algorithm paired with a cross-validation procedure. Through a Monte Carlo simulation study, the SaS-Funclust method is shown to outperform other methods already appeared in the literature, both in terms of clustering performance and interpretability. Finally, three real-data examples are presented to demonstrate the favourable performance of the proposed method. The SaS-Funclust method is implemented in the $\textsf{R}$ package $\textsf{sasfunclust}$, available online at https://github.com/unina-sfere/sasfunclust.


Introduction
In the last years, due to advances in technology and computational power, most of the data collected by practitioners and scientists in many fields bring information about curves or surfaces that are apt to be modelled as functional data, i.e., continuous random functions defined on a compact domain.A thorough overview of functional data analysis (FDA) techniques can be found in Ramsay and Silverman (2005); Ramsay et al. (2009); Horv'ath and Kokoszka (2012); Hsing and Eubank (2015) and Kokoszka and Reimherr (2017).As in the classical (non-functional) statistical literature, cluster analysis is an important topic in FDA, with many applications in various fields.The primary concern of functional clustering techniques is to classify a sample of data into homogeneous groups of curves, without having any prior knowledge about the true underlying clustering structure.The clustering of functional data is generally a difficult task because of the infinite dimensionality of the problem.For this reason, methods for functional data clustering have received a lot of attention in recent years, and different approaches have been proposed and discussed in the last decade.To the best of authors' knowledge, the most used approach is the filtering approach (Jacques and Preda, 2014), which relies on the reduction of the infinite dimensional problem by approximating functional data in a finite dimensional space and, then, uses traditional clustering tools on the basis expansion coefficients.Along this line, Abraham et al. (2003) propose an advanced version of the k-means algorithm to the coefficients obtained by projecting the original profiles onto a lower-dimensional subspace spanned by B-spline basis functions.A similar method is proposed by Rossi et al. (2004) who apply a Self-Organizing Map (SOM) on the resulting coefficient instead of the k-means algorithm.Elaborating on this path, Serban and Wasserman (2005) present a technique for the nonparametric estimation and clustering of a large number of functional data that is still based on the k-means algorithm applied to the basis expansion coefficients obtained through smoothing techniques.A step forward is moved by Chiou and Li (2007), who introduce the k-centers functional clustering method to account, differently from the previous methods, for both the means and the mode of variation differentials between clusters by predicting cluster membership with a reclassification step.
Instead of considering the basis expansion coefficients as parameters, a different idea is that of using a model-based approach where coefficients are treated as random variables themselves with a cluster-specific probability distribution.The seminal work of James and Sugar (2003) is the first one to develop a flexible model-based procedure to cluster functional data based on a random effects model for the coefficients.This allows for borrowing strength across curves and, thus, for superior results when data contain a large number of sparsely sampled curves.More recently, Bouveyron and Jacques (2011) propose a new functional clustering method, which is referred to as funHDDC and based on a functional latent Gaussian mixture model, to fit the functional data in group-specific functional subspaces.
By constraining model parameters within and between groups, they obtain a family of parsimonious models that allow for more flexibility.Analogously, Jacques and Preda (2013) assume cluster-specific Gaussian distribution on the principal components resulting from the Karhunen-Loeve expansion of the curves, and Giacofci et al. (2013) propose to use a Gaussian mixture model on the wavelet decomposition of the curves, which turns out to be particularly appropriate for peak-like data, as opposed to methods based on splines.
In the multivariate cluster analysis, some attributes could be, however, completely noninformative for uncovering the clustering structure of interest.As an example, this often happens in high-dimensional problems, i.e., where the number of variables is considerably larger than the number of observations.In this setting, the task of identifying the features in which respect true clusters differ the most is of great interest to achieve (a) a more accurate identification of the groups, as noninformative features may hide the true clustering structure, and (b) an higher interpretability of the analysis, by imputing the presence of the clustering structure to a small number of features.More in general, the methods capable of selecting informative features and eliminating noninformative ones are referred to as sparse.Such class of methods can be usually reconducted and regarded as variable selection methods.Sparse clustering has received increasing attention in the recent literature.Based on conventional heuristic clustering algorithms, Friedman and Meulman (2004) develop a new procedure to automatically detect subgroups of objects, which preferentially cluster on subsets of features.Witten and Tibshirani (2010) elaborate a novel clustering framework based on an adaptively chosen subset of features that are selected by means of a lasso-type penalty.In terms of model-based approaches, the method introduced by Raftery and Dean ( 2006) is able to sequentially compare nested models through approximate Bayes factor and to select the informative features.Maugis et al. (2009) improve this method by considering the noninformative features as independent from some informative ones.
It is moreover worth mentioning quite promising variable selection approaches that make use of a regularization framework.The seminal work in this direction is that of Pan and Shen (2007), who introduce a penalized likelihood approach with an L 1 penalty function, which is able to automatically achieve variable selection via thresholding and delivering a sparse solution.Similarly, Wang and Zhu (2008) suggest a solution by replacing the L 1 penalty with either the L ∞ penalty or the hierarchical penalization function, which take into account the fact that cluster means, corresponding to the same feature, can be treated as grouped.Xie et al. (2008) also account for grouped parameters through the use of two planes of grouping, named vertical and horizontal grouping.In all sparse clustering methods just mentioned, a feature is selected if it is informative for at least one pair of clusters and eliminated otherwise, i.e., if it is noninformative for all clusters.However, some variables could be informative only for specific pairs of clusters.For this reason, Guo et al. (2010) propose a pairwise fusion penalty that penalizes, for each feature, the differences between all pairs of cluster means and fuses only the non separated clusters.
Only recently, the notion of sparseness has been translated into a functional data clustering framework.Specifically, sparse functional clustering methods aim to cluster the curves while jointly detecting the most informative portion of domain to the clustering in order to improve both the accuracy and the interpretability of the analysis.Based on the idea of Chen et al. (2014), Floriello and Vitelli (2017) propose a sparse functional clustering method based on the estimation of a suitable weight function that is capable of identifying the informative part of the domain.Vitelli (2019) proposes a novel framework for sparse functional clustering that also embeds an alignment step.To the best of the authors' knowledge, these are the only works that propose sparse functional clustering methods so far.
In this article, we present a model-based procedure for the sparse clustering of functional data, named sparse and smooth functional clustering (SaS-Funclust) method, where the basic idea is to provide accurate and interpretable cluster analysis.Specifically, the parameters of a general functional Gaussian mixture model are estimated by maximizing a penalized version of the log-likelihood function, where a functional adaptive pairwise fusion penalty, the functional extension of the penalty proposed by Guo et al. (2010), is introduced.Firstly, it penalizes the pointwise differences between all pairs of cluster functional means and locally shrinks the means of cluster pairs to some common values.Secondly, a roughness penalty on cluster functional means is considered to further improve the interpretability of the cluster analysis.Therefore, the SaS-Funclust method gains the ability to detect, for each cluster pair, the informative portion of domain to the clustering, hereinafter always intended in terms of mean differences.If a specific mean pair is fused over a portion of the domain, it is labelled as noninformative to the clustering of that pair.
Otherwise, it is labelled as informative.In other words, the proposed method is able to detect portions of domain that are noninformative pairwise, i.e., for at least a specific cluster pair, differently from the method proposed by Floriello and Vitelli (2017) that is only able to detect portions of domain that are noninformative overall, i.e., for all the cluster pairs simultaneously.Moreover, the model-based fashion of the proposed method provides greater flexibility than the latter, which basically relies on k-means clustering.A specific expectation-conditional maximization (ECM) algorithm is designed to perform the maximization of the penalized log-likelihood function, which is a non-trivial problem, and a cross-validation based procedure is proposed to select the appropriate model.The method presented in this article is implemented in the R package sasfunclust, openly available online at https://github.com/unina-sfere/sasfunclust.
To give a general idea of the sparseness property of the proposed method, Figure 1 shows the cluster means estimated by the latter for three different simulated data sets with (a) two, (b) three and (c) four clusters.Data are generated as described in Section 3. In Figure 1(a), the estimated means are correctly fused over t ∈ (0.2, 1.0].Hence, the proposed method is shown to be able to identify the informative portion of domain [0.0, 0.2], for the unique pair of clusters and not for all.In Figure 1(b) and Figure 1(c), several cluster pairs are available, because the number of clusters is larger than 2, and, thus, a given portion of domain could be informative for a specific pair of clusters.In Figure 1(b), the informative portion of domain for each pair of clusters is correctly recovered.The estimated cluster means are indeed pairwise fused over approximately the same portion of domain as the true cluster means pairs.Note that, for the clusters whose true means are equal over t ∈ (0.2, 1.0], the SaS-Funclust method identifies the informative portion of domain roughly in [0.0, 0.2].In Figure 1(c), the sparseness property of the SaS-Funclust method is even more striking.In this case, in the face of many cluster pairs, the proposed method is still able to successfully detect the informative portion of domain.The properties of the proposed method will be deepened in Section 3.
The remainder of this article is organized as follows.Section 2 presents the proposed methodology.Specifically, Section 2.1 and 2.2 introduce the general functional Gaussian mixture model and the penalized maximum likelihood estimator, respectively.Whereas, the optimization algorithm and model selection considerations are discussed in Section 2.3 and Section 2.4, respectively.Properties and performance of the SaS-Funclust method are assessed through a wide Monte Carlo simulation study presented in Section 3. Section 4 illustrates the potentiality of the SaS-Funclust method by means of of three real datasets: the Berkeley Growth Study, the Canadian weather, and the ICOSAF project data.
2 The SaS-Funclust method for functional clustering 2.1 A general functional clustering model The SaS-Funclust method is based on the general functional clustering model introduced by James and Sugar (2003).Suppose that N observations are spread among G unknown clusters, and that the probability of each observation to belong to the gth cluster is π g , G g=1 π g = 1.Moreover, let us denote with Z i = (Z 1i , . . ., Z Gi ) T the unknown componentlabel vector corresponding to the ith observation, where Z gi equals 1 if the ith observation is in the gth cluster and 0 otherwise.Then, let us assume that for each i observation, i = 1, . . ., N in the cluster g = 1, . . ., G, it is available a vector Y i = (y i1 , . . ., y in i ) T of size n i , which can differ across observations, of observed values of a function g i over the time points t i1 , . . ., t in i .The function g i is assumed a Gaussian random process with mean µ g , covariance ω g , and values in L 2 (T ), the separable Hilbert space of square integrable functions defined on the compact domain T .We assume that, conditionally on that Z gi = 1, Y i is modelled as where g i = (g i (t i1 ) , . . ., g i (t in i )) T contains the values of the function g i at t i1 , . . ., t in i and i is a vector of measurement errors that are mutually independent and normally distributed with mean 0 and constant variance σ 2 e .Let us suppose also that the unknown componentlabel vector Z i has a multinomial distribution, which consists of one draw on g categories with probabilities π 1 , . . ., π G .Then, for every i, the unconditional density function f (•) of where As discussed in James and Sugar (2003), it is necessary to impose some structure curves g i because both the curves could be observed at different time domain points and the dimensionality of the model in Equation ( 2) could be too high in comparison to the sample size.Therefore, similarly to the filtering approach for clustering, we assume that each function g i , for i = 1, . . ., N , may be represented in terms of a q-dimensional set of basis functions Φ = (φ 1 , . . ., φ q ) T , that is where η i = (η i1 , . . ., η iq ) T are vectors of basis coefficients.Then, η i are modelled as Gaussian random vectors, that is, given that Z gi = 1, where µ g = (µ g1 , . . ., µ gq ) T are q-dimensional vectors and γ ig are Gaussian random vectors with zero mean and covariance Γ g .
With these assumption the unconditional density function where T is the basis matrix for the ith curve and e .Therefore, the log-likelihood function corresponding to Y 1 , . . ., Y N is given by where Θ = {π g , µ g , Γ g , σ 2 e } g=1,...,G is the parameter set of interest.Based on an estimate Θ = {π g , μg , Γg , σ2 e } g=1,...,G , an observation Y * is assigned to the cluster g that achieves the largest posterior probability estimate πg ψ Y * ; S i μg , Σig , with Σig = S i Γg S T i + I σ2 e .

The penalized maximum likelihood estimator
work, we propose a different estimator ΘPMLE of Θ that is the maximizer of the following penalized log-likelihood where P (•) is a penalty function defined as where λ L , λ s ≥ 0 are tuning parameters, and τ g,g are prespecified weight functions.The symbol f (s) (•) denotes the sth-order derivative of f if the latter a function or the entries of f if it is a vector.Note that in Equation ( 8) each function g i may be represented as in Equation ( 3), then it follows that The first element of the right-hand side of Equation ( 8) is the functional extension of the penalty introduced by Guo et al. (2010) and is referred to as functional adaptive pairwise fusion penalty (FAPFP).The aim of the FAPFP is to shrink the differences between every pair of cluster means for each value of t ∈ T .Due to the property of the absolute value function of being singular at zero, some of these differences are shrunken exactly to zero.
In particular, the FAPFP allows pair of cluster means to be equal over specific portion of domain that are, thus, noninformative for separating the means of that pair of clusters.
The choice of the weight function τ g,g in Equation ( 8) and Equation ( 8) should be based on the idea that if a given portion of domain is informative for separating the means of the corresponding pair of clusters, then, the values of τ g,g over that portion should be small.
In this way, the absolute difference |µ g (•) − µ g (•) | is penalized more over noninforvative portions of domain than over informative ones.Following the standard practice for the adaptive penalties (Zou, 2006), we propose to use where μg are initial estimates of the cluster means.
Finally, the term λ s G g=1 T µ (s) 2 dt is a smoothness penalty that penalizes the sth derivative of the cluster means.This term aims to further improve the interpretability of the results by constraining, of a magnitude quantified by λ s , the cluster means to own a certain degree of smoothness, measured by the derivative order s.Following the common practice in FDA (Ramsay and Silverman, 2005), we suggest to penalize the cluster mean curvature by setting s = 2, which implies that the basis functions chosen are differentiable at least s times.As a remark, the penalization in Equation ( 7) is applied only to the parameter vectors corresponding to the cluster means, i.e., to µ 1 , . . ., µ G .The reason is that, in this work, we consider the case where a portion domain is defined as informative only in terms of cluster mean differences.However, portions of domain could be informative also in terms of differences in cluster covariances, which together with the means uniquely identify each cluster.

The penalty approximation and the optimization algorithm
To perform the maximization of the penalized log-likelihood in Equation ( 7), the penalty P (•), defined as in Equation ( 8), can be written as where the weight functions τ g,g (t) are expressed as in Equation ( 10), and the initial estimates of the cluster means are represented through the set of basis functions Φ as μg A great simplification of the optimization problem can be achieved if the first element on the right-hand side of Equation ( 11) can be expressed as linear function of |µ T g − µ g |.The following theorem provides a practical way to rewrite the first term of the right-hand side of Equation ( 11) as linear function of |µ T g − µ g |, when Φ is a set of B-splines (De Boor et al., 1978;Schumaker, 2007).
Theorem 1.Let Φ = (φ 1 , . . ., φ q ) T be the set of B-splines of order k and non-decreasing knots sequences {x 0 , x 1 , . . ., x M j , x M +1 } defined on the compact set T = [x 0 , x M +1 ], with q = M +k, and {τ j } q+1 j=1 being a sequence with and zero elsewhere, is such that where δ = max i |x i+1 − x i |, that is f (t) − f (t) converges uniformly to the zero function.
Theorem 1, whose proof is deferred to the Supplementary material, basically states that when δ is small, f is well approximated by f .In other words, the approximation error |f − f | can be made arbitrarily small by increasing the number of knots.If we further assume the knots sequence evenly spaced, δ turns out to be equal to 1/M .These considerations allow where T .Thus, Equation ( 11) can be rewritten as where M = diag is the diagonal matrix with diagonal entries The goodness of the approximations in Equation ( 13) depends on the cardinality q of the set of B-splines Φ, which, thus, should be as large as possible.However, the number of parameters in Equation ( 2), which depends quadratically on q, becomes very large even for moderate values of q.To mitigate this issue, one may further assume equal and diagonal coefficient covariance matrices across all clusters, that is As a remark, with this assumption, we are implicitly assuming that clusters are separated only by their mean values, and, thus, informative portion of domain are identified only by cluster mean differences and not in terms of covariances.
The penalized log-likelihood function in Equation ( 7) becomes with The maximization of this objective function is a nontrival problem.A specifically designed algorithm is proposed, which is a modification of the expectation maximization (EM) algorithm proposed by James and Sugar (2003).By treating the component-label vectors Z i (defined at the beginning of Section 2.1) and γ ig in Equation (4) as missing data, the complete penalized log-likelihood is given by The EM algorithm consists in the maximization, at each iteration t = 0, 1, 2, . . ., of the expected value of L cp , calculated with respect the joint distribution of Z i and γ ig , given Y 1 , . . ., Y N and the current parameter estimates Θ(t) = {π q , σ2(t) e } g=1,...,G The algorithm stops when a pre specified stopping condition is met.At each t, the expected value of L cp as a function the probability of membership π 1 , . . ., π G is then maximized by setting π(t+1) with π(t+1 where γ 2 ig(j) indicates the jth entry of γ 2 ig .The value of E γ 2 ig(j) |Y i , Z gi = 1, , Θ(t) can be calculated by using the property that the (conditional where γ(t) ig = E γ ig |Y i , Z gi = 1, Θ(t) .The mean vectors µ 1 , . . ., µ G that maximize the conditional expectation of L cp are the solution of the following optimization problem The optimization problem in Equation ( 20) is a difficult task of the non differentiability of the absolute value function in zero, and, it has not a closed form solution.However, following the idea of Fan and Li (2001), it can be solved by means of the standard local quadratic approximation method, i.e., by iteratively solving the following quadratic optimization problem for s = 0, 1, 2, . . . where

|
, and μ(t+1,0) 21) is based on the following approximation (Fan and Li, 2001) The solution to the original problem in Equation ( 20) can be satisfactorily approximated by the solution at iteration s * of the optimization problem in Equation ( 22) when a pre specified stopping condition is met, i.e., μ(t+1) For numerical stability, a reasonable suggestion is to set a lower bound on |μ and then to shrink to zero all the estimates below the lower bound.It is worth noting that the proposed modification to the algorithm of James and Sugar (2003) falls within the class of the expectation conditional maximization (ECM) algorithms (Meng and Rubin, 1993).Based on the convergence property of the ECM algorithms, which also holds for the local quadratic approximation in variable selection problems (Fan and Li, 2001;Hunter and Li, 2005), the proposed algorithm can be proved to converge to a stationary point of the penalized log-likelihood in Equation ( 15).

Model selection
The proposed SaS-Funclust method requires the choice several hyper-parameters viz., the number of clusters G, tuning parameters λ s , λ L , dimension q and the order k of the set of B-spline functions Φ as well as the knot locations.A standard choice for Φ is the cubic B-splines (i.e., k = 4) with equally spaced knot sequence, which enjoy the optimal property of interpolation (De Boor et al., 1978).Moreover, the dimension q should be set as large as possible to reduce, to the greatest possible extent, the approximation error in Equation ( 13).This facilitates the estimated cluster means to successfully capture the local feature of the true cluster means.Unfortunately, the larger the value of q, the higher the complexity of the model in Equation (2), i.e., the number of parameters to estimate.The presence of the smoothness penalty on µ g , as well as the constraint imposed on Γ g , allows one to control the complexity of the model and, thus, to prevent over-fitting issues.The choice of G, λ s , and λ L may be based on a K-fold cross-validation procedure.Based on observations divided into K equal-sized disjoint subsets f 1 , . . ., f k , . . ., f K , G, λ s , and λ L are chosen as the maximizers of the following function λ L , with m = m 2 ; thirdly, to choose λ L at fixed λ s and G, by using m = m 3 .In this way, the estimated model is not unnecessarily complex and achieves predictive performance that is comparable to that of the best model (i.e., the one that maximizes the CV function in Equation ( 23)).
As a remark, although the component-wise procedure proposed to choose λ s , λ L and G proves itself to be very effective in the simulation study of Section 3, we recommand whenever possible to directly plot and inspect the CV curve as a function of G, λ s , and λ L and to use any information available from the specific application.

Simulation study
In this section, the performance of the SaS-Funclust method is assessed by means of an extensive Monte Carlo simulation study.The SaS-Funclust method, implemented through the R package sasfunclust, is compared with the following methods that have already appeared in the literature before.In particular, we refer to the method proposed by Giacofci et al.
These methods are implemented through the homonymous R packages curvclust (Giacofci et al., 2012) and funHDDC (Schmutz and Bouveyron, 2019).In addition, we consider also filtering approaches based on two main steps.The first step consists in the estimation of the functions g i by means of either smoothing B-splines or functional principal component analysis (Ramsay and Silverman, 2005); whereas the second step aims to apply standard clustering algorithms, viz.hierarchical, k-means and finite mixture model clustering methods (Everitt et al., 2011), on either the resulting B-spline coefficients or the functional principal components scores.Filtering approaches based on the smoothing B-splines and the hierarchical, k-means and finite mixture model clustering methods will be hereinafter referred to as B-HC, B-KM and B-FMM, respectively, whereas methods based on the functional principal component analysis and the hierarchical, k-means and finite mixture model clustering methods are referred to as FPCA-HC, FPCA-KM and FPCA-FMM.Finally, we evaluate also the method presented by Ieva et al. (2013), which is referred to as DIS-KM and it basically consists in the application of the k-means clustering to the L 2 distances among the observed curves.
The number of clusters is selected through the Bayesian information criterion (BIC) for the curvclust and funHDDC methods, as suggested by Giacofci et al. (2013) and Bouveyron and Jacques (2011), respectively; whereas the silhouette index (Rousseeuw, 1987) is used for the DIS-KM method.The majority rule applied to several validity indices (Charrad et al., 2014) is used to determine the number of clusters for all the filtering approaches.
The number of clusters and the tuning parameters needed to implement the SaS-Funclust method are determined through the CV based procedure described in Section 2.4 with q = 30, K = 5, m 1 = m 3 = 0.5, and m 2 = 0.The values of m 1 and m 3 ensure parsimony in the choice of λ L and G, whereas for picking λ s the m-standard deviation is not applied.
The initial values of the parameters for the ECM algorithm are chosen by applying the k-means algorithm on the coefficients estimated through smoothing B-spline.
The performance of the clustering procedures in selecting the proper number of clusters and identifying the clustering structure, when the true number of cluster is known, is assessed separately.In particular, the former is measured through the mean number of selected clusters, whereas the latter is compared through the adjusted Rand index (Hubert and Arabie, 1985) denoted by aRand.This index accounts for the agreement between true data partitions and clustering results corrected by chance, based on the number of paired objects that are either in the same group or in different groups in both partitions.
The aRand yields values between 0 and 1.The larger its value, the higher the similarity between the two partitions.
Three different scenarios are analysed where data are generated from G t = 2, 3, 4 clusters  The true functions are obtained as g i = η T i Θ B , with Θ B representing a set of 30 evenly spaced knot cubic B-splines.The coefficients η i = (η i1 , . . ., η i30 ) T are Gaussian random coefficients with mean vectors (depending on both the scenario and the cluster considered) reported in Table 1, and covariance matrix Γ = σ 2 c I, where σ c = 0.5 and I denotes the identity matrix.Then, to obtain the contaminated-with-error observations Y i , the true functions g i are evaluated over a grid of 50 points, with the addition of measurement errors independently generated as Gaussian random variables with zero mean and five different values of standard error σ e = 1, 1.5, 2, 2.5, 3. Note that, from Scenario I to Scenario III, the portion of domain that is noninformative for all cluster pairs decreases, whereas, the portions of domain that are informative for specific cluster pairs increases.
Figure 2 shows the average aRand index values for Scenario I, through III as a function of the standard error σ e .In Scenario I, at small values of σ e , all methods perform comparably and provide clustering partitions with aRand values very close to 1, which corresponds to the perfect cluster identification.However, as σ e increases, the SaS-Funclust method turns out to be the best method, closely followed by the curvclust method.Also the B-FMM performs very well, except when σ e = 3.0.In Scenario II and III, the SaS-Funclust method is still the best, followed by the curvclust and B-FMM case in Scenario II and only by the curvclust method in Scenario III.Note that in these scenarios, the DIS-KM underperforms also in the most favourable cases as a consequence of the lesser capacity of the L 2 distance to recover the true clustering structure.
Figure 3 shows the mean number of selected clusters in all scenarios.It is clear that the SaS-Funclust method is able to identify the true number of clusters much better than the competitors in all the considered scenarios.In particular, Scenario II highlights that, especially for large measurement error σ e , the competing methods reduce their complexity and select, on average, a number of clusters smaller than the true number of clusters G t = 3.This is evident in Scenario III, where the competing methods select, on average, a number of clusters G = 2 for σ e = 2.5, 3.0, which is smaller than G t = 4. Figure 4 and Table 2 highlight the ability of the SaS-Funclust method in recovering the true cluster means and detecting the informative portions of domain.The average root mean squared error calculated as , with μg (t) = μT g Φ (t), t ∈ T , is plotted in Figure 4 for each method as a function of σ e in all three scenarios.By this figure, the SaS-Funclust method outperforms the competitors in each scenario, especially for large measurement errors, even though the curvclust method shows comparable performance.Table 2 reports, for each σ e and scenario, the average fractions of correctly identified noninformative portions of domain by the SaS-Funclust method, which can be regarded as a measure of the interpretability (i.e., sparseness) of the proposed solution.In more detail, each entry of the table is obtained as the mean of the average fraction of correctly identified noninformative portions of domain, over the 100 generated datasets, for each pair of clusters, weighted by the size of the corresponding true noninformative portions of domain.In Scenario I, it trivially coincides with the average, because the true number of clusters is G t = 2.The proposed method is clearly able to provide an interpretable clustering.The fraction of correctly identified noninformative portions of domain is almost larger than or equal to 0.90 for σ e ≤ 2.5 and decreases to 0.80 for σ e = 3.5.It is worth noting that when σ e = 1.0, the pairs of clusters in each scenario are correctly fused over almost all the noninformative portion of domain in terms of mean differences.This confirms what is shown in Figure 1 of Section 1.  4 Real-data Examples

Berkeley Growth Study Data
In this section, the SaS-Funclust method is applied to the growth dataset from the Berkeley growth study (Tuddenham, 1954), which is available in the R package fda (Ramsay et al., 2020).In this study, the heights of 54 girls and 39 boys were measured 31 times at age 1 through 18.The aim of the analysis is to cluster the growth curves and compare the results with the partition based on the gender difference.This problem has been already addressed by Chiou and Li (2007); Jacques and Preda (2013); Floriello and Vitelli (2017).
In particular, we focus on the growth velocities from age 2 to 17, whose discrete values are estimated through the central differences method applied to the growth curves.Figure 5 shows the interpolating growth velocity curves for all the individuals.In view of the analysis objective, all clustering methods described in Section 3 are applied by setting G = 2.As shown in the first row of Table 3, all clustering methods, excluded the B-HC, perform similarly in terms of the aRand index with respect to the gender difference partition.Moreover, by looking at the second row of Table 3, which shows the aRand index with respect to the SaS-Funclust partition, the competing methods provide partitions very similar to the SaS-Funclust one.
As expected, the SaS-Funclust method allows for a more interpretable analysis.Figure 6 shows the estimated cluster means and the clustered growth curves for the SaS-Funclust method.The estimated cluster means are fused over the first portion of the domain, whereas they are separated over the remaining portions.This implies that the two identified clusters are not different on average over the first portion of domain which can be, thus, regarded as noninformative.Separation between the two group arises over the remaining informative portion of domain, where two sharp peaks of growth velocity arise, instead.
The latter peaks are known in the medical literature as pubertal spurts, in which respect the attained results indicate two main timing/duration groups.In particular, male pubertal spurt happens later and lasts longer than female one.Nevertheless, some individuals show unusual growth patterns that are not captured by the cluster analysis.Additionally, the estimated cluster means from the competing methods, not shown here, do not allow for a similar straightforward interpretation.

Canadian Weather Data
The Canadian weather dataset contains the daily mean temperature curves, measured in Celsius degree, recorded at 35 cities in Canada.The temperature profiles are obtained by averaging over the years 1960 through 1994.This is a benchmark dataset available in the R package fda (Ramsay et al., 2020) that has been already studied by Ramsay and Silverman (2005); Centofanti et al. (2020).Figure 7 displays the interpolating profiles, where, for computational reasons, temperature curves are sampled each five days.
The ultimate goal of the cluster analysis applied to these curves is the geographical  4 reports the aRand index for all the competing methods calculated with respect to the SaS-Funclust method.As expected, the proposed clustering agrees with filtering methods based on B-splines, while mostly disagrees with the others.
In terms of interpretability, Figure 8 shows the estimated cluster means and geographical distribution of the curves in the clusters obtained by the SaS-Funclust method.From Figure 8(a), the estimated means for clusters 1, 2 and 4 are shown to fuse approximately from day 100 through 250.This is a strong evidence that the mean temperature in this period of the year is not significantly different among zones in cluster 1, 2 and 4. Hence, this portion of domain turns out to be noninformative for the separation of these clusters, whereas the mean temperature is different for the rest of the year.A different pattern is followed by the curves in cluster 3, which shows significantly smaller mean temperature all over the year.The geographical displacement of the temperature profiles coloured by the clusters identified through the SaS-Funclust method is reported in Figure 8 cluster 3, which correspond to northern stations, show lower mean temperature.This nice and plausible interpretation of this well-known real-data example is not possible by means of any competing method.

ICOSAF Project Data
The ICOSAF project dataset contains 538 dynamic resistance curves (DRCs), collected during resistance spot welding lab tests at Centro Ricerche Fiat in 2019.The DRCs are collected over a regular grid of 238 points equally spaced by 1 ms.Further details on this dataset can be found in Capezza et al. (2020) and the data are publicly available online at https://github.com/unina-sfere/funclustRSW/.In this example, we focus on the first derivative of the DRCs, estimated by means of the central differences method applied to the DRC values sampled each 2 ms. Figure 9 shows the first derivative of the DRCs defined, without loss of generality, on the domain [0, 1].In this setting, the aim of the analysis is to cluster DRCs to identify homogenous groups of spot welds that share common mechanical and metallurgical properties.Differently from the previous datasets, no information are available about a reasonable partition of the DRCs.Therefore, based on the considerations provided by Capezza et al. (2020) and on cluster number selection methods described for the SaS-Funclust and competing methods in Section 2.4 and 3, respectively, we set G = 3.Table 5 shows aRand values obtained for all method pairs with respect to the SaS-Funclust partition.In this case, the SaS-Funclust method provides partitions that are more similar to those obtained through the FPCA-based methods than those obtained with the B-splines filtering approaches.However, the clusters identified by the SaS-Funclust method do not resemble those of the other methods.It is worth noting that, for this dataset, even if results are not reported here, the partition obtained by curvclust differs dramatically from the others and does not provide meaningful clusters.
Also in this case, the SaS-Funclust method allows for an insightful interpretation of the results.The estimated cluster means and the corresponding clustered curves obtained through the SaS-Funclust method displayed in Figure 10 confirm the ability of the proposed method to fuse cluster means, as it is clear over the second part of the domain.In particular, the mean of cluster 1 and 3 are fused from 0.5 to 1, which accounts for the comparable decreasing rate of the DRCs over these clusters.Differently, the mean of cluster 2 is fused with other cluster means between 0.8 and 1, only.This indicates that between 0.5 and

Conclusions and discussions
This article presented the SaS-Funclust method, a new approach to the sparse clustering of functional data.Differently from methods that have already appeared in the literature before, it was shown to be capable of successfully detecting where cluster pairs are separated.In many applications, this involves limited portions of domain, which are referred to as informative, and thus, the proposed method allows for a more accurate and interpretable cluster analysis.The SaS-Funclust method can be considered as belonging to the model-based clustering procedures with parameters of a general functional Gaussian mixture model estimated by maximizing a penalized version of the log-likelihood function.The key element is the functional adaptive pairwise fusion penalty that, by locally shrinking mean differences, allows pairs of cluster means to be exactly equal over portions of domain where cluster pairs are not well separated, referred to as noninformative.
In addition, a smoothness penalty is introduced to further improve cluster interpretability.The penalized log-likelihood function was maximized by means of a specifically designed expectation-conditional expectation algorithm, and model selection was addressed through a cross-validation technique.An extensive Monte Carlo simulation study showed the favourable performance of the proposed method over several competing methods both in terms of clustering accuracy and interpretability.Lastly, real-data examples further demonstrated the practical advantages of the proposed method, which provided, thanks to its sparseness property, new insightful and interpretable solutions to cluster analysis.In the Berkeley growth study example, the SaS-Funclust method highlighted that growth velocity curves of boys and girls show different pubertal spurt, which happens later and last longer for male than female.Whereas, in the Canadian weather example, the mean temperatures over the Pacific, Atlantic and southern continental regions were found to be equal over the middle days of the year and different otherwise.Moreover, the proposed method was applied to the ICOSAF project dataset, where, differently from the previous datsets, no information are available about a reasonable partition.The SaS-Funclust method also in this case identified homogenous groups of spot welds that showed, only during the first part of the process, differences in the rate of change of dynamic resistance curves, which are likely to be responsible of distinct mechanical and metallurgical properties of the spot welds.
As closing remarks, we can envisage several important extensions to refine the proposed method.Regarding the structure of the functional clustering model, the assumption of a common diagonal coefficient covariance matrix across all clusters may be too restrictive and may result in a poor fit.Unfortunately, more flexible covariance structures dramatically increase the number of parameters to be estimated, already enlarged to achieve sparseness, in the SaS-Funclust method.For this reason, regularization framework shall necessarily be addressed to avoid overfitting, possibly either by constraining the covariance structure, as done in this article, or by means of shrinkage estimators.However, the choice of the best approach still remains not straightforward.Furthermore, the covariance structure of the measurement errors could be modified to include more complex relationships, and the model can be extended also by including covariates (James and Sugar, 2003).
Another natural extension of the SaS-Funclust method that in worth considering is the integration of a proper pairwise penalty applied to the covariance functions, useful in those settings where portions of domain are informative for the clustering also in terms of covariance functions.Unfortunately, the choice of such penalty and the resulting computational issues are non-trivial and need for additional careful investigation.

Figure 1 :
Figure 1: True and estimated cluster means obtained through the SaS-Funclust method for three different simulated data sets with (a) two, (b) three and (c) four clusters generated as described in Section 3.
I is the identity matrix, and ψ (•; µ, Σ) is the multivariate Gaussian density distribution with mean µ and covariance Σ.The model in Equation (2) is the classical G-component Gaussian mixture model(McLachlan and Peel, 2004).
-Funclust estimates of π g , µ g and Σ i obtained by leaving out the observations in the k-th subset f k .Usually, the CV function is numerically calculated over a finite grid of values.As in the multivariate regression setting, the uncertainty of the CV function in estimating the log-likelihood observed for an out-ofsample observation is taken into account by means of the so called m-standard deviation rule.This heuristic rule suggests to pick up the most parsimonious model among those achieving values of the CV function that are no more than m standard errors below the maximum of Equation (23).Note that, in this problem, parsimony is reflected into large λ s , λ L and small G.By elaborating on the m-standard deviation rule, we propose to firstly choose G for each value of λ s , λ L , with m = m 1 ; secondly, at fixed G, choose λ s for each as Scenario I, II and III, respectively.For each scenario, the considered methods are evaluated by assessing the performance over 100 independently generated datasets.From each cluster, 200 observations are generated over the domain T = [0, 1].

Figure 2 :Figure 3 :
Figure 2: Average aRand index for (a) Scenario I, (b) Scenario II, and (c) Scenario III as a function of σ e when the true number of clusters is known.

Figure 4 :
Figure 4: Average selected number of clusters G for Scenario I (a), Scenario II (b), and, Scenario III (c) as a function of σ e .

Figure 5 :
Figure 5: Growth velocities of 54 girls and 39 boys in the Berkeley growth study dataset.

Figure 6 :Figure 7 :
Figure 6: (a) Estimated cluster curve means and (b) curve clusters for the SaS-Funclust method in the Berkeley growth study dataset.

Figure 8 :
Figure 8: (a) Estimated cluster curve mean and (b) geographical displacement of the curves pertaining to clusters obtained through SaS-Funclust method.

Figure 9 :
Figure 9: First derivatives of the 538 DRCs in the ICOSAF project dataset.

Figure 10 :
Figure 10: (a) Estimated cluster curve means and (b) curve clusters for the SaS-Funclust method in the ICOSAF project dataset.

Table 1 :
Coefficient mean vectors for each scenario and cluster.

Table 2 :
Average fractions of correctly identified noninformative portions of domain by the SaS-Funclust method for each scenario.

Table 3 :
The values of the aRand index for all the clustering methods with respect to gender difference based grouping and the SaS-Funclust partition for the Berkeley growth

Table 4 :
Jacques and Preda, 2013)index for all the clustering methods with respect to climate zones grouping and the SaS-Funclust partition for the Canadian weather datasetJacques and Preda, 2013).The first row of Table4shows the aRand index values of the resulting clusters calculated with respect to the 4-climatezone grouping.Although the SaS-Funclust and the B-HC methods achieve the largest aRand in this case, aRand values are in all cases inadequately low, which indicates the clustering structure disagrees with such grouping.That is, different method performance cannot properly evaluated by using the 4-climate-zone grouping.The second row of Table

Table 5 :
aRand index calculated on the ICOSAF project dataset for all competing method partitions with respect to the SaS-Funclust one.