# Interpretable discriminant analysis for functional data supported on random
non-linear domains

Eardi Lila  Department of Biostatistics, University of Washington  Wenbo Zhang
Department of Biostatistics, University of Washington  Department of
Statistics, University of California, Irvine  Swati Rane  Department of
Radiology, University of Washington
for the Alzheimer’s Disease Neuroimaging Initiative111Data used in preparation
of this article were obtained from the Alzheimer’s Disease Neuroimaging
Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators
within the ADNI contributed to the design and implementation of ADNI and/or
provided data but did not participate in analysis or writing of this report. A
complete listing of ADNI investigators can be found at:
http://adni.loni.usc.edu/wp-
content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf

###### Abstract

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.

##  1 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) [46, 20, 7, 38] and has motivated
a large body of literature [[, see
e.g.,]]James2001,Muller2005b,Preda2007a,Delaigle2012,Dai2017,Berrendero2018,Kraus2019a,Park2021.
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.,]]Su2014,Dai2018,Lin2018,Dubey2020,Kim2018.

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) [72]. 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 [4, 8, 66, 47].

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 [41, 53, 60, 68, 77], 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.

Figure 1: Panel A: FoSs of three subjects in the training sample, where
$g_{i}\in\\{\text{`C'},\text{`AD'}\\}$ denotes the disease state of the $i$th
individual (C: Control, AD: Alzheimer’s Disease), $\mathcal{M}_{i}$ is a two-
dimensional manifold encoding the geometry of the cerebral cortex, and
$z_{i}:\mathcal{M}_{i}\rightarrow\mathbb{R}$ is a real function, supported on
$\mathcal{M}_{i}$, describing cortical thickness (in mm). Panel B: Linear
representation $\left(v_{i},x_{i}\right)$ of each FoS
$\left(\mathcal{M}_{i},z_{i}\right)$ shown in Panel A. Here
$v_{i}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ is a vector-valued function
encoding the geometry of the $i$th individual. This is depicted as a
collection of 3D vectors $\\{v_{i}(p_{j})\\}$ for a dense set of points
$\\{p_{j}\\}\subset\mathbb{R}^{3}$. For clarity, the function $v_{i}$ is
displayed only on half of its domain $\mathbb{R}^{3}$. The function
$x_{i}:\mathcal{M}\rightarrow\mathbb{R}$ describes the spatially normalized
thickness map of the $i$th individual on the fixed template $\mathcal{M}$.
Panel C: FoS $(\varphi_{v_{i}}(\mathcal{M}),x_{i}\circ\varphi^{-1}_{v_{i}})$
parameterized by the associated functions $\left(v_{i},x_{i}\right)$ in Panel
B. These is a close approximation of the FoS
$\left(\mathcal{M}_{i},z_{i}\right)$ in Panel A.

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 [13]. 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.]]Chung2014, Lila2016, Mejia2019. 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 [40, 69].

Little attention has been given to the most general setting of functional data
supported on random domains, with the exception of [50, 72], where the authors
focus on unsupervised models, and [52], where the authors leverage the fshape
framework [50] 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.

##  2 Functional data supported on general random domains

The data motivating this work are a sample of triplets

$\left\\{(g_{i},\mathcal{M}_{i},z_{i}),i=1,\ldots,n\right\\},$ (1)

where $g_{i}$ is a binary label, $\mathcal{M}_{i}\subset\mathbb{R}^{3}$ is a
sample-specific closed two-dimensional manifold embedded in $\mathbb{R}^{3}$,
and $z_{i}:\mathcal{M}_{i}\rightarrow\mathbb{R}$ is a scalar function
supported on the geometric object $\mathcal{M}_{i}$. We moreover assume that
the points on the geometries $\\{\mathcal{M}_{i}\\}$ of the observed FoSs are
in one-to-one correspondence across subjects. We refer to the pairs
$\\{(\mathcal{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,
$\mathcal{M}_{i}$ encodes the geometry of the cerebral cortex, and $z_{i}$
encodes the cortical thickness map supported on $\mathcal{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 $(\mathcal{M}^{*},z^{*})$.

###  2.1 Linear functional representations

In our motivating application, the geometric objects
$\\{\mathcal{M}_{i}:i=1,\ldots,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
$\\{(\mathcal{M}_{i},z_{i}):i=1,\ldots,n\\}$ in terms of objects belonging to
linear function spaces.

Let $\mathcal{M}$ be a template two-dimensional manifold embedded in
$\mathbb{R}^{3}$ that is topologically equivalent to a sphere. We denote by
$\mathcal{L}^{2}(\mathcal{M})$ the space of square integrable functions over
$\mathcal{M}$, equipped with the standard inner product
$\langle\cdot,\cdot\rangle_{\mathcal{L}^{2}(\mathcal{M})}$ and norm
$\|\cdot\|_{\mathcal{L}^{2}(\mathcal{M})}$ and denote by
$\mathcal{L}^{2}(\mathbb{R}^{3},\mathbb{R}^{3})$ the space of square
integrable vector-valued functions from $\mathbb{R}^{3}$ to $\mathbb{R}^{3}$,
with inner product
$\langle\cdot,\cdot\rangle_{\mathcal{L}^{2}(\mathbb{R}^{3},\mathbb{R}^{3})}$
and norm $\|\cdot\|_{\mathcal{L}^{2}(\mathbb{R}^{3},\mathbb{R}^{3})}$. Let
$\mathcal{V}(\mathbb{R}^{3})\subset\mathcal{L}^{2}(\mathbb{R}^{3},\mathbb{R}^{3})$
be a Reproducing Kernel Hilbert Space (RKHS) of smooth functions with compact
support in $\mathbb{R}^{3}$. We then introduce a diffeomorphic operator
$\varphi$ such that $\varphi_{v}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ is a
diffeomorphic function for every choice of $v\in\mathcal{V}(\mathbb{R}^{3})$
[16]. We denote by $\varphi_{v}(\mathcal{M})$ the geometric object given by
displacing every point $p\in\mathcal{M}\subset\mathbb{R}^{3}$ to the new
location $\varphi_{v}(p)$. A direct consequence of $\varphi_{v}$ being
diffeomorphic is that $\varphi_{v}(\mathcal{M})$ is topologically equivalent
to a sphere for every choice of $v\in\mathcal{V}(\mathbb{R}^{3})$. The
construction of the diffeomorphic operator adopted in this paper and the
computation of $v_{i}\in\mathcal{V}(\mathbb{R}^{3})$ such that
$\varphi_{v_{i}}(\mathcal{M})\approx\mathcal{M}_{i}$ are detailed in Appendix
A.1.

We can then represent each FoS $(\mathcal{M}_{i},z_{i})$ with a pair of
functions $(v_{i},x_{i})$ such that

$(\mathcal{M}_{i},z_{i})=(\varphi_{v_{i}}(\mathcal{M}),x_{i}\circ\varphi^{-1}_{v_{i}}),$
(2)

where $v_{i}\in\mathcal{V}(\mathbb{R}^{3})$ and
$x_{i}\in\mathcal{L}^{2}(\mathcal{M})$. Intuitively, this means that every
geometric object $\mathcal{M}_{i}$ is modeled as a (diffeomorphic) deformation
of the template, i.e. $\varphi_{v_{i}}(\mathcal{M})$, while the associated
function $z_{i}$ is given by ‘transporting’ the signal $x_{i}$ onto
$\mathcal{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},\mathcal{M}_{i},z_{i})|i=1,\ldots,n\\}$ into the problem of
learning a classifier from

$\left\\{(g_{i},v_{i},x_{i})|i=1,\ldots,n\right\\},$ (3)

where
$v_{i}\in\mathcal{V}(\mathbb{R}^{3})\subset\mathcal{L}^{2}(\mathbb{R}^{3},\mathbb{R}^{3})$
and $x_{i}\in\mathcal{L}^{2}(\mathcal{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
$\beta^{G}\in\mathcal{V}(\mathbb{R}^{3}),\beta^{F}\in\mathcal{L}^{2}(\mathcal{M})$,
such as a ‘direction’ that best discriminates between two classes, can be
mapped onto the original space, defining the associated trajectories of FoSs

$\left\\{\left(\varphi_{c_{1}\beta^{G}}({\mathcal{M}_{0}}),c_{2}\beta^{F}\circ\varphi_{c_{1}\beta^{G}}^{-1}\right),c_{1},c_{2}\in\mathbb{R}\right\\},$

where $\varphi_{c_{1}\beta^{G}}({\mathcal{M}_{0}})$ is guaranteed to be
topologically equivalent to a sphere, and therefore, to satisfy the physical
constraints of the problem considered.

###  2.2 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,\ldots,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.

##  3 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
$\mathcal{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 $\mathcal{L}^{2}(\mathcal{M})$ and $G$ a binary
random variable such that
$P\left(G=g_{1}\right)=\pi_{1},\left(G=g_{2}\right)=\pi_{2}$. Let
$\mu_{1}=\mathbb{E}\left[X|G=g_{1}\right]$ and
$\mu_{2}=\mathbb{E}\left[X|G=g_{2}\right]$ denote the conditional means of $X$
and assume $\mu_{1}\neq\mu_{2}$. Moreover, let
$C(p,q)=\mathbb{E}\left[X(p)X(q)\right],p,q\in\mathcal{M}$ denote the
covariance function of $X$ and assume this is square integrable, i.e.
$\int_{\mathcal{M}}\int_{\mathcal{M}}C(p,q)^{2}dpdq<\infty$. For a real,
symmetric, square integrable, and non-negative function
$R\in\mathcal{L}^{2}(\mathcal{M}\times\mathcal{M})$ let the integral operator
$L_{R}:\mathcal{L}^{2}(\mathcal{M})\rightarrow\mathcal{L}^{2}(\mathcal{M})$ be
defined as

$L_{R}(\beta)(\cdot)=\int_{\mathcal{M}}R(p,\cdot)\beta(p)dp,\,\forall\beta\in\mathcal{L}^{2}(\mathcal{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

$L_{C}(\beta)=\sum_{k=1}^{\infty}\theta_{k}\langle\beta,e_{k}\rangle e_{k},$
(4)

in terms of the eigenvalues $\theta_{1}\geq\theta_{2}\geq\ldots\geq 0$ and
associated eigenfunctions
$e_{1},e_{2},\ldots\subset\mathcal{L}^{2}(\mathcal{M})$ of $L_{C}$.

Assume that
$\|L_{C}^{-1}(\mu_{2}-\mu_{1})\|_{\mathcal{L}^{2}(\mathcal{M})}<\infty$ and
define the population quantity $\beta^{0}\in{\mathcal{L}^{2}(\mathcal{M})}$
such that

$L_{C}\beta^{0}=\mu_{2}-\mu_{1}.$

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 $\beta^{0}$ to study convergence properties of the
proposed models. For a discussion on the case
$\|L_{C}^{-1}(\mu_{2}-\mu_{1})\|_{\mathcal{L}^{2}(\mathcal{M})}=\infty$, which
is related to the perfect classification phenomenon see [34, 54, 56, 64].

The function $\beta_{0}$ can be understood as a functional analog of the
multivariate discriminant vector of a linear discriminant analysis [25] and,
for a new observation with predictor $x^{*}\in\mathcal{L}^{2}(\mathcal{M})$,
can be used to predict the associated label $g^{*}$ with the linear
classification rule
$\langle\beta^{0},x^{*}\rangle_{\mathcal{L}^{2}(\mathcal{M})}>c^{\text{th}}$,
with $c^{\text{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 $\beta^{0}$ defines the linear classifier that minimizes the
misclassification error rate and it is therefore optimal [34]. The
discriminant direction $\beta^{0}$ can also be equivalently defined as the
minimizer of the functional

$\frac{1}{2}\langle\beta,L_{C}\beta\rangle-\langle\mu_{2}-\mu_{1},\beta\rangle.$
(5)

In practice, the population quantities $C$, $\mu_{1}$, and $\mu_{2}$ are
unknown and need to be estimated from the data. The goal of a classification
model is therefore to recover $\beta_{0}$ from the training sample
$\\{(g_{i},x_{i}):i=1,\ldots,n\\}$ of $n$ independent copies of $(G,X)$. This
can be achieved by using the sample covariance $L_{\hat{C}}$, with

$\hat{C}(p,q)=\frac{1}{n}\sum_{i=1}^{n}x_{i}(p)x_{i}(q),\qquad
p,q\in\mathcal{M}$

and the sample conditional means $\hat{\mu}_{1}$ and $\hat{\mu}_{2}$ to
replace the population counterparts in (5). An estimate $\hat{\beta}$ of
$\beta^{0}$ can then be defined as the minimizer of

$\frac{1}{2}\langle\beta,L_{\hat{C}}\beta\rangle-\langle\hat{\mu}_{2}-\hat{\mu}_{1},\beta\rangle+\mathcal{P}(\beta),$
(6)

where a penalty term $\mathcal{P}(\beta)$ is generally added to overcome the
ill-posedenss of the minimization problem, which is due to the low-rank
structure of $L_{\hat{C}}$. For instance, [76] define a penalty $\mathcal{P}$
that encorages the estimate $\hat{\beta}$ to be smooth and sparse, while [64]
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 $\beta^{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 [34], the discriminant direction
$\beta^{0}$ can be equivalently characterized as the solution to the
minimization problem

$\beta^{0}=\operatorname*{arg\,min}_{\beta}\mathbb{E}\left[Y-\langle
X,\beta\rangle_{\mathcal{L}^{2}(\mathcal{M})}\right]^{2},$ (7)

where $Y$ is an auxiliary scalar random variable such that
$Y=-\frac{1}{\pi_{1}}$ if $G=g_{1}$ and $Y_{2}=\frac{1}{\pi_{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 $\beta^{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 $\mathcal{M}$ and overcome the ill-posedness of the problem. In the
multivariate setting, analogous least-squares formulations have also been
adopted, for instance, in [35, 70].

###  3.1 Regularized estimation

Given the training sample $(g_{i},x_{i})$ introduce a scalar variable $y_{i}$
such that $y_{i}=-\frac{n}{n_{1}}$ if $g_{i}=1$ and $y_{i}=\frac{n}{n_{2}}$
otherwise. Note that $-\frac{n}{n_{1}}$ and $\frac{n}{n_{2}}$ are an estimate
of $Y_{1}=-\frac{1}{\pi_{1}}$ and $Y_{2}=\frac{1}{\pi_{2}}$, respectively. Let
$\mathcal{W}^{2}(\mathcal{M})$ be the Sobolev space of functions in
$\mathcal{L}^{2}(\mathcal{M})$ with first and second distributional
derivatives in $\mathcal{L}^{2}(\mathcal{M})$. We define an estimate
$\hat{\beta}\in\mathcal{W}^{2}(\mathcal{M})$ of the population quantity
$\beta^{0}$ as the solution to the following minimization problem

$\displaystyle\hat{\beta}=\operatorname*{arg\,min}_{\beta\in\mathcal{W}^{2}(\mathcal{M})}\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}-\langle
x_{i},\beta\rangle_{\mathcal{L}^{2}(\mathcal{M})}\right)^{2}+\lambda
J(\beta),$ (8)

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
$\lambda$ controls the trade-off between the least-squares term of the
objective function and the penalty term. Our choice of the regularization term
is

$J(\beta)=\|\Delta_{\mathcal{M}}\beta\|^{2}_{\mathcal{L}^{2}(\mathcal{M})}+\varepsilon\|\beta\|^{2}_{\mathcal{L}^{2}(\mathcal{M})},$
(9)

with $\varepsilon\geq 0$. This is a linear combination of two terms. The first
one is based on the Laplace-Beltrami operator
$\Delta_{\mathcal{M}}:\mathcal{W}^{2}\subset\mathcal{L}^{2}(\mathcal{M})\rightarrow\mathcal{L}^{2}(\mathcal{M})$
and describes the local smoothness of the function
$\beta:\mathcal{M}\rightarrow\mathbb{R}$ on the non-linear manifold domain
$\mathcal{M}$. Specifically, it allows the model estimate $\hat{\beta}$, at
any fixed point $p\in\mathcal{M}$, to borrow strength from the other points on
$\mathcal{M}$ while constraining the ‘information’ to propagate coherently
with the non-linear manifold structure of the anatomical object $\mathcal{M}$.
The second term is a generic shrinkage-type regularization.

###  3.2 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

$\mathbb{E}^{*}\left[\langle X^{*},\beta^{0}-\hat{\beta}\rangle\right]^{2},$
(10)

where $X^{*}$ is an independent copy of $X$ that is independent of the
training data and $\mathbb{E}^{*}$ is the expectation taken over $X^{*}$.
Equation (10) measures how close the prediction made with the estimated
parameter $\hat{\beta}$ is to the ‘optimal’ one made with the unknown
population quantity $\beta^{0}$.

Assume for simplicity that $\varepsilon>0$. Then, thanks to the Sobolev
embedding theorem [31], $\exists M\geq 0$ such that for any $p\in\mathcal{M}$

$f(p)\leq\sup_{q}|f(q)|\leq
M\left(\|\Delta_{\mathcal{M}}f\|^{2}_{\mathcal{L}^{2}(\mathcal{M})}+\varepsilon\|f\|^{2}_{\mathcal{L}^{2}(\mathcal{M})}\right)^{1/2},\qquad\forall
f\in\mathcal{W}^{2}(\mathcal{M}),$

that is, the evaluation operator is a continuous functional. A direct
consequence is that the space $\mathcal{W}^{2}(\mathcal{M})$ equipped with the
norm
$J^{1/2}(\cdot)=\left(\|\Delta_{\mathcal{M}}\cdot\|^{2}_{\mathcal{L}^{2}(\mathcal{M})}+\varepsilon\|\cdot\|^{2}_{\mathcal{L}^{2}(\mathcal{M})}\right)^{1/2}$
is an RKHS with a symmetric, positive definite kernel function
$K_{\mathcal{M}}:\mathcal{M}\times\mathcal{M}\rightarrow\mathbb{R}$ [15]. The
kernel function $K_{\mathcal{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_{K_{\mathcal{M}}}^{\frac{1}{2}}(\mathcal{L}^{2}(\mathcal{M}))=\mathcal{W}^{2}(\mathcal{M})$
[14]. For $\varepsilon=0$, $J^{1/2}$ defines a semi-norm rather than a norm
and similar arguments hold by restricting ourselves to the subspace of
$\mathcal{L}^{2}(\mathcal{M})$ that is orthogonal to the null space of
$J^{1/2}$.

Next, we define the sandwich operator
$T=L_{K_{\mathcal{M}}^{\frac{1}{2}}CK_{\mathcal{M}}^{\frac{1}{2}}}=L_{K_{\mathcal{M}}}^{\frac{1}{2}}\circ
L_{C}\circ L_{K}^{\frac{1}{2}}$ [36] 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\|_{\mathcal{L}^{2}(\mathcal{M})}\leq M_{2}$ a.s.

######  Assumption 3.2.

There exist a smooth function $\beta^{0}\in\mathcal{W}^{2}(\mathcal{M})$ such
that $\beta^{0}=L_{C}^{-1}(\mu_{2}-\mu_{1})$.

######  Assumption 3.3.

The effective dimension of $T$ satisfies
$D(\lambda)=\operatorname{Tr}((T+\lambda I)^{-1}T)\leq c\lambda^{-\theta}$ for
constants $c,\theta>0$. Here $\operatorname{Tr}$ denotes the trace operator.

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
$\beta^{0}$ to be estimated is well-defined and belongs to the space of smooth
functions $\mathcal{W}^{2}(\mathcal{M})$. Assumption 3.3 is expressed in terms
of properties of the effective dimension $D(\lambda)$. For our final choice of
$\lambda$, it is straightforward to check that this assumption holds by
assuming that the eigenvalues $\\{\tau_{k}\\}$ of $T$ decay as $\tau_{k}\asymp
k^{-2r}$, with $r>\frac{1}{2}$, which is a typical assumption in the
functional linear model literature [36] and is related to the rate of decay of
the eigenvalues of $L_{K_{\mathcal{M}}}$ and $L_{C}$, and their alignment.

The following theorem provides an upper bound for the excess prediction risk.

######  Theorem 3.1.

Under the Assumptions 3.1-3.3, if $\lambda\asymp n^{-\frac{1}{1+\theta}}$, the
estimate $\hat{\beta}$ in (8) is such that

$\displaystyle\mathbb{E}^{*}\left[\langle
X^{*},\beta^{0}-\hat{\beta}\rangle_{\mathcal{L}^{2}(\mathcal{M})}\right]^{2}=\mathcal{O}_{p}\left(n^{-\frac{1}{1+\theta}}\right).$
(11)

Similar rates of convergence have been shown to hold for regularized estimates
in the functional linear regression setting [[, see e.g.
]]Cai2012a,Tong2018,Sun2018,Reimherr2018. However, a key difference in our
model is that the residual random variable $\varepsilon=Y-\langle
X,\beta^{0}\rangle_{\mathcal{L}^{2}(\mathcal{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.

##  4 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 $\\{(g_{i},v_{i},x_{i})\\}$, where
$v_{i}\in\mathcal{V}(\mathbb{R}^{3})$ is a vector field representing the
subject-specific geometry of the $i$th subject. Recall that
$\left(\mathcal{V}(\mathbb{R}^{3}),\|\cdot\|_{\mathcal{V}(\mathbb{R}^{3})}\right)$
is an RKHS of smooth functions with compact support in $\mathbb{R}^{3}$.
Moreover, denote with $K_{\mathcal{V}}(\mathbb{R}^{3})$ its reproducing
kernel.

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 $\mathcal{V}(\mathbb{R}^{3})$, $X$
a zero-mean random function taking values in $\mathcal{L}^{2}(\mathcal{M})$,
and $G$ a binary random variable such that
$P\left(G=g_{1}\right)=\pi_{1},\left(G=g_{2}\right)=\pi_{2}$.

We now adopt the multivariate functional data notation in [59], and define
$\bm{Y}(p)=\left(V(p_{1}),X(p_{2})\right)$ with $p=(p_{1},p_{2})\in
D=D_{1}\times D_{2}=\mathbb{R}^{3}\times\mathcal{M}$. The multivariate random
function $\bm{Y}$ takes values in a Hilbert space
$\mathcal{H}=\mathcal{L}^{2}(\mathbb{R}^{3},\mathbb{R}^{3})\times\mathcal{L}^{2}(\mathcal{M})$
with inner product $\langle\bm{f},\bm{g}\rangle_{\mathcal{H}}=\langle
f^{(1)},g^{(1)}\rangle_{\mathcal{L}^{2}(\mathbb{R}^{3},\mathbb{R}^{3})}+\langle
f^{(2)},g^{(2)}\rangle_{\mathcal{L}^{2}(\mathcal{M})}$ for
$f,g\in\mathcal{H}$. Here $f^{(j)}$, with $j\in\\{1,2\\}$, denotes the $j$th
functional component of the multivariate function $\bm{f}$. For $p,q\in D$,
define the matrix of covariances
$\bm{C}(p,q)=\mathbb{E}(\bm{Y}(p)\otimes\bm{Y}(q))$ with elements
$C_{lj}(p_{l},q_{j})=\mathbb{E}[Y^{(l)}(p_{l})Y^{(j)}(q_{j})]$ with $p_{l}\in
D_{l},q_{j}\in D_{j}$.

Denote the conditional means of $\bm{Y}$ by
$\bm{\mu}_{1}=\left(\mu_{1}^{(1)},\mu_{1}^{(2)}\right):=\left(\mathbb{E}\left[V|G=g_{1}\right],\mathbb{E}\left[X|G=g_{1}\right]\right)$
and
$\bm{\mu}_{2}=\left(\mu_{2}^{(1)},\mu_{2}^{(2)}\right):=\left(\mathbb{E}\left[V|G=g_{2}\right],\mathbb{E}\left[X|G=g_{2}\right]\right)$,
and assume $\bm{\mu}_{1}\neq\bm{\mu}_{2}$. Next, we define the covariance
operator $L_{\bm{C}}:\mathcal{H}\rightarrow\mathcal{H}$, with the $j$th
element of $L_{\bm{C}}\bm{f},\bm{f}\in\mathcal{H}$ given by

$(L_{\bm{C}}\bm{f})^{(j)}(p_{j})=\sum_{i=1}^{2}\int_{D_{i}}C_{ij}(q_{i},p_{j})f^{(i)}(q_{i}).$
(12)

As in the univariate case, we then assume that the population quantity
$\bm{\beta}^{0}\in\mathcal{H}$ is well defined and satisfies the equation

$L_{\bm{C}}\bm{\beta}^{0}=\bm{\mu}_{2}-\bm{\mu}_{1}.$

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 $\bm{\beta}^{0}$.

###  4.1 Regularized estimation

Let the variable $y_{i}$ be such that $y_{i}=-\frac{n}{n_{1}}$ if $g_{i}=1$
and $y_{i}=\frac{n}{n_{2}}$ otherwise. We define the multivariate functional
estimate $\bm{\hat{\beta}}=\left(\hat{\beta}^{G},\hat{\beta}^{F}\right)$ of
the population quantity $\bm{\beta^{0}}$ to be the solution to the following
minimization problem

$\left(\hat{\beta}^{G},\hat{\beta}^{F}\right)=\operatorname*{arg\,min}_{\begin{subarray}{c}\beta^{G}\in\mathcal{V}(R^{3})\\\
\beta^{F}\in\mathcal{W}^{2}(\mathcal{M})\end{subarray}}\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}-\langle
v_{i},\beta^{G}\rangle_{\mathcal{L}^{2}(\mathbb{R}^{3},\mathbb{R}^{3})}-\langle
x_{i},\beta^{F}\rangle_{\mathcal{L}^{2}(\mathcal{M})}\right)^{2}+\lambda_{1}\|\beta^{G}\|_{\mathcal{V}(\mathbb{R}^{3})}^{2}+\lambda_{2}J(\beta^{F}).$
(13)

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 $\hat{\beta}^{G}$ and
$\hat{\beta}^{F}$ in their respective function spaces.

###  4.2 Differential regularization and kernel penalty: a unified modeling
perspective

In Section 4.1, we have adopted two different approaches to produce smooth
estimates $\hat{\beta}^{G}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ and
$\hat{\beta}^{F}:\mathcal{M}\rightarrow\mathbb{R}$. The smoothness of
$\hat{\beta}^{F}$ is enforced by means of a penalty $J(\cdot)$ defined in
terms of the $\mathcal{L}^{2}(\mathcal{M})$ norm of a differential operator,
which implicitly defines a kernel $K_{\mathcal{M}}$. On the other hand, the
smoothness of $\hat{\beta}^{G}$ is enforced by means of a norm
$\|\cdot\|_{\mathcal{V}(\mathbb{R}^{3})}$, defined implicitly through the
direct definition of a kernel $K_{\mathbb{R}^{3}}$. For clarity, we summarize
the relevant function spaces and associated norms and kernels in Table 1.

Estimate | Function space | Norm | Kernel
---|---|---|---
$\hat{\beta}^{F}:\mathcal{M}\rightarrow\mathbb{R}$ | $\mathcal{W}^{2}(\mathcal{M})$ | $J^{1/2}(\cdot)=\left(\|\Delta_{\mathcal{M}}\cdot\|^{2}_{\mathcal{L}^{2}(\mathcal{M})}+\varepsilon\|\cdot\|^{2}_{\mathcal{L}^{2}(\mathcal{M})}\right)^{1/2}$ |  $K_{\mathcal{M}}$ (implicit)
$\hat{\beta}^{G}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ | $\mathcal{V}(\mathbb{R}^{3})$ |  $\|\cdot\|_{\mathcal{V}(\mathbb{R}^{3})}$ (implicit) | $K_{\mathbb{R}^{3}}$
Table 1: Table summarizing estimates and associated function spaces, norms and
kernels.

From a methodological perspective, the reproducing kernel
$K_{\mathbb{R}^{3}}(p,q),p,q\in\mathbb{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 [6, 30], 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
$\hat{\beta}^{F}:\mathcal{M}\rightarrow\mathbb{R}$. In other words, can we
define a smooth real function space on $\mathcal{M}$ by explicitly defining a
kernel $K_{\mathcal{M}}:\mathcal{M}\times\mathcal{M}\rightarrow\mathbb{R}$
that encodes a measure of influence that is coherent with the non-linear
geometry of $\mathcal{M}$? This turns out to be a non-trivial task due to the
positive definiteness property that $K_{\mathcal{M}}$ needs to satisfy.
Consider, for instance, the popular exponential kernel. Its natural extension
to the manifold setting is
$K_{\mathcal{M}}(p,q)=\exp(-c\,\text{d}_{\mathcal{M}}(p,q)^{2})$, where
$\text{d}_{\mathcal{M}}(p,q)$ is the geodesic distance between
$p\in\mathcal{M}$ and $q\in\mathcal{M}$. Unfortunately, this kernel is not
guaranteed to be positive definite for a generic non-linear manifold
$\mathcal{M}$ [44, 45].

Alternatively, we could try to compute an explicit form of the kernel
$K_{\mathcal{M}}$ from $J(\cdot)$ by employing the following identity

$f(p)=\langle K(p,\cdot),f\rangle_{\mathcal{W}^{2}(\mathcal{M})},\forall
f\in\mathcal{W}^{2}(\mathcal{M}),$

that is, for any fixed $p\in\mathcal{M}$,
$K_{\mathcal{M}}(p,\cdot)\in\mathcal{W}^{2}(\mathcal{M})$ is given by the
solution to the partial differential equation

$(\Delta\Delta+\varepsilon\text{I})K_{\mathcal{M}}(p,\cdot)=\delta_{p},\qquad
p\in\mathcal{M},$ (14)

where here $\delta_{p}$ denotes the Dirach delta [6, 37]. However, closed form
solutions of (14) are not available in our setting. Approximate solutions
could be computed by Finite Elements Analysis (FEA) [28]. Yet, we would face
the issue of having to store the dense object $K_{\mathcal{M}}(\cdot,\cdot)$.
As described in Section 5, we will instead leverage FEA to directly discretize
the function $\hat{\beta}^{F}:\mathcal{M}\rightarrow\mathbb{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.

###  4.3 Theory

Define the diagonal matrix of reproducing kernels $\bm{K}(p,q)$ with entries
$\bm{K}_{11}(p_{1},q_{1})=K_{\mathbb{R}^{3}}(p_{1},q_{1})$ and
$\bm{K}_{22}(p_{2},q_{2})=K_{\mathcal{M}}(p_{2},q_{2})$ with $p_{i}\in
D_{i},p_{j}\in D_{j}$. Let $L_{\bm{K}}:\mathcal{H}\rightarrow\mathcal{H}$ be
the associated integral operator and, analogously to the univariate setting,
define the sandwich operator
$T=L_{\bm{K}^{\frac{1}{2}}\bm{C}\bm{K}^{\frac{1}{2}}}$. We make the following
assumptions.

######  Assumption 4.1.

The norm of $\bm{Y}$ is a.s. bounded by a positive real number, i.e., there
exists $M>0$ such that $\|\bm{Y}\|_{\mathcal{H}}\leq M$ (a.s.).

######  Assumption 4.2.

There exist a smooth function
$\bm{\beta}^{0}\in\mathcal{V}(\mathbb{R}^{3})\times\mathcal{W}^{2}(\mathcal{M})$
such that $\bm{\beta}^{0}=L_{\bm{C}}^{-1}(\bm{\mu}_{2}-\bm{\mu}_{1})$.

######  Assumption 4.3.

The penalty coefficients $\lambda_{1}=\lambda_{2}=:\lambda$ and the effective
dimension of $T$ satisfies $D(\lambda)=\operatorname{Tr}((T+\lambda
I)^{-1}T)\leq c\lambda^{-\theta}$ for constants $c,\theta>0$.

The following theorem, which is an extension of Theorem 3.1, provides an upper
bound for the excess prediction risk.

######  Theorem 4.1.

Under Assumptions 4.1-4.3,, if $\lambda_{1}=\lambda_{2}\asymp
n^{-\frac{1}{1+\theta}}$, the estimate
$\bm{\hat{\beta}}=\left(\hat{\beta}^{G},\hat{\beta}^{F}\right)$ in (13) is
such that

$\displaystyle\mathbb{E}^{*}\left[\langle\bm{Y}^{*},\bm{\beta}^{0}-\bm{\hat{\beta}}\rangle_{\mathcal{H}}\right]^{2}=\mathcal{O}_{p}\left(n^{-\frac{1}{1+\theta}}\right),$
(15)

where $\bm{Y}^{*}$ is an independent copy of $\bm{Y}$ that is independent of
the training data and $\mathbb{E}^{*}$ is the expectation taken over
$\bm{Y}^{*}$.

##  5 Discretization

Consider now a triangular mesh $\mathcal{M}_{\mathcal{T}}$ that is the union
of the finite set of triangles $\mathcal{T}$ defined over a set of $s$ common
vertices $\xi_{1},\ldots,\xi_{s}$. Let this be an approximate representation
of the manifold $\mathcal{M}$. We then introduce the linear finite element
space $\mathcal{W}_{\mathcal{T}}$ consisting of a set of globally continuous
functions over $\mathcal{M}_{\mathcal{T}}$ that are affine within each
triangle $\tau$ in $\mathcal{T}$, i.e.

$\mathcal{W}_{\mathcal{T}}=\\{w\in
C^{0}(\mathcal{M}_{\mathcal{T}}):w|_{\tau}\text{ is affine
}\forall\tau\in\mathcal{T}\\}.$

The space $\mathcal{W}_{\mathcal{T}}$ is spanned by the Finite Elements (FE)
basis $\psi_{1},\ldots,\psi_{s}$, where $\psi_{l}(\xi_{j})=1$, if $l=j$, and
$\psi_{l}(\xi_{j})=0$ otherwise. In Figure 2, we show one element of such a
basis. Moreover, define $\bm{\psi}$ as the vector-valued function
$\bm{\psi}(\cdot)=\left(\psi_{1}(\cdot),\ldots,\psi_{s}(\cdot)\right)^{T}$. We
search for an approximate solution $\hat{\beta}^{F}_{\mathcal{T}}$ of the form

$\beta_{\mathcal{T}}^{F}(\cdot)=\sum_{l=1}^{s}c^{F}_{l}\psi_{l}(\cdot)=(\bm{c}^{F})^{T}\bm{\psi}(\cdot),$
(16)

where $\bm{c}^{F}=\left(c^{F}_{1},\cdots,c^{F}_{s}\right)^{T}$.

Let now $M$ and $S$ be the sparse mass and stiffness $s\times s$ matrices
defined as
$(M)_{jj^{\prime}}=\int_{\mathcal{M}_{\mathcal{T}}}\psi_{j}\psi_{j^{\prime}}$
and
$(S)_{jj^{\prime}}=\int_{\mathcal{M}_{\mathcal{T}}}\nabla_{\mathcal{M}_{\mathcal{T}}}\psi_{j}\cdot\nabla_{\mathcal{M}_{\mathcal{T}}}\psi_{j^{\prime}}$,
where $\nabla_{\mathcal{M}_{\mathcal{T}}}$ is the gradient operator on the
manifold $\mathcal{M}$. For $\beta_{\mathcal{T}}^{F}$ of the form given in
equation (16), we have that the associated penalty term $J(\cdot)$ can be
approximated as $(\bm{c}^{F})^{T}D_{\mathcal{M}_{\mathcal{T}}}\bm{c}^{F}$,
with $D_{\mathcal{M}_{\mathcal{T}}}=SM^{-1}S+\varepsilon M$ [49]. Further,
following an approach often adopted in FEA [1, 2], we replace the dense matrix
$M^{-1}$ with the sparse matrix $\tilde{M}^{-1}$, where $\tilde{M}$ is the
diagonal matrix such that $\tilde{M}_{jj}=\sum_{j^{\prime}}M_{jj^{\prime}}$.
This leads to the sparse penalty matrix
$D_{\mathcal{M}_{\mathcal{T}}}=S\tilde{M}^{-1}S+\varepsilon M$. In practice,
each functional observation $x_{i}$ is also of the form given in equation
(16). Therefore, denoting by $\mathbb{X}$ the $n\times s$ matrix, where each
row contains the basis coefficients of $x_{i}$, the term $\langle
x_{i},\beta^{F}\rangle_{\mathcal{L}^{2}(\mathcal{M})}$ can be discretized by
$\mathbb{X}M\bm{c}^{F}$.

We now turn our attention to the estimate
$\hat{\beta}^{G}\in\mathcal{V}(\mathbb{R}^{3})$. Since for
$\mathcal{V}(\mathbb{R}^{3})$ we have an explicit form of the associated
reproducing kernel, we employ the representer theorem [6, 30] and take
$\hat{\beta}^{G}$ of the form

$\hat{\beta}^{G}(\cdot)=\sum_{i=1}^{n}c^{G}_{i}\int_{\mathbb{R}^{3}}K_{\mathbb{R}^{3}}(p,\cdot)v_{i}(p)dp.$
(17)

On the right hand side of Figure 2, we show an example of a function
$\int_{\mathbb{R}^{3}}K_{\mathbb{R}^{3}}(p,\cdot)v_{i}(p)dp$. Then, we have
that $\|\hat{\beta}^{G}\|^{2}_{\mathcal{V}}=(\bm{c}^{G})^{T}\Sigma\bm{c}^{G}$,
where $\bm{c}^{G}=\left(c_{1}^{G},\ldots,c_{n}^{G}\right)^{T}$ and $\Sigma$ is
a $n\times n$ matrix with entries

$\Sigma_{ij}=\int\int v_{i}(p)^{T}K_{\mathbb{R}^{3}}(p,q)v_{j}(q)\,dp\,dq.$

Figure 2: On the left hand side, an element of the FE basis
$\left\\{\psi_{l}:\mathcal{M}_{\mathcal{T}}\rightarrow\mathbb{R},l=1,\ldots,s\right\\}$.
This is a scalar linear function within each triangle of the mesh
$\mathcal{M}_{\mathcal{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
$\left\\{\int_{\mathbb{R}^{3}}K_{\mathbb{R}^{3}}(p,\cdot)v_{i}(p)dp,i=1,\ldots,n\right\\}$.
This is a smooth vector-valued function from $\mathbb{R}^{3}$ to
$\mathbb{R}^{3}$.

As a result, the coefficients of the approximate solution of model (13) are
given by

$\left(\bm{\hat{c}}^{G},\bm{\hat{c}}^{F}\right)=\operatorname*{arg\,min}_{\bm{c}^{G}\in\mathbb{R}^{n},\bm{c}^{F}\in\mathbb{R}^{s}}\|\bm{y}-\Sigma\bm{c}^{G}-\mathbb{X}M\bm{c}^{F}\|^{2}_{2}+\lambda_{1}(\bm{c}^{G})^{T}\Sigma\bm{c}^{G}+\lambda_{2}(\bm{c}^{F})^{T}D_{\mathcal{M}_{\mathcal{T}}}\bm{c}^{F},$
(18)

where $\bm{y}=\left(y_{1},\ldots,y_{n}\right)^{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

$\left(\bm{\hat{c}}^{G},\bm{\hat{c}}^{F}\right)=\operatorname*{arg\,min}_{\bm{c}^{G}\in\mathbb{R}^{n},\bm{c}^{F}\in\mathbb{R}^{s}}\left\|\begin{bmatrix}\bm{y}\\\
\bm{0}\end{bmatrix}-A\begin{bmatrix}\bm{c}^{G}\\\
\bm{c}^{F}\end{bmatrix}\right\|^{2}_{2},$ (19)

with

$A=\begin{bmatrix}\Sigma&\mathbb{X}M\\\
0&\lambda_{2}^{\frac{1}{2}}S\tilde{M}^{-\frac{1}{2}}\\\
0&\lambda_{2}^{\frac{1}{2}}\varepsilon^{\frac{1}{2}}M^{\frac{1}{2}}\\\
\lambda_{1}^{\frac{1}{2}}\Sigma^{\frac{1}{2}}&0\end{bmatrix}.$

Note that for $n\ll 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 [3], 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.

##  6 Applications

###  6.1 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 [48]), 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 [10, 9].
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 [11], fMRI signals, or
connectivity maps [39, 33]. 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
$\mathcal{M}^{\mathcal{T}}_{i}$ consisting of $s=64$K 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 [47], translation, rigid rotation, and scale were removed
from the data, while jointly estimating a template
$\mathcal{M}^{\mathcal{T}}$, which is also a triangular mesh with $s=64$K
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^{\mathcal{T}}_{i}$ over the mesh $\mathcal{M}^{\mathcal{T}}_{i}$.

The preprocessing stage results in a set FoSs
$\left\\{\left(\mathcal{M}^{\mathcal{T}}_{i},z^{\mathcal{T}}_{i}\right),i=1,\ldots,n\right\\}$,
which are discretized versions of the continuous idealized objects
$\left\\{\left(\mathcal{M}_{i},z_{i}\right),i=1,\ldots,n\right\\}$ introduced
in Section 2. To simplify the notation, we drop the superscript $\mathcal{T}$.
Moreover, we denote the diagnostic labels by
$\left\\{g_{i}\in\\{\text{C},\text{AD}\\},i=1,\ldots,n\right\\}$. 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 data-driven hypotheses about the pathophysiology
of AD.

###  6.2 Linear functional representations

For each FoS, we compute a function $v_{i}\in\mathcal{V}(\mathbb{R}^{3})$ such
that $\varphi_{v_{i}}(\mathcal{M})$ is a close approximation of
$\mathcal{M}_{i}$, where closeness is measured as the Euclidean distance
between the corresponding vertices of
$\varphi_{v_{i}}\left(\mathcal{M}\right)$ and $\mathcal{M}_{i}$. We can then
transport the function $z_{i}:\mathcal{M}_{i}\rightarrow\mathbb{R}$ onto the
template defining a continuous piecewise linear function
$x_{i}=z_{i}\circ\varphi_{v_{i}}$. This leads to the definition of the
bivariate functional data $\left(v_{i},x_{i}\right)$ that are linear
representations of the FoSs
$(\mathcal{M}_{i},z_{i})\approx(\varphi_{v_{i}}(\mathcal{M}),x_{i}\circ\varphi^{-1}_{v_{i}})$.
Further details on the computation of $v_{i}$ can be found in the appendices.

###  6.3 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}-\bar{v},x_{i}-\bar{x}\\}$, where $\bar{v}$ and $\bar{x}$ 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 estimates
$\hat{\beta}^{G}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ and
$\hat{\beta}^{F}:\mathcal{M}\rightarrow\mathbb{R}$. Given a new subject with
predictors $(v_{*},x_{*})$, these estimates can be used in conjunction with
the classification rule $\langle v_{*}-\bar{v},\hat{\beta}^{G}\rangle+\langle
x_{*}-\bar{x},\hat{\beta}^{F}\rangle>c^{\text{th}}$ to predict the diagnostic
label of a new subject. The cut-off level $c^{\text{th}}$ can be chosen by
computing sensitivity and specificity on a test set, for different values
$c^{\text{th}}$, and selecting the desired level and type of accuracy.

Figure 3: On the left hand side, we show the most discriminant geometric and
thickness ‘directions’ as estimated from the linear representations
$\\{(v_{i}-\bar{v},x_{i}-\bar{x})\\}$. These are a vector field
$\hat{\beta}^{G}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$, representing the
most predictive geometric pattern of AD, and a function
$\hat{\beta}^{F}:\mathcal{M}\rightarrow\mathbb{R}$, representing the most
predictive cortical thickness pattern of AD. For a new FoS, with linear
representation $(v_{*},x_{*})$, we compute the score $\langle
v_{*}-\bar{v},\hat{\beta}^{G}\rangle+\langle
x_{*}-\bar{x},\hat{\beta}^{F}\rangle$ and predict whether the subject has AD
by comparing the score value with a predetermined threshold $c^{\text{th}}$.
On the right hand side, we depict the process of mapping back the estimates
$\hat{\beta}^{G}$ and $\hat{\beta}^{F}$ onto the space of FoSs. On the same
space, we also pictorially map the classification rule adopted. The blue areas
in the $\hat{\beta}^{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
$\hat{\beta}^{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
$\varphi_{\bar{v}-c_{1}\hat{\beta}^{G}}(\mathcal{M})$ and
$\varphi_{\bar{v}+c_{1}\hat{\beta}^{G}}(\mathcal{M})$, on the right hand side
diagram.

These estimates effectively identify linear directions
$\\{c_{1}\hat{\beta}^{G},c_{1}\in\mathbb{R}\\}$ and
$\\{c_{2}\hat{\beta}^{F},c_{2}\in\mathbb{R}\\}$, in their respective spaces,
that can be interpreted as the most discriminant geometric and thickness
directions. Specifically, large values of $c_{1}\in\mathbb{R}$ and
$c_{2}\in\mathbb{R}$ describe configurations $c_{1}\hat{\beta}^{G}$ and
$c_{2}\hat{\beta}^{F}$ that are predictive of AD. Low values of
$c_{1}\in\mathbb{R}$ and $c_{2}\in\mathbb{R}$ describe configurations that are
instead predictive of the subject being healthy. Moreover, given the additive
modeling assumption on the geometric and thickness components, for every
configuration $\left(c_{1}\hat{\beta}^{G},c_{2}\hat{\beta}^{F}\right)$, 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

$\left(\varphi_{\bar{v}+c_{1}\hat{\beta}^{G}}(\mathcal{M}),c_{2}\hat{\beta}^{F}\circ\varphi^{-1}_{\bar{v}+c_{1}\hat{\beta}^{G}}\right),\qquad
c_{1},c_{2}\in\mathbb{R}.$

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
[19].

We fit the proposed model for different choices of the parameters
$\lambda_{1}$ and $\lambda_{2}$. Recall that $\lambda_{1}$ controls the
regularity of the geometric estimate and $\lambda_{2}$ that of the thickness,
with high values virtually constraining the solution to be the zero function.
The final choice of $\lambda_{1}$ and $\lambda_{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
$\lambda_{2}$ and with $\lambda_{1}=+\infty$, i.e., discarding the geometric
component.

###  6.4 Results

On the left hand side of Figure 3, we show the estimated most discriminant
geometric and thickness directions, i.e.,
$\hat{\beta}^{G}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ and
$\hat{\beta}^{F}:\mathcal{M}\rightarrow\mathbb{R}$. These have been estimated
by applying the model in (13) to the linear representations
$\\{(v_{i}-\bar{v},x_{i}-\bar{x})\\}$. The colormap describing
$\hat{\beta}^{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.,
]]Bondareff1989,Dickerson2009,Sabuncu2011. The geometric component
$\hat{\beta}^{G}$ is instead a vector field in $\mathbb{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 $\hat{\beta}^{G}$ back to the space of FoSs, i.e.,
exploring $\varphi_{\bar{v}+c_{1}\hat{\beta}^{G}}(\mathcal{M})$ for different
choices of $c_{1}\in\mathbb{R}$, the magnitude of the vector field
$\hat{\beta}^{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 $\hat{\beta}^{G}$ and $\hat{\beta}^{F}$, that is,
$\left(\varphi_{\bar{v}+c_{1}\hat{\beta}^{G}}(\mathcal{M}),c_{2}\hat{\beta}^{F}\circ\varphi^{-1}_{\bar{v}+c_{1}\hat{\beta}^{G}}\right)$
with $c_{1},c_{2}\in\mathbb{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.

###  6.5 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.,
]]Bondareff1989,Dickerson2009,Sabuncu2011. However, there is some variability
in the estimated regions. For instance, [32] 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 [71]. 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.

##  7 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 [58] or cardiac
imaging [55]. 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;
Alzheimer’s Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.;
Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Cogstate; Eisai Inc.;
Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La
Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE
Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research &
Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development
LLC.; Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale Diagnostics, LLC.;
NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals
Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical
Company; and Transition Therapeutics. The Canadian Institutes of Health
Research is providing funds to support ADNI clinical sites in Canada. Private
sector contributions are facilitated by the Foundation for the National
Institutes of Health (www.fnih.org). The grantee organization is the Northern
California Institute for Research and Education, and the study is coordinated
by the Alzheimer’s Therapeutic Research Institute at the University of
Southern California. ADNI data are disseminated by the Laboratory for Neuro
Imaging at the University of Southern California.

Part of the data used in the preparation of this article was obtained from the
Parkinson’s Progression Markers Initiative (PPMI) database (www.ppmi-
info.org/data). For up-to-date information on the study, visit www.ppmi-
info.org. PPMI — a public-private partnership — is funded by the Michael J.
Fox Foundation for Parkinson’s Research and funding partners, including [list
the full names of all of the PPMI funding partners found at www.ppmi-
info.org/fundingpartners].

##  Appendix A Linear functional representation model

###  A.1 Diffeomorphic deformation operator

The diffeomorphic operator $\varphi$ can be constructed as follows. Let
$\\{v_{t}\in\mathcal{V}(\mathbb{R}^{3}):t\in[0,1]\\}$ be a time-variant vector
field such that
$\int_{0}^{1}\|v_{t}\|^{2}_{\mathcal{V}(\mathbb{R}^{3})}\,dt<\infty$. Then,
the solution $\phi_{v}:[0,1]\times\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$, at
time $t=1$, to the Ordinary Differential Equation (ODE)

$\displaystyle\begin{cases}\frac{\partial\phi_{v}}{\partial
t}(t,x)=v_{t}\circ\phi_{v}(t,x)\ \qquad&t\in[0,1],x\in\mathbb{R}^{3},\\\
\phi_{v}(0,x)=x\qquad&x\in\mathbb{R}^{3},\\\ \end{cases}$ (20)

is a diffeomorphic deformation of $\mathbb{R}^{3}$ [[, see, e.g.,]]Younes2010.
We then model $\\{v_{t}:t\in[0,1]\\}$ to be a minimizer of the quantity
$\int_{0}^{1}\|v_{t}\|^{2}_{\mathcal{V}}dt$, for a given initial vector field
$v_{0}\in\mathcal{V}(\mathbb{R}^{3})$ [21]. Finally, the diffeomorphic
operator is defined as $\varphi_{v_{0}}(x)=\phi_{v}(1,x)$, where
$v_{0}\in\mathcal{V}(\mathbb{R}^{3})$ is the initial vector field generating
$\\{v_{t}:t\in[0,1]\\}$, and $\phi_{v}$ is the solution of the ODE (20) for
the computed $\\{v_{t}:t\in[0,1]\\}$.

###  A.2 Computation

In practice, each surface $\mathcal{M}_{i}$ has a computational representation
$\mathcal{M}^{\mathcal{T}}_{i}$ that is a triangular mesh with $s$ vertices
$\xi^{i}_{1},\ldots,\xi^{i}_{s}$ in correspondence across the $n$ subjects. We
use these vertices to perform Procrustes analysis [47] and remove translation,
size, and rigid rotations from the surfaces $\mathcal{M}^{\mathcal{T}}_{i}$.
In addition, Procrustes analysis yields a template $\mathcal{M}^{\mathcal{T}}$
with vertices $\xi_{1},\ldots,\xi_{s}$.

The linear representation functions
$\\{v_{i}:v_{i}\in\mathcal{V}(\mathbb{R}^{3})\\}$, associated with the
surfaces $\mathcal{M}^{\mathcal{T}}_{i}$, are then computed by solving the
minimization problem

${v}_{i}=\operatorname*{arg\,min}_{v\in\mathcal{V}}\sum_{l=1}^{s}\|\varphi_{v}(\xi_{l})-\xi^{i}_{l}\|^{2}_{\mathbb{R}^{3}}+\lambda\|v\|^{2}_{\mathcal{V}},\qquad
i=1,\ldots,n,$ (21)

where the least-squares term ensures that the deformed template
$\varphi_{v_{i}}(\mathcal{M}^{\mathcal{T}})$ is a close approximation of
$\mathcal{M}^{\mathcal{T}}_{i}$. The term
$\|v\|^{2}_{\mathcal{V}(\mathbb{R}^{3})}$ is a regularizing term that
encourages the solution to achieve such a close approximation with a ‘minimal’
deformation. The constant $\lambda$ 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 [43, 50]. 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 [18, 24].

##  Appendix B Proofs

Recall that we denote by $C$ the covariance function of the functional
predictors and by $\hat{C}$ its empirical counterpart. Moreover,
$K_{\mathcal{M}}$ denotes the evaluation kernel of the Sobolev space
$\mathcal{W}^{2}(\mathcal{M})$ endowed with the norm
$\left(\|\Delta_{\mathcal{M}}\beta\|^{2}_{\mathcal{L}^{2}(\mathcal{M})}+\varepsilon\|\beta\|^{2}_{\mathcal{L}^{2}(\mathcal{M})}\right)^{1/2}$.
The sandwich operator $T$ is defined as
$T=L_{K_{\mathcal{M}}}^{\frac{1}{2}}L_{C}L_{K_{\mathcal{M}}}^{\frac{1}{2}}$.
To simplify the notation, we drop the subscript from
$\|\cdot\|_{\mathcal{L}^{2}(\mathcal{M})}$ and denote it by $\|\cdot\|$. We
define $\|A\|_{\operatorname{op}}=\sup_{f:\|f\|=1}\|Af\|$ to be the operator
norm of a linear operator
$A:\mathcal{L}^{2}(\mathcal{M})\rightarrow\mathcal{L}^{2}(\mathcal{M})$.

Thanks to the fact that
$L_{K_{\mathcal{M}}}^{\frac{1}{2}}(\mathcal{L}^{2}(\mathcal{M}))=\mathcal{W}^{2}(\mathcal{M})$
[14], it is clear that we can reformulate problem (8) as

$\text{minimize}_{f}\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}-\langle
x_{i},L_{K_{\mathcal{M}}}^{\frac{1}{2}}f\rangle_{\mathcal{L}^{2}(\mathcal{M})}\right)^{2}+\lambda\|f\|^{2}_{\mathcal{L}^{2}(\mathcal{M})},$
(22)

whose solution $\hat{f}_{\lambda}$ is given by

$\hat{f}_{\lambda}=(T_{n}+\lambda
I)^{-1}\frac{1}{n}\sum_{i=1}^{n}y_{i}L_{K_{\mathcal{M}}}^{\frac{1}{2}}x_{i},$
(23)

and where
$T_{n}=L_{K_{\mathcal{M}}}^{\frac{1}{2}}L_{\hat{C}}L_{K_{\mathcal{M}}}^{\frac{1}{2}}$.
Moreover, the out-of-sample risk $\mathbb{E}^{*}\left[\langle
X^{*},\beta^{0}-\hat{\beta}\rangle_{\mathcal{L}^{2}(\mathcal{M})}\right]^{2}$
can be easily rewritten as
$\left\|T^{\frac{1}{2}}(\hat{f}_{\lambda}-f_{0})\right\|^{2}$, with
$f_{0}\in\mathcal{L}^{2}(\mathcal{M})$ such that
$f_{0}=T^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}(\mu_{2}-\mu_{1})$

Observe that

$\displaystyle T^{\frac{1}{2}}(\hat{f}_{\lambda}-f_{0})$ (24)
$\displaystyle=T^{\frac{1}{2}}\left[(T_{n}+\lambda
I)^{-1}\left(\frac{1}{n}\sum_{i=1}^{n}y_{i}L_{K_{\mathcal{M}}}^{\frac{1}{2}}x_{i}\right)-f_{0}\right]$
(25) $\displaystyle=T^{\frac{1}{2}}\left[(T_{n}+\lambda
I)^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}\left(\frac{1}{n}\sum_{i=1}^{n}y_{i}x_{i}-(\mu_{2}-\mu_{1})\right)+(T_{n}+\lambda
I)^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}(\mu_{2}-\mu_{1})-T^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}(\mu_{2}-\mu_{1})\right]$
(26)

Let moreover $\hat{d}=\frac{1}{n}\sum_{i=1}^{n}y_{i}x_{i}$ and
$d=\mu_{2}-\mu_{1}$, and notice that $\mathbb{E}[\hat{d}]=d$.

Then, we have

$\displaystyle\left\|T^{\frac{1}{2}}(\hat{f}_{\lambda}-f_{0})\right\|$
$\displaystyle\leq\left\|T^{\frac{1}{2}}(T_{n}+\lambda
I)^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}\left(\hat{d}-d\right)\right\|+\left\|T^{\frac{1}{2}}\left[(T_{n}+\lambda
I)^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}d-T^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}d\right]\right\|$
(27) $\displaystyle=I_{1}+I_{2}$

We first derive a bound for the term $I_{1}$ in (B), 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

$\displaystyle I_{1}$ $\displaystyle=\left\|T^{\frac{1}{2}}(T_{n}+\lambda
I)^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}\left(\hat{d}-d\right)\right\|$
$\displaystyle=\left\|T^{\frac{1}{2}}(T_{n}+\lambda
I)^{-\frac{1}{2}}(T_{n}+\lambda I)^{-\frac{1}{2}}(T+\lambda
I)^{\frac{1}{2}}(T+\lambda
I)^{-\frac{1}{2}}L_{K_{\mathcal{M}}}^{\frac{1}{2}}\left(\hat{d}-d\right)\right\|$
$\displaystyle\leq\|T^{\frac{1}{2}}(T_{n}+\lambda
I)^{-\frac{1}{2}}\|_{\text{op}}\|(T_{n}+\lambda I)^{-\frac{1}{2}}(T+\lambda
I)^{\frac{1}{2}}\|_{\text{op}}\left\|(T+\lambda
I)^{-\frac{1}{2}}L^{\frac{1}{2}}_{K_{\mathcal{M}}}\left(\hat{d}-d\right)\right\|$
$\displaystyle\leq\|(T+\lambda I)^{\frac{1}{2}}(T_{n}+\lambda
I)^{-\frac{1}{2}}\|_{\text{op}}\|(T_{n}+\lambda I)^{-\frac{1}{2}}(T+\lambda
I)^{\frac{1}{2}}\|_{\text{op}}\left\|(T+\lambda
I)^{-\frac{1}{2}}L^{\frac{1}{2}}_{K_{\mathcal{M}}}\left(\hat{d}-d\right)\right\|$
$\displaystyle\leq\|(T+\lambda I)(T_{n}+\lambda
I)^{-1}\|_{\text{op}}\left\|(T+\lambda
I)^{-\frac{1}{2}}L^{\frac{1}{2}}_{K_{\mathcal{M}}}\left(\hat{d}-d\right)\right\|,$

where in the last inequality we have used
$\|A^{\gamma}B^{\gamma}\|_{\operatorname{op}}\leq\|AB\|^{\gamma}_{\operatorname{op}}$,
for any $0<\gamma<1$ [29]. We can then bound the first term thanks to Theorem
B.3. We therefore turn to the second term.

First observe that

$\displaystyle\mathbb{E}\left[\hat{d}\right]=d,\qquad\mathbb{E}\left[(\hat{d}-d)\otimes(\hat{d}-d)\right]=\frac{1}{n}C-\pi_{1}\mu_{1}\otimes\mu_{1}-\pi_{2}\mu_{2}\otimes\mu_{2}$

and, therefore,

$\displaystyle\langle
L_{K_{\mathcal{M}}}^{\frac{1}{2}}(\hat{d}-d),\eta_{k}\rangle^{2}$
$\displaystyle=\langle
L_{K_{\mathcal{M}}}^{\frac{1}{2}}L_{\mathbb{E}\left[(\hat{d}-d)\otimes(\hat{d}-d)\right]}L_{K_{\mathcal{M}}}^{\frac{1}{2}}\eta_{k},\eta_{k}\rangle$
$\displaystyle=\frac{1}{n}\langle
T\eta_{k},\eta_{k}\rangle-\frac{\pi_{1}}{n}\langle
L_{K_{\mathcal{M}}}^{\frac{1}{2}}L_{\mu_{1}\otimes\mu_{1}}L_{K_{\mathcal{M}}}^{\frac{1}{2}}\eta_{k},\eta_{k}\rangle-\frac{\pi_{2}}{n}\langle
L_{K_{\mathcal{M}}}^{\frac{1}{2}}L_{\mu_{2}\otimes\mu_{2}}L_{K_{\mathcal{M}}}^{\frac{1}{2}}\eta_{k},\eta_{k}\rangle$
$\displaystyle=\frac{1}{n}\left(\tau_{k}-\|L_{K_{\mathcal{M}}}^{\frac{1}{2}}\mu_{1}\|^{2}-\|L_{K_{\mathcal{M}}}^{\frac{1}{2}}\mu_{2}\|^{2}\right),$

with $\\{\tau_{k}\\}$ and $\\{\eta_{k}\\}$ denoting the eigenvalues and
eigenfunctions of $T$, respectively.

We then have

$\displaystyle\mathbb{E}\left\|(T+\lambda
I)^{-\frac{1}{2}}L^{\frac{1}{2}}_{K_{\mathcal{M}}}\left(\hat{d}-d\right)\right\|^{2}$
$\displaystyle=\mathbb{E}\left\|\sum_{k}(T+\lambda
I)^{-\frac{1}{2}}\eta_{k}\langle
L_{K_{\mathcal{M}}}^{\frac{1}{2}}(\hat{d}-d),\eta_{k}\rangle\right\|^{2}$
$\displaystyle\leq\frac{1}{n}\sum_{k}\left(\frac{\tau_{k}}{\lambda+\tau_{k}}-\pi_{1}\frac{\|L_{K_{\mathcal{M}}}^{\frac{1}{2}}\mu_{1}\|^{2}}{\lambda+\tau_{l}}-\pi_{2}\frac{\|L_{K_{\mathcal{M}}}^{\frac{1}{2}}\mu_{2}\|^{2}}{\lambda+\tau_{l}}\right)$
$\displaystyle\leq\frac{1}{n}\operatorname{D}(\lambda).$

Moreover, by appealing to the Markov inequality, we have that with confidence
at least $1-\delta/2$

$\|(T+\lambda
I)^{-\frac{1}{2}}L^{\frac{1}{2}}_{K_{\mathcal{M}}}\left(\hat{d}-d\right)\|\leq\frac{2}{\delta}\sqrt{\frac{\operatorname{D}(\lambda)}{n}}\leq\frac{1}{\kappa\delta}B_{n,\lambda},$
(28)

with
$B_{n,\lambda}=\frac{2\kappa}{\sqrt{n}}\left(\frac{\kappa}{\sqrt{n\lambda}}+\sqrt{D(\lambda)}\right)$
and $\kappa=M_{2}\|L^{\frac{1}{2}}_{K_{\mathcal{M}}}\|_{\operatorname{op}}$.

Thanks to the inequality in (28) and Theorem B.3, we have that with
probability at least $1-\delta$

$\displaystyle I_{1}=\left\|T^{\frac{1}{2}}(T_{n}+\lambda
I)^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}\left(\frac{1}{n}\sum_{i=1}^{n}y_{i}x_{i}-d\right)\right\|\leq\left(\frac{B_{n,\lambda}\log(2/\delta)}{\lambda}+1\right)^{2}\left(\frac{1}{2\kappa\delta}B_{n,\lambda}\right).$
(29)

We now turn to the bias term $I_{2}$ in (B).

###  B.2 Bias

By simple calculations, we have

$\displaystyle I_{2}$
$\displaystyle=\left\|T^{\frac{1}{2}}\left[(T_{n}+\lambda
I)^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}d-T^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}d\right]\right\|$
$\displaystyle=\left\|T^{\frac{1}{2}}(T_{n}+\lambda I)^{-1}[T-T_{n}+\lambda
I]\right\|_{\text{op}}\|f_{0}\|$ $\displaystyle\leq\left\|(T+\lambda
I)^{\frac{1}{2}}(T_{n}+\lambda
I)^{-\frac{1}{2}}\right\|_{\text{op}}\left\|(T_{n}+\lambda
I)^{-\frac{1}{2}}(T+\lambda
I)^{\frac{1}{2}}\right\|_{\text{op}}\left\|(T+\lambda
I)^{-\frac{1}{2}}[T-T_{n}+\lambda I]\right\|_{\text{op}}\|f_{0}\|$
$\displaystyle\leq\left\|(T+\lambda I)(T_{n}+\lambda
I)^{-1}\right\|_{\text{op}}\left\|(T+\lambda I)^{-\frac{1}{2}}[T-T_{n}+\lambda
I]\right\|_{\text{op}}\|f_{0}\|$ $\displaystyle\leq\left\|(T+\lambda
I)(T_{n}+\lambda I)^{-1}\right\|_{\text{op}}\left\|(T+\lambda
I)^{-\frac{1}{2}}[T-T_{n}]\right\|_{\text{op}}\|f_{0}\|$
$\displaystyle+\left\|(T+\lambda I)(T_{n}+\lambda
I)^{-1}\right\|_{\text{op}}\lambda\left\|(T+\lambda
I)^{-\frac{1}{2}}\right\|_{\text{op}}\|f_{0}\|.$

Moreover, by using Theorems B.2-B.3, and the inequality

$\lambda\left\|(T+\lambda
I)^{-\frac{1}{2}}\right\|_{\text{op}}\leq\frac{\lambda}{\sqrt{\lambda}}=\sqrt{\lambda},$
(30)

we have that with probability at least $1-\delta$

$\displaystyle\left\|T^{\frac{1}{2}}\left[(T_{n}+\lambda
I)^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}d-T^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}d\right]\right\|\leq\left(\frac{B_{n,\lambda}\log(4/\delta)}{\sqrt{\lambda}}+1\right)^{2}(B_{n,\lambda}\log(4/\delta)+\sqrt{\lambda})\|f_{0}\|$
(31)

###  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:

$\displaystyle\left\|T^{\frac{1}{2}}(\hat{f}_{\lambda}-f_{0})\right\|^{2}$
(32) $\displaystyle\leq 2\left\|T^{\frac{1}{2}}(T_{n}+\lambda
I)^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}\left(\hat{d}-d\right)\right\|^{2}+2\left\|T^{\frac{1}{2}}\left[(T_{n}+\lambda
I)^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}d-T^{-1}L_{K_{\mathcal{M}}}^{\frac{1}{2}}d\right]\right\|^{2}$
(33) $\displaystyle\leq
2\left(\frac{B_{n,\lambda}\log(2/\delta)}{\sqrt{\lambda}}+1\right)^{4}\left(\frac{1}{2\kappa\delta}B_{n,\lambda}\right)^{2}$
(34)
$\displaystyle+2\left(\frac{B_{n,\lambda}\log(4/\delta)}{\sqrt{\lambda}}+1\right)^{4}\left(B_{n,\lambda}\log(4/\delta)+\sqrt{\lambda}\right)^{2}\|f_{0}\|^{2}$
(35) $\displaystyle\leq
2\frac{\lambda}{\delta^{2}}\left(\frac{B_{n,\lambda}\log(2/\delta)}{\sqrt{\lambda}}+1\right)^{4}\left(\frac{1}{2\kappa}\frac{B_{n,\lambda}}{\sqrt{\lambda}}\right)^{2}$
(36)
$\displaystyle+2\lambda\left(\frac{B_{n,\lambda}\log(4/\delta)}{\sqrt{\lambda}}+1\right)^{4}\left(\frac{B_{n,\lambda}}{\sqrt{\lambda}}\log(4/\delta)+1\right)^{2}\|f_{0}\|^{2}$
(37) $\displaystyle\leq
C\frac{\left(\log(4/\delta)\right)^{6}}{\delta^{2}}n^{-\frac{1}{1+\theta}},$
(38)

where in the last inequality we have chosen $\lambda=n^{-\frac{1}{1+\theta}}$
and used the simple inequality $B_{n,\lambda}\leq
2\kappa(\kappa+\sqrt{c})\sqrt{\lambda}$. 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 $\mathcal{H}$ be a Hilbert space endowed with a norm
$\|\cdot\|_{\mathcal{H}}$ and let $X$ be a random variable taking values in
$\mathcal{H}$. Let $\xi_{1},\ldots,\xi_{n}$ be a sequence of $n$ independent
copies of $X$. Assume that $\|\xi\|_{\mathcal{H}}\leq M$ (a.s.), then for
$0<\delta<1$

$\left\|\frac{1}{n}\sum_{i=1}^{n}\left(\xi_{i}-\mathbb{E}[\xi]\right)\right\|^{2}_{\mathcal{H}}\leq\frac{2M\log(2/\delta)}{n}+\sqrt{\frac{2\mathbb{E}\left[\|\xi\|_{\mathcal{H}}^{2}\right]\log(2/\delta)}{n}}$

with probability at least $1-\delta$.

###### Proof.

A proof can be found in [22]. ∎

Using Theorem B.1 it can be shown that the following two theorems hold.

######  Theorem B.2.

Under Assumption 3.1, for any $0<\delta<1$, the inequality

$\|(T+\lambda I)^{-\frac{1}{2}}(T-T_{n})\|_{\operatorname{op}}\leq
B_{n,\lambda}\log(2/\delta),$ (39)

holds with confidence at least $1-\delta$.

###### Proof.

A proof can be found in [63]. ∎

######  Theorem B.3.

Under Assumption 3.1, for any $0<\delta<1$ with confidence at least
$1-\delta$,

$\|(T+\lambda I)(T_{n}+\lambda
I)^{-1}\|_{\operatorname{op}}\leq\left(\frac{B_{n,\lambda}\log(2/\delta)}{\sqrt{\lambda}}+1\right)^{2}.$
(40)

Moreover, the confidence set is the same as that in Theorem B.2.

###### Proof.

A proof can be found in [63]. ∎

###  B.5 Multivariate model

Thanks to the multivariate notation introduced, we can reformulate the
multivariate model (13) as

$\text{minimize}_{\bm{f}}\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}-\langle\bm{x}_{i},L_{\bm{K}}^{\frac{1}{2}}\bm{f}\rangle_{\mathcal{H}}\right)^{2}+\lambda\|\bm{f}\|^{2}_{\mathcal{H}},$
(41)

with $\bm{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.

##  Appendix 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 [27] 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 [49] to perform dimension reduction in
the FPCA+LDA approach.

We use a triangulated surface $\mathcal{M}_{\mathcal{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,\ldots,40\\}$ consisting of
$40$ eigenfunctions of the Laplace-Beltrami operator. Then, we generate two
sets of smooth functional data supported on $\mathcal{M}_{\mathcal{T}}$ as
follows:

$\begin{aligned} &\text{Group 1:}\quad
x_{i1}=w_{i1}v_{1}+w_{i2}v_{2}+\ldots+w_{i40}v_{40}\quad i=1,\ldots,n\\\
&\text{Group 2:}\quad
x_{i2}=\alpha\mu+u_{i1}v_{1}+u_{i2}v_{2}+\ldots+u_{i40}v_{40}\quad
i=1,\ldots,n\end{aligned},$ (42)

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
$\sigma_{j}^{2}$ decreasing in $j$. The function $\mu$ is the group
difference, chosen to be the 10th eigenfunction, and its magnitude is
controlled by a parameter $\alpha>0$. The parameter $\alpha$ controls the
‘difficulty’ of the problem. We then identify three regimes, with
$\alpha\in{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 $20$K 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 $\alpha=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 $\alpha=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 $\alpha=0.2$ and $\alpha=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 $\alpha=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\in\\{128,256,512,1024\\}$, $50$ datasets were generated.

## References

  * [1] Isaac Fried and David S. Malkus  “Finite element mass matrix lumping by numerical integration with no convergence rate loss”  In _International Journal of Solids and Structures_ 11.4, 1975, pp. 461–466  DOI: 10.1016/0020-7683(75)90081-5
  * [2] E. Hinton, T. Rock and O. C. Zienkiewicz  “A note on mass lumping and related processes in the finite element method”  In _Earthquake Engineering & Structural Dynamics_ 4.3, 1976, pp. 245–249  DOI: 10.1002/eqe.4290040305
  * [3] Christopher C. Paige and Michael A. Saunders  “Algorithm 583: LSQR: Sparse Linear Equations and Least Squares Problems”  In _ACM Transactions on Mathematical Software (TOMS)_ 8.2, 1982, pp. 195–209  DOI: 10.1145/355993.356000
  * [4] David G. Kendall  “Shape manifolds, procrustean metrics, and complex projective spaces”  In _Bulletin of the London Mathematical Society_ 16.2, 1984, pp. 81–121  DOI: 10.1112/blms/16.2.81
  * [5] William Bondareff, C. Q. Mountjoy, M. Roth and D. L. Hauser  “Neurofibrillary degeneration and neuronal loss in alzheimer’s disease”  In _Neurobiology of Aging_ 10.6, 1989, pp. 709–715  DOI: 10.1016/0197-4580(89)90007-9
  * [6] Grace Wahba  “Spline Models for Observational Data”  In _Spline Models for Observational Data_ Society for IndustrialApplied Mathematics, 1990  DOI: 10.1137/1.9781611970128
  * [7] Mike West and Jeff Harrison  “Inference for Functional Data with Applications”  In _Bayesian Forecasting and Dynamic Models_ 200, Springer Series in Statistics  New York, NY: Springer New York, 1997, pp. 1–694  DOI: 10.1007/978-1-4614-3655-3
  * [8] Ulf Grenander and Michael I. Miller  “Computational anatomy: an emerging discipline”  In _Quarterly of Applied Mathematics_ 56.4, 1998, pp. 617–694  DOI: 10.1090/qam/1668732
  * [9] Bruce Fischl, Martin I. Sereno and Anders M. Dale  “Cortical surface-based analysis: II. Inflation, flattening, and a surface-based coordinate system”  In _NeuroImage_ 9.2, 1999, pp. 195–207  DOI: 10.1006/nimg.1998.0396
  * [10] Bruce Fischl, Martin I. Sereno, Roger B.H. Tootell and Anders M. Dale  “High-resolution intersubject averaging and a coordinate system for the cortical surface”  In _Human Brain Mapping_ 8.4, 1999, pp. 272–284  DOI: 10.1002/(SICI)1097-0193(1999)8:4¡272::AID-HBM10¿3.0.CO;2-4
  * [11] Bruce Fischl and Anders M. Dale  “Measuring the thickness of the human cerebral cortex from magnetic resonance images”  In _Proceedings of the National Academy of Sciences of the United States of America_ 97.20 National Academy of Sciences, 2000, pp. 11050–11055  DOI: 10.1073/pnas.200033797
  * [12] Gareth M. James and Trevor J. Hastie  “Functional linear discriminant analysis for irregularly sampled curves”  In _Journal of the Royal Statistical Society. Series B: Statistical Methodology_ 63.3, 2001, pp. 533–550  DOI: 10.1111/1467-9868.00297
  * [13] Paul M. Thompson et al.  “Cortical change in Alzheimer’s disease detected with a disease-specific population-based brain atlas”  In _Cerebral Cortex_ 11.1, 2001, pp. 1–16  DOI: 10.1093/cercor/11.1.1
  * [14] Felipe Cucker and Steve Smale  “On the mathematical foundations of learning”  In _Bulletin of the American Mathematical Society_ 39.1, 2002, pp. 1–49  DOI: 10.1090/S0273-0979-01-00923-5
  * [15] Alain Berlinet and Christine Thomas-Agnan  “Reproducing Kernel Hilbert Spaces in Probability and Statistics”  In _Reproducing Kernel Hilbert Spaces in Probability and Statistics_ Boston, MA: Springer US, 2004  DOI: 10.1007/978-1-4419-9096-9
  * [16] M. Vaillant, M. I. Miller, L. Younes and A. Trouvé  “Statistics on diffeomorphisms via tangent space representations”  In _NeuroImage_ 23.SUPPL. 1, 2004, pp. S161–S169  DOI: 10.1016/j.neuroimage.2004.07.023
  * [17] Hans Georg Müller  “Functional modelling and classification of longitudinal data”  In _Scandinavian Journal of Statistics_ 32.2, 2005, pp. 223–240  DOI: 10.1111/j.1467-9469.2005.00429.x
  * [18] Marc Vaillant and Joan Glaunès  “Surface matching via currents”  In _Lecture Notes in Computer Science_ 3565 Berlin, Heidelberg: Springer, 2005, pp. 381–392  DOI: 10.1007/11505730˙32
  * [19] Heiko Braak et al.  “Staging of Alzheimer disease-associated neurofibrillary pathology using paraffin sections and immunocytochemistry”  In _Acta Neuropathologica_ 112.4, 2006, pp. 389–404  DOI: 10.1007/s00401-006-0127-z
  * [20] Frédéric Ferraty and Philippe Vieu  “Nonparametric Functional Data Analysis”  In _Nonparametric Functional Data Analysis_ , Springer Series in Statistics  Springer New York, 2006  DOI: 10.1007/0-387-36620-2
  * [21] Michael I. Miller, Alain Trouvé and Laurent Younes  “Geodesic shooting for computational anatomy”  In _Journal of Mathematical Imaging and Vision_ 24.2, 2006, pp. 209–228  DOI: 10.1007/s10851-005-3624-0
  * [22] Iosif Pinelis  “Optimum Bounds for the Distributions of Martingales in Banach Spaces”  In _The Annals of Probability_ 22.4, 2007, pp. 347–370  DOI: 10.1214/aop/1176988477
  * [23] Cristian Preda, Gilbert Saporta and Caroline Lévéder  “PLS classification of functional data”  In _Computational Statistics_ 22.2, 2007, pp. 223–235  DOI: 10.1007/s00180-007-0041-4
  * [24] Marc Vaillant, Anqi Qiu, Joan Glaunès and Michael I. Miller  “Diffeomorphic metric surface mapping in subregion of the superior temporal gyrus”  In _NeuroImage_ 34.3, 2007, pp. 1149–1159  DOI: 10.1016/j.neuroimage.2006.08.053
  * [25] Hyejin Shin  “An extension of Fisher’s discriminant analysis for stochastic processes”  In _Journal of Multivariate Analysis_ 99.6, 2008, pp. 1191–1216  DOI: 10.1016/j.jmva.2007.08.001
  * [26] Bradford C Dickerson et al.  “The cortical signature of Alzheimer’s disease: Regionally specific cortical thinning relates to symptom severity in very mild to mild AD dementia and is detectable in asymptomatic amyloid-positive individuals”  In _Cerebral Cortex_ 19.3, 2009, pp. 497–510  DOI: 10.1093/cercor/bhn113
  * [27] Trevor Hastie, Robert Tibshirani and Jerome Friedman  “The Elements of Statistical Learning”, Springer Series in Statistics  New York, NY: Springer New York, 2009  DOI: 10.1007/978-0-387-84858-7
  * [28] Alfio Quarteroni  “Numerical Models for Differential Problems”  In _Modeling, Simulation and Applications_ 2, 2009, pp. 1–593  DOI: 10.1007/978-88-470-1071-0
  * [29] Gilles Blanchard and Nicole Krämer  “Optimal learning rates for Kernel Conjugate Gradient regression”  In _Advances in Neural Information Processing Systems 23: 24th Annual Conference on Neural Information Processing Systems 2010, NIPS 2010_ , 2010  DOI: 10.1.1.231.648
  * [30] Ming Yuan and T. Tony Cai  “A reproducing Kernel Hilbert space approach to functional linear regression”  In _Annals of Statistics_ 38.6, 2010, pp. 3412–3444  DOI: 10.1214/09-AOS772
  * [31] Haim Brezis  “Functional Analysis, Sobolev Spaces and Partial Differential Equations”  In _Functional Analysis, Sobolev Spaces and Partial Differential Equations_ New York, NY: Springer New York, 2011  DOI: 10.1007/978-0-387-70914-7
  * [32] Mert R. Sabuncu et al.  “The dynamics of cortical and hippocampal atrophy in Alzheimer disease”  In _Archives of Neurology_ 68.8, 2011, pp. 1040–1048  DOI: 10.1001/archneurol.2011.167
  * [33] B. T. Thomas Yeo et al.  “The organization of the human cerebral cortex estimated by intrinsic functional connectivity”  In _Journal of Neurophysiology_ 106.3, 2011, pp. 1125–1165  DOI: 10.1152/jn.00338.2011
  * [34] Aurore Delaigle and Peter Hall  “Achieving near perfect classification for functional data”  In _Journal of the Royal Statistical Society. Series B: Statistical Methodology_ 74.2, 2012, pp. 267–286  DOI: 10.1111/j.1467-9868.2011.01003.x
  * [35] Qing Mai, Hui Zou and Ming Yuan  “A direct approach to sparse discriminant analysis in ultra-high dimensions”  In _Biometrika_ 99.1, 2012, pp. 29–42  DOI: 10.1093/biomet/asr066
  * [36] T. Tony Cai and Ming Yuan  “Minimax and adaptive prediction for functional linear regression”  In _Journal of the American Statistical Association_ 107.499, 2012, pp. 1201–1216  DOI: 10.1080/01621459.2012.716337
  * [37] Gregory E. Fasshauer and Qi Ye  “Reproducing kernels of Sobolev spaces via a green kernel approach with differential operators and boundary operators”  In _Advances in Computational Mathematics_ 38.4, 2013, pp. 891–921  DOI: 10.1007/s10444-011-9264-6
  * [38] Tailen Hsing and Randall Eubank  “Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators”  In _Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators_ , Wiley Series in Probability and Statistics  Chichester, UK: John Wiley & Sons, Ltd, 2013, pp. 1–333  DOI: 10.1002/9781118762547
  * [39] Stephen M. Smith et al.  “Resting-state fMRI in the Human Connectome Project”  In _NeuroImage_ 80, 2013, pp. 144–168  DOI: 10.1016/j.neuroimage.2013.05.039
  * [40] Irene Epifanio and Noelia Ventura-Campos  “Hippocampal shape analysis in Alzheimer’s disease using functional data analysis”  In _Statistics in Medicine_ 33.5, 2014, pp. 867–880  DOI: 10.1002/sim.5968
  * [41] Jeff Goldsmith, Lei Huang and Ciprian M Crainiceanu  “Smooth scalar-on-image regression via spatial Bayesian variable selection”  In _Journal of Computational and Graphical Statistics_ 23.1, 2014, pp. 46–64  DOI: 10.1080/10618600.2012.743437
  * [42] Jingyong Su, Sebastian Kurtek, Eric Klassen and Anuj Srivastava  “Statistical analysis of trajectories on riemannian manifolds: Bird migration, hurricane tracking and video surveillance”  In _Annals of Applied Statistics_ 8.1, 2014, pp. 530–552  DOI: 10.1214/13-AOAS701
  * [43] B. Charlier, G. Nardi and A. Trouvé  “The matching problem between functional shapes via a BV-penalty term: a $\Gamma$-convergence result”, 2015, pp. 1–31  arXiv: http://arxiv.org/abs/1503.07685
  * [44] Aasa Feragen, François Lauze and Søren Hauberg  “Geodesic exponential kernels: When curvature and linearity conflict”  In _Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ 07-12-June.1 IEEE, 2015, pp. 3032–3042  DOI: 10.1109/CVPR.2015.7298922
  * [45] Sadeep Jayasumana et al.  “Kernel Methods on Riemannian Manifolds with Gaussian RBF Kernels”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 37.12, 2015, pp. 2464–2477  DOI: 10.1109/TPAMI.2015.2414422
  * [46] James O. Ramsay and Bernard W. Silverman  “Functional Data Analysis”  In _International Encyclopedia of the Social & Behavioral Sciences: Second Edition_, Springer Series in Statistics  New York: Springer-Verlag, 2015, pp. 514–518  DOI: 10.1016/B978-0-08-097086-8.42046-5
  * [47] Ian L. Dryden and Kanti V. Mardia  “Statistical shape analysis, with applications in R: Second edition”  In _Statistical Shape Analysis, with Applications in R: Second Edition_ , Wiley Series in Probability and Statistics  Chichester, UK: John Wiley & Sons, Ltd, 2016, pp. 1–454  DOI: 10.1002/9781119072492
  * [48] Clifford R. Jack et al.  “A/T/N: An unbiased descriptive classification scheme for Alzheimer disease biomarkers”  In _Neurology_ 87.5, 2016, pp. 539–547  DOI: 10.1212/WNL.0000000000002923
  * [49] Eardi Lila, John A.D. Aston and Laura M. Sangalli  “Smooth principal component analysis over two-dimensional manifolds with an application to neuroimaging”  In _Annals of Applied Statistics_ 10.4, 2016, pp. 1854–1879  DOI: 10.1214/16-AOAS975
  * [50] B. Charlier, N. Charon and A. Trouvé  “The Fshape Framework for the Variability Analysis of Functional Shapes”  In _Foundations of Computational Mathematics_ 17.2, 2017, pp. 287–357  DOI: 10.1007/s10208-015-9288-2
  * [51] Xiongtao Dai, Hans Georg Müller and Fang Yao  “Optimal Bayes classifiers for functional data and density ratios”  In _Biometrika_ 104.3, 2017, pp. 545–560  DOI: 10.1093/biomet/asx024
  * [52] Sieun Lee et al.  “Atlas-based shape analysis and classification of retinal optical coherence tomography images using the functional shape (fshape) framework”  In _Medical Image Analysis_ 35 Elsevier B.V., 2017, pp. 570–581  DOI: 10.1016/j.media.2016.08.012
  * [53] Xiao Wang and Hongtu Zhu  “Generalized Scalar-on-Image Regression Models via Total Variation”  In _Journal of the American Statistical Association_ 112.519 Taylor & Francis, 2017, pp. 1156–1168  DOI: 10.1080/01621459.2016.1194846
  * [54] José R. Berrendero, Antonio Cuevas and José L. Torrecilla  “On the Use of Reproducing Kernel Hilbert Spaces in Functional Classification”  In _Journal of the American Statistical Association_ 113.523 American Statistical Association, 2018, pp. 1210–1218  DOI: 10.1080/01621459.2017.1320287
  * [55] Carlo Biffi et al.  “Three-dimensional cardiovascular imaging-genetics: A mass univariate framework”  In _Bioinformatics_ 34.1, 2018, pp. 97–103  DOI: 10.1093/bioinformatics/btx552
  * [56] Lu Hung Chen and Ci Ren Jiang  “Sensible functional linear discriminant analysis”  In _Computational Statistics and Data Analysis_ 126 Elsevier B.V., 2018, pp. 39–52  DOI: 10.1016/j.csda.2018.04.005
  * [57] Xiongtao Dai and Hans Georg Müller  “Principal component analysis for functional data on riemannian manifolds and spheres”  In _Annals of Statistics_ 46.6B, 2018, pp. 3334–3361  DOI: 10.1214/17-AOS1660
  * [58] A. H. Gee, G. M. Treece and K. E.S. Poole  “How does the femoral cortex depend on bone shape? A methodology for the joint analysis of surface texture and shape”  In _Medical Image Analysis_ 45 Elsevier B.V., 2018, pp. 55–67  DOI: 10.1016/j.media.2018.01.001
  * [59] Clara Happ and Sonja Greven  “Multivariate Functional Principal Component Analysis for Data Observed on Different (Dimensional) Domains”  In _Journal of the American Statistical Association_ 113.522, 2018, pp. 649–659  DOI: 10.1080/01621459.2016.1273115
  * [60] Jian Kang, Brian J Reich and Ana Maria Staicu  “Scalar-on-image regression via the soft-thresholded Gaussian process”  In _Biometrika_ 105.1, 2018, pp. 165–184  DOI: 10.1093/biomet/asx075
  * [61] Matthew Reimherr, Bharath Sriperumbudur and Bahaeddine Taoufik  “Optimal prediction for additive function-on-function regression”  In _Electronic Journal of Statistics_ 12.2, 2018, pp. 4571–4601  DOI: 10.1214/18-EJS1505
  * [62] Xiaoxiao Sun, Pang Du, Xiao Wang and Ping Ma  “Optimal Penalized Function-on-Function Regression Under a Reproducing Kernel Hilbert Space Framework”  In _Journal of the American Statistical Association_ 113.524, 2018, pp. 1601–1611  DOI: 10.1080/01621459.2017.1356320
  * [63] Hongzhi Tong and Michael Ng  “Analysis of regularized least squares for functional linear regression model”  In _Journal of Complexity_ 49 Elsevier Inc., 2018, pp. 85–94  DOI: 10.1016/j.jco.2018.08.001
  * [64] David Kraus and Marco Stefanucci  “Classification of functional fragments by regularized linear classifiers with domain selection”  In _Biometrika_ 106.1, 2019, pp. 161–180  DOI: 10.1093/biomet/asy060
  * [65] Zhenhua Lin and Fang Yao  “Intrinsic riemannian functional data analysis1”  In _Annals of Statistics_ 47.6, 2019, pp. 3533–3577  DOI: 10.1214/18-AOS1787
  * [66] Laurent Younes  “Shapes and diffeomorphisms”  In _Applied Mathematical Sciences (Switzerland)_ 171, Applied Mathematical Sciences  Berlin, Heidelberg: Springer, 2019, pp. 1–316  DOI: 10.1007/978-3-642-12055-8
  * [67] Paromita Dubey and Hans Georg Müller  “Functional models for time-varying random objects”  In _Journal of the Royal Statistical Society. Series B: Statistical Methodology_ 82.2, 2020, pp. 275–327  DOI: 10.1111/rssb.12337
  * [68] Xiangnan Feng, Tengfei Li, Xinyuan Song and Hongtu Zhu  “Bayesian Scalar on Image Regression With Nonignorable Nonresponse”  In _Journal of the American Statistical Association_ 115.532 ApplicationsCase Studies, 2020, pp. 1574–1597  DOI: 10.1080/01621459.2019.1686391
  * [69] L. Ferrando, N. Ventura-Campos and I. Epifanio  “Detecting and visualizing differences in brain structures with SPHARM and functional data analysis”  In _NeuroImage_ 222.August Elsevier Inc., 2020, pp. 117209  DOI: 10.1016/j.neuroimage.2020.117209
  * [70] Irina Gaynanova  “Prediction and estimation consistency of sparse multi-class penalized optimal scoring”  In _Bernoulli_ 26.1, 2020, pp. 286–322  DOI: 10.3150/19-BEJ1126
  * [71] Jingyi Jessica Li and Xin Tong  “Statistical Hypothesis Testing versus Machine Learning Binary Classification: Distinctions and Guidelines”  In _Patterns_ 1.7, 2020, pp. 100115  DOI: 10.1016/j.patter.2020.100115
  * [72] Eardi Lila and John A.D. Aston  “Statistical Analysis of Functions on Surfaces, With an Application to Medical Imaging”  In _Journal of the American Statistical Association_ 115.531, 2020, pp. 1420–1434  DOI: 10.1080/01621459.2019.1635479
  * [73] Amanda F. Mejia et al.  “A Bayesian General Linear Modeling Approach to Cortical Surface fMRI Data Analysis”  In _Journal of the American Statistical Association_ 115.530, 2020, pp. 501–520  DOI: 10.1080/01621459.2019.1611582
  * [74] Moo K Chung, Jamie L Hanson and Seth D Pollak  “Statistical Analysis on Brain Surfaces”  In _Handbook of Neuroimaging Data Analysis_ , 2021, pp. 273–302  DOI: 10.1201/9781315373652-20
  * [75] Kwang Rae Kim, Ian L. Dryden, Huiling Le and Katie E. Severn  “Smoothing splines on Riemannian manifolds, with applications to 3D shape space”  In _Journal of the Royal Statistical Society. Series B: Statistical Methodology_ 83.1, 2021, pp. 108–132  DOI: 10.1111/rssb.12402
  * [76] Juhyun Park, Jeongyoun Ahn and Yongho Jeon  “Sparse functional linear discriminant analysis”  In _Biometrika_ , 2021  DOI: 10.1093/biomet/asaa107
  * [77] Shan Yu, Guannan Wang, Li Wang and Lijian Yang  “Multivariate spline estimation and inference for image-on-scalar regression”  In _Statistica Sinica_ 31.3, 2021, pp. 1463–1487  DOI: 10.5705/ss.202019.0188