A Method of L1-Norm Principal Component Analysis for Functional Data

Abstract

Recently, with the popularization of intelligent terminals, research on intelligent big data has been paid more attention. Among these data, a kind of intelligent big data with functional characteristics, which is called functional data, has attracted attention. Functional data principal component analysis (FPCA), as an unsupervised machine learning method, plays a vital role in the analysis of functional data. FPCA is the primary step for functional data exploration, and the reliability of FPCA plays an important role in subsequent analysis. However, classical L2-norm functional data principal component analysis (L2-norm FPCA) is sensitive to outliers. Inspired by the multivariate data L1-norm principal component analysis methods, we propose an L1-norm functional data principal component analysis method (L1-norm FPCA). Because the proposed method utilizes L1-norm, the L1-norm FPCs are less sensitive to the outliers than L2-norm FPCs which are the characteristic functions of symmetric covariance operator. A corresponding algorithm for solving the L1-norm maximized optimization model is extended to functional data based on the idea of the multivariate data L1-norm principal component analysis method. Numerical experiments show that L1-norm FPCA proposed in this paper has a better robustness than L2-norm FPCA, and the reconstruction ability of the L1-norm principal component analysis to the original uncontaminated functional data is as good as that of the L2-norm principal component analysis.

1. Introduction

In recent years, with the rapid popularization of intelligent terminals and sensors, massive data have been rapidly accumulated, and the processing technology of intelligent big data has attracted more and more attention. Among these data, kinds of intelligent big data with function characteristics, such as physiological indicator data, growth curve data, air quality data, and temperature data, has also attracted people’s attention. In fact, these data are discrete samples of a continuous function, so such data are known in the literature as functional data [1,2,3,4,5,6,7,8,9,10]. The difference between functional data and traditional multivariate data is that the former regards the observed discrete data as a whole and as a realization of a random process. Therefore, the first step of statistical analysis is to fit the discrete data into smooth curves; this can solve the problems of missing data and inconsistent sampling intervals, which are difficult issues for multivariate data. Moreover, if the fitting curve is smooth enough, we can get more information from its derivatives, which is impossible for traditional multivariate data. As a nonparametric statistical method, functional data analysis is not limited by a model and its parameters, so it can better reflect real laws in nature. At present, statistical analysis methods for functional data have been widely used in the fields of biology, medicine, economics and meteorology [11,12,13,14,15,16,17,18].

Functional data principal component analysis (FPCA), as an unsupervised machine learning method, plays a vital role in the analysis of functional data. The central idea of FPCA is to use a few orthogonal dimensions to express most of the information of the original functional data. Through dimensionality reduction, the analysis of the original functional data can be transformed into the analysis of the characteristic functions of a few dimensions, thus greatly reducing the complexity of the functional data and allowing for the better interpretation of the function data. Since J.Q. Ramsay proposed the idea of functional principal component analysis in 1991 [19], various pieces of research on functional principal component analysis have emerged one after another. Classical functional principal components are the characteristic functions of the symmetric empirical covariance operator [20]. As early as 1982, Pousse and Romain studied the asymptotic properties of the characteristic functions of the empirical covariance operator: the empirical functional principal components [21]. In order to avoid the violent oscillation of the obtained principal component weight function, Rice and Silverman (1991) proposed a smooth functional principal component estimation method that smoothed the principal component weight function by adding penalties to the variance after projection [22]. The consistency of the estimate of the smooth functional principal component was then confirmed by Pezzulli and Silverman (1993) [23]. Silverman (1996) proposed another method of smooth functional principal components. Unlike the methods of Rice and Silverman (1991), the new method achieved the smoothness of the principal component function by penalizing the norm of the projected variance [24]. Gareth (2000) studied principal component analysis for sparse function data [25]. Boente (2000) studied the functional principal components-based kernel [26] Hall (2006) studied the properties of functional principal components [27]. Benko (2009) studied common functional principal components [28], and Hormann (2015) studied dynamic functional principal components [29].

Functional data principal component analysis (FPCA) is an important research subject of machine learning and artificial intelligence, and it is the primary step for functional data exploration. Therefore, the reliability of FPCA plays an important role in subsequent analysis. The aforementioned principal component methods for functional data were established in L2-norm framework. However, because the L2-norm enlarges the influence of outliers, the traditional functional principal components analysis method is sensitive to outliers. On the other hand, in regard to multivariate data, relevant research of principal component analysis methods [30,31,32,33,34,35,36,37] has shown that the principal component analysis method of L1-norm for multivariate data has a better robustness than that of the L2-norm. In [30], Kwak (2008) proposed an L1-PCA optimization model based on L1-norm maximization for multivariable data, i.e., $W_{L1}=\underset{W\in R^{D\times K},W^{T}W=I}{\arg \max }{\Vert W^{T}X\Vert }_1$. The algorithm in [30] gives an approximate solver for $W_{L1}=\underset{W\in R^{D\times K},W^{T}W=I}{\arg \max }{\Vert W^{T}X\Vert }_1$ through a sequence of deflating nullspace projections with cost $O(N^{2}DM)$, and it is robust to outliers and invariant to rotations. In [31], Nie et al. (2011) simultaneously approximated all M L1-PCs of X with complexity $O(N^{2}DM+N{M}^3)$; however, the principal components obtained by [31] were highly dependent on the the finding of the dimension M of a subspace. For example, the projection vector obtained when M = 1 may not be in a subspace obtained when M = 2. The optimal algorithm in [33] introduced a bit-flipping-based approximate solver for $W_{L1}=\underset{W\in R^{D\times K},W^{T}W=I}{\arg \max }{\Vert W^{T}X\Vert }_1$ with complexity $O(ND\min \{N,D\}+N^2(M^4+d{M}^2)+Nd{M}^3)$, where $d=rank(X)$; this solution has a low performance degradation, and is close to L2-PCA, but the cost is that it is not as robust as that in [30]. The work in [32] offered an algorithm for exact calculation $W_{L1}=\underset{W\in R^{D\times K},W^{T}W=I}{\arg \max }{\Vert W^{T}X\Vert }_1$ with complexity $O(2^{NM})$; however, when X is big data of large N and/or large dimension D, the cost is prohibitive. The authors of [34] studied the relationship of independent component analysis (ICA) and L1-PCA, and they proved that independent component analysis (ICA) can be performed by L1-norm PCA under the assumption of whitening. The authors of [36] computed L1-PCA by an incremental algorithm, in which only one measurement was processed at a time, and the changes in the nominal signal subspace could be tracked. Instead of maximizing the L1-norm deviation of the projected data, the authors of [35,37] focused on minimizing the L1-norm reconstruction error. However, in contrast to the conventional L2-PCA, the solutions of the minimization of the L1-norm reconstruction error might not be same as the solutions of the maximization of the L1-norm deviation of projected data.

Inspired by these pieces of research on L1-PCA for multivariable data, in this paper, we try to construct a robust L1-norm principal component analysis method for functional data (L1-norm FPCA). Firstly, we build a functional data L1-norm maximized principal component optimization model, and then a corresponding algorithm for solving the L1-norm maximized optimization model is extended to functional data based on the idea of a multivariate data L1-norm principal component analysis method [30]. Numerical experiments show that the L1-norm functional principal component analysis method provides a more robust estimation of principal components than the traditional L2-norm functional principal component analysis method (L2-norm FPCA). Finally, by comparing the reconstruction errors of the L1-norm FPCA and L2-norm FPCA, it is found that the reconstruction ability of the L1-norm principal components to the original uncontaminated functional data is as good as that of the L2-norm functional principal components.

2. Problem Description

2.1. L2-Norm Functional Principal Component Analysis (L2-Norm FPCA)

Suppose $x_1(t),x_2(t),\cdots ,x_n(t),t\in \tau \subset R$ are implementations of the square integrable random process $X(t)\in L_2(\tau )$. Without a loss of generality, we assume that $x_1(t),x_2(t),\cdots ,x_n(t),t\in \tau \subset R$ are centralized. The purpose of functional principal component analysis (FPCA) is to express as much information as possible of the original functional data with as few dimensions as possible. Firstly, the case of only one principal component is considered. At this point, the task of FPCA is to find a “projection direction” in infinite dimensional space so that the variance of projection of $x_1(t),x_2(t),\cdots ,x_n(t),t\in \tau$ to that direction is maximum. Assuming that the projection direction is ${\xi }_1(t)$, which is called the first functional principal component weight function of functional data $x_1(t),x_2(t),\cdots ,x_n(t),t\in \tau$, then ${\xi }_1(t)$ should be the solution of the following optimization problem:

$$\underset{}{\underset{{\xi }_1(t)}{\max }}\frac{1}{n}{\displaystyle \sum _{i=1}^n({\displaystyle \int {\xi }_1(t)x_i(t)dt{)}^2}}$$

(1)

$$s.t.\hspace{1em}{\displaystyle \int {\xi }_1{}^2(t)dt}=1$$

If the information that is expressed by one principal component is insufficient, a second projection direction ${\xi }_2(t)$, which is orthogonal to the first principal component direction ${\xi }_1(t)$ and maximizes the variance of the functional data $x_1(t),x_2(t),\cdots ,x_n(t), t\in \tau$ under the orthogonality condition, is necessary. This is the second functional principal component weight function. And so on, this process continues until the obtained principal components can express enough information. Therefore, the subsequent principal component weight functions need to satisfy the following optimization model:

$$\begin{array}{l}\underset{{\xi }_j(t)}{\max }\frac{1}{n}{\displaystyle \sum _{i=1}^n({\displaystyle \int {\xi }_j(t)x_i(t)dt{)}^2}}\\ \begin{array}{l}s.t.\hspace{1em}{\displaystyle \int {\xi }_j{}^2(t)dt}=1,j=2,3,\cdots m\\ {\displaystyle \hspace{1em}\int {\xi }_j(t){\xi }_k(t)}dt=0,k=1,2,\cdots ,j-1\end{array}\end{array}$$

(2)

J.Q. Ramsay proved that the principal component weight functions ${\xi }_1(t),{\xi }_2(t),\cdots ,{\xi }_m(t)$ of functional data $x_1(t),x_2(t),\cdots ,x_n(t),t\in \tau$ are the eigenfunctions that correspond to the first m largest eigenvalues of sample covariance function of functional data $x_1(t),x_2(t),\cdots ,x_n(t),t\in \tau$, i.e., ${\displaystyle \int \stackrel{\wedge }{C}(t,s){\xi }_i(t)dt={\rho }_i{\xi }_i(s)},i=1,2,\cdots ,m$, where ${\rho }_1\ge {\rho }_2\ge \cdots \ge {\rho }_m$ are the eigenvalues of the covariance function $\widehat{C}(t,s)$.

From the optimization Formulas (1) and (2), it is easy to find that the above L2-norm functional principal components enlarge the influence of outliers and are sensitive to outliers. Therefore, L1-norm functional principal components are constructed in this paper. Compared with the traditional L2-norm, the L1-norm weakens the influence of outliers. It can be expected that the L1-norm functional principal components have a good anti-noise ability.

2.2. L1-Norm Functional Principal Component Analysis (L1-Norm FPCA)

Suppose $x_1(t),x_2(t),\cdots ,x_n(t)$ are the implementations of square integrable stochastic process $x(\cdot )$. Without a loss of generality, suppose that $x_1(t),x_2(t),\cdots ,x_n(t)$ have been centralized. Now we want to find an m-dimensional linear subspace so that the L1-norm dispersion of the projection of $x_1(t),x_2(t),\cdots ,x_n(t)$ in this subspace is the largest. Assume that the subspace is spanned by ${\beta }_1(t),{\beta }_2(t),\cdots ,{\beta }_m(t)$, and the optimization problem corresponding to Formulas (1) and (2) can be obtained:

$$\underset{{\beta }_1,{\beta }_2,\cdots {\beta }_m}{\max }{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m\left|{\displaystyle \int {\beta }_j(t)x_i(t)dt}\right|}}$$

(3)

$$\begin{array}{l}s.t.\hspace{1em}\hspace{1em}{\displaystyle \int {\beta }_j{}^2(t)}dt=1,j=1,2,\cdots ,m\\ \hspace{1em}\hspace{1em}\hspace{1em}{\displaystyle \int {\beta }_j(t){\beta }_k(t)dt=0}j,k=1,2,\cdots ,m;j\ne k\end{array}$$

where ${\beta }_1(t),{\beta }_2(t),\cdots ,{\beta }_m(t)$ are called L1-norm principal component weight functions for $x_1(t),x_2(t),\cdots ,x_n(t)$.

It is not easy to solve the Optimization Problem (3) because the objective function is non-differentiable, non-convex, and contains an absolute value operation. Next, we try to find the solution of Optimization Problem (3) from the perspective of orthogonal basis expansion.

Assuming that $x_1(t),x_2(t),\cdots ,x_n(t);{\beta }_1(t),{\beta }_2(t),\cdots ,{\beta }_m(t),t\in \tau$ can be linearly represented by the same standard orthogonal basis functions ${\varphi }_1(t),{\varphi }_2(t),\cdots$ with the same number of basis functions $K$, i.e.,

$$x_i(t)={\displaystyle \sum _{v=1}^K{c}_{iv}{\varphi }_v(t)=c_i^T\varphi (t)}i=1,2,\cdots ,n$$

$${\beta }_j(t)={\displaystyle \sum _{v=1}^K{b}_{jv}{\varphi }_v(t)=b_j^T\varphi (t)}j=1,2,\cdots ,m$$

where $K$ is a positive integer. $c_i={(c_{i1},c_{i2},\cdots ,c_{iK})}^T;b_i={(b_{j1},b_{j2},\cdots ,b_{jK})}^T;$ $\varphi (t)={({\varphi }_1(t),{\varphi }_2(t),\cdots ,{\varphi }_K(t))}^T;$

Under the above assumptions, we get:

$$\underset{{\beta }_1,{\beta }_2,\cdots {\beta }_m}{\max }{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m\left|{\displaystyle \int {\beta }_j(t)x_i(t)dt}\right|}}=\underset{b_1,b_2,\cdots ,b_m}{\max }{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m\left|{\displaystyle \int b_j^T\varphi (t){\varphi }^T(t)c_{i}dt}\right|}}=\underset{b_1,b_2,\cdots ,b_m}{\max }{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m\left|b_j^T{c}_i\right|}}$$

Since ${\displaystyle \int {\beta }_j{}^2(t)}dt={\displaystyle \int b_j^T\varphi (t){\varphi }^T(t)b_j}dt=b_j^T{b}_j$, the constraints ${\displaystyle \int {\beta }_j{}^2(t)}dt=1,j=1,2,\cdots ,m$ can be expressed as $b_j^T{b}_j=1,j=1,2,\cdots ,m$, and the constrains ${\displaystyle \int {\beta }_j(t){\beta }_k(t)dt=0}$ $j,k=1,2,\cdots ,m;$ $j\ne k$ can be expressed as $b_j^T{b}_k=0j,k=1,2,\cdots ,m;j\ne k$. Therefore, the Optimization Problem (3) can be transformed into the following Optimization Problem (4):

$$\underset{b_1,b_2,\cdots ,b_m}{\max }{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m\left|b_j^T{c}_i\right|}}=\underset{b_1,b_2,\cdots ,b_m}{\max }{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m\left|{\displaystyle \sum _{k=1}^K{b}_{jk}{c}_{ik}}\right|}}$$

(4)

$$\begin{array}{l}s.t.\hspace{1em}\hspace{1em}{b}_j^T{b}_j=1j=1,2,\cdots ,m\\ \hspace{1em}\hspace{1em}\hspace{1em} b_j^T{b}_k=0j,k=1,2,\cdots ,m;j\ne k\end{array}$$

If we can get the solution of Optimization Problem (4), according to ${\beta }_j(t)=b_j^T\varphi (t),$ $j=1,2,\cdots ,m$, we can get the solution of Optimization Problem (3). There are several algorithms to solve Optimization Problem (4), such as those in [30,31,33], each of which has its own advantages. According to the goal of building robust principal components for functional data, we finally choose the algorithm in [30], because the principal components calculated in [30] are more robust to outliers, and this algorithm is relatively low-complexity when the data number, data dimension, and the principal components number are large.

Next, based on the orthogonal basis expansion of functional data, we employ the L1-norm PCA algorithm of multivariate data [30] to get the solving algorithm of the L1-norm functional principal component weight functions (Abbreviation: L1-FPCA algorithm). The algorithm is rewritten in the next section.

3. The Solving Algorithm of L1-Norm Functional Principal Component Weight Functions (L1-FPCA Algorithm)

3.1. Only One Principal Component

First, we discuss the case where there is only one principal component, namely m = 1. In this case, the Optimization Problems (3) and (4) are, respectively, simplified as follows:

$$\underset{\beta (t)}{\max }{\displaystyle \sum _{i=1}^n\left|{\displaystyle \int \beta (t)x_i(t)dt}\right|}$$

(5)

$$s.t.\hspace{1em}{\displaystyle \int {\beta }^2(t)}dt=1$$

and

$$\underset{b}{\max }{\displaystyle \sum _{i=1}^n\left|b_{}^T{c}_i\right|}$$

(6)

$$s.t.\hspace{1em}{b}_{}^{T}b=1$$

Next, we construct L1-FPCA algorithm to solve the Optimization Problems (5) and (6).

L1-FPCA Algorithm:

Step 1: Arbitrarily choose the initial projection direction ${\beta }^0(t)$, get $b^0$ by ${\beta }^0(t)={(b^0)}^T\varphi (t)$, normalize $b^0$:$b^0=\frac{b^0}{{\Vert b^0\Vert }_2}$, and set the iteration number $k$ to be 0.

Step 2: For all $i\in (1,2,\cdots ,n)$, if ${\displaystyle \int {\beta }^k(t)x_i(t)dt<0}$, i.e., ${(b^k)}^T{c}_i<0$, let $p_i^k=-1$; otherwise $p_i^k=1$.

Step 3: Let $b^k={\displaystyle \sum _{i=1}^n{p}_i^{k-1}{c}_i}$, normalize $b^k$: $b^k=\frac{b^k}{{\Vert b^k\Vert }_2}$, and get the corresponding ${\beta }^k(t)$ by ${\beta }^k(t)={(b^k)}^T\varphi (t)$.

Step 4: If ${\beta }^k(t)\ne {\beta }^{k-1}(t)$, return to step 2. If there is $i$ such that ${\displaystyle \int {\beta }^k(t)x_i(t)dt=0}$, i.e., ${(b^k)}^T{c}_i=0$, then let $b^k=\frac{(b^K+\Delta b)}{{\Vert b^K+\Delta b\Vert }_2}$ and get the corresponding ${\beta }^k(t)$, then return to step 2, where $\Delta b$ is a small non-zero vector. Otherwise, let ${\beta }^{\ast }(t)={\beta }^k(t)$, $b^{*}=b^k$ and ${\beta }^{\ast }(t)={(b^{\ast })}^T\varphi (t)$, stop.

Theorem 1.

The L1-FPCA algorithm is convergent, and its convergence point$b^{*}$is the local maximum point of the Optimization Problem (6) and${\beta }^{\ast }(t)$is the local maximum point of Optimization Problem (5).

Proof. 

First, we prove that the objective function ${\displaystyle \sum _{i=1}^n\left|{\displaystyle \int {\beta }^k(t)x_i(t)dt}\right|}$ and ${\displaystyle \sum _{i=1}^n\left|b_{}^T{c}_i\right|}$ are nondecreasing in the iteration process of L1-FPCA algorithm, i.e.,

$$\begin{array}{l}{\displaystyle \sum _{i=1}^n\left|{\displaystyle \int {\beta }^k(t)x_i(t)dt}\right|}={\displaystyle \sum _{i=1}^n\left|{(b^k)}^T{c}_i\right|}={\displaystyle \sum _{i=1}^n{(b^k)}^T{p}_i^k{c}_i}={(b^k)}^T{\displaystyle \sum _{i=1}^n{p}_i^k{c}_i}\ge {(b^k)}^T{\displaystyle \sum _{i=1}^n{p}_i^{k-1}{c}_i}\\ \ge {(b^{k-1})}^T{\displaystyle \sum _{i=1}^n{p}_i^{k-1}{c}_i}={\displaystyle \sum _{i=1}^n{(b^{k-1})}^T{p}_i^{k-1}{c}_i}={\displaystyle \sum _{i=1}^n\left|{(b^{k-1})}^T{c}_i\right|}\\ ={\displaystyle \sum _{i=1}^n\left|{\displaystyle \int {\beta }^{k-1}(t)x_i(t)dt}\right|}\end{array}$$

Therefore, the objective function ${\displaystyle \sum _{i=1}^n\left|{\displaystyle \int {\beta }^k(t)x_i(t)dt}\right|}$ and ${\displaystyle \sum _{i=1}^n\left|b_{}^T{c}_i\right|}$ are nondecreasing. Additionally, because there are only finite number of data points, the convergence points ${\beta }^{\ast }(t)$ and $b^{*}$ of the L1-FPCA algorithm exist.

Next, we prove that $b^{*}$ and ${\beta }^{\ast }(t)$ are the local maxima of the corresponding optimization problem.

Suppose that $b^{*}=b^k$, that is the convergence point $b^{*}$ is found after $k$ iterations. Because for any $i\in (1,2,\cdots ,n)$, $b^T{p}_i^{k}c\to {(b^{*})}^T{p}_i^k{c}_i$, there is a neighborhood $N(b^{*})$ of $b^{*}$ so that for $b\in N(b^{*})$, $b^T{p}_i^{k}c\ge 0$ and

$${\displaystyle \sum _{i=1}^n\left|{\displaystyle \int {\beta }^{\ast }(t)x_i(t)}\right|}={\displaystyle \sum _{i=1}^n\left|{(b^{\ast })}^T{c}_i\right|}={\displaystyle \sum _{i=1}^n{(b^{\ast })}^T{p}_i^k{c}_i}={(b^{\ast })}^T{\displaystyle \sum _{i=1}^n{p}_i^k{c}_i}$$

Because $b^{*}$ is the convergence point, $b^{*}$ is parallel to ${\displaystyle \sum _{i=1}^n{p}_i^k{c}_i}$; therefore, ${(b^{\ast })}^T{\displaystyle \sum _{i=1}^n{p}_i^k{c}_i}\ge b^T{\displaystyle \sum _{i=1}^n{p}_i^k{c}_i}={\displaystyle \sum _{i=1}^n\left|b^T{c}_i\right|}$, so for $b\in N(b^{*})$, ${\displaystyle \sum _{i=1}^n\left|{\displaystyle \int {\beta }^{\ast }(t)x_i(t)}\right|}={\displaystyle \sum _{i=1}^n\left|{(b^{\ast })}^T{c}_i\right|}\ge {\displaystyle \sum _{i=1}^n\left|b^T{c}_i\right|}={\displaystyle \sum _{i=1}^n\left|\beta (t)x_i(t)\right|}$; that is, $b^{*}$ is the local maximum of ${\displaystyle \sum _{i=1}^n\left|b_{}^T{c}_i\right|}$ and ${\beta }^{\ast }(t)$ is the local maximum of ${\displaystyle \sum _{i=1}^n\left|{\displaystyle \int \beta (t)x_i(t)dt}\right|}$.

Therefore, the L1-FPCA procedure finds a local maximum point $b^{*}$ of ${\displaystyle \sum _{i=1}^n\left|b_{}^T{c}_i\right|}$ and ${\beta }^{\ast }(t)$ of ${\displaystyle \sum _{i=1}^n\left|{\displaystyle \int \beta (t)x_i(t)dt}\right|}$. □

Since the L1-FPCA algorithm obtains a local optimal solution, we expect to find the global optimal solution with great probability by appropriately setting the initial projection direction ${\beta }^0(t)$, e.g., by setting ${\beta }^0(t)=\arg \underset{x_i(t)}{\max }{\displaystyle \int x_i^2(t)dt}$ or by setting it to be the solution of L2-FPCA. In practice, we usually select several different initial projection directions ${\beta }^0(t)$ and calculate the respective local optimal solutions, and the solution with maximized the objective function ${\displaystyle \sum _{i=1}^n\left|{\displaystyle \int \beta (t)x_i(t)dt}\right|}$ is selected as the optimal solution.

3.2. Multiple Principal Components

Suppose that m principal components (m > 1) are needed, and the L1-FPCA algorithm needs to sequentially find m principal component projection directions $b_1,b_2,\cdots ,b_m$ and corresponding ${\beta }_1(t),{\beta }_2(t),\cdots ,{\beta }_m(t)$. The specific algorithm is as follows:

Step 1: Let ${\beta }_0(t)=0$, i.e.,$b_0=0$,${\{c_i^0=c_i\}}_{i=1}^n$.

Step 2: For all $i\in (1,2,\cdots ,n)$, let $c_i^1=c_i^0-b_0(b_0^T{c}_i^0)$ and apply the L1-FPCA algorithm to $c^1=(c_1^1,c_2^1,\cdots ,c_n^1)$ to obtain the projection vector $b_1$ and the corresponding ${\beta }_1(t)$.

Step 3: For all $i\in (1,2,\cdots ,n)$, let $c_i^j=c_i^{j-1}-b_{j-1}(b_{j-1}^T{c}_i^{j-1})$ and apply the L1-FPCA algorithm to $c^j=(c_1^j,c_2^j,\cdots ,c_n^j)$ to obtain the projection vector $b_j$ and the corresponding ${\beta }_j(t)$.

Step 4: Repeat Step 3 until m projection vectors $b_1,b_2,\cdots ,b_m$ and corresponding ${\beta }_1(t),{\beta }_2(t),\cdots ,{\beta }_m(t)$ are obtained.

Since $b_1,b_2,\cdots ,b_m$ are standard orthogonal dimensions in $R^K$ space [38], the principal component weight functions ${\beta }_1(t),{\beta }_2(t),\cdots ,{\beta }_m(t)t\in \tau$ are also standard orthogonal dimensions because:

$$<{\beta }_j(t),{\beta }_k(t)>={\displaystyle \int {\beta }_j(t){\beta }_k(t)}dt={\displaystyle \int b_j^T\varphi (t)}{\varphi }^T(t)b_{k}dt=b_j^T{b}_k=0j,k=1,2,\cdots ,m;j\ne k$$

$$<{\beta }_j(t),{\beta }_j(t)>={\displaystyle \int {\beta }_j(t){\beta }_j(t)}dt={\displaystyle \int b_j^T\varphi (t)}{\varphi }^T(t)b_{j}dt=b_j^T{b}_j=1j=1,2,\cdots ,m$$

As with the L2-norm functional principal component analysis, it is necessary to consider how many principal components are appropriate. This problem needs to be determined by the cumulative variance contribution rate. That is, according to the variance of the j projection direction, $v_j=\frac{1}{n}{\displaystyle \sum _{i=1}^n({\displaystyle \int {\beta }_j(t)x_i(t)dt{)}^2}},j=1,2,\cdots ,m$, the total variance of the first k projection directions can be calculated as $S(k)={\displaystyle \sum _{j=1}^k{v}_j}$, and the total variance of the original functional data is $S=\frac{1}{n}{\displaystyle \sum _{i=1}^n({\displaystyle \int x_i^2(t)dt{)}^2}}$. Thus, for the actual problems, the number of final principal component weight functions can be determined when $\frac{S(k)}{S}$ is more than 80% or 85%.

4. Numerical Examples

4.1. Simulation

In order to compare the robustness to outliers of L1-norm functional principal components (L1-FPCs) that are proposed in this paper and the classical L2-norm functional principal components (L2-FPCs), we performed this simulation. We referred to the simulation setting given by Fraiman and Muniz (2001) [38]. Here, we considered that functional data $x_1(t),x_2(t),\cdots ,x_n(t)$ are the implementations of squared integrable stochastic process $X(\cdot )$, and the function curves were generated from different model. There was no contamination in Model 1, and several other models suffered from different types contamination based on Model 1.

Model 1 (no contamination): $x_i(t)=m(t)+{\epsilon }_i(t),i=1,2,\cdots ,n$, where error term ${\epsilon }_i(t)$ is a stochastic Gaussian process with zero mean and covariance function $\mathrm{cov}(s,t)=(1/2){(1/2)}^{0.9\left|t-s\right|}$ and $m(t)=4t$,$t\in [0,1]$.

Model 2 (asymmetric contamination): $y_i(t)=x_i(t)+c_{i}M,i=1,2,\cdots ,n$, where $c_i$ is the sample of the 0–1 distribution with the parameter $q$, and $M$ is the contamination constant.

Model 3 (symmetric contamination): $y_i(t)=x_i(t)+c_i{\sigma }_{i}M,i=1,2,\cdots ,n$, where $c_i$ and $M$ are defined as in Model 2 and ${\sigma }_i$ is a sequence of random variables with values of 1 and −1 with a probability of 1/2 that is independent of $c_i$.

Model 4 (partially contaminated): $y_i(t)=\{\begin{array}{l}{x}_i(t)+c_i{\sigma }_{i}M,t\ge T_i\\ x_i(t),\hspace{1em}\hspace{1em}\hspace{1em} t<T_i\end{array},i=1,2,\cdots ,n$, where $T_i$ is a random number generated from a uniform distribution on [0,1].

Model 5 (peak contamination): $y_i(t)=\{\begin{array}{l}{x}_i(t)+c_i{\sigma }_{i}M,T\le t\le T_i+l\\ x_i(t),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}t\notin [T_i,T_i+l]\end{array},i=1,2,\cdots ,n$, where $l=1/15$ and $T_i$ is a random number generated from a uniform distribution on $[0,1-l]$.

Figure 1 shows the simulated curves of these five models. For each model, we set 100 equal-interval sampling points in [0,1] and generated 200 replications. For Model 1, the parameter $q$ was 0 and the contamination constant $M$ was 0. For several other contaminated models, we considered several levels of contamination, with q = 5% and 10% and contamination constants M = 5 and 10. When fitting function curves, we use generalized cross validation (GCV) to obtain the number of bases. The results showed that the number of bases of Model 1–3 were the same, while those of Models 4 and 5 were different. However, due to the need of calculating the change of principal component coefficient, we had to calculate it on the same basis. Therefore, for comparison purposes, in Models 4 and 5, we selected the same number of bases as that of Model 1.

Classical L2-norm FPCA and L1-norm FPCA were used for the simulated functional data corresponding to these five models. We focused on their robust to various abnormal disturbances. When implementing L1-norm FPCA on Model 1, by comparing the value of objective function, the initial value was chosen as the first L2-norm functional principal component weight function, i.e., ${\beta }^0(t)=\xi (t)$, where $\xi (t)$ is the eigenfunction corresponding to the largest eigenvalue of the sample covariance function of the functional data in Model 1. Because the L1-norm FPCA of the following several disturbance models should be compared with Model 1, in order to ensure the consistency of conditions when calculating the L1-norm FPCA of the following several disturbance models, the initialization values also adopted the eigenfunction corresponding to the largest eigenvalue of the sample covariance function of the corresponding functional data.

The sums of absolute values of the coefficient differences of several principal components under non-contamination and contamination were compared and analyzed. The sum of the absolute values of the corresponding coefficient changes are given in

Table 1. The sum of the absolute values of the first principal component weight function coefficient changes for no contamination and asymmetric contamination (5% and 10%).

	M = 10		M = 5
q	L1-norm FPC	L2-norm FPC	L1-norm FPC	L2-norm FPC
	1st FPC	1st FPC	1st FPC	1st FPC
5%	0.17	1.13	0.13	0.91
10%	0.22	1.24	0.18	1.16

,

Table 2. The sum of the absolute values of the first principal component weight function coefficient changes for no contamination and symmetric contamination (5% and 10%).

	M = 10		M = 5
q	L1-norm FPC	L2-norm FPC	L1-norm FPC	L2-norm FPC
	1st FPC	1st FPC	1st FPC	1st FPC
5%	0.2	1.13	0.1	0.96
10%	0.25	1.17	0.23	1.04

,

Table 3. The sum of the absolute values of the first principal component weight function coefficient changes for no contamination and partial contamination (5% and 10%).

	M = 10		M = 5
q	L1-norm FPC	L2-norm FPC	L1-norm FPC	L2-norm FPC
	1st FPC	1st FPC	1st FPC	1st FPC
5%	1.17	13.47	0.81	10.67
10%	1.70	14.77	1.24	12.29

and

Table 4. The sum of the absolute values of the first four principal component weight functions. Coefficient changes for no contamination and peak contamination (5% and 10%).

	M = 10								M = 5
q	L1-norm FPC				L2-norm FPC				L1-norm FPC				L2-norm FPC
	1st	2nd	3rd	4th	1st	2nd	3rd	4th	1st	2nd	3rd	4th	1st	2nd	3rd	4th
5%	0.5	1.5	3.9	9.8	3.9	66.3	39.2	45.2	0.2	0.7	1.3	10.1	3.9	66.8	39.9	45.7
10%	0.8	1.8	4.1	8.4	5.4	24.2	48.4	54.3	0.4	1.2	2.2	3.4	2.3	10.9	22.3	42.2

. Since the variance contribution rate of the first principal component reached 80%, only the first principal component function was taken in Models 1, 2, 3 and 4. However, in order to achieve a similar variance contribution rate, at least four principal components were needed in Model 5. Thus, for Models 1, 2, 3 and 4, we only show the changes of the first principal component function. For Model 5, we show the changes of the first four principal component functions.

It can be seen from Table 1, Table 2, Table 3 and Table 4 that under the same contamination ratio and contamination size, the coefficient changes of the principal component weight functions of the L1-norm were significantly smaller than those of the L2-norm, which shows that the functional principal components of the L1-norm were more stable than those of the L2-norm, no matter which form of contamination was received. This conclusion can also be confirmed from the boxplots of the coefficient changes of the principal component weight functions.

As can be seen from Figure 2, Figure 3, Figure 4 and Figure 5, in the same contamination ratio and size, the changes of L1-norm principal component coefficient are more concentrated near zero compared with the changes of the L2-norm principal component coefficient, which shows that under the same contamination mode, L1-norm functional principal components were more robust to outliers and more reliable.

From the above research, we found that the L1-norm functional principal components were more robust than L2-norm functional principal components. Thus, how can one reconstruct the original functional data with these two types of principal components? In order to study this problem, we reconstructed the original uncontaminated functional data with the same number of functional principal components of L1-norm and L2-norm under each model. The scatter plots of the coefficients of the two types of reconstructed error curves are shown in Figure 6, Figure 7, Figure 8 and Figure 9.

In Figure 6, Figure 7, Figure 8 and Figure 9, we can see that the scatter plots of the reconstruction error curve coefficients of L1-norm and L2-norm were always near the line y = x under the first three pollution models, and under peak pollution, the reconstruction error of the L1-norm was smaller than that of the L2-norm. When using the paired one-sided T-test, the p-values were found to all be close to 1, indicating that the reconstruction error curve coefficients of the L1-norm were not greater than those of the L2-norm. Thus, the reconstruction ability of the L1-norm principal components to the original uncontaminated functional data was not worse than that of the L2-norm principal components. The results of the paired one-sided T-test are shown in

Table 5. The table of the one-sided paired T-test of the coefficients of the reconstruction error curves of L1-norm and L2-norm under asymmetric contamination (Alternative hypothesis: The true difference of the reconstruction error curve coefficients of L1-norm and L2-norm in means was greater than 0.).

q	M = 5			M = 10
	t	df	p-Value	t	df	p-Value
5%	−2.8447	199	0.9975	−2.1651	199	0.9842
10%	−2.2484	199	0.9872	−2.5843	199	0.9948

,

Table 6. The table of the one-sided paired T-test of the coefficients of the reconstruction error curves of L1-norm and L2-norm under symmetric contamination (Alternative hypothesis: The true difference of the reconstruction error curve coefficients of L1-norm and L2-norm in means was greater than 0.).

q	M = 5			M = 10
	t	df	p-Value	t	df	p-Value
5%	−3.8761	199	0.9999	−3.34	199	0.9995
10%	−4.7628	199	1	−3.5293	199	0.9997

,

Table 7. The table of the one-sided paired T-test of the coefficients of the reconstruction error curve of L1-norm and L2-norm under partial contamination (Alternative hypothesis: The true difference of the reconstruction error curve coefficients of L1-norm and L2-norm in means was greater than 0.).

q	M = 5			M = 10
	t	df	p-Value	t	df	p-Value
5%	−5.2373	199	1	−4.9371	199	1
10%	−7.7896	199	1	−5.033	199	1

and

Table 8. The table of the one-sided paired T-test of the coefficients of the reconstruction error curves of L1-norm and L2-norm under peak contamination (Alternative hypothesis: The true difference of the reconstruction error curve coefficients of L1-norm and L2-norm in means was greater than 0.).

q	M = 5			M = 10
	t	df	p-Value	t	df	p-Value
5%	−6.6502	199	1	−6.6212	199	1
10%	−7.6313	199	1	−6.8564	199	1

.

The above experiments showed that the functional principal component of the L1-norm was not just stable and reliable, it also had the same reconstruction ability as the L2-norm.

4.2. Canadian Weather Data

We used Canadian weather data, which provide daily temperatures at 35 different locations in Canada averaged over 1960–1994, in order to compare the robust to outliers of the L1-norm functional principal components and L2-norm functional principal components when the functional data were contaminated by abnormal data. Firstly, by considering the periodic characteristics of the data, the discrete temperature observation data were fitted into 35 functional curves by a Fourier basis function, and the number of the basis functions was 65. The fitting curves are shown in Figure 10a. As can be seen when using the function data outlier detection method [39], the temperature modes of the four stations of Vancouver, Victoria, Pr. Rupert and Resolute were different from those of the other stations. Figure 10b shows this function after removing the data from these four observatories. The functional data of the 35 observatories were called the whole data, and the functional data after removing Vancouver, Victoria, Pr. Rupert and Resolute were normal data, so the whole data can be understood as the addition of abnormalities to the normal data.

In order to compare the robustness between the L2-norm functional principal component weighting functions and the L1-norm functional principal component weighting functions to outliers, the L2-norm functional principal components and L1-norm functional principal components were, respectively, used for normal data and data added with outliers. For each method, the results of the two cases were compared, because the variance contribution rate of the first two principal components reached 90%, though the latter analysis only focused on the first two functional principal components.

Figure 11 shows the change of the first principal component weight function before and after adding outliers by using two functional principal component analysis methods. Figure 11a is a graph of the first principal component weight function that was obtained by using the L2-norm functional principal component method. The solid line is the result of normal data, and the dashed line is the result of adding four abnormal curves. Figure 11b is a graph of the first principal component weight function that was obtained by using the proposed L1-norm functional principal component method. After comparing the objective function, the initial value was chosen as the first L2-norm functional principal component weight function, i.e., ${\beta }^0(t)=\xi (t)$, where $\xi (t)$ is the eigenfunction corresponding to the largest eigenvalue of the sample covariance function of normal functional data and the same method for whole functional data. The solid line is the result of normal data, and the dashed line is the result of adding four abnormal curves. By comparing the coefficients of the two first functional principal component weighting functions, it was found that the sum of the absolute change of the coefficients of the first principal component weighting functions that were obtained by the L1-norm method before and after adding abnormal values was 0.16, which was less than the 0.18 corresponding to the L2-norm. Next, the performance of the second principal component weight function is discussed.

Figure 12 shows the change of the second principal component weight function before and after the addition of outliers by using two functional principal component analysis methods. Figure 12a is a graph of the second principal component weight function that was obtained by using the L2-norm function principal component method. The solid line is the result of normal data, and the dashed line is the result of adding four abnormal curves. Figure 12b is a graph of the second principal component weight function that was obtained by using the proposed L1-norm function principal component method. The solid line is the result of normal data, and the dashed line is the result of adding four abnormal curves. By comparing the coefficients of the two second function principal component weighting functions, it was found that the sum of absolute change of the coefficients of the second principal component weighting functions that were obtained by the L1-norm method before and after adding abnormal values was 0.33, which was less than the 0.76 corresponding to the L2-norm. the sums of the absolute values of the coefficient change of the principal component weight functions under the two methods are shown in Table 5.

Table 9. The sum of absolute change of the coefficients of the first two principal component weighting functions.

The Sum of Absolute Change of the Coefficients	The 1st Function Principal Component Weighting Function	The 2nd Function Principal Component Weighting Function
L2-norm	0.18	0.76
L1-norm	0.16	0.33

shows that the classical L2-norm principal components weight functions greatly changed before and after removing outliers, reflecting its sensitivity to outliers. However, the L1-norm functional principal components weight functions presented in this paper had little change before and after adding abnormal values. Therefore, that the L1-norm principal component weight function proposed in this paper has a strong anti-noise ability and a good stability.

We also compared the reconstruction ability of two types of principal components to normal data. The scatter plots of the coefficients of the two types of reconstructed error curves are shown in Figure 13.

From Figure 13, it can be seen the scatter plots of the reconstruction error curve coefficients of L1-norm and L2-norm were always near the line y = x; When we performed a paired one-sided T-test for the two groups of reconstruction error curve coefficients, the t value was found to 1.0323, the degree of freedom for the t-statistic was 33, and the p-value was 0.1547, which indicates that the reconstruction error curve coefficients of L1-norm were not greater than those of L2-norm. Thus, the reconstruction ability of the L1-norm principal components to the original uncontaminated functional data was not worse than the L2-norm principal components.

5. Concluding Remarks

FPCA is a primary step for functional data exploration, and the reliability of FPCA plays an important role in subsequent analysis. The existing principal component methods of functional data were established in an L2-norm framework. However, because the L2-norm enlarges the influence of outliers, the traditional functional principal components analysis method is sensitive to outliers. On the other hand, in regard to multivariate data, the relevant research on the principal component analysis method [30,31,32,33,34,35,36,37] have shown that the principal component analysis method of L1-norm for multivariate data has a better robustness than that of the L2-norm. Motivated by this research, in this paper, we tried to construct an L1-norm principal component analysis method for functional data. Firstly, we built a functional data L1-norm maximized principal component optimization model. Then, a corresponding algorithm for solving the L1-norm maximized optimization model was constructed based on the idea of multivariate data L1-norm principal component analysis method [30]. An extensive simulation study was conducted, and a real dataset of Canadian weather was employed to assess the robustness of the L1-norm functional principal component analysis. From the simulation study that considered different contamination configurations (symmetric, asymmetric, partial and peak), we found that the L1-norm functional principal component analysis method provides a more robust estimation of principal components than the traditional L2-norm principal component analysis method. Finally, by comparing the reconstruction errors of the L1-norm FPCA and L2-norm FPCA, it was found that the reconstruction ability of the L1-norm principal components to the original uncontaminated functional data is as good as that of the L2-norm principal components. Therefore, when functional data contain outliers, the estimation given by the L1-norm functional principal component analysis method is more reliable. The proposed L1-norm FPCA may prove to be an useful addition to functional data analysis.