Shape-Informed Clustering of Multi-Dimensional Functional Data via Deep Functional Autoencoders
Authors: Samuel Singh
Shirley Coyle
Mimi Zhang
Venue: NeurIPS 2025 poster ()
Keywords: autoencoder
cluster analysis
functional data
manifold
multi-dimensional
neural network
phase variation
universal approximation
Primary area: general_machine_learning
TL;DR: We developed a functional autoencoder framework for cluster analysis of multi-dimensional functional data, data that are random realizations of vector-valued random functions.

Abstract: We introduce FAEclust, a novel functional autoencoder framework for cluster analysis of multi-dimensional functional data, data that are random realizations of vector-valued random functions. Our framework features a universal-approximator encoder that captures complex nonlinear interdependencies among component functions, and a universal-approximator decoder capable of accurately reconstructing both Euclidean and manifold-valued functional data. Stability and robustness are enhanced through innovative regularization strategies applied to functional weights and biases.  Additionally, we incorporate a clustering loss into the network's training objective, promoting the learning of latent representations that are conducive to effective clustering. A key innovation is our shape-informed clustering objective, ensuring that the clustering results are resistant to phase variations in the functions.  We establish the universal approximation property of our non-linear decoder and validate the effectiveness of our model through extensive experiments.

Shape-Informed Clustering of Multi-Dimeensinall Functional Data via Deep Functional Autoencoders
Samuel Singh
Shirley Coyle
School of Computer Science and Statistics Trinity College Dublin
School of Electronic Engineering
Dublin City University
Mimi Zhang
School of Computer Science and Statistics
Trinity College Dublin
Abstract
We introduce FAEclust, novel functional autoencoder framework for cluster analysis of multi celtiiddmeensonall functional data, data that are random realizations of vector-valued random functions. Our framework features a universal- -pprroxiappo encoder that captures complex nonlinear interdependencies among component func- decoder capable of accurately reconstructing
tions, and a aniversal-approxing both Euclidean and manifold-valued functional data. Stability and robustness are enhanced through innovative regularization strategies applied functional weights and biases. Additionally, we incorporate clustering loss into the network's 's train- ing objective, promoting the learning of latent representations that are conducive to effective clustering. A key innovation is our shape-informed clustering objective, ensuring that the clustering results are resistant to phase variations in the functions. We establish the universal approximation property of our non-linear decoder and validate the effectiveness of our model through extensive experiments.
1 Introduction
Advancements in information technology have enabled collection of multi-dimensional functional data (FD) across various fields. Analyzing multi-dimensional FD is challenging due to its complexity and the intricate relationships between variables that change over time and space. The large volume of FD, especially in fields like medical imaging and climate science, further complicates efficient data processing. Zhang and Parnell [37] recently conducted a comprehensive review clustering methods for FD. The review highlights that current approaches to clustering multi-dimensional FD typically follow one of two strategies: (1) performing multivariate functional principal component analysis S and subsequently applying traditional multivariate clustering techniques to the functional principal component score vectors, or (2) defining a (dis) similariyy measure for multi diltii-dimessonall FD and then using a (dis)similarity- based clustering method. Although few other works apply alternative basis systems rather than the eigen-functions of the covariance operator, they are still limited to learning linear representations of FD. Given that the relationships between dimensions can be nonlinear and complex, linear methods, including those derived from multivariate functional principal component analysis, often fall short in effectively learning FD.
Neural networks, known for their ability to model nonlinear relationships, offer a compelling alter- native for multi-dimensional FD analysis. Rossi et al. (33 was among the first to apply multilayer perceptrons (MLPs) to FD by introducing functional weights in the first hidden layer. They demon- strated that functional MLPs approximate continuous mappings from compact subset of functional space to with arbitrary precision. In the context of cenettionoonfuucttii regression,
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Wang et al. [3] converted FD into vectorial form using multivariate functional principal component analysis, subsequently building traditional neural network for regression where both inputs and outputs were multivariate (i.e., functional principal component score vectors). Unlike [23 which modeled functional weights by cubic B-splines, Yao et al. 355 approximated each functional weight (namely their discrete evaluations over a grid) by neural network. Thind et al.
[22] extended
the functional MLP framework to handle inputs that include both functional and multivariate data. Heinrichs et al. [ defined each functional neuron to be convolution (w * y), where the functional weight w is a translatooninnvaiian
kernel: w(s,t) t) = w(s - f). A few recent works have developed neural network architectures for operator learning, namely learning mappings between two function spaces [6 6 24 Notably, all these studies focus on supervised learning, designing neural networks for tasks like regression classification and training neural networks with labeled data. As result, they are not suited for cluster analysis, which requires identifying patterns in unlabeled data.
For cluster analysis of FD, a logical progression is to adapt the autoencoder architecture to FD. Hsieh et al. [00 made strides in this direction by developing functional autoencoder (FAE), where the functional weights are modeled as integral kernels. However, this approach poses substantial computational challenges, as it requires training a separate integral kernel for each connection between neurons, making the training process exceedingly complex. Seidman et al. [25] applied variational autoencoders to FD. However, their approach has a critical limitation: the inputs (and outputs) are not continuous functions but rather function evaluations over fixed grid. Therefore, the model is not ceaeeizzation-nvariinnt and new architectures with new parameters may be needed to achieve the same error for data with varying discretization.
In this work, we introduce FAEclust, Ia comprehensive Python framework for clustering multi- dimensional FD. FAEclust is capable of handling FD in both linear spaces and Riemannian mani- folds, combining computational efficiency with robustness to phase variation (time warping) and high stability. Moreover, we prove that the functional decoder is a universal approximator. The paper is organized as follows: objective, and Section E against statc-of-the-an methods, and Section discusses implications and future directions.
Section details the FAEclust architecture, Section 3 introduces the training presents the clustering methodology. Section S benchmarks FAEclust
2 Functional autoencoder
In FD analysis, subjects are represented by random functions rather than traditional (multivariate) random variables. Below, we introduce the definitions and assumptions used throughout this work.
Let (y: = 1 1....,77 n) be a set of n independent p-dimensional sample functions: (t) = (t)..... v (t)) E TRP for any & TT where 'T is a compact interval of R; the superscript T is the transpose operator. Each sample function yı is a random realization of a p-dimensional random function y = -Y'..... yP.T. For d = ...... P, the n miedimensional sample functions (v/...., y are independent realizations of the component random function yd yd(t) = (t, defined on a probability space (S2, J, Pr) and taking values in H(T, ]R.. Here, H(T, R) is the separable Hilbert space of all square-integrable measurable functions that are defined on T and taking values in R. That is, yd is measurable map from l to H(T, R). Alternatively. we can view the function value yf (4)
realization of the random variable yd(,,-), mapping from (S, to (R, BR), where BR is the Borel r-algebra of R. In the following, we will suppress the dependence of a random function Y y(,ww) on Ω and simply write Y. Appendix A provides additional preliminaries on multi-dimensiona FD.
To formally define the FD clustering problem, assume there are K(> 2) clusters in the population. If a sample function y; belongs to the kth cluster (1 < k S K), then it is a realization of the kth random function YE, defined on the probability space (Ω, Prx). In other words, there are K different probability measures defined on the r-algebra (22, Given access only the discrete evaluations (ta)............. (tir.) } 1. a clustering method is to identify the underlying random function for each sample function y,. Table 4 in Appendix A summarizes the notations we will use throughout the work. All vectors are column vectors.
2

---

2.1 Network architecture
FAE consists of an encoder & and a decoder D: the encoder & maps a sample function y H(T, into latent representation E R,, and the decoder D maps the representation to reconstruction of y. denoted by E H(T IR.). That is, the encoder a mapping from the function space to latent finite-dimensional space & H(T, IRP) + R and the decoder is mapping from the latent space back to the function space D IR + H(T RP). Cluster analysis performed on the embedded data in the latent space. Figure presents a vanilla FAE architecture,
where the functional weights [wit) 4 01 S d S p) and { q.d
d S p, S S @ are continuous functions, distinguishing them from typical scalar weights. Each node in the input
encoder
decoder
Figure An illustrative FAE architecture, where we have five hidden layers.
layer accepts function in H(T, RR, and in the output layer delivers function in H(T, R). The FAE architecture depicted Figure can be succinctly described by the following flow of information:
(1) EH(T.R)
(1) (X((T
Here, MLP stands for Multi-Layer Perceptron. This representation clearly demonstrates that the flow information through the hidden layers, namely (1) MLP MLP essentially forms classical MLP autoencoder. Define weight matrices of appropriate dimensions:
(() [w ((() E RRIPP, and (())()) (1) e R*x
The feedforward equations for the first hidden layer and the output layer are given by:
x (2) al (t) (t) (t)
(1)
where b is a bias vector, and is non-linear activation function. It is important to note that for the output layer, the activation function is linear, and there no bias vector.
The output layer in Eq. T involves only matrix multiplication, which is inherently linear operation. Therefore, the decoder lacks the capacity to learn nonlinear mappings from a Euclidean space to function space. The reconstruction error Σp= lly - gdd 1 is fundamentally constrained by how well the FD can represented finite-dimensiona linear subspace. When the FD lie on low-dimensional submanifold, the FAE architecture requires a nonlinear decoder for efficient reconstruction. To incorporate nonlinearity, we extend the FAE architecture depicted Tigure[bb introducing additional hidden layers, where each hidden layer is composition of linear operations and nonlinear activation function. For example, with two more hidden layers, the information flow is follows:
( EH(T.R)
& (1) W EH(T.R)
Wa.e EH(T (R)
EN((TB) (4)
x (1)
(()
(2)
(2)
(3)

---

Note that the weights for the additional hidden layers are all continuous functions. The corresponding feedforward equations for the additional hidden layers and the output layer are:
(1) = ((W(()E (b) y g(a)(t) = (+)
(2)
(3)
(t) = (tg
where b1 (t) and (t) are bias functions. Again, for the output layer, the activation function is linear, and there is no bias function.
While the encoder & is generally not an injective mapping, Stinchcombe 29 proved that & acts as a universal approximator when the activation function a is continuous and non-polynomial. We establish below that the decoder structured as above is also universal approximator in the Euclidean setting and achieves patchwise universal approximation when the output lies on a compact, connected Riemannian manifold. Theorem D ensures that the decoder can accurately reconstruct any latent representation, allowing the encoder to fully utilize the latent space without avoiding regions where the decoder might otherwise fail. The construction of the readout map p is provided in Appendix B.2
Theorem 1. Fet : HM) denote a family of continuous mappings. When M is a P 1-ddmenniooall Euclidean space. for any f F and e E (0. 1), there exists functionall network D R' H(T, RP) with the structure given by Eqs. (2-4) such that supzeE E (D(x) - f fz))ll f any compact set F. C TK. When M is compact, connected Riemannian manifold that isometrically embeds in TRP, for any f E F and E (0, 1), there exists a functional network DD TR -> H(T, RR) with the structure given by Eqs. (2-4) and readout map p. such that sup E
Dlpc f()) l| <
for controlled compact sets E C R whose maximal diameter depends on the curvature of M.
2.2 The joint training and clustering framework
Given the unsupervised nature of our problem, network training and cluster analysis must be per- formed jointly in each Nowaard-backward loop. Specifically, during the forward phase, we update the learned latent representations {E()) ) which necessitates concurrent update of the clustering results. In the backward phase, we update the network parameters by minimizing unified objective function that incorporates both the network training objective and the clustering regularizer.
Our methodological framework is depicted in Figure 2 where network update and cluster up- date carried out jointly in an iterative manner. In the for- ward phase, clusters are updated by minimizing a clustering ob- jective function L.; in the back- ward phase, network parameters are updated by minimizing an in- tegrated objective function L = L, + C. + XeLe. Here, L, is the reconstruction loss, L. is regularization term on the func- tional weights and biases, and L. is measure of clustering va- lidity. The tuning parameter Aw controls the amount of regular- ization. The combined term L, + Ar w. is interpreted penalized reconstruction loss. Incorporat- ing the clustering loss Le into the training objective encourages the encoder to learn representations that are conducive to effective clustering.
Network Updare
Clustering Loss Penalized Reconstruction Loss
Input
Reconstraction
MM
MWM
Encoder
Bacoder To
Fintte-iiensionall
Cluster Update
Figure 2: The joint network training and clustering framework.
an initial clustering of the FD available, then training can performed through backpropagation in the traditional way. Otherwise, if no initial clustering is available, we need to pre-train the network by minimizing only the penalized reconstruction loss L, + Lw After obtaining a partition of the embedded data, then fine-tune the network by minimizing the integrated loss function L.

---

Backward: network update
3.1 Clustering loss Le
To encourage the FAE to learn specific loss, denoted by Le, into the network training objective function. Let the generic notation X E IR denote the data matrix in the embedding space: X =
representations, we need to incorporate a clustering-
.,,,, and let { ) denote the current partition of the data X, where C2 is the kth cluster with ny members. The fundamental concepts of clustering validity are compactness and separation. A natural measure of compactness is the within-cluster variance, and a natural measure of separation the total distance from the cluster centroids to the overall sample average. To encourage both compactness within clusters and separation between clusters, we formulate the clustering loss penalty to be:
(x
E ZIECL
|
Ms
Σ
Th | a
(2 Σ
||=
1x
|| x 12
B.EC,
where μı Ea,ec, 3; is the centroid the cluster Cg. and 2 = Ei= H; is the overall sample average. The formulation of the clustering loss function index, clustering validation measure. The main difference is that we additionally include the dimension parameter 8, to make the partition [C= comparable over different embedding spaces.
accordance with the
3.2 Penalized reconstruction loss L, + AwLw
The network parameters are weights and biases and are learned through minimizing a training objective function in the backpropagation manner. Let denote the collection of all network parameters. The reconstruction loss of the n input-output pairs (yi 9, D is
1,, 1 - -" 1.
(,)
71
In the relevant works reviewed before, the network parameters are learned by minimizing the reconstruction loss (0) only. However, for each functional dimension d, without appropriately regularizing the functional weights (w
g.d S q S 91), the training process might learn similar functional weights and hence extract redundant information from the functional input y² Likewise, the absence of constraints on the functional weights {w S S )1) can render the training process unstable and prone to overfitting. We therefore introduce a few regularization terms into the reconstruction loss function to improve the stability, robustness, and generalization performance of the trained model.
Orthogonality Penalty To encourage different functional weights in the encoder to learn different (uncorrelated) information about the input function, we can regularize them be orthogonal. Note that we only require the within-component functional weights (w 1 I q S 91 } to be orthogonal, while the choice of functional weights for one component random function independent of that for any other component random function. To encourage orthogonality, natural choice of penalty is the l -tyype regularization:
En However, the above functional is neither differentiable nor convex. In fact, even in vector calculus, an absolute inner product (W1, W2)l is neither differentiable nor a convex function of the argument pair (W1, W2). Therefore, computing a functional variant of "subgradient" (to perform gradient descent) requires heavy machinery of the notion variation the calculus of variations. We are looking for exactly orthogonal functional weights, and hence we here adopt the f2-type regularization (rather than the & -type regularization):
d=1 E (w (1) 9.ds 9.d/ Rel +
Σ 9=1
w C2 1[].
L (w =
E Ei
(w w q.d:
ge + Σ.r
21
(1) g.d
1)
d=1
To perform gradient backpropagation and directly optimize over the functional weights, definition of functional derivative is required. A natural choice is the Fréchet derivative in the calculus of
5

---

variations 22 Chapter 4]. However, implementing optimization over functional weight w in computer code entails optimizing the function's 's values over fine grid, leading to significant computational overhead. Furthermore, while enforcing the orthogonal penalty Lw L. (W(1)) can aid in regularization, the optimized functional weights could exhibit violent fluctuations. A natural expectation for the functional weights (w/td 1 S S Q1, - d S p} is their continuity and differentiability. Therefore, to uphold the continuity and differentiability of the functional weights (word S q 41, s d S p} throughout the training process, we formulate the functional weights into linear combinations of predetermined (continuously differentiable) basis functions. This reformulation reduces the optimization problem from one over the functional space to one over the coefficient vectors of the basis expansion. This approach not only preserves the smoothness and differentiability of the functional weights but also significantly reduces computational complexity.
q.d
Roughness Penalty The functional weights (and functional biases) in the decoder are for approxi- mating the FD, not for learning latent representations. Therefore, to prevent the functional weights and biases from exhibiting sharp bends or excessive curvature, we introduce a second-order roughness penalty, namely,
(w(())== - (w (wh)) dt,
to penalize functional weights and functional biases with high curvature. Again, to maintain the continuity and differentiability of the functional weights and functional biases throughout the training process, we formulate them into linear combinations of predetermined (continuously differentiable) basis functions. Consequently, in the code implementation, the roughness penalty on, e.g., the functional weight wa.q is effectively replaced by the & penalty on its basis-expansion coefficient vector. Here, we adopt the L regularization to encourages sparsity in the coefficient vectors.
Incorporating all the regularization terms on the functional parameters, the penalized reconstruction loss function is given by:
(0) + Am [Cw (W(() + C. + (ubl),
where the penalty terms on the weights and biases are functions of their basis-expansion coefficient vectors. In other words, by representing the functional weights and functional biases linear expansions of basis functions, all the trainable parameters in the FAE become scalar values, allowing us apply standard scalar backpropagation techniques.
The scalar weights in the two MLPs are regularized using batch normalization and dropout. Details of these techniques, as implemented in our Python framework FAEclust, provided in Appendix The optimization algorithm used to train the FAE is described in Appendix D
Forward: cluster update
4.1 Clustering objective function
Once the latent representations are available, we then apply a clustering method on the embedded data to obtain clusters. We want the clustering method to be able to adapt to the original FD, rather than completely ignoring the original FD. One appropriate choice is the technique of convex clustering [36,, with the clustering objective function being:
(U. x) 1 || - Σ 1sisJsn s(y. B,) uylle,
(5)
where || || is the Frobenius norm, and Uy): ) is similarity/affiiity measure hetween the sample functions y: and zj- Here, UT / un,, and E R (i (i= 1,... n) is the centroid of the cluster that object z; belongs to. Therefore, by forcing u, Uy. the two data points z; and Φ, will belong to the same cluster. The second term in Eq. (5) is regularizer, putting a constraint on the number of distinct cluster centroids. If X O, the minimum is attained when point occupies a unique cluster u; As À increases, the cluster centroids begin to coalesce. For sufficiently large X, all related points will coalesce into single cluster. Section 4.2 we develop an efficient algorithm for minimizing the loss function C, (U, X) w.r.t. U, for any value of X(> 0).
X, and each

---

For traditional multivariate data, the similarity measure (-.-)) is usually defined to be distance- dependent. e.g. (() exp(- Iz z; ()llea), order to make the estimates of the centroids enjoy asymptotic consistency. Let al(.;)) denote a distance metric, and Nm denote the set of m-nearest "neighbors" ye. where the neighbors are determined by the distance metric cl(.,)) -) In analogy to [36,, the similarity measure (W).V,) v2) can be formulated as
syy: U3) 5(y E (y or yı E (w))) x exp(-al(s:9))).
The introduction of the neighbor set N, (y:) aims to preserve the locality of the clusters by pushing nearby data points together. In our Python framework FAEclust, offer two options for determining the optimal neighborhood size m: the distance knee method and graph connectivity analysis.
Unlike traditional multivariate data, FD may differ in two types of variation: amplitude variation in curve height and phase variation in lateral displacements of curve features (e.g. peaks, points of inflection, and threshold crossings). If no phase variation is presented in the sample functions, or if it the joint variation between amplitude and phase that determines the clusters, then al is the standard
distance metric for Hilbert spaces: d yj) V 1l 254 l If amplitude variation the
main focus, with phase variation being a nuisance, then we can utilize the square-root velocity (SRV) framework [27]. In particular, we require that any sample function y: absolutely continuous on T; that is, the component functions [ d== are all absolutely continuous. Let H = {h T T denote the set of all orientation preserving diffcomorphisms. Then for any h E H, the composition y, h is parameterizzaion of y:- The SRV representation of is
SRV (() - {
(1 #0;
0,
otherwise.
Then the SRV representation of the re- re-parameterized function y; o h is
SRV o (t) (RV())(((t))
The distance between y: and yj is defined be dly. vı) inf d o h,y,) inf SRV o h) SRV SRV keH
(6)
hell
where FR is the Fisher-Rao Riemannian distance metric.
The distance metric defined in Eq. (6) has the property that it is invariant to both translation and reparameterization. Its computation relies on dynamic programming, which, however, has a time compiexity of O(N2), where N is the number of nodes of grid on the interval T [26, Appendix BJ. Therefore, calculating the pairwise distances (yl By) : S < j S n) requires substantial computational resources for datasets of medium to large size. To address this, our Python framework FAEclust additionally implements fast and ultra-fast dynamic time warping (DTW) techniques to approximate the similarity s (gy, y,). The fast DTW algorithm achieves a time complexity of O(N),
while the ultra-fast DTW [21] further reduces complexity to -O(N), enabling scalable and rapid computation for large datasets.
4.2 Optimization algorithm
The objective function C. (U, X) is separable on dimensions, and hence the minimization can be carried out separately in parallel for each embedding dimension. We might let the generic vector x represent an arbitrary column X, and u the corresponding column of U. Minimizing (U X) amounts to solving minimization problems of the following form: u(>) = arg mill Li (u), where
(u) } | ulli, +
(,, Σ s
d4
We here develop a path-folli hollowollowing homotopy algorithm with a complexity of O(n log(n)) for finding the solution (\)) for any value of 0, utilizing the property that each element u. (A( is piecewise
linear function the parameter À 9
According to Eq. ( for given embedding X, the convex clustering objective admits a unique solution for any value of X. When O, the solution simply u(0) = æ. assigning each point its own cluster. As À increases, pairs of cluster centroids merge sequentially, giving rise to a complete
7

---

agglomerative clustering hierarchy. Crucially, these merges happen only at specific values of A, referred to as breakpoints. Let À take the value Ag. the breakpoint at which two clusters are being merged into one, and now the data are partitioned into K clusters; that is, there are K unique values in the solution u u(Ax). ), denoted by {ū.....; üx} Right after the merge, the algorithm computes the next breakpoint XK-1, at which the next merge will happen, producing a partition with (K - 1) clusters. When À lies between the two consecutive breakpoints, namely AK S À < AK-1, the clustering remains unchanged. Our path-following homotopy algorithm efficiently identifies all the breakpoints in (n log(n)) time, thereby constructing the full clustering hierarchy. The computational details of this procedure presented in Appendix E A high-level outline of the algorithm is provided below. Input: latent representations X and similarity measures (ss (s(. vy) S < S n). while number of clusters >> 2 do
Compute the next breakpoint Aσ where a pair of centroids merge. Update the clustering by merging the corresponding pair of clusters.
end while Output: breakpoints {1 }2) and corresponding hierarchical clustering.
Once the hierarchy is constructed, an internal validation index (such as the silhouette score or Davies- Bouldin index) is used select the optimal clustering result, thereby determining the number of clusters in a data-driven manner. This selection is entirely independent of any particular À value; that is, the À parameter in Eq. (5) functions as an algorithmic variable for path tracing, not as a tunable hyperparameter.
The network training process consists of several iterations of the forward pass, we repeat the hierarchy-construetie and clustering-selection procedure described above. The final clustering result from the last iteration is reported as the output of the model.
loop. In each
5 Experiments
We selected four FD clustering algorithms available on CRAN - funHDDC, funclust, FADPclust (FADP1 & FADP2) - for benchmarking because they natively support multi-dimensional functional datasets, where the input is a list of fd objects (one per dimension). Other FD clustering meth- ods available CRAN either do not accept fd objects as input (c.g.. kmeans __lggn)__ reduce multi-dimensional functions to a one-dimensional form, e.g., by concatenating the basis-cxpansion coefficients. On the Python side, include the VANO model (25), whose reference implementation available on GitHub. Additionally, we re- re-implemented the methods from g (FNN) and [0] (FAE) ourselves no official implementations are currently available. To ensure consistency, our FNN uses the architecture depicted in Figure of [8, and our FAE follows Algorithm I of [TO]. In Appendix F.1 detail the input arguments of our FAEclust algorithm and summarize the configuration settings used for every algorithm in this benchmark.
When the FD are in Euclidean space, applied all eight FD clustering algorithms to nine one- dimensional and eight muuttiddimessionll functional datasets from the UEA & UCR Time Series Classification Repository. Clustering performance, evaluated using the Adjusted Mutual Information (AMI), is summarized in Table D. Results on the Adjusted Rand Index (ARI) are provided in Appendix H.3 Dataset abbreviations used in the table include: DSR (DiatomSizeReduuttion) BM (BasicMotions), SWJ (SkkndWalJJump), EOS (EyesOpenShut), FM (FingerMovemmntts) and JV (Japanese Vowels). FAEclust is the most consistent performer: it achieves the best AMI on 12 of 17 datasets and ranks top-2 on 15 of 17. Table | highlights the robustness and generalizability of FAEclust in clustering both univariate and multivariate FD.
For manifold-valued FD, we investigate five types of manifolds: Hypersphere, Hyperbolic, Swiss roll, Lorenz, and Pendulum. Details on the simulated FD are given in Appendix F22 For each simulation scenario, we apply each model on 100 simulated datasets and report the mean and standard deviation of AMI (Table 2 and ARI (Table 77. Figure T0 in Appendix of clusters identified across the 100 repetitions. FAEclust attains the best performance on all five manifolds, with especially strong gains nonlinear dynamics (e.g., Hyperbolic and Lorenz) and near-perfest accuracy on Pendulum, while exhibiting low variability across repetitions.
F3 shows box plots of the number
We evaluate the robustness of FAEclust to phase variation by extending the 12 simulation scenarios from [I. For each scenario, we begin by generating an unwarped functional dataset and applying FAEclust. We then introduce phase variation by composing each function with randomly gener-
8

---

Table 1: AMI scores for the 17 Euclidean functional datasets.
Dataset
funHDDC
funclust
FADPI
FADP2
FNN
FAE
VANO
FAEclust
BirdChicken
0.055
0.080
0.019
0.055
0.290
0.259
0.302
0.339
CBF
0.362
0,066
0,435
0.318
0.363
0.370
0.743
0.724
Chinatown
0.147
0.147
0.246
0.055
0.236
0.110
0.339
0.343
DSR
0.786
0.017
0.692
0.876
0.679
0.645
0.713
0.887
ECG200
0.173
0.101
0.201
0.143
0.214
0.213
0.217
0.261
Fungi
0.773
0.186
0.477
0.501
0.244
0.624
0.798
0.925
Plane
0.819
0.013
0.725
0.741
0,841
0.825
0.846
0.907
Rock
0.373
0.228
0.216
0.335
0.184
0.089
0.355
0.447
Symbols
0.767
0.001
0.435
0.712
0.748
0.800
0.817
0.824
Blink
0.410
0.048
0.189
0.177
0.453
0.506
0.522
0.633
BM
0.377
0.031
0.191
0.422
0.401
0.676
0,592
0.539
EOS
0.209
0.018
0.104
0.102
0.152
0.206
0.180
0.266
Epilepsy
0.143
0,028
0.077
0.225
0,274
0.209
0.297
0.485
ERing
0.735
0.012
0.288
0.714
0.643
0.743
0.733
0.664
FM
0.001
0.001
0.002
0.002
0.174
0.138
0.174
0.228
IV
0.840
0.069
0.294
0.466
0,236
0.854
0.899
0.893
SWJ
0.268
0.040
0.248
0.046
0.174
0.170
0.344
0.324
Table 2: AMI scores for the five manifold-valued functional datasets. The table reports the mean (top row) and standard deviation (bottom row) of the scores over 100 repetitions.
Dataset
funHDDC
funclust
FADPI
FADP2
FNN
FAE
VANO
FAEclust
Hypersphere
0.016
0.478
0.137
0.067
0.089
0.307
0.443
0.737
0.041
0.036
0.127
0.071
0.038
0,052
0.063
0.026
Hyperbolic
0.005
0.013
0.001
0.001
0.004
0.047
0.410
0.798
0.013
0.044
0.006
0.006
0.003
0.016
0.042
0.034
Swiss roll
0.127
0.114
0.382
0.189
0.016
0.125
0.242
0.432
0.059
0.035
0.046
0.065
0.026
0.035
0.062
0.040
Lorenz
0.109
0.389
0.092
0.144
0.023
0.246
0,251
0.457
0.057
0.049
0.060
0.081
0.047
0.022
0.034
0.038
Pendulum
0.887
0.376
0.797
0.808
0.253
0.794
0.905
0.986
0.029
0.058
0.049
0.063
0.044
0.074
0.038
0.006
ated time-warping function, created using the Lagmeeter wtteeeerrtimeeeeeeee time_waar module in Python. FAEclust is subsequently reapplied to the time- warped dataset. This entire procedure repeated 100 times, and (3plortreports the mean and standard deviation of ARI and AMI scores across 100 repetitions for both the original (unwarped) and time warped datasets. Figure [3 in Appendix [F4]]] shows box plots of the number of clusters identified across the 100 repetitions. Tables 8 and in Appendix F:4 summarize the AMI and ARI performance, respectively, of the seven baseline methods on the time-warped datasets.
The results in Tables and 9 show that FAEclust exhibits strong robustness to phase variation: introducing random time warping changes the mean AMI/ARI by only 0,013 on average. Per- formance remains essentially unchanged and near-perfect in structured settings, while the largest degradations are modest. Against seven baselines on the warped datasets, FAEclust attains the highest scores in 1/12 scenarios and does so by substantial margins; for example, relative the best baseline it improves by +0.10/+.. 199 (AMI/ARI in A, +0.3 +0.315/+0.3 315/+0.31 311 in B, +0.093/+0.1 121 in E, and +0.501/+0.465 (0.728 vs. 0.699), but FAEclust still yields the best ARI (0.709 vs. 0.667). Overall, these results demonstrate that FAEclust highly robust to phase variation and provides ofte-of-the-- clustering accuracy time-warped functional data.
K. The only exception is scenario L, where FNN attains slightly higher AMI

---

Table 3: AMI and ARI scores of FAEclust for the 12 simulation scenarios. For both the original (unwarped) and time-warped datasets, the table reports the mean (top row) and standard deviation (bottom row) of the scores over 100 repetitions.
Scenario
Original
Warped
Scenario
Original
Warped
AMI
ARI
AMI
ARI
AMI
ARI
AMI
ARI
A
0.782
0.730
0.754
0.723
0.582
0.484
0.580
0.492
0.026
0.024
0.020
0.022
0.106
0.110
0.098
0.116
B
0.428
0.441
0.412
0.385
1.000
1.000
1.000
1.000
0.060
0.085
0,074
0.092
0.000
0.000
0.000
0.000
C
0.820
0.751
0.782
0.710
I
0.997
0.993
0.996
0.992
0.042
0.057
0.036
0.045
0.007
0.004
0.009
0.006
D
0.581
0.514
0.556
0.518
J
0.
0.086
0.145
0.110
0.085
0.069
0.072
0.058
0.065
0.061
0.074
0.066
E
0.985
0.990
0.983
0.990
0.871
0.819
0.867
0.813
0.016
0.010
0.019
0.012
0.030
0.026
0.036
0.030
F
0.888
0,872
0.869
0.819
L
0.712
0.734
0.699
0.709
0.031
0,050
0.057
0.079
0.049
0.067
0.056
0.081
6
Conclusion
We introduced FAEclust, novel deep functional autoencoder framework for clustering multi- dimensional functional data. By combining universal universal-approximmai approximation guarantees for both encoder and decoder with chare clustering objective. FAEclust remains robust to phase variation and can capture complex, nonlinear (and even non-Euclidean) data geometries. We operationalize clustering through a convex objective and derive path-following homotopy algorithm that constructs the full clustering hierarchy in On log(n)), with model selection performed along the path via internal validation. Extensive benchmarking on real-world and simulated datasets demonstrated the clear advantages of FAEclust over existing methods. Overall, FAEclust sets a new standard for clustering complex functional data by integrating geometric structure, regularized functional representations, and a shape-aware clustering loss into a unified deep learning framework.
Acknowledgments
This work was conducted with the financial support of the Research Ireland Centre for Research Training in Digitally-Enhanced Reality (d-real) under Grant No. 18/CRT/6224. We are grateful the four anonymous reviewers for their valuable suggestions.
References
[1] E. Akeweje and M. Zhang. Learning mixtures of Gaussian processes through random projection. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings Machine Learning Research, pages 720-739. PMLR, 21-27 Jul 2024.
[2] A. Beck and M. Teboulle. A fast iterative problems. SIAM Journal on Imaging Sciences, 2(1):183-22) 2022 2009. doi: 10.1137/08071653 16542.
algorithm for linear inverse
[3] C. Bouveyron and J. Jacques. Model-based clustering of time series in group-specific functional subspaces. Advances in Data Analysis and Classification, 5(4):281 300, 2011.
[4] Chen and H. Chen. Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems. JEEE Transactions on Neural Networks, 6(4):911-91) 1995. doi: 10.1109/72.922553.
[5] J.-M. Chiou, Y.T. Chen, and Y.-F Yang. Multivariate functional principal component analysis: normalization approach. Statistica Sinica, 1714)1157 OCT 2014. doi: 10.5705/ss. 2013.305.
10

---

[6] G. Gupta, X. Xiao, and P. Bogdan. Multiwavelet-assed operator learning for differential equations. Advances in Neural Information Processing Systems, 29:24048 24062, 2021.
[7] W. H. Guss and R. Salakhutdinov. On universal approximation by neural networks with uniform guarantees on approximation of infinite dimensional maps, 2019. URL httphsts:///xxx org/ abs/1910. 01545) arViv.
[8] Heinrichs, M. Heim, and C. Weber. Functional neural networks: Shift invariant models for functional data with applications to EEG classification. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 12866-12282 12881. PMLR, 23-29 Jul 2023.
[9] H. Hoefling. A path algorithm for the fused lasso signal approximator. Journal of Computational and Graphical Statistics, 19(4):94--0066 2010. doi: 10.1 1998/cgs.20100092000
[10] T.-.. Hsieh, Y. Sun, S. Wang, and V. Honavar. Functional autoencoders for functional data representation learning. In Proceedings of the 2021 SIAM International Conference on Data Mining (SDM), pages 666-674, 2021, doi: 10. 10.1137/1.99 0177/197711111 077816611767777 1976700.75.
[11] loffe and C. Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In Proceedings of the 32nd International Conference on Machine Learning- Volume 37, page 448 456. JMLR.org, 2015.
[12] J. Jacques and C. Preda. Funclust: A curves clustering method using functional random variables density approximation. Neurocomputingg 12:1664 171, 2013.
[13] A. Kratsios and I. Bilokopytov. Non-euclidean universal approximation. In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 10635-10646. Curran Associates, Inc., 2020.
[14] A. Kratsios and Papon. Universal approximation theorems for differentiable geometric deep learning. Journal of Machine Learning Research, 23(196):1-73, 2022.
[15] M. Leshno, V. Y. Lin, A. Pinkus, and S. Schocken. Multilayer feedforward networks with a nonpolynomial activation function can approximate any function. Neural Networks, 6(6):
861-867, 1993. doi: 10. 1016/50893360080055)))00
[16] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anand- kumar. Fourier neural operator for parametric partial differential equations. ICLR 2021 9th International Conference on Learning Representations, 2021.
[17] E. N. Lorenz. Deterministic nonperiodic flow. Journal of Atmospheric Sciences, 20(2):130- 141. 1963. doi: 10.1175/152004400 0469 04656//9660000000000000222000..0
[18] I. Loshchilov and F. Hutter. Decoupled weight decay regularization. arXiv, 2017. doi: 10. 485550//RXIV.1111005000
[19] K. V. Mardia. Statistics of directional data. Journal of the Royal Statistical Society: Series
B (Methodological), 37(3):39--37, 1975. doi: tb01550.x.
[20] M. Nickel and D. Kiela. Poincaré embeddings for learning hierarchical representations. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
[21] T. Rakthanmanon, B. Campana, A. Mueen, G. Batista, B. Westover, Q. Zhu, J. Zakaria, and E. Keogh. Searching and mining trillions of time series subsequences under dynamic time warping. In Proceedings of the 18th ACM SIGKDD International Conference Knowledge Discovery and Data Mining, page 262-270 2012. doi: 10.1145/2339530.23397777
[22] D. H. Richard Courant. Methods of Mathematical Physics, volume 1. John Wiley & Sons, Ltd, 1989.
[23] F. Rossi, B. Conan-Guez, and F. Fleuret. Functional data analysis with multi layer perceptrons. In 2002 International Joint Conference on Neural Networks (ICNN), volume pages 2843-
2848 vol.3, 2002. doi: 10. 10.1109/ICNN.2002100000

---

[24] J. Seidman, G. Kissas, P. Perdikaris, and G. J. Pappas. Nomad: Nonlinear manifold decoders for operator learning. Advances in Neural Information Processing Systems, volume 35, pages 5601-5613. Curran Associates, Inc., 2022.
[25] J. H. Seidman, G. Kissas, G. Pappas. and P. Perdikaris. Variational autoencoding neural operators, 2023. URL https:t/arxii....... orgg abs/ abs/2302. 2302. 10351 arXiv.
[26] A. Srivastava and E. P. Klassen. Functional and Shape Data Analysis. Springer New York, NY, 2016.
[27] A. Srivastava, E. Klassen, S. H. Joshi, and I. H. Jermyn. Shape analysis of clastic curves in Euclidean spaces. IEEE Transactions on Pattern Analysis and Machine Intelligence, 33(7): 1415-1428, 2011.
doi: 10.1109/TPAM..200..11
[28] N. Srivastava, G. Hinton, Krizhevsky, Sutskever, and Salakhutdinov. Dropout: A simple way to prevent neural networks from overfitting. Journal of Machine Learning Research, 15 (56):1929-1195)) 1958, 2014.
[29] M. Stinchcombe. Neural network approximation continuous functionals and continu- functions compactifications. Neural Networks, 12(3):467-477, 1999. doi: 10.1016/ 50893-6080098301088-..
[30] S. H. Strogatz. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering (2nd ed.). CRC Press, 2015.
[31] B. Tenenbaum, V. de Silva, and J. Langford. A global geometric framework for nonlinear dimensionality reduction. Science, 290(5500):2319, 2000.
[32] B. Thind, K. Multani, and J. Cao. Deep learning with functional inputs. Journal of Computa- tional and Graphical Statistics, 0(0):1-10, 2022. doi: 10. 1080/ 2000122000002
[33] Q. Wang, H. Wang, C. Gupta, A. Rao, and H. Khorasgani. A non-linear function-on-funetio model for regression with time series data. In 2020 IEEE International Conference Big Data (Big Data), pages 232-239. IEEE Computer Society, 2020. doi: 10. 109/BigDData5002222220 9378087.
[34] X.-F. Wang and Y. Xu. Fast clustering using adaptive density peak detection. Statistical Methods in Medical Research, 26(6))2800-288 201228000211111- 2017.
[35] J. Yao, J. Mueller, and J.-L. Wang. Deep learning for functional data analysis with adaptive basis layers. In Proceedings of the 38th International Conference on Machine Learning, volume 139, pages 11898-11999 11908, 18-24 Jul 2021.
[36] M. Zhang. sorward-stagewisse clustering: An algorithm for convex clustering. Pattern Recogni- tion Letters, 128:283-2888 2899, 2019. doi: 10. 1a..006//..atree2000
[37] M. Zhang and A. Parnell. Review clustering methods for functional data, ACM Transactions on Knowledge Discovery from Data, 17(7):34, 2023. doi: 10. 1053849/3581789.
12

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: All the contributions and claims are properly explained in the introduction and abstract.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
A discussion of the limitations of the current model is provided in the Ap-
pendix.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiitiinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
13

---

Answer: [Yes]
Justification: Assumptions for Theorem are clearly stated in the statement. Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
The linked repository includes all necessary scripts and detailed instructions to replicate the experimental setup and results. It is well-structured and clearly documented, enabling users to reproduce both the simulated and real-data experiments. This ensures that the main claims and conclusions of the paper are independently verifiable.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data provided or not.
the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access a hosted model (e.g.. in the of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
14

---

5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
The real datasets sourced from the UCR Time Series Archive, while the simulated datasets can be reproduced using the provided scripts. Detailed experimental settings are available in the Appendix and the associated Python package.
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
//nips.cc/ lttipss.ccppublicuuu.eeee
for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The training and testing procedures thoroughly documented both in the paper and in the accompanying Python package documentation. Users can directly run the provided scripts to reproduce the experimental results reported in this study.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, or as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer:
[Yes]
Justification: Error bars (standard deviations) are reported Appendix F
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emoor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
15

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative еπτог rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information the com- puter resources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Compute resources are reported in Appendix F
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics :tttp:://neuri: wttphhewwwss..ccpppoeeeeeeeeeeee
Answer: [Yes]
Justification: This paper has no foreseeable ethical issues.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction..
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification: This paper presents work whose goal is advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
16

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
We have properly cited all existing assets: code, data and prior works.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
17

---

If assets are released, the license, copyright information, and terms of use in the thcode. has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, batasets.papersstt paperswi
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors encouraged to reach out the asset's s creators.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: The developed Python package is well documented and will be uploaded to GitHub after the reviewing process.
Guidelines:
The answer NA means that the paper does not release assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
18

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in this research does involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM for what should or should not be described.
19

---

A Preliminaries
While the n sample functions (y: :1 L....... ......n] are called functional data, in real applications, we only have access to their discrete evaluations. Moreover, the observation of y, at any time point t may come with an additive error. Hence, we write (;(t) (t()) + <((,, where (.( is the noise process with (t]] O, ((e] = oz and (t)
for & + t and S d, S p. The sampling scheme for the ith subject is denoted by the column vector t; (4a1,..., lin. where laj E T for = ....., ri, and 4/1 <... < lir.. Here, to avoid overly complex notation, we assume that all component random functions have an identical function domain T, and the p component functions {v/ ... v,') of each subject have the same sampling scheme. The sequence of observations {út; (4:1)...... We (fir,) } called a sample path the sample function y,. The sample path of mie-dmmnsioona sample function y written as V (t;). where we employ compact notation for functions applied to collections of input points.
If sample function y; belongs the kth cluster S k S K), then it is realization of the kth random function Yg, defined on the probability space F,, Prę). In other words, there are K' different probability measures defined on the o-algebra F).. Given the sample paths (ta).... ye (tar, ))=1- clustering method is to identify the underlying random function for each sample function y. Table 4 summarizes the notations we will use throughout the work. All vectors
Table 4: Notations adopted throughout the paper.
number of subjects
i
index for subject, 1 S. i S number of component random functions index for dimension, S d S number of clusters in the population index for cluster, < k < K
p
d
K
k
Ck
kth cluster of size n2: Σ =k=1 Rig
[yi Wi:i= ......n]
n independent p-dimensional sample functions compact function domain E T
(1())
white noise with variance oî for = 1, .,....P observation of (1): (1) y. (4) (L))
(1(()
t.
sampling scheme for the ith subject
number of sampling points in the ith sampling scheme sample path of the p-dimensional sample function yı sample path of the component sample function vi kth random function Yú = (Y (t, w).... (t (t, probability space for YE, inner product and norm for the Hilbert space H(T, RP) & and l2 norm on an Euclidean space first-order derivative of the function y() 0-1 indicator function signum function
..... ... ((r))
(t,)
Ys
(2, Præ)
1,...,,k K
H
16/
y()
8(.)
sgn
are column vectors.
The network model developed in Section 2 only accepts input continuous functions {y/ 1,..., n.. In general, however, we do not have access to any sample function yi, but only its noise- contaminated discrete evaluations (t.))...... (f..r)))) Therefore, in practice, need to perform smoothing before cluster analysis.
Smoothing Given the sample path (ti)......... ........... (tir.), }}. need to recover the underlying continuous function (t), which is commonly achieved through the smoothing technique. Let (b,, dm,. be m element functions of preedeter basis (e.g., the cubic B-spline basis) defined on the interval T. Then we can approximate each component sample function v by Σ be. where the coefficient vector c¹ (ch.... ...,.................. & TT is commonly obtained through minimizing
20

---

regularized residual sum of squares:
min CR
r=1 (tirr) - m
C (tr)) regularizer).
The choice of the smoothing technique depends on fields of application: wavelets are effective in capturing discrete jumps or edges and especially useful for modeling signals and images; splines are more often applied for approximating a function of given order. Kernel smoothing is another common method of estimating the mean function in the regression model 07 ()) (t) + (t). we formulate our model in terms of the sample functions {y/ with the implication that they will be replaced by their smooth estimates in real applications.
From Section 2
,....,,)) n) only,
Standardization As with the analysis of traditional multivariate data, where a proccessssing step is to standardize the data to have unit standard deviation in each dimension, Chiou et al. [5] proposed to standardize each component sample function to account for differences in degrees of variability and in units of measurements among the component random functions. The formula is in analogy to standardizing multivariate data: each component sample function subtracted by its mean function and then divided by the square root of the variance function. Standardization is necessary if the component random functions have quite different ranges or if they exhibit different amounts of variation. If the sample functions
scaled, algorithms that compute the overall
distance Ei= 1o0f 11 are biased towards component random functions with numerically larger values. Therefore, standardization a crucial step for distance- based clustering algorithms. Moreover, standardization helps machine learning algorithms train and converge faster. Henceforth, we assume that standardization has been done before applying the network model.
B Proof of Theorem I
B.1 Euclidean Universal Approximation
Existing universal approximation theorems primarily focus on mappings between finite-dimensional topological vector spaces. While some works have explored the universal approximation property for nonlinear functionals (29 and nonlinear operators 4 only single study has established universality for mappings from a finite-dimensional space to an dnfiniteddimensiona space (7]. Unlike previous approaches, which rely on integral kernels, our network architecture employs multiplicative parameters defined by continuous univariate functions. This key distinction necessitates a tailored universality proof for our specific framework, which we present here.
When M is the p-dimensional Euclidean space RP, it suffices only consider the case
1. This is because uniform convergence in H(T, IRP) reduces to uniform convergence in each component space H(T, RR). We now prove that neu- ral network ' with a single hidden layer can uniformly ap- proximate any continuous mapping f R² H(T, R) over any compact domain E' C RR; see Figure 3 Given that is continuous mapping from compact
( ad
P LL C (()
9
E' R²
into M(T, R). it is straightforward to prove that the range f(E) = (f(z) x E E} is a compact set H(TRR)
With the notation given in Figure 3 define (
(1),... (1) (t)T,, where is the number of nodes in
the hidden layer. The output function is
m
((x)(t)
Σ q=1
wi wi.g (f)a(w
(2) C (1)" (()
Figure 3: A functional network D IR + H(T, R) with only one hidden layer.
2 (t) aw (t)z b(t)),
where (1) = (w) (1)
w (T)) and (L) = [w (2) (I),... .... w (1). The activation function is chosen as Tauber- Wiener function (c.g., sigmoid) [4: that is, the set of all linear
21

---

combinations E cya(wnz + by) over a closed interval I C R is dense in C(T), where cg. Ug. and by are real numbers. The feedforward equation (7) differs from that of the traditional three-layer feedforward network only that the parameters are now functions, varying continuously with the value of the output function. One can immediately claim that, for any continuous mapping f RR* - H(T, IR) and any fixed / E T. there exist (constant) weight parameters {w (w()()), (I) (I))(()) such that E E (D)x)(() [()) < e for any compact set E C IR: see, e.g.. 15 This leads the following property:
Property Let (T, C H(T,R) be the set of all piecewise constant functions over the compact
interval T. For any continuous mapping Se R + (T, R., there exist piecewise constant weight
parameters (w) (t) (t), such that supzEE [D(xx) fa(xx)llx < 6 forr any compact set E C.
Given Property it remains to prove that the set of all piecewise constant functions CT, , R) dense in H(T RR. It well-known that the set of continuous functions C(T. IR. is dense in H(T, R). Hence, for any f(x)= E H(T, R), there exists a continuous function E C(T, IR) such that [x)( yllx < e/2. Since T is compact and y is continuous, y is uniformly continuous. Hence, there exists an h > O, such that for all &1, l2 E T with |61 12l < h, we have (ı((1) 22(2) < e/(2vITT) where [T] is the measure (length) of T. Let [T be partition of the domain T with max [ITy] : 1 1 j S M) < h. Define (t) = Σ == (tt)) E Tz), where ty is an arbitrary point in Tj. Then we have lyy fa( l < e/2, and it follows that [z) - fc(x)]l Since e >> 0 is arbitrary, the set of piecewise constant functions is dense H (T,
B.2 Non-Euclidean Universal Approximation
denote family of continuous mappings, where M is a complete connected Riemannian manifold with metric dm- We consider the class of feedforward neural networks Jm with J + layers and at most m neurons per layer, defined as:
F== [f
NN J.m (T) o0............... o T2o Ti R + H(T,
where T are affine maps of the form in Eq. 77.. Let NN = {NNN J.m : J 2, m 2 1} denote the class of neural networks of arbitrary depth and width. The notion of convergence on F is still that of uniform convergence on compact sets of R In analogy with Section B.1 we have the following property:
Property 2. Let (T, C H(T,M, be the set of all piecewise constant functions over the compact interval T. Under certain regularity conditions on the manifold M. for any con- tinuous mapping fc : R² + (T,M), there exist piecewise constant weight parameters = L....,JJ
such that SuPaeE lp o NN - ff) < for any e > 0. Here, p is a readout map that projects the output function from H(T RP) onto (((TM)); E C RR is controlled compact set whose maximal diameter depends on the curvature of M.
(E) (U)
Proof. If M is Cartan-Hadamare manifold (i.e., a simply connected, complete Riemannian man- ifold with non-positive sectional curvature), the exponential map Exp,, Ty.M + M is global diffcomorphism. and its inverse (the logarithmic map Exp, : + R²) is well- everywhere. This allows us to linearize manifold-valle Exp,y and process them using the neural network class NN. By Corollary 3.14 of [13], since Exp, is smooth and NN is dense in the space of continuous mappings from R to H(T, RR), the class (Exp, AN J.ml JJ 2 m 2 1} is universal in F. That is, for any compact set E C R and any continuous mapping F, there exists neural network NW J.m. such that
functions by mapping them into Euclidean space via
sup df Exp, o NN J. (z)(t))))) <,, #eE
(8)
for any T. The readout map Exp,, acts a fixed tnonntrainable) layer, projecting Euclidean outputs back onto M. The reference point y is typically chosen as the Fréchet mean of the data, minimizing intrinsic variance.
For an arbitrary complete connected Riemannian manifold M, the exponential map Exp, y only locally diffcomorphic near 0 E My If M is compact, it admits finite atlas (a finite collection of
22

---

manifolds, aligning with our reliance on smooth coordinate charts covering the manifold), In particular, we can cover it with finitely many logarithmic charts {(U,, Exp, ))== where each U; is a geodesic ball of radius ri S inj (the global injectivity radius) and Exp; is the locally defined logarithmic map on Ug. For any f E F. 1,.... .N. Because each logarithmic map
define E, (L) {ar E R (1) U.} for Exp;" is only valid the open subset U; M, the approximation (8) holds solely for compact subsets of (t)),, not arbitrary compact subsets of R². In particular, there exists network NNN J,m such that for any compact subset S (t))()),.
sup (f))()))t) Exp, J. (z) <c =EE
This reflects the intrinsic local nature of exponential charts general manifolds, contrasting with the global diffcomorphism property in the Cartan- Hartan-Hddamard case.
The proof Property provides an explicit construction for the readout map p. Specifically, since M is compact, there exists a finite set of anchor points {yI...., .yN C .M such that the geodesic balls B,, (y:) cover M, where each radius r; upper-bounded by the injectivity radius inj (3:). On each ball, the logarithmic map Expyy defines a smooth coordinate chart, and the transition maps Exp,, EXPy, are smooth diffeomorphisms overlaps Br, (B:) n Br, (Bj). The readout map is then constructed the collection of exponential maps {Exp,, 1==1 which project Euclidean outputs from H(T, IRP) back H(T.M) patchwise.
Our proof of Property 2 generalizes the framework of [13] to arbitrary compact, connected Rieman- nian manifolds, leveraging local exponential charts to construct patchwise universal approximator. Unlike the global diffeomorphism property of Cartan-Hadamard manifolds, our approach accounts for the intrinsic locality of general manifolds by partitioning the input space into regions where logarithmic maps are well-defined. Notably, Theorem 6 of [14) underscores the impossibility of global universal approximation for local coordinate systems.
C Scalar weight regularization
The scalar weights in the two MLPs are regularized through batch normalization and dropout tech- niques, either together or alternatively. We here briefly explain the dropout and batch normalization techniques; more information be found in the original publications
(28
The feedforward operation each hidden layer in Figure can be broken down into two layers: a fully connected layer (with the formula, e.g., _ = (z + b) and then an activation layer (with
the formula (1+-1) ("). With batch normalization only, the feedforward operation (over mini-batch) becomes
x (!)
W x (1)
+ b,
z 7 Norm (z) (0) + n.
(1)
x (++1) a( 2
where o denotes the element-wise product, and Norm ()) is the normalizing transform (by the batch mean and batch variance). With dropout only, the feedforward operation becomes
W (1) + b
x (1) a(x) (1)
1 Bernouli(r7)
x
(1+1) 1)
*p x
where 1 - is the dropout rate. If both the dropout and batch normalization techniques are imple- mented, eexplains the order of the different layers. Different from the batch normalization
CCmpaatnness or bounded geometry ensures inj(M) > 0, allowing finite atlases with geodesic balls of uniform size.
23

---

Linear
Activation Function
x
x ((++1)
W,
Linear
Batch Norm
Activation Function
Dropout
x1(
x((+1)
W, b
3
Figure 4: The order of the different layers: the fully connected layer is followed by batch normaliza- tion, then non-linear activation, and finally dropout.
technique, where y and r are trainable network parameters, the parameter in the dropout technique is a cyper-parameter
D Optimization algorithm: network training
Gradient descent optimizers with adaptive learning rates, e.g., Adam, have become a default method of choice for training feedforward and recurrent neural networks. However, adaptive gradient methods do not generalize well when a regularization term is added to the loss function. Loshchilov and Hutter [18] improved Adam by decaying the weight by little bit, rather than implementing the gradient of the penalized loss function. However, our objective of regularizing functional weights to encourage orthogonality, not shrinking their magnitudes. We therefore apply the technique of mini-batch gradient descent with momentum.
Write C y.); the loss function evaluated on the ith batch of the data, where 0 now may include the parameters from bath normalization. The algorithmic framework the gradient-descent technique is:
for 2, ... n_epochs do
randomly partition the data into mini-batches
for i 1.2.... n batches do
mi 3mo + (1 Ve 0o:yn.)
0,- am1
m. = m
0u 01
end for
end for
Here, > is the step size, and 0 < < is the momentum weight. The momentum starts at O. Note that we shuffle the data after every epoch. Within each epoch, after updating mo and 6o- we pick the next mini-batch, feed it to the FAE, and calculate the mean gradient of the mini-batch V L ).
Optimization algorithm: clustering
To obtain Ax-1, define the index set I₁ = {i (;(Ax) = üu iiliinn for & = 1...., K. We then re-formulate the objective function L, to be
(i(S)
K Σ Σx Rー
(si (xi-Cx12
+
EkolGa Gol.
24

---

where <
(G1....,Sxt is the vector of K optimization variables, and Skv Eiez, EJEz, B)). Let (()) arg mrg (i( denote the solution, again function of the parameter X. Apparently, when X = Ax. the vector i = (i1..... ux)² is the minimizer of the function L1 (C). and hence we have
0 =
ic, 05x
( ما
LEI xz) EVEI+.)++ AK
Σ
v=1 Sku x sgn (Üx üy);
where sgn(O)=O.. In fact, the above equation is valid for any value of À between the current breakpoint Ax and the next breakpoint -1- and therefore we can write
(A)
1
Ti
T 2
D v=1
&ko x sgn( ü.)l. AK S. < AK-1-
(9)
IE
We now can conclude that, when Ax S < A and i E Ig we have uz(A) = (x()), and that (A) is line segment, starting at üg and with the constant slope
8Ck 0A
1
K
Sku x sgn(üx- - -..)
(10)
2 I
==1
Equations (9) and (10) are valid for any value of À as long as the index sets {I+] remain unchanged. However, two things could happen at a breakpoint: two index sets get merged into one, or an index set I& splits into two. That is, if currently uy (axA))
u,(Ak) then could happen that
ua (A) F (A(A) for certain A(> )x). We assume that, at the breakpoints of the solution path, the index sets can never split
To determine the value of the next breakpoint Ax-1, for any pair of
(üge ú,), define
(Gn(A) " aCa(A)
ū.)) Ak.= (üg (ūx)=ū.)]
Or
Then &&v +A is the value for \ at which the two paths (A) and (a( will be merged together, assuming that no other merge occurs before that. A negative value of A indicates that the two paths (A) and (A) are actually moving apart for increasing À and hence will be ignored. Therefore, the next breakpoint is
AK- min 4a 0} AK-
(11)
If 412= min < 0}. then the minimizer u(Ax-1) of (u) for À = Ax-1 will have K 1 unique values, and the two clusters (a, } and {2, 10:28 i E I2} will be merged into one.
We have explained above how the solution (A)) evolves from the current breakpoint to the next. Now given any (warm) start X+ > O, let (X+)) denote the minimizer of L1 (u), obtained by any convex optimization algorithm. (In our code implementation, we adopt the FISTA algorithm developed by Beck and Teboulle (21.) Let the unique values in u uX++) be denoted by {ū) {......; üx). Exploiting the piecewise linearity of the evolving paths, we can quickly determine the next breakpoint: (1) determine
the index sets: I (i: uu(A+) üg, < i S n): (2) calculate the affinity weights: &
EiEL ΣjEz, syy.,y);; (3) calculate the derivatives: ACa. 7X
2/ Σx= $kv X ogn(ü, ü.).
"
Then the next breakpoint is Ax-1 min Abr:DA>>>>>>>
0} + X+ Apparently, trivial additional
computations, our math-following homotopy algorithm will produce a hierarchy, with the clusters being repeatedly merged.
F
Supplementary materials for the experiments
F11
Details on algorithm configuration
The same basis family and evaluation grid used for smoothing across all algorithms to maintain consistency. For all baseline methods, we run each method with the number of clusters ranging from 2 to 10 (extended to 20 for the Fungi dataset), ensuring that the true number of clusters is included in this range. We then report the best clustering performance based on ARI and AMI. In contrast, FAEclust determines the number of clusters in fully data-driven manner using the silhouette score within each forward phase of the joint training and clustering framework.
'If this assumption is valid, then our holl-following homotopy algorithm will produce the exact solution; otherwise, we will interpret our solution approximation. Splitting events infrequent.
25

---

1. funHDDC (from R package FunHDDC BD: All six supported models were evaluated: "AkjBkOkDK", "AkjBQkDk", "AKKB "AKBLOKDK", "ABkQRDK", "AkBQRDk", and "ABQkDk", Additional parameters included initialization via init "means". a convergence threshold of threshold 0. 1, model selection based on criterion "bic", and maximum of ittemax = 100 iterations. For each dataset, the final clustering result is given by the model that achieved the lowest BIC value.
2. funclust (from R package funclust 121::
The relevant arguments are set to thd = 05, increaseDimension = FALSE, hard = FALSE, fixedDimension integer (0), nbInit 20.
3. FADPclust (from R package FADPclust 34 34 The package consists of two clustering methods method "FADP1" " and method "FADP2". We evaluated both methods using the default parameters proportion NULL and cut 0.15.
4. FNN (18): The FNN architecture is (T) AO) (ATER((XX)) (22)
((1) ( G. (S, and (S, 7) are Legendre- basis convolution kernels. The first hidden layer is h (s)
(()
(2) x
+ (s)) where is ELU activation. Channel sizes 32 for U and 16 for U (() is Fourier- basis functional producing 16-dimensional vector (1)), The block x11) - x + (1)) is a 3-layer MLP, each layer with 16 neurons. The output layer is linear without bias: x(() (()())) Learned embeddings are
clustered by k-means and evaluated via AMI/RRI.
Settings: latent_din = 16, conv sonv_channels = c32,, 16), conv maav_saiis manse_aaiis "Fourier", activation "ELU", padding "same", optimizer "Adam", Ir = 1e-3, loss "MSE", , epochs 200,
"Legendre".
kmeans _ init 20.
5. FAE ([10):: The FAE architecture is y y. The functional layer uses tanh activation, and MLP layers use ELU activation: the output layer is linear, bias. Functional weights are represented via basis expansions (B-splines). Learned embeddings are clustered by k-means and evaluated via AMI/ARI.
w (MEN(T.R)
(())
2 (OE
Settings: basis "B-spline", basis_dia 10. fateet_dinn 16, layers 32, 16, activation "tanh", optimizer "ELU", Ir 1e-3, epochs 200, loss
= "MSE", kmeans irmeans__niitt = 20.
6. VANO (VAND 25]):
VANO is trained with Adam on reconstruction MSE plus 3 --wehtteed KL term. Layer widths are [64, 64, 16, 64, 64). Encoder means are clustered by k-means and evaluated with AMI/ARI. Settings:
latent_ dim 16, hidden_sizes c64, 64), activation "GELU", reparantterization = TRUE, optimizer = "Adam", Ir = 1e-3, beta = 1e-4. epochs 200, kmeans_n_aaaaeeeeeeeeaaaaaeaaa_m init 20.
7. FAFclust (FAEclust) FAEclust shares hyperparameters with standard autoencoders, including the number of training epochs (epochs), learning rate (Ir), batch size (batch_size), and dropout rate (tau). The network architecture is defined by the list layers, where the first entry and last three entries specify functional layers, and the middle entries form an MLP autoencoder. Functional weights and biases are parameterized using a basis family network, basis of size 1. The network's 's objective function includes two reg- ularization terms: an orthogonality penalty on the encoder's functional weights (lambda_e e) and sparsity penalty on the decoder's functional weights and biases (lambda_dd d). All hyperparameters epochs, 1r, batch_siza tau, layers, lambda_ _ee, and lambda_d -d are optimized via Bayesian optimization over predefined search space. The objective function in the Bayesian optimization problem is still the integrated objective function L C= -, + Cw X.C.. However, in the context of Bayesian optimization, the variables being optimized are the hyperparameters, not the model weights. In our implementation, we utilized the Optuna package to perform the hyperparameter tuning task.
In all experiments, both real and simulated, we fixed the FAEclust architecture to a seven- layer structure: one functional layer in the encoder, three functional layers in the decoder,
26

---

and three-layer MLP in between. The third layer in FAEclust serves as the bottleneck, producing the latent embedding. The maximum number of nodes allowed per layer was constrained to (64, 32, 16, 32, 64, 64, 64).
F.2 Details on the simulated manifold-valued functional datasets
Table 5 summarizes the number of samples, dimensions, time steps, and ground-truth clusters for each simulated functional dataset.
Table 5: Dataset parameters.
Dataset
# Samples
Dimensions
Time steps
# Clusters
Hypersphere
100
100
Hyperbolic
200
2
50
2
Swiss roll
300
200
Lorenz
100
100
Pendulum
200
100
Hypersphere [9]
Figure visualizes the clusters on a hypersphere S consisting of all unit vectors in Rt. We generate trajectories along great circles on the sphere, using different directions for the two clusters. The two clusters differ in period and are phase-shifted randomly, making this a challenging dataset for clustering algorithms that do not perform curve registration.
hypersphere-
Feature vs Time
hypersphere Feature vs Time
1.00
0.75
0.75
0.25
0.55
60
0.00
6.25
0.5
0.75
0.75
80
100
40
100
Time stop
Time stop
hypersphere Feature vs Time
1.00
0.75
0.25
00.55
0.50
0.75
1.00
Time stap
Figure 5: Hypersphere functional data. Trajectories on the surface of a hypersphere S with clusters defined by distinct great-cirele paths and phase shifts.
illustrates trajectories in hyperbolic space using the Poincaré ball model, which represents a space of constant negative curvature. We simulate geodesic motion within this ball. One class stays near the center, and the other ventures closer to the boundary.
Hyperbolic 220
Figure 6
Lorenz Figure shows time series generated from the Lorenz system, well-known chaotic system in three dimensions. For certain parameter values (sigma=10, 10, rho=28, beta=8/3), the solutions
27

---

hyperbolic Feature 1 vs Time
hyperbolic Feature vs Time
0.6
0.6
0.4
0.-
0.2
Value
-0.2
-0.2
0.4
004
0.66
-0.6
10
20
30
40
50
10
20
30
40
50
step
hyperbolic Trajectories features time)
40
30
0.5
0.5 Feature y
Feature
0.0
005
Figure 6: Hyperbolie functional data. Trajectories simulated in the Poincaré ball model of hyperbolic space, with clusters differentiated by proximity to the center or boundary.
approach the butterfly-shaped Lorenz attractor. We simulate three classes: non-chaotie class with r4o=14, and two chaotic classes with rho=21 and rho=28.
Swiss roll 31 Figure 8 depicts the Swiss roll, a classic example of 2D manifold embedded in 3D. We generate trajectories along the rolled surface, with different clusters corresponding to different vertical ranges.
Pendulum [30] Figure shows pendulum trajectories in cylindrical phase space (R x St), defined by angular position and angular velocity. We simulate both oscillatory motion (low energy) and full rotations (high energy), resulting in four distinct clusters.
F.3 ARI tables for the Euclidean and manifold-valuee functional datasets
Tables 6 and present the ARI scores for our benchmarks. Table covers the 17 Euclidean functional datasets (see Table I for the corresponding AMI results), while Table summarizes the ARI scores for the five manifold-based scenarios (cf'. Table
For the manifold-valued FD simulation, Figure 10 shows boxplots of the number of clusters identified over 100 repetitions for each manifold type. The boxplots indicate that FAEclust consistently finds the true number of clusters with low run-to-run variability across 100 repetitions.
F.4 Experiments on time-warped functional data
Figures II and 12 respectively illustrate the original and time-warped functional dataset in one repetition, with different clusters shown in different colors.
Tables and summarize the AMI and ARI performance, respectively, of the seven baseline methods on the timewwarped datasets. Overall, FAEclust delivers strong, robust clustering time- warped functional data, outperforming baselines in most scenarios. These results highlight the advantages of using a shape-informed metric for clustering problems where phase variation is nuisance.
For the 12 simulation scenarios, Figure shows boxplots of the number of clusters identified over 100 repetitions for each simulation scenario. The dashed line in each panel indicates the true numbers
28

---

lorenz - Feature 1 vs Time
lorenz - Feature 2 vs Time
30
20
6
0-
-10
-40
20
60
100
20
60
100
Time step
Time stop
Time
618
20
-10
60
100
Time step
Figure 7: Lorenz functional data. Trajectories generated from the Lorenz system, illustrating chaotic dynamics across different parameter settings.
swiss_roll Feature vs Time
swiss_rol colss_rroll- Feature vs Time
10
Value
-5
2-
0
50
100
150
200
50
100
150
200
stop
Wolls_roll- 3D Trajectories (2 features with time)
150
Time
10 5
2 4
Peature
Feature
10
Figure 8: Swiss roll functional data. Trajectories on a 2D manifold embedded in 3D space, with different colors indicating different clusters.
29

---

pendulum Feature vs Time
pendulum Feature vs Time
1.00
1.00
0.50
0.50
0.00
0.00
0.25
0.25
00.00
-0.50
-0.75
-0.75
100
100
Time step
Time step
pendulum 3D Trajectories features with time)
mi
1.0
0.5 posture?
-10_05
0.0
0.0
-0.5
Frature 1
0.5
Figure 9: Pendulum functional data. Trajectories in cylindrical phase space (IR 8t), with colors indicating different clusters.
Table 6: ARI scores for the 17 Euclidean functional datasets.
Dataset
funHDDC
funclust
FADP1
FADP2
FNN
FAE
VANO
FAEclust
BirdChicken
0.068
0.101
0.026
0.068
0.232
0.169
0.148
0.227
CBF
0.343
0.030
0.438
0.234
0.329
0.339
0.643
0.709
Chinatown
0.198
0.198
0.080
0.124
0.174
0.093
0,256
0.282
DSR
0.750
0.025
0.708
0.822
0.563
0.462
0.644
0.862
ECG200
0.186
0.158
0.129
0.144
0.134
0.137
0.119
0.164
Fungi
0.593
0.091
0.258
0.317
0.126
0.495
0.702
0.853
Plane
0.764
0.006
0.586
0.662
0.783
0.740
0.728
0.825
Rock
0.321
0.267
0.104
0.291
0.137
0.101
0.263
0.297
Symbols
0.664
0,000
0.362
0.576
0.630
0.638
0.676
0.682
Blink
0.369
0.039
0.164
0.153
0.418
0.449
0.442
0.563
BM
0.306
0.017
0.146
0.344
0.334
0.626
0.528
0.497
EOS
0.160
0.004
0.067
0.072
0.056
0.175
0.131
0.193
Epilepsy
0.095
0.014
0.049
0.166
0.180
0.137
0.193
0.328
ERing
0.641
0.031
0.172
0.651
0.576
0.677
0.669
0.591
FM
0.001
0.004
0.003
0.002
0.105
0.143
0,077
0.139
IV
0.810
0.054
0.182
0.344
0.085
0.843
0.887
0.879
SWJ
0.191
0.029
0.183
0.048
0.135
0.122
0.288
0.261
30

---

Cuoketng muthods FAt
Hypergpeer
Pyperbo
Swia
Penculuu
Figure 10: Per manifold scenario, 100 FD datasets generated and analyzed using the eight clustering methods; the boxplots (one per method) summarize the number clusters identified across the 100 runs. Only FAEclust and VANO yield consistent results, with FAEclust showing lower run-to-run variability.
Grwnarto
Connanio
Scenwato
16
-10
-40
0.0
22
0.0
0.2
3.4
6.5
0.
3.0
0.0
62
0.4
0.5
0.8
1,3
2.5
1.8
Sconario
Scenaroo
Sconarla
122
22
2.
6.
2.4
6.4
Figure/ 11: Example of original unwarped functional data from repetition. Different colors represent different clusters.
31

---

Table 7: ARI scores for the five coanifold-valued functional datasets. The table reports the mean (top
row) and standard deviation (bottom row) of the scores over 100 repetitions.
Dataset
funHDDC
funclust
FADP1
FADP2
FNN
FAE
VANO
FAEclust
Hypersphere
0.014
0.461
0.107
0.038
0.067
0.228
0.411
0.696
0.032
0.039
0.124
0.046
0.058
0.041
0.060
0.032
Hyperbolic
0.005
0.007
0.001
0.001
0.004
0.058
0.213
0.744
0.014
0.030
0.008
0.008
0.003
0.010
0.035
0.037
Swiss roll
0.069
0.097
0.271
0.109
0.016
0.057
0.142
0.205
0.034
0.032
0.054
0.071
0.026
0.018
0.045
0.049
Lorenz
0.103
0.348
0.063
0.126
0.036
0.180
0.189
0.416
0.062
0.048
0.042
0.079
0.039
0.025
0.018
0.031
Pendulum
0.847
0.316
0.735
0.744
0,195
0.680
0.874
0.985
0.031
0.053
0.076
0.068
0.033
0.050
0.047
0.009
Sconario
Sconario
Sconario
40
-10
-30
-40
0.0
32
2.5
0.0
9.5
0.8
142
Scenario
Scenadoo
Sconadi
2.5
130
1.38
0.5
Scamario
Scamanio
Scanario
-10
7.
1.5
60
72
0.4
05
1.5
0.0
0.4
0.8
10
24
Figure 12: Example of time-warped functional data from one repetition. Different colors represent different clusters.
32

---

Table 8: AMI scores of seven baseline methods on time-warped datasets. For each simulation scenario, the table shows the mean (top row) and standard deviation (bottom row) of AMI scores across 100 independent runs.
Scenario
funHDDC
funclust
FADP1
FADP2
FNN
FAE
VANO
A 0.506
0,407
0.411
0.402
0,470
0.436
0.654
0.101
0.011
0.012
0.026
0.064
0.019
0.029
B
0.096
0,013
0.009
0.097
0.001
0.000
0.005
0.029
0.024
0.025
0.069
0.006
0.005
0.017
C
0.542
0.199
0.303
0.402
0.500
0.457
0.571
0.078
0.103
0.094
0.013
0.006
0.025
0.029
D
0.299
0.069
0.290
0.321
0.058
0.057
0.489
0.051
0.093
0.073
0.071
0.029
0.025
0.080
E
0.469
0.687
0.578
0.878
0.868
0.890
0.871
0.094
0.084
0.010
0.024
0.029
0.015
0.030
F
0.664
0.277
0.338
0.324
0.308
0.459
0.664
0.049
0.078
0.061
0.066
0.051
0.044
0.045
G
0.285
0.005
0.202
0.149
0.280
0.253
0.206
0.036
0.011
0.045
0.061
0.043
0.044
0.062
0.592
0.579
0.366
0.579
0.641
0.724
0.861
0.047
0.013
0.117
0.032
0.032
0.026
0.058
0.592
0.579
0.336
0.560
0.635
0.658
0.672
0.052
0.022
0.109
0.019
0.030
0.015
0.057
0.001
0.001
0.000
0.000
0.000
0.000
0.001
0.003
0.008
0.005
0.005
0.004
0.005
0.006
0.366
0.003
0.200
0.049
0.129
0.097
0.057
0,085
0.004
0.097
0.061
0.104
0.101
0.065
L
0.381
0.001
0.514
0.494
0.728
0.714
0.648
0.047
0.006
0.076
0.060
0.064
0.061
0.129
of clusters. Again, FAEclust performs best: consistently recovers the true cluster number and exhibits less variability in the identified counts than competitors. Even in Scenario where variability increases, it remains the only method to find the true cluster count.
E5 Limitations
In all the experiments above, dynamic time warping computations performed on two 2.90 GHz CPUs (each with 16 cores), and the resulting distance matrices were used to train the FAE network on an NVIDIA Quadro RTX 5000 GPU with 16 GB of VRAM.
The main limitations of our model are sensitivity to hyperparameter settings and its higher computational cost. Even with fixed hyperparameter values, FAEclust takes significantly longer to run than most existing FD clustering methods. For example, on the "Plane" dataset from the UCR repository comprising 210 samples, each of length 144, and 7 clusters funHDDC completes in 24 seconds, FADPI 31 VANO in 51 while FAEclust requires 220 seconds. This makes it roughly ten times slower than the fastest baselines (funHDDC. FADPI, FADP2, and FAE). In practical applications, users can first apply fast baseline method such as FADP2 to obtain initial clustering results. These results can then serve as a warm start for our joint network training and clustering framework (Figure 2
seconds, FADP2 in seconds, funclust in 284 seconds, FNN 48 s. FAE in
To enhance clustering performance, we optionally include Bayesian optimization step for hyperpa- rameter tuning. While more sample-efficient than grid or random search, this step further increases the overall runtime. In practical applications, users must weigh FAEclust's improved modeling flexibility and accuracy - achieved through its deep learning foundation and shape-aware objective- against this added computational burden.
33

---

Table ARI scores of seven baseline methods on randomly generated time- warped datasets. For each simulation scenario, the table shows the mean (top row) and standard deviation (bottom row) of ARI scores across 100 independent runs.
Scenario
funHDDC
funclust
FADP1
FADP2
FNN
FAE
VANO
A
0.423
0.362
0.365
0.372
0.334
0.335
0.524
0.108
0,010
0.024
0.030
0.074
0.018
0.027
B
0.063
0.016
0.009
0.074
0.001
0.001
0.005
0.029
0.032
0.019
0.059
0.008
0.006
0.013
C
0.449
0.151
0.281
0.330
0.436
0.344
0.502
0.086
0.095
0.088
0.018
0.028
0.028
0.033
D
0.284
0.080
0.267
0.291
0.053
0.050
0,456
0.063
0.101
0.086
0.076
0.030
0.023
0.087
E
0.438
0.660
0.571
0.869
0.828
0.864
0,850
0.122
0.096
0.016
0.031
0.036
0.029
0.017
F
0,610
0.208
0.282
0.243
0.188
0.339
0,580
0,048
0.070
0.057
0.052
0.046
0.051
0.043
G
0,294
0.003
0.210
0.155
0.228
0.203
0.161
0.059
0.008
0.063
0.074
0.036
0.033
0.053
H
0.524
0.551
0.302
0.551
0.554
0.701
0.826
0.053
0.019
0.101
0.026
0.023
0.060
0.059
0.497
0.571
0.265
0.547
0.551
0.554
0.631
0.081
0.025
0.120
0.020
0.017
0.048
0.080
0.001
0.001
0.003
0.005
0.001
0.001
0.002
0.004
0.009
0.005
0.004
0.006
0.006
0.005
K
0.346
0.001
0.195
0.047
0.165
0.127
0.061
0.134
0.005
0.119
0.070
0.130
0.129
0.065
L
0.439
0.001
0.560
0.532
0.667
0.652
0.622
0.040
0.007
0.069
0.059
0.067
0.062
0.079
34

---

7MB
WANG
Anduv
30-
Scenati
Scanario
Sconario (4)
Somarii Ht Ci)
Scanari LC))
Scesaro
Figure 13: Per simulation scenario, 100 FD datasets were generated and analyzed using the eight clustering methods; the boxplots (one per method) summarize the number of clusters identified across the 100 runs. The numbers in the parentheses are the numbers of clusters.
35


Reviews:
[{"review_id": "Qyqu3G3kid", "reviewer": "Reviewer_CGH6", "summary": "This paper proposes a new unsupervised method to learn functional autoencoders for clustering multi-dimensional functional data. The new training objective applies to autoencoders, that involves a reconstruction loss, a clustering loss, and a regularization term on the functional weights and biases. The paper also discusses the forward and backward update of the clustering objective. The proposed method can handle functional data in both linear spaces and Riemannian manifolds. Experiments show that the proposed method outperforms existing state-of-the-art methods on learning time series data.", "questions": "1.\tIs there an ablation study on the influence of the architecture on the learning? Namely, is your method sensitive to the choice of architecture?", "limitations": "yes", "rating": 5, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths:**\n\n1.\tThe paper is well-written and easy-to-follow. The notation was clear, which makes for great readability. \n2.\tI am not an expert of the area, so I am not in the position to comment on its novelty. But suppose the limitations of the previous methods for unsupervised learning of functional data were true, then the method is novel and sufficiently addresses the shortcomings of previous methods. \n3.\tEmpirically, the proposed method seem to greatly outperforms existing methods. \n\n **Weaknesses:**\n\n1.\tPerhaps I missed the discussion in the paper, but it’s not clear how one can select the hyperparameter on the number of clusters. For example, L152 says “If an initial clustering of the FD is available”. In what cases will initial clustering of FD be available?", "was_revised": true, "final_justification": "The paper is well-written and easy to follow. I had some minor concerns to which the authors have addressed in their rebuttal. I am not an expert in the area as I have mentioned, so I cannot comment on the novelty and technique details. Suppose the novelty and technique is sound, then I would argue this is a fair contribution and deserves acceptance.", "rebuttal": "**Response to Weaknesses:**\n\nThank you for this observation. We agree that the discussion on determining the number of clusters should be made more prominent in the paper. Currently, this is only briefly mentioned in lines 314–315. We will revise the manuscript to clarify this point in multiple relevant sections. Below we provide a detailed explanation of our procedure:\n\nDuring the forward phase, we first update the embeddings and then apply our path-following homotopy algorithm to the new embeddings to obtain a hierarchy of agglomerative clustering solutions:\n1. According to Eq. (5), for a given embedding $\\pmb{X}$, the convex clustering objective has a unique solution for any value of $\\lambda$. \n2. As $\\lambda$ increases, clusters merge one pair at a time, tracing out a full agglomerative hierarchy. Crucially, these merges happen only at specific values of $\\lambda$ (termed breakpoints), such as $\\lambda_K$, at which two clusters merge, yielding a partition of exactly K clusters (see lines 278–283). Between any two consecutive breakpoints (for example, when $\\lambda_K\\leq \\lambda<\\lambda_{K-1}$), the clustering remains exactly the same.\n3. Our path-following homotopy algorithm efficiently identifies these breakpoints in O(n log(n)) time, thereby constructing the entire clustering hierarchy. This process can be further accelerated by initializing $\\lambda$ with a warm-start value (lines 305-308). \n\nOnce the hierarchy is constructed, we apply an internal validation index (e.g., silhouette score or Davies–Bouldin index) to select the optimal clustering result (and hence the number of clusters) from the hierarchy. Because this decision is driven entirely by the data, it does not depend on any particular $\\lambda$ value; hence, the $\\lambda$ in Eq. (5) is not a tunable hyperparameter.\n\nThe network training process consists of several iterations of the forward-backward loop. In each iteration, we repeat the clustering and selection procedure described above. The final clustering result from the last iteration is reported as the output of the model.\n\nRegarding the mention of \"an initial clustering of the FD,\" our intention was to convey that domain knowledge or simple and fast clustering algorithms (such as funHDDC) can be applied directly to the original functional data to obtain a rough initial estimate of cluster structure. While not required, such an initialization can help stabilize or accelerate training in practice. We will revise the text to better explain this and avoid potential confusion.\n\n**Response to Q1:**\n\nThank you for the insightful question.  FAEclust extends a classical autoencoder architecture into the functional domain by introducing functional layers at both the input and output ends. Specifically:\n\n- Encoder: The first layer is a functional layer where functional weights (continuous functions) map the multi-dimensional functional input into a finite-dimensional latent space. After this, one or more fully connected MLP layers (with scalar weights) further process the latent representation.\n- Decoder: The decoder consists of three or more functional layers, each involving matrix multiplications with functional weights and biases followed by nonlinear activation functions (except for the final layer which is linear and has no bias). These layers map the latent code back into the space of functions, enabling reconstruction.\n\nBecause FAEclust retains the same MLP core—depth, width, latent dimension, activation functions—as standard autoencoders, any architectural ablation study on these shared components would lead to similar conclusions to those already established in classical settings. (For example, Telgarsky 2016 [1] shows that additional depth and width exponentially expands representational power.) For this reason and to empirically validate the necessity of using multiple functional layers in the decoder, we focus our ablation study on the distinctive architectural component of FAEclust: the depth of the functional decoder. To this end, we compare two architectures:\n\n- a shallow variant with a single functional layer at the output end (architecture explained by the flow of information between lines 94 and 95 in the paper), and\n\n- the full version with three functional layers in the decoder (architecture explained by the flow of information between lines 109 and 110 in the paper).\n\nBoth models share identical hyperparameter settings, where applicable. We ran each model 20 times per dataset and report the mean and standard deviation of AMI and ARI. The variation is due to the stochastic nature of Bayesian optimization and the mini-batch gradient descent algorithm. The results across the 17 real datasets are summarized below. It is evident that the three-layer decoder consistently improves clustering performance in terms of both AMI and ARI, often substantially. This highlights that deeper functional decoders not only align with theoretical guarantees but also yield tangible performance gains in practice.\n\n| Datasets | One-layer FAEclust |  &emsp;  Three-layer FAEclust |\n|--|--|--|\n|  | &emsp;&emsp; **AMI** &emsp;&emsp;&emsp;&emsp;&emsp; **ARI** | &emsp;&emsp; **AMI** &emsp;&emsp;&emsp; &emsp; &emsp; **ARI** |\n| BirdChicken | $0.295 \\\\pm 0.053$ &emsp; $0.239 \\\\pm 0.061$ | &emsp; $0.329 \\\\pm 0.037$ &emsp; $0.260 \\\\pm 0.052$ |\n| CBF | $0.483 \\\\pm 0.071$ &emsp; $0.465 \\\\pm 0.064$ | &emsp; $0.719 \\\\pm 0.043$ &emsp; $0.689 \\\\pm 0.054$ |\n| Chinatown | $0.194 \\\\pm 0.056$ &emsp; $0.170 \\\\pm 0.045$ | &emsp; $0.319 \\\\pm 0.024$ &emsp; $0.252 \\\\pm 0.020$ |\n| DSR | $0.769 \\\\pm 0.021$ &emsp; $0.752 \\\\pm 0.027$ | &emsp; $0.881 \\\\pm 0.029$ &emsp; $0.865 \\\\pm 0.041$ |\n| ECG200 | $0.117 \\\\pm 0.050$ &emsp; $0.136 \\\\pm 0.104$ | &emsp; $0.240 \\\\pm 0.035$ &emsp; $0.160 \\\\pm 0.065$ |\n| Fungi | $0.840 \\\\pm 0.042$ &emsp; $0.712 \\\\pm 0.063$ | &emsp; $0.926 \\\\pm 0.019$ &emsp; $0.847 \\\\pm 0.048$ |\n| Plane | $0.841 \\\\pm 0.043$ &emsp; $0.792 \\\\pm 0.067$ | &emsp; $0.904 \\\\pm 0.016$ &emsp; $0.826 \\\\pm 0.039$ |\n| Rock | $0.405 \\\\pm 0.019$ &emsp; $0.286 \\\\pm 0.025$ | &emsp; $0.455 \\\\pm 0.023$ &emsp; $0.309 \\\\pm 0.026$ |\n| Symbols | $0.567 \\\\pm 0.049$ &emsp; $0.410 \\\\pm 0.061$ | &emsp; $0.839 \\\\pm 0.031$ &emsp; $0.674 \\\\pm 0.039$ |\n| Blink | $0.516 \\\\pm 0.057$ &emsp; $0.470 \\\\pm 0.074$ | &emsp; $0.637 \\\\pm 0.038$ &emsp; $0.572 \\\\pm 0.043$ |\n| BM | $0.389 \\\\pm 0.102$ &emsp; $0.326 \\\\pm 0.116$ | &emsp; $0.542 \\\\pm 0.037$ &emsp; $0.493 \\\\pm 0.046$ |\n| EOS | $0.182 \\\\pm 0.019$ &emsp; $0.137 \\\\pm 0.027$ | &emsp; $0.261 \\\\pm 0.014$ &emsp; $0.206 \\\\pm 0.019$ |\n| Epilepsy | $0.255 \\\\pm 0.051$ &emsp; $0.192 \\\\pm 0.025$ | &emsp; $0.464 \\\\pm 0.039$ &emsp; $0.288 \\\\pm 0.034$ |\n| ERing | $0.586 \\\\pm 0.073$ &emsp; $0.454 \\\\pm 0.088$ | &emsp; $0.657 \\\\pm 0.047$ &emsp; $0.582 \\\\pm 0.053$ |\n| FM | $0.259 \\\\pm 0.045$ &emsp; $0.164 \\\\pm 0.056$ | &emsp; $0.261 \\\\pm 0.034$ &emsp; $0.180 \\\\pm 0.051$ |\n| JV | $0.723 \\\\pm 0.038$ &emsp; $0.651 \\\\pm 0.063$ | &emsp; $0.883 \\\\pm 0.016$ &emsp; $0.870 \\\\pm 0.020$ |\n| SWJ | $0.159 \\\\pm 0.064$ &emsp; $0.096 \\\\pm 0.080$ | &emsp; $0.322 \\\\pm 0.063$ &emsp; $0.257 \\\\pm 0.051$ |\n|  |  |  |\n\n[1] Matus Telgarsky. Benefits of depth in neural networks. 29th Annual Conference on Learning Theory, PMLR 49:1517-1539, 2016."}, {"review_id": "qhDpgW0hyK", "reviewer": "Reviewer_9SuC", "summary": "The paper presents an autoencoder-based approach (FAEclust) to clustering multidimensional functional data. It is based on embeddings created by training a functional autoencoder minimizing a combined loss related to reconstruction error and clustering error. Moreover, an additional clustering objective (“shape-informed” objective) is optimized that takes into account the pairwise similarity of functional input data. Weight functions are represented as weighted linear sums of basis functions (qubic B-splines).", "questions": "Q1. It is not easy to judge the novelty of the proposed approach. The novel aspects with respect to its competitors are not sufficiently highlighted.\nQ2. The use of two clustering objective functions creates some confusion. The notation U is used in both L_c and L_s, however their dimensions are different. Why not integrating the regularization term of (5) into L_c, instead of performing separate optimizations?\nQ3. A pseudocode is needed to describe in detail the steps of the full method. Figure 2 is by no means sufficient. The algorithm in lines 515-530 refer to a loss function L, but it is not clear what L represents.\nQ4. How do the results change if the cluster update is omitted?\nQ5. There too many hyperparameters that need to be set in order to apply the method (including the number m of basis functions). It is not clear how all those parameters are set. If they are set taking into account the ground truth labels, then this procedure is not acceptable in a clustering setting. Typically, default values should be provided for all hyperparameters.\nQ6. The experiments with competitors are made using their default values, while for the proposed method a search for optimal hyperparameter values is conducted.\nQ7. In several datasets with AMI/ARI values are small (less than 0.35). This means that the ground truth clustering structure is not recovered, thus comparison is meaningless.\nQ8. There is no information on whether the ground truth number of clusters is recovered or not. It is just mentioned that an internal validation metric is used.", "limitations": "The main limitations of the method concern hyperparameter tuning and execution time.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths\nS1. The use of the shape informed objective introduces some novelty in the method.\nS2. Theoretical results are provided on the approximation capability of the functional decoder.\nS3. Better results have been obtained compared to other functional clustering approaches.\n\nWeaknesses\nW1. There are several concerns related to novelty, methodology, presentation and clarity (see questions section).\nW2. The method involves too many hyperparameters. It is mentioned that they are determined “through Bayesian optimization over a predefined search space”. It is not clear if such an optimization requires the ground truth cluster labels or not.", "was_revised": true, "final_justification": "Based on the clarifications provided, the additional experimental results and the positive opinion of other reviewers I have increased my score.", "rebuttal": "**Q1:  \"It is not easy to judge the novelty of the proposed approach.\"**\n\nWe apologize for the brevity in our original submission; space constraints limited how much detail we could include. The main contributions of the work are as follows.\n\n1. FAEclust is the first autoencoder architecture built specifically to multi-dimensional functional data *cluster analysis*. By using functional weights, it operates directly on functions rather than pointwise values, and its deep structure automatically captures complex dependencies among the component random functions.\n2. We establish that our decoder—which maps from a finite-dimensional space into an infinite‐dimensional function space—is a universal approximator in the Euclidean sense, and furthermore achieves patchwise universal approximation on any compact, connected Riemannian manifold. (This theoretical guarantee ensures the encoder never “avoid” regions of latent space for fear that the decoder cannot reconstruct them, a point we’ve further detailed in our response to Reviewer 5zSQ.)\n3. To address phase variation in functional datasets, we introduce a clustering objective that simultaneously captures shape information and remains separable across dimensions. We devise a fully data-driven method for selecting the neighborhood parameter m, and deliver an O(nlog(n)) optimization algorithm.\n4. FAEclust is available as a fully documented Python package (anonymous GitHub link in the paper), with sensible default hyperparameter settings.  Unlike related methods without released code, our package enables practitioners to apply, reproduce, and extend FAEclust immediately.\n\n\n**Q2:  \"The use of two clustering objective functions creates some confusion. The notation U is used in both L_c and L_s, ... Why not integrating the regularization term of (5) into Lc, instead of performing separate optimizations?\"**\n\nThank you for catching that ambiguity. We will rename $\\pmb{u}_k$ in $\\mathcal{L}_c$ to avoid confusion.\n\n$\\mathcal{L}_c$  is a regularization term in our training objective; it helps shape the latent space for clustering but is not itself a clustering criterion. Both  $\\mathcal{L}_c$ and  $\\mathcal{L}_s$ include a first term that measures within-cluster dispersion. Their key difference lies in the second term. In  $\\mathcal{L}_c$, the second term measures between-cluster dispersion, so adding  $\\mathcal{L}_c$ to the training objective encourages clusters to spread apart in the embedding space. In $\\mathcal{L}_s$, the second term captures shape information in the original functional data. Therefore, using $\\mathcal{L}_c$ in the loss prevents the latent representations from collapsing, a problem that can arise if one replaces  $\\mathcal{L}_c$ with $\\mathcal{L}_s$ in the network-training objective.\n\nWhy performing separate optimizations: network training and clustering serve distinct purposes and involve different objective functions. Once the network parameters (and hence the embeddings) are updated, the number of clusters and cluster centroids need be updated. In other words, both the number of clusters and the cluster centroids are unknown and dynamically determined. As such, they cannot be directly incorporated into the network training objective, which requires a fixed structure.\n\n**Q3:  \"A pseudocode is needed ... The algorithm in lines 515-530 refer to a loss function L, but it is not clear what L represents.\"**\n\nIn the backward network-training phase (see Figure 2),  we use stochastic gradient descent (Appendix D) to minimize the training objective function. In the forward cluster-update phase, we use the homotopy algorithm (Section 4.2) to minimize the clustering objective function. In the Python package (anonymous GitHub link in the paper), we explained the modular pipeline of the full method.\n\nAppendix D is for explaining and justifying the optimization algorithm for training the network, and hence the loss function $\\mathcal{L}$ there is the training objective function ($\\mathcal{L}=\\mathcal{L}_r+\\lambda_w\\mathcal{L}_w+\\lambda_c\\mathcal{L}_c$ in lines 141 & 142). \n\nA high-level logic of the  path-following homotopy algorithm is below:\n\n**Input**: latent embeddings $\\pmb{X}$, pairwise similarities $s(\\pmb{y}_i, \\pmb{y}_j)$, and $\\lambda=0$ (or a warm start).\n\n**While** number of clusters > 2:\n\n- Calculate the next breakpoint $\\lambda_{K-1}$ where a pair of centroids fuse (line 302).\n- Merge the corresponding clusters.\n\n**Return**: breakpoints {$\\lambda_{n-1}$, ..., $\\lambda_2$} and corresponding hierarchical clusterings.\n\nWe will integrate the above logic, Appendix D, and Figure 2 to provide a clearer explanation of the full method. Thanks for the suggestion.\n\n**Q4:  \"How do the results change if the cluster update is omitted?”**\n\nOur training process follows the usual pattern in “deep clustering” (as explained in lines 152-155): (1) pre-train an autoencoder on penalized reconstruction-loss only; (2) initialize clusters; (3) joint fine-tuning. If cluster update is omitted, then the final clusters are just the initial clusters from step (2). The pre-training only minimizes the reconstruction loss, and therefore the latent representations are agnostic to clustering structure (e.g., shape information). Iterative updating aligns latent representations with the clustering objective, leading to better-defined, more meaningful clusters.\n\n**Q5 and Weakness 2: \"There are too many hyperparameters. … If they are set taking into account the ground truth cluster labels, then this procedure is not acceptable in a clustering setting.”**\n\nWe will provide more information in Appendix E.1. As explained in lines 549-588, FAEclust shares hyperparameters with standard autoencoders. The additional hyperparameters are related to the functional nature of the network and the two penalty parameters ($\\lambda_c$ and $\\lambda_w$). Note that $\\lambda$ in Eq. (5) is not a hyperparameter; please see our first response to Reviewer CGH6 for further clarification. Tuning hyperparameters using Bayesian optimization is a standard and well-established approach for machine learning hyperparameter tuning (see, e.g., [1]). There’s no need for ground-truth cluster labels. Users simply define each parameter’s search domain. FAEclust is available as a fully documented Python package (anonymous GitHub link in the paper). Its README lists sensible default values for every hyperparameter.\n\n[1]: Turner, R. et al. “Bayesian Optimization is Superior to Random Search for Machine Learning Hyperparameter Tuning: Analysis of the Black-Box Optimization Challenge 2020”. Proceedings of the NeurIPS 2020 Competition and Demonstration Track, PMLR 133:3-26, 2021.\n\n**Q6:  \"The experiments with competitors are made using their default values, while for the proposed method a search for optimal hyperparameter values is conducted.\"**\n\nWe would like to clarify that for all baseline methods, we configured each algorithm to achieve its best possible performance, rather than simply using default settings.\n\n- funHDDC: We evaluated all six supported model types (line 538) and reported the one that gave the best performance for each dataset. Other parameters were related to algorithm convergence.\n- funclust: The parameter thd is the significance level in a hypothesis test, and setting thd=0.05 is a well-accepted choice. The parameters increaseDimension = FALSE and hard = FALSE were also tested and selected based on performance. Other arguments were kept at default as they pertain only to algorithmic stability or convergence.\n- FADPclust: We set proportion = NULL, which triggers the algorithm's internal grid search over 10 candidate values, as recommended by the authors.\n- Across all baseline methods, we conducted a grid search over the number of clusters and reported the result that achieved the best clustering performance. In contrast, FAEclust determines the number of clusters in a fully data-driven manner, as explained in our response to Reviewer CGH6.\n\n**Q7:  \"In several datasets with AMI/ARI values are small (less than 0.35). This means that the ground truth clustering structure is not recovered, thus comparison is meaningless.\"**\n\nThank you for the comment. We should have clarified that FAEclust was evaluated using a fixed architecture across all experiments: 7 layers in total (1 encoder + 3 decoder functional layers, and a 3-layer MLP), with maximum layer widths of {32, 16, 8, 16, 32, 32, 32}. Even with this conservative setup, FAEclust consistently outperformed all baseline methods. The key message we aim to convey is that FAEclust is substantially more effective than existing approaches.\n\nIf we were to increase the model capacity, FAEclust could achieve even better performance. For example, on the challenging hypersphere dataset, expanding the search space to {64, 32, 16, 32, 64, 64, 64} enables Bayesian optimization to find an architecture {64, 32, 16, 32, 64, 32, 32}, which achieves ARI = 0.5194 and AMI = 0.6489. Moreover, increasing the sample size from 100 to 500 yields ARI = 1 and AMI = 1, a perfect recovery of the true clustering. Please see our response to Reviewer WLrd for further details and additional supporting experiments.\n\n**Q8: \"There is no information on whether the ground truth number of clusters is recovered or not. It is just mentioned that an internal validation metric is used.''**\n\nWe use an internal validation metric to determine the number K of clusters in each forward phase, much like using classification accuracy in a supervised setting; please see our first response to Reviewer CGH6 for further clarification. To evaluate how well FAEclust recovers the true cluster structure, in our experiments, we reported AMI and ARI, both of which are computed against the ground-truth labels. We will include the number of identified clusters in our revised version.\n\nWe hope our responses have addressed the reviewer’s concerns. We would be pleased to provide any further clarification during the discussion period."}, {"review_id": "4pnM5yciqJ", "reviewer": "Reviewer_WLrd", "summary": "The authors propose FAEclust, a clustering algorithm for functional data, based on a functional autoencoder. It jointly learns a  nonlinear, shape-aware embeddings and performs convex clustering for multi-dimensional functional data. Across many Euclidean and five manifold datasets FAEclust consistently outperforms classic functional-data baselines, though performance drops on extremely low-signal manifold cases.", "questions": "-Please answer the question regarding the limitation below. **And**\nAs noted, the experimental evaluation would be more convincing if it included methods [9] and [23], which are closely related and represent strong baselines. Is it possible to include them in experiments?\n\nIf the authors adequately address the points raised above—particularly by discussing the observed limitations and including the suggested baselines—I would be inclined to raise my score. These additions would significantly strengthen the paper’s empirical foundation and overall clarity.", "limitations": "The paper demonstrates commendable breadth by evaluating on a diverse set of datasets—both Euclidean and manifold—and consistently outperforms selected baselines. However, key insights into why the method performs better, and under what conditions it fails, are either missing or only superficially addressed. This weakens the overall empirical evidence chain. For example, on many manifold datasets (Table 2), the clustering performance remains below 0.4, which raises concerns about the method's robustness in such settings. Could the authors provide an explanation for this behavior? It would also strengthen the paper to explicitly acknowledge this limitation in the discussion section.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "The paper proposes a new algorithm and strong results over some competitors (funHDDC, funclust, FADP). Unfortunately, the paper fails in giving the comprehensive details of these related works, does not even cite all of them. Instead, a few other related works are cited, e.g. [9] and [23] (see their paper) which appear to be more related as they apply also autoencoders to functional data. They exclude them from experimentation by the hint \"substantial computational challenges\". This is too general and does not convince me.\n\nWhile I believe that the algorithm is indeed new, the paper needs to discuss the novelty in the light of related work, not just pointing out \"my algorithm yields better numbers than something...\"\nOverall many components are known from the paper. The paper combines them in unique way therefore the novelty is medium. \n\n\nThe results in Table 1 are promising and suggest strong overall performance. However, the outcomes in Table 2 warrant more careful consideration. In particular, the consistently low performance on several manifold datasets calls for a deeper explanation. Clarifying the underlying factors contributing to these results—such as dataset characteristics, model limitations, or failure modes—would provide a more balanced and informative assessment of the method's strengths and weaknesses.\n\nAs noted, the experimental evaluation would be more convincing if it included methods [9] and [23], which are closely related and represent strong baselines. Including these would provide a more comprehensive comparison and help better contextualize the reported performance gains.", "was_revised": true, "final_justification": "In their rebuttal, the authors provided detailed results for the previously omitted methods, thereby enhancing the overall quality of the paper.", "rebuttal": "Thank you for these valuable suggestions. We have explained the novelty of our work in our response to Reviewer 9SuC. Before presenting additional comparative studies, we would like to clarify our choice of benchmarks: We compared FAEclust to every method surveyed in [34] that (a) is explicitly designed for clustering multi-dimensional functional data, and (b) provides publicly available code. While the approaches in [7] and [9] address similar problems, their implementations are not accessible. We did reach out to the authors of the two works when preparing the benchmark study; however, both informed us that they were unable to share their implementations. We excluded [23] because its autoencoder processes input and output as plain vectors of function evaluations, hence not fundamentally a functional network.\n\nTo ensure a fair and comprehensive comparison, we re-implemented the methods from [7] **and** [9] ourselves and included them in our comparative study, along with [23]. The network architecture of [7] follows Figure 1 in that work, while the architecture of [9] is based on Algorithm 1. For each method, we used Bayesian optimization for hyperparameter tuning. \n\nIn the table below, we compare the four methods across the 17 real datasets, with each model run 20 times per dataset. The reported variation is due to the stochastic nature of Bayesian optimization and the network training process. FAEclust consistently outperforms the other three methods on the majority of datasets. \n\n|Datasets|&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;FAEclust|&emsp;&emsp;&emsp;&emsp;&emsp;FNN[7]+K-means|&emsp;&emsp;&emsp;&emsp;&emsp;FAE[9]+K-means|&emsp;&emsp;&emsp;  VANO[23]+K-means|\n|--|--|--|--|--|\n| |&emsp;&emsp;&emsp;**AMI**&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;**ARI**|&emsp;&emsp;&emsp; **AMI**&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;**ARI**|&emsp;&emsp;&emsp;**AMI**&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;**ARI**|&emsp;&emsp;&emsp;**AMI**&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;**ARI**|\n|BirdChicken|&emsp;$0.329\\\\pm0.037$&emsp;$0.260\\\\pm0.052$|&emsp;$0.282\\\\pm0.026$&emsp;$0.219\\\\pm0.032$|&emsp;$0.266\\\\pm0.023$&emsp;$0.225\\\\pm0.033$|&emsp;$0.287\\\\pm0.030$&emsp;$0.164\\\\pm0.042$|\n|CBF|&emsp;$0.719\\\\pm0.043$&emsp;$0.689\\\\pm0.054$|&emsp;$0.367\\\\pm0.013$&emsp;$0.338\\\\pm0.020$|&emsp;$0.361\\\\pm0.015$&emsp;$0.343\\\\pm0.025$|&emsp;$0.723\\\\pm0.032$&emsp;$0.653\\\\pm0.040$|\n|Chinatown|&emsp;$0.319\\\\pm0.024$&emsp;$0.252\\\\pm0.020$|&emsp;$0.264\\\\pm0.035$&emsp;$0.198\\\\pm0.087$|&emsp;$0.132\\\\pm0.021$&emsp;$0.116\\\\pm0.027$|&emsp;$0.313\\\\pm0.038$&emsp;$0.239\\\\pm0.057$|\n|DSR|&emsp;$0.881\\\\pm0.029$&emsp;$0.865\\\\pm0.041$|&emsp;$0.688\\\\pm0.095$&emsp;$0.579\\\\pm0.112$|&emsp;$0.636\\\\pm0.030$&emsp;$0.462\\\\pm0.052$|&emsp;$0.704\\\\pm0.051$&emsp;$0.637\\\\pm0.063$|\n|ECG200|&emsp;$0.240\\\\pm0.035$&emsp;$0.160\\\\pm0.065$|&emsp;$0.230\\\\pm0.039$&emsp;$0.134\\\\pm0.032$|&emsp;$0.215\\\\pm0.035$&emsp;$0.122\\\\pm0.041$|&emsp;$0.219\\\\pm0.028$&emsp;$0.125\\\\pm0.023$|\n|Fungi|&emsp;$0.926\\\\pm0.019$&emsp;$0.847\\\\pm0.048$|&emsp;$0.261\\\\pm0.018$&emsp;$0.148\\\\pm0.013$|&emsp;$0.604\\\\pm0.025$&emsp;$0.475\\\\pm0.030$|&emsp;$0.821\\\\pm0.027$&emsp;$0.687\\\\pm0.042$|\n|Plane|&emsp;$0.904\\\\pm0.016$&emsp;$0.826\\\\pm0.039$|&emsp;$0.842\\\\pm0.031$&emsp;$0.787\\\\pm0.036$|&emsp;$0.832\\\\pm0.013$&emsp;$0.752\\\\pm0.016$|&emsp;$0.855\\\\pm0.029$&emsp;$0.702\\\\pm0.041$|\n|Rock|&emsp;$0.455\\\\pm0.023$&emsp;$0.309\\\\pm0.026$|&emsp;$0.198\\\\pm0.037$&emsp;$0.160\\\\pm0.033$|&emsp;$0.104\\\\pm0.028$&emsp;$0.116\\\\pm0.021$|&emsp;$0.354\\\\pm0.027$&emsp;$0.256\\\\pm0.032$|\n|Symbols|&emsp;$0.839\\\\pm0.031$&emsp;$0.674\\\\pm0.039$|&emsp;$0.748\\\\pm0.005$&emsp;$0.631\\\\pm0.026$|&emsp;$0.793\\\\pm0.008$&emsp;$0.645\\\\pm0.010$|&emsp;$0.794\\\\pm0.020$&emsp;$0.662\\\\pm0.041$|\n|Blink|&emsp;$0.637\\\\pm0.038$&emsp;$0.572\\\\pm0.043$|&emsp;$0.483\\\\pm0.039$&emsp;$0.430\\\\pm0.048$|&emsp;$0.510\\\\pm0.031$&emsp;$0.455\\\\pm0.014$|&emsp;$0.497\\\\pm0.045$&emsp;$0.419\\\\pm0.051$|\n|BM|&emsp;$0.542\\\\pm0.037$&emsp;$0.493\\\\pm0.046$|&emsp;$0.393\\\\pm0.044$&emsp;$0.326\\\\pm0.031$|&emsp;$0.671\\\\pm0.061$&emsp;$0.655\\\\pm0.050$|&emsp;$0.596\\\\pm0.042$&emsp;$0.530\\\\pm0.058$|\n|EOS|&emsp;$0.261\\\\pm0.014$&emsp;$0.206\\\\pm0.019$|&emsp;$0.161\\\\pm0.036$&emsp;$0.068\\\\pm0.049$|&emsp;$0.206\\\\pm0.024$&emsp;$0.163\\\\pm0.020$|&emsp;$0.179\\\\pm0.036$&emsp;$0.136\\\\pm0.056$|\n|Epilepsy|&emsp;$0.464\\\\pm0.039$&emsp;$0.288\\\\pm0.034$|&emsp;$0.290\\\\pm0.027$&emsp;$0.181\\\\pm0.013$|&emsp;$0.169\\\\pm0.033$&emsp;$0.122\\\\pm0.025$|&emsp;$0.176\\\\pm0.023$&emsp;$0.118\\\\pm0.018$|\n|ERing|&emsp;$0.657\\\\pm0.047$&emsp;$0.582\\\\pm0.053$|&emsp;$0.662\\\\pm0.052$&emsp;$0.593\\\\pm0.011$|&emsp;$0.762\\\\pm0.027$&emsp;$0.736\\\\pm0.036$|&emsp;$0.734\\\\pm0.040$&emsp;$0.684\\\\pm0.063$|\n|FM|&emsp;$0.261\\\\pm0.034$&emsp;$0.180\\\\pm0.051$|&emsp;$0.187\\\\pm0.035$&emsp;$0.094\\\\pm0.047$|&emsp;$0.159\\\\pm0.016$&emsp;$0.117\\\\pm0.023$|&emsp;$0.192\\\\pm0.029$&emsp;$0.086\\\\pm0.058$|\n|JV|&emsp;$0.883\\\\pm0.016$&emsp;$0.870\\\\pm0.020$|&emsp;$0.258\\\\pm0.072$&emsp;$0.098\\\\pm0.032$|&emsp;$0.857\\\\pm0.006$&emsp;$0.841\\\\pm0.008$|&emsp;$0.892\\\\pm0.006$&emsp;$0.883\\\\pm0.009$|\n|SWJ|&emsp;$0.322\\\\pm0.063$&emsp;$0.257\\\\pm0.051$|&emsp;$0.146\\\\pm0.039$&emsp;$0.127\\\\pm0.018$|&emsp;$0.154\\\\pm0.019$&emsp;$0.132\\\\pm0.030$|&emsp;$0.337\\\\pm0.068$&emsp;$0.286\\\\pm0.055$|\n| | | | | |\n\nIn the table below, we compare the four methods on the 12 simulated time-warped datasets. We generate 100 datasets for each scenario. FAEclust demonstrates robust and consistently superior performance across most scenarios. The drop of performance in Scenario L has been explained in Appendix E.3.\n\n||Original|Warped|Warped|Warped|Warped|\n|--|--|--|--|--|--|\n|Scenario|FAEclust|FAEclust|FNN|FAE|VANO|\n| |AMI ARI|AMI ARI|AMI ARI|AMI ARI|AMI ARI|\n|A|0.756 0.725|0.751 0.724|0.443 0.384|0.437 0.334|0.616 0.560|\n| |0.010 0.006|0.016 0.009|0.062 0.068|0.017 0.013|0.034 0.034|\n|B|0.426 0.446|0.405 0.383|0.003 0.004|0.002 0.003|0.004 0.005|\n| |0.069 0.084|0.092 0.114|0.008 0.011|0.006 0.009|0.010 0.013|\n|C|0.671 0.608|0.665 0.597|0.499 0.236|0.450 0.416|0.561 0.484|\n| |0.096 0.081|0.109 0.073|0.004 0.007|0.018 0.022|0.031 0.048|\n|D|0.573 0.529|0.559 0.526|0.041 0.038|0.042 0.035|0.294 0.257|\n| |0.094 0.075|0.088 0.062|0.019 0.021|0.018 0.013|0.038 0.039|\n|E|0.846 0.776|0.851 0.789|0.763 0.773|0.743 0.696|0.728 0.692|\n| |0.054 0.051|0.038 0.026|0.035 0.029|0.015 0.020|0.123 0.171|\n|F|0.864 0.793|0.839 0.781|0.309 0.185|0.477 0.358|0.658 0.548|\n| |0.046 0.039|0.032 0.045|0.035 0.034|0.035 0.046|0.047 0.050|\n|G|0.579 0.478|0.585 0.487|0.285 0.222|0.284 0.226|0.216 0.170|\n| |0.127 0.142|0.103 0.158|0.034 0.031|0.047 0.038|0.055 0.041|\n|H|0.998 0.993|0.997 0.991|0.629 0.542|0.723 0.709|0.858 0.848|\n| |0.008 0.005|0.011 0.006|0.040 0.022|0.026 0.032|0.055 0.068|\n|I|0.971 0.960|0.968 0.955|0.632 0.550|0.658 0.552|0.687 0.652|\n| |0.018 0.017|0.009 0.008|0.034 0.022|0.019 0.090|0.049 0.062|\n|J|0.151 0.080|0.142 0.116|0.001 0.001|0.000 0.001|0.002 0.002|\n| |0.076 0.053|0.092 0.089|0.004 0.005|0.005 0.007|0.003 0.003|\n|K|0.873 0.825|0.876 0.824|0.082 0.106|0.083 0.110|0.271 0.260|\n| |0.022 0.015|0.032 0.024|0.082 0.103|0.084 0.110|0.049 0.050|\n|L|0.704 0.729|0.359 0.370|0.656 0.678|0.636 0.702|0.632 0.640|\n| |0.087 0.064|0.108 0.123|0.048 0.042|0.035 0.038|0.066 0.086|\n| | | | | | |\n\nThe lower performance in Table 2 is due to the narrow architecture of FAEclust. We should have provided more detail regarding the architecture choices of FAEclust. In all experiments—both real and simulated—we fixed the architecture to a seven-layer structure: one functional layer in the encoder, three functional layers in the decoder, and a three-layer MLP in between. The third layer in FAEclust serves as the bottleneck, producing the latent embedding. Please refer to our response to Reviewer CGH6 for an explanation of this choice. The maximum number of nodes allowed per layer was constrained to {32, 16, 8, 16, 32, 32, 32}. For the five datasets in Table 2, the final architectures were {16, 8, 2, 8, 16, 16, 16} for Hypersphere and Lorenz, and {32, 16, 8, 16, 32, 32, 32} for Hyperbolic, Swiss Roll, and Pendulum. Even under this restricted architecture, FAEclust consistently outperforms all benchmark methods. \n\nWhile a small network may guard against overfitting with only 100-300 samples, it underfits when tasked with capturing intricate, non‑linear patterns. If we were to increase the model capacity, FAEclust could achieve even better performance. The below table shows the new clustering results, where we expand the width-search space to {64, 32, 16, 32, 64, 64, 64}. These results confirm that FAEclust can achieve better performance when given additional capacity, especially on more complex datasets.\n\n|Datasets|FAEclust|||FNN||FAE||VANO||\n|:--|:--|:--|:--|:--|:--|:--|:--|:--|:--|\n| |layer width|AMI|ARI|AMI|ARI|AMI|ARI|AMI|ARI|\n|Hypersphere|{64, 32, 16, 32, 64, 32, 32}|0.649|0.519|0.016|0.023|0.050|0.037|0.343|0.327|\n|Hyperbolic|{64, 32, 16, 32, 32, 32, 32}|0.842|0.814|0.001|0.003|0.043|0.060|0.028|0.039|\n|Swiss Roll|{64, 32, 16, 32, 64, 64, 64}|0.459|0.376|0.014|0.018|0.120|0.106|0.104|0.068|\n|Lorenz|{32, 16, 16, 16, 32, 32, 32}|0.488|0.440|0.023|0.028|0.154|0.106|0.229|0.180|\n|Pendulum|{32, 16, 16, 16, 16, 16, 32}|1|1|0.254|0.195|0.781|0.621|0.958|0.941|\n| | | | | | | | | | |"}, {"review_id": "W9qD6EOXnK", "reviewer": "Reviewer_5zSQ", "summary": "This work develops a clustering method for multivariate functional data based on deep functional autoencoder. It incorporates layers with functional weights in the network and clustering loss in the objective function to adapt classic autoencoder for the functional clustering task. Additionally, a shape-informed SRV metric is used in the clustering objective function to account for phase variation in functional data. A universal approximation theorem is presented for reconstruction accuracy and numerical experiments confirm the effectiveness of the proposed method.", "questions": "1. Could you explain briefly the methods used as the benchmarks? For example, which category of the two types of clustering methods (line 21-24) do they fall into? Are they deterministic or probabilistic?\n2. Could you explain how do you explore the number of clusters (line 535-536)? Does it recover the true number of clusters (if available)? Do any of the benchmarks include choosing the number of clusters in their framework?\n3. Could you further clarify your theoretical result based on the questions in Limitations, first point (see below)?\n4. Could you provide additional comparisons based on the questions in Limitations, second point (see below)?", "limitations": "1. It may be better to provide readers some intuitions and clarifications before and after Theorem 1. For example, what question it tries to answer?  Why it is useful? How does it connect and compare to classic universal approximation theorem? Also, the readout map $\\rho$ is not clearly defined in the theorem.\n2. $\\textbf{Shape-informed clustering}$ is the first phrase in the title, but only the clustering objective function in Section 4.1 is relevant to it. Is this one of the designs in your method that you would like to highlight? If it is the case, maybe it would be better to put more efforts in motivating and demonstrating the merits of introducing a shape-informed metric, both throughout writing and in experiments. For example, if a classic $\\mathbb{L}^2$ metric is used, how would the results change in Table 1? If phase variation is introduced (Appendix E.3), how would the benchmark methods perform?\n3. In general it would be better to report uncertainty for evaluation metrics (Table 1 and 2), e.g., repetition over data splits.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strength:\n\n1. The main paper has very clear structure and notation and is very well-written. The appendix is very comprehensive and include sufficient amount of details that are not included in the main paper.\n2. The idea of using autoencoder for functional clustering is very interesting and novel. The method is flexible and well-designed for this task.\n\nWeaknesses:\n\n1. The connection between theoretical result (existence of a good functional network) and practical design of the method (how to find it) is relatively weak.\n2. There are a lot of modeling choices, e.g., penalty tuning, number of basis for functional weights, number of clusters, which can potentially affect the robustness of the method.\n3. As pointed out by the authors, the method is computationally expensive.", "was_revised": true, "final_justification": "The rebuttal answers my questions. The paper is well-written with sufficient amount of novelty, supported by comprehensive theoretical and empirical results.", "rebuttal": "Thank you for these constructive comments, which, together with those from other reviewers, have helped us strengthen the paper’s contribution and clarify key aspects.\n\n**Response to Weakness:** In our response to Reviewer CGH6 (Q1), we presented additional experiments demonstrating the importance of using three functional layers in the decoder (rather than a single layer). These results empirically validate the necessity of multiple functional layers and support Theorem 1. Across all experiments, we fixed the network architecture to 7 layers, which we found sufficient in practice. We then apply Bayesian optimization to select the width of each layer, with the number of nodes per layer constrained to the set {32, 16, 8, 16, 32, 32, 32} as the maximum allowed configuration. In the revised manuscript, we will provide layer-wise information of the trained FAEclust for each dataset.  For more challenging datasets (e.g., manifold-valued data), users may increase the layer widths to further enhance FAEclust’s representational power. Please see our response to Reviewer WLrd for further details and supporting experiments.\n\nRegarding modeling choices, the two penalty parameters $\\lambda_w$ and $\\lambda_c$,  along with the number of basis functions for the four functional layers, are additional hyperparameters introduced beyond those in traditional autoencoders.  The parameter $\\lambda$ in the clustering objective function and the number of clusters are not hyperparameters (please see our first response to Reviewer CGH6 for clarification). To assess robustness, we ran FAEclust 20 times per dataset in all experiments reported above, presenting the mean and standard deviation of ARI and AMI scores. The results (summarized in the tables above) demonstrate that FAEclust is stable and robust across datasets.\n\n**Response to Q1:** Thank you for the comment. In our revision, we will also cite the publications corresponding to each R package in the introduction.\n\n- funHDDC [1]: This method projects functional data onto a large set of basis functions and models the projection coefficient vectors within each cluster as lying in a low-dimensional subspace. A parametric model is then applied to these coefficient vectors for clustering. funHDDC falls under the third category (lines 24-26) and is probabilistic in nature.\n- funclust [2]: This method also employs a parametric model but uses functional PCA scores (instead of general basis coefficients) as inputs for clustering. It falls under the first category and is also probabilistic.\n- FADPclust [3]: This method applies the density peak clustering algorithm directly to the pairwise distance matrix of functional data, making it a deterministic method in the second category.\n\n[1] Bouveyron and Jacques. 2011. Model-based clustering of time series in group-specific functional subspaces. Advances in Data Analysis and Classification.\n\n[2] Jacques and Preda. 2013. Funclust: A curves clustering method using functional random variables density approximation. Neurocomputing.\n\n[3] Wang and Xu 2015. Fast clustering using adaptive density peak detection. Statistical Methods in Medical Research.\n\n**Response to Q2:** All baseline methods (including FNN, FAE and VANO) require the number of clusters to be specified as an input. In contrast, FAEclust automatically determines the number of clusters in a fully data-driven manner using an internal validation index.\n\nFor all baseline methods, we run each method with the number of clusters ranging from 2 to 10, ensuring that the true number of clusters is included in this range. We then report the best clustering performance based on ARI and AMI.\n\nFor FAEclust, the number of clusters is adaptively determined within each forward phase of the joint training and clustering framework. Please refer to our first response to Reviewer CGH6 for a detailed explanation of this procedure.\n\n**Response  to Q3 and Limitation 1:** Thanks for the suggestion. We will give more explanation after Theorem 1 in the revision.\n1.  Why Theorem 1 matters: **(a)** When a decoder cannot faithfully reconstruct every possible latent code, the encoder quickly discovers which regions of the latent space will lead to plausible outputs and which will not. In effect, the encoder “steers clear” of those latent directions where the decoder’s mapping breaks down. As a result, the geometry of the learned representation is warped: some dimensions or combinations of dimensions are never explored, because they would produce meaningless or distorted reconstructions. Consequently, any clustering performed on this latent representation yields groupings that conform to the decoder’s permissible manifold rather than to the data’s intrinsic structure.  **(b)** Additionally, limited decoder capacity leads to consistent failures on certain patterns or classes. The network compensates by balancing reconstruction loss against clustering objectives, resulting in unpredictable trade-offs that undermine both reconstruction quality and cluster integrity. **(c)** Finally, even small differences between the initial clusters (after pre-training) and the final clusters (after fine-tuning) become magnified by a weak decoder, making the learned clusters fragile to noise. Theorem 1 guarantees that our decoder has the universal-approximation property, so none of these pathologies can arise.\n2. Connection to classic universal approximation theorem: This was briefly explained in lines 399-405. While the universal-approximation property of finite-dimensional neural networks is well understood, far less is known about infinite-dimensional mappings ([6], [11]). Unlike existing functional networks that rely on integral kernels, our framework uses univariate functions as weights. This key distinction requires a bespoke universality proof. In our paper, we not only establish Euclidean universal approximation but also extend it to patchwise approximation on compact, connected Riemannian manifolds, substantially strengthening and generalizing the classical result.\n3. The full construction of the readout map $\\rho$ appears later in lines 479-485. In the revised manuscript, we will add a concise explanation of $\\rho$ directly in the main text.\n\n**Response to Q4 and Limitation 2:** Phase variation (a.k.a., time warping) in functional data refers to variability in the timing of key features across curves. Addressing phase variation in functional data is crucial because it distinguishes true shape differences from misalignments in time or space. In most functional data studies, curve alignment is treated as a separate pre-processing step, independent of subsequent analysis. As noted in [34] (Section 7), only a small number of works explicitly integrate phase variation handling within clustering, highlighting the challenge of this problem.\n\nUnlike classical symmetry groups such as the Euclidean group E(2), the group of warping functions (acting on functions via composition) is nonlinear and does not admit a standard (left) group representation. This nonlinear structure means we cannot easily design network layers that are equivariant or invariant to phase variation in the same way as for transformations such as rotation. In our current work, we handle phase variation through the clustering objective that incorporates a shape-informed elastic metric, which captures phase-amplitude separation. Developing neural network architectures explicitly designed to identify and disentangle phase variation remains an important future direction.\n\nIn Table 1, most datasets there do not exhibit phase variation. To investigate the effect of ignoring phase variation, we applied the L2-variant of FAEclust (replacing the elastic distance in Eq. (6) with the standard L2 distance) across the 12 simulation scenarios. The clustering results are shown below. For comparison with other baseline methods, please refer to the second table in our response to Reviewer WLrd. These results highlight the advantages of using a shape-informed metric for clustering problems where phase variation is a nuisance.\n\n| |Warped | | Warped| |\n|:--|:--|:--|:--|:--|\n|Scenario|Elastic||$L^2$||\n||AMI|ARI|AMI|ARI|\n|A|0.751|0.724|0.582|0.465|\n||0.016|0.009|0.021|0.038|\n|B|0.405|0.383|0.009|0.002|\n||0.092|0.114|0.004|0.001|\n|C|0.665|0.597|0.574|0.460|\n||0.109|0.073|0.042|0.059|\n|D|0.559|0.526|0.274|0.216|\n||0.088|0.062|0.011|0.019|\n|E|0.851|0.789|0.717|0.570|\n||0.038|0.026|0.034|0.028|\n|F|0.839|0.781|0.735|0.613|\n||0.032|0.045|0.024|0.036|\n|G|0.585|0.487|0.086|0.051|\n||0.103|0.158|0.009|0.012|\n|H|0.997|0.991|0.734|0.672|\n||0.011|0.006|0.037|0.029|\n|I|0.968|0.955|0.733|0.591|\n||0.009|0.008|0.022|0.045|\n|J|0.142|0.116|0.003|0.002|\n||0.092|0.089|0.001|0.002|\n|K|0.876|0.824|0.163|0.094|\n||0.032|0.024|0.038|0.037|\n|L|0.359|0.370|0.419|0.361|\n||0.108|0.123|0.064|0.080|\n||||||\n\n**Response to Limitation 3:** Thank you for the suggestion. In Table 1, most datasets are small (around 50 curves) to medium (around 200 curves) in size. For example, BirdChicken has 40 curves, Rock has 70, and BM has 80. Bootstrapping is less suitable for unsupervised learning because clustering depends on global data structure and lacks ground truth labels to anchor resamples. Resampling can distort the underlying relationships, especially in small datasets, and lead to unstable or non-comparable cluster solutions. To assess the robustness of FAEclust, we instead ran FAEclust 20 times per dataset and report the mean and standard deviation of AMI and ARI in the tables above (in responses to Reviewer CGH6 and to Reviewer WLrd). The observed variation arises from the stochastic nature of Bayesian optimization and the mini-batch gradient descent algorithm.\n\nFor the 12 simulation scenarios, we generated 100 independent datasets for each scenario and report the mean (top row) and standard deviation (bottom row) of ARI and AMI, thereby quantifying uncertainty due to both data variability and algorithmic randomness."}]

Author rebuttal:

Meta-review: ### Summary of scientific claims and findings

The paper introduces a functional autoencoder-based method for clustering multi-dimensional functional data in both Euclidean and Riemannian manifolds. The approach incorporates functional layers with basis function representations, a shape-informed clustering objective to address phase variation, and a convex clustering procedure. The authors prove a universal approximation theorem for the functional decoder and propose a homotopy-based algorithm for efficient clustering updates. 

### Strengths

* **Clear contribution**: Extends deep autoencoders to functional data clustering with theoretical guarantees.
* **Novelty**: Introduces functional layers and a shape-informed clustering loss, which address limitations of existing functional clustering methods.
* **Breadth of evaluation**: Covers real, simulated, and manifold datasets, with systematic reporting of AMI and ARI across multiple repetitions.
* **Rebuttal improvements**: Added ablation studies, clarified the cluster selection procedure, re-implemented missing baselines, and reported robustness results.
* **Reproducibility**: Code package available with sensible defaults and documentation.

### Weaknesses and missing elements

* **Presentation issues**: Initial submission lacked pseudocode and clear explanations of some objectives, leading to confusion.
* **Hyperparameter complexity**: Method involves many hyperparameters, raising robustness and practicality concerns.
* **Novelty positioning**: While components are individually known (autoencoders, convex clustering, shape metrics), the distinct contribution relative to prior autoencoder-based works was under-emphasized initially.
* **Comparative evaluation**: Important baselines were missing in the original submission and only included after reviewer prompting.

### Key reasons for acceptance

1. The paper establishes a sound, well-motivated framework for clustering functional data, with both theoretical and empirical contributions.
2. The rebuttal successfully clarified earlier ambiguities and substantially strengthened the empirical evaluation by including missing baselines.
3. The work is novel in integrating functional layers, shape-informed clustering, and convex optimization into a unified deep clustering pipeline.

### Summary of discussion and rebuttal

**Before rebuttal**:

* `CGH6`: Positive on clarity and empirical performance but lacked domain expertise; questioned cluster selection procedure and architectural sensitivity.
* `9SuC`: Raised multiple concerns about novelty, clarity, hyperparameter tuning, and missing cluster-recovery details; leaned negative.
* `WLrd`: Questioned novelty positioning, missing baselines, and low manifold performance.
* `5zSQ`: Highlighting clarity, novelty, and design, but noted weak theory-to-practice connection and robustness concerns.

**Authors’ rebuttal**:

* Clarified cluster selection and hyperparameter tuning (no ground truth labels needed).
* Provided ablation on functional decoder depth.
* Re-implemented and added comparisons with missing baselines.
* Reported robustness via repeated runs and variance statistics.
* Addressed phase variation via elastic shape metric, showing clear benefits.

**Post-rebuttal positions**:

* `CGH6`: Maintained positive stance but deferred on novelty due to lack of expertise.
* `9SuC`: Raised to borderline accept but retained concerns about hyperparameters and reporting cluster numbers.
* `WLrd`: Increased score after additional experiments and baseline inclusion.
* `5zSQ`: Raised rating, satisfied with clarifications.

### Weighing these points

While initial reviews expressed substantial concerns, the rebuttal was strong and persuasive, addressing methodological clarity, missing baselines, and robustness. Remaining weaknesses relate mainly to presentation and complexity rather than to core soundness. 

### Additional AC review

The integration into the current functional data literature appears a bit weak. Several papers on functional data were presented at last year's NeurIPS, but are not cited, nor are other functional machine learning approaches beyond clustering. I would encourage the authors to add these references.
Decision: Accept (poster)