Sparse Principal Component Analysis

Hui Zou∗, Trevor Hastie†, Robert Tibshirani‡

April 26, 2004

Abstract

Principal component analysis (PCA) is widely used in data processing and dimensionality

reduction. However, PCA suffers from the fact that each principal component is a linear combi-

nation of all the original variables, thus it is often difficult to interpret the results. We introduce

a new method called sparse principal component analysis (SPCA) using the lasso (elastic net)

to produce modified principal components with sparse loadings. We show that PCA can be

formulated as a regression-type optimization problem, then sparse loadings are obtained by im-

posing the lasso (elastic net) constraint on the regression coefficients. Efficient algorithms are

proposed to realize SPCA for both regular multivariate data and gene expression arrays. We

also give a new formula to compute the total variance of modified principal components. As

illustrations, SPCA is applied to real and simulated data, and the results are encouraging.

Keywords: multivariate analysis, gene expression arrays, elastic net, lasso, singular value de-

composition, thresholding

∗
Hui Zou is a Ph.D student in the Department of Statistics at Stanford University, Stanford, CA 94305. Email:
hzou@[Link].
†
Trevor Hastie is Professor, Department of Statistics and Department of Health Research & Policy, Stanford
University, Stanford, CA 94305. Email: hastie@[Link].
‡
Robert Tibshirani is Professor, Department of Health Research & Policy and Department of Statistics, Stanford
University, Stanford, CA 94305. Email: tibs@[Link].

1
1 Introduction

Principal component analysis (PCA) (Jolliffe 1986) is a popular data processing and dimension

reduction technique . As an un-supervised learning method, PCA has numerous applications such

as handwritten zip code classification (Hastie et al. 2001) and human face recognition (Hancock

et al. 1996). Recently PCA has been used in gene expression data analysis (Misra et al. 2002).

Hastie et al. (2000) propose the so-called Gene Shaving techniques using PCA to cluster high

variable and coherent genes in microarray data.

PCA seeks the linear combinations of the original variables such that the derived variables

capture maximal variance. PCA can be done via the singular value decomposition (SVD) of the

data matrix. In detail, let the data X be a n × p matrix, where n and p are the number of

observations and the number of variables, respectively. Without loss of generality, assume the

column means of X are all 0. Suppose we have the SVD of X as

X = UDVT (1)

where T means transpose. U are the principal components (PCs) of unit length, and the columns

of V are the corresponding loadings of the principal components. The variance of the ith PC is

2 . In gene expression data the PCs U are called the eigen-arrays and V are the eigen-genes
Di,i

(Alter et al. 2000). Usually the first q (q " p) PCs are chosen to represent the data, thus a great

dimensionality reduction is achieved.

The success of PCA is due to the following two important optimal properties:

1. principal components sequentially capture the maximum variability among X, thus guaran-

teeing minimal information loss;

2. principal components are uncorrelated, so we can talk about one principal component without

2
 referring to others.

However, PCA also has an obvious drawback, i.e., each PC is a linear combination of all p variables

and the loadings are typically nonzero. This makes it often difficult to interpret the derived PCs.

Rotation techniques are commonly used to help practitioners to interpret principal components

(Jolliffe 1995). Vines (2000) considered simple principal components by restricting the loadings to

take values from a small set of allowable integers such as 0, 1 and -1.

We feel it is desirable not only to achieve the dimensionality reduction but also to reduce the

size of explicitly used variables. An ad hoc way is to artificially set the loadings with absolute

values smaller than a threshold to zero. This informal thresholding approach is frequently used in

practice but can be potentially misleading in various respects (Cadima & Jolliffe 1995). McCabe

(1984) presented an alternative to PCA which found a subset of principal variables. Jolliffe & Uddin

(2003) introduced SCoTLASS to get modified principal components with possible zero loadings.

Recall the same interpretation issue arising in multiple linear regression, where the response is

predicted by a linear combination of the predictors. Interpretable models are obtained via variable

selection. The lasso (Tibshirani 1996) is a promising variable selection technique, simultaneously

producing accurate and sparse models. Zou & Hastie (2003) propose the elastic net, a generalization

of the lasso, to further improve upon the lasso. In this paper we introduce a new approach to get

modified PCs with sparse loadings, which we call sparse principal component analysis (SPCA).

SPCA is built on the fact that PCA can be written as a regression-type optimization problem, thus

the lasso (elastic net) can be directly integrated into the regression criterion such that the resulting

modified PCA produces sparse loadings.

In the next section we briefly review the lasso and the elastic net. The method details of SPCA

are presented in Section 3. We first discuss a direct sparse approximation approach via the elas-

tic net, which is a useful exploratory tool. We then show that finding the loadings of principal

3
components can be reformulated as estimating coefficients in a regression-type optimization prob-

lem. Thus by imposing the lasso (elastic net) constraint on the coefficients, we derive the modified

principal components with sparse loadings. An efficient algorithm is proposed to realize SPCA.

We also give a new formula, which justifies the correlation effects, to calculate the total variance

of modified principal components. In Section 4 we consider a special case of the SPCA algorithm

to efficiently handle gene expression arrays. The proposed methodology is illustrated by using real

data and simulation examples in Section 5. Discussions are in Section 6. The paper ends up with

an appendix summarizing technical details.

2 The Lasso and The Elastic Net

Consider the linear regression model. Suppose the data set has n observations with p predictors.

Let Y = (y1 , . . . , yn )T be the response and Xj = (x1j , . . . , xnj )T , i = 1, . . . , p are the predictors.

After a location transformation we can assume all Xj and Y are centered.

The lasso is a penalized least squares method, imposing a constraint on the L 1 norm of the

regression coefficients. Thus the lasso estimates β̂lasso are obtained by minimizing the lasso criterion

! !2
! "p ! "p
! !
β̂lasso !
= arg min !Y − !
Xj βj ! + λ |βj | , (2)
β ! !
j=1 j=1

where λ is a non-negative value. The lasso was originally solved by quadratic programming

(Tibshirani 1996). Efron et al. (2004) proved that the lasso estimates as a function of λ are piece-

wise linear, and proposed an algorithm called LARS to efficiently solve the whole lasso solution

path in the same order of computations as a single least squares fit.

The lasso continuously shrinks the coefficients toward zero, thus gaining its prediction accuracy

via the bias variance trade-off. Moreover, due to the nature of the L1 penalty, some coefficients

4
will be shrunk to exact zero if λ1 is large enough. Therefore the lasso simultaneously produces an

accurate and sparse model, which makes it a favorable variable selection method. However, the

lasso has several limitations as pointed out in Zou & Hastie (2003). The most relevant one to this

work is that the number of selected variables by the lasso is limited by the number of observations.

For example, if applied to the microarray data where there are thousands of predictors (genes)

(p > 1000) with less than 100 samples (n < 100), the lasso can only select at most n genes, which

is clearly unsatisfactory.

The elastic net (Zou & Hastie 2003) generalizes the lasso to overcome its drawbacks, while

enjoying the similar optimal properties. For any non-negative λ1 and λ2 , the elastic net estimates

β̂en are given as follows

! !2
! p ! p p
! " ! " "
β̂en !
= (1 + λ2 ) arg min !Y − !
Xj βj ! + λ2 2
|βj | + λ1 |βj | . (3)
β ! !
j=1 j=1 j=1

Hence the elastic net penalty is a convex combination of ridge penalty and the lasso penalty .

Obviously, the lasso is a special case of the elastic net with λ2 = 0. Given a fixed λ2 , the LARS-

EN algorithm (Zou & Hastie 2003) efficiently solves the elastic net problem for all λ 1 with the

computation cost as a single least squares fit. When p > n, we choose some λ 2 > 0. Then the

elastic net can potentially include all variables in the fitted model, so the limitation of the lasso is

removed. An additional benefit offered by the elastic net is its grouping effect, that is, the elastic

net tends to select a group of highly correlated variables once one variable among them is selected.

In contrast, the lasso tends to select only one out of the grouped variables and does not care which

one is in the final model. Zou & Hastie (2003) compare the elastic net with the lasso and discuss

the application of the elastic net as a gene selection method in microarray analysis.

5
3 Motivation and Method Details

In both lasso and elastic net, the sparse coefficients are a direct consequence of the L 1 penalty,

not depending on the squared error loss function. Jolliffe & Uddin (2003) proposed SCoTLASS

by directly putting the L1 constraint in PCA to get sparse loadings. SCoTLASS successively

maximizes the variance

aTk (XT X)ak (4)

subject to

aTk ak = 1 and (for k ≥ 2) aTh ak = 0, h < k; (5)

and the extra constraints

p
"
|ak,j | ≤ t (6)
j=1

for some tuning parameter t. Although sufficiently small t yields some exact zero loadings, SCoT-

LASS seems to lack of a guidance to choose an appropriate t value. One might try several t values,

but the high computational cost of SCoTLASS makes it an impractical solution. The high compu-

tational cost is due to the fact that SCoTLASS is not a convex optimization problem. Moreover,

the examples in Jolliffe & Uddin (2003) show that the obtained loadings by SCoTLASS are not

sparse enough when requiring a high percentage of explained variance.

We consider a different approach to modify PCA, which can more directly make good use of

the lasso. In light of the success of the lasso (elastic net) in regression, we state our strategy

We seek a regression optimization framework in which PCA is done exactly. In addition,

the regression framework should allow a direct modification by using the lasso (elastic

net) penalty such that the derived loadings are sparse.

6
3.1 Direct sparse approximations

We first discuss a simple regression approach to PCA. Observe that each PC is a linear combination

of the p variables, thus its loadings can be recovered by regressing the PC on the p variables.

Theorem 1 ∀i, denote Yi = Ui Di . Yi is the i-th principal component. ∀ λ > 0, suppose β̂ridge is

the ridge estimates given by

β̂ridge = arg min |Yi − Xβ|2 + λ |β|2 . (7)

β

β̂ridge
Let v̂ = , then v̂ = Vi .
|β̂ridge |

The theme of this simple theorem is to show the connection between PCA and a regression

method is possible. Regressing PCs on variables was discussed in Cadima & Jolliffe (1995), where

they focused on approximating PCs by a subset of k variables. We extend it to a more general

ridge regression in order to handle all kinds of data, especially the gene expression data. Obviously

when n > p and X is a full rank matrix, the theorem does not require a positive λ. Note that if

p > n and λ = 0, ordinary multiple regression has no unique solution that is exactly V i . The same

story happens when n > p and X is not a full rank matrix. However, PCA always gives a unique

solution in all situations. As shown in theorem 1, this discrepancy is eliminated by the positive

ridge penalty (λ |β|2 ). Note that after normalization the coefficients are independent of λ, therefore

the ridge penalty is not used to penalize the regression coefficients but to ensure the reconstruction

of principal components. Hence we keep the ridge penalty term throughout this paper.

Now let us add the L1 penalty to (7) and consider the following optimization problem

β̂ = arg min |Yi − Xβ|2 + λ |β|2 + λ1 |β|1 . (8)

β

7
 β̂
We call V̂i = an approximation to Vi , and XV̂i the ith approximated principal component. (8)
|β̂ |
is called naive elastic net (Zou & Hastie 2003) which differs from the elastic net by a scaling factor

(1 + λ). Since we are using the normalized fitted coefficients, the scaling factor does not affect V̂i .

Clearly, large enough λ1 gives a sparse β̂, hence a sparse V̂i . Given a fixed λ, (8) is efficiently solved

for all λ1 by using the LARS-EN algorithm (Zou & Hastie 2003). Thus we can flexibly choose a

sparse approximation to the ith principal component.

3.2 Sparse principal components based on SPCA criterion

Theorem 1 depends on the results of PCA, so it is not a genuine alternative. However, it can be

used in a two-stage exploratory analysis: first perform PCA, then use (8) to find suitable sparse

approximations.

We now present a “self-contained” regression-type criterion to derive PCs. We first consider

the leading principal component.

Theorem 2 Let Xi denote the ith row vector of the matrix X. For any λ > 0, let

n
" ! !
(α̂, β̂) = arg min !Xi − αβ T Xi !2 + λ |β|2 (9)
α,β
i=1

subject to |α|2 = 1.

Then β̂ ∝ V1 .

The next theorem extends theorem 2 to derive the whole sequence of PCs.

Theorem 3 Suppose we are considering the first k principal components. Let α and β be p × k

matrices. Xi denote the i-th row vector of the matrix X. For any λ > 0, let

n
" k
"
! !
(α̂, β̂) = arg min !Xi − αβ T Xi !2 + λ |βj |2 (10)
α,β
i=1 j=1

8
 subject to α T α = Ik .

Then β̂i ∝ Vi for i = 1, 2, . . . , k.

Theorem 3 effectively transforms the PCA problem to a regression-type problem. The critical el-
#n !! !2
! #n !! !2
!
ement is the object function i=1 Xi − αβ Xi . If we restrict β = α, then
T
i=1 Xi − αβ Xi =
T

#n !! !2
i=1 Xi − ααT Xi ! , whose minimizer under the orthonormal constraint on α is exactly the first k

loading vectors of ordinary PCA. This is actually an alternative derivation of PCA other than the

maximizing variance approach, e.g. Hastie et al. (2001). Theorem 3 shows that we can still have

exact PCA while relaxing the restriction β = α and adding the ridge penalty term. As can be seen

later, these generalizations enable us to flexibly modify PCA.

To obtain sparse loadings, we add the lasso penalty into the criterion (10) and consider the

following optimization problem

n
" k
" k
"
! !
(α̂, β̂) = arg min !Xi − αβ T Xi !2 + λ |βj |2 + λ1,j |βj |
α,β 1
i=1 j=1 j=1

(11)

subject to αT α = Ik .

Whereas the same λ is used for all k components, different λ1,j s are allowed for penalizing the

loadings of different principal components. Again, if p > n, a positive λ is required in order to get

exact PCA when the sparsity constraint (the lasso penalty) vanishes (λ1,j = 0). (11) is called the

SPCA criterion hereafter.

9
3.3 Numerical solution

We propose an alternatively minimization algorithm to minimize the SPCA criterion. From the

proof of theorem 3 (see appendix for details) we get

#n !! !2 # #k
i=1Xi − αβ T Xi ! + λ k j=1 |βj |
2
+ j=1 λ1,j |βj |1

#k $ %
= TrXT X + j=1 βjT (XT X + λ)βj − 2αjT XT Xβj + λ1,j |βj |1 . (12)

Hence if given α, it amounts to solve k independent elastic net problems to get β̂j for j = 1, 2, . . . , k.

On the other hand, we also have (details in appendix)

#n !! !2 # #k
i=1Xi − αβ T Xi ! + λ k j=1 |βj |
2
+ j=1 λ1,j |βj |1

#k
= TrXT X − 2TrαT XT Xβ + Trβ T (XT X + λ)β + j=1 λ1,j |βj |1 . (13)

Thus if β is fixed, we should maximize TrαT (XT X)β subject to αT α = Ik , whose solution is given

by the following theorem.

Theorem 4 Let α and β be m × k matrices and β has rank k. Consider the constrained maxi-

mization problem
& '
α̂ = arg max Tr αT β subject to α T α = Ik . (14)
α

Suppose the SVD of β is β = U DV T , then α̂ = U V T .

Here are the steps of our numerical algorithm to derive the first k sparse PCs.

General SPCA Algorithm

1. Let α start at V[, 1 : k], the loadings of first k ordinary principal components.

10
2. Given fixed α, solve the following naive elastic net problem for j = 1, 2, . . . , k

βj = arg min
∗
β ∗T (XT X + λ)β ∗ − 2αjT XT Xβ ∗ + λ1,j |β ∗ |1 . (15)
β

3. For each fixed β, do the SVD of XT Xβ = U DV T , then update α = U V T .

4. Repeat steps 2-3, until β converges.

βj
5. Normalization: V̂j = |βj | , j = 1, . . . , k.

Some remarks:

1. Empirical evidence indicates that the outputs of the above algorithm vary slowly as λ changes.

For n > p data, the default choice of λ can be zero. Practically λ is a small positive number

to overcome potential collinearity problems of X. Section 4 discusses the default choice of λ

for the data with thousands of variables, such as gene expression arrays.

2. In principle, we can try several combinations of {λ1,j } to figure out a good choice of the

tunning parameters, since the above algorithm converges quite fast. There is a shortcut

provided by the direct sparse approximation (8). The LARS-EN algorithm efficiently deliver

a whole sequence of sparse approximations for each PC and the corresponding values of λ 1,j .

Hence we can pick a λ1,j which gives a good compromise between variance and sparsity. In

this selection, variance has a higher priority than sparsity, thus we tend to be conservative in

pursuing sparsity.

XT X
3. Both PCA and SPCA depend on X only through XT X. Note that n is actually the

sample covariance matrix of variables (Xi ). Therefore if Σ, the covariance matrix of (Xi ), is

known, we can replace XT X with Σ and have a population version of PCA or SPCA. If X is

11
 standardized beforehand, then PCA or SPCA uses the (sample) correlation matrix, which is

preferred when the scales of the variables are different.

3.4 Adjusted total variance

The ordinary principal components are uncorrelated and their loadings are orthogonal. Let Σ̂ =

XT X, then VT V = Ik and VT Σ̂V is diagonal. It is easy to check that only the loadings of

ordinary principal components can satisfy both conditions. In Jolliffe & Uddin (2003) the loadings

were forced to be orthogonal, so the uncorrelated property was sacrificed. SPCA does not explicitly

impose the uncorrelated components condition too.

Let Û be the modified PCs. Usually the total variance explained by Û is calculated by

trace(ÛT Û). This is unquestionable when Û are uncorrelated. However, if they are correlated, the

computed total variance is too optimistic. Here we propose a new formula to compute the total

variance explained by Û, which takes into account the correlations among Û.

Suppose (Ûi , i = 1, 2, . . . , k) are the first k modified PCs by any method. Denote Ûj·1,...,j−1 the

reminder of Ûj after adjusting the effects of Û1 , . . . , Ûj−1 , that is

Ûj·1,...,j−1 = Ûj − H1,...,j−1 Ûj , (16)

where H1,...,j−1 is the projection matrix on Ûi i = 1, 2, . . . , j − 1. Then the adjusted variance
! !2 #k !! !2
! ! !
of Ûj is !Ûj·1,...,j−1 ! , and the total explained variance is given by j=1 !Ûj·1,...,j−1 ! . When the

modified PCs Û are uncorrelated, then the new formula agrees with trace(ÛT Û). Note that the

above computations depend on the order of Ûi . However, since we have a natural order in PCA,

ordering is not an issue here.

Using the QR decomposition, we can easily compute the adjusted variance. Suppose Û = QR,

12
where Q is orthonormal and R is upper triangular. Then it is straightforward to see that

! !2
! !
!Ûj·1,...,j−1 ! = R2j,j . (17)

#k
Hence the explained total variance is equal to j=1 Rj,j .
2

3.5 Computation complexity

PCA is computationally efficient for both n > p or p ( n data. We separately discuss the

computational cost of the general SPCA algorithm for n > p and p ( n.

1. n > p. Traditional multivariate data fit in this category. Note that although the SPCA

criterion is defined using X, it only depends on X via XT X. A trick is to first compute the

p × p matrix Σ̂ = XT X once for all, which requires np2 operations. Then the same Σ̂ is

used at each step within the loop. Computing XT Xβ costs p2 k and the SVD of XT Xβ is of

order O(pk 2 ). Each elastic net solution requires at most O(p3 ) operations. Since k ≤ p, the

total computation cost is at most np2 + mO(p3 ), where m is the number of iterations before

convergence. Therefore the SPCA algorithm is able to efficiently handle data with huge n, as

long as p is small (say p < 100).

2. p ( n. Gene expression arrays are typical examples of this p ( n category. The trick of Σ̂

is no longer applicable, because Σ̂ is a huge matrix (p × p) in this case. The most consuming

step is solving each elastic net, whose cost is of order O(pJ 2 ) for a positive finite λ, where J is

the number of nonzero coefficients. Generally speaking the total cost is of order mO(pJ 2 k),

which is expensive for a large J. Fortunately, as shown in Section 4, there exits a special

SPCA algorithm for efficiently dealing with p ( n data.

13
4 SPCA for p ! n and Gene Expression Arrays

Gene expression arrays are a new type of data where the number of variables (genes) are much

bigger than the number of samples. Our general SPCA algorithm still fits this situation using a

positive λ. However the computation cost is expensive when requiring a large number of nonzero

loadings. It is desirable to simplify the general SPCA algorithm to boost the computation.

Observe that theorem 3 is valid for all λ > 0, so in principle we can use any positive λ. It turns

out that a thrifty solution emerges if λ → ∞. Precisely, we have the following theorem.

β̂i
Theorem 5 Let V̂i (λ) = be the loadings derived from criterion (11). Define (α̂∗ , β̂ ∗ ) as the
|β̂i |
solution of the optimization problem

k
" k
"
∗ ∗
(α̂ , β̂ ) = arg min −2Trα X Xβ + T T
βj2 + λ1,j |βj |1 (18)
α,β
j=1 j=1

subject to α T α = Ik .

β̂i∗
When λ → ∞, V̂i (λ) → .
|β̂i∗ |

By the same statements in Section 3.3, criterion (18) is solved by the following algorithm, which

is a special case of the general SPCA algorithm with λ = ∞.

Gene Expression Arrays SPCA Algorithm

Replacing step 2 in the general SPCA algorithm with

Step 2∗ : Given fixed α, for j = 1, 2, . . . , k

( )
! T T ! λ1,j
! !
βj = αj X X − Sign(αjT XT X). (19)
2 +

The operation in (19) is called soft-thresholding. Figure 1 gives an illustration of how the

soft-thresholding rule operates. Recently soft-thresholding has become increasingly popular in

14
 y

x
(0,0)

Figure 1: An illstration of soft-thresholding rule y = (|x| − ∆)+ Sign(x) with ∆ = 1.

the literature. For example, nearest shrunken centroids (Tibshirani et al. 2002) adopts the soft-

thresholding rule to simultaneously classify samples and select important genes in microarrays.

5 Examples

5.1 Pitprops data

The pitprops data first introduced in Jeffers (1967) has 180 observations and 13 measured variables.

It is the classic example showing the difficulty of interpreting principal components. Jeffers (1967)

tried to interpret the first 6 PCs. Jolliffe & Uddin (2003) used their SCoTLASS to find the modified

PCs. Table 1 presents the results of PCA, while Table 2 presents the modified PCs loadings by

SCoTLASS and the adjusted variance computed using (17).

As a demonstration, we also considered the first 6 principal components. Since this is a usual

15
n ( p data set, we set λ = 0. λ1 = (0.06, 0.16, 0.1, 0.5, 0.5, 0.5) were chosen according to Figure 2

such that each sparse approximation explained almost the same amount of variance as the ordinary

PC did. Table 3 shows the obtained sparse loadings and the corresponding adjusted variance.

Compared with the modified PCs by SCoTLASS, PCs by SPCA account for nearly the same

amount of variance (75.8% vs. 78.2%) but with a much sparser loading structure. The important

variables associated with the 6 PCs do not overlap, which further makes the interpretations easier

and clearer. It is interesting to note that in Table 3 even though the variance does not strictly

monotonously decrease, the adjusted variance follows the right order. However, Table 2 shows this

is not true in SCoTLASS. It is also worthy to mention that the whole computation of SPCA was

done in seconds in R, while the implementation of SCoTLASS for each t was expensive (Jolliffe &

Uddin 2003). Optimizing SCoTLASS over several values of t is even a more difficult computational

challenge.

Although the informal thresholding method, which is referred to as simple thresholding hence-

forth, has various drawbacks, it may serve as the benchmark for testing sparse PCs methods. An

variant of simple thresholding is soft-thresholding. We found that used in PCA, soft-thresholding

performs very similarly to simple thresholding. Thus we omitted the results of soft-thresholding

in this paper. Both SCoTLASS and SPCA were compared with simple thresholding. Table 4

presents the loadings and the corresponding explained variance by simple thresholding. To make

fair comparisons, we let the numbers of nonzero loadings by simple thresholding match the results

of SCoTLASS and SPCA. In terms of variance, it seems that simple thresholding is better than

SCoTLASS and worse than SPCA. Moreover, the variables with non-zero loadings by SPCA are

very different to that chosen by simple thresholding for the first three PCs; while SCoTLASS seems

to create a similar sparseness pattern as simple thresholding does, especially in the leading PC.

16
 PC 1 PC 2

0.30

0.15
0.25
0.20

0.10
PEV

PEV
0.15
0.10

0.05
0.05
0.00

0.00
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.5 1.0 1.5 2.0 2.5

λ1 λ1

PC 3 PC 4
0.15

0.08
0.06
0.10
PEV

PEV

0.04
0.05

0.02
0.00

0.00

0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5

λ1 λ1

PC 5 PC 6
0.08

0.08
0.06

0.06
0.04

0.04
PEV

PEV
0.02

0.02
0.00

0.00

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

λ1 λ1

Figure 2: Pitprops data: The sequences of sparse approximations to the first 6 principal components.
Plots show the percentage of explained variance (PEV) as a function of λ1 .

17
 Table 1: Pitprops data: loadings of the first 6 principal components
Variable PC1 PC2 PC3 PC4 PC5 PC6
topdiam -0.404 0.218 -0.207 0.091 -0.083 0.120
length -0.406 0.186 -0.235 0.103 -0.113 0.163
moist -0.124 0.541 0.141 -0.078 0.350 -0.276
testsg -0.173 0.456 0.352 -0.055 0.356 -0.054
ovensg -0.057 -0.170 0.481 -0.049 0.176 0.626
ringtop -0.284 -0.014 0.475 0.063 -0.316 0.052
ringbut -0.400 -0.190 0.253 0.065 -0.215 0.003
bowmax -0.294 -0.189 -0.243 -0.286 0.185 -0.055
bowdist -0.357 0.017 -0.208 -0.097 -0.106 0.034
whorls -0.379 -0.248 -0.119 0.205 0.156 -0.173
clear 0.011 0.205 -0.070 -0.804 -0.343 0.175
knots 0.115 0.343 0.092 0.301 -0.600 -0.170
diaknot 0.113 0.309 -0.326 0.303 0.080 0.626
Variance (%) 32.4 18.3 14.4 8.5 7.0 6.3
Cumulative Variance (%) 32.4 50.7 65.1 73.6 80.6 86.9

Table 2: Pitprops data: loadings of the first 6 modified PCs by SCoTLASS

t = 1.75
Variable PC1 PC2 PC3 PC4 PC5 PC6
topdiam 0.546 0.047 -0.087 0.066 -0.046 0.000
length 0.568 0.000 -0.076 0.117 -0.081 0.000
moist 0.000 0.641 -0.187 -0.127 0.009 0.017
testsg 0.000 0.641 0.000 -0.139 0.000 0.000
ovensg 0.000 0.000 0.457 0.000 -0.614 -0.562
ringtop 0.000 0.356 0.348 0.000 0.000 -0.045
ringbut 0.279 0.000 0.325 0.000 0.000 0.000
bowmax 0.132 -0.007 0.000 -0.589 0.000 0.000
bowdist 0.376 0.000 0.000 0.000 0.000 0.065
whorls 0.376 -0.065 0.000 -0.067 0.189 -0.065
clear 0.000 0.000 0.000 0.000 -0.659 0.725
knots 0.000 0.206 0.000 0.771 0.040 0.003
diaknot 0.000 0.000 -0.718 0.013 -0.379 -0.384
Number of nonzero loadings 6 7 7 8 8 8
Variance (%) 27.2 16.4 14.8 9.4 7.1 7.9
Adjusted Variance (%) 27.2 15.3 14.4 7.1 6.7 7.5
Cumulative Adjusted Variance (%) 27.2 42.5 56.9 64.0 70.7 78.2

18
 Table 3: Pitprops data: loadings of the first 6 sparse PCs by SPCA
Variable PC1 PC2 PC3 PC4 PC5 PC6
topdiam -0.477 0.000 0.000 0 0 0
length -0.476 0.000 0.000 0 0 0
moist 0.000 0.785 0.000 0 0 0
testsg 0.000 0.620 0.000 0 0 0
ovensg 0.177 0.000 0.640 0 0 0
ringtop 0.000 0.000 0.589 0 0 0
ringbut -0.250 0.000 0.492 0 0 0
bowmax -0.344 -0.021 0.000 0 0 0
bowdist -0.416 0.000 0.000 0 0 0
whorls -0.400 0.000 0.000 0 0 0
clear 0.000 0.000 0.000 -1 0 0
knots 0.000 0.013 0.000 0 -1 0
diaknot 0.000 0.000 -0.015 0 0 1
Number of nonzero loadings 7 4 4 1 1 1
Variance (%) 28.0 14.4 15.0 7.7 7.7 7.7
Adjusted Variance (%) 28.0 14.0 13.3 7.4 6.8 6.2
Cumulative Adjusted Variance (%) 28.0 42.0 55.3 62.7 69.5 75.8

5.2 A simulation example

We first created three hidden factors

V1 ∼ N (0, 290), V2 ∼ N (0, 300)

V3 = −0.3V1 + 0.925V2 + $, $ ∼ N (0, 1)

V1 , V2 and $ are independent.

Then 10 observed variables were generated as the follows

Xi = V1 + $1i , $1i ∼ N (0, 1), i = 1, 2, 3, 4,

Xi = V2 + $2i , $2i ∼ N (0, 1), i = 5, 6, 7, 8,

Xi = V3 + $3i , $3i ∼ N (0, 1), i = 9, 10,

{$ji } are independent, j = 1, 2, 3 i = 1, · · · , 10.

19
 Table 4: Pitprops data: loadings of the first 6 modified PCs by simple thresholding
Variable PC1 PC2 PC3 PC4 PC5 PC6
topdiam -0.439 0.234 0.000 0.092 0.000 0.120
length -0.441 0.000 -0.253 0.104 0.000 0.164
moist 0.000 0.582 0.000 0.000 0.361 -0.277
testsg 0.000 0.490 0.379 0.000 0.367 0.000
ovensg 0.000 0.000 0.517 0.000 0.182 0.629
ringtop 0.000 0.000 0.511 0.000 -0.326 0.000
ringbut -0.435 0.000 0.272 0.000 -0.222 0.000
bowmax -0.319 0.000 -0.261 -0.288 0.191 0.000
bowdist -0.388 0.000 0.000 -0.098 0.000 0.000
whorls -0.412 -0.267 0.000 0.207 0.000 -0.174
clear 0.000 0.221 0.000 -0.812 -0.354 0.176
knots 0.000 0.369 0.000 0.304 -0.620 -0.171
diaknot 0.000 0.332 -0.350 0.306 0.000 0.629
Number of nonzero loadings 6 7 7 8 8 8
Variance (%) 28.9 16.6 14.2 8.6 6.9 6.3
Adjusted Variance (%) 28.9 16.5 14.0 8.5 6.7 6.2
Cumulative Adjusted Variance (%) 28.9 45.4 59.4 67.9 74.6 80.8
Variable PC1 PC2 PC3 PC4 PC5 PC6
topdiam -0.420 0.000 0.000 0 0 0
length -0.422 0.000 0.000 0 0 0
moist 0.000 0.640 0.000 0 0 0
testsg 0.000 0.540 0.425 0 0 0
ovensg 0.000 0.000 0.580 0 0 0
ringtop -0.296 0.000 0.573 0 0 0
ringbut -0.416 0.000 0.000 0 0 0
bowmax -0.305 0.000 0.000 0 0 0
bowdist -0.370 0.000 0.000 0 0 0
whorls -0.394 0.000 0.000 0 0 0
clear 0.000 0.000 0.000 -1 0 0
knots 0.000 0.406 0.000 0 -1 0
diaknot 0.000 0.365 -0.393 0 0 1
Number of nonzero loadings 7 4 4 1 1 1
Variance (%) 30.7 14.8 13.6 7.7 7.7 7.7
Adjusted Variance (%) 30.7 14.7 11.1 7.6 5.2 3.6
Cumulative Adjusted Variance (%) 30.7 45.4 56.5 64.1 68.3 71.9

20
To avoid the simulation randomness, we used the exact covariance matrix of (X 1 , . . . , X10 ) to

perform PCA, SPCA and simple thresholding. In other words, we compared their performances

using an infinity amount of data generated from the above model.

The variance of the three underlying factors is 290, 300 and 283.8, respectively. The numbers of

variables associated with the three factors are 4, 4 and 2. Therefore V 2 and V1 are almost equally

important, and they are much more important than V3 . The first two PCs together explain 99.6%

of the total variance. These facts suggest that we only need to consider two derived variables with

‘right’ sparse representations. Ideally, the first derived variable should recover the factor V 2 only

using (X5 , X6 , X7 , X8 ), and the second derived variable should recover the factor V1 only using

(X1 , X2 , X3 , X4 ). In fact, if we sequentially maximize the variance of the first two derived variables

under the orthonormal constraint, while restricting the numbers of nonzero loadings to four, then

the first derived variable uniformly assigns nonzero loadings on (X5 , X6 , X7 , X8 ); and the second

derived variable uniformly assigns nonzero loadings on (X1 , X2 , X3 , X4 ).

Both SPCA (λ = 0) and simple thresholding were carried out by using the oracle information

that the ideal sparse representations use only four variables. Table 5 summarizes the comparison

results. Clearly, SPCA correctly identifies the sets of important variables. As a matter of fact,

SPCA delivers the ideal sparse representations of the first two principal components. Mathemat-

ically, it is easy to show that if t = 2 is used, SCoTLASS is also able to find the same sparse

solution. In this example, both SPCA and SCoTLASS produce the ideal sparse PCs, which may

be explained by the fact that both methods explicitly use the lasso penalty.

In contrast, simple thresholding wrongly includes X9 , X10 in the most important variables. The

explained variance by simple thresholding is also lower than that by SPCA, although the relative

difference is small (less than 5%). Due to the high correlation between V2 and V3 , variables X9 , X10

gain loadings which are even higher than that of the true important varaibles (X 5 , X6 , X7 , X8 ). Thus

21
 Table 5: Results of the simulation example: loadings and variance
PCA SPCA (λ = 0) Simple Thresholding
PC1 PC2 PC3 PC1 PC2 PC1 PC2
X1 0.116 -0.478 -0.087 0.0 0.5 0.000 -0.5
X2 0.116 -0.478 -0.087 0.0 0.5 0.000 -0.5
X3 0.116 -0.478 -0.087 0.0 0.5 0.000 -0.5
X4 0.116 -0.478 -0.087 0.0 0.5 0.000 -0.5
X5 -0.395 -0.145 0.270 0.5 0.0 0.000 0.0
X6 -0.395 -0.145 0.270 0.5 0.0 0.000 0.0
X7 -0.395 -0.145 0.270 0.5 0.0 -0.497 0.0
X8 -0.395 -0.145 0.270 0.5 0.0 -0.497 0.0
X9 -0.401 0.010 -0.582 0.0 0.0 -0.503 0.0
X10 -0.401 0.010 -0.582 0.0 0.0 -0.503 0.0
Adjusted
Variance (%) 60.0 39.6 0.08 40.9 39.5 38.8 38.6

the truth is disguised by the high correlation. On the other hand, simple thresholding correctly

discovers the second factor, because V1 has a low correlation with V3 .

5.3 Ramaswamy data

Ramaswamy data (Ramaswamy et al. 2001) has 16063 (p = 16063) genes and 144 (n = 144) samples.

Its first principal component explains 46% of the total variance. In a typical microarray data like

this, it appears that SCoTLASS cannot be practically useful. We applied SPCA (λ = ∞) to find

the sparse leading PC. A sequence of λ1 were used such that the number of nonzero loadings varied

in a rather wide range. As displayed in Figure 3, the percentage of explained variance decreases at

a slow rate, as the sparsity increase. As few as 2.5% of these 16063 genes can sufficiently construct

the leading principal component with little loss of explained variance (from 46% to 40%). Simple

thresholding was also applied to this data. It seems that when using the same number of genes,

simple thresholding always explains slightly higher variance than SPCA does. Among the same

number of selected genes by SPCA and simple thresholding, there are about 2% different genes,

and this difference rate is quite consistent.

22
 Ramaswamy data

0.46
0.45
0.44
0.43
0.42
PEV

0.41
0.40

SPCA
simple thresholding
0.37

200 1600 4000 8000 16063

number of nonzero loadings

Figure 3: The sparse leading principal component: percentage of explained variance versus sparsity.
Simple thresholding and SPCA have similar performances. However, there still exists consistent
difference in the selected genes (the ones with nonzero loadings).

6 Discussion

It has been a long standing interest to have a formal approach to derive principal components with

sparse loadings. From a practical point of view, a good method to achieve the sparseness goal

should (at least) possess the following properties.

• Without any sparsity constraint, the method should reduce to PCA.

• It should be computationally efficient for both small p and big p data.

• It should avoid mis-identifying the important variables.

The frequently used simple thresholding is not criterion based. However, this informal ad hoc

method seems to have the first two of the good properties listed above. If the explained variance

and sparsity are the only concerns, simple thresholding is not such a bad choice, and it is extremely

convenient. We have shown that simple thresholding can work pretty well in gene expression

23
arrays. The serious problem with simple thresholding is that it can mis-identify the real important

variables. Nevertheless, simple thresholding is regarded as a benchmark for any potentially better

method.

Using the lasso constraint in PCA, SCoTLASS successfully derives sparse loadings. However,

SCoTLASS is not computationally efficient, and it lacks a good rule to pick its tunning parameter.

In addition, it is not feasible to apply SCoTLASS to gene expression arrays, while in which PCA

is a quite popular tool.

In this work we have developed SPCA using the SPCA criterion. The new SPCA criterion

gives exact PCA results when its sparsity (lasso) penalty term vanishes. SPCA allows a quite

flexible control on the sparse structure of the resulting loadings. Unified efficient algorithms have

been proposed to realize SPCA for both regular multivariate data and gene expression arrays.

As a principled procedure, SPCA enjoys advantages in several aspects, including computational

efficiency, high explained variance and ability of identifying important variables.

7 Appendix: proofs

Theorem 1 proof: Using XT X = VD2 VT and VT V = I, we have

& T '−1 T
β̂ridge = X X + λI X (XVi )
( )
D2
= V V T Vi
D2 + λI
D2
= Vi 2 i . (20)
Di + λ

24
Theorem 2 proof: Note that

n
" ! ! #n
!Xi − αβ T Xi !2 =
i=1 TrXi (I − βαT )(I − αβ T )Xi
T

i=1
#n
= i=1 Tr(I − βαT )(I − αβ T )Xi XTi

#
= Tr(I − βαT )(I − αβ T )( ni=1 Xi XTi )

= Tr(I − βαT − αβ T + βαT αβ T )XT X

= TrXT X + Trβ T XT Xβ − 2TrαT XT Xβ. (21)

Since αT XT Xβ and β T XT Xβ are both scalars, we get

#n !! !2
i=1 Xi − αβ T Xi ! + λ |β|2

= TrXT X − 2αT XT Xβ + β T (XT X + λ)β. (22)

For a fixed α, the above quantity is minimized at

& '−1 T
β = XT X + λ X Xα. (23)

Substituting (23) into (22) gives

#n !! !2
i=1 Xi − αβ T Xi ! + λ |β|2

= TrXT X − 2αT XT X(XT X + λ)−1 XT Xα. (24)

Therefore

α̂ = arg max αT XT X(XT X + λ)−1 XT Xα (25)

α

25
 subject to αT α = 1.

& '−1 T
And β̂ = XT X + λ X Xα̂.

By X = UDVT , we have

D4
XT X(XT X + λ)−1 XT X = V VT . (26)
D2 + λ

D2
Hence α̂ = sV1 with s=1 or -1. Then β̂ = s D2 +λ
1
V1 .
1

Theorem 3 proof: By the same steps in the proof of theorem 2 we derive (22) as long as α T α = Ik .

Hence we have

#n !! !2 #
i=1 Xi − αβ T Xi ! + λ k j=1 |βj |
2

= TrXT X − 2TrαT XT Xβ + Trβ T (XT X + λ)β (27)

#k $ %
= TrXT X + j=1 βjT (XT X + λ)βj − 2αjT XT Xβj (28)

& '−1
Thus given a fixed α, the above quantity is minimized at βj = XT X + λ XT Xαj for

j = 1, 2, . . . , k; or equivalently
& '−1 T
β = XT X + λ X Xα. (29)

Therefore

α̂ = arg max TrαT XT X(XT X + λ)−1 XT Xα (30)

α

subject to αT α = Ik .

This is an eigen-analysis problem whose solution is α̂j = sj Vj with sj =1 or -1 for j = 1, 2, . . . , k,

26
 D2
because the eigenvectors of XT X(XT X + λ)−1 XT X are V. Hence (29) gives β̂j = sj D2 +λ
j
Vj for
j

j = 1, 2, . . . , k.

Theorem 4 proof: By assumption β = U DV T with U T U = Ik and V V T = V T V = Ik . The

k(k+1)
constraint αT α = Ik is equivalent to 2 constraints

αiT αi = 1, i = 1, 2 . . . , k (31)

αiT αj = 0, j > i. (32)

Using Lagrangian multipliers method, we define

k
" k
" k
"
1
L=− βiT αi + λi,i (αiT αi − 1) + λi,j (αiT αj ). (33)
2
i=1 i=1 j>i

Setting ∂L
∂αi = 0 gives βi = λi,i α̂i + λi,j α̂j ; or in a matrix form β = α̂Λ, where Λi,j = λj,i . Both β

and α are full rank, so Λ is invertible and α = βΛ−1 . We have

& '
Trα̂T β = TrΛ−1 β T β = Tr Λ−1,T V D2 V T , (34)

Ik = α̂T α̂ = Λ−1,T β T βΛ−1 = Λ−1,T V D2 V T Λ−1 . (35)

Let A = V T Λ−1 V , observe

k
"
& ' & '
Tr Λ−1 V D2 V T = Tr V T Λ−1 V D2 = TrAT D2 = 2
Ajj Djj , (36)
j=1

AT D2 A = Ik . (37)

27
 Since A2jj Djj
2 ≤ 1,

k
" k
"
2
Ajj Djj ≤ Djj . (38)
j=1 j=1

−1
The “=” is taken if only if A is diagonal and Ajj = Djj . Therefore Λ−1 = V AV T = V D−1 V T ,

and α̂ = βΛ = U DV T V D−1 V T = U V T .

!
β̂i∗ β̂ ∗
Theorem 5 proof: Let β̂ ∗ = (1 + λ)β̂, then we observe V̂i (λ) = . On the other hand, β̂ = 1+λ
|β̂i∗ |
means

n !
" !2 k !
" ! k ! !
! βT ! ! βj !2 " ! βj !
∗
(α̂, β̂ ) = arg min ! ! ! ! ! !
α,β !X i − α 1 + λ X i ! + λ !1 + λ! + λ1,j !
1 + λ!
i=1 j=1 j=1 1

(39)

subject to αT α = Ik .

Then by (12), we have

#n !! βT
!2
! #k !! βj !!2 #k ! !
! βj !
i=1 ! X i − α X
1+λ i ! + λ !
j=1 1+λ ! + λ
j=1 1,j 1+λ !
!
1
$# $ T
%%
k
= TrXT X + 1
1+λ j=1 βjT X 1+λ
X+λ
βj − 2αjT XT Xβj + λ1,j |βj |1
$# $ T X+λ
% %
k
= TrXT X + 1
1+λ j=1 βjT X 1+λ βj + λ1,j |βj |1 − 2TrαT XT Xβ . (40)

k
" " k
XT X + λ
(α̂, β̂ ∗ ) = arg min −2TrαT XT Xβ + βjT βj + λ1,j |βj |1 (41)
α,β 1+λ
j=1 j=1

subject to α α = Ik .T

As λ → ∞, (41) approaches (18). Thus the conclusion follows.

28
References

Alter, O., Brown, P. & Botstein, D. (2000), ‘Singular value decomposition for genome-wide ex-

pression data processing and modeling’, Proceedings of the National Academy of Sciences

97, 10101–10106.

Cadima, J. & Jolliffe, I. (1995), ‘Loadings and correlations in the interpretation of principal com-

ponents’, Journal of Applied Statistics 22, 203–214.

Efron, B., Hastie, T., Johnstone, I. & Tibshirani, R. (2004), ‘Least angle regression’, The Annals

of Statistics 32, In press.

Hancock, P., Burton, A. & Bruce, V. (1996), ‘Face processing: human perception and principal

components analysis’, Memory and Cognition 24, 26–40.

Hastie, T., Tibshirani, R., Eisen, M., Brown, P., Ross, D., Scherf, U., Weinstein, J., Alizadeh, A.,

Staudt, L. & Botstein, D. (2000), ‘’gene shaving’ as a method for identifying distinct sets of

genes with similar expression patterns’, Genome Biology 1, 1–21.

Hastie, T., Tibshirani, R. & Friedman, J. (2001), The Elements of Statistical Learning; Data

mining, Inference and Prediction, Springer Verlag, New York.

Jeffers, J. (1967), ‘Two case studies in the application of principal component’, Applied Statistics

16, 225–236.

Jolliffe, I. (1986), Principal component analysis, Springer Verlag, New York.

Jolliffe, I. (1995), ‘Rotation of principal components: choice of normalization constraints’, Journal

of Applied Statistics 22, 29–35.

29
Jolliffe, I. T. & Uddin, M. (2003), ‘A modified principal component technique based on the lasso’,

Journal of Computational and Graphical Statistics 12, 531–547.

McCabe, G. (1984), ‘Principal variables’, Technometrics 26, 137–144.

Misra, J., Schmitt, W., Hwang, D., Hsiao, L., Gullans, S., Stephanopoulos, G. & Stephanopou-

los, G. (2002), ‘Interactive exploration of microarray gene expression patterns in a reduced

dimensional space’, Genome Research 12, 1112–1120.

Ramaswamy, S., Tamayo, P., Rifkin, R., Mukheriee, S., Yeang, C., Angelo, M., Ladd, C., Reich,

M., Latulippe, E., Mesirov, J., Poggio, T., Gerald, W., Loda, M., Lander, E. & Golub, T.

(2001), ‘Multiclass cancer diagnosis using tumor gene expression signature’, Proceedings of the

National Academy of Sciences 98, 15149–15154.

Tibshirani, R. (1996), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal

Statistical Society, Series B 58, 267–288.

Tibshirani, R., Hastie, T., Narasimhan, B. & Chu, G. (2002), ‘Diagnosis of multiple cancer types by

shrunken centroids of gene’, Proceedings of the National Academy of Sciences 99, 6567–6572.

Vines, S. (2000), ‘Simple principal components’, Applied Statistics 49, 441–451.

Zou, H. & Hastie, T. (2003), Regression shrinkage and selection via the elastic net, with applications

to microarrays, Technical report, Department of Statistics, Stanford University. Available at

[Link]

30