Theory of functional principal component analysis
for discretely observed data

Modern data collection technologies have rapidly evolved, leading to the widespread emergence of functional data that have been extensively studied over the past few decades. Generally, functional data are considered stochastic processes that satisfy certain smoothness conditions or realizations of random elements valued in Hilbert space. These two perspectives highlight the essential natures of functional data, namely, their smoothness and infinite dimensionality, which distinguish them from high-dimensional and vector-valued data. For a comprehensive treatment of functional data, we recommend the monographs by Ramsay and Silverman (2006), Ferraty and Vieu (2006), Horváth and Kokoszka (2012), and Hsing and Eubank (2015), among others.

Although functional data provide information over a continuum, which is often time or spatial location, in reality, data are collected or observed discretely with measurement errors. For instance, we usually use $n$ to denote the sample size, which is the number of subjects corresponding to random functions, and $N_{i}$ to denote the number of observations for the $i$th subject. Thanks to the smooth nature of functional data, having a large number of observations per subject is more of a blessing than a curse, in contrast to high-dimensional data (Hall, Müller and Wang, 2006; Zhang and Wang, 2016). There is extensive literature on nonparametric methods that address the smoothness of functional data, including kernel or local polynomial methods (Yao, Müller and Wang, 2005a; Hall, Müller and Wang, 2006; Zhang and Wang, 2016), and various types of spline methods (Rice and Wu, 2001; Yao and Lee, 2006; Paul and Peng, 2009; Cai and Yuan, 2011).

When employing a smoothing method, there are two typical strategies to be considered. If the observed time points per subject are relatively dense, it is recommended to pre-smooth each curve before further analysis, as suggested by Ramsay and Silverman (2006) and Zhang and Chen (2007). However, if the sampling scheme is rather sparse, it is preferred to pool observations together from all subjects (Yao, Müller and Wang, 2005a). The choice of individual versus pooled information affects the convergence rates and phase transitions in estimating population quantities, such as mean and covariance functions. When $N_{i}\gtrsim O(n^{5/4})$ and the tuning parameter is optimally chosen per subject, the estimated mean and covariance functions based on the reconstructed curves through pre-smoothing are $\sqrt{n}$-consistent, the so-called optimal parametric rate. On the other hand, by borrowing information from all subjects, the pooling method only requires $N_{i}\gtrsim O(n^{1/4})$ for mean and covariance estimation to reach optimal (Cai and Yuan, 2010, 2011; Zhang and Wang, 2016), providing theoretical insight into the advantage of the pooling strategy.

However, estimating the mean and covariance functions does not account for the infinite dimensionality of functional data. Due to the non-invertibility of covariance operators for functional random objects, regularization is required in models that involve inverse issues with functional covariates, such as the functional linear model (Yao, Müller and Wang, 2005b; Hall and Horowitz, 2007; Yuan and Cai, 2010), functional generalized linear model (Müller and Stadtmüller, 2005; Dou, Pollard and Zhou, 2012), and functional Cox model (Qu, Wang and Wang, 2016). Truncation of the leading functional principal components (FPC) is a well-developed approach to addressing this inverse issue (Hall and Horowitz, 2007; Dou, Pollard and Zhou, 2012). In order to suppress the model bias, the number of principal components used in truncation should grow slowly with sample size. As a result, the convergence rate of the estimated eigenfunctions with diverging indices becomes a fundamental issue, which is not only important in its own right but also crucial for most models and methods involving functional principal components regularization.

For fully observed functional data, Hall and Horowitz (2007) obtained the optimal convergence rate $j^{2}/n$ for the $j$th eigenfunction, which served as a cornerstone in establishing the optimal convergence rate in functional linear model (Hall and Horowitz, 2007) and functional generalized linear model (Dou, Pollard and Zhou, 2012). In the discretely observed case, stochastic bounds for a fixed number of eigenfunctions have been obtained by different methods. Using a local linear smoother, Hall, Müller and Wang (2006) showed that the $\mathcal{L}^{2}$ rate of a fixed eigenfunction for finite $N_{i}$ is $O_{P}(n^{-4/5})$. Under the reproducing kernel Hilbert space framework, Cai and Yuan (2010) claimed that eigenfunctions with fixed indices admit the same convergence rate as the covariance function, which is $O_{P}((n/\log n)^{-4/5})$. It is important to note that, although both results are one-dimensional nonparametric rates (differing at most by a factor of $(\log n)^{4/5}$), the methodologies and techniques used are completely disparate, and a detailed discussion can be found in Section 2. Additionally, Paul and Peng (2009) proposed a reduced rank model and studied its asymptotic properties under a particular setting. In Zhou, Yao and Zhang (2023), the authors studied the convergence rate for the functional linear model and obtained an improved bound for the eigenfunctions with diverging indices. However, this rate will not reach the optimal rate of $j^{2}/n$ for any sampling rate $N_{i}$. As explained in Section 2, while some bounds can be obtained for eigenfunctions with diverging indices, attaining an optimal bound presents a substantially greater challenge. Lack of such an optimal bound for eigenfunctions poses considerable challenge in analyzing the standard and efficient plug-in estimator in functional linear model (Hall and Horowitz, 2007). Consequently, Zhou, Yao and Zhang (2023) resorted to a complex sample-splitting strategy, which results in lower estimation efficiency. To the best of our knowledge, there has been no progress in obtaining the optimal convergence rate of eigenfunctions with diverging indices when the data are discretely observed with noise contamination.

The distinction between estimating a diverging number and a fixed number of eigenfunctions is rooted in the infinite-dimensional nature of functional data. Analyzing eigenfunctions with diverging indices presents challenges due to the decaying eigenvalues. For fully observed data, the cross-sectional sample covariance based on the true functions facilitates the application of perturbation results, as shown in prior work (Hall and Horowitz, 2007; Dou, Pollard and Zhou, 2012). This approach simplifies each term in the perturbation series to the principal component scores. However, when the trajectories are observed at discrete time points, this virtue no longer exists, leading to a summability issue arising from the estimation bias and decaying eigenvalues. This renders existing techniques invalid and remains an unresolved problem; see Section 2 for further elaboration.

This paper addresses this significant yet challenging task of estimating an increasing number of eigenfunctions from discretely observed functional data, and presents a unified theory. The contributions of this paper are at least threefold. First, we establish an $\mathcal{L}^{2}$ bound for the eigenfunctions and the asymptotic normality of the eigenvalues with increasing indices, reflecting a transition from nonparametric to parametric regimes and encompassing a wide range from sparse to dense sampling. We show that when $N_{i}$ reaches a magnitude of $n^{1/4+\delta}$, where $\delta$ depends on the smoothness parameters of the underlying curves, the convergence rate becomes optimal as if the curves are fully observed. Second, we introduce a novel double truncation method that yields uniform convergence across the time domain, surmounting theoretical barriers in the existing literature. Through this approach, uniform convergence rates for the covariance and eigenfunctions are achieved under mild conditions across various sampling schemes. Notably, this includes the uniform convergence of eigenfunctions with increasing indices, which is new even in scenarios where data are fully observed. Third, we provide a new technical route for addressing the perturbation series of the functional covariance operator, bridging the gap between the “ideal” fully observed scenario and the noisy, discrete “real-world” context. These advanced techniques pave the way for their application in downstream FPCA-related analyses, and the achieved optimal rate of eigenpairs facilitates the extension of existing theoretical results for “fully” observed functional data to discreetly observed case.

The rest of the paper is organized as follows. In Section 2, we give a synopsis of covariance and eigencomponents estimation in functional data. We present the $\mathcal{L}^{2}$ convergence of eigenfunctions in Section 3, and discuss the uniform convergence problem of functional data in Section 5. Asymptotic normality of eigenvalues is presented in Section 4. Section 6 provides an illustration of the phase transition phenomenom in eigenfunctions with synthetic data. The proofs of Theorem 1 can be found in Appendix, while the proofs of other theorems and lemmas are collected in the Supplementary Material.

In what follows, we denote by $A_{n}=O_{P}(B_{n})$ the relation $\mathbb{P}(A_{n}\leq MB_{n})\geq 1-\epsilon$, and by $A_{n}=o_{P}(B_{n})$ the relation $\mathbb{P}(A_{n}\leq\epsilon B_{n})\rightarrow 0$ as $n\rightarrow\infty$, for each $\epsilon>0$ and a positive constant $M$. A non-random sequence $a_{n}$ is said to be $O(1)$ if it is bounded. For any non-random sequence $b_{n}$, we say $b_{n}=O(a_{n})$ if $b_{n}/a_{n}=O(1)$, and $b_{n}=o(a_{n})$ if $b_{n}/a_{n}\rightarrow 0$. The notation $a_{n}\lesssim b_{n}$ indicates $a_{n}\leq Cb_{n}$ for sufficiently large $n$ and postive constant $C$, and the relation $\gtrsim$ is defined similarly. We write $a_{n}\asymp b_{n}$ if $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$. For $a\in\mathbb{R}$, $\lfloor a\rfloor$ denotes the largest integer less than or equal to $a$. For a function $f\in\mathcal{L}^{2}[0,1]$, where $\mathcal{L}^{2}[0,1]$ denotes the space of square-integrable functions on $[0,1]$, $\|f\|^{2}$ denotes $\int_{[0,1]}f(s)^{2}\mathrm{d}s$, and $\|f\|_{\infty}$ denotes $\sup_{s\in[0,1]}|f(s)|$. For a function $A(s,t)\in\mathcal{L}^{2}[0,1]^{2}$, define $\|A\|_{\mathrm{HS}}^{2}=\iint_{[0,1]^{2}}A(s,t)^{2}\mathrm{d}s\mathrm{d}t$ and $\|A\|_{(j)}^{2}=\int_{[0,1]}\{\int_{[0,1]}A(s,t)\phi_{j}(s)\mathrm{d}s\}^{2}\mathrm{d}t$, where $\{\phi_{j}\}_{j=1}^{\infty}$ are the eigenfunctions of interest. We write $\int pq$ and $\int Apq$ for $\int p(u)q(u)\mathrm{d}u$ and $\iint A(u,v)p(u)q(v)\mathrm{d}u\mathrm{d}v$ occasionally for brevity.

Let $X(t)$ be a square integrable stochastic process on $[0,1]$, and let $X_{i}(t)$ be independent and identically distributed copies of $X(t)$. The mean and covariance functions of $X(t)$ are denoted by $\mu(t)=\mathbb{E}[{X(t)}]$ and $C(s,t)=\mathbb{E}[{\{X(s)-\mu(s)\}\{X(t)-\mu(t)\}}]$, respectively. According to Mercer’s Theorem (Indritz, 1963), $C(s,t)$ has the spectral decomposition

where $\lambda_{1}>\lambda_{2}>\ldots>0$ are eigenvalues and $\{\phi_{j}\}_{j=1}^{\infty}$ are the corresponding eigenfunctions, which form a complete orthonormal system on $\mathcal{L}^{2}[0,1]$. For each $i$, the process $X_{i}$ admits the so-called Karhunen-Loève expansion

where $\xi_{ik}=\int_{0}^{1}\{X_{i}(t)-\mu(t)\}\phi_{k}(t)\mathrm{d}t$ are functional principal component scores with zero mean and variances $\lambda_{k}$.

However, in practice, it is only an idealization to have each $X_{i}(t)$ for all $t\in[0,1]$ to simplify theoretical analysis. Measurements are typically taken at $N_{i}$ discrete time points with noise contamination. Specifically, the actual observations for each $X_{i}$ are given by

where $\varepsilon_{ij}$ are random copies of $\varepsilon$, with $\mathbb{E}(\varepsilon)=0$ and $\text{Var}(\varepsilon)=\sigma_{X}^{2}$. We further assume the measurements errors $\{\varepsilon_{ij}\}_{i,j}$ are independent of $X_{i}$.

Local linear regression is a popular smoothing technique in functional data analysis due to its attractive theoretical properties (Yao, Müller and Wang, 2005a; Hall, Müller and Wang, 2006; Li and Hsing, 2010; Zhang and Wang, 2016). The primary goal of this paper is to develop a unified theory for estimating a larger number of eigenfunctions from discretely observed functional data. To maintain focus and avoid distractions, we assume that the mean function $\mu(t)$ is known, and set $\mu(t)=0$ without loss of generality. The scenario involving an unknown mean function is discussed later in Section 3. We denote by $\delta_{ijl}=X_{ij}X_{il}$ the raw covariance, and define $v_{i}=\{nN_{i}(N_{i}-1)\}^{-1}$. The local linear estimator of the covariance function is given by $\hat{C}(s,t)=\hat{\beta}_{0}$,

where $\mathrm{K}$ is a symmetric, Lipschitz continuous density kernel on $[-1,1]$ and $h$ is the tuning parameter. The estimated covariance function $\hat{C}(s,t)$ can be expressed as an empirical version of the spectral decomposition in (1), i.e.,

where $\hat{\lambda}_{k}$ and $\hat{\phi}_{k}$ are estimators for $\lambda_{k}$ and $\phi_{k}$, respectively. We assume that $\langle\hat{\phi}_{k},\phi_{k}\rangle\geq 0$ for specificity.

Before delving into the theoretical details, we provide an overview of eigenfunction estimation in functional data analysis. We start with the resolvent series shown in Equation (6) and illustrate its application in statistical analysis,

Such expansions can be found in Bosq (2000), Dou, Pollard and Zhou (2012), and Li and Hsing (2010); see Chapter 5 in Hsing and Eubank (2015) for details. Denote by $\eta_{j}$ the eigengap of $\lambda_{j}$, that is, $\eta_{j}:=\min_{k\neq j}|\lambda_{k}-\lambda_{j}|$. An basic rough bound for $\mathbb{E}(\|\hat{\phi}_{j}-\phi_{j}\|^{2})$ can be derived from Equation (6) and Bessel’s inequality,

However, this bound is suboptimal for two reasons. First, while $\eta_{j}$ is bounded away from $0$ for a fixed $j$, the bound implies that the eigenfunctions converge at the same rate as the covariance function. This is counterintuitive since integration usually brings extra smoothness (Cai and Hall, 2006), which typically results in the eigenfunction estimates converging at a faster rate than the two-dimensional kernel smoothing rate of $\|\hat{C}-C\|_{\mathrm{HS}}^{2}$. Second, for $j$ that diverges with the sample size, $\eta_{j}^{-2}\rightarrow\infty$ in the bound cannot be improved. Assuming $\lambda_{j}\asymp j^{-a}$, the rate obtained by (7) is slower than $j^{2a+2}/n$, which is known to be suboptimal (Wahl, 2022). Therefore, to obtain a sharp rate for $\mathbb{E}(\|\hat{\phi}_{j}-\phi_{j}\|^{2})$, one should use the original perturbation series given by (6), rather than its approximation given by (7).

When each trajectory $X_{i}(t)$ is fully observed for all $t\in[0,1]$, the cross-sectional sample covariance $\hat{C}(s,t)=n^{-1}\sum_{i=1}^{n}X_{i}(s)X_{i}(t)$ is a canonical estimator of $C(s,t)$. Then, the numerators in each term of (6) can be reduced to the principal components scores under some mild assumptions, for example, $\mathbb{E}[\{\iint(\hat{C}-C)\phi_{j}\phi_{k}\}^{2}]\lesssim n^{-1}\lambda_{j}\lambda_{k}$ (Hall and Horowitz, 2007; Dou, Pollard and Zhou, 2012). Subsequently, $\mathbb{E}(\|\hat{\phi}_{j}-\phi_{j}\|^{2})$ is bounded by $(\lambda_{j}/n)\sum_{k\neq j}\lambda_{k}/(\lambda_{k}-\lambda_{j})^{2}$. With the common assumption of the polynomial decay of eigenvalues, the aforementioned summation is dominated by $\lambda_{j}/n\sum_{\lfloor j/2\rfloor\leq k\leq 2j}\lambda_{k}/(\lambda_{k}-\lambda_{j})^{2}$, which is $O(j^{2}/n)$ and optimal in the minimax sense (Wahl, 2022). See Lemma 7 in Dou, Pollard and Zhou (2012) for a detailed elaboration. This suggests that the convergence rate caused by the inverse issue can be captured by the summation over the set $\{k\leq 2j\}$, and the tail sum on $\{k>2j\}$ can be treated as a unity.

However, we would like to emphasize that when it comes to discretely observed functional data, all the existing literature utilizing a bound similar to (7) excludes the case of diverging indices. For instance, the result in Cai and Yuan (2010) is simply a direct application of the bound in (7). Moreover, their one-dimensional rate is inherited from the covariance estimator, which is assumed to be in a tensor product space that is smaller than the space $\mathcal{L}^{2}[0,1]^{2}$. On the other hand, the one-dimensional rates obtained by Hall, Müller and Wang (2006) and Li and Hsing (2010) utilize detailed calculations based on the approximation of the perturbation series in (6). However, these results are based on the assumption that $\eta_{j}$ is bounded away from zero, which implies that $j$ must be a fixed constant. This is inconsistent with the nonparametric nature of functional data models, which aim to approximate or regularize an infinite-dimensional process. Therefore, when dealing with discretely observed functional data, the key to obtaining a sharp bound for estimated eigenfunctions with diverging indices lies in effectively utilizing the perturbation series (6).

The main challenges arise from quantifying the summation in (6) without the fully observed sample covariance. For the pre-smoothing method, the reconstructed $\hat{X}_{i}$ achieves a $\sqrt{n}$ convergence in the $\mathcal{L}^{2}$ sense when each $N_{i}$ reaches a magnitude of $n^{5/4}$, and then the estimated covariance function $\hat{C}(s,t)=n^{-1}\sum_{i=1}^{n}\hat{X}_{i}(s)\hat{X}_{i}(t)$ has an optimal rate $\|\hat{C}-C\|_{\mathrm{HS}}=O_{P}(n^{-1/2})$. However, this does not guarantee optimal convergence of a diverging number of eigenfunctions. The numerators in each term of (6) are no longer the principal component scores, and the complex form of this infinite summation makes it difficult to quantify when $|\lambda_{k}-\lambda_{j}|\rightarrow 0$. Similarly, the pooling method also encounters significant challenges in summing all $\mathbb{E}[\{\iint(\hat{C}-C)\phi_{j}\phi_{k}\}^{2}]$ with respect to $j$ and $k$. Specifically, the convergence rate of $\|\hat{C}-C\|_{\mathrm{HS}}^{2}$ should be a two-dimensional kernel smoothing rate with variance $n^{-1}\{1+(Nh)^{-2}\}$ (Zhang and Wang, 2016). However, after being integrated twice, $\mathbb{E}[\{\iint(\hat{C}-C)\phi_{j}\phi_{k}\}^{2}]$ has a degenerated kernel smoothing rate with variance $n^{-1}$. According to Bessel’s inequality, $\mathbb{E}(\|\hat{C}-C\|_{\mathrm{HS}}^{2})$ can be expressed as $\sum_{j,k}\mathbb{E}[\{\iint(\hat{C}-C)\phi_{j}\phi_{k}\}^{2}]$. However, due to estimation bias, one cannot directly sum all $\mathbb{E}[\{\iint(\hat{C}-C)\phi_{j}\phi_{k}\}^{2}]$ with respect to all $j,k$.

Based on the issues discussed above, we propose a novel technique for addressing the perturbation series (6) when dealing with discretely observed functional data. To this end, we make the following regularity assumptions.

Assumptions 1 and 2 are widely adopted in the functional data literature related to smoothing (Yao, Müller and Wang, 2005a; Cai and Yuan, 2010; Zhang and Wang, 2016). The decay rate assumption on the eigenvalues provides a natural theoretical framework for justifying and analyzing the properties of functional principal components (Hall and Horowitz, 2007; Dou, Pollard and Zhou, 2012; Zhou, Yao and Zhang, 2023). The number of eigenfunctions that can be well estimated from exponentially decaying eigenvalues is limited to the order of $\log n$, which lacks practical interest. Consequently, we adopt the assumption of polynomial decay in eigenvalues. To achieve quality estimates for a specific eigenfunction, its smoothness should be taken into account. Generally, the frequency of $\phi_{j}$ is higher for larger $j$, which requires a smaller bandwidth to capture its local variation. Assumption 4 characterizes the frequency increment of a specific eigenfunction via the smoothness of its derivatives. For some commonly used bases, such as the Fourier, Legendre, and wavelet bases, $c=2$. In Hall, Müller and Wang (2006), the authors assumed that $\max_{1\leq j\leq r}\max_{s=0,1,2}\sup_{t\in[0,1]}|\phi_{j}^{(s)}(t)|\leq C$, which is only achievable for a fixed $r$. Therefore, Assumption 4 can be interpreted as a generalization of this assumption. To analyze the convergence of eigenfunctions effectively, a uniform convergence rate of the covariance function is needed to handle the local linear estimator (4). Assumption 5 is the moment assumption required for uniform convergence of covariance function and adopted in Li and Hsing (2010) and Zhang and Wang (2016).

For the observation time points $\{t_{ij}\}_{i,j}$, there are two typical types of designs: the random design, in which the design points are random samples from a distribution, and the regular design, where the observation points are predetermined mesh grid. For the random design, the following assumption is commonly adopted (Yao, Müller and Wang, 2005a; Li and Hsing, 2010; Cai and Yuan, 2011; Zhang and Wang, 2016):

For the regular design, each sample path is observed on an equally spaced grid $\{t_{j}\}_{j=1}^{N}$, where $N_{i}=\cdots=N_{n}=N$ for all subjects. This longitudinal design is frequently encountered in a various scientific experiments and has been studied in Cai and Yuan (2011). Assumption 7 guarantees a sufficient number of observations within the local window for the kernel smoothing method. Furthermore, Assumption 8 is needed for the Riemann sum approximation in the fixed regular design.

The following theorem is one of our main results. It characterizes the $\mathcal{L}^{2}$ convergence of the estimated eigenfunctions with diverging indices for both random design (Assumption 6) and fixed regular design (Assumption 7).

The integer $m$ in Theorem 1 represents the maximum number of eigenfunctions that can be well estimated from the observed data using appropriate tuning parameters. It is important to note that in Theorem 1, $m$ could diverge to infinity. The upper bound of $m$ is a function of the sample size $n$, the sampling rate $\bar{N}_{2}$, the smoothing parameter $h$, and the decaying eigengap $\eta_{j}$ or $a$. As the frequency of $\phi_{j}$ is higher for larger $j$, smaller $h$ is required to capture its local variations. If $a$ is large, the eigengap $\eta_{j}$ diminishes rapidly, posing a greater challenge in distinguishing between adjacent eigenvalues. Note that the assumptions of $m$ in Theorem 1 could encompass most downstream analyses that involve a functional covariate, such as functional linear regression as discussed in (Hall and Horowitz, 2007).

Theorem 1 provides a good illustration of both the infinite dimensionality and smoothness nature of functional data. To better understand this result, note that $j^{2}/n$ is the optimal rate in the fully observed case. The additional terms on the right-hand side of Equation (8) represent contamination introduced by discrete observations and measurement errors. In particular, the term with $h^{4}$ corresponds to the smoothing bias, while the term $(n\bar{N}_{2}h)^{-1}$ reflects the variance typically associated with one-dimensional kernel smoothing. Terms including $\bar{N}_{2}^{-1}$ arise from the discrete approximation, and the terms that involve $j$ with positive powers are due to the decaying eigengaps associated with an increasing number of eigencomponents.

The first part of Theorem 1 provides a unified result for eigenfunction estimates under random design without imposing any restrictions on the sampling rate $N_{i}$. Similar to the phase transitions of mean and covariance functions studied in Cai and Yuan (2011) and Zhang and Wang (2016), Corollary 1 presents a systematic partition that ranges from “sparse” to “dense” schemes for eigenfunction estimation under the random design, which is meaningful for FPCA-based models and methods.

Note that $h_{opt}(j)$ is the optimal bandwidth for estimating a specific eigenfunction $\phi_{j}$. However, in practice, it suffices to estimate the covariance function just once using the optimal bandwidth associated with the largest eigenfunction $\phi_{m_{\text{max}}}$. This ensures that both the subspace spanned by the first $m_{\text{max}}$ eigenfunctions and their corresponding projections are well estimated. When conducting downstream analysis with FPCA, the evaluation of error rates typically involves the term $\sum_{j=1}^{m_{\max}}j^{-\beta}\|\hat{\phi}_{j}-\phi_{j}\|^{2}$. Here $m_{\max}$ denotes the maximum number of principal components used, and $\beta$ varies across different scenarios, reflecting how the influence of each principal component on the error rate is weighted. For example, in the functional linear model, $\beta>0$ represents the rate at which the Fourier coefficients of the regression function decay. By the first statement of Theorem 1,

It can be found that the optimal bandwidth $h$ in the above equation aligns with $h_{opt}(m_{\max})$ for all $\beta$. This implies that our theoretical framework can be seamlessly adapted to scenarios where a single bandwidth is employed to attain the optimal convergence rate in the downstream analysis.

For the commonly used bases where $c=2$, the convergence rate for the $j$th eigenfunction achieves optimality as if the curves are fully observed when $\bar{N}_{2}>n^{1/4}j^{a-1}$. Keeping $j$ fixed, the phase transition occurs at $n^{1/4}$, aligning with results in Hall, Müller and Wang (2006) and Cai and Yuan (2010), as well as mean and covariance functions discussed in Zhang and Wang (2016). For $n$ subjects, the maximum index of the eigenfunction that can be well estimated is less than $m_{\mathrm{max}}:=n^{1/(2a+2)}$, which directly follows from the assumption $m^{2a+2}/n\rightarrow 0$.The phase transition in estimating $\phi_{m_{\mathrm{max}}}$ occurs at $n^{1/4+(a-1)/(2a+2)}$. This can be interpreted from two aspects. On one hand, compared to mean and covariance estimation, more observations per subject are required to obtain optimal convergence for eigenfunctions with increasing indices, showing the challenges tied to infinite dimensionality and decaying eigenvalues. On the other hand, the fact that $n^{1/4+(a-1)/(2a+2)}$ is only slightly larger than $n^{1/4}$ justifies the merits of the pooling method and supports our intuition. When $j$ is fixed and $\bar{N}_{2}$ is finite, the convergence rate obtained by Corollary 1 is $(nh)^{-1}+h^{4}$, which corresponds to a typical one-dimensional kernel smoothing rate and achieves optimal at $h\asymp n^{-1/5}$. This result aligns with those in Hall, Müller and Wang (2006) and is optimal in the minimax sense. When allowing $j$ to diverge, the known lower bound $j^{2}/n$ for fully observed data is attained by applying van Trees’ inequality to the special orthogonal group (Wahl, 2022). However, the argument presented in Wahl (2022) does not directly extend to the discrete case and there are currently no available lower bounds for the eigenfunctions with diverging indices based on discrete observations, which remains an open problem that requires further investigation.

Comparing the results of this work with those in Zhou, Yao and Zhang (2023) is also of interest. The $\mathcal{L}^{2}$ convergence rate for the $j$th eigenfunction obtained in (Zhou, Yao and Zhang, 2023) is

It is evident that the results in Zhou, Yao and Zhang (2023) will never reach the optimal rate $j^{2}/n$ for any sampling rate $N$. In contrast, Corollary 1 provide a systematic partition ranging from “sparse” to “dense” schemes for eigenfunction estimation, and the optimal rate $j^{2}/n$ can be achieved when the sampling rate $N_{i}$ exceeds the phase transition point. The optimal rate achieved here represents more than just a theoretical improvement over previous findings; it also carries substantial implications for downstream analysis. Further discussion can be found in Section 7.

The following Corollary presents the phase transition of eigenfunctions for fixed regular design. Note that for fixed regular design, the number of distinct observed time points in an interval of length $h$ is on the order of $Nh$, so $Nh\gtrsim 1$ is required to ensure there is at least one observation in the bandwidths of kernel smoothing (Shao, Lin and Yao, 2022). Moreover, $Nh^{a}\gtrsim 1$ is required to eliminate the Riemann sum approximation bias. The condition $Nh^{a}\gtrsim 1$ is parallel of part (a) in Corollary 1 in the scenario of the random design, which is similar as the mean and covariance estimation where consistency can only be obtained under the dense case for regular design (Shao, Lin and Yao, 2022).

If the mean function $\mu(t)$ is unknown, one could use local linear smoother to fit $\hat{\mu}(t)=\hat{\alpha}_{0}$ with

Then the covariance estimator $\hat{C}$ is obtained by replacing the raw covariance $\delta_{ijl}$ in (4) by $\hat{\delta}_{ijl}=\{X_{ij}-\hat{\mu}(t_{ij})\}\{X_{il}-\hat{\mu}(t_{il})\}$. The following corollary presents the convergence rate and phase transition for the case where $\mu(t)$ is unknown. It should be noted that the Fourier coefficients of $\mu(t)$ with respect to the eigenfunction ${\phi_{j}}$ generally do not exhibit a decaying trend. Therefore, to eliminate the estimation error caused by the mean estimation, an additional lower bound on $N$ is necessary. This lower bound, denoted as $N\gtrsim j^{a}$, aligns with the partition described in Corollary 1 for the random design case.

The distribution of eigenvalues plays a crucial role in statistical learning and is of significant interest in the high-dimensional setting. Random matrix theory provides a systematic tool for deriving the distribution of eigenvalues of a squared matrix (Anderson, Guionnet and Zeitouni, 2010; Pastur and Shcherbina, 2011), and has been successfully applied in various statistical problems, such as signal detection (Nadler, Penna and Garello, 2011; Onatski, 2009; Bianchi et al., 2011), spiked covariance models (Johnstone, 2001; Paul, 2007; El Karoui, 2007; Ding and Yang, 2021; Bao et al., 2022), and hypothesis testing (Bai et al., 2009; Chen and Qin, 2010; Zheng, 2012). For a comprehensive treatment of random matrix theory in statistics, we recommend the monograph by Bai and Silverstein (2010) and the review paper by Paul and Aue (2014).

Despite the success of random matrix theory in high-dimensional statistics, its application to functional data analysis is not straightforward due to the different structures of functional spaces. If the observations are taken at the same time points $\{t_{j}\}_{j=1}^{N}$ for all $i$, one can obtain an estimator for $\Sigma_{T}=\operatorname{Cov}(\tilde{X}_{i},\tilde{X}_{i})+\sigma_{X}^{2}I_{N}$, where $\tilde{X}_{i}=(X_{i}(t_{1}),\ldots,X_{i}(t_{N}))^{\top}$. Note that $\tilde{X}_{i}$ can be regarded as a random vector; however the adjacent elements in $\tilde{X}_{i}$ are highly correlated as $N$ increases due to the smooth nature of functional data. This correlation violates the independence assumption required in most random matrix theory settings.

In the context of functional data, variables of interest become the eigenvalues of the covariance operator. However, the rough bound obtained by Weyl’s inequality, $|\hat{\lambda}_{k}-\lambda_{k}|\leq\|\Delta\|_{\mathrm{HS}}$, is suboptimal from two respects. First, $|\hat{\lambda}_{k}-\lambda_{k}|$ should have a degenerated kernel smoothing rate with variance $n^{-1/2}$, whereas $\|\Delta\|_{\mathrm{HS}}$ has a two-dimensional kernel smoothing rate with variance $(n\bar{N}_{2}^{2}h^{2})^{-1/2}$. Second, due to the infinite dimensionality of functional data, the eigenvalues $\{\lambda_{k}\}_{k}$ tend to zero as $k\rightarrow\infty$, so a general bound for all eigenvalues provides little information for those with larger indices. Although expansions and asymptotic normality have been studied for a fixed number of eigenvalues, as well as for those with diverging indices for fully observed functional data, the study of eigenvalues with a diverging index under the discrete sampling scheme remains limited.

In light of the aforementioned issues, we employ the perturbation technique outlined in previous sections to establish the asymptotic normality of eigenvalues with diverging indices, which holds broad implications for inference problems in functional data analysis. Before presenting our results, we introduce the following assumption, which is standard in FPCA for establishing asymptotic normality.

Since $\lambda_{j}$ approaches zero as $j$ approaches infinity, we need to regularize the eigenvalues so that they can be compared on the same scale of variability. For a fixed $j$, $(\hat{\lambda}_{j}-\lambda_{j})/\lambda_{j}$ is $\sqrt{n}$-consistent when $h$ is small. For diverging $j$, Corollary 4 presents three different types of asymptotic normalities, depending on the value of $P_{2}(N)$.