# Sparse Principal Component Analysis via Fractional Function Regularity

Authors: Xuanli Han; Jigen Peng; Angang Cui; Fujun Zhao
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, we describe a novel approach to sparse principal component analysis (SPCA) via a nonconvex sparsity-inducing fraction penalty function SPCA (FP-SPCA). Firstly, SPCA is reformulated as a fraction penalty regression problem model. Secondly, an algorithm corresponding to the model is proposed and the convergence of the algorithm is guaranteed. Finally, numerical experiments were carried out on a synthetic data set, and the experimental results show that the FP-SPCA method is more adaptable and has a better performance in the tradeoff between sparsity and explainable variance than SPCA.


---

## Body

## 1. Introduction

Principal component analysis (PCA) [1] has become more popular in data analysis, dimension reduction, image processing, and feature extraction [2]. It seeks the linear combinations of the original variables such that the derived variables capture maximal variance of the total original variables. We can obtain PCA(s) of an original variable by the singular value decomposition (SVD) of the data matrix which is composed of the observed variables. Let X be an n×p matrix, where n and p are the number of observations and the number of variables, respectively. Suppose the average value of columns of X are all zero and the singular value decomposition of X is(1)X=UDV⊤.Then, the columns of matrixUD are the principal components of X, and the columns of V are the corresponding loadings of the principal components. The PCA method is popular due to the following two properties: first, principal components capture the maximum variability from the columns of matrix X, which guarantees minimal information loss; second, principal components are uncorrelated, so we can use one of them without considering the others.At the same time, however, PCA has its own drawback, i.e., each principal component of matrixX is a linear combination of all p variables and the elements of the loading vectors (the columns of matrix V) are usually nonzero, which make them difficult to interpret the gained PCs.Based on the abovementioned priority properties and drawback, some scholars consider it will be a wise option to keep the dimension reduction property at the same time to reduce the number of explainable variables. For this purpose, one method is to set the loading whose absolute value is smaller than a threshold to zero.The same problem also arises in multiple linear regression, where the response variable is explained by a linear combination of the explainable variables. But, which are the important explainable variables that account for response variable information most? To answer this question, researchers have proposed several approaches, and the first type method is the lasso-based approaches [3–6]. Others, for example, Efron et al. [7] proposed the LARS which is sensitive to noise of samples in 2004, and Friedman et al. [8] put forward a pathwise coordinate optimization approach in 2007. The second type is the lasso-bootstrap methods [9, 10]. Other methods including overview type articles [11–13] and Radchenko et al. [14] proposed a variable inclusion and shrinkage method, Fan and Li [15] proposed a nonconcave penalized likelihood method, and Candès and Tao [16] proposed a dantzig selector method. Among them, the lasso method [5] is a well-known variable selection technique which produces an accurate and sparse model. As time went on, Zou and Hastie [4] proposed the elastic net in 2005, an adaptive method, which has the advantages of both lasso and ridge regression. The elastic net approach can be considered as a generalization of lasso and ridge regression. In 2006, Zou et al. [3] proposed sparse principal analysis (SPCA) using the lasso and ridge regression to produce a modified principal component with sparse loadings. They show that PCA can be formulated as a regression-type optimization problem and sparse loadings are obtained by imposing l1-norm and l2-norm regularization on the regression coefficients of the model. In 2009, Leng and Wang [17] proposed a method of simple adaptive sparse principal component analysis, which uses the adaptive lasso penalty instead of the lasso penalty in SPCA. Moreover, they replace the least-squares objective function in SPCA by a general least-squares objective function, which leads to study many related SPCA estimators under a unified theoretical framework and get the method of general adaptive sparse principal component analysis. SPCA can produce a sparse model that has only fewer active coefficients, and the other coefficients are all set to zero. The proposed sparse model SPCA is more interpretable than PCA. In a word, SPCA can specifically identify structure information of the data matrix.The idea of making the model’s coefficients sparse is not a new one. Jollife et al. [18] first proposed a SPCA method using l1-norm regularization which leads to a variety of sparsity-inducing methods. The success of SPCA in gaining more interpretable model motivates us to propose the method in this paper, where we proposed a fractional function penalty method with respect to sparse PCA. The method is to exploit the fractional function [19] to penalize the model coefficients in order to obtain sparse coefficients, which makes the PCA to have more interpretable ability.The rest of the paper is organized as follows. In Section2, definition and derivation of principal components are reviewed. The sparse principal component is presented in Section 3. In Section 4, we propose the fractional functional penalized principal components. Numerical experiments and conclusions are formulated in Sections 5 and 6, respectively.In this paper, without specification, scalar is denoted as lower case letterx, and vector is denoted as bold lower case letter x=x1,x2,…,xn⊤. The matrix is denoted by bold capital letter M, and I denotes the identity matrix. The transpose of a real matrix M is denoted by M⊤. The Frobenius norm of · is denoted by ⋅, where · denotes column vector or matrix, and the spectral norm of a matrix M is denoted by M2.


## 2. Principal Components and Sparse Principal Components

### 2.1. Definition of Derivation of Principal Components

Letx=x1,x2,…,xp be a vector of p random variables, and the variances of the p random variables and the covariance between the p variables will be of importance. If p is large, an alternative key approach is to seek for fewer variables that preserve well most of the information provided by variances and covariances of these p variables. One of the alternative techniques is principle component analysis (PCA) which concerns more variance of variables than covariance. PCA first looks for a linear function xα1 owning to maximum variance of x, where α1=α11,α12,…,α1p⊤, that is,(2)xα1=α11x1+α12x2+⋯+α1pxp.Then, in the same way, we can find a linear functionxα2, in which xα2 and xα1 are linear independence, so that at the k-th step, a linear function xαk can be found which has maximum variance and is linear independence with the former k−1 linear functions. The xαk is the k-th principal component (PC) [1]. Many of the estimation problems in array signal processing can also start with the same issue in (2); however, there may be exist some perturbations [20–22]. Having known what are the PCs, the remaining question is how to find them. If the vector of random variables x has a known covariance matrix Σ whose i,jth entry is the known covariance of random variables xi and xj when i≠j, and variance of the element xi of x=x1,x2,…,xp when i=j. When Σ is unknown, the more realistic method is to substitute Σ for sample covariance matrix S.In fact, becauseΣ is often unknown, we consider the case S in the following only.To obtain the PCs of the sampling matrixXn×p, where n and p are the number of observations and variables of random vector x, respectively. Consider the coefficient vector α1∈ℝp, which is forced to maximize the sample variance of Xα1, where α1 is normalized. Let xi denote the i-th column of sampling matrix X, and xi is zero-centered, i=1,2,…,p. The problem can be reformulated as the follows:(3)maxα1α1⊤X⊤Xα1subject toα1⊤α1=1,and solving (3), we can obtain the relationship X⊤Xα1=λ1α1; furthermore, we have(4)α1⊤X⊤Xα1=α1⊤λ1α1=λ1α1⊤α1=λ1,where λ1 is the largest eigenvalue of matrix X⊤X and α1 is the eigenvector of X⊤X corresponding to λ1. Xα1 is called the 1-th PC of sampling matrix X.In the same way, thek-th PC of X is αk and the sample variance of Xαk is λk, where λk is the k-th largest eigenvalue of X⊤X and αk is the corresponding eigenvector.In general, ifn×p matrix X is known, xi, the i-th column of X, is zero-centered. The singular value decomposition of X is UDV⊤. Then, the principal components of matrix X, which are denoted by Y, are UDr, where UDr is the matrix composed of the first r columns of the matrix UD and r is the rank of X. So, we can get yi=Xvi [1, 3], where yi and vi are ith columns of matrix X and V, respectively.Although principal component analysis is widely used in various technical fields, it has an obvious drawback, that is, each principal component is a linear combination of all original variables, but it is often difficult to interpret the result. To solve this problem, Zou et al. [3] proposed sparse principal component analysis. They first consider PCA as a linear regression question.


### 2.2. Preparation of Sparse Principal Components

#### 2.2.1. Linear Regression and LASSO

Consider matrixXn×p with n observations of the p variables, and let y=y1,y2,…,yn⊤ be the response vector and X=x1,x2,…,xp the predictor variables, where xj=x1j,x2j,…,xnj⊤ is the j-th column of matrix Xj=1,2,…,p. Assume y and xj have zero means. The question of evaluating of PCA of matrix X can be considered as a linear regression problem:(5)minβ=β1,β2,…,βp⊤y−∑j=1pβjxj2=minβ=β1,β2,…,βp⊤y−Xβ2.This model is simple, and the result is easy to get in some cases; however, the result is difficult to interpret [3]. Based on model (5), Tibshirani [5] proposed the lasso method:(6)β^lasso=argminβy−∑j=1pβjxj2+λ∑j=1pβj,where β=β1,β2,…,βp⊤ and λ≥0.Model (6) can produce accurate and sparse result at the same time, which strengths the result’s interpretable ability. However, the lasso method has some limitations, and the most relevant one is that the number of variables produced by the lasso is limited by the number of observations. Given Xn×p, when the number of rows n and number of columns p of X satisfy p>n, the lasso can select n explainable variables at most.To overcome the abovementioned drawback, Zou and Hastie [4] proposed the elastic net model:(7)β^en=1+λ2argminβy−∑j=1pβjxj2+λ2∑j=1pβj2+λ1∑j=1pβj,where β=β1,β2,…,βp⊤, λ1≥0, and λ2≥0.The elastic net model shares the properties of lasso and ridge regression at the same time because whenλ2=0 the elastic is the lasso, and when λ1=0, it becomes the ridge regression. When p>n, by choosing proper λ2, the elastic net model can potentially include all variables.


#### 2.2.2. Ridge Regression and the Naive Elastic Net

In model (6) and (7), the sparse coefficients are all l1-norm-induced and the sparsity is none of l2-norm’s business. Jolliffe et al. [18] proposed the SCoTLASS method:(8)αi⊤X⊤Xαis.t.αi⊤αi=1,αl⊤αi=0i≥2,l<i,∑j=1pαij≤t,where αi=αi1,αi2,…,αip⊤ and t is a parameter. This method has the ability to obtain sparse loadings by imposing l1-norm penalty on principle component analysis (PCA).The SCoTLASS method has two drawbacks. One is difficulty of choosing parametert and another problem is the high computational cost. So, Zou et al. [3] consider a modified method. They first show that principle component analysis (PCA) can be expressed as a ridge regression problem.Suppose the singular value decomposition ofXn×p is UDV⊤. Let Ui and Dii be the i-th column of matrix U and the i-th element of diagonal matrix D, respectively, and then, DiiUi is the i-th principal component of matrix Xn×p,i=1,2,…,r, r=rankXn×p. The ridge estimation of β is given by the following expression:(9)β^ridge=argminβDiiUi−Xn×pβ2+λβ2.Letv^=β^ridge/β^ridge, and then, v^=vi, where λ>0 and vi is the i-th column of V.This conclusion shows the relationship between principal component analysis and the regression method. The termλβ2 in (9) is to guarantee the unique of β^ridge when the inverse of X⊤X does not exist. Furthermore, Zou and Hastie [4] proposed a method to add l1-norm penalty to model (9) and got the following optimization problem:(10)β^nen=argminβDiiUi−Xβ2+λβ2+λ1β1,They call thev^i=β^nen/β^nen is an approximation of vi, and Xv^i is the ith approximated principal component. To distinguish this model from model (7), model (10) is called the naive elastic net.


### 2.3. Sparse Principle Component Based on PCA

Based on the general conclusion in the end of Section2.1, Zou et al. [3] proposed a two-stage analysis: perform PCs first, and then, find suitable spare approximation of PCs.Consider the firstr principal components of Xn×p, where r=rankXn×p. Let Ap×r=α1,α2,…,αr, Bp×r=β1,β2,…,βr, and(11)A^,B^=argminA,BX⊤−AB⊤X⊤2+λB2s.t.A⊤A=Ir,for any λ>0. Then, βi^∝vi,i=1,2,…,r, where βi^ is the estimation of βi.Then, we can perform SPCA, which is based on the PCA and regression method, by addingl1-norm penalty to the regression coefficients to induce sparse loadings. Thus, the following result is obtained:(12)A^,B^=argminA,BX−XBA⊤2+λB2+∑j=1rλ1,jβj1,=argminA,BX⊤−AB⊤X⊤2+λB2+∑j=1rλ1,jβj1,=argminA,B∑i=1nxi−AB⊤xi2+λ∑j=1rβj2+∑j=1rλ1,jβj1,subject to A⊤A=Ir. Where xi is the ith column of matrix X⊤ and βj is the jth column of matrix B, r=rankX.


## 2.1. Definition of Derivation of Principal Components

Letx=x1,x2,…,xp be a vector of p random variables, and the variances of the p random variables and the covariance between the p variables will be of importance. If p is large, an alternative key approach is to seek for fewer variables that preserve well most of the information provided by variances and covariances of these p variables. One of the alternative techniques is principle component analysis (PCA) which concerns more variance of variables than covariance. PCA first looks for a linear function xα1 owning to maximum variance of x, where α1=α11,α12,…,α1p⊤, that is,(2)xα1=α11x1+α12x2+⋯+α1pxp.Then, in the same way, we can find a linear functionxα2, in which xα2 and xα1 are linear independence, so that at the k-th step, a linear function xαk can be found which has maximum variance and is linear independence with the former k−1 linear functions. The xαk is the k-th principal component (PC) [1]. Many of the estimation problems in array signal processing can also start with the same issue in (2); however, there may be exist some perturbations [20–22]. Having known what are the PCs, the remaining question is how to find them. If the vector of random variables x has a known covariance matrix Σ whose i,jth entry is the known covariance of random variables xi and xj when i≠j, and variance of the element xi of x=x1,x2,…,xp when i=j. When Σ is unknown, the more realistic method is to substitute Σ for sample covariance matrix S.In fact, becauseΣ is often unknown, we consider the case S in the following only.To obtain the PCs of the sampling matrixXn×p, where n and p are the number of observations and variables of random vector x, respectively. Consider the coefficient vector α1∈ℝp, which is forced to maximize the sample variance of Xα1, where α1 is normalized. Let xi denote the i-th column of sampling matrix X, and xi is zero-centered, i=1,2,…,p. The problem can be reformulated as the follows:(3)maxα1α1⊤X⊤Xα1subject toα1⊤α1=1,and solving (3), we can obtain the relationship X⊤Xα1=λ1α1; furthermore, we have(4)α1⊤X⊤Xα1=α1⊤λ1α1=λ1α1⊤α1=λ1,where λ1 is the largest eigenvalue of matrix X⊤X and α1 is the eigenvector of X⊤X corresponding to λ1. Xα1 is called the 1-th PC of sampling matrix X.In the same way, thek-th PC of X is αk and the sample variance of Xαk is λk, where λk is the k-th largest eigenvalue of X⊤X and αk is the corresponding eigenvector.In general, ifn×p matrix X is known, xi, the i-th column of X, is zero-centered. The singular value decomposition of X is UDV⊤. Then, the principal components of matrix X, which are denoted by Y, are UDr, where UDr is the matrix composed of the first r columns of the matrix UD and r is the rank of X. So, we can get yi=Xvi [1, 3], where yi and vi are ith columns of matrix X and V, respectively.Although principal component analysis is widely used in various technical fields, it has an obvious drawback, that is, each principal component is a linear combination of all original variables, but it is often difficult to interpret the result. To solve this problem, Zou et al. [3] proposed sparse principal component analysis. They first consider PCA as a linear regression question.


## 2.2. Preparation of Sparse Principal Components

### 2.2.1. Linear Regression and LASSO

Consider matrixXn×p with n observations of the p variables, and let y=y1,y2,…,yn⊤ be the response vector and X=x1,x2,…,xp the predictor variables, where xj=x1j,x2j,…,xnj⊤ is the j-th column of matrix Xj=1,2,…,p. Assume y and xj have zero means. The question of evaluating of PCA of matrix X can be considered as a linear regression problem:(5)minβ=β1,β2,…,βp⊤y−∑j=1pβjxj2=minβ=β1,β2,…,βp⊤y−Xβ2.This model is simple, and the result is easy to get in some cases; however, the result is difficult to interpret [3]. Based on model (5), Tibshirani [5] proposed the lasso method:(6)β^lasso=argminβy−∑j=1pβjxj2+λ∑j=1pβj,where β=β1,β2,…,βp⊤ and λ≥0.Model (6) can produce accurate and sparse result at the same time, which strengths the result’s interpretable ability. However, the lasso method has some limitations, and the most relevant one is that the number of variables produced by the lasso is limited by the number of observations. Given Xn×p, when the number of rows n and number of columns p of X satisfy p>n, the lasso can select n explainable variables at most.To overcome the abovementioned drawback, Zou and Hastie [4] proposed the elastic net model:(7)β^en=1+λ2argminβy−∑j=1pβjxj2+λ2∑j=1pβj2+λ1∑j=1pβj,where β=β1,β2,…,βp⊤, λ1≥0, and λ2≥0.The elastic net model shares the properties of lasso and ridge regression at the same time because whenλ2=0 the elastic is the lasso, and when λ1=0, it becomes the ridge regression. When p>n, by choosing proper λ2, the elastic net model can potentially include all variables.


### 2.2.2. Ridge Regression and the Naive Elastic Net

In model (6) and (7), the sparse coefficients are all l1-norm-induced and the sparsity is none of l2-norm’s business. Jolliffe et al. [18] proposed the SCoTLASS method:(8)αi⊤X⊤Xαis.t.αi⊤αi=1,αl⊤αi=0i≥2,l<i,∑j=1pαij≤t,where αi=αi1,αi2,…,αip⊤ and t is a parameter. This method has the ability to obtain sparse loadings by imposing l1-norm penalty on principle component analysis (PCA).The SCoTLASS method has two drawbacks. One is difficulty of choosing parametert and another problem is the high computational cost. So, Zou et al. [3] consider a modified method. They first show that principle component analysis (PCA) can be expressed as a ridge regression problem.Suppose the singular value decomposition ofXn×p is UDV⊤. Let Ui and Dii be the i-th column of matrix U and the i-th element of diagonal matrix D, respectively, and then, DiiUi is the i-th principal component of matrix Xn×p,i=1,2,…,r, r=rankXn×p. The ridge estimation of β is given by the following expression:(9)β^ridge=argminβDiiUi−Xn×pβ2+λβ2.Letv^=β^ridge/β^ridge, and then, v^=vi, where λ>0 and vi is the i-th column of V.This conclusion shows the relationship between principal component analysis and the regression method. The termλβ2 in (9) is to guarantee the unique of β^ridge when the inverse of X⊤X does not exist. Furthermore, Zou and Hastie [4] proposed a method to add l1-norm penalty to model (9) and got the following optimization problem:(10)β^nen=argminβDiiUi−Xβ2+λβ2+λ1β1,They call thev^i=β^nen/β^nen is an approximation of vi, and Xv^i is the ith approximated principal component. To distinguish this model from model (7), model (10) is called the naive elastic net.


## 2.2.1. Linear Regression and LASSO

Consider matrixXn×p with n observations of the p variables, and let y=y1,y2,…,yn⊤ be the response vector and X=x1,x2,…,xp the predictor variables, where xj=x1j,x2j,…,xnj⊤ is the j-th column of matrix Xj=1,2,…,p. Assume y and xj have zero means. The question of evaluating of PCA of matrix X can be considered as a linear regression problem:(5)minβ=β1,β2,…,βp⊤y−∑j=1pβjxj2=minβ=β1,β2,…,βp⊤y−Xβ2.This model is simple, and the result is easy to get in some cases; however, the result is difficult to interpret [3]. Based on model (5), Tibshirani [5] proposed the lasso method:(6)β^lasso=argminβy−∑j=1pβjxj2+λ∑j=1pβj,where β=β1,β2,…,βp⊤ and λ≥0.Model (6) can produce accurate and sparse result at the same time, which strengths the result’s interpretable ability. However, the lasso method has some limitations, and the most relevant one is that the number of variables produced by the lasso is limited by the number of observations. Given Xn×p, when the number of rows n and number of columns p of X satisfy p>n, the lasso can select n explainable variables at most.To overcome the abovementioned drawback, Zou and Hastie [4] proposed the elastic net model:(7)β^en=1+λ2argminβy−∑j=1pβjxj2+λ2∑j=1pβj2+λ1∑j=1pβj,where β=β1,β2,…,βp⊤, λ1≥0, and λ2≥0.The elastic net model shares the properties of lasso and ridge regression at the same time because whenλ2=0 the elastic is the lasso, and when λ1=0, it becomes the ridge regression. When p>n, by choosing proper λ2, the elastic net model can potentially include all variables.


## 2.2.2. Ridge Regression and the Naive Elastic Net

In model (6) and (7), the sparse coefficients are all l1-norm-induced and the sparsity is none of l2-norm’s business. Jolliffe et al. [18] proposed the SCoTLASS method:(8)αi⊤X⊤Xαis.t.αi⊤αi=1,αl⊤αi=0i≥2,l<i,∑j=1pαij≤t,where αi=αi1,αi2,…,αip⊤ and t is a parameter. This method has the ability to obtain sparse loadings by imposing l1-norm penalty on principle component analysis (PCA).The SCoTLASS method has two drawbacks. One is difficulty of choosing parametert and another problem is the high computational cost. So, Zou et al. [3] consider a modified method. They first show that principle component analysis (PCA) can be expressed as a ridge regression problem.Suppose the singular value decomposition ofXn×p is UDV⊤. Let Ui and Dii be the i-th column of matrix U and the i-th element of diagonal matrix D, respectively, and then, DiiUi is the i-th principal component of matrix Xn×p,i=1,2,…,r, r=rankXn×p. The ridge estimation of β is given by the following expression:(9)β^ridge=argminβDiiUi−Xn×pβ2+λβ2.Letv^=β^ridge/β^ridge, and then, v^=vi, where λ>0 and vi is the i-th column of V.This conclusion shows the relationship between principal component analysis and the regression method. The termλβ2 in (9) is to guarantee the unique of β^ridge when the inverse of X⊤X does not exist. Furthermore, Zou and Hastie [4] proposed a method to add l1-norm penalty to model (9) and got the following optimization problem:(10)β^nen=argminβDiiUi−Xβ2+λβ2+λ1β1,They call thev^i=β^nen/β^nen is an approximation of vi, and Xv^i is the ith approximated principal component. To distinguish this model from model (7), model (10) is called the naive elastic net.


## 2.3. Sparse Principle Component Based on PCA

Based on the general conclusion in the end of Section2.1, Zou et al. [3] proposed a two-stage analysis: perform PCs first, and then, find suitable spare approximation of PCs.Consider the firstr principal components of Xn×p, where r=rankXn×p. Let Ap×r=α1,α2,…,αr, Bp×r=β1,β2,…,βr, and(11)A^,B^=argminA,BX⊤−AB⊤X⊤2+λB2s.t.A⊤A=Ir,for any λ>0. Then, βi^∝vi,i=1,2,…,r, where βi^ is the estimation of βi.Then, we can perform SPCA, which is based on the PCA and regression method, by addingl1-norm penalty to the regression coefficients to induce sparse loadings. Thus, the following result is obtained:(12)A^,B^=argminA,BX−XBA⊤2+λB2+∑j=1rλ1,jβj1,=argminA,BX⊤−AB⊤X⊤2+λB2+∑j=1rλ1,jβj1,=argminA,B∑i=1nxi−AB⊤xi2+λ∑j=1rβj2+∑j=1rλ1,jβj1,subject to A⊤A=Ir. Where xi is the ith column of matrix X⊤ and βj is the jth column of matrix B, r=rankX.


## 3. Fractional Function-Penalized Sparse Principal Components

### 3.1. Fraction Function-Penalized Method

Recently, Li et al. [19] proposed a fraction-penalized function. They study the problem of a nonconvex sparsity promoting penalty function,(13)Pax=∑i=1npaxi=∑i=1naxi1+axi,to relax the l0-norm in compressed sensing [15, 16, 23, 24] to get the sparsity signal, where a>0 is a parameter and x=x1,x2,…,xn⊤. Pax is a fraction function which is nonconvex and sparse promoting. They studied the properties of this fractional function and derived a closed form expression of the optimal solution to the fractional function penalty problem and proposed an iterative thresholding algorithm to solve the fraction function penalty problem. A series of experiments have been conducted and the results show, compared with the soft thresholding algorithm and the half thresholding algorithms, that the proposed algorithm has the property of better sparse signal recovery with and without measurement noise.Based on the good performance of the fraction function in sparse promotion in compressed sensing [23, 25], in this paper, we shall adopt this method to the sparse principal component analysis (SPCA), i.e., replace the l1-norm penalty in SPCA with the fraction function penalty.Combiningfd12(12) and (13)fd13, we can get(14)A^,B^=argminA,B∑i=1nxi−AB⊤xi2+λ∑j=1rβj2+∑j=1rλ1,jPaβj,subject toA⊤A=Ir.We now introduce the following fact for the purpose of solving the problem (14).(15)∑i=1nxi−AB⊤xi2=X⊤−AB⊤X⊤2=X−XBA⊤2.Owing toA⊤A=Ir, let A⊥ be the matrix satisfying AA⊥ a p×p orthonormal matrix, which AA⊥ is a block matrix. Let Q:=AA⊥. Then,(16)X−XBA⊤Q2=trX−XBA⊤Q⊤X−XBA⊤Q,=trQ⊤X−XBA⊤⊤X−XBA⊤Q,=trX−XBA⊤⊤X−XBA⊤,=X−XBA⊤2,while(17)X−XBA⊤Q2=X−XBA⊤AA⊥2,=X−XBA⊤AX−XBA⊤A⊥2,=XA−XBXA⊥2,=XA−XB2+XA⊥2,so(18)X−XBA⊤2=XA−XB2+XA⊥2,=∑j=1rXαj−Xβj2+XA⊥2.So, problem (14) can be transformed into(19)A^,B^=argminA,B∑j=1rXαj−Xβj2+XA⊥2+λ∑j=1rβj2+∑j=1rλ1,jPaβj,subject toA⊤A=Ir.Problem (19) can be solved by an alternating algorithm. If A is given, B^=β^1,β^2,…,β^r can be estimated by(20)β^j=argminβjXαj−Xβj2+λβj2+λ1,jPaβj,and if B is given, problem (14) is to solve(21)A^=argminA∑i=1nxi−AB⊤xi2,=argminAX−XBA⊤2subject toA⊤A=Ir,where xi is the i-th column of matrix X⊤. According to the procrustes rotation [26–28], the estimation of A can be obtained; suppose the SVD of matrix X⊤XB is U1D1V1⊤, by U1V1⊤ [3], that is, A^=U1V1⊤.


### 3.2. Transformation of Problem (16)

In (20), let A=α1,α2,…,αr be given by αj=vj for j=1,2,…,r, where vi is the i-th column of V-the matrix in the UDV⊤ which is the singular decomposition of the matrix X. Let yj=Xvj, and (22) can be reformulated as(22)β^j=argminβjyj−Xβj2+λβj2+λ1,jPaβj.Let(23)X∗=11+λXλIn+p×p,yj∗=yj0n+p×1,γj=λ1,j/1+λ and βj∗=1+λβj. Then, (22) can be rewritten as(24)β^j∗=argminβj∗yj∗−X∗βj∗2+γjPaβj∗.The relationship ofβj∗^ and βj^ can be explained as the following. Let(25)Fβj=yj−Xβj2+λβj2+λ1,jPaβj,Hβj=Gβj∗=yj∗−X∗βj∗2+γjPaβj∗.According to the expressions ofX∗,yj∗,γj,βj∗ and the definition of Paβj, we can get(26)Fβj=yj−Xβj2+λβj2+λ1,j∑i=1paβji1+aβji,Hβj=yj−Xβj2+λβj2+λ1,j1+λ∑i=1pa1+λβji1+a1+λβji,=yj−Xβj2+λβj2+λ1,j∑i=1paβji1+a1+λβji.Owing toλ>0, it is easy to obtain that Hβj<Fβj, so that(27)minβjHβj<minβjFβj,that is,(28)minβj∗Gβj∗<minβjFβj0.BothGβj∗ and Hβj are continuous functions with respect to βj∗ and βj, respectively, so as λ⟶0, we have βj∗⟶βj, and further, we get Gβj∗⟶Gβj⟶Fβj. That is, Gβj∗⟶Fβj as λ⟶0.So, we conclude thatβ^j can be approximated well by βj∗^ when positive λ is small enough.Remark 1.
The reason why we substituteFβj with Gβj∗ is model (24) is easier to solve than model (22).
Now, we will show that the optimal solution to the minimization problem,(29)minβj∗yj∗−X∗βj∗2+γjPaβj∗,can be expressed as a thresholding operation.


### 3.3. The Algorithm of the Problem (29)

For anyγj,μ∈0,+∞, and z∈ℝp, let(30)Cγjβj∗=yj∗−X∗βj∗2+γjPaβj∗,Cμβj∗,z=μCγjβj∗−X∗βj∗−X∗z2+βj∗−z2,Bμβj∗=βj∗+μX∗⊤yj∗−X∗βj∗.We have the following results.Lemma 1.
For any fixed parameterμ,a,γj∈0,+∞ and z∈ℝp, ifβ¯j∗=β¯j1∗,β¯j2∗,…,β¯jp∗⊤ is a local minimizer to Cμβj∗,z, then(31)β¯ji∗=0⟺Bμzi≤t,β¯ji∗=gγjμBμzi⟺Bμzi>t,where i=1,2,…,p,(32)t=γj2a,ifγj≤1a2,γj−12a,ifγj>1a2,and according to [19],(33)ϕs=arccos27γjμa241+as3−1,gγjμs=sgns1+as/31+2cosϕs/3−π/3−1a.Proof.
We can see thatCμβj∗,z can be expressed as(34)Cμβj∗,z=βj∗−z−μX∗⊤X∗z+μX∗⊤yj∗2+γjμPaβj∗+μyj∗2+z2−μX∗z2−z−μX∗⊤X∗z+μX∗⊤yj∗2,=βj∗−Bμz2+γjμPaβj∗+μyj∗2+z2−μX∗z2−Bμz2,that is to say, minimizing Cμβj∗,z for z, given μ, γj, and a, is equivalent to(35)minβj∗∈Rp∑i=1pβji∗−Bμzi2+γjμ∑i=1ppaβji∗.
So,β¯j∗=β¯j1∗,β¯j2∗,…,β¯jp∗⊤ is a local minimizer of Cμβj∗,z if and only if β¯ji∗ solves the following problem:(36)minβji∗∈Rβji∗−Bμzi2+γjμpaβji∗.
According to Lemma 10 in [19], we complete the proof.
Furthermore, we can get Theorem1.Theorem 1.
Ifβ˜j∗=β˜j1∗,β˜j2∗,…,β˜jp∗⊤ is an optimal solution to (29), and γj is positive value and μ satisfies 0<μ≤X∗2−2, and then, the optimal solution β˜ji∗ can be given by(37)β˜ji∗=gγjμBμβ˜j∗i,ifBμβ˜j∗i>t,0,ifBμβ˜j∗i≤t,where i=1,2,…,p,(38)t=γjμa2,ifγj≤1a2μ,γjμ−12a,ifγj>1a2μ.Proof.
According to0<μ≤X∗2−2, we have(39)Cμβj∗,β˜j∗=μyj∗−X∗βj∗2+γjPaβj∗−X∗βj∗−X∗β˜j∗2+βj∗−β˜j∗2,=μyj∗−X∗βj∗2+γjPaβj∗−μX∗βj∗−β˜j∗2+βj∗−β˜j∗2≥μyj∗−X∗βj∗2+γjPaβj∗≥μyj∗−X∗β˜j∗2+γjPaβ˜j∗,=Cμβ˜j∗,β˜j∗.
That is to say, forβj∗∈Rp, β˜j∗ is a local minimizer of Cμβj∗,β˜j∗ as long as β˜j∗ is a solution to (29). According to Lemma 10 in [19] and Lemma 1, we complete the proof.
With the expression ofβ˜j∗ in Theorem 1, we can get an alternating thresholding algorithm for problem (29), that is,(40)βj∗n+1i=gγjμBμβj∗ni,for i=1,2,…,p, where Bμβj∗=βj∗+μX∗⊤yj∗−X∗βj∗ and gγjμ· is the thresholding operator in Lemma 1. We call this algorithm the alternating iterative FP-SPCA thresholding algorithm.Remark 2.
The dot· in gγjμ· may be a real number as in Lemma 1 or a vector whose components are all real numbers. If the dot · represents a vector, gγjμ· means the result vector which components are the result that gγjμ· acts on every component of the vector · sequentially.
It is well known that the solution to problem (29) depends on the choice of the regularization parameter γj. Before we give out the proper choice for the regularization parameter γj, we need the following definition.Definition 1.
The nonincreasing rearrangement of the vectorβ¯j∈ℝp is the vector β¯j∈ℝp for which(41)β¯j1≥β¯j2≥⋯≥β¯jp,and there is a permutation π:1,2,…,p⟶1,2,…,p with β¯ji=β¯jπi for all i∈1,2,…,p, and β¯jπi is the πi-th component of the vector β¯j.
Similar to the choice of the regularization parameter in FP algorithm (see Scheme 2 in [19]), here, we can approximate β˜j∗ by using the current step value βj∗n and the γj can be chosen as(42)γj=2aμBμβj∗nk+1,if2aμBμβj∗nk+1≤1a2μ,1−ε¯4a2μ2aBμβj∗nk+12,if2aμBμβj∗nk+1>1a2μ,where ε¯ is a small positive number, such as 0.1, 0.01, or 0.001, and k is sparsity degree of β˜j∗. By this way, the iterative algorithm is adaptive and free from the regular parameter choice. Noticing that 0<μ<X∗2−2 in the (2) of Theorem 2, we take μ=1−ε/X∗22, where ε∈0,1.
With abovementioned parameter-choosing method, according to the iterative method of (40), we have the following alternating iterative FP-SPCA thresholding algorithm. If A is given, for each j=1,2,…,r, yj∗ is the same as in formula (24), by using the iterative formula (40), βj∗n+1 is obtained. If B is fixed, according to (21), it is only try to minimize ∑i=1nxi−AB⊤xi2 = X−XBA⊤2 with respect to A, where A⊤A=Ir. According to the reduced rank form of the Procrustes rotation in [26], if the singular value decomposition of matrix X⊤XB is U1D1V1⊤, the solution is A^=U1V1⊤ [3]. The alternating iterative FP-SPCA thresholding algorithm is described in Algorithm 1.
At last, we discuss the convergence of the iterative FP-SPCA thresholding algorithm which converges to a stationary point of the iteration (40).Algorithm 1: Iterative FP-SPCA thresholding algorithm.
Initialize: Choose βj∗0, k, μ=1−ε/X∗22, ϵ¯, W, a, ϵ, and λ.while not converged doforj=1:rBμβj∗n=βj∗n+μX∗⊤yj∗−X∗βj∗nγj1=2/aμBμβj∗k+1; γj2n=1−ε¯/4a2μ2aBμβj∗nk+12ifγj1≤1/a2μ thenγj=γj1; t=γjμa/2elseγj=γj2; t=γjμ−1/2aendλ1,j=γj1+λ;fori=1:pifBμβj∗ni>tβj∗n+1i=gγjμBμβj∗nielseβj∗n+1i=0endendW:,j=βj∗n+1end;C:=X∗⊤X∗WU1,S1,V1=svdCV=U1V1⊤n⟶n+1end whileform=1:rβ^m=W:,m/normW:,mendoutput:β^m,m=1,2,…,r.Theorem 2.
Letβj∗n be the sequence generated by the iterative FP-SPCA thresholding algorithm (40) with 0<μ<X∗2−2. Then,(i)
The sequenceCγjβj∗n=yj∗−X∗βj∗n2+γjPaβj∗n is monotonically decreasing(ii)
βj∗n is asymptotically regular, i.e., βj∗n+1−βj∗n⟶0asn⟶∞(iii)
βj∗n converges to a stationary point of the iteration (40)Proof.
The proof is similar to the proof of Theorem 4.1 in [29] and Lemma 2 in [30], so it is omitted.


## 3.1. Fraction Function-Penalized Method

Recently, Li et al. [19] proposed a fraction-penalized function. They study the problem of a nonconvex sparsity promoting penalty function,(13)Pax=∑i=1npaxi=∑i=1naxi1+axi,to relax the l0-norm in compressed sensing [15, 16, 23, 24] to get the sparsity signal, where a>0 is a parameter and x=x1,x2,…,xn⊤. Pax is a fraction function which is nonconvex and sparse promoting. They studied the properties of this fractional function and derived a closed form expression of the optimal solution to the fractional function penalty problem and proposed an iterative thresholding algorithm to solve the fraction function penalty problem. A series of experiments have been conducted and the results show, compared with the soft thresholding algorithm and the half thresholding algorithms, that the proposed algorithm has the property of better sparse signal recovery with and without measurement noise.Based on the good performance of the fraction function in sparse promotion in compressed sensing [23, 25], in this paper, we shall adopt this method to the sparse principal component analysis (SPCA), i.e., replace the l1-norm penalty in SPCA with the fraction function penalty.Combiningfd12(12) and (13)fd13, we can get(14)A^,B^=argminA,B∑i=1nxi−AB⊤xi2+λ∑j=1rβj2+∑j=1rλ1,jPaβj,subject toA⊤A=Ir.We now introduce the following fact for the purpose of solving the problem (14).(15)∑i=1nxi−AB⊤xi2=X⊤−AB⊤X⊤2=X−XBA⊤2.Owing toA⊤A=Ir, let A⊥ be the matrix satisfying AA⊥ a p×p orthonormal matrix, which AA⊥ is a block matrix. Let Q:=AA⊥. Then,(16)X−XBA⊤Q2=trX−XBA⊤Q⊤X−XBA⊤Q,=trQ⊤X−XBA⊤⊤X−XBA⊤Q,=trX−XBA⊤⊤X−XBA⊤,=X−XBA⊤2,while(17)X−XBA⊤Q2=X−XBA⊤AA⊥2,=X−XBA⊤AX−XBA⊤A⊥2,=XA−XBXA⊥2,=XA−XB2+XA⊥2,so(18)X−XBA⊤2=XA−XB2+XA⊥2,=∑j=1rXαj−Xβj2+XA⊥2.So, problem (14) can be transformed into(19)A^,B^=argminA,B∑j=1rXαj−Xβj2+XA⊥2+λ∑j=1rβj2+∑j=1rλ1,jPaβj,subject toA⊤A=Ir.Problem (19) can be solved by an alternating algorithm. If A is given, B^=β^1,β^2,…,β^r can be estimated by(20)β^j=argminβjXαj−Xβj2+λβj2+λ1,jPaβj,and if B is given, problem (14) is to solve(21)A^=argminA∑i=1nxi−AB⊤xi2,=argminAX−XBA⊤2subject toA⊤A=Ir,where xi is the i-th column of matrix X⊤. According to the procrustes rotation [26–28], the estimation of A can be obtained; suppose the SVD of matrix X⊤XB is U1D1V1⊤, by U1V1⊤ [3], that is, A^=U1V1⊤.


## 3.2. Transformation of Problem (16)

In (20), let A=α1,α2,…,αr be given by αj=vj for j=1,2,…,r, where vi is the i-th column of V-the matrix in the UDV⊤ which is the singular decomposition of the matrix X. Let yj=Xvj, and (22) can be reformulated as(22)β^j=argminβjyj−Xβj2+λβj2+λ1,jPaβj.Let(23)X∗=11+λXλIn+p×p,yj∗=yj0n+p×1,γj=λ1,j/1+λ and βj∗=1+λβj. Then, (22) can be rewritten as(24)β^j∗=argminβj∗yj∗−X∗βj∗2+γjPaβj∗.The relationship ofβj∗^ and βj^ can be explained as the following. Let(25)Fβj=yj−Xβj2+λβj2+λ1,jPaβj,Hβj=Gβj∗=yj∗−X∗βj∗2+γjPaβj∗.According to the expressions ofX∗,yj∗,γj,βj∗ and the definition of Paβj, we can get(26)Fβj=yj−Xβj2+λβj2+λ1,j∑i=1paβji1+aβji,Hβj=yj−Xβj2+λβj2+λ1,j1+λ∑i=1pa1+λβji1+a1+λβji,=yj−Xβj2+λβj2+λ1,j∑i=1paβji1+a1+λβji.Owing toλ>0, it is easy to obtain that Hβj<Fβj, so that(27)minβjHβj<minβjFβj,that is,(28)minβj∗Gβj∗<minβjFβj0.BothGβj∗ and Hβj are continuous functions with respect to βj∗ and βj, respectively, so as λ⟶0, we have βj∗⟶βj, and further, we get Gβj∗⟶Gβj⟶Fβj. That is, Gβj∗⟶Fβj as λ⟶0.So, we conclude thatβ^j can be approximated well by βj∗^ when positive λ is small enough.Remark 1.
The reason why we substituteFβj with Gβj∗ is model (24) is easier to solve than model (22).
Now, we will show that the optimal solution to the minimization problem,(29)minβj∗yj∗−X∗βj∗2+γjPaβj∗,can be expressed as a thresholding operation.


## 3.3. The Algorithm of the Problem (29)

For anyγj,μ∈0,+∞, and z∈ℝp, let(30)Cγjβj∗=yj∗−X∗βj∗2+γjPaβj∗,Cμβj∗,z=μCγjβj∗−X∗βj∗−X∗z2+βj∗−z2,Bμβj∗=βj∗+μX∗⊤yj∗−X∗βj∗.We have the following results.Lemma 1.
For any fixed parameterμ,a,γj∈0,+∞ and z∈ℝp, ifβ¯j∗=β¯j1∗,β¯j2∗,…,β¯jp∗⊤ is a local minimizer to Cμβj∗,z, then(31)β¯ji∗=0⟺Bμzi≤t,β¯ji∗=gγjμBμzi⟺Bμzi>t,where i=1,2,…,p,(32)t=γj2a,ifγj≤1a2,γj−12a,ifγj>1a2,and according to [19],(33)ϕs=arccos27γjμa241+as3−1,gγjμs=sgns1+as/31+2cosϕs/3−π/3−1a.Proof.
We can see thatCμβj∗,z can be expressed as(34)Cμβj∗,z=βj∗−z−μX∗⊤X∗z+μX∗⊤yj∗2+γjμPaβj∗+μyj∗2+z2−μX∗z2−z−μX∗⊤X∗z+μX∗⊤yj∗2,=βj∗−Bμz2+γjμPaβj∗+μyj∗2+z2−μX∗z2−Bμz2,that is to say, minimizing Cμβj∗,z for z, given μ, γj, and a, is equivalent to(35)minβj∗∈Rp∑i=1pβji∗−Bμzi2+γjμ∑i=1ppaβji∗.
So,β¯j∗=β¯j1∗,β¯j2∗,…,β¯jp∗⊤ is a local minimizer of Cμβj∗,z if and only if β¯ji∗ solves the following problem:(36)minβji∗∈Rβji∗−Bμzi2+γjμpaβji∗.
According to Lemma 10 in [19], we complete the proof.
Furthermore, we can get Theorem1.Theorem 1.
Ifβ˜j∗=β˜j1∗,β˜j2∗,…,β˜jp∗⊤ is an optimal solution to (29), and γj is positive value and μ satisfies 0<μ≤X∗2−2, and then, the optimal solution β˜ji∗ can be given by(37)β˜ji∗=gγjμBμβ˜j∗i,ifBμβ˜j∗i>t,0,ifBμβ˜j∗i≤t,where i=1,2,…,p,(38)t=γjμa2,ifγj≤1a2μ,γjμ−12a,ifγj>1a2μ.Proof.
According to0<μ≤X∗2−2, we have(39)Cμβj∗,β˜j∗=μyj∗−X∗βj∗2+γjPaβj∗−X∗βj∗−X∗β˜j∗2+βj∗−β˜j∗2,=μyj∗−X∗βj∗2+γjPaβj∗−μX∗βj∗−β˜j∗2+βj∗−β˜j∗2≥μyj∗−X∗βj∗2+γjPaβj∗≥μyj∗−X∗β˜j∗2+γjPaβ˜j∗,=Cμβ˜j∗,β˜j∗.
That is to say, forβj∗∈Rp, β˜j∗ is a local minimizer of Cμβj∗,β˜j∗ as long as β˜j∗ is a solution to (29). According to Lemma 10 in [19] and Lemma 1, we complete the proof.
With the expression ofβ˜j∗ in Theorem 1, we can get an alternating thresholding algorithm for problem (29), that is,(40)βj∗n+1i=gγjμBμβj∗ni,for i=1,2,…,p, where Bμβj∗=βj∗+μX∗⊤yj∗−X∗βj∗ and gγjμ· is the thresholding operator in Lemma 1. We call this algorithm the alternating iterative FP-SPCA thresholding algorithm.Remark 2.
The dot· in gγjμ· may be a real number as in Lemma 1 or a vector whose components are all real numbers. If the dot · represents a vector, gγjμ· means the result vector which components are the result that gγjμ· acts on every component of the vector · sequentially.
It is well known that the solution to problem (29) depends on the choice of the regularization parameter γj. Before we give out the proper choice for the regularization parameter γj, we need the following definition.Definition 1.
The nonincreasing rearrangement of the vectorβ¯j∈ℝp is the vector β¯j∈ℝp for which(41)β¯j1≥β¯j2≥⋯≥β¯jp,and there is a permutation π:1,2,…,p⟶1,2,…,p with β¯ji=β¯jπi for all i∈1,2,…,p, and β¯jπi is the πi-th component of the vector β¯j.
Similar to the choice of the regularization parameter in FP algorithm (see Scheme 2 in [19]), here, we can approximate β˜j∗ by using the current step value βj∗n and the γj can be chosen as(42)γj=2aμBμβj∗nk+1,if2aμBμβj∗nk+1≤1a2μ,1−ε¯4a2μ2aBμβj∗nk+12,if2aμBμβj∗nk+1>1a2μ,where ε¯ is a small positive number, such as 0.1, 0.01, or 0.001, and k is sparsity degree of β˜j∗. By this way, the iterative algorithm is adaptive and free from the regular parameter choice. Noticing that 0<μ<X∗2−2 in the (2) of Theorem 2, we take μ=1−ε/X∗22, where ε∈0,1.
With abovementioned parameter-choosing method, according to the iterative method of (40), we have the following alternating iterative FP-SPCA thresholding algorithm. If A is given, for each j=1,2,…,r, yj∗ is the same as in formula (24), by using the iterative formula (40), βj∗n+1 is obtained. If B is fixed, according to (21), it is only try to minimize ∑i=1nxi−AB⊤xi2 = X−XBA⊤2 with respect to A, where A⊤A=Ir. According to the reduced rank form of the Procrustes rotation in [26], if the singular value decomposition of matrix X⊤XB is U1D1V1⊤, the solution is A^=U1V1⊤ [3]. The alternating iterative FP-SPCA thresholding algorithm is described in Algorithm 1.
At last, we discuss the convergence of the iterative FP-SPCA thresholding algorithm which converges to a stationary point of the iteration (40).Algorithm 1: Iterative FP-SPCA thresholding algorithm.
Initialize: Choose βj∗0, k, μ=1−ε/X∗22, ϵ¯, W, a, ϵ, and λ.while not converged doforj=1:rBμβj∗n=βj∗n+μX∗⊤yj∗−X∗βj∗nγj1=2/aμBμβj∗k+1; γj2n=1−ε¯/4a2μ2aBμβj∗nk+12ifγj1≤1/a2μ thenγj=γj1; t=γjμa/2elseγj=γj2; t=γjμ−1/2aendλ1,j=γj1+λ;fori=1:pifBμβj∗ni>tβj∗n+1i=gγjμBμβj∗nielseβj∗n+1i=0endendW:,j=βj∗n+1end;C:=X∗⊤X∗WU1,S1,V1=svdCV=U1V1⊤n⟶n+1end whileform=1:rβ^m=W:,m/normW:,mendoutput:β^m,m=1,2,…,r.Theorem 2.
Letβj∗n be the sequence generated by the iterative FP-SPCA thresholding algorithm (40) with 0<μ<X∗2−2. Then,(i)
The sequenceCγjβj∗n=yj∗−X∗βj∗n2+γjPaβj∗n is monotonically decreasing(ii)
βj∗n is asymptotically regular, i.e., βj∗n+1−βj∗n⟶0asn⟶∞(iii)
βj∗n converges to a stationary point of the iteration (40)Proof.
The proof is similar to the proof of Theorem 4.1 in [29] and Lemma 2 in [30], so it is omitted.


## 4. Numerical Experiments

In this section, we present numerical results of the FP-SPCA method and compare the result with the SPCA method [3]. The simulations are conducted on a personal computer (Intel(R) core(TM) i7-6700 cpu@3.40GHZ 3.40GHZ, RAM 16.0 GB) with MATLAB R2014a programming platform.


### 4.1. Experimental Preparation

In order to compare with the SPCA method, we use the synthetic data as in [3]. The synthetic data has three hidden factors,(43)V1∼N0,290,V2∼N0,300,V3=−0.3V1+0.925V2+ε,ε∼N0,1,where V1, V2, and ε are independent.The 10 observable variables are generated as follows:(44)Xi=V1+εi1,εi1∼N0,1,i=1,2,3,4,Xi=V2+εi2,εi2∼N0,1,i=5,6,7,8,Xi=V3+εi3,εi3∼N0,1,i=9,10,where εij are independent, j=1,2,3 and i=1,2,…,10. We use the matrix X=X1,X2,…,X10 to perform PCA, SPCA, and FP-SPCA.The three hidden factors have three different variance 290, 300, and 283.7875. The number of variables with respect to the three factors are 4, 4, and 2, soV1 and V2 are nearly equally important, and they are more important than V3. These facts suggest that we need to consider the relatively important variables with proper sparse representations only. The two derived variables should recover the factors V1 and V2 well using X1,X2,X3,X4 and X5,X6,X7,X8 accordingly.The SPCA and FP-SPCA methods were carried out by the method proposed by Zou et al. [3] and the fractional penalty method proposed in [19]. We compared the performance of the two methods using the synthetic data matrix X and summarize the results in two cases below.


### 4.2. Parameter Selection and the Result of Experiments

Case 1.
In this situation, the method proposed by Zou et al. [3] which is called SPCA and our proposed method which is called FP-SPCA have the same performance for the almost equally important hidden factors. That is, for most of the generated data matrices, the SPCA and FP-SPCA have the same loadings but different adjusted variance.
Table1 reports a numerical result of PCA, SPCA, and FP-SPCA and summarizes the results in loadings and adjusted variance. The parameters are λ,λ1,j,ε=0.01,3.7,0.01 in SPCA and λ,a,ε,ε¯=0.01,1.2,0.01,0.001 in FP-SPCA. We can see, from Table 1, that both SPCA and FP-SPCA have same orthogonal loadings of PC1 and PC2 after normalization, but the adjusted variance obtained from FP-SPCA is bigger than that from SPCA. That is to say, PC1 and PC2 obtained by FP-SPCA method have more information than PC1 and PC2 by SPCA. Numerical simulations using other generated data matrices have the same result also.Table 1
Loadings and variance of numerical results of PCA, SPCA, and FP-SPCA methods in Case1, where the SPCA and FP-SPCA methods have the same performance in obtaining the sparse loadings while FP-SPCA performs better than SPCA in adjusted variance. AV† (%) denotes the adjusted variance.
PCASPCAFP-SPCAPC1PC2PC3PC1PC2PC1PC2X10.37080.1993−0.26980.500.50X20.37080.1993−0.26980.500.50X30.37080.1993−0.26980.500.50X40.37080.1993−0.26980.500.50X50.2625−0.42280.04830−0.50−0.5X60.2625−0.42280.04830−0.50−0.5X70.2625−0.42280.04830−0.50−0.5X80.2625−0.42280.04830−0.50−0.5X90.29540.25100.59140000X100.29540.25100.59140000AV† (%)58.9933.107.9132.4929.8640.0036.76Case 2.
For a part of the generated matrices, according to our experimental results, the FP-SPCA does work, but SPCA does not work well. That is, by adjusting parametersλ and λ1,j in the SPCA method, we cannot get the sparse loadings of PC1 or PC2. The parameters and the corresponding numerical results are reported in Table 2, where capital letter N represents no sparse loadings of PC1 and PC2 have been obtained when parameters are the value above the N and the value(s) on the left-hand side of the N. NaN means not a number when parameters acquire the value(s) in the same way as capital letter N in Table 2. Also, we can see, from Table 2, that the loadings of PC1 or PC2 show a meaningless result when λ1,j is between 3.5 and 4.0. Numerical simulations show the same meaningless result when λ1,j is increasing by a large degree. We give, in Table 3, a result of loadings of PC1 and PC2 when λ1,j,λ,ε=1.8,0.5,0.01 using the SPCA method and the result of the FP-SPCA method when λ,a,ε,ε¯=0.001,1.3,0.01,0.001. In a word, the FP-SPCA method can find the sparse loadings of PC1 and PC2 while the SPCA method is at least difficult to obtain the sparse result in this case.
From Case1 and Case 2, it can be seen that the FP-SPCA method is more adaptable than the SPCA method. The reason why the FP-SPCA method is more preferable is that it has a parameter a in the penalty function which makes FP-SPCA more flexible and easier to adjust to get the desirable result.Table 2
The parameters selection in experiments and the results, where N represents no sparse loadings of PC1 or PC2 have been obtained and NaN means not a number as the MATLAB experimental result shows. 3.5∼4.0∗ denotes the numbers changing from the number 3.5 to 4.0 increasing by 0.1 each time.
λ1,jλ0.0010.0100.1001.0000.001NNNN0.01NNNN0.1∼3.1NNNN3.2∼3.4NNNNaN3.5∼4.0∗NaNNaNNaNNaNTable 3
Loadings and variance of numerical results of PCA, SPCA, and FP-SPCA methods in Case2, where the FP-SPCA method works better than SPCA in obtaining the sparse loadings. AV† (%) denotes the adjusted variance.
PCASPCAFP-SPCAPC1PC2PC3PC1PC2PC1PC2X1−0.3472−0.35030.0819−0.3565−0.3466−0.50X2−0.3472−0.35030.0819−0.3565−0.3466−0.50X3−0.3472−0.35030.0819−0.3565−0.3466−0.50X4−0.3472−0.35030.0819−0.3565−0.3466−0.50X5−0.34440.35660.0651−0.35050.360400.5X6−0.34440.35660.0651−0.35050.360400.5X7−0.34440.35660.0651−0.35050.360400.5X8−0.34440.35660.0651−0.35050.360400.5X9−0.1472−0.0159−0.69140000X10−0.1472−0.0159−0.69140000AV† (%)49.2332.1018.6744.6129.8840.0038.44


## 4.1. Experimental Preparation

In order to compare with the SPCA method, we use the synthetic data as in [3]. The synthetic data has three hidden factors,(43)V1∼N0,290,V2∼N0,300,V3=−0.3V1+0.925V2+ε,ε∼N0,1,where V1, V2, and ε are independent.The 10 observable variables are generated as follows:(44)Xi=V1+εi1,εi1∼N0,1,i=1,2,3,4,Xi=V2+εi2,εi2∼N0,1,i=5,6,7,8,Xi=V3+εi3,εi3∼N0,1,i=9,10,where εij are independent, j=1,2,3 and i=1,2,…,10. We use the matrix X=X1,X2,…,X10 to perform PCA, SPCA, and FP-SPCA.The three hidden factors have three different variance 290, 300, and 283.7875. The number of variables with respect to the three factors are 4, 4, and 2, soV1 and V2 are nearly equally important, and they are more important than V3. These facts suggest that we need to consider the relatively important variables with proper sparse representations only. The two derived variables should recover the factors V1 and V2 well using X1,X2,X3,X4 and X5,X6,X7,X8 accordingly.The SPCA and FP-SPCA methods were carried out by the method proposed by Zou et al. [3] and the fractional penalty method proposed in [19]. We compared the performance of the two methods using the synthetic data matrix X and summarize the results in two cases below.


## 4.2. Parameter Selection and the Result of Experiments

Case 1.
In this situation, the method proposed by Zou et al. [3] which is called SPCA and our proposed method which is called FP-SPCA have the same performance for the almost equally important hidden factors. That is, for most of the generated data matrices, the SPCA and FP-SPCA have the same loadings but different adjusted variance.
Table1 reports a numerical result of PCA, SPCA, and FP-SPCA and summarizes the results in loadings and adjusted variance. The parameters are λ,λ1,j,ε=0.01,3.7,0.01 in SPCA and λ,a,ε,ε¯=0.01,1.2,0.01,0.001 in FP-SPCA. We can see, from Table 1, that both SPCA and FP-SPCA have same orthogonal loadings of PC1 and PC2 after normalization, but the adjusted variance obtained from FP-SPCA is bigger than that from SPCA. That is to say, PC1 and PC2 obtained by FP-SPCA method have more information than PC1 and PC2 by SPCA. Numerical simulations using other generated data matrices have the same result also.Table 1
Loadings and variance of numerical results of PCA, SPCA, and FP-SPCA methods in Case1, where the SPCA and FP-SPCA methods have the same performance in obtaining the sparse loadings while FP-SPCA performs better than SPCA in adjusted variance. AV† (%) denotes the adjusted variance.
PCASPCAFP-SPCAPC1PC2PC3PC1PC2PC1PC2X10.37080.1993−0.26980.500.50X20.37080.1993−0.26980.500.50X30.37080.1993−0.26980.500.50X40.37080.1993−0.26980.500.50X50.2625−0.42280.04830−0.50−0.5X60.2625−0.42280.04830−0.50−0.5X70.2625−0.42280.04830−0.50−0.5X80.2625−0.42280.04830−0.50−0.5X90.29540.25100.59140000X100.29540.25100.59140000AV† (%)58.9933.107.9132.4929.8640.0036.76Case 2.
For a part of the generated matrices, according to our experimental results, the FP-SPCA does work, but SPCA does not work well. That is, by adjusting parametersλ and λ1,j in the SPCA method, we cannot get the sparse loadings of PC1 or PC2. The parameters and the corresponding numerical results are reported in Table 2, where capital letter N represents no sparse loadings of PC1 and PC2 have been obtained when parameters are the value above the N and the value(s) on the left-hand side of the N. NaN means not a number when parameters acquire the value(s) in the same way as capital letter N in Table 2. Also, we can see, from Table 2, that the loadings of PC1 or PC2 show a meaningless result when λ1,j is between 3.5 and 4.0. Numerical simulations show the same meaningless result when λ1,j is increasing by a large degree. We give, in Table 3, a result of loadings of PC1 and PC2 when λ1,j,λ,ε=1.8,0.5,0.01 using the SPCA method and the result of the FP-SPCA method when λ,a,ε,ε¯=0.001,1.3,0.01,0.001. In a word, the FP-SPCA method can find the sparse loadings of PC1 and PC2 while the SPCA method is at least difficult to obtain the sparse result in this case.
From Case1 and Case 2, it can be seen that the FP-SPCA method is more adaptable than the SPCA method. The reason why the FP-SPCA method is more preferable is that it has a parameter a in the penalty function which makes FP-SPCA more flexible and easier to adjust to get the desirable result.Table 2
The parameters selection in experiments and the results, where N represents no sparse loadings of PC1 or PC2 have been obtained and NaN means not a number as the MATLAB experimental result shows. 3.5∼4.0∗ denotes the numbers changing from the number 3.5 to 4.0 increasing by 0.1 each time.
λ1,jλ0.0010.0100.1001.0000.001NNNN0.01NNNN0.1∼3.1NNNN3.2∼3.4NNNNaN3.5∼4.0∗NaNNaNNaNNaNTable 3
Loadings and variance of numerical results of PCA, SPCA, and FP-SPCA methods in Case2, where the FP-SPCA method works better than SPCA in obtaining the sparse loadings. AV† (%) denotes the adjusted variance.
PCASPCAFP-SPCAPC1PC2PC3PC1PC2PC1PC2X1−0.3472−0.35030.0819−0.3565−0.3466−0.50X2−0.3472−0.35030.0819−0.3565−0.3466−0.50X3−0.3472−0.35030.0819−0.3565−0.3466−0.50X4−0.3472−0.35030.0819−0.3565−0.3466−0.50X5−0.34440.35660.0651−0.35050.360400.5X6−0.34440.35660.0651−0.35050.360400.5X7−0.34440.35660.0651−0.35050.360400.5X8−0.34440.35660.0651−0.35050.360400.5X9−0.1472−0.0159−0.69140000X10−0.1472−0.0159−0.69140000AV† (%)49.2332.1018.6744.6129.8840.0038.44


## 5. Conclusions

In this paper, the FP-SPCA method is proposed based on the SPCA method [3, 4], where a fractional penalty function [19] is proposed to replace the l1-norm in the SPCA method. Numerical simulations show that the proposed FP-SPCA method is more adaptable and flexible than the SPCA method. However, the FP-SPCA has its limitation, and during the numerical simulations, we find a few generated matrices that FP-SPCA fails to produce sparse loadings of principal components and the SPCA method does not work also. How to deal with this problem and how to compare the method we proposed with other existing penalty methods are remaining questions that will be worth investigating further.



---

*Source: 7874140-2020-08-19.xml*