Nonlinear Functional Principal Component Analysis Using Neural Networks

Let $X(t)$ be a smooth random function in $L^{2}(\mathcal{T})$, where $\mathcal{T}$ is a bounded and closed interval. In this paper, $\mathcal{T}$ is set as $[0,1]$ if there is no specific explanation. Let $\mu(t)$ and $\Sigma(s,t)$ denote the mean function and covariance function of $X(t)$, respectively. For linear FPCA, according to the Karhunen-Loève Theorem, $X(t)$ admits the following expansion

$\displaystyle X(t)=\mu(t)+\sum_{k=1}^{\infty}\xi_{k}\phi_{k}(t),$

where $\xi_{k}$ is the $k$-th functional principal component score, $\phi_{k}(t)$ is the $k$-th eigenfunction of $\Sigma(s,t)$ and satisfies $\int_{\mathcal{T}}\phi_{k}(t)^{2}dt=1$ and $\int_{\mathcal{T}}\phi_{k}(t)\phi_{l}(t)dt=0$ for $l\neq k$. Moreover, the functional principal component scores are uncorrelated random variables with mean zero and $E\xi_{k}^{2}=\lambda_{k}$, where $\lambda_{k}$ is the $k$-th eigenvalue of $\Sigma(s,t)$. In practice, as only the first several functional principal component scores dominate the variation, a truncated expansion

$\displaystyle X(t)=\mu(t)+\sum_{k=1}^{K}\xi_{k}\phi_{k}(t)$

(1)

is often applied, where $K$ is the number of functional principal components that used. Furthermore, we also have that

$\displaystyle\xi_{k}=\int_{\mathcal{T}}\{X(t)-\mu(t)\}\phi_{k}(t)dt.$

(2)

It is obvious that $X(t)$ is mapped into a lower dimensional space through a linear transformation. However, with the constraint of linear map, the nonlinear structure is ignored, which may lower the efficiency of dimension reduction.

For nonlinear FPCA, we extend the linear map in (2) to arbitrary nonlinear map. That is

$\displaystyle\xi_{k}=G_{k}(X),$

(3)

where $G_{k}$ is a nonlinear function that maps function in square integrable space $L^{2}(\mathcal{T})$ to scalar in $\mathbb{R}$. Similarly, (1) can be generalized into nonlinear version as

$\displaystyle X(t)=H(\bm{\xi})(t),$

where $\bm{\xi}=(\xi_{1},\ldots,\xi_{K})^{\top}$ and $H$ is a nonlinear function that maps from a vector space to a square integrable space. The scores obtained from nonlinear FPCA also contain the main information of $X(t)$, but the dimension of the feature space can be lower, as the nonlinear structure is taken into account in the process of dimension reduction. To estimate the nonlinear functional principal component scores, the nonlinear functions $G_{k}$ and $H$ have to be learnt, and a neural network method is employed.

In this section, we construct a functional autoassociative neural network for nonlinear FPCA. The structure of the proposed neural network is shown in Figure 1. The output is the reconstruction of the input data, thus both input and output are functions in our network, which brings challenge to the optimization of the neural network. Furthermore, dimension reduction is realized by the second hidden layer, which can also be called “bottleneck” layer as in (Kramer, 1991). More details about the computation of the proposed functional autoassociative neural network are given below.

To be specific, the left two hidden layers are designed to learn the nonlinear functions $G_{k}$’s that map the functional input to the scores, which can be viewed as a dimension reduction process. For the $j$th neuron in the first hidden layer, we define that

$\displaystyle H_{j}^{(1)}=\sigma\Big{\{}b_{j}+\int_{\mathcal{T}}\beta_{j}(t)X(t)dt\Big{\}},\ j=1,\ldots,J,$

where $b_{j}\in\mathbb{R}$ is the intercept, $\beta_{j}(t)\in L^{2}(\mathcal{T})$ is the weight function, $\sigma(\cdot)$ is a nonlinear activation function, and $J$ is the number of neurons in the first hidden layer. This is a natural generalization of the first two layers of a multilayer perceptron to adapt to the functional input. As $H_{j}^{(1)}\in\mathbb{R}$, computation of the score in the second hidden layer can be promoted naturally. The $k$th neuron in the second hidden layer is

$\displaystyle H_{k}^{(2)}=\sum_{j=1}^{J}w_{jk}H_{j}^{(1)},\ k=1,\ldots,K,$

where $w_{jk}\in\mathbb{R}$ is the weight. Moreover, let $\mathcal{S}(\sigma,L^{2}(\mathcal{T}))$ be the set of functions from $L^{2}(\mathcal{T})$ to $\mathbb{R}$ of the form

$\displaystyle X\mapsto\sum_{j=1}^{J}w_{jk}\sigma\Big{\{}b_{j}+\int_{\mathcal{T}}\beta_{j}(t)X(t)dt\Big{\}}.$

According to Corollary 5.1.2 in (Stinchcombe, 1999) and the proof in Section 6.1.2 of (Rossi and Conan-Guez, 2006), $\mathcal{S}(\sigma,L^{2}(\mathcal{T}))$ has the universal approximation property for $L^{2}(\mathcal{T})$. That means any nonlinear function from $L^{2}(\mathcal{T})$ to $\mathbb{R}$ can be approximated up to arbitrary degree of precision by functions in $\mathcal{S}(\sigma,L^{2}(\mathcal{T}))$.

The right two layers in Figure 1, which correspond to the estimation of the nonlinear function $H$, are used to reconstruct $X(t)$ by the low-dimensional scores $H_{k}^{(2)}$ in the second hidden layer. The procedure can be more challenging, since we have to get functional output from scalars in the second hidden layer. To this end, the $r$th neuron in the third layer is defined as

$\displaystyle H_{r}^{(3)}(t)=\sigma\Big{\{}a_{r}(t)+\sum_{k=1}^{K}\gamma_{kr}(t)H_{k}^{(2)}\Big{\}},\ t\in\mathcal{T},\ r=1,\ldots,R,$

where $a_{r}(t)\in L^{2}(\mathcal{T})$ is the intercept function, $\gamma_{kr}(t)\in L^{2}(\mathcal{T})$ is the weight function, and $R$ is the number of the neurons in the third hidden layer. It can be observed that each neuron in the third hidden layer is a function. Then, $X(t)$ is reconstructed by

$\displaystyle\widehat{X}(t)=\sum_{r=1}^{R}u_{r}H_{r}^{(3)}(t),\ t\in\mathcal{T},$

(4)

where $u_{r}\in\mathbb{R}$ is the weight. The whole network is trained by minimizing the following reconstruction error

$\displaystyle RE=\int_{\mathcal{T}}\{X(t)-\widehat{X}(t)\}^{2}dt.$

As functional data is involved in the proposed network, some of the parameters need to be estimated are functions, which makes the optimization of the network nontrivial. In Section 2.3, we introduce the optimization algorithm for practical implementation.

Note that $H_{k}^{(2)}$ can be viewed as the estimation of $\xi_{k}$ in (3). Therefore, the dimension of $X$ is reduced to $K$ through the functional autoassociative neural network. We can use the low-dimensional vector to complete further inference, such as curve reconstruction, regression and clustering. In this paper, we mainly focus on the curve reconstruction by the low-dimensional representation.

As discussed in Section 2.2, it can be hard to optimize the proposed functional autoassociative neural network in Figure 1, since many parameters appear as a function. Here, we employ the B-spline basis functions to transform the estimation of functions to the estimation of their coefficients. Let $\textbf{B}_{L}=\{B_{l}(t),l=1,\ldots,L\}$ be a set of B-spline basis functions with degree $d$, where $L$ is the number of basis functions. Then, we have

$\displaystyle\beta_{j}(t)=\sum_{l=1}^{L}c_{jl}B_{l}(t),\ a_{r}(t)=\sum_{l=1}^{L}\alpha_{rl}B_{l}(t),\ \gamma_{kr}(t)=\sum_{l=1}^{L}v_{krl}B_{l}(t),$

for $j=1,\ldots,J$, $k=1,\ldots,K$, and $r=1,\ldots,R$, where $c_{jl}$, $\alpha_{rl}$ and $v_{krl}$ are the basis expansion coefficients of $\beta_{j}(t)$, $a_{r}(t)$ and $\gamma_{kr}(t)$, respectively.

With the use of B-spline basis functions, the computation of the first two layers for the proposed functional autoassociative neural network turns out to be

$\displaystyle H_{k}^{(2)}=\sum_{j=1}^{J}w_{jk}\sigma\Big{\{}b_{j}+\sum_{l=1}^{L}c_{jl}\int_{\mathcal{T}}X(t)B_{l}(t)dt\Big{\}}\triangleq\sum_{j=1}^{J}w_{jk}\sigma\Big{(}b_{j}+\sum_{l=1}^{L}c_{jl}\widetilde{X}_{l}\Big{)},$

(5)

where $\widetilde{X}_{l}=\int_{\mathcal{T}}X(t)B_{l}(t)dt$. In practice, $\widetilde{X}_{l}$ is calculated through the B-spline expansion of $X$, that is $\widetilde{X}_{l}=\sum_{h=1}^{L}x_{h}\{\int_{\mathcal{T}}B_{h}(t)B_{l}(t)dt\}$, where $x_{h}$ is the basis expansion coefficient of $X(t)$. Then, we have

$\displaystyle H_{k}^{(2)}=\sum_{j=1}^{J}w_{jk}\sigma\Big{[}b_{j}+\sum_{h=1}^{L}\Big{\{}\sum_{l=1}^{L}c_{jl}\int_{\mathcal{T}}B_{h}(t)B_{l}(t)dt\Big{\}}x_{h}\Big{]}\triangleq\sum_{j=1}^{J}w_{jk}\sigma\Big{(}b_{j}+\sum_{h=1}^{L}d_{jh}x_{h}\Big{)},$

where $d_{jh}=\sum_{l=1}^{L}c_{jl}\int_{\mathcal{T}}B_{h}(t)B_{l}(t)dt$. The following theorem discusses the universal approximation property of the first two layers based on B-spline basis functions. The proof is provided in the Appendix.

For the computation of the third hidden layer of the proposed functional autoassociative neural network, we have

	$\displaystyle a_{r}(t)+\sum_{k=1}^{K}\gamma_{kr}(t)H_{k}^{(2)}$	$\displaystyle=\sum_{l=1}^{L}\alpha_{rl}B_{l}(t)+\sum_{k=1}^{K}\sum_{l=1}^{L}v_{krl}B_{l}(t)H_{k}^{(2)}=A_{r}^{\top}B(t)+\sum_{k=1}^{K}V_{kr}^{\top}B(t)H_{k}^{(2)}$
		$\displaystyle=\Big{(}A_{r}+\sum_{k=1}^{K}V_{kr}H_{k}^{(2)}\Big{)}^{\top}B(t),$

where $A_{r}=(\alpha_{r1},\ldots,\alpha_{rL})^{\top}$, $V_{kr}=(v_{kr1},\ldots,v_{krL})^{\top}$, and $B(t)=(B_{1}(t),\ldots,B_{L}(t))^{\top}$. For the term $A_{r}+\sum_{k=1}^{K}V_{kr}H_{k}^{(2)}$ above, it can be viewed as weighted sums of $H_{1}^{(2)},\ldots,H_{K}^{(2)}$. Thus, the computation from the second hidden layer to the third hidden layer of the functional autoassociative neural network can be transformed accordingly, that is

$\displaystyle H_{r}^{(3)}(t)=\sigma\Big{\{}\Big{(}A_{r}+\sum_{k=1}^{K}V_{kr}H_{k}^{(2)}\Big{)}^{\top}B(t)\Big{\}},t\in\mathcal{T},$

(6)

the structure of which is shown in Figure 2. The middle layer in Figure 2 represents the elements of the $L$-dimensional vector $A_{r}+\sum_{k=1}^{K}V_{kr}H_{k}^{(2)}$. Furthermore, the reconstruction of $X(t)$ can be obtained by (4) as discussed in Section 2.2.

In practice, suppose that $X_{1},\ldots,X_{n}$ are $n$ independent realizations of $X$, where $n$ is the sample size. Then, the loss function for the network can be represented as

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\int_{\mathcal{T}}\{X_{i}(t)-\widehat{X}_{i}(t)\}^{2}dt.$

However, integral appeared in the loss function can increase the difficulty of optimization. Hence, right-hand Riemann sum is employed for the computation. Specifically, let $0=t_{1}<t_{2}<\ldots<t_{M}=1$ be some equally spaced times points on $\mathcal{T}$. Moreover, denote $s_{1},\ldots,s_{T}$ be the observation time points of the random curves, where $T$ is the observation size. Then, we consider the following loss function

$\displaystyle\widetilde{RE}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{M-1}\sum_{m=2}^{M}\{\widetilde{X}_{i}(t_{m})-\widehat{X}_{i}(t_{m})\}^{2},$

(7)

where $\widetilde{X}_{i}(t_{m})$ is the estimation of $X_{i}(t_{m})$ using the observed data by smoothing. Note that only values of the curves at $t_{1},\ldots,t_{M}$ involved in the loss function. Therefore, we just need to consider discrete values of $X(t)$ in the last two layers of the network. Let $\textbf{B}_{l}=(B_{l}(t_{0}),\ldots,B_{l}(t_{M}))^{\top}$ for $l=1,\ldots,L$, and $\textbf{B}=(\textbf{B}_{1},\ldots,\textbf{B}_{L})^{\top}$. The $r$th neuron of the third hidden layer can be computed by

$\displaystyle\widetilde{H}_{r}^{(3)}=\sigma\Big{\{}\Big{(}A_{r}+\sum_{k=1}^{K}V_{kr}H_{k}^{(2)}\Big{)}^{\top}\textbf{B}\Big{\}}.$

(8)

Then the output is given by

$\displaystyle\widehat{\textbf{X}}=\sum_{r=1}^{R}u_{r}\widetilde{H}_{r}^{(3)},$

(9)

where $\widehat{\textbf{X}}=(\widehat{X}(t_{0}),\ldots,\widehat{X}(t_{M}))$. To be clear, the transformed functional autoassociative neural network is shown in Figure 3. The computation involved in the transformed autoassociative neural network corresponds to (5), (8) and (9) respectively.

By turning the proposed functional autoassociative neural network in Section 2.2 into the transformed functional autoassociative neural network, the Adam algorithm (Kingma and Ba, 2014) can be employed in the optimization. This algorithm is popularly-used, and it can be realized by the Python package torch. Moreover, we also provide the Python package nFunNN for the implementation specific to our method.

In this section, we conduct a simulation study to explore the performance of the proposed nFunNN method. For the simulated data, we take into account the measurement error to make it more consistent with the practical cases. Specifically, the observed data is generated by

$\displaystyle Y_{ij}=X_{i}(s_{j})+\epsilon_{ij},\ i=1,\ldots,n,\ j=1,\ldots,T,$

where $Y_{ij}$ is the $j$th observation for the $i$th subject, and $\epsilon_{ij}$’s are the independent measurement errors, which we obtain from the normal distribution $\mathcal{N}(0,\delta^{2})$. We set $\delta=0.1$, and the observation size is set as $T=51$. The observation grids are equally spaced on $\mathcal{T}=[0,1]$. For the setting of $X_{i}$, we consider the following cases.

• •

  Case 1: $X_{i}(t)=\xi_{i1}\sin(2\pi t)+\xi_{i2}\cos(2\pi t),\ t\in\mathcal{T}$, where $\xi_{i1}$’s and $\xi_{i2}$’s are simulated from $\mathcal{N}(0,3^{2})$ and $\mathcal{N}(0,2^{2})$, respectively.

• •

  Case 2: $X_{i}(t)=\xi_{i2}\sin(\xi_{i1}t),\ t\in\mathcal{T}$, where both $\xi_{i1}$’s and $\xi_{i2}$’s are simulated from $\mathcal{N}(0,2^{2})$.

• •

  Case 3: $X_{i}(t)=\xi_{i2}\cos(\xi_{i1}t),\ t\in\mathcal{T}$, where both $\xi_{i1}$’s and $\xi_{i2}$’s are simulated in the same way as that in Case 2.

• •

  Case 4: $X_{i}(t)=\xi_{i1}\sin(2\pi t)+\xi_{i2}\cos(2\pi t)+\xi_{i2}\sin(\xi_{i1}t),\ t\in\mathcal{T}$, where both $\xi_{i1}$’s and $\xi_{i2}$’s are simulated in the same way as that in Case 2.

• •

  Case 5: $X_{i}(t)=\xi_{i1}\sin(2\pi t)+\xi_{i2}\cos(2\pi t)+\xi_{i2}\cos(\xi_{i1}t),\ t\in\mathcal{T}$, where both $\xi_{i1}$’s and $\xi_{i2}$’s are simulated in the same way as that in Case 2.

The above setups include various structures of $X_{i}(t)$. In Case 1, $X_{i}(t)$ is actually generated through the Karhunen-Loève expansion, with zero mean, $\lambda_{1}=3^{2}$, $\lambda_{2}=2^{2}$, and $\lambda_{k}=0$ for $k\geq 3$. Thus, $X_{i}(t)$ in Case 1 has a linear structure, and a linear method may be suitable enough for this case. Moreover, the other four cases consider the nonlinear structure of $X_{i}(t)$. Case 2 and Case 3 impose only one nonlinear term in the setup, while Case 4 and Case 5 combine the linear terms in Case 1 with nonlinear terms in Case 2 and Case 3 respectively. Further, whether a linear structure or a nonlinear structure is considered, all the five cases are set to contain two principal components.

The proposed nFunNN method is compared with the classical linear FPCA method (Ramsay and Silverman, 2005). Specifically, the numbers of neurons in different layers for our transformed functional autoassociative network are set as $L=10$, $J=20$, $K=2$, and $R=20$. And the number of principal components for the linear FPCA method is selected as $2$. To evaluate the performance of curve reconstruction, we consider the following criteria:

	RMSE	$\displaystyle=\sqrt{\frac{1}{nM}\sum_{i=1}^{n}\sum_{m=1}^{M}\{\widehat{X}_{i}(t_{m})-X_{i}(t_{m})\}^{2}},$
	RRMSE	$\displaystyle=\sqrt{\sum_{i=1}^{n}\sum_{m=1}^{M}\{\widehat{X}_{i}(t_{m})-X_{i}(t_{m})\}^{2}}\Bigg{/}\sqrt{\sum_{i=1}^{n}\sum_{m=1}^{M}X_{i}(t_{m})^{2}},$

where RRMSE means the relative RMSE. Note that the prediction at $t_{1},\ldots,t_{M}$ is assessed, and these time points can be different from the observation time points. In our simulation study, $t_{1},\ldots,t_{M}$ are also equally spaced on $[0,1]$ with $M$ being set as $101$. We consider the performance for both training set and test set in the evaluation, and the sizes of both sets are set as $1000$. Moreover, $100$ Monte Carlo runs are conducted for each considered case.

Table 1 lists the simulation results of the proposed nFunNN method and FPCA method in all the five cases. In Case 1, though FPCA method yields slightly better RMSE and RRMSE than our nFunNN method, both methods perform well for training set and test set. The result is not surprising, since $X_{i}(t)$ is generated by the Karhunen-Loève expansion in Case 1, and classical linear FPCA is good enough to handle such case. For Case 2 and Case 3, where nonlinear structure is considered, it is evident that the proposed nFunNN method outperforms the FPCA method and gives more accurate prediction. That implies the advantage of our nFunNN method when solving nonlinear cases. Furthermore, FPCA method is almost invalid for Case 4 and Case 5, while the proposed nFunNN method still provides encouraging curve reconstruction results. In the setting of $X_{i}(t)$ in Case 4 and Case 5, we add a nonlinear term besides two linear terms, which makes the linear FPCA method cannot fulfill the prediction with only two principal components. It can be observed that the prediction error of our nFunNN method is much lower than that of the FPCA method in Case 4 and Case 5. That indicates our method can achieve more effective dimension reduction results in nonlinear settings.

To sum up, the proposed nFunNN method gives great predicting results for all the five cases. Although the error can be a bit larger than that of the linear method in linear case, the predicting results of our nFunNN method is still reasonable. Furthermore, the nFunNN method shows obvious superiority for the nonlinear cases. Therefore, our nFunNN method can be a good choice, when we have no idea whether the data at hand has a linear structure or a nonlinear structure.

In this section, we discuss the effect of $L$, $J$, $K$, and $R$ on the performance of the proposed nFunNN method. Though the values of $J$ and $R$ can be different for our method, we set $J=R$ for simplicity here. We consider the settings of Case 1, Case 2 and Case 4 for explanation in this section. For the discussion of the effect of $L$, we fix $J=20$ and $K=2$, and the value of $L$ can be selected by $10$, $15$, and $20$. Moreover, we set $J=10,15,20$, and $L=10$, $K=2$ for the exploration of the influence of $J$. When discussing the effect of $K$, we set $K=2,3,4$, $L=10$ and $J=20$. Note that the number of parameters in the network is $J(KL+K+2L+2)$ when $J=R$. As the sample size should be larger than the number of parameters, we increase the sample size of training set to $2000$, and the size of test set is still $1000$ as in Section 3.1.

Tables 2–4 present the simulation results of our nFunNN method with the use of various $L$, $J$ and $K$. From Table 2, it can be observed that with the rise of $L$, both RMSE and RRMSE increase slightly. As shown in Table 3, there is no obvious difference in the prediction errors with the use of various $J$, which implies that the effect of $J$ is not significant. For results in Table 4, the prediction errors for Case 1 are similar when different values of $K$ are considered. However, for Case 2 and Case 4, which are both nonlinear cases, various values of $K$ lead to different performance of the network. In Case 2, the prediction error first decreases with the increase of $K$, and then shows a minor growth. Furthermore, the performance of the nFunNN method gets much better when larger $K$ is used in Case 4. We conjecture the reason is related to the complex setting of Case 4.

To summarize, according to the simulation results, the effects of $L$ and $J$ are not very obvious, while different values of $K$ can bring large changes for the prediction in some complex cases. Moreover, the selection of tuning parameters in the neural network can be completed through the validation set.