{'fraction_non_alphanumeric': 0.042441321657556784, 'fraction_numerical': 0.018326934352126796, 'mean_word_length': 4.222123456790124, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 35, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interpretable discriminant analysis for functional data supported on random non-linear domains

We introduce a novel framework for the classification of functional data supported on non-linear, and possibly random, manifold domains. The motivating application is the identification of subjects with Alzheimer's disease from their cortical surface geometry and associated cortical thickness map. The proposed model is based upon a reformulation of the classification problem into a regularized multivariate functional linear regression model. This allows us to adopt a direct approach to the estimation of the most discriminant direction while controlling for its complexity with appropriate differential regularization. Our approach does not require prior estimation of the covariance structure of the functional predictors, which is computationally not feasible in our application setting. We provide a theoretical analysis of the out-of-sample prediction error of the proposed model and explore the finite sample performance in a simulation setting. We apply the proposed method to a pooled dataset from the Alzheimer's Disease Neuroimaging Initiative and the Parkinson's Progression Markers Initiative, and are able to estimate discriminant directions that capture both cortical geometric and thickness predictive features of Alzheimer's Disease, which are consistent with the existing neuroscience literature.


Introduction
Functional discriminant analysis, which is the problem of predicting categorical outcomes from functional predictors, has been well studied within the field of Functional Data Analysis (FDA) (Ramsay and Silverman 2015;Ferraty and Vieu 2006;West and Harrison 1997;Hsing and Eubank 2013) and has motivated a large body of literature (see e.g., James and Hastie 2001;Müller 2005;Preda, Saporta, and Lévéder 2007;Delaigle and Hall 2012;Dai, Müller, and Yao 2017;Berrendero, Cuevas, and Torrecilla 2018;Kraus and Stefanucci 2019;Park, Ahn, and Jeon 2021). However, most of the proposed methods are concerned with the classification of functional data supported on one-dimensional linear domains, which is a limiting assumption in many relevant applications. Recent work on the analysis of functional data with manifold structure has mostly focused on vector-valued functional data with non-Euclidean constraints in the image space (see e.g., Su et al. 2014;Dai and Müller 2018;Lin and Yao 2019;Dubey and Müller 2020;Kim et al. 2021).
In this paper, motivated by the analysis of modern multi-modal imaging data, we propose a classification model for functional predictors supported on non-linear sample-specific manifold domains, which we call Functions on Surfaces (FoSs) (Lila and Aston 2020). An example of such data is provided in Figure 1A. The statistical analysis of FoSs comes with challenges that are largely unaddressed. This is mainly due to the non-Euclidean structure of each individual domain, which poses challenges to the definition of appropriate spatial regularizations, and due to the more abstract non-Euclidean structure of the latent space where the random domains are supported, which further invalidates traditional linear statistical models (Kendall 1984;Grenander and Miller 1998;Younes 2019;Dryden and Mardia 2016).
The approach proposed in this paper tackles the non-Euclidean nature of the latent space of random domains by defining appropriate linear functional representations of FoSs, effectively recasting the problem of classifying FoSs into the problem of classifying bivariate functional data supported on general, but fixed, domains. We then define a direct model for the estimation of a functional discriminant direction in the parametrizing space exploiting a convenient reformulation of the problem into a regularized functional linear regression model, where appropriate differential spatial regularization is introduced to produce interpretable and well-defined estimates. While several regularized regression models for functional data supported on multidimensional domains have been proposed (Goldsmith, Huang, and Crainiceanu 2014;Wang and Zhu 2017;Kang, Reich, and Staicu 2018;Feng et al. 2020;Yu et al. 2021), these are limited to flat domains as they rely on a global Cartesian coordinate system, which is in general not available on non-linear domains. Moreover, some of these approaches rely on pre-computations of the covariance structure of the functional predictors, which is ultimately not possible in our final application due to the high dimensionality of the data.
Our motivating application is the identification of subjects with Alzheimer's disease from FoSs that are subject-specific cortical surfaces coupled with cortical thickness measurements. In fact, it is well known that Alzheimer's disease and related neurodegenerative diseases affect both the geometry of the cortical surface of an individual and their cortical thickness (Thompson et al. 2001). However, current statistical approaches to the analysis of such data have mostly focused on the case of spatially normalized functional data, i.e. where the geometric domain is fixed (see e.g. Chung, Hanson, and Pollak 2021;Lila, Aston, and Sangalli 2016;Mejia et al. 2020). On the other hand, most of the methods for the analysis of subject-specific geometries are not able to handle coupled functional components, e.g. neuronal activation or thickness measurements (Epifanio and Ventura-Campos 2014;Ferrando, Ventura-Campos, and Epifanio 2020).
Little attention has been given to the most general setting of functional data supported on random domains, with the exception of Charlier, Charon, and Trouvé (2017) and Lila and Aston (2020), where the authors focus on unsupervised models, and Lee et al. (2017), where the authors leverage the fshape framework (Charlier, Charon, and Trouvé 2017) to represent such data and propose a linear discriminant model on finite-dimensional representations of the predictors. The framework proposed here bypasses any dimension reduction steps applied to the predictors. We instead leverage the spatial nature of the data allowing for direct control of the complexity of the discriminant directions through efficient spatial regularization.
The rest of the paper is organized as follows. In Section 2, we describe the representation model adopted to parametrize the non-Euclidean space of FoSs in terms of linear function spaces. In Section 3 and 4, we develop a novel discriminant analysis model on the parametrizing linear function spaces and provide theoretical guarantees for the prediction performance of the proposed approach. We propose an efficient discretization approach in Section 5 and apply the proposed model to the Alzheimer's Disease Neuroimaging Initiative (ADNI) and Parkinson's Progression Markers Initiative (PPMI) data in Section 6. Proofs and simulations are left to the appendices.
Here v i : R 3 → R 3 is a vector-valued function encoding the geometry of the ith individual. This is depicted as a collection of 3D vectors {v i (p j )} for a dense set of points {p j } ⊂ R 3 . For clarity, the function v i is displayed only on half of its domain R 3 . The function x i : M → R describes the spatially normalized thickness map of the ith individual on the fixed template M.
These is a close approximation of the FoS (M i , z i ) in Panel A.

Functional data supported on general random domains
The data motivating this work are a sample of triplets where g i is a binary label, M i ⊂ R 3 is a sample-specific closed two-dimensional manifold embedded in R 3 , and z i : M i → R is a scalar function supported on the geometric object M i . We moreover assume that the points on the geometries {M i } of the observed FoSs are in one-to-one correspondence across subjects. We refer to the pairs {(M i , z i )} as FoSs.
In Figure 1A, we display three observations of the training sample of our final application, where g i encodes the disease state of the subject, M i encodes the geometry of the cerebral cortex, and z i encodes the cortical thickness map supported on M i . Our goal is to build a classifier from the given training sample that can predict the binary label g * of a previously unseen FoS (M * , z * ).

Linear functional representations
In our motivating application, the geometric objects {M i : i = 1, . . . , n} are topologically equivalent to a sphere, and therefore, they do not display holes or self-intersections. To inform our model of such physical non-Euclidean constraints, we define a convenient unconstrained representation model for the FoSs {(M i , z i ) : i = 1, . . . , n} in terms of objects belonging to linear function spaces.
Let M be a template two-dimensional manifold embedded in R 3 that is topologically equivalent to a sphere. We denote by L 2 (M) the space of square integrable functions over M, equipped with the standard inner product ·, · L 2 (M) and norm · L 2 (M) and denote by L 2 (R 3 , R 3 ) the space of square integrable vector-valued functions from R 3 to R 3 , with inner product ·, · L 2 (R 3 ,R 3 ) and norm ) be a Reproducing Kernel Hilbert Space (RKHS) of smooth functions with compact support in R 3 . We then introduce a diffeomorphic operator ϕ such that ϕ v : R 3 → R 3 is a diffeomorphic function for every choice of v ∈ V(R 3 ) (Vaillant et al. 2004). We denote by ϕ v (M) the geometric object given by displacing every point p ∈ M ⊂ R 3 to the new location ϕ v (p). A direct consequence of ϕ v being diffeomorphic is that ϕ v (M) is topologically equivalent to a sphere for every choice of v ∈ V(R 3 ). The construction of the diffeomorphic operator adopted in this paper and the computation of where v i ∈ V(R 3 ) and x i ∈ L 2 (M). Intuitively, this means that every geometric object M i is modeled as a (diffeomorphic) deformation of the template, i.e. ϕ v i (M), while the associated function z i is given by 'transporting' the signal x i onto M i with such a deformation. The approach described allows us to recast the original non-Euclidean problem of learning a classifier from the training sample {(g i , M i , z i )|i = 1, . . . , n} into the problem of learning a classifier from where v i ∈ V(R 3 ) ⊂ L 2 (R 3 , R 3 ) and x i ∈ L 2 (M) are two functional predictors both belonging to linear function spaces. In Figure 1B, we show the functional linear representations associated with the three FoSs in Figure 1A. Most importantly, the representation mapping adopted here is 'invertible', i.e. any pair of estimates β G ∈ V(R 3 ), β F ∈ L 2 (M), such as a 'direction' that best discriminates between two classes, can be mapped onto the original space, defining the associated trajectories of FoSs where ϕ c 1 β G (M 0 ) is guaranteed to be topologically equivalent to a sphere, and therefore, to satisfy the physical constraints of the problem considered.

Discriminant analysis on the parametrizing linear function spaces
The aim of Section 3 and 4 is to provide methodology for learning a linear classifier, from the training data {(g i , v i , x i )|i = 1, . . . , n} displayed in Figure 1B, by introducing a novel functional classification model that has the following crucial characteristics: • Does not rely on Functional Principal Components Analysis (FPCA), or related dimension reduction methods, to reduce the dimension of the functional predictors, bypassing the intrinsic assumption that the discriminant direction is well represented by the space spanned by the first few (unsupervised) PC functions; • Can be applied to bivariate, and possibly multivariate, functional predictors each supported on different domains; • Allows for explicit spatial regularization of the estimates on possibly non-linear manifold domains, yielding well-defined and interpretable estimates; • Provides a direct estimation of the discriminant directions without relying on prior computation of the covariance structure, which is prohibitive in our application setting.

Linear discriminant analysis over general domains
We start with considering the sub-problem of defining a classifier for the training data {(g i , x i )}, i.e. for the spatially normalized functional predictors supported on the fixed non-linear manifold M. In Section 4, we will then extend the proposed model to account for the geometric component v i in an additive fashion. Assume the training set {(g i , x i )} consists of n independent copies of (G, X), a pair of random variables with X a zero-mean random function taking values in L 2 (M) and G a binary random variable such that P (G = g 1 ) = π 1 , (G = g 2 ) = π 2 . Let µ 1 = E [X|G = g 1 ] and µ 2 = E [X|G = g 2 ] denote the conditional means of X and assume µ 1 = µ 2 . Moreover, let C(p, q) = E [X(p)X(q)] , p, q ∈ M denote the covariance function of X and assume this is square integrable, i.e. M M C(p, q) 2 dpdq < ∞. For a real, symmetric, square integrable, and non-negative function R ∈ L 2 (M × M) let the integral operator L R : L 2 (M) → L 2 (M) be defined as L R (β)(·) = M R(p, ·)β(p)dp, ∀β ∈ L 2 (M). Therefore, with the notation introduced, L C denotes the covariance operator of X, which is a compact self-adjoint operator and therefore admits the following spectral representation in terms of the eigenvalues θ 1 ≥ θ 2 ≥ . . . ≥ 0 and associated eigenfunctions e 1 , e 2 , . . . ⊂ L 2 (M) of L C . Assume that L −1 C (µ 2 − µ 1 ) L 2 (M) < ∞ and define the population quantity β 0 ∈ L 2 (M) such that Note that this is an assumption on the underlying population quantities and will not have practical implications. However, it allows us to have a well-defined population quantity β 0 to study convergence properties of the proposed models. For a discussion on the case L −1 C (µ 2 − µ 1 ) L 2 (M) = ∞, which is related to the perfect classification phenomenon see Delaigle and Hall (2012), Berrendero, Cuevas, and Torrecilla (2018), Chen and Jiang (2018), and Kraus and Stefanucci (2019).
The function β 0 can be understood as a functional analog of the multivariate discriminant vector of a linear discriminant analysis (Shin 2008) and, for a new observation with predictor x * ∈ L 2 (M), can be used to predict the associated label g * with the linear classification rule β 0 , x * L 2 (M) > c th , with c th an appropriately chosen threshold. Moreover, if X is a Gaussian random function within each group in G, it can be shown that the function β 0 defines the linear classifier that minimizes the misclassification error rate and it is therefore optimal (Delaigle and Hall 2012). The discriminant direction β 0 can also be equivalently defined as the minimizer of the functional In practice, the population quantities C, µ 1 , and µ 2 are unknown and need to be estimated from the data. The goal of a classification model is therefore to recover β 0 from the training sample {(g i , x i ) : i = 1, . . . , n} of n independent copies of (G, X). This can be achieved by using the sample covariance LĈ, withĈ (p, q) = 1 n n i=1 x i (p)x i (q), p, q ∈ M and the sample conditional meansμ 1 andμ 2 to replace the population counterparts in (5). An estimatê β of β 0 can then be defined as the minimizer of where a penalty term P(β) is generally added to overcome the ill-posedenss of the minimization problem, which is due to the low-rank structure of LĈ. For instance, Park, Ahn, and Jeon (2021) define a penalty P that encorages the estimateβ to be smooth and sparse, while Kraus and Stefanucci (2019) consider a ridge-type penalty. The functional discriminant model in (6) relies on the precomputation of the empirical covariance function, which is general not possible for dense functional data supported on multidimensional domains and is ultimately not possible in our application setting. We therefore propose a direct regularized approach to estimating β 0 that does not rely on the precomputation of the empirical covariance. This will be possible thanks to the following simple observation. As noted for instance in Delaigle and Hall (2012), the discriminant direction β 0 can be equivalently characterized as the solution to the minimization problem where Y is an auxiliary scalar random variable such that Y = − 1 π 1 if G = g 1 and Y 2 = 1 π 2 otherwise. In other words, the classification problem considered can be reformulated into a functional regression problem. This motivates the adoption of a least-squares approach to estimating β 0 , based on the empirical counterpart of the objective function in (7), where an additional differential regularization term is introduced to incorporate information on the geometric domain M and overcome the illposedness of the problem. In the multivariate setting, analogous least-squares formulations have also been adopted, for instance, in Mai, Zou, and Yuan (2012) and Gaynanova (2020).

Regularized estimation
Given the training sample (g i , x i ) introduce a scalar variable y i such that y i = − n n 1 if g i = 1 and y i = n n 2 otherwise. Note that − n n 1 and n n 2 are an estimate of Y 1 = − 1 π 1 and Y 2 = 1 π 2 , respectively. Let W 2 (M) be the Sobolev space of functions in L 2 (M) with first and second distributional derivatives in L 2 (M). We define an estimateβ ∈ W 2 (M) of the population quantity β 0 as the solution to the following minimization problem β = arg min where the first term is a least-squares estimate of the objective function in (7) and the second term is a differential regularization term. The parameter λ controls the trade-off between the least-squares term of the objective function and the penalty term. Our choice of the regularization term is with ε ≥ 0. This is a linear combination of two terms. The first one is based on the Laplace-Beltrami operator ∆ M : W 2 ⊂ L 2 (M) → L 2 (M) and describes the local smoothness of the function β : M → R on the non-linear manifold domain M. Specifically, it allows the model estimateβ, at any fixed point p ∈ M, to borrow strength from the other points on M while constraining the 'information' to propagate coherently with the non-linear manifold structure of the anatomical object M. The second term is a generic shrinkage-type regularization.

Theory
The aim of this section is to provide theoretical guarantees for the performance of the proposed model. Specifically, we provide a probability bound for the out-of-sample risk, i.e. the random variable where X * is an independent copy of X that is independent of the training data and E * is the expectation taken over X * . Equation (10) measures how close the prediction made with the estimated parameter β is to the 'optimal' one made with the unknown population quantity β 0 . Assume for simplicity that ε > 0. Then, thanks to the Sobolev embedding theorem (Brezis 2011), that is, the evaluation operator is a continuous functional. A direct consequence is that the space is an RKHS with a symmetric, positive definite kernel function K M : M × M → R (Berlinet and Thomas-Agnan 2004). The kernel function K M is used only for theoretical investigation here and obtaining its explicit form is, in general, not computationally feasible and will not be necessary. We will however leverage the fact that L 1 2 K M (L 2 (M)) = W 2 (M) (Cucker and Smale 2002). For ε = 0, J 1/2 defines a semi-norm rather than a norm and similar arguments hold by restricting ourselves to the subspace of L 2 (M) that is orthogonal to the null space of J 1/2 . Next, we define the sandwich operator Cai and Yuan 2012) and make the following assumptions.
Assumption 3.1. The norm of X is a.s. bounded by a positive real number, i.e., there exists M 2 > 0 such that X L 2 (M) ≤ M 2 a.s.
Assumption 3.1 allows us to use a Hoeffding type inequality for Hilbert space valued random elements, i.e., Theorem B.1. Such an assumption does not have practical implications. Assumption 3.2 guarantees that the population quantity β 0 to be estimated is well-defined and belongs to the space of smooth functions W 2 (M). Assumption 3.3 is expressed in terms of properties of the effective dimension D(λ). For our final choice of λ, it is straightforward to check that this assumption holds by assuming that the eigenvalues {τ k } of T decay as τ k k −2r , with r > 1 2 , which is a typical assumption in the functional linear model literature (Tony Cai and Yuan 2012) and is related to the rate of decay of the eigenvalues of L K M and L C , and their alignment.
The following theorem provides an upper bound for the excess prediction risk.
Similar rates of convergence have been shown to hold for regularized estimates in the functional linear regression setting (see e.g. Tony Cai and Yuan 2012;Tong and Ng 2018;Sun et al. 2018;Reimherr, Sriperumbudur, and Taoufik 2018). However, a key difference in our model is that the residual random variable ε = Y − X, β 0 L 2 (M) and the functional predictor X are not independent, preventing the direct application of such results. Therefore, Theorem 3.1 shows that in spite of such a dependence structure we are nevertheless able to recover the functional linear model rates of convergence. The proof is provided in the appendix.

Additive multivariate generalizations
We now consider a bivariate extension of the functional model proposed in the previous section to include the geometric component of the original data. We therefore consider the training sample Assume the training set {(g i , v i , x i )} consists of n independent copies of (G, V, X), a triplet of random variables with V a zero-mean random function taking values in V(R 3 ), X a zero-mean random function taking values in L 2 (M), and G a binary random variable such that P (G = g 1 ) = π 1 , (G = g 2 ) = π 2 .
We now adopt the multivariate functional data notation in Happ and Greven (2018), and define , and assume µ 1 = µ 2 . Next, we define the covariance operator L C : H → H, with the jth element of L C f , f ∈ H given by As in the univariate case, we then assume that the population quantity β 0 ∈ H is well defined and satisfies the equation This can be regarded as a multivariate generalization of the linear discriminant direction defined in the previous section. We now turn to the problem of defining an estimator for β 0 .

Regularized estimation
Let the variable y i be such that y i = − n n 1 if g i = 1 and y i = n n 2 otherwise. We define the multivariate functional estimateβ = β G ,β F of the population quantity β 0 to be the solution to the following minimization problem Model (13) is a multivariate extension of the one proposed in Section 3, where the geometric functional descriptor v i and the functions x i are accounted for in an additive fashion. The remaining terms are regularizations that enforce smoothness on the functional estimatesβ G andβ F in their respective function spaces.

Differential regularization and kernel penalty: a unified modeling perspective
In Section 4.1, we have adopted two different approaches to produce smooth estimatesβ G : R 3 → R 3 andβ F : M → R. The smoothness ofβ F is enforced by means of a penalty J(·) defined in terms of the L 2 (M) norm of a differential operator, which implicitly defines a kernel K M . On the other hand, the smoothness ofβ G is enforced by means of a norm · V(R 3 ) , defined implicitly through the Estimate Function space Norm Kernel direct definition of a kernel K R 3 . For clarity, we summarize the relevant function spaces and associated norms and kernels in Table 1. From a methodological perspective, the reproducing kernel K R 3 (p, q), p, q ∈ R 3 can be understood as a measure of the influence of the function value at p on the function value at q and vice-versa. Defining a smooth function space through its kernel has arguably an advantage when it comes to discretizing an infinite-dimensional minimization problem over that function space. In fact, thanks to the well-known representer theorem (Wahba 1990;Yuan and Cai 2010), under mild conditions, its exact solution can be expressed as a linear combination of the elements of a n-dimensional basis, which involves the explicit expression of the kernel.
Then, it is natural to ask whether a similar approach could be adopted forβ F : M → R. In other words, can we define a smooth real function space on M by explicitly defining a kernel K M : M × M → R that encodes a measure of influence that is coherent with the non-linear geometry of M? This turns out to be a non-trivial task due to the positive definiteness property that K M needs to satisfy. Consider, for instance, the popular exponential kernel. Its natural extension to the manifold setting is K M (p, q) = exp(−c d M (p, q) 2 ), where d M (p, q) is the geodesic distance between p ∈ M and q ∈ M. Unfortunately, this kernel is not guaranteed to be positive definite for a generic non-linear manifold M (Feragen, Lauze, and Hauberg 2015;Jayasumana et al. 2015).
Alternatively, we could try to compute an explicit form of the kernel K M from J(·) by employing the following identity f (p) = K(p, ·), f W 2 (M) , ∀f ∈ W 2 (M), that is, for any fixed p ∈ M, K M (p, ·) ∈ W 2 (M) is given by the solution to the partial differential equation where here δ p denotes the Dirach delta (Wahba 1990;Fasshauer and Ye 2013). However, closed form solutions of (14) are not available in our setting. Approximate solutions could be computed by Finite Elements Analysis (FEA) (Quarteroni 2009). Yet, we would face the issue of having to store the dense object K M (·, ·). As described in Section 5, we will instead leverage FEA to directly discretize the functionβ F : M → R in the objective function (13). This highlights that the choice of the two modeling approaches is not arbitrary and that, arguably, for functional estimates supported on Euclidean spaces, defining explicitly a reproducing kernel is likely the preferred choice, while for general non-Euclidean domains, where defining a reproducing kernel is not trivial, the differential penalization approach is preferable.

Theory
Define the diagonal matrix of reproducing kernels K(p, q) with entries K 11 (p 1 , q 1 ) = K R 3 (p 1 , q 1 ) and K 22 (p 2 , q 2 ) = K M (p 2 , q 2 ) with p i ∈ D i , p j ∈ D j . Let L K : H → H be the associated integral operator and, analogously to the univariate setting, define the sandwich operator T = L K 1 2 CK 1 2 . We make the following assumptions. Assumption 4.2. There exist a smooth function β 0 ∈ V(R 3 ) × W 2 (M) such that β 0 = L −1 C (µ 2 − µ 1 ).
The following theorem, which is an extension of Theorem 3.1, provides an upper bound for the excess prediction risk.
where Y * is an independent copy of Y that is independent of the training data and E * is the expectation taken over Y * .

Discretization
Consider now a triangular mesh M T that is the union of the finite set of triangles T defined over a set of s common vertices ξ 1 , . . . , ξ s . Let this be an approximate representation of the manifold M. We then introduce the linear finite element space W T consisting of a set of globally continuous functions over M T that are affine within each triangle τ in T , i.e.
Let now M and S be the sparse mass and stiffness s × s matrices defined as (M ) jj = M T ψ j ψ j and (S) jj = M T ∇ M T ψ j · ∇ M T ψ j , where ∇ M T is the gradient operator on the manifold M. For β F T of the form given in equation (16), we have that the associated penalty term J(·) can be approximated as (c F ) T D M T c F , with D M T = SM −1 S + εM (Lila, Aston, and Sangalli 2016). Further, following an approach often adopted in FEA (Fried and Malkus 1975;Hinton, Rock, and Zienkiewicz 1976), we replace the dense matrix M −1 with the sparse matrixM −1 , whereM is the diagonal matrix such that M jj = j M jj . This leads to the sparse penalty matrix D M T = SM −1 S + εM . In practice, each functional observation x i is also of the form given in equation (16). Therefore, denoting by X the n × s matrix, where each row contains the basis coefficients of x i , the term x i , β F L 2 (M) can be discretized by XM c F .
We now turn our attention to the estimateβ G ∈ V(R 3 ). Since for V(R 3 ) we have an explicit form of the associated reproducing kernel, we employ the representer theorem (Wahba 1990;Yuan and Cai 2010) and takeβ G of the formβ On the right hand side of Figure 2, we show an example of a function R 3 K R 3 (p, ·)v i (p)dp. Then, we have that β G 2 V = (c G ) T Σc G , where c G = c G 1 , . . . , c G n T and Σ is a n × n matrix with entries . . , s}. This is a scalar linear function within each triangle of the mesh M T that takes value 1 on a fixed vertex and value 0 on every other vertex. On the right hand side, an element of the basis R 3 K R 3 (p, ·)v i (p)dp, i = 1, . . . , n . This is a smooth vector-valued function from R 3 to R 3 .
As a result, the coefficients of the approximate solution of model (13) are given by where y = (y 1 , . . . , y n ) T is the vector of auxiliary response variables. It easy to check that this minimization problem can be equivalently written as the following augmented quadratic least-squares problem Note that for n s, which is the setting of our application, the matrix A is sparse. Therefore, the minimization problem (19) can be efficiently solved by conjugate gradients, or its variations, e.g. LSQR (Paige and Saunders 1982), without requiring explicit computation of the high-dimensional normal matrix A T A -a quantity related to the covariance structure of the functional predictors.
An approximate solution to the univariate model in (8) follows as a special case of the multivariate case considered here.

Data and preprocessing
We analyze a cohort of n = 484 subjects from the ADNI and PPMI studies. On the basis of three molecular biomarkers (the ATN classification scheme (Jack et al. 2016)), each subject was assigned to one of the two diagnostic categories -C: Control (n 1 = 100) and AD: Alzheimer's Disease (n 2 = 384). Here we focus on data collected at the baseline visit, which include, among other imaging modalities, structural T1-weighted MRI.
The T1-weighted images were preprocessed using FreeSurfer Fischl, Sereno, and Dale 1999). Specifically, white matter, grey matter, and cerebrospinal fluid were segmented and used to extract the outer and inner surface of the cerebral cortex. From these two surfaces, we generated a central surface interpolating the midpoints between the outer and inner layers, which offers an accurate representation of the two-dimensional anatomical structure of the cerebral cortex. This representation has the benefit of encoding a notion of distance between brain regions that is neurologically more relevant than the original volumetric representation. The cortical surface can moreover be coupled with one or more maps describing complementary structural or functional properties of the cortex, such as cortical thickness measurements (Fischl and Dale 2000), fMRI signals, or connectivity maps (Smith et al. 2013;Thomas Yeo et al. 2011). In this study, we focus on cortical thickness, which is estimated from the distances between the outer and inner surface of the cerebral cortex. Next, the n surfaces were registered and sub-sampled.
Hence, from the preprocessing stage, we obtain n = 484 triangular meshes M T i consisting of s = 64K vertices in correspondence across subjects and a set of triangular faces defining how these vertices are connected to define the surfaces. By classical Generalized Procrustes Analysis (Dryden and Mardia 2016), translation, rigid rotation, and scale were removed from the data, while jointly estimating a template M T , which is also a triangular mesh with s = 64K vertices in correspondence with those of the individual subjects. Each surface has been moreover coupled with cortical thickness measurements at the mesh vertices, which we model as a real piecewise linear continuous function z T i over the mesh M T i . The preprocessing stage results in a set FoSs M T i , z T i , i = 1, . . . , n , which are discretized versions of the continuous idealized objects {(M i , z i ) , i = 1, . . . , n} introduced in Section 2. To simplify the notation, we drop the superscript T . Moreover, we denote the diagnostic labels by {g i ∈ {C, AD}, i = 1, . . . , n}. Three examples of the such FoSs, and associated diagnostic labels, are given in Figure 1A.
Here, we are interested in identifying subjects with AD using the FoSs extracted from the MRI images. A key feature of the proposed analysis is interpretability. While it is important to build models with good classification accuracy, it is also important to describe what the estimated relationship between the predictors and the outcome variable is, in order to inform studies and generate datadriven hypotheses about the pathophysiology of AD.

Linear functional representations
For each FoS, we compute a function v i ∈ V(R 3 ) such that ϕ v i (M) is a close approximation of M i , where closeness is measured as the Euclidean distance between the corresponding vertices of ϕ v i (M) and M i . We can then transport the function z i : M i → R onto the template defining a continuous piecewise linear function x i = z i • ϕ v i . This leads to the definition of the bivariate functional data (v i , x i ) that are linear representations of the FoSs ( Further details on the computation of v i can be found in the appendices.

Discriminant analysis
Our aim is to estimate directions in the parameterizing geometric and thickness spaces that are most predictive of AD. To this end, we apply the model introduced in Section 4.1 to the training data {g i , v i −v, x i −x}, wherev andx are the sample means of {v i } and {x i }. The training data are a subset of the entire dataset containing 60% of the observations. This leads to the functional estimateŝ β G : R 3 → R 3 andβ F : M → R. Given a new subject with predictors (v * , x * ), these estimates can be used in conjunction with the classification rule v * −v,β G + x * −x,β F > c th to predict the diagnostic label of a new subject. The cut-off level c th can be chosen by computing sensitivity and specificity on a test set, for different values c th , and selecting the desired level and type of accuracy.
These estimates effectively identify linear directions {c 1β G , c 1 ∈ R} and {c 2β F , c 2 ∈ R}, in their respective spaces, that can be interpreted as the most discriminant geometric and thickness directions. Specifically, large values of c 1 ∈ R and c 2 ∈ R describe configurations c 1β G and c 2β F that are predictive of AD. Low values of c 1 ∈ R and c 2 ∈ R describe configurations that are instead predictive These are a vector fieldβ G : R 3 → R 3 , representing the most predictive geometric pattern of AD, and a functionβ F : M → R, representing the most predictive cortical thickness pattern of AD. For a new FoS, with linear representation (v * , x * ), we compute the score v * −v,β G + x * −x,β F and predict whether the subject has AD by comparing the score value with a predetermined threshold c th . On the right hand side, we depict the process of mapping back the estimatesβ G andβ F onto the space of FoSs. On the same space, we also pictorially map the classification rule adopted. The blue areas in theβ F figure indicate what are the parts of the cortical surface where a thinner cortex, than the average population, is predictive of AD. These are mostly localized in the lateral temporal, entorhinal, inferior parietal, precuneus, and posterior cingulate cortices. The red arrows in theβ G are an indication of what are the areas where differences in the morphological configuration of the cerebral cortex, with respect to the average population, are most predictive of AD. The specific types of morphological change can be inspected by comparing the surfaces ϕv −c 1β G (M) and ϕv +c 1β G (M), on the right hand side diagram. of the subject being healthy. Moreover, given the additive modeling assumption on the geometric and thickness components, for every configuration c 1β G , c 2β F , an increase of c 1 (c 2 ) for a fixed c 2 (c 1 ) describes a configuration that is more strongly associated with AD. These linear trajectories on the parametrizing space can be mapped back onto the original space of FoSs by using equation (2), i.e., by defining the curved space We fit the proposed model for different choices of the parameters λ 1 and λ 2 . Recall that λ 1 controls the regularity of the geometric estimate and λ 2 that of the thickness, with high values virtually constraining the solution to be the zero function. The final choice of λ 1 and λ 2 are the results of a compromise between classification accuracy on the test set (final test AUC = 0.7006) and compatibility of the estimated discriminant directions with the neurodegenerative nature of the disease (see also Discussion). For comparison, we also fit the proposed model for different choices of the parameters λ 2 Figure 4: A comparison of the cortical thickness discriminant direction obtained by fitting Model (8) (i.e., without modeling geometric differences) with that obtained by fitting Model (13) (i.e. by explicitly modeling inter-subject geometric differences). The figure on the left indicates that a thicker cortex in the red areas is associated with a higher risk of AD, which is not physiologically feasible. Once we explicitly model geometric differences, this evidence seems to disappear. This likely means that there is a non-negligible dependence structure between the predictors modeling geometry and those modeling thickness. What appeared to be differences in the thickness of the cerebral cortex in the model without geometric component are now captured by the term modeling geometric differences. Furthermore, it is only when we model inter-subject geometric differences that entorhinal cortex atrophy in the medial temporal lobe is identified as the strongest predictor of AD. This is consistent with pathological findings and staging of early AD (Braak et al. 2006). and with λ 1 = +∞, i.e., discarding the geometric component.

Results
On the left hand side of Figure 3, we show the estimated most discriminant geometric and thickness directions, i.e.,β G : R 3 → R 3 andβ F : M → R. These have been estimated by applying the model in (13) to the linear representations {(v i −v, x i −x)}. The colormap describingβ F indicates what types of variations with respect to the population average cortical thickness are most predictive of AD. Specifically, a thinner cerebral cortex in the blue areas (i.e., lateral temporal, entorhinal, inferior parietal, precuneus, and posterior cingulate cortices) is associated with a higher risk of AD. Therefore, our study identifies the typical thickness signature of AD observed to date (see, e.g., Bondareff et al. 1989;Dickerson et al. 2009;Sabuncu et al. 2011). The geometric componentβ G is instead a vector field in R 3 . This is a linear representation of the types of morphological variations, with respect to the population average cortex geometry, that are associated with a higher risk of AD. While a full understanding of its information content is only possible by mappingβ G back to the space of FoSs, i.e., exploring ϕv +c 1β G (M) for different choices of c 1 ∈ R, the magnitude of the vector fieldβ G , at any fixed point, gives a rough indication of what are the parts of the cortex whose morphological variations are most relevant to the classification problem.
On the right hand side of Figure 3, we show the FoSs associated with the linear representationsβ G andβ F , that is, ϕv +c 1β G (M), c 2β F • ϕ −1 v+c 1β G with c 1 , c 2 ∈ R. These describe the most predictive patterns in terms of the original neurobiological objects. We have circled a specific area of the brain to ease comparison and appreciate the specific morphological patterns that we estimate being relevant to the classification problem.

Discussion
The methodological approach here introduced is an important addition to current neuroimaging methods. In fact, these are mostly based on cortical thickness analyses on pre-specified parcellations to reduce the dimension of the cortical thickness maps, which unavoidably constrains the final results to the specific choice of the parcellation. Moreover, such models usually discard the morphological information of the cerebral cortex, which we have shown could lead to misleading results. Our approach aims to overcome both limitations.
As mentioned, our analysis identified the typical AD thickness pattern. Several studies that focus on identifying AD-vulnerable areas include the regions found in our analysis (see, e.g., Bondareff et al. 1989;Dickerson et al. 2009;Sabuncu et al. 2011). However, there is some variability in the estimated regions. For instance, Sabuncu et al. (2011) used a dynamic model and found strongest changes in the inferior parietal regions and the posterior cingulate. It should be noted that these studies typically rely on massive univariate analyses (between the cortical thickness at each vertex, or each parcel, and the diagnostic label) and show the resulting p-value maps. They are therefore taking a feature-centric perspective on the problem (the features here being the thickness at each location). It is not clear how these findings would generalize to the prediction problem on previously unobserved data (Li and Tong 2020). In contrast, we formulate the problem as a linear classification problem, which is tested against an unobserved test set, yet yielding interpretable classification rules thanks to a linear regularized approach.
Moreover, to demonstrate the importance of modeling cortical geometry, we compare our results to those obtained by fitting model (8), i.e., by discarding inter-subject geometric differences. We compare the estimates of the two models in Figure 4. What we observe is that, with the exception of the entorhinal cortex, the simpler model is able to capture the main areas whose cortex thinning is associated with AD. However, the model appears also to estimate an association between an increase in cortical thickness in certain areas (dark red) with higher risk of AD. Given the neurodegenerative nature of AD this is not physiologically feasible. Interestingly, introducing the geometric component in the model reduces such effects. This may also be caused by the geometric component now capturing systematic misregistration. In order to verify such a hypothesis, further validation of the estimated geometric component is required in controlled settings where the registration is more reliable, e.g., in the longitudinal setting.

Conclusions
We introduce a framework for the discriminant analysis of functional data supported on random manifold domains, i.e., FoSs. To this aim, we adopt linear representations of these objects that are bivariate functional data belonging to linear spaces. We then define a functional linear classification model on the parametrizing space. Thanks to a penalized least-squares formulation, the proposed model is able to estimate the most discriminant direction in the data without requiring the explicit computation of the covariance function of the predictors or low-rank approximations thereof. This allows us to reduce the memory requirements by five orders of magnitude and ultimately be able to run our model on a standard workstation. The complexity of the solution is controlled by means of differential penalties that are aware of the geometry of the domain where the functional data are supported.
We apply the proposed model to the analysis of modern multi-modal neuroimaging data. Specifically, we estimate interpretable discriminant directions that are able to exploit both geometric and thickness features of the cerebral cortex to identify subjects with AD. We show that our results are consistent with those in the neuroscience literature.
The model proposed can be applied as-is to different imaging settings that lead to FoSs representations, such as musculoskeletal imaging (Gee, Treece, and Poole 2018) or cardiac imaging (Biffi et al. 2018). It is also important to highlight that the proposed model is not a mere generalization of existing models for functional data supported on one-dimensional domains to multidimensional domains. We indeed believe that its application to one-dimensional functional data, where the bivariate representation is given by the registered functions and associated registration maps, leads to a novel classification approach in such a simplified setting.

Acknowledgments
Part of the data collection and sharing for this project was funded by the Alzheimer's Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: AbbVie, Alzheimer's Association; is a diffeomorphic deformation of R 3 (see, e.g., Younes 2019). We then model {v t : t ∈ [0, 1]} to be a minimizer of the quantity 1 0 v t 2 V dt, for a given initial vector field v 0 ∈ V(R 3 ) (Miller, Trouvé, and Younes 2006). Finally, the diffeomorphic operator is defined as is the initial vector field generating {v t : t ∈ [0, 1]}, and φ v is the solution of the ODE (20) for the computed {v t : t ∈ [0, 1]}.

A.2 Computation
In practice, each surface M i has a computational representation M T i that is a triangular mesh with s vertices ξ i 1 , . . . , ξ i s in correspondence across the n subjects. We use these vertices to perform Procrustes analysis (Dryden and Mardia 2016) and remove translation, size, and rigid rotations from the surfaces M T i . In addition, Procrustes analysis yields a template M T with vertices ξ 1 , . . . , ξ s . The linear representation functions {v i : v i ∈ V(R 3 )}, associated with the surfaces M T i , are then computed by solving the minimization problem where the least-squares term ensures that the deformed template is a regularizing term that encourages the solution to achieve such a close approximation with a 'minimal' deformation. The constant λ is selected by inspecting the solutions on a small subset of the full cohort. To perform the actual computations, we use the MATLAB implementation fshapetk (Charlier, Nardi, and Trouvé 2015;Charlier, Charon, and Trouvé 2017). Note that if the vertices were not in correspondence, that is, the surfaces had not been registered beforehand, the proposed framework would still be applicable by replacing the least-squares term in (21) with a more general shape similarity measure, an example of which is given in Vaillant and Glaunès (2005) and Vaillant et al. (2007).

B Proofs
Recall that we denote by C the covariance function of the functional predictors and byĈ its empirical counterpart. Moreover, K M denotes the evaluation kernel of the Sobolev space W 2 (M) endowed with K M . To simplify the notation, we drop the subscript from · L 2 (M) and denote it by · . We define A op = sup f : f =1 Af to be the operator norm of a linear operator A : L 2 (M) → L 2 (M).
Thanks to the fact that L 1 2 K M (L 2 (M)) = W 2 (M) (Cucker and Smale 2002), it is clear that we can reformulate problem (8) as whose solutionf λ is given byf and where T n = L 1 2 can be easily Let moreoverd = 1 n n i=1 y i x i and d = µ 2 − µ 1 , and notice that E[d] = d.
Then, we have = I 1 + I 2 We first derive a bound for the term I 1 in (27), i.e., the variance term, and then proceed with bounding the term I 2 , i.e., the bias term. To do so, we will also use Theorems B.2-B.3, stated in Section B.4.

B.1 Variance
By simple calculations we have that where in the last inequality we have used A γ B γ op ≤ AB γ op , for any 0 < γ < 1 (Blanchard and Krämer 2010). We can then bound the first term thanks to Theorem B.3. We therefore turn to the second term.
First observe that and, therefore, with {τ k } and {η k } denoting the eigenvalues and eigenfunctions of T , respectively.
We then have Moreover, by appealing to the Markov inequality, we have that with confidence at least 1 − δ/2 with B n,λ = 2κ √ n κ √ nλ + D(λ) and κ = M 2 L 1 2 K M op . Thanks to the inequality in (28) and Theorem B.3, we have that with probability at least 1 − δ We now turn to the bias term I 2 in (27).

B.2 Bias
By simple calculations, we have Moreover, by using Theorems B.2-B.3, and the inequality we have that with probability at least 1 − δ T 1 2

B.3 Final rates
Using the variance and bias bounds in equations (29) and (31), we get the following bound for the out-of-sample risk: where in the last inequality we have chosen λ = n − 1 1+θ and used the simple inequality B n,λ ≤ 2κ(κ + √ c) √ λ. The constant c here is from Assumption 3.3. This implies the result stated in Theorem 3.1.

B.4 Auxiliary results
Theorem B.1. Let H be a Hilbert space endowed with a norm · H and let X be a random variable taking values in H. Let ξ 1 , . . . , ξ n be a sequence of n independent copies of X. Assume that ξ H ≤ M (a.s.), then for 0 < δ < 1 2M log(2/δ) n + 2E ξ 2 H log(2/δ) n with probability at least 1 − δ.
Proof. A proof can be found in Pinelis (2007).
Using Theorem B.1 it can be shown that the following two theorems hold.
Proof. A proof can be found in Tong and Ng (2018).
Proof. A proof can be found in Tong and Ng (2018).

B.5 Multivariate model
Thanks to the multivariate notation introduced, we can reformulate the multivariate model (13) as with x i := (v i , x i ). Therefore, the proof of Theorem 4.1 follows the same lines as that of Theorem 3.1 and is therefore omitted.

C Simulations
In this section, we conduct simulations to assess the finite sample classification performance of the model proposed (FLDA) when compared with FPCA followed by multivariate LDA (Hastie, Tibshirani, and Friedman 2009) on the PC scores (FPCA+LDA). The main purpose is to highlight the importance of developing classification models that do not rely on unsupervised dimension reductions of the functional predictors, which could remove relevant information. Specifically, we focus on the univariate model proposed in Section 3. This will allow us to adopt the geometry-aware FPCA model proposed in Lila, Aston, and Sangalli (2016) to perform dimension reduction in the FPCA+LDA approach.
We use a triangulated surface M T with 642 nodes that is an approximation of the brainstem. On this triangulated surface we generate the orthonormal functions {v l : l = 1, 2, . . . , 40} consisting of 40 eigenfunctions of the Laplace-Beltrami operator. Then, we generate two sets of smooth functional data supported on M T as follows: Group 1: x i1 = w i1 v 1 + w i2 v 2 + . . . + w i40 v 40 i = 1, . . . , n Group 2: x i2 = αµ + u i1 v 1 + u i2 v 2 + . . . + u i40 v 40 i = 1, . . . , n , where u ij and w ij are zero-mean independent random variables that represent the scores, and are distributed according to a normal with variance σ 2 j decreasing in j. The function µ is the group difference, chosen to be the 10th eigenfunction, and its magnitude is controlled by a parameter α > 0. The parameter α controls the 'difficulty' of the problem. We then identify three regimes, with α ∈ 0.2, 0.4, 0.6. Given a new function x * , the aim is to recover the group this belongs to.
We compare the proposed FLDA technique with FPCA+LDA. To choose the models' hyperparameters, we use a validation set approach. We evaluate the performance for different sample sizes of the training data: n = 128, 256, 512 and 1024. For every n, besides the training set, we generate a validation set of size n and a test set including 20K samples, and repeat the experiment 50 times. We summarize the classification performance on the test set in Figure 5.
In spite of the simulation setting being 'fair' to FPCA+LDA, as the discriminant direction was chosen to be the 10th eigenfunction -and therefore easy to capture with a small number of eigenfunctions with respect to the number of observations -the results show that for α = 0.2 and 0.4, FLDA tends to perform better than FPCA+LDA. FPCA+LDA performs rather poorly when the difference in the means is small, i.e. the discrimination task is more difficult. When α = 0.6, the two approaches achieve similar performances.  Figure 5: Area Under the Curve (AUC) scores of the two classification approaches considered when applied to a test set. In the left and middle panels, we show the results for α = 0.2 and α = 0.4, i.e., small and medium differences in the group means. These two settings highlight the fact that the proposed approach performs better than FPCA+LDA when the classification task is challenging. In the rightmost panel, we show the results for α = 0.6, i.e. large differences in the means of the two groups. This setting shows that when the classification task becomes simpler the two approaches give comparable results. For each setting, and for every choice of n ∈ {128, 256, 512, 1024}, 50 datasets were generated.