Multiclass classification for multidimensional functional data through
  deep neural networks
The intrinsically infinite-dimensional features of the functional
observations over multidimensional domains render the standard classification
methods effectively inapplicable. To address this problem, we introduce a novel
multiclass functional deep neural network (mfDNN) classifier as an innovative
data mining and classification tool. Specifically, we consider sparse deep
neural network architecture with rectifier linear unit (ReLU) activation
function and minimize the cross-entropy loss in the multiclass classification
setup. This neural network architecture allows us to employ modern
computational tools in the implementation. The convergence rates of the
misclassification risk functions are also derived for both fully observed and
discretely observed multidimensional functional data. We demonstrate the
performance of mfDNN on simulated data and several benchmark datasets from
different application domains.
Shuoyang Wang, Guanqun Cao
# Multiclass classification for multidimensional functional data through deep neural networks

###### Abstract

The intrinsically infinite-dimensional features of the functional observations over multidimensional domains render the standard classification methods effectively inapplicable. To address this problem, we introduce a novel multiclass functional deep neural network (mfDNN) classifier as an innovative data mining and classification tool. Specifically, we consider sparse deep neural network architecture with rectifier linear unit (ReLU) activation function and minimize the cross-entropy loss in the multiclass classification setup. This neural network architecture allows us to employ modern computational tools in the implementation. The convergence rates of the misclassification risk functions are also derived for both fully observed and discretely observed multidimensional functional data. We demonstrate the performance of mfDNN on simulated data and several benchmark datasets from different application domains.

P 62G05, 62G08; secondary 62G35. [table]capposition=top

**MSC2020 subject classifications:** Primary 62G05, 62G08; secondary 62G35.

Functional data analysis, Deep Neural networks, Multiclass classification, Rate of convergence, Multidimensional functional data.

## 1 Introduction

Functional data classification has wide applications in many areas such as machine learning and artificial intelligence [19, 10, 17, 5]. While the majority of the work on functional data classification focuses on one-dimensional functional data cases, for example, temperature data over a certain period [16], and speech recognition data (log-periodogram data) [11, 23], and there are few results obtained for multidimensional functional data classification, for example, 2D or 3D images classification. To the best of our knowledge, there is only one recent work [22] where binary classification for general non-Gaussian multidimensional functional data was considered. However, their proposed classifier is designed for fully observed functional data and only able to conduct binary classification. There is also a lack of literature working on multiclass classification for functional data. [11] investigated a Bayesian approach to estimate parameters in multiclass functional models. However, their method is still only suitable for fully observed one-dimensional functional data classification.

The availability of massive databases resulted in the development of scalable machine learning methods to process these data. In this paper, we are interested in multiclass classification for multidimensional functional data, including but not limited to
2D and 3D imaging data classification, in the framework of DNN. The existing statistical and machine learning methods belong to the general framework of multivariate analysis. That is, data are treated as vectors of discrete samples and permutation of components will not affect the analysis results, hence the ordering of the pixels or voxels is irrelevant in these analyses and the imaging data were treated as high-dimensional data. In recent years, many sparse discriminant analysis methods have been proposed for i.i.d. high-dimensional data classification and variable selection. Most of these proposals focus on binary classification and they are not directly applicable to multiclass classification problems. A popular multiclass sparse discriminant analysis proposal is the \(\ell_{1}\) penalized Fisher's discriminant [24]. However, [24] does not have theoretical justifications. It is generally unknown how close their estimated discriminant directions are to the true directions, and whether the final classifier will work similarly as the Bayes rule. More recently, [14] proposed a multiclass sparse discriminant analysis method that estimates all discriminant directions simultaneously with theoretical justification for high-dimensional data setting. We compare our method with both [24] and [14] carefully on numerical performance in Sections 5 and 6.

An efficient way to handle multi-dimensional functional data, for example, imaging data, is to incorporate the information that is inherent in order and smoothness of processes over pixels or voxels. In this paper, we propose a novel approach, called as multi-class functional deep neural network (mfDNN, multiclass-functional+DNN), for multidimensional functional data multiclass classification. We extract the functional principal components scores (FPCs) of the data functions, and then train a DNN based classifier on these FPCs as well as their corresponding class memberships. Given these FPCs as inputs and their corresponding class labels as outputs, we consider sparse deep ReLU network architecture and minimize the cross-entropy (CE) loss in the multiclass classification setup. CE loss is a standard tool in machine learning classification, but the technique has not been utilized to its fullest potential in the context of functional data classification, especially regarding the study of theoretical properties of DNN classifiers. A recent work [2], has an example of multiclassification by minimizing CE loss in sparse DNN for cross-sectional data setting. While the powerful capabilities of DNN and CE loss have been clearly demonstrated for i.i.d. data, how to best adapt these models to functional data remains under-explored. We apply the same loss function in [2] but develop it in the context of multidimensional functional multiclassification. Furthermore, convergence rates and the conditional class probabilities are developed under both fully observed and discretely observed multi-dimensional functional data setting, respectively.

The rest of this article is organized as follows. In Section 2, we introduce the \(K\)-class functional data classification models and propose the multiclass functional deep neural network classifiers. In Section 3, we establish theoretical properties of mfDNN under suitable technical assumptions. Section 4 provides two progressive examples to demonstrate the validity of these technical assumptions. In Section 5, performances of mfDNN and its competitors are demonstrated through simulation studies. The corresponding R programs for implementing our method are provided on GitHub. In Section 6, we apply mfDNN to the zip code dataset and Alzeheimer's Disease data. Section 7 summarizes the conclusions. Technical proofs are provided in Appendix.
## 2 Methodology

### \(K\)-class functional data classification models

Suppose we observe \(n\) i.i.d. functional training samples \(\{(X_{i}(s),\mathbf{Y}_{i}):1\leq i\leq n,s\in\left[0,1\right]^{d}\}\), which are independent of \(X(s)\) to be classified, and class label \(\mathbf{Y}_{i}=(Y_{i1},\ldots,Y_{iK})^{\mathsf{T}}\), such that \(Y_{ik}=1\) if the sample is in the \(k\)-th group, \(k=1,2,\ldots,K\), and \(0\) otherwise. Throughout the entire paper, we consider the number of class \(K\) is a universal constant not depending on sample size \(n\). In multiclass classification with \(K\geq 2\) classes, we are interested in grouping a newly observed continuous curves \(X(s)\) to one of the \(K\) classes given the training data.

If \(Y_{ik}=1\), \(X_{i}(s)\) presumably has some unknown mean function \(\mathbb{E}X_{i}(s)=\mu_{k}(s)\) and unknown covariance function \(\Omega_{k}(s,s^{\prime})=\mathbb{E}\left[\left(X_{i}(s)-\mu_{k}(s)\right)(X_{ i}(s^{\prime})-\mu_{k}(s^{\prime}))\right]\), for \(s,s^{\prime}\in\left[0,1\right]^{d}\) and \(1\leq k\leq K\). Suppose the covariance function \(\Omega_{k}(\cdot,\cdot)\) satisfies the Karhunen-Loeve decomposition:

\[\Omega_{k}(s,s^{\prime})=\sum_{j=1}^{\infty}\lambda_{kj}\psi_{kj}(s)\psi_{kj} (s^{\prime}),s,s^{\prime}\in\left[0,1\right]^{d}, \tag{1}\]

where \(\psi_{kj}\), \(j\geq 1\) is an orthonormal basis of \(L^{2}(\left[0,1\right]^{d})\) with respect to the usual \(L^{2}\) inner product, and \(\lambda_{k1}\geq\lambda_{k2}\geq\cdots>0\) are nonincreasing positive eigenvalues. For any \(X(s)\in L^{2}(\left[0,1\right]^{d})\) under the \(k\)-th class, write \(X(s)=\sum_{j=1}^{\infty}\xi_{j}\psi_{kj}(s)\), where \(\xi_{j}\)'s are pairwise uncorrelated random coefficients. By projecting the function space to vector space, we define the conditional probabilities

\[\pi_{k}(\mathbf{x})=\mathrm{P}(\mathbf{Y}=\mathbf{e}_{k}|\mathbf{\xi}=\mathbf{x}),\ \ k\in\{1,\ldots,K\},\]

where \(\mathbf{e}_{k}=(0,\ldots,0,1,0,\ldots,0)^{\mathsf{T}}\) is a \(K\)-dimensional standard basis vector denoting the label of the k-th class is observed and \(\mathbf{\xi}=(\xi_{1},\xi_{2},\ldots)^{\mathsf{T}}\) is an infinite dimension random vector. Notably, \(\sum_{k=1}^{K}\pi_{k}(\mathbf{x})=1\) for any \(\mathbf{x}\in\mathbb{R}^{\infty}\).

### Multiclass functional deep neural network classifier

For \(k\in\{1,\ldots,K\}\), define sample covariance function

\[\widehat{\Omega}_{k}(s,s^{\prime})=\frac{1}{n_{k}}\sum_{i\in I_{k}}(X_{i}(s)- \bar{X}_{k}(s))(X_{i}(s^{\prime})-\bar{X}_{k}(s^{\prime})),\ \ s,s^{\prime}\in\left[0,1\right]^{d},\]

where \(I_{k}\) is the collection of \(i\) such that \(Y_{ik}=1\), \(n_{k}=|I_{k}|\) is the sample size, and \(\bar{X}_{k}(s)=n_{k}^{-1}\sum_{i\in I_{k}}X_{i}(s)\) is the sample mean function of class \(k\). The Karhunen-Loeve decomposition for \(\widehat{\Omega}_{k}\) is

\[\widehat{\Omega}_{k}(s,s^{\prime})=\sum_{j=1}^{\infty}\widehat{\lambda}_{kj} \widehat{\psi}_{kj}(s)\widehat{\psi}_{kj}(s^{\prime}),s,s^{\prime}\in\left[0,1 \right]^{d}.\]
The sample data function \(X_{i}(s)\), \(i=1,\ldots,n\), under \(Y_{ik}=1\), can be represented as \(X_{i}(s)=\sum_{j=1}^{\infty}\widehat{\xi}_{ij}\widehat{\psi}_{kj}(s)\). Intuitively, \(\widehat{\boldsymbol{\xi}}^{(i)}:=(\widehat{\xi}_{i1},\widehat{\xi}_{i2},\ldots)\) is an estimator of \(\boldsymbol{\xi}^{(i)}:=(\xi_{i1},\xi_{i2},\ldots)\), in which \(\xi_{ij}\)'s are unobservable random coefficients of \(X_{i}\) with respect to the population basis \(\psi_{kj}\). Hence, it is natural to design classifiers based on \(\widehat{\boldsymbol{\xi}}^{(i)}\)'s.

Let \(\widehat{\boldsymbol{\xi}}_{J}^{(i)}=(\widehat{\xi}_{i1},\ldots,\widehat{\xi }_{iJ})^{\intercal}\) be the \(J\)-dimensional truncation of \(\widehat{\boldsymbol{\xi}}^{(i)}\) for \(i=1,\ldots,n\). Denote \(h_{k}(\cdot)\) as the conditional probability densities of \(\boldsymbol{\xi}\) given \(\boldsymbol{Y}=\boldsymbol{e}_{k}\), such that \(\pi_{k}(\boldsymbol{x})=\frac{h_{k}\mathbf{P}(\boldsymbol{Y}=\boldsymbol{e}_ {k})}{\sum_{\ell=1}^{K}h_{\ell}\mathbf{P}(\boldsymbol{Y}=\boldsymbol{e}_{\ell })}\). When \(X_{i}\)'s are some random processes with complex structures, such as non-Gaussian processes, one major challenge is the underlying complicated form of \(\{h_{k}\}_{k=1}^{K}\) so that estimation of \(\{\pi_{k}\}_{k=1}^{K}\) is typically difficult. In this section, inspired by the rich approximation power of DNN, we propose a new classifier so-called mfDNN, which can still approach Bayes classifiers closely even when \(h_{k}\)'s are non-Gaussian and complicated.

Let \(\sigma(x)=(x)_{+}\) be the ReLU activation function, and \(K\)-dimensional softmax activation function \(\boldsymbol{\sigma}^{*}(\boldsymbol{y})=\left(\frac{\exp{(y_{1})}}{\sum_{k=1} ^{K}\exp{(y_{k})}},\ldots,\frac{\exp{(y_{K})}}{\sum_{k=1}^{K}\exp{(y_{k})}}\right)\) for \(\boldsymbol{y}=(y_{1},\ldots,y_{K})\in\mathbb{R}^{K}\). For any real vectors \(\boldsymbol{V}=(v_{1},\ldots,v_{w})^{\intercal}\) and \(\boldsymbol{z}=(z_{1},\ldots,z_{w})^{\intercal}\), define the shift activation function \(\sigma_{\boldsymbol{V}}(\boldsymbol{z})=(\sigma(z_{1}-v_{1}),\ldots,\sigma(z _{w}-v_{w}))^{\intercal}\). For \(L\geq 1\), \(\boldsymbol{p}=(p_{1},\ldots,p_{L})\in\mathbb{N}^{L}\), let \(\mathcal{F}(L,J,\boldsymbol{p})\) denote the class of fully connected feedforward DNN with \(J\) inputs, \(L\) hidden layers and, for \(l=1,\ldots,L\), \(p_{l}\) nodes on the \(l\)-th hidden layer. Equivalently, any \(f\in\mathcal{F}(L,J,\boldsymbol{p})\) has an expression

\[f(\boldsymbol{x})=\boldsymbol{\sigma}^{*}\mathbf{W}_{L}\sigma_{\mathbf{V}_{L} }\mathbf{W}_{L-1}\sigma_{\mathbf{V}_{L-1}}\ldots\mathbf{W}_{1}\sigma_{\mathbf{ V}_{1}}\mathbf{W}_{0}\boldsymbol{x},\ \ \ \boldsymbol{x}\in\mathbb{R}^{J}, \tag{2.2}\]

where \(\mathbf{W}_{l}\in\mathbb{R}^{p_{l+1}\times p_{l}}\), for \(l=0,\ldots,L\), are weight matrices, \(\boldsymbol{V}_{l}\in\mathbb{R}^{p_{l}}\), for \(l=1,\ldots,L\), are shift vectors. Here, we adopt the convention that \(p_{0}=J\) and \(p_{L+1}=1\).

Due to the large capacity of the fully connected DNN class and to avoid overparameterization, we consider the following sparse DNN class:

\[\mathcal{F}(L,J,\boldsymbol{p},s)= \tag{2.3}\] \[\left\{f\in\mathcal{F}(L,J,\boldsymbol{p}):\sum_{l=0}^{L}\| \mathbf{W}_{l}\|_{0}+\|\mathbf{V}_{l}\|_{0}\leq s,\max_{0\leq l\leq L}\|\mathbf{ W}_{l}\|_{\infty}\vee\|\mathbf{V}_{l}\|_{\infty}\leq 1\right\},\]

where \(\|\cdot\|_{\infty}\) denotes the maximum-entry norm of a matrix/vector, which is controlled by \(1\), \(\|\cdot\|_{0}\) denotes the number of non-zero elements of a matrix/vector, and \(s\) controls the total number of active neurons.

Given the training data \((\boldsymbol{\xi}_{J}^{(1)},\boldsymbol{Y}_{1}),\ldots,(\boldsymbol{\xi}_{J}^ {(n)},\boldsymbol{Y}_{n})\), let

\[\widehat{\boldsymbol{f}}_{\phi}(\cdot)=\left(\widehat{f}_{\phi}^{(1)},\ldots, \widehat{f}_{\phi}^{(K)}\right)^{\intercal}=\arg\min_{\boldsymbol{f}\in \mathcal{F}(L,J,\boldsymbol{p},s)}-\frac{1}{n}\sum_{i=1}^{n}\phi(\boldsymbol{Y }_{i},\boldsymbol{f}(\boldsymbol{\tilde{\xi}}_{J}^{(i)})), \tag{2.4}\]

where \(\phi(\boldsymbol{x}_{1},\boldsymbol{x}_{2})=\boldsymbol{x}_{1}^{\intercal}\log {(\boldsymbol{x}_{2})}\) denotes the CE loss function.

We then propose the following mfDNN classifier: for \(X\in L^{2}([0,1]^{d})\),

\[\widehat{G}^{mfDNN}(X)=k,\ \ \text{if}\ \widehat{f}_{\phi}^{(k)}(\boldsymbol{\xi}_{J} )\geq\widehat{f}_{\phi}^{(k^{\prime})}(\boldsymbol{\xi}_{J}),\ \ \text{for all}\ \ 1\leq k^{\prime}\neq k\leq K. \tag{2.5}\]

The detailed implementation for the proposed mfDNN classifier is explained in Algorithm 1. The tuning parameters include \(J,L,\boldsymbol{p}\), and \(s\), which theoretically contribute
to the performance of classification. We discuss the related theories in Section 3. In practice, we suggest the following data-splitting approach to select \((J,L,\mathbf{p},s)\). Note that the sparsity is not directly applicable by assigning the number of active neurons. Instead, we use dropout technique to achieve sparse network.

```
Input: Collection of samples \(\left\{\left\{X_{i}(s_{j})\right\}_{j=1}^{m},\mathbf{Y}_{i}\right\}_{i=1}^{n}\), training index set \(\mathcal{I}_{1}\), testing index set \(\mathcal{I}_{2}\), \(|\mathcal{I}_{1}|=\left\lfloor 0.7n\right\rfloor\), \(|\mathcal{I}_{2}|=\left\lfloor 0.3n\right\rfloor\), candidates for number of projection scores \(\left\{J_{0},\ldots,J_{N_{1}}\right\}\), candidates for depth \(\left\{L_{0},\ldots,L_{N_{2}}\right\}\), candidates for width \(\left\{p_{0},\ldots,p_{N_{3}}\right\}\), candidates for dropout rates \(\left\{s_{0},\ldots,s_{N_{4}}\right\}\), number of epochs, batch size, learning rate Output:\(\left\{\widehat{G}^{mHDNN}(X_{i})\right\}_{i=1}^{n}\)
1 Extract \(\widehat{\mathbf{\xi}}_{J}^{(i)}\) from \(\left\{X_{i}(s_{j})\right\}_{i,j=1}^{n,m}\) by integration, and \(J=\max\{J_{0},\ldots,J_{N_{1}}\}\)
2 For all candidates, (i) Train \(\widehat{\mathbf{f}}_{J\ell_{1},L_{\ell_{2}},p_{\ell_{3}},s_{\ell_{4}}}\) by Equation (2.4) with \(\left(\{X_{i}(s_{j})\}_{j=1}^{m},\mathbf{Y}_{i}\right)_{i\in\mathcal{I}_{1}}\). (ii) Compute \(\text{err}(\ell_{1},\ell_{2},\ell_{3},\ell_{4}):=\frac{1}{|\mathcal{I}_{2}|} \sum_{i\in\mathcal{I}_{2}}I(\widehat{\mathbf{f}}_{J\ell_{1},L_{\ell_{2}},p_{\ell_ {3}},s_{\ell_{4}}}\neq\mathbf{Y}_{i})\) Obtain \(\left(J_{\ell_{1}^{\star}},L_{\ell_{2}^{\star}},p_{\ell_{3}^{\star}},s_{\ell_ {4}^{\star}}\right)=\arg\min_{\ell_{1},\ell_{2},\ell_{3},\ell_{4}}\text{ err}(\ell_{1},\ell_{2},\ell_{3},\ell_{4})\) Train \(\widehat{\mathbf{f}}_{J\ell_{1}^{\star},L_{\ell_{2}^{\star}},p_{\ell_{3}^{\star}},s _{\ell_{4}^{\star}}}\) by Equation (2.4) with \(\left(\{X_{i}(s_{j})\}_{j=1}^{m},\mathbf{Y}_{i}\right)_{i=1}^{n}\)
```

**Algorithm 1**Training algorithm for mfDNN

## 3 Theoretical properties

### Functional data space

#### 3.1.1 Assumptions for functional multi-classification setup

To formulate our classification task, we first introduce two mild assumptions for \(\pi_{k}(\cdot)\). Without further notice, \(c\), \(C\), \(C_{1}\), \(C_{2}\),..., represent some positive constants and can vary from line to line.

**Assumption 1**.: _(Boundary condition) There exist an absolute constant \(C>0\) and some positive vector \(\mathbf{\alpha}=(\alpha_{1},\ldots,\alpha_{K})^{\intercal}\) such that_

\[\mathrm{P}\left(\pi_{k}(\mathbf{\xi})\leq x\right)\leq Cx^{\alpha_{k}},\qquad \forall x\in\left(0,1\right],\forall k\in\{1,\ldots,K\}. \tag{3.1}\]

Assumption 1 describes the behaviour of \(\pi_{k}\) around \(0\). Trivially, when \(\alpha_{k}=0\) for \(k=1,\ldots,K\), the condition holds universally only if \(C\geq 1\). Similar conditions can be found for multi-class classification for multivariate data in [2] and binary classification in [15] and [20].

**Assumption 2**.: _(Uniform lower bound) There exists a relatively small \(\epsilon\in(0,1)\), such that_

\[\mathrm{P}\left(\pi_{k}(\mathbf{\xi})>\epsilon\right)=1\]

_for all \(k=1,\ldots,K\)._
Assumption 2 provides a uniform lower bound of \(\pi_{k}\), indicating that \(\pi_{k}\) is bounded away from \(\epsilon\) in probability.

**Remark 1**.: _Assumptions 1 and 2 both depict the characteristic of \(\pi_{k}\) around zero in different aspects. Specifically, Assumption 1 controls the decay rate of the probability measure of \(\{\pi_{k}(\boldsymbol{\xi}):0\leq\pi_{k}(\boldsymbol{\xi})\leq x\}\). Assumption 2 truncates the probability measure below \(\epsilon\) to be zero. It can be seen that Assumptions 1 and 2 are closely related in certain circumstances. For example, given some \(\epsilon>0\), Assumption 2 implies Assumption 1 for arbitrary \(\boldsymbol{\alpha}\) and \(C=\max_{k}\epsilon^{-\alpha_{k}}\), and Assumption 1 implies Assumption 2 when \(C=0\) in a trivial case. However, they are not equivalent in most scenarios when \(C\), \(\{\boldsymbol{\alpha}_{k}\}_{k=1}^{K}\), and \(\epsilon\) are all universally positive constants._

Essentially, functional data have intrinsically infinite dimension. Owing to this unique phenomenon, we introduce the following finite approximation condition for \(\boldsymbol{\pi}\). For any positive integer \(J\), define \(\boldsymbol{\xi}^{(J)}=(\xi_{1},\ldots,\xi_{J})^{\intercal}\) and \(\boldsymbol{\pi}_{J}(\boldsymbol{x}_{J})=\left(\pi_{1}^{(J)}(\boldsymbol{x}_{ J}),\ldots,\pi_{K}^{(J)}(\boldsymbol{x}_{J})\right)^{\intercal}\), where \(\pi_{k}^{(J)}(\boldsymbol{x}_{J})=\mathrm{P}(\boldsymbol{Y}_{k}=\boldsymbol{ e}_{k}|\boldsymbol{\xi}_{J}=\boldsymbol{x}_{J})\), \(k=1,\ldots,K\).

**Assumption 3**.: _(Approximation error of \(\pi_{k}^{(J)}\)) There exist a constant \(J_{0}\geq 1\) and decreasing functions \(\zeta(\cdot):[1,\infty)\to\mathbb{R}_{+}\) and \(\Gamma(\cdot):[0,\infty)\to\mathbb{R}_{+}\), with \(\sup_{J\geq 1}J^{\varrho}\zeta(J)<\infty\) for some \(\varrho>0\) and \(\int_{0}^{\infty}\Gamma(x)dx<\infty\), such that for any \(J\geq J_{0}\), \(k=1\ldots,K\) and \(x>0\),_

\[\mathrm{P}\left(|\pi_{k}(\boldsymbol{\xi})-\pi_{k}^{(J)}(\boldsymbol{\xi}_{J}) |\geq x\right)\leq\zeta(J)\Gamma(x). \tag{3.2}\]

Assumption 3 provides a uniform upper bound on the probability of \(\pi_{k}\) differing from \(\pi_{k}^{(J)}\) by at least \(x\), which approaches zero if either \(J\) or \(x\) tends to infinity, implying that when dimension is reduced, there exists some large enough \(J\), such that \(\pi_{k}^{(J)}\) is an accurate approximation of \(\pi_{k}\). All three assumptions can be verified in concrete examples included in Section 4.

#### 3.1.2 Conditional probability \(\pi_{k}\) with complex structure

In the following, we impose composition structures on \(\pi_{k}(\cdot)\). For \(t\geq 1\), a measurable subset \(D\subset\mathbb{R}^{t}\), and constants \(\beta,R>0\), define

\[\mathcal{C}^{\beta}(D,R) = \left\{f:D\mapsto\mathbb{R}\big{|}\sum_{\boldsymbol{\alpha}:| \boldsymbol{\alpha}|<\beta}\|\partial^{\boldsymbol{\alpha}}f\|_{\infty}f\right.\] \[\left.+\sum_{\boldsymbol{\alpha}:|\boldsymbol{\alpha}|=\lfloor \beta\rfloor}\sup_{\boldsymbol{z},\boldsymbol{z}^{\prime}\in D,\boldsymbol{z }\neq\boldsymbol{z}^{\prime}}\frac{|\partial^{\boldsymbol{\alpha}}f( \boldsymbol{z})-\partial^{\alpha}f(\boldsymbol{z}^{\prime})|}{\|\boldsymbol{z }-\boldsymbol{z}^{\prime}\|_{\infty}^{\beta-\lfloor\beta\rfloor}}\leq R\right\},\]

where \(\partial^{\boldsymbol{\alpha}}\) = \(\partial^{\alpha_{1}}\ldots\partial^{\alpha_{t}}\) denotes the partial differential operator with multi-index \(\boldsymbol{\alpha}\) = \((\alpha_{1},\ldots,\alpha_{t})\in\mathbb{N}^{t}\), \(|\boldsymbol{\alpha}|=\alpha_{1}+\cdots+\alpha_{t}\). Equivalently, \(\mathcal{C}^{\beta}(D,C)\) is the ball of \(\beta\)-Holder smooth functions on \(D\) with radius \(R\). A function \(f:\mathbb{R}^{t}\to\mathbb{R}\) is said to be locally \(\beta\)-Holder smooth if for any \(a,b\in\mathbb{R}\), there exists a constant \(C\) (possibly depending on \(a,b\)) such that \(f\in\mathcal{C}^{\beta}([a,b]^{t},C)\).
For \(q\geq 0,J\geq 1\), let \(d_{0}=J\) and \(d_{q+1}=1\). For \(\mathbf{d}=(d_{1},\ldots,d_{q})\in\mathbb{N}_{+}^{q}\), \(\mathbf{t}=(t_{0},\ldots,t_{q})\in\mathbb{N}_{+}^{q+1}\) with \(t_{u}\leq d_{u}\) for \(u=0,\ldots,q\), \(\mathbf{\beta}:=(\beta_{0},\ldots,\beta_{q})\in\mathbb{R}_{+}^{q+1}\), let \(\mathcal{G}(q,J,\mathbf{d},\mathbf{t},\mathbf{\beta})\) be the class of functions \(g\) satisfying a modular expression

\[g(\mathbf{z})=g_{q}\circ\cdots\circ g_{0}(\mathbf{z})\in(0,1)\,,\ \forall\mathbf{z}\in \mathbb{R}^{d_{0}}, \tag{3.3}\]

where \(g_{u}=(g_{u1},\ldots,g_{ud_{u+1}}):\mathbb{R}^{d_{u}}\mapsto\mathbb{R}^{d_{u +1}}\) and \(g_{uv}:\mathbb{R}^{t_{u}}\mapsto\mathbb{R}\) are locally \(\beta_{u}\)-Holder smooth. The \(d_{u}\) arguments of \(g_{u}\) are locally connected in the sense that each component \(g_{uv}\) only relies on \(t_{u}(\leq d_{u})\) arguments. Similar structures have been considered by [18, 1, 13, 21, 2, 9, 8] in multivariate regression or classification to overcome high-dimensionality. Generalized additive model [7] and tensor product space ANOVA model [12] are special cases; see [13].

We define the class of \(\mathbf{h}=\{h_{1},\ldots,h_{K}\}\) as

\[\mathcal{H}\equiv \mathcal{H}\left(K,\{q^{(k)}\}_{k=1}^{K},\{\mathbf{d}^{(k)}\}_{k=1}^{ K},\{\mathbf{t}^{(k)}\}_{k=1}^{K},\{\mathbf{\beta}^{(k)}\}_{k=1}^{K},\{\alpha_{k}\}_{k= 1}^{K},\{\pi_{k}(\cdot)\}_{k=1}^{K},\right.\] \[\left.\zeta(\cdot),\Gamma(\cdot),C,\rho,\epsilon\right),\]

such that for any \(J\geq 1\), \(\pi_{k}^{(J)}\in\mathcal{G}\left(q^{(k)},J,\mathbf{d}^{(k)},\mathbf{t}^{(k)},\mathbf{ \beta}^{(k)}\right)\), where \(\mathbf{d}^{(k)}=(d_{1}^{(k)},\ldots,d_{q^{(k)}}^{(k)})\)\(\in\mathbb{N}_{+}^{q^{(k)}}\), \(\mathbf{t}^{(k)}=(t_{0}^{(k)},\ldots,t_{q^{(k)}}^{(k)})\in\mathbb{N}_{+}^{q^{(k)}+1}\) with \(t_{u}^{(k)}\leq d_{u}^{(k)}\) for \(u=0,\ldots,q^{(k)}\), \(\mathbf{\beta}^{(k)}:=(\beta_{0},\ldots,\beta_{q^{(k)}}^{(k)})\in\mathbb{R}_{+}^{q ^{(k)}+1}\). Note that this density class \(\mathcal{H}\) includes many popular models studied in literature, both Gaussian and non-Gaussian; see Section 4.

Throughout the paper, we explore \(\pi_{k}\) in some complicated \(\mathcal{G}\) with group-specific parameters \(q^{(k)},\mathbf{d}^{(k)},\mathbf{t}^{(k)}\) and \(\mathbf{\beta}^{(k)}\). The selection range of the truncation parameter \(J\) is based on the asymptotic order provided in Assumptions 4 and 5 in the next section. Although \(\pi_{k}^{(J)}\) has \(J\) arguments, it involves at most \(t_{0}^{(k)}d_{1}^{(k)}\) effective arguments, implying that the two population densities differ by a small number of variables. Relevant conditions are necessary for high-dimensional classification. For instance, in high-dimensional Gaussian data classification, [4, 3] show that, to consistently estimate Bayes classifier, it is necessary that the mean vectors differ at a small number of components. The modular structure holds for arbitrary \(J\), which may be viewed as a extension of [18] in the functional data analysis setting.

### Convergence of Kullback-Leibler divergence

#### 3.2.1 Kullback-Leibler divergence

For the true probability distribution \(\mathbf{\pi}(\mathbf{x})=(\pi_{1}(\mathbf{x}),\ldots,\pi_{K}(\mathbf{x}))^{\mathsf{T}}\) and any generic estimation \(\widehat{\mathbf{\pi}}(\mathbf{x})=(\widehat{\pi}_{1}(\mathbf{x}),\ldots,\widehat{\pi}_{K }(\mathbf{x}))^{\mathsf{T}}\), define the corresponding discrete version Kullback-Leibler divergence

\[\text{KL}\left(\mathbf{\pi}(\mathbf{x}),\widehat{\mathbf{\pi}}(\mathbf{x})\right)=\sum_{k=1}^{ K}\pi_{k}(\mathbf{x})\log\left(\frac{\pi_{k}(\mathbf{x})}{\widehat{\pi}_{k}(\mathbf{x})} \right).\]
For any \(\widehat{\pi}\) trained with \(\left\{(X_{i}(s),\mathbf{Y}_{i})\right\}_{i=1}^{n}\), we evaluate its performance by the log-likelihood ratio

\[\mathbb{E}\left[\sum_{k=1}^{K}Y_{k}\log\left(\frac{\pi_{k}(\mathbf{\xi})}{\widehat{ \pi}_{k}(\mathbf{\xi})}\right)\right]=\mathbb{E}\left[\mathit{KL}\left(\pi(\mathbf{\xi} ),\widehat{\pi}(\mathbf{\xi})\right)\right].\]

Note that the estimation risk is associate with the well-known CE loss in Section 2.2.

Different from the popular least square loss, the CE loss is the expectation with respect to the input distribution of the Kullback-Leibler divergence of the conditional class probabilities. If anyone of the conditional class probabilities has zero estimation while the underlying conditional class probability is positive, the risk can even become infinite. To avoid the infinite risk, we truncate the CE loss function and derive convergence rates without assuming either the true conditional class probabilities or the estimators away from zero. Instead, our misclassification risks depend on an index quantifying the behaviour of the conditional class probabilities near zero.

Given an absolute constant \(C_{0}\geq 2\), for any classifier \(\widehat{\pi}\), define the truncated Kullback-Leibler risk for some density \(\mathbf{h}\) as

\[R_{\mathbf{h},C_{0}}\left(\widehat{\pi}\right)=\mathrm{E}_{\mathbf{h}}\left\{\sum_{k =1}^{K}\pi_{k}(\mathbf{\xi})\left[C_{0}\wedge\log\left(\frac{\pi_{k}(\mathbf{\xi})}{ \widehat{\pi}_{k}(\mathbf{\xi})}\right)\right]\right\}. \tag{3.4}\]

**Remark 2**.: \(C_{0}\) _is commonly introduced to avoid infinity value of the ordinary Kullback-Leibler risk. It is trivial to show that \(C_{0}\) can only be abandoned when \(\widehat{\pi}_{k}\) is lower bounded by \(\exp\left(-C_{0}\right)\) for all \(k\). When \(\widehat{\pi}_{k}\) has a large deviation from \(\pi_{k}\) for some \(k\), the Kullback-Leibler risk explodes to infinity, which makes it is infeasible to evaluate the performance of \(\widehat{\pi}\)._

#### 3.2.2 Convergence rate for fully observed functional data

In this section, we provide the non-asymptotic Kullback-Leibler risk of mfDNN classifier. Let

\[(\widehat{k},\widehat{u})=\arg\max_{\begin{subarray}{c}(k,u),k=1\ldots,K,\\ u=0,\ldots,q^{(k)}\end{subarray}}\frac{t_{u}^{(k)}}{\widehat{\beta}_{u}^{(k )}}, \tag{3.5}\]

where \(\widehat{\beta}_{u}^{(k)}=\beta_{u}^{(k)}\prod_{\ell=u+1}^{q^{(k)}}\beta_{ \ell}^{(k)}\wedge 1\), and

\[(\widehat{l},\widehat{v})=\arg\max_{\begin{subarray}{c}(k,u),k=1\ldots,K,\\ u=0,\ldots,q^{(k)}\end{subarray}}\frac{t_{u}^{(k)}}{\widehat{\beta}_{u}^{(k )}(1+\alpha_{k})}.\]

Let \(\theta=\frac{(1+\alpha_{\widehat{l}})\widehat{\beta}_{\widehat{v}}^{(\widehat {l})}}{(1+\alpha_{\widehat{l}})\widehat{\beta}_{\widehat{v}}^{(l)}+t_{\widehat {v}}^{(l)}}\) and \(\nu=\frac{\theta t_{\widehat{v}}^{(\widehat{k})}}{\widehat{\beta}_{\widehat{u} }^{(k)}(1+\widehat{\alpha})}\), where \(\widetilde{\alpha}=\min_{k}\alpha_{k}\wedge 1\).

**Assumption 4**.: _There exist some constants \(C_{1},C_{2}^{\prime},C_{2},C_{3}\), only depending on \(\mathcal{H}\) and \(C_{0}\), such that the DNN class \(\mathcal{F}(L,J,\mathbf{p},s)\) satisfies_
1. \(L\leq C_{1}\log n\)_;_
2. \(C_{2}^{\prime}n^{\theta/\rho}\leq J\leq C_{2}n^{\nu}\)_;_
3. \(\max_{1\leq\ell\leq L}p_{\ell}\leq C_{2}n^{\nu}\)_;_
4. \(s\leq C_{3}n^{\nu}\log n\)_._

Assumption 4 provides exact orders on \((L,\mathbf{p},s)\) for network, respectively. Assumption 4 provides the precise range on \(J\). It is worth mentioning this condition implies \(\varrho\geq\theta\nu^{-1}\), i.e., the function \(\zeta(J)\) converges to zero in a relatively fast rate when \(J\to\infty\).

In the following, we provide the convergence rate in the ideal case when the entire functional curve is fully observed.

**Theorem 3.1**.: _There exist a positive constant \(\omega_{1}\), only depending on \(\mathcal{H}\) and \(C_{0}\), such that_

\[\sup_{\mathbf{h}\in\mathcal{H}}R_{\mathbf{h},C_{0}}\left(\widehat{\mathbf{\pi}}\right) \leq\omega_{1}n^{-\theta}\log^{3}n,\]

_where network classifier \(\widehat{\mathbf{\pi}}\) belongs to \(\mathcal{F}(L,J,\mathbf{p},s)\) in Assumption 4._

**Remark 3**.: _Theorem 3.1 provides the upper bound of the KL misclassification risk of the proposed mfDNN. When \(K=2\), the multiclass classification downgrades to the binary classification problem. Compared with the minimax excess misclassification risk derived in [22] for the binary classification, the upper bound rate in Theorem 3.1 provides a slightly larger order, Specifically, the leading term of the upper bound rate in Theorem 3.1 is \(n^{-\theta}\) with \(\theta=\frac{(1+\alpha)\beta}{(1+\alpha)\beta+t}\), while the leading terms in Theorem 1 in [22] is \(n^{-S_{0}}\) with \(S_{0}=\frac{(\alpha+1)\beta}{(\alpha+2)\beta+t}\). Hence, when \(\beta\)'s, i.e. the degrees of Holder smooth functions in the class \(\mathcal{H}\), are large enough, the discrepancy between these two bounds are rarely negligible. Hence, to reduce potential slightly larger risks, we recommend [22] for the binary classification problems and when the distribution functions are not smooth enough. Meanwhile, we recommend the mfDNN classifier for multiclassification problems regardless the smoothness of the conditional distribution functions._

#### 3.2.3 Convergence rate for discretely observed functional data

Practically, it is usually unrealistic to observe the full trajectory of each individual, thus the rate in Theorem 3.1 can only be reached if sampling frequency is dense enough. Hence, it is interesting to discuss the upper bound of the risk of mfDNN classifier when functional data are discretely observed at \(m\) occasions for each subject. Let \(\widetilde{\beta}=\max_{k=1\ldots,K}\left(\widetilde{\beta}_{0}^{(k)}\wedge 1\right)\), \(\theta^{\prime}=\tau\widetilde{\beta}\), and \(\nu^{\prime}=\frac{\theta^{\prime}t_{\widehat{\sigma}}^{(\widehat{k})}}{ \widetilde{\beta}_{\widehat{u}}^{(\widehat{k})}(1+\widetilde{\alpha})}\), where \(\widehat{k}\) is defined in (3.5) and \(\tau\) is a positive universal constant.

**Assumption 5**.: _There exist some constants \(C_{1},C_{2}^{\prime},C_{2},C_{3}\), \(\widetilde{C}_{1},\widetilde{C}_{2}^{\prime},\widetilde{C}_{2},\) and \(\widetilde{C}_{3}\) only depending on \(\mathcal{H}\), \(C_{0}\) and \(\tau\), and a phase transition point \(m^{*}\in\mathbb{N}^{+}\), such that the DNN class \(\mathcal{F}(L,J,\mathbf{p},s)\) satisfies_

1. \(L\leq C_{1}\log n\mathbb{I}(m\geq m^{*})+\widetilde{C}_{1}\log m\mathbb{I}(m<m ^{*})\)
_;_
2. \(C_{2}^{\prime}n^{\theta/\rho}\mathbb{I}(m\geq m^{*})+\widetilde{C}_{2}^{\prime}m^{ \theta^{\prime}/\rho}\mathbb{I}(m<m^{*})\leq J\leq C_{2}n^{\nu}\mathbb{I}(m\geq m ^{*})+\widetilde{C}_{2}m^{\nu^{\prime}}\mathbb{I}(m<m^{*})\)_;_
3. \(\max_{1\leq\ell\leq L}p_{\ell}\leq C_{2}n^{\nu}\mathbb{I}(m\geq m^{*})+ \widetilde{C}_{2}m^{\nu^{\prime}}\mathbb{I}(m<m^{*})\)_;_
4. \(s\leq C_{3}n^{\nu}\log n\mathbb{I}(m\geq m^{*})+\widetilde{C}_{3}m^{\nu^{ \prime}}\log m\mathbb{I}(m<m^{*})\)_._

Similar to Assumption 4, Assumption 5 provides exact orders on \(L,\mathbf{p},s\), and range of \(J\) when sampling frequency \(m\) is involved. When \(m\geq m^{*}\), Assumption 5 coincides with Assumption 4 for dense functional data.

The following theorem provides the phase transition rate when functional data are discretely observed on \(m\) locations at a certain rate with respect to \(\tau\).

**Theorem 3.2**.: _When \(E|\xi_{j}-\widehat{\xi}_{j}|\lesssim m^{-\tau}\) for all \(j=1,\ldots,J\), there exists positive constants \(\omega_{1},\omega_{2}\) and \(\omega_{3}\) only depending on \(\mathcal{H}\), \(C_{0}\) and \(\tau\), such that_

\[\sup_{\mathbf{h}\in\mathcal{H}}R_{\mathbf{h},C_{0}}\left(\widehat{\mathbf{\pi}}\right) \leq\omega_{1}n^{-\theta}\log^{3}n\mathbb{I}(m\geq m^{*})+\omega_{2}m^{-\theta^ {\prime}}\mathbb{I}(m<m^{*}),\]

_where \(m^{*}=\lfloor\left(\omega_{3}n^{\theta}/\log^{3}n\right)^{1/\theta^{\prime}}\rfloor\), and \(\widehat{\mathbf{\pi}}\in\mathcal{F}(L,J,\mathbf{p},s)\) defined in Assumption 5._

**Remark 4**.: _When functional curves are discretely observed, Theorem 3.2 provides the convergence rate of truncated KL risk when the biases of projection scores are uniformly bounded by \(m^{-\tau}\). This assumption does not hold universally for all scenarios. Nevertheless, it can be satisfied in various examples._

_For any empirical process on \([0,1]\), if we use Fourier basis with \(m\) terms to decompose the curve, we can show that \(E|\xi_{j}-\widehat{\xi}_{j}|\leq O\left(\max_{k=1,\ldots,K}\sqrt{\sum_{j=m+1}^ {\infty}\left(\lambda_{kj}+\mu_{kj}^{2}\right)}\right)\), where \(\mu_{kj}=\operatorname{E}\xi_{j}\) in the \(k\)-th group. As \(\left\{\lambda_{kj}\right\}_{j=1}^{\infty}\) and \(\left\{\mu_{kj}^{2}\right\}_{j=1}^{\infty}\) are both convergent, when \(\lambda_{kj}\) and \(\mu_{kj}^{2}\) are decreasing no faster than some polynomial order, the assumption easily follows._

_Another well-known example is the FPCA provided by [6]. According to Theorem 1 in [6], for all \(j=1,\ldots,J\), the estimators of projection scores satisfy \(\operatorname{E}|\xi_{j}-\widehat{\xi}_{j}|\leq\max_{k}\left(\int_{[0,1]^{d} \times[0,1]^{d}}\left(\Omega_{k}-\widehat{\Omega}_{k}\right)^{2}(s,s^{\prime}) dsds^{\prime}\right)^{1/2}\). Hence, \(E|\xi_{j}-\widehat{\xi}_{j}|\lesssim m^{-\tau}\) holds easily when \(\max_{k}\left(\int_{[0,1]^{d}\times[0,1]^{d}}\left(\Omega_{k}-\widehat{\Omega} _{k}\right)^{2}(s,s^{\prime})dsds^{\prime}\right)\) is bounded properly._

## 4 Examples

In this section, we provide two examples of exponential families to justify the validation of our model assumptions, and emphasize the necessity of applying DNN approach owing to the complicated structure of data population. For simplicity, we assume the prior probability satisfies \(\mathbb{P}(\mathbf{Y}=\mathbf{e}_{k})=k^{-1}\) for all \(k\) throughout the section.

### Independent exponential family

We first consider independent projection scores, which are from exponential families. For some collection of unknown parameters \(\left\{\theta_{kj}\right\}_{k=1,j=1}^{K,\infty}\), and unknown collections of
functions \(\left\{\eta_{kj}\right\}_{k=1,j=1}^{K,\infty}\), \(\left\{U_{kj}\right\}_{k=1,j=1}^{K,\infty}\), and \(\left\{W_{kj}\right\}_{k=1,j=1}^{K,\infty}\), we consider the \(k\)-th class conditional density, such that

\[\boldsymbol{\xi}|\left\{\theta_{kj}\right\}_{j=1}^{\infty},\boldsymbol{Y}= \boldsymbol{e}_{k}\sim h_{k}(\boldsymbol{x})=\exp\left(\sum_{j=1}^{\infty} \eta_{kj}(\theta_{kj})U_{kj}(x_{j})+W_{kj}(x_{j})\right).\]

For all \(1\leq k,k^{\prime}\leq K\), let \(A_{kk^{\prime}}=\left\{j:\eta_{kj}(\theta_{kj})U_{kj}(x_{j})\neq\eta_{k^{\prime }j}(\theta_{k^{\prime}j})U_{k^{\prime}j}(x_{j}),\forall x_{j}\in\mathbb{R}\right\}\) and \(B_{kk^{\prime}}=\left\{j:W_{kj}\neq W_{k^{\prime}j}\right\}\) be the two sets identifying the difference between \(h_{k}\) and \(h_{k^{\prime}}\). Therefore, we have the pairwise log likelihood

\[\log\left(h_{k}/h_{k^{\prime}}\right) = \sum_{j\in(A_{kk^{\prime}}\bigcup B_{kk^{\prime}})}\left\{[\eta_{ kj}(\theta_{kj})U_{kj}(x_{j})-\eta_{k^{\prime}j}(\theta_{k^{\prime}j})U_{k^{ \prime}j}(x_{j})]\right.\] \[\left.+\left[W_{kj}(x_{j})-W_{k^{\prime}j}(x_{j})\right]\right\}.\]

Given some universal constant \(N_{kk^{\prime}}\), when \(|A_{kk^{\prime}}\bigcup B_{kk^{\prime}}|\leq N_{kk^{\prime}}\), there exists a positive integer \(J_{max}=\max\bigcup_{k,k^{\prime}}\left(A_{kk^{\prime}}\bigcup B_{kk^{\prime}}\right)\), such that \(\pi_{k}=\pi_{k}^{(J)}\) for all \(J\geq J_{max}\).

By definition, Assumption 1 holds for \(\alpha=1\). Assumption 2 holds when \(h_{k}/h_{k^{\prime}}\) is bounded for all pairs. Note that it is trivial when \(\{h_{k}\}_{k=1}^{K}\) share the same \(\left\{U_{kj}\right\}_{k=1,j=1}^{K,\infty}\) and \(\left\{W_{kj}\right\}_{k=1,j=1}^{K,\infty}\), such as Gaussian distribution, student's t distribution and exponential distribution, whose density ratio of their kind is always bounded. Assumption 3 holds for \(J_{0}=J_{max}\), and for arbitrary function \(e(\cdot)\) with exponential tails and density \(\pi(\cdot)\). Since \(\pi_{k}=\left[1+\sum_{k^{\prime}\neq k}\log(h_{k}/h_{k^{\prime}})\right]^{-1}\), the smoothness is determined by \(\left\{\eta_{kj}\right\}_{k=1,j=1}^{K,\infty}\), \(\left\{U_{kj}\right\}_{k=1,j=1}^{K,\infty}\), and \(\left\{W_{kj}\right\}_{k=1,j=1}^{K,\infty}\), thus \(\boldsymbol{h}\) is trivially in some \(\mathcal{H}\).

### Exponential family with in-block interaction

In this example, we consider \(\xi_{j}\) are dependent with each other in a block, but independent across blocks, which is an extension of the example in Section 4.1. Given a sequence of positive integers \(\left\{\ell_{p}\right\}_{j=1}^{\infty}\), such that \(0=\ell_{1}<\ell_{2}<\ldots\), we define the \(p\)-th group index set \(\mathcal{E}_{p}=\left\{\ell_{p}+1,\ldots,\ell_{p+1}\right\}\), such that the cardinality \(|\mathcal{E}_{p}|=\ell_{p+1}-\ell_{p}\), therein grouping are based on adjacent members for simplicity. For a collection of unknown parameters \(\left\{\boldsymbol{\theta}_{kj}\right\}_{k=1,p=1}^{K,\infty}\), and unknown collections of functions \(\left\{\widetilde{\eta}_{kp}\right\}_{k=1,p=1}^{K,\infty}\), \(\left\{\widetilde{U}_{kp}\right\}_{k=1,p=1}^{K,\infty}\), and \(\left\{\widetilde{W}_{kp}\right\}_{k=1,p=1}^{K,\infty}\), where \(\widetilde{U}_{kp}\) and \(\widetilde{W}_{kp}\) are functions from \(\mathbb{R}^{|\mathcal{E}_{p}|}\) to \(\mathbb{R}\).

Consider the joint conditional density

\[h_{k}(\boldsymbol{x})=\exp\left(\sum_{p=1}^{\infty}\eta_{kp}(\boldsymbol{ \theta}_{kp})\widetilde{U}_{kp}(\boldsymbol{x}_{p})+\widetilde{W}_{kp}( \boldsymbol{x}_{p})\right)\]

for class \(k\), where \(\boldsymbol{x}_{p}=(x_{j})_{j=\ell_{p}+1,\ldots,\ell_{p+1}}\).

For any \(1\leq k,k^{\prime}\leq K\), define the density difference sets

\[\widetilde{A}_{kk^{\prime}}=\left\{p:\widetilde{\eta}_{kp}(\boldsymbol{\theta }_{kp})\widetilde{U}_{kp}(\boldsymbol{x}_{p})\neq\eta_{k^{\prime}p}(\boldsymbol{ \theta}_{k^{\prime}p})\widetilde{U}_{k^{\prime}p}(\boldsymbol{x}_{p}),\forall \boldsymbol{x}_{p}\in\mathbb{R}^{|\mathcal{E}_{p}|}\right\}\]
and \(\widetilde{B}_{kk^{\prime}}=\left\{p:\widetilde{W}_{kp}\neq\widetilde{W}_{k^{ \prime}p}\right\}\), and the pairwise log likelihood is thus given by

\[\log\left(h_{k}/h_{k^{\prime}}\right)\] \[= \sum_{p\in\widetilde{A}_{kk^{\prime}}\bigcup\widetilde{B}_{kk^{ \prime}}}\left\{\left[\widetilde{\eta}_{kp}(\boldsymbol{\theta}_{kp}) \widetilde{U}_{kp}(\boldsymbol{x})-\widetilde{\eta}_{k^{\prime}p}(\boldsymbol {\theta}_{k^{\prime}p})\widetilde{U}_{k^{\prime}p}(\boldsymbol{x})\right]+ \left[\widetilde{W}_{kp}(\boldsymbol{x})-\widetilde{W}_{k^{\prime}p}( \boldsymbol{x})\right]\right\}.\]

Given some finite positive number \(N_{kk^{\prime}}\), such that \(|A_{kk^{\prime}}\bigcup B_{kk^{\prime}}|\leq N_{kk^{\prime}}\), the verification can be similarly derived from Section 4.1.

## 5 Simulation studies

In this section, we provide numerical evidences to demonstrate the superior performance of mfDNN. In all simulations, we generated \(n_{k}=200,350,700\) training samples for each class, and testing sample sizes \(100,150,300\), respectively. Based on the fact that there is no existing multi-classification method specifically designed for multidimensional functional data, for comparison, we include the sparse discriminate analysis and \(\ell_{1}\) penalized Fisher's discriminant analysis (MSDA) approach introduced in [14] and penalized linear discriminant analysis (PLDA) classifier in [24]. In fact, MSDA and PLDA are efficient classifiers designed for high-dimensional i.i.d. observations. To make these two methods directly applicable to functional data, we first pre-processed 2D or 3D functional data by vectorization. The realization is via the R packages msda and PenalizedLDA, where the tuning parameter candidates for the penalty term are generated by default. We use the default five-fold and six-fold cross-validation to tune MSDA and PLDA, respectively. For mfDNN, we use tensor of Fourier basis to extract projection scores by integration. The structure parameters \((L,J,\boldsymbol{p},s)\) are selected by Algorithm 1, where the candidates are given based on Theorem 3.2. We summarize R codes and examples for the proposed mfDNN algorithms on GitHub ([https://github.com/FDASTATAUBURN/mfdnn](https://github.com/FDASTATAUBURN/mfdnn)).

### 2D functional data

For \(k=1,2,3\), we generated functional data \(X_{i}^{(k)}(s,s^{\prime})=\sum_{j=1}^{5}\xi_{ij}^{(k)}\psi_{j}(s,s^{\prime}),s,s^{\prime}\in[0,1]\). Let \(\psi_{1}(s,s^{\prime})=s\), \(\psi_{2}(s,s^{\prime})=s^{\prime}\), \(\psi_{3}(s,s^{\prime})=ss^{\prime}\), \(\psi_{4}(s,s^{\prime})=s^{2}s^{\prime}\), \(\psi_{5}(s,s^{\prime})=ss^{\prime}^{2}\). Define \(\mathbf{1}_{k}\) be a \(k\times 1\) vector with all the elements one. We specify the distribution of \(\xi_{ij}^{(k)}\)'s as following.

_Model 1_ (2D Gaussian): Let \(\left(\xi_{i1}^{(k)},\ldots,\xi_{i5}^{(k)}\right)^{\intercal}\sim N( \boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k})\), where \(\boldsymbol{\mu}_{1}=(4,4,3,3,3)^{\intercal}\), \(\boldsymbol{\Sigma}_{1}^{1/2}=\text{diag}\left(8,7,6,5,4\right)\), \(\boldsymbol{\mu}_{2}=-\mathbf{1}_{5}\), \(\boldsymbol{\Sigma}_{2}^{1/2}=\text{diag}\left(5,4,3,2,1\right)\), \(\boldsymbol{\mu}_{3}=\mathbf{0}_{5}\), \(\boldsymbol{\Sigma}_{3}^{1/2}=\text{diag}\left(2.5,2,1.5,1,0.5\right)\).

_Model 2_ (2D Mixed 1): Let \(\left(\xi_{i1}^{(k)},\ldots,\xi_{i5}^{(k)}\right)^{\intercal}\sim N( \boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k})\) for \(k=1,2\), and \(\xi_{ij}^{(3)}\sim t_{2j+1}(\nu_{j})\), where \(\boldsymbol{\mu}_{1}=-\mathbf{1}_{5}\), \(\boldsymbol{\Sigma}_{1}^{1/2}=\text{diag}\left(5,4,3,2,1\right)\), \(\boldsymbol{\mu}_{2}=\mathbf{0}_{5}\), \(\boldsymbol{\Sigma}_{2}^{1/2}=\text{diag}\left(\frac{5}{2},2,\frac{3}{2},1, \frac{1}{2}\right)\), \((\nu_{1},\ldots,\nu_{5})^{\intercal}=3\cdot\mathbf{1}_{5}\).
_Model 3_ (2D Mixed 2): Let \(\left(\xi_{i1}^{(1)},\ldots,\xi_{i5}^{(1)}\right)^{\intercal}\sim N(\mathbf{\mu}_{k}, \mathbf{\Sigma}_{k})\), \(\xi_{ij}^{(2)}\sim t_{j+1}(\nu_{2j})\), \(\xi_{ij}^{(2)}\sim t_{2j+1}(\nu_{3j})\), where \(\mathbf{\mu}_{1}=\mathbf{0}_{5}\), \(\mathbf{\Sigma}_{1}^{1/2}=\text{diag}\left(\frac{5}{2},2,\frac{3}{2},1,\frac{1}{2}\right)\), \((\nu_{21},\ldots,\nu_{25})^{\intercal}=\mathbf{1}_{5}\), \((\nu_{31},\ldots,\nu_{35})^{\intercal}=3\cdot\mathbf{1}_{5}\).

_Model 4_ (2D Mixed 3): Let \(\xi_{ij}^{(1)}\sim\text{{Exp}}(r_{1j})\), \(\xi_{ij}^{(2)}\sim t_{2j+1}(\nu_{2j})\), and \(\left(\xi_{i1}^{(3)},\ldots,\xi_{i5}^{(3)}\right)^{\intercal}\)\(\sim N(\mathbf{\mu}_{3},\mathbf{\Sigma}_{3})\), where \((r_{11},\ldots,r_{15})^{\intercal}=(0.1,0.3,0.5,0.7,0.9)^{\intercal}\), \((\nu_{21},\ldots,\nu_{25})^{\intercal}=3\times\mathbf{1}_{5}\), \(\mathbf{\mu}_{3}=\mathbf{0}_{5}\), \(\mathbf{\Sigma}_{3}^{1/2}=\text{diag}\left(2.5,2,1.5,1,0.5\right)\).

For each model, we observe the functional data on \(3\times 3\), \(5\times 5\), \(10\times 10\), and \(20\times 20\) grid points over \(\left[0,1\right]^{2}\), respectively. As a result, the sampling frequency \(m=9,25,100,400\), which indicates that the functional observations are from sparse to dense. Tables 1 and 2 demonstrate the results of \(100\) simulations. For mfDNN, it can be seen that the misclassification risks decrease as the sample size \(n\) increasing, as well as the increase of the sampling frequency \(m\). This founding further confirms Theorem 3.2. Given the relatively sparse sampling frequency, i.e., \(m=9\), MSDA has slightly better performance than mfDNN does. However, despite the increase of \(m\) in Table 2, there is no improvement of both MSDA and PLDA methods in terms of misclassification risks. This finding indicates that MSDA and PLDA classifiers can not be improved with more gathered information. In summary, the simulation results illustrate that the proposed mfDNN method outperforms the existing sparse and penalized discriminant analysis when classifying 2D dense functional data.

### 3D functional data

For \(k=1,2,3\) and \(s_{1},s_{2},s_{3}\in[0,1]\), we generate 3D functional data \(X_{i}^{(k)}(s_{1},s_{2},s_{3})=\sum_{j=1}^{9}\xi_{ij}^{(k)}\psi_{j}(s_{1},s_{ 2},s_{3})\) where \(\psi_{1}(s_{1},s_{2},s_{3})=s_{1}\), \(\psi_{2}(s_{1},s_{2},s_{3})=s_{2}\), \(\psi_{3}(s_{1},s_{2},s_{3})=s_{3},\psi_{4}(s_{1},s_{2},s_{3})=s_{1}s_{2}\), \(\psi_{5}(s_{1},s_{2},s_{3})=s_{1}s_{3}\), \(\psi_{6}(s_{1},s_{2},s_{3})=s_{2}s_{3}\), \(\psi_{7}(s_{1},s_{2},s_{3})=s_{1}^{2}\), \(\psi_{8}(s_{1},s_{2},s_{3})=s_{2}^{2}\), \(\psi_{9}(s_{1},s_{2},s_{3})=s_{3}^{2}\), and the distribution of \(\xi_{ij}^{(k)}\)'s are specified as below.

_Model 5_ (3D Gaussian): Let \(\left(\xi_{i1}^{(k)},\ldots,\xi_{i9}^{(k)}\right)^{\intercal}\sim N(\mathbf{\mu}_{ k},\mathbf{\Sigma}_{k})\), where \(\mathbf{\mu}_{1}=2\times\mathbf{1}_{9}\), \(\mathbf{\Sigma}_{1}^{1/2}=\mathbf{\Sigma}_{2}^{1/2}=\text{diag}\left(9,8,7,6,5,4,3,2,1\right)\), \(\mathbf{\mu}_{2}=\mathbf{\mu}_{3}=\mathbf{0}_{9}\), \(\mathbf{\Sigma}_{3}^{1/2}=1/3\times\mathbf{\Sigma}_{1}^{1/2}\).

_Model 6_ (3D Mixed 1): Let \(\left(\xi_{i1}^{(k)},\ldots,\xi_{i9}^{(k)}\right)^{\intercal}\sim N(\mathbf{\mu}_{ k},\mathbf{\Sigma}_{k})\) for \(k=1,2\), and \(\xi_{ij}^{(3)}\sim t_{j+1}(\nu_{j})\), where \((\nu_{1},\ldots,\nu_{5})^{\intercal}=3\times\mathbf{1}_{5}\), \(\mathbf{\mu}_{1}=-\mathbf{1}_{9}\), \(\mathbf{\Sigma}_{1}^{1/2}=\text{diag}(5.5,5,4.5,4,\)\(3.5,3.2,5.2,1.5)\), \(\mathbf{\mu}_{2}=\mathbf{0}_{9}\), \(\mathbf{\Sigma}_{2}^{1/2}=\text{diag}\left(4.5,4,3.5,3,2.5,2,1.5,1,0.5\right)\).

_Model 7_ (3D Mixed 2): Let \(\left(\xi_{i1}^{(1)},\ldots,\xi_{i9}^{(1)}\right)^{\intercal}\sim N(\mathbf{\mu}_{ k},\mathbf{\Sigma}_{k})\), \(\xi_{ij}^{(2)}\sim t_{j+1}(\nu_{2j})\), where \(\mathbf{\mu}_{1}=\mathbf{0}_{9}\), \(\mathbf{\Sigma}_{1}^{1/2}=\text{diag}\left(4.5,4,3.5,3,2.5,2,1.5,1,0.5\right)\), \((\nu_{21},\ldots,\nu_{29})^{\intercal}=-\mathbf{1}_{9}\), \((\nu_{31},\ldots,\nu_{39})^{\intercal}=0.5\times\mathbf{1}_{9}\).

_Model 8_ (3D Mixed 3): Let \(\xi_{ij}^{(1)}\sim\text{{Exp}}(r_{1j})\), \(\xi_{ij}^{(2)}\sim t_{j+1}(\nu_{2j})\), and \(\left(\xi_{i1}^{(3)},\ldots,\xi_{i9}^{(3)}\right)^{\intercal}\)\(\sim N(\mathbf{\mu}_{3},\mathbf{\Sigma}_{3})\), where \((r_{11},\ldots,r_{19})^{\intercal}=0.1\times(1,3,5,7,9,11,13,15,17)^{\intercal}\), \((\nu_{21},\ldots,\nu_{29})^{\intercal}=0.6\times\mathbf{1}_{9}\), \(\mathbf{\mu}_{3}=\mathbf{0}_{9}\), \(\mathbf{\Sigma}_{3}^{1/2}=\text{diag}\left(4.5,4,3.5,3,2.5,2,1.5,1,0.5\right)\).

For the 3D functional data, we apply similar setups as 2D cases. We observe the functional data on \(2\times 2\times 2\), \(3\times 3\times 3\), \(4\times 4\times 4\), and \(5\times 5\times 5\) grid points over
\(\left[0,1\right]^{3}\), respectively, and the sampling frequency \(m=8,27,64,125\). Tables 3 and 4 demonstrate the results of \(100\) simulations. The proposed mfDNN classifier is superior to its counterparts for all 3D functional data cases. Meanwhile, there also exists the phase transition patterns for mfDNN method. However, the performance of MSDA and PLDA methods lacks of improvement with the increase of \(m\). It can be seen that when \(m=125\), the misclassification error rates of mfDNN are almost one third of MSDA's and one forth of PLDA's in Gaussian case, and almost a half of either MSDA's or PLDA's error rates in Models 7 and 8. A plausible reason is that given the functional data framework, our proposed mfDNN can properly accommodate the repeatedly observed data over pixels or voxels, while other competitors only treat those information as common high-dimensional covariates and ignore the underlining smoothing structures. By efficiently extracting the projection scores of the continuum, the proposed mfDNN has full potential to discover the underlying distributions of the functional data clusters. This again demonstrates our proposed classifier has a distinct advantage over these competitors in complex imaging data classification problems.

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \multirow{2}{*}{Model} & \multirow{2}{*}{\(n_{k}\)} & \multicolumn{4}{c}{\(m=9\)} & \multicolumn{4}{c}{\(m=25\)} \\ \cline{3-8}  & & mfDNN & MSDA & PLDA & mfDNN & MSDA & PLDA \\ \hline \multirow{4}{*}{2D Gaussian} & 200 & 0.259 & **0.246** & 0.323 & **0.205** & 0.245 & 0.325 \\  & & (0.031) & (0.027) & (0.028) & (0.026) & (0.028) & (0.028) \\ \cline{2-8}  & 350 & **0.232** & 0.243 & 0.325 & **0.200** & 0.243 & 0.328 \\  & & (0.020) & (0.022) & (0.022) & (0.020) & (0.021) & (0.022) \\ \cline{2-8}  & 700 & **0.227** & 0.241 & 0.323 & **0.196** & 0.241 & 0.325 \\  & & (0.014) & (0.014) & (0.016) & (0.014) & (0.015) & (0.016) \\ \hline \multirow{4}{*}{2D Mixed 1} & 200 & 0.158 & **0.150** & 0.227 & **0.100** & 0.150 & 0.229 \\  & & (0.020) & (0.021) & (0.026) & (0.018) & (0.022) & (0.025) \\ \cline{2-8}  & 350 & 0.153 & **0.144** & 0.227 & **0.085** & 0.144 & 0.229 \\  & & (0.015) & (0.016) & (0.020) & (0.012) & (0.016) & (0.019) \\ \cline{2-8}  & 700 & **0.152** & 0.156 & 0.229 & **0.081** & 0.123 & 0.231 \\  & & (0.011) & (0.014) & (0.014) & (0.010) & (0.014) & (0.014) \\ \hline \multirow{4}{*}{2D Mixed 2} & 200 & 0.166 & **0.152** & 0.198 & **0.140** & 0.153 & 0.200 \\  & & (0.023) & (0.022) & (0.025) & (0.020) & (0.020) & (0.025) \\ \cline{2-8}  & 350 & 0.165 & **0.152** & 0.200 & **0.135** & 0.152 & 0.202 \\  & & (0.016) & (0.016) & (0.022) & (0.016) & (0.016) & (0.022) \\ \cline{2-8}  & 700 & 0.163 & **0.148** & 0.197 & **0.123** & 0.148 & 0.199 \\  & & (0.011) & (0.012) & (0.013) & (0.010) & (0.012) & (0.013) \\ \hline \multirow{4}{*}{2D Mixed 3} & 200 & 0.168 & **0.115** & 0.264 & 0.136 & **0.116** & 0.264 \\  & & (0.030) & (0.018) & (0.023) & (0.022) & (0.018) & (0.024) \\ \cline{1-1} \cline{2-8}  & 350 & 0.164 & **0.115** & 0.265 & 0.132 & **0.115** & 0.265 \\ \cline{1-1}  & & (0.024) & (0.014) & (0.021) & (0.019) & (0.014) & (0.021) \\ \cline{1-1} \cline{2-8}  & 700 & 0.156 & **0.114** & 0.264 & 0.123 & **0.114** & 0.264 \\ \cline{1-1} \cline{2-8}  & & (0.014) & (0.011) & (0.016) & (0.014) & (0.011) & (0.016) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Averaged misclassification rates with standard errors in brackets for 2D simulations when \(m=9\) and \(m=25\) over \(100\) replicates.
\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \multirow{2}{*}{Model} & \multirow{2}{*}{\(n_{k}\)} & \multicolumn{3}{c}{\(m=100\)} & \multicolumn{3}{c}{\(m=400\)} \\ \cline{3-8}  & & mfDNN & MSDA & PLDA & mfDNN & MSDA & PLDA \\ \hline \multirow{4}{*}{2D Gaussian} & 200 & **0.147** & 0.244 & 0.326 & **0.145** & 0.245 & 0.329 \\  & & (0.024) & (0.027) & (0.029) & (0.025) & (0.027) & (0.029) \\ \cline{2-8}  & \multirow{2}{*}{350} & **0.139** & 0.242 & 0.328 & **0.139** & 0.242 & 0.331 \\  & & (0.017) & (0.021) & (0.023) & (0.017) & (0.021) & (0.022) \\ \cline{2-8}  & \multirow{2}{*}{700} & **0.132** & 0.241 & 0.327 & **0.131** & 0.241 & 0.329 \\  & & (0.011) & (0.014) & (0.016) & (0.011) & (0.014) & (0.016) \\ \hline \multirow{4}{*}{2D Mixed 1} & 200 & **0.090** & 0.150 & 0.229 & **0.089** & 0.150 & 0.232 \\  & & (0.018) & (0.021) & (0.025) & (0.017) & (0.022) & (0.026) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{350} & **0.076** & 0.144 & 0.229 & **0.075** & 0.144 & 0.232 \\  & & (0.010) & (0.016) & (0.020) & (0.010) & (0.016) & (0.019) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{700} & **0.070** & 0.142 & 0.232 & **0.069** & 0.142 & 0.235 \\  & & (0.009) & (0.011) & (0.014) & (0.008) & (0.011) & (0.015) \\ \hline \multirow{4}{*}{2D Mixed 2} & 200 & **0.121** & 0.152 & 0.201 & **0.119** & 0.152 & 0.204 \\  & & (0.019) & (0.022) & (0.025) & (0.018) & (0.022) & (0.025) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{350} & **0.116** & 0.152 & 0.202 & **0.114** & 0.151 & 0.206 \\  & & (0.014) & (0.016) & (0.021) & (0.013) & (0.016) & (0.021) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{700} & **0.108** & 0.148 & 0.199 & **0.069** & 0.148 & 0.203 \\  & & (0.009) & (0.012) & (0.012) & (0.009) & (0.012) & (0.013) \\ \hline \multirow{4}{*}{2D Mixed 3} & 200 & **0.107** & 0.116 & 0.264 & **0.099** & 0.116 & 0.264 \\  & & (0.024) & (0.018) & (0.023) & (0.020) & (0.019) & (0.024) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{350} & **0.098** & 0.116 & 0.265 & **0.090** & 0.115 & 0.267 \\  & & (0.012) & (0.014) & (0.021) & (0.012) & (0.014) & (0.021) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{700} & **0.097** & 0.114 & 0.265 & **0.085** & 0.114 & 0.265 \\  & & (0.012) & (0.011) & (0.016) & (0.010) & (0.011) & (0.016) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Averaged misclassification rates with standard errors in brackets for 2D simulations when \(m=100\) and \(m=400\) over \(100\) replicates.
\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \multirow{2}{*}{Model} & \multirow{2}{*}{\(n_{k}\)} & \multicolumn{3}{c}{\(m=8\)} & \multicolumn{3}{c}{\(m=27\)} \\ \cline{3-8}  & & mfDNN & MSDA & PLDA & mfDNN & MSDA & PLDA \\ \hline \multirow{8}{*}{3D Gaussian} & 200 & **0.374** & 0.465 & 0.485 & **0.234** & 0.367 & 0.491 \\  & & (0.024) & (0.037) & (0.045) & (0.020) & (0.030) & (0.046) \\ \cline{2-8}  & \multirow{2}{*}{350} & **0.373** & 0.474 & 0.495 & **0.233** & 0.369 & 0.500 \\  & & (0.020) & (0.029) & (0.042) & (0.020) & (0.023) & (0.041) \\ \cline{2-8}  & \multirow{2}{*}{700} & **0.363** & 0.472 & 0.488 & **0.232** & 0.374 & 0.493 \\  & & (0.014) & (0.025) & (0.041) & (0.013) & (0.018) & (0.041) \\ \hline \multirow{8}{*}{3D Mixed 1} & 200 & **0.215** & 0.224 & 0.238 & **0.168** & 0.164 & 0.242 \\  & & (0.022) & (0.024) & (0.022) & (0.021) & (0.022) & (0.022) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{350} & **0.200** & 0.221 & 0.236 & **0.151** & 0.156 & 0.240 \\  & & (0.018) & (0.019) & (0.020) & (0.017) & (0.016) & (0.020) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{700} & **0.188** & 0.220 & 0.237 & **0.149** & 0.152 & 0.241 \\  & & (0.011) & (0.014) & (0.012) & (0.010) & (0.011) & (0.013) \\ \hline \multirow{8}{*}{3D Mixed 2} & 200 & **0.311** & 0.325 & 0.346 & **0.277** & 0.294 & 0.348 \\  & & (0.018) & (0.025) & (0.026) & (0.020) & (0.022) & (0.025) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{350} & **0.307** & 0.328 & 0.345 & **0.268** & 0.289 & 0.348 \\  & & (0.018) & (0.021) & (0.022) & (0.018) & (0.018) & (0.022) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{700} & **0.253** & 0.326 & 0.347 & **0.234** & 0.286 & 0.350 \\  & & (0.014) & (0.015) & (0.015) & (0.013) & (0.014) & (0.015) \\ \hline \multirow{8}{*}{3D Mixed 3} & 200 & **0.272** & 0.308 & 0.315 & **0.268** & 0.292 & 0.313 \\  & & (0.026) & (0.031) & (0.030) & (0.026) & (0.026) & (0.029) & (0.030) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{350} & **0.268** & 0.303 & 0.309 & **0.239** & 0.286 & 0.310 \\  & & (0.019) & (0.023) & (0.024) & (0.018) & (0.024) & (0.023) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{700} & **0.253** & 0.302 & 0.311 & **0.234** & 0.286 & 0.312 \\  & & (0.014) & (0.017) & (0.017) & (0.013) & (0.015) & (0.017) \\ \hline \hline \end{tabular}
\end{table}
Table 3: Averaged misclassification rates with standard errors in brackets for 3D simulations when \(m=8\) and \(m=27\) over \(100\) replicates.
\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \multirow{2}{*}{Model} & \multirow{2}{*}{\(n_{k}\)} & \multicolumn{3}{c}{\(m=64\)} & \multicolumn{3}{c}{\(m=125\)} \\ \cline{3-8}  & & mfDNN & MSDA & PLDA & mfDNN & MSDA & PLDA \\ \hline \multirow{4}{*}{3D Gaussian} & 200 & **0.135** & 0.369 & 0.493 & **0.134** & 0.369 & 0.494 \\  & & (0.018) & (0.029) & (0.048) & (0.017) & (0.029) & (0.048) \\ \cline{2-8}  & \multirow{2}{*}{350} & **0.123** & 0.369 & 0.502 & **0.118** & 0.370 & 0.501 \\  & & (0.015) & (0.022) & (0.042) & (0.012) & (0.022) & (0.042) \\ \cline{2-8}  & \multirow{2}{*}{700} & **0.114** & 0.374 & 0.495 & **0.108** & 0.374 & 0.495 \\  & & (0.010) & (0.019) & (0.04) & (0.009) & (0.019) & (0.042) \\ \hline \multirow{4}{*}{3D Mixed 1} & 200 & **0.130** & 0.163 & 0.243 & **0.127** & 0.162 & 0.244 \\  & & (0.021) & (0.022) & (0.022) & (0.020) & (0.022) & (0.023) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{350} & **0.109** & 0.156 & 0.241 & **0.106** & 0.155 & 0.241 \\  & & (0.014) & (0.016) & (0.020) & (0.014) & (0.016) & (0.020) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{700} & **0.098** & 0.153 & 0.243 & **0.098** & 0.152 & 0.242 \\  & & (0.010) & (0.011) & (0.013) & (0.009) & (0.011) & (0.013) \\ \hline \multirow{4}{*}{3D Mixed 2} & 200 & **0.166** & 0.294 & 0.350 & **0.154** & 0.295 & 0.348 \\  & & (0.020) & (0.022) & (0.025) & (0.019) & (0.024) & (0.026) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{350} & **0.147** & 0.289 & 0.349 & **0.135** & 0.289 & 0.349 \\  & & (0.015) & (0.018) & (0.023) & (0.013) & (0.018) & (0.023) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{700} & **0.141** & 0.286 & 0.351 & **0.126** & 0.286 & 0.351 \\  & & (0.012) & (0.015) & (0.014) & (0.011) & (0.014) & (0.014) \\ \hline \multirow{4}{*}{3D Mixed 3} & 200 & **0.184** & 0.292 & 0.312 & **0.176** & 0.292 & 0.312 \\  & & (0.022) & (0.029) & (0.030) & (0.021) & (0.028) & (0.030) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{350} & **0.171** & 0.287 & 0.310 & **0.162** & 0.286 & 0.310 \\  & & (0.016) & (0.023) & (0.024) & (0.016) & (0.024) & (0.023) \\ \cline{1-1} \cline{2-8}  & \multirow{2}{*}{700} & **0.160** & 0.286 & 0.312 & **0.151** & 0.285 & 0.312 \\  & & (0.011) & (0.015) & (0.017) & (0.012) & (0.015) & (0.017) \\ \hline \hline \end{tabular}
\end{table}
Table 4: Averaged misclassification rates with standard errors in brackets for 3D simulations when \(m=64\) and \(m=125\) over \(100\) replicates.
## 6 Real data analysis

### Handwritten digits

The first benchmark data example was extracted from the MNIST database ([http://yann.lecun.com/exdb/mnist/](http://yann.lecun.com/exdb/mnist/)). This classical MNIST database contains 60,000 training images and 10,000 testing images of handwritten digits (\(0,1,\ldots,9\)), and the black and white images were normalized to fit into a \(28\times 28\) pixel bounding box and anti-aliased. We used tensor of Fourier basis for data processing. According to our numerical experience, we choose candidates for \((L,J,\mathbf{p},s)\), such that \(L=(2,3,4)^{\intercal}\), \(J=(300,500,800)^{\intercal}\), \(\|\mathbf{p}\|_{\infty}~{}=~{}(500,1000,2000)^{\intercal}\), and \(s=(0.01,0.1,0.5)\) for dropout rate. Here we abuse the notation of \(s\), as the dropout is the technique we choose to sparsify the neural network. With the optimal parameters \(L_{opt}=3\), \(J_{opt}=500\), \(\|\mathbf{p}_{opt}\|_{\infty}~{}=1000\), \(s_{opt}=0.01\) through validation, we demonstrate the misclassification risk in Table 5. We estimate the rules given by MSDA, PLDA and our proposal on the training set. As most observations for each subject are zeros, PenalizedLDA reports errors and does not work any more. It can be seen that our proposal achieves the higher accuracy with the sparsest classification rule. This again shows that our method is a competitive classifier and has more broader applications.

### ADNI database

The dataset used in the preparation of this article were obtained from the ADNI database ([http://adni.loni.usc.edu](http://adni.loni.usc.edu)). The ADNI is a longitudinal multicenter study designed to develop clinical, imaging, genetic, and biochemical biomarkers for the

Figure 1: Samples from MNIST data.
early detection and tracking of Alzheimer's Disease (AD). From this database, we collect PET data from \(79\) patients in AD group, \(45\) patients in Early Mild Cognitive Impairment (EMCI) group, and \(101\) people in Control (CN) group. This PET dataset has been spatially normalized and post-processed. These AD patients have three to six times doctor visits and we select the PET scans obtained in the third visits. People in EMCI group only have the second visit, and we select the PET scans obtained in the second visits. For AD group, patients' age ranges from \(59\) to \(88\) and average age is \(76.49\), and there are \(33\) females and \(46\) males among these \(79\) subjects. For EMCI group, patients' age ranges from \(57\) to \(89\) and average age is \(72.33\), and there are \(26\) females and \(19\) males among these \(45\) subjects. For CN group, patients' age ranges from \(62\) to \(87\) and average age is \(75.98\), and there are \(40\) females and \(61\) males among these \(101\) subjects. All scans were reoriented into \(79\times 95\times 68\) voxels, which means each patient has \(68\) sliced 2D images with \(79\times 95\) pixels. For 2D case, it means each subject has \(N=79\times 95=7,505\) observed pixels for each selected image slice. For 3D case, the observed number of voxels for each patient's brain image observation is \(N=79\times 95\times 68=510,340\). We randomly split the datasets with a \(7:3\) ratio in a balanced manner to form the training set and the testing set, with \(100\) repetitions.

We choose candidates for \((L,J,\boldsymbol{p},s)\), such that \(L=(2,3)^{\intercal}\), \(J=(50,100,200)^{\intercal}\), \(\|\boldsymbol{p}\|_{\infty}\ =(200,500,800)^{\intercal}\), and \(s=(0.01,0.1,0.5)\) for dropout rate. We still compare our method with MSDA. For 2D case, it means each subject has \(N=79\times 95=7,505\) observed pixels for each selected image slice. Table 6 displays the misclassification rates for 2D brain imaging data of AD, EMCI and CN. For 3D case, the observed number of voxels for each patient's brain sample is \(N=79\times 95\times 68=510,340\). Unfortunately, given more than a half million covariates, MSDA method crashed down as such gigantic covariance matrix (almost million by million dimension) needs around 2TB RAM to store. It easily exceeds the memory limit of the common supercomputer. Hence, MSDA for 3D data results are unavailable. Meanwhile, as PLDA recasts Fisher's discriminant problem as a biconvex problem that can be optimized using a simple iterative algorithm, PLDA avoids the heavy computation burden of the covariance matrix and it still works for this 3D case. Table 7 presents the empirical misclassification risk for mfDNN and PLDA.

There are several interesting findings in Tables 6 and 7. First, our proposed classifier has better performance than other competitors in any 2D slice data or 3D data cases. Second, from Table 7, we can conclude that given a single slice of 2D imaging data, the misclassification rates are consistently larger than using the 3D data. It indicates that 3D data contains more helpful information to label the brain images among three stages of the disease. Third, the \(10\)-th and the \(20\)-th slices provide the lowest ones among all 2D data. As a matter of fact, it is well known that Alzheimer's disease destroys neurons and their connections in hippocampus, the entorhinal cortex, and the cerebral cortex. The related regions are corresponding to the first \(25\) slices. This is a promising finding for neurologists, as this smallest risk indicates this particular slice presents useful information to distinguish the CN, EMCI and AD groups. Further medical checkups are meaningful for this special location in the brain.
\begin{table}
\begin{tabular}{l c c c c c c} \hline \hline methods & \(10\)-th & \(20\)-th & \(30\)-th & \(40\)-th & \(50\)-th & \(60\)-th \\ \hline mfDNN & **0.341** & **0.325** & **0.409** & **0.359** & **0.365** & **0.452** \\  & (0.057) & (0.049) & (0.049) & (0.050) & (0.041) & (0.055) \\ \hline MSDA & 0.347 & 0.374 & 0.444 & 0.366 & 0.385 & 0.453 \\  & (0.038) & (0.034) & (0.035) & (0.034) & (0.039) & (0.040) \\ \hline PLDA & 0.354 & 0.339 & 0.440 & 0.376 & 0.424 & 0.474 \\  & (0.057) & (0.059) & (0.067) & (0.060) & (0.068) & (0.064) \\ \hline \end{tabular}
\end{table}
Table 6: Averaged misclassification rates with standard errors in brackets for ADNI 2D brain images.

\begin{table}
\begin{tabular}{l c} \hline \hline methods & Misclassification rates \\ \hline mfDNN & **0.274 (0.044)** \\ \hline MSDA & – (–) \\ \hline PLDA & 0.295 (0.056) \\ \hline \end{tabular}
\end{table}
Table 7: Averaged misclassification rates with standard errors in brackets for ADNI 3D brain images.

Figure 2: Averaged images of the the \(20\)-th, the \(40\)-th and the \(60\)-th slices of EMCI (left column) group and AD group (right column).
## 7 Summary

In this paper, we propose a new formulation to derive multiclass classifiers for multidimensional functional data. We show that our proposal has a solid theoretical foundation and can be solved by a very efficient computational algorithm. Our proposal actually gives a unified treatment of both one-dimensional and multidimensional classification problems. In light of this evidence, our proposal is regarded as an efficient multiclass and generalization of the multiclassification methods from i.i.d. data to multidimensional functional data cases. To the best of our knowledge, the present work is the first work on multiclassification for multidimensional functional data with theoretical justification.

## Acknowledgements

ADNI data 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](http://adni.loni.usc.edu/wp-content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf).
