Factor-guided functional PCA
for high-dimensional functional data

Suppose that the observations ${\bf X}_{i}(t),i=1,\cdots,n$ are independent identically distributed (iid), with ${\bf X}_{i}(t)=(X_{i1}(t),\cdots,X_{ip}(t))^{\prime}$ being the high-dimensional functional variables for $p\gg n$. For simplicity, we assume that the mean of $X_{ij}(t)$ has already been removed; namely, $\mathbb{E}\big{(}X_{ij}(t)\big{)}=0$ for any $j$ and $t\in[0,\ 1]$. To reflect the situation of irregular and possibly subject-specific time points, we assume that ${\bf X}_{i}(\cdot)$ are measured at ${\bf t}_{i}=(t_{i1},\cdots,t_{in_{i}})^{\prime}$.

Illustrated by the idea of the factor analysis (Bai and Ng,, 2002; Bai,, 2003; Bai and Ng,, 2006, 2013; Bai and Liao,, 2013), we assume the ${\bf X}_{i}(t)$ are correlated because they share a vector of factor processes ${\bf h}_{i}(t)$. The following model is considered:

$\displaystyle{\bf X}_{i}(t)={\bf B}{\bf h}_{i}(t)+{\bf u}_{i}(t),$

(1)

where ${\bf h}_{i}(t)=(h_{i1}(t),\cdots,h_{iq}(t))^{\prime}$ is a $q$-dimensional vector of factor processes with $q\ll p$, ${\bf B}=({\bf b}_{1},\cdots,{\bf b}_{p})^{\prime}=(b_{jk})_{p\times q}$ is a deterministic matrix, and ${\bf u}_{i}(t)=(u_{i1}(t),\cdots,u_{ip}(t))^{\prime}$ denotes measurement error independent of ${\bf h}_{i}(\cdot)$ with $\mathbb{E}(u_{ij}(t))=0$ and $\hbox{var}(u_{ij}(t))=\sigma^{2}$. Furthermore, to consider the dependence of functional data over time, we apply the Kauhunen-Lo$\rm\grave{e}$ve expansion (Ash and Gardner,, 1975) to the factor process $h_{ij}(t)$, and have,

$\displaystyle h_{ij}(t)=\sum_{k=1}^{\infty}\xi_{ijk}\phi_{jk}(t),$

(2)

where $\phi_{jk}(t)$ is the $k$th orthonormal eigenfunction of the covariance function $C_{j}(s,t)=\mbox{cov}(h_{ij}(s),h_{ij}(t))$ for factor process $j$, which satisfies $\int\phi_{jk}(t)\phi_{jk^{\prime}}(t)dt=1$ if $k=k^{\prime}$, and it is $0$ otherwise; $\xi_{ijk},k=1,2,\cdots$ are the functional principal component scores for the stochastic process $h_{ij}(t)$ with $\mathbb{E}(\xi_{ijk})=0$, $\hbox{var}(\xi_{ijk})=\rho_{jk}$ and $\mathrm{cov}(\xi_{ijk},\xi_{ijk^{\prime}})=0$ if $k\neq k^{\prime}$. $\rho_{jk}$ is the eigenvalue corresponding to the eigenfunction $\phi_{jk}(\cdot)$, with the constraint of $\rho_{j1}\geq\rho_{j2}\geq\cdots>0$ and $\sum_{k=1}^{\infty}\rho_{jk}<\infty$ for any $j=1,\cdots,q$, which implies that $\sup_{t\in[0,1]}\mathbb{E}\big{(}\sum_{k=1}^{K}\xi_{ijk}\phi_{jk}(t)-\sum_{k=1}^{\infty}\xi_{ijk}\phi_{jk}(t)\big{)}^{2}\to 0$ as $K\rightarrow\infty$, we hence suppose that

$\displaystyle h_{ij}(t)=\sum_{k=1}^{K}\xi_{ijk}\phi_{jk}(t).$

(3)

Model (3) has been extensively studied in the literature of FPCA for the case when $K$ is fixed (James et al.,, 2000; Yao et al.,, 2005; Hall et al.,, 2008; Zhou et al.,, 2018). To improve flexibility, we allow $K\rightarrow\infty$ (Hall and Hosseini-Nasab,, 2006; Kong et al.,, 2016; Lin et al.,, 2018). In addition, $K$ can vary with $j$, and we use the same $K$ for $j=1,\cdots,q$ for simple notation.

Let $\bm{\zeta}_{ij}=(\xi_{ij1},\cdots,\xi_{ijK})^{\prime}$, $\bm{\zeta}_{i}=(\bm{\zeta}^{\prime}_{i1},\cdots,\bm{\zeta}^{\prime}_{iq})^{\prime}$ and ${\bm{\Phi}}(t)=\mbox{diag}({\bm{\Phi}}_{1}(t),\cdots,{\bm{\Phi}}_{q}(t))$, where ${\bm{\Phi}}(t)$ is a $Kq\times q$ block diagonal matrix with block $j$ being ${\bm{\Phi}}_{j}(t)=(\phi_{j1}(t),\cdots,\phi_{jK}(t))^{\prime}$. Combined with the expression in (3), the model in (1) can be written as the following factor-guided functional principal component analysis (FaFPCA) model:

$\displaystyle{\bf X}_{i}(t)={\bf B}{\bm{\Phi}}^{\prime}(t)\bm{\zeta}_{i}+{\bf u}_{i}(t).$

(4)

As mentioned above, model (1) highlights the correlation among the functional variables, and the model (3) expresses the temporal dependence of the functional data. As a result, the proposed FaFPCA model in (4) simultaneously organizes the two kinds of dependence structures so that a simple and sufficient expression for the high-dimensional functional may be obtained. Different from the framework of uFPCA (uFPCA-HD), the proposed FaFPCA utilizes factor processes to extract features from functional variables, hence can efficiently avoid overlapping features and organize the correlation information among functional variables. In addition, compared with the framework of factor models (FM-HD), FaFPCA captures the information of temporal dependence by imposing a low-rank decomposition ${\bf B}{\bm{\Phi}}^{\prime}(t)$ on the loading function in FM-HD, and hence is more structured and efficient than FM-HD for dimension reduction.

Obviously, model (4) is not identifiable. We impose three constraints for identifiability:

(I1)

$p^{-1}{\bf B}^{\prime}{\bf B}={\bf I}_{q}$ and the first nonzero element of each column of ${\bf B}$ is positive.

(I2)

$\bm{\zeta}^{\prime}\bm{\zeta}$ is diagonal with decreasing diagonal entries, where $\bm{\zeta}=(\bm{\zeta}_{1},\cdots,\bm{\zeta}_{n})^{\prime}$.

(I3)

$\int{\bm{\Phi}}(t){\bm{\Phi}}(t)^{\prime}dt={\bf I}_{Kq}$ and $\phi_{jk}(0)>0$.

Conditions (I1) and (I2), which are restrictions on loadings and factors, are commonly used in the literature of factor analysis (Bai and Ng,, 2013; Fan et al.,, 2013; Li et al.,, 2018; Jiang et al.,, 2019; Liu et al.,, 2022). Condition (I3), which restricts the eigenfunctions, is widely used in the literature of FPCA (James et al.,, 2000; Yao et al.,, 2005; Hall et al.,, 2008; Zhou et al.,, 2018). We state the identifiability in the following proposition, and its proof is detailed in the Supplementary Materials B.

The estimation involves three terms ${\bf B}$, $\bm{\zeta}$ and ${\bm{\Phi}}(\cdot)$. Although model (4) is a kind of factor model, the existing methods for factor analysis are not feasible for the estimation of model (4) because the matrix-structured data required for factor analysis is not available for functional data due to their individual-specific observation times. A possible method is to iterate among the three terms ${\bf B}$, $\bm{\zeta}$ and ${\bm{\Phi}}(\cdot)$, which is considerably computationally intensive because all of them are infinite. Moreover, the iteration usually fails to converge without appropriate initial values, which are almost impossible to obtaine in the current situation. Here, we avoid the iteration method and solve the computational problem by building some moment equations. Concretely, we denote the covariance matrices of ${\bf X}_{i}(t)$ and ${\bf u}_{i}(t)$ by ${\bf\Sigma}_{{\bf X}}(t)$ and ${\bf\Sigma}_{{\bf u}}(t)$, respectively. Let $\widetilde{\bf\Sigma}_{{\bf X}}=\int{\bf\Sigma}_{{\bf X}}(t)dt$, $\widetilde{\bf\Sigma}_{{\bf u}}=\int{\bf\Sigma}_{{\bf u}}(t)dt$ and ${\bf\Lambda}_{\bm{\zeta}}=\mbox{diag}\big{(}\sum_{k=1}^{K}\hbox{var}(\xi_{i1k}),\cdots,\sum_{k=1}^{K}\hbox{var}(\xi_{iqk})\big{)}$. By calculating the covariance matrix and integrating over $t$ on both sides of (4), we have

$\displaystyle\widetilde{\bf\Sigma}_{{\bf X}}={\bf B}{\bf\Lambda}_{\bm{\zeta}}{\bf B}^{\prime}+\widetilde{\bf\Sigma}_{{\bf u}},$

(5)

which implies that $p^{-1/2}{\bf B}$ is the eigenvector of $\widetilde{\bf\Sigma}_{{\bf X}}$. Considering the identification condition $p^{-1}{\bf B}^{\prime}{\bf B}={\bf I}_{q}$, we estimate $p^{-1/2}\widehat{\bf B}$ for $p^{-1/2}{\bf B}$ by the orthogonal eigenvectors of the numerical approximation version of $\widetilde{\bf\Sigma}_{{\bf X}}$; i.e.,

$\displaystyle\widehat{\bf B}=p^{1/2}{\bf E}_{eigen}\Big{(}n^{-1}\sum_{i=1}^{n}n^{-1}_{i}\sum_{l=1}^{n_{i}}{\bf X}_{i}(t_{il}){\bf X}^{\prime}_{i}(t_{il});q\Big{)},$

(6)

where ${\bf E}_{eigen}(A;q)$ is a matrix composed of the orthogonal eigenvectors corresponding to the $q$ largest eigenvalues of matrix $A$.

Now we consider ${\bm{\Phi}}$ and $\bm{\zeta}$. Without loss of generality, we assume that the support of $t_{ij}$ is $[0,\ 1]$ by proper scaling. ${\bf M}(\cdot)=(M_{1}(\cdot),\ \cdots,\ M_{\tau_{n}}(\cdot))^{\prime}$ denote a vector of B-spline basis functions on $[0,\ 1]$ with $\tau_{n}=O(n^{v})$; then we have $\phi_{jk}(t)\approx\bm{\Theta}^{\prime}_{jk}{\bf M}(t)$. Let $\bm{\Theta}_{j}=(\bm{\Theta}_{j1},\cdots,\bm{\Theta}_{jK})^{\prime}\in\mathbb{R}^{K\times\tau_{n}}.$ With the B-spline approximation, identification condition (I2) on ${\bm{\Phi}}(\cdot)$ can be written as $\tau^{-1}_{n}\bm{\Theta}_{j}\bm{\Theta}^{\prime}_{j}={\bf I}_{K},j=1,\cdots,q$ and $\tau_{n}\int{\bf M}(t){\bf M}^{\prime}(t)dt={\bf I}_{\tau_{n}}$. Then by (3), we have

$\displaystyle h_{ik}(t_{il})={\bm{\Phi}}^{\prime}_{k}(t_{il})\bm{\zeta}_{ik}\approx{\bf M}^{\prime}(t_{il})\bm{\Theta}^{\prime}_{k}\bm{\zeta}_{ik}.$

(7)

On the other hand, by multiplying $p^{-1}{\bf B}^{\prime}$ on both sides of (1), we have $p^{-1}{\bf B}^{\prime}{\bf X}_{i}(t)\approx p^{-1}{\bf B}^{\prime}{\bf B}{\bf h}_{i}(t).$ Since $p^{-1}{\bf B}^{\prime}{\bf B}={\bf I}_{q}$, we obtain

$\displaystyle h_{ik}(t_{il})\approx p^{-1}\sum_{j=1}^{p}b_{jk}X_{ij}(t_{il}).$

(8)

Combining (7) and (8), we obtain,

$\displaystyle p^{-1}\sum_{j=1}^{p}b_{jk}X_{ij}(t_{il})\approx{\bf M}^{\prime}(t_{il})\bm{\Theta}^{\prime}_{k}\bm{\zeta}_{ik}.$

(9)

Then, multiplyling ${\bf M}(t_{il})$ on both sides of (9) and making the summation over the observation times of individual $i$, we have

$$\big{\{}\sum_{l=1}^{n_{i}}{\bf M}(t_{il}){\bf M}^{\prime}(t_{il})\big{\}}^{-1}\big{\{}\sum_{l=1}^{n_{i}}\sum_{j=1}^{p}p^{-1}{\bf M}(t_{il})b_{jk}X_{ij}(t_{il})\big{\}}\approx\bm{\Theta}^{\prime}_{k}\bm{\zeta}_{ik},$$

which is a factor model with response ${\bf w}_{ik}=\big{\{}\sum_{l=1}^{n_{i}}{\bf M}(t_{il}){\bf M}^{\prime}(t_{il})\big{\}}^{-1}\big{\{}\sum_{l=1}^{n_{i}}\sum_{j=1}^{p}p^{-1}{\bf M}(t_{il})\widehat{b}_{jk}X_{ij}(t_{il})\big{\}}$, factor $\bm{\zeta}_{ik}$ and loading $\bm{\Theta}_{k}$, where $\hat{b}_{jk}$ is the estimator from (6). Here ${\bf W}_{k}=({\bf w}_{1k},\cdots,{\bf w}_{nk})^{\prime}\in\mathbb{R}^{n\times\tau_{n}}$ and $\bm{\zeta}_{[k]}=(\bm{\zeta}_{1k},\cdots,\bm{\zeta}_{nk})^{\prime}\in\mathbb{R}^{n\times K}$; hence, we estimate $(\bm{\zeta}_{[k]},\bm{\Theta}_{k})$ by

$\displaystyle(\widehat{\bm{\zeta}}_{[k]},\widehat{\bm{\Theta}}_{k})=\mathop{\arg\min}_{\bm{\zeta}_{[k]},\bm{\Theta}_{k}}\Big{\|}{\bf W}_{k}-\bm{\zeta}_{[k]}\bm{\Theta}_{k}\Big{\|}^{2}_{F},$

(10)

with $\|\cdot\|_{F}$ being the Frobenius-norm of a matrix. Following Bai and Ng, (2013), which states that the factors can be estimated asymptotically by principal component analysis (PCA), we obtain estimators for $\bm{\Theta}_{k}$ and $\bm{\zeta}_{[k]}$, respectively by

	$\displaystyle\widehat{\bm{\Theta}}_{k}^{\prime}={\bf E}_{eigen}({\bf W}_{k}^{\prime}{\bf W}_{k};K),$		(11)
	$\displaystyle\widehat{\bm{\zeta}}_{[k]}=\tau_{n}^{-1/2}{\bf W}_{k}\times{\bf E}_{eigen}({\bf W}_{k}^{\prime}{\bf W}_{k};K),$		(12)

for $k=1,\cdots,q$. Finally, we estimate $\phi_{jk}(t)$ by $\widehat{\phi}_{jk}(t)=\widehat{\bm{\Theta}}^{\prime}_{jk}{\bf M}(t)$.

The estimators (6), (11) and (12) for ${\bf B}$, $\bm{\Theta}$ and $\bm{\zeta}$, respectively, have closed forms based on eigenvalue decomposition without any rotation, and hence, the computation and programming are very simple.

We need to determine the dimension of the factor components $q$ and the number of eigenfunctions $K$. Since we transform the estimation for ${\bf B}$, $\bm{\zeta}$ and ${\bm{\Phi}}(\cdot)$ to the problem of principal component analysis, $q$ and $K$ can be selected by calculating the proportion of variability explained by the principal components (James et al.,, 2000; Happ and Greven,, 2018). Particularly, we choose $q$ and $K$ such that

$$\sum_{j=1}^{q}\lambda_{j}\left(\sum_{i=1}^{n}n^{-1}_{i}\sum_{l=1}^{n_{i}}{\bf X}_{i}(t_{il}){\bf X}_{i}^{\prime}(t_{il})\right)/\sum_{j=1}^{p}\lambda_{j}\left(\sum_{i=1}^{n}n^{-1}_{i}\sum_{l=1}^{n_{i}}{\bf X}_{i}(t_{il}){\bf X}_{i}^{\prime}(t_{il})\right)>95\%,$$

$$\min_{k\in\{1,\cdots,q\}}\sum_{j=1}^{K}\lambda_{j}({\bf W}^{\prime}_{k}{\bf W}_{k})/\sum_{j=1}^{\tau_{n}}\lambda_{j}({\bf W}^{\prime}_{k}{\bf W}_{k})>95\%,$$

respectively, where $\lambda_{k}(A)$ is the $k$th eigenvalue of $A$. In practice, we can choose a different $K_{k}$ for $k=1,\cdots,q$ so that $\sum_{j=1}^{K_{k}}\lambda_{j}({\bf W}^{\prime}_{k}{\bf W}_{k})/\sum_{j=1}^{\tau_{n}}\lambda_{j}({\bf W}^{\prime}_{k}{\bf W}_{k})>95\%$ is more flexible.

For each trajectory ${\bf X}_{i}(t)$, $20$ observation time points are randomly sampled from $U(0,10)$ and 200 replicates are applied in each scenario unless otherwise stated.

Scenario 1: To be fair, we generate data from the following model, which does not follow the assumptions of the methods we consider,

$${\bf X}_{i}(t)=\sum^{K}_{k=1}\bm{\zeta}_{ik}\phi_{k}(t)+\bm{\epsilon_{i}}(t),\ i=1,\cdots,n,$$

where $\bm{\zeta}_{ik}\in\mathbb{R}^{p}\sim N({\bf 0},{\bf\Sigma}_{\bm{\zeta},k})$, $\bm{\epsilon}_{i}(t)\sim N({\bf 0},{\bf\Sigma})$ for each $t$, $({\bf\Sigma}_{\bm{\zeta},k})_{ij}=0.5^{|i-j|}$ if $i\neq j$ and $({\bf\Sigma}_{\bm{\zeta},k})_{ii}=3-(k-1)/(K-1)$ for each $i=1,\cdots,p$ and $k=1,\cdots,K$. We take ${\bf\Sigma}=(\sigma_{ij})_{p\times p}$ as ${\bf Case\ I}:{\bf\Sigma}={\bf I}$ and ${\bf Case\ II}:\sigma_{ii}=1$ and $\sigma_{ij}=0.3$ if $i\neq j$, and let $\phi_{k}(t)=\sin[(2k-1)\pi t/10]$, $n=p=100$ and $K=2,10$.

Scenario 2: To investigate the performance of the estimation for loadings, factors and eigenfunctions, we generate functional data by ${\bf X}_{i}(t)={\bf B}{\bf h}_{i}(t)+{\bf u}_{i}(t),\ i=1,\cdots,n,$ which satisfie the assumption of the proposed method, where random error ${\bf u}_{i}(t)\sim N({\bf 0},{\bf I})$. To construct ${\bf B}$, we first generate $n$ samples of the $p$-dimensional ${\bf k}_{i}$ from $N(0,(\sigma_{ij})_{p\times p})$ with $\sigma_{ij}=0.5^{|i-j|}$. Here, ${\bf K}=({\bf k}_{1},\cdots,{\bf K})^{\prime}.$ We apply eigendecomposition to the matrix ${\bf K}{\bf K}^{\prime}$ to obtain ${\bf B}_{q}$ and ${\bf K}_{q}$, where ${\bf B}_{q}$ is a diagonal matrix whose diagonal entries are the first $q$ largest eigenvalues of ${\bf K}{\bf K}^{\prime}$, and ${\bf K}_{q}\in\mathbb{R}^{n\times q}$ consists of the corresponding eigenvectors. ${\bf B}_{n}={\bf K}^{\prime}{\bf K}_{q}$, and QR decomposition of ${\bf B}_{n}={\bf Q}_{n}{\bf R}_{n}$ is performed. Then, ${\bf B}=p^{1/2}{\bf Q}_{n}$ so that the identification condition on ${\bf B}$ can be satisfied. Similarly, to construct $\bm{\zeta}_{i}$ for ${\bf h}_{i}(t)$, we independently generate $\xi^{*}_{ik}\sim N(0,Kqk^{-1})$, denote $\bm{\zeta}^{*}_{i}=(\xi^{*}_{i1},\cdots,\xi^{*}_{i,Kq})^{\prime}$ and $\bm{\zeta}^{*}=(\bm{\zeta}^{*}_{1},\cdots,\bm{\zeta}^{*}_{n})$, and then perform eigen-decomposition of the matrix $\bm{\zeta}^{*\prime}\bm{\zeta}^{*}$ for matrices ${\bf\Lambda}_{Kq}$ and ${\bf M}_{Kq}$, which consist of the first $Kq$ largest eigenvalues and the corresponding eigenvectors, respectively. We take $\bm{\zeta}\hat{=}(\bm{\zeta}_{1},\cdots,\bm{\zeta}_{n})^{\prime}={\bf\Lambda}^{1/2}_{Kq}{\bf M}^{\prime}_{Kq}.$ It is easy to check that $\bm{\zeta}$ satisfies Identification Condition (I2), as described in Section 2.1. Given $\bm{\zeta}_{i},$ the factor process ${\bf h}_{i}(t)=(h_{i1}(t),\cdots,h_{iq}(t))^{\prime}$ is generated by $h_{ij}(t)=\sum_{k=1}^{K}\xi_{ijk}\phi_{jk}(t)$, where $\phi_{jk}(t)=\sqrt{2}\sin[\{2(j\neq 2)+(j=2)\}k\pi t/10]$ if $k$ is odd and $\phi_{jk}(t)=\sqrt{2}\cos[\{2(j\neq 2)k+(j=2)(2k+1)\}\pi t/10]$ if $k$ is even. The setting is used to ensure the orthogonality of eigenfunctions.

Scenario 1: We first compare FaFPCA with the existing methods in terms of prediction accuracy PE. Figure 1 shows the PE of the proposed FaFPCA, FM-HD, DPCA, sFPCA-HD and uPFCA-HD methods with the corresponding optimal tuning parameters under the cases listed in Scenario 1. From Figure 1, we can see that (1) When ${\bf u}_{i}(t)$ becomes correlated, that is, when ${\bf\Sigma}$ changes from ${\bf I}$ to ${\bf\Sigma}\neq{\bf I}$ in Figures 1(a) and 1(b), respectively, FaFPCA, FM-HD, and DPCA perform better, while the PE performance of sFPCA-HD and uPFCA-HD worsens. This indicates that FaFPCA, FM-HD and DPCA are efficient in making use of the correlations among functional variables while sFPCA-HD and uPFCA-HD are not. (2) On the other hand, when $K$ changes from 10 to 2 in Figures 1(c) and 1(d), which means the dependence over time gets stronger, sFPCA-HD and uPFCA-HD perform better, while FM-HD and DPCA perform worse. This is not surprising because sFPCA-HD and uPFCA-HD focus on dependence over time while FM-HD and DPCA do not. (3) The proposed FaFPCA approach performs the best and is stable in all the cases in terms of PE; particularly, either the dependence over time or correlations among variables increase, the accuracy of FaFPCA is improved. Compared to the other methods, the obtained by FaFPCA is due to the usage of the two kind of correlations.

We carried out the computation on a single 4-core personal laptop with 16 GB of RAM. Table 1 shows the average running time of 200 repetitions for Case I with $K=2$. Table 1 shows that the proposed FaFPCA approach is much faster than the other methods except FM-HD, which only considers model (1) and does not estimate eigenfunctions and factor scores. uFPCA-HD requires an extremely large amount of time because uFPCA is applied to each of the $p$ curves separately.

Scenario 2: Figure 2 displays RMSE_l for loadings, RMSE_f for factors and RMSE_e for eigenfunctions using the proposed FaFPCA method for Scenario 2 with $(q,K)=(5,2)$ when $(n,p)$ varies. The results in Figure 2 indicate that the precision of ${\bf B},\bm{\zeta}$ and ${\bm{\Phi}}(t)$ increases as $n$ or $p$ increases, which is consistent with the theoretical results. In particular, the precision of $\bm{\zeta}$ is more sensitive to $p$ than to $n$, while that of ${\bf B}$ is more sensitive to $n$. This is because the estimator for $\bm{\zeta}_{i}$ is based on $p$ functions of individual $i$; however, the estimator for ${\bf B}$ is mainly based on $n$ individuals.

To show the effect of the number of observations, we take $n_{i}$ from 5 to 50 given $n=p=100$ and $(q,K)=(5,2)$. The results in Figure 3 show that the estimators are quite stable with respect to $n_{i}$ except for small $n_{i}$ values. When $n_{i}$ is small, the matrix $\sum_{l=1}^{n_{i}}{\bf M}(t_{il}){\bf M}^{\prime}(t_{il})$ in ${\bf w}_{ik}$ might be irreversible. To avoid this problem, we add a term $\lambda{\bf I}_{\tau_{n}}$, where $\lambda$ is small, to $\sum_{l=1}^{n_{i}}{\bf M}(t_{il}){\bf M}^{\prime}(t_{il})$ to ensure invertibility. This may cause extra bias for the estimations of factors and eigenfunctions when $n_{i}$ is small. From Figure 3, we also observe an interesting phenomenon as $n_{i}$ increases, RMSE_l for loadings increases, while RMSE_f and RMSE_e for factors and eigenfunctions, respectively, decrease. This can be attributed the following two aspects. On the one hand, the loadings are estimated mainly based on $n$ samples, rather than $n_{i}$ or $p$; hence, the increasing $n_{i}$ values may introduce extra noise due to the integration of all observations using (6). On the other hand, the factors $\bm{\zeta}_{i}$ are estimated mainly based on the observations from individual $i$; hence, the estimation of the factors becomes better as $n_{i}$ increases, although becomes stable when $n_{i}$ is large enough, for example, when $n_{i}>10$ in the cases. The estimation of eigenfunctions has a similar performance as that of the factors since they are simultaneously estimated.