Home | About | Journals | Submit | Contact Us | Français |

**|**Biostatistics**|**PMC2697346

Formats

Article sections

- Abstract
- 1. INTRODUCTION
- 2. THE PMD
- 3. SPARSE PCA VIA PMD
- 4. PENALIZED CCA VIA PMD
- 5. DISCUSSION
- FUNDING
- References

Authors

Related links

Biostatistics. 2009 July; 10(3): 515–534.

Published online 2009 April 17. doi: 10.1093/biostatistics/kxp008

PMCID: PMC2697346

Daniela M. Witten^{*}

Department of Statistics, Stanford University, Stanford, CA 94305, USA

*ude.drofnats@nettiwd*

Robert Tibshirani

Department of Health Research & Policy and Department of Statistics, Stanford University, Stanford, CA 94305, USA

Department of Statistics, Stanford University, Stanford, CA 94305, USA

Received 2008 July 28; Revised 2008 December 23; Revised 2009 February 4; Accepted 2009 February 24.

Copyright © The Author 2009. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oxfordjournals.org.

This article has been cited by other articles in PMC.

We present a penalized matrix decomposition (PMD), a new framework for computing a rank-*K* approximation for a matrix. We approximate the matrix **X** as , where *d*_{k}, **u**_{k}, and **v**_{k} minimize the squared Frobenius norm of **X**, subject to penalties on **u**_{k} and **v**_{k}. This results in a regularized version of the singular value decomposition. Of particular interest is the use of *L*_{1}-penalties on **u**_{k} and **v**_{k}, which yields a decomposition of **X** using sparse vectors. We show that when the PMD is applied using an *L*_{1}-penalty on **v**_{k} but not on **u**_{k}, a method for sparse principal components results. In fact, this yields an efficient algorithm for the “SCoTLASS” proposal (Jolliffe *and others* 2003) for obtaining sparse principal components. This method is demonstrated on a publicly available gene expression data set. We also establish connections between the SCoTLASS method for sparse principal component analysis and the method of Zou *and others* (2006). In addition, we show that when the PMD is applied to a cross-products matrix, it results in a method for penalized canonical correlation analysis (CCA). We apply this penalized CCA method to simulated data and to a genomic data set consisting of gene expression and DNA copy number measurements on the same set of samples.

Consider a matrix **X** with *n* rows and *p* columns. In this paper, we present a new method for computing a rank-*K* approximation for **X**:

(1.1)

where **u**_{k} and **v**_{k} are unit vectors in ^{n} and ^{p}, respectively, and *d*_{k} are nonnegative constants. We estimate **u**_{k} and **v**_{k} subject to a penalty on their elements; as a result, we call this the “penalized matrix decomposition” (PMD) of **X**.

In this paper, we will show that this decomposition has many uses:

- Applying PMD to a data matrix can yield interpretable factors that provide insight into the data.
- Applying PMD to a data matrix with
*L*_{1}-constraints on the columns but not the rows yields an efficient algorithm for the SCoTLASS method for finding sparse principal components. This is similar to a method of Shen and Huang (2008). - Applying PMD to a cross-product matrix yields a new method for penalized CCA.

The main area of application in this paper relates to (3) above. In recent years, it has become increasingly common for biologists to perform 2 different assays on the same set of samples. For instance, both gene expression and DNA copy number measurements often are available on a set of patient samples. In this situation, an integrative analysis of both the gene expression and the copy number data is desired. If **X** and **Y** are *n* × *p* and *n* × *q* matrices with standardized columns, then PMD applied to the matrix of cross-products **X**^{T}**Y** results in an efficient method for performing penalized CCA. This method can be applied to gene expression and copy number data in order to identify sets of genes that are correlated with regions of copy number change. We will demonstrate the use of our penalized CCA method for this purpose on a publicly available breast cancer data set.

In Section 2, we present the PMD. We show that the PMD can be used to identify shared regions of gain and loss in simulated DNA copy number data. In Section 3, we use the PMD to arrive at an efficient algorithm for finding sparse principal components, and we use PMD to unify preexisting methods for sparse principal component analysis (PCA). In Section 4, we extend the PMD framework in order to develop a method for penalized CCA, and we demonstrate its use on a breast cancer data set consisting of both gene expression and DNA copy number measurements on the same set of patients. Section 5 contains the discussion.

Let **X** denote an *n* × *p* matrix of data with rank *K* ≤ min(*n*, *p*). Without loss of generality, assume that the overall mean of **X** is 0. The singular value decomposition (SVD) of the data can be written as follows:

(2.1)

Let **u**_{k} denote column *k* of **U**, let **v**_{k} denote column *k* of **V**, and note that *d*_{k} denotes the *k*th diagonal element of the diagonal matrix **D**. Then, it is a well-known fact (see e.g. Eckart and Young, 1936) that for any *r* ≤ *K*,

(2.2)

where *M*(*r*) is the set of rank-*r* *n* × *p* matrices and ‖·‖_{F}^{2} indicates the squared Frobenius norm (the sum of squared elements of the matrix). In other words, the first *r* components of the SVD give the best rank-*r* approximation to a matrix, in the sense of the Frobenius norm.

In this paper, we develop generalizations of this decomposition by imposing additional constraints on the elements of **U** and **V**. We start with a rank-1 approximation.

Consider the following optimization problem:

(2.3)

Here, *P*_{1} and *P*_{2} are convex penalty functions, which can take on a variety of forms. Useful examples are

- lasso:
*P*_{1}(**u**) = ∑_{i = 1}^{n}|*u*_{i}| and - fused lasso:
*P*_{1}(**u**) = ∑_{i = 1}^{n}|*u*_{i}| +*λ*∑_{i = 2}^{n}|*u*_{i}−*u*_{i − 1}|, where*λ*>0

Only certain ranges of *c*_{1} and *c*_{2} will lead to feasible solutions, as discussed in Section 2.3. (Note that throughout this paper, the notation ‖**u**‖_{p} indicates the *L*_{p}-norm of the vector **u**, i.e. .) We now derive a more convenient form for this criterion.

The following decomposition holds.

THEOREM 2.1

Let **U** and **V** be *n* × *K* and *p* × *K* orthogonal matrices and **D** a diagonal matrix with diagonal elements *d*_{k}. Then,

(2.4)

The theorem's proof is given in the Appendix. Hence, using the case *K* = 1, we have that the values of **u** and **v** that solve (2.3) also solve the following problem:

(2.5)

and the value of *d* solving (2.3) is **u**^{T}**X****v**. The objective function **u**^{T}**X****v** in (2.5) is bilinear in **u** and **v**: that is, with **u** fixed, it is linear in **v**, and vice versa. In fact, with **v** fixed, the criterion in (2.5) takes the following form:

(2.6)

This criterion is not convex due to the *L*_{2}-equality penalty on **u**.

We can finesse this as follows. We define the (rank-1) PMD by

(2.7)

With **v** fixed, this criterion takes the form

(2.8)

which is convex. This means that (2.7) is biconvex, and this suggests an iterative algorithm for optimizing it. Moreover, it turns out that the solution to (2.8) also satisfies ‖**u**‖_{2}^{2} = 1, provided that *c*_{1} is chosen so that (for fixed **v**) the **u** that maximizes

(2.9)

has *L*_{2}-norm greater than or equal to 1. This follows from the Karush–Kuhn–Tucker conditions in convex optimization (see, e.g. Boyd and Vandenberghe, 2004). Therefore, for *c*_{1} chosen appropriately, the solution to (2.8) solves (2.6).

The following iterative algorithm is used to optimize the criterion for the rank-1 PMD.

Algorithm 1: Computation of single-factor PMD model

- Initialize
**v**to have*L*_{2}-norm 1. 2. - Iterate until convergence:
- (a)
**u**←argmax_{u}**u**^{T}**X****v**subject to*P*_{1}(**u**) ≤*c*_{1}and ‖**u**‖_{2}^{2}≤ 1. - (b)
**v**←argmax_{v}**u**^{T}**X****v**subject to*P*_{2}(**v**) ≤*c*_{2}and ‖**v**‖_{2}^{2}≤ 1.

*d*←**u**^{T}**X****v**.

In Section 2.2, we present an algorithm for obtaining multiple-factor solutions for the PMD. When *P*_{1} and *P*_{2} both are *L*_{1}-penalties, maximizations in Steps 2(a) and 2(b) are simple. This is explained in Algorithm 3 in Section 2.3.

It can be seen that without the *P*_{1}- and *P*_{2}-constraints, the algorithm above leads to the usual rank-1 SVD. Starting with **v**^{(0)}, one can show that at the end of iteration *i*,

(2.10)

This is the well-known “power method” for computing the largest eigenvector of **X**^{T}**X**, which is the leading singular vector of **X**.

In practice, we suggest using the first right singular vector of **X** as the initial value **v**. In general, Algorithm 1 does not necessarily converge to a global optimum for (2.7); however, our empirical studies indicate that the algorithm does converge to interpretable factors for appropriate choices of the penalty terms. Note that each iteration of Step 2 in Algorithm 1 results in a decrease in the objective in (2.7).

The PMD is similar to a method of Shen and Huang (2008) for identifying sparse principal components; we will elaborate on the relationship between the 2 methods in Section 3.

In order to obtain multiple factors of the PMD, we minimize the single-factor criterion (2.7) repeatedly, each time using as the **X** matrix the residuals obtained by subtracting from the data matrix the previous factors found. The algorithm is as follows.

Algorithm 2: Computation of *K* factors of PMD

- Let
**X**^{1}←**X**. - For
*k*1, …,*K*:- (a)Find
**u**_{k},**v**_{k}, and*d*_{k}by applying the single-factor PMD algorithm (Algorithm 1) to data**X**^{k}. - (b)
**X**^{k + 1}←**X**^{k}−*d*_{k}**u**_{k}**v**_{k}^{T}.

Without the *P*_{1}- and *P*_{2}-penalty constraints, it can be shown that the *K*-factor PMD algorithm leads to the rank-*K* SVD of **X**. In particular, the successive solutions are orthogonal. This can be seen since the solutions **u**_{k} and **v**_{k} are in the column and row spaces of **X**^{k}, which has been orthogonalized with respect to **u**_{j}, **v**_{j} for *j* 1,…, *k* − 1. With *P*_{1} and/or *P*_{2} present, the solutions are no longer in the column and row spaces, and so the orthogonality does not hold. In Section 3.2, we discuss an alternative multifactor approach, in the setting where PMD is specialized to sparse principal components.

We are most interested in 2 specific forms of the PMD, which we call the “PMD(*L*_{1}, *L*_{1})” and “PMD(*L*_{1}, FL)” methods. The PMD(*L*_{1}, *L*_{1}) criterion is as follows:

(2.11)

This method results in factors **u** and **v** that are sparse for *c*_{1} and *c*_{2} chosen appropriately. As shown in Figure 1, we restrict *c*_{1} and *c*_{2} to the ranges .

A graphical representation of the *L*_{1}- and *L*_{2}-constraints on **u** in the PMD(*L*_{1}, *L*_{1}) criterion. The constraints are as follows: ‖**u**‖_{2}^{2} ≤ 1 and ‖**u**‖_{1} ≤ *c*. Here, **u** is two-dimensional, and the grey lines indicate **...**

Let *S* denote the soft thresholding operator; that is, *S*(*a*,*c*) = sgn(*a*)(|*a*| − *c*)_{+}, where *c*>0 is a constant and where *x*_{+} is defined to equal *x* if *x*>0 and 0 if *x* ≤ 0. We have the following lemma.

LEMMA 2.2

Consider the optimization problem

(2.12)

The solution satisfies , with Δ = 0 if this results in ‖**u**‖_{1} ≤ *c*; otherwise, Δ is chosen so that ‖**u**‖_{1} = *c*.

The proof is given in the Appendix. We solve the PMD criterion in (2.11) using Algorithm 1, with Steps 2(a) and 2(b) adjusted as follows.

Algorithm 3: Computation of single-factor PMD(*L*_{1}, *L*_{1}) model

- Initialize
**v**to have*L*_{2}-norm 1. - Iterate until convergence:
- (a), where Δ
_{1}= 0 if this results in ‖**u**‖_{1}≤*c*_{1}; otherwise, Δ_{1}is chosen to be a positive constant such that ‖**u**‖_{1}=*c*_{1}. - (b), where Δ
_{2}= 0 if this results in ‖**v**‖_{1}≤*c*_{2}; otherwise, Δ_{2}is chosen to be a positive constant such that ‖**v**‖_{1}=*c*_{2}.

*d*←**u**^{T}**X****v**.

If one desires that **u** and **v** be equally sparse, one can simply fix a constant *c* and set . For each update of **u** and **v**, Δ_{1} and Δ_{2} are chosen by a binary search.

Figure 1 shows a graphical representation of the *L*_{1}- and *L*_{2}-constraints on **u** that are present in the PMD(*L*_{1}, *L*_{1}) criterion: namely, ‖**u**‖_{2}^{2} ≤ 1 and ‖**u**‖_{1} ≤ *c*_{1}. From the figure, it is clear that in two dimensions, the intersection of the *L*_{1}- and *L*_{2}-constraints results in both *u*_{1} and *u*_{2} nonzero. However, when *n* = 2, the dimension of **u**, is at least 3, then the right panel of Figure 1 can be thought of as the hyperplane {*u*_{i} = 0,∀*i* > 2}. In this case, the small circles indicate regions where both constraints are active and the solution is sparse (since *u*_{i} = 0 for *i*>2).

The PMD(*L*_{1}, FL) criterion is as follows (where “FL” stands for the “fused lasso” penalty, proposed in Tibshirani *and others*, 2005):

(2.13)

This method results in **u** sparse and **v** sparse and somewhat smooth (depending on the value of *λ* ≥ 0). However, for simplicity, rather than solving (2.13), we solve a slightly different criterion which results from using the Lagrange form, rather than the bound form, of the constraints on **v**:

(2.14)

We can solve this by replacing Steps 2(a) and 2(b) in Algorithm 1 with the appropriate updates:

Algorithm 4: Computation of single-factor PMD(*L*_{1}, FL) model

- Initialize
**v**to have*L*_{2}-norm 1. - Iterate until convergence:
- (a), where Δ = 0 if this results in ‖
**u**‖_{1}≤*c*; otherwise, Δ is chosen to be a positive constant such that ‖**u**‖_{1}=*c*. - (b).

*d*←**u**^{T}**X****v**.

Step 2(b) can be performed using fast software implementing fused lasso regression, as described in Friedman *and others* (2007), Tibshirani and Wang (2008), and Hoefling (2009).

The algorithm for computing the PMD works even in the case of missing data. When some elements of the data matrix **X** are missing, those elements can simply be excluded from all computations. Let *C* denote the set of indices of nonmissing elements in **X**. The criterion is as follows:

(2.15)

The PMD can therefore be used as a method for missing data imputation. This is related to SVD-based data imputation methods proposed in the literature (see, e.g. Troyanskaya *and others*, 2001).

The possibility of computing the PMD in the presence of missing data leads to a simple and automated method for the selection of the constants *c*_{1} and *c*_{2} in the PMD criterion. We can treat *c*_{1} and *c*_{2} as tuning parameters and can take an approach similar to cross-validation in order to select their values. For simplicity, we demonstrate this method for the rank-1 case here.

Algorithm 5: Selection of tuning parameters for PMD

- From the original data matrix
**X**, construct 10 data matrices**X**_{1},…,**X**_{10}, each of which is missing a nonoverlapping one-tenth of the elements of**X**, sampled at random from the rows and columns. - For each candidate value of
*c*_{1}and*c*_{2}:- (a)For
*i*1, …, 10:- i)Fit the PMD to
**X**_{i}with tuning parameters*c*_{1}and*c*_{2}and calculate_{i}=*d***uv**^{T}, the resulting estimate of**X**_{i}. - ii)Record the mean squared error of the estimate
_{i}. This mean squared error is obtained by computing the mean of the squared differences between elements of**X**and the corresponding elements of_{i}, where the mean is taken only over elements that are missing from**X**_{i}.

- (b)Record the average mean squared error across
**X**_{1},…,**X**_{10}for tuning parameters*c*_{1}and*c*_{2}.

- The optimal values of
*c*_{1}and*c*_{2}are those which correspond to the lowest mean squared error.

Note that in Step 1 of this method, we construct each **X**_{i} by randomly removing scattered elements of the matrix **X**. That is, we are not removing entire rows of **X** or entire columns of **X**, but rather individual elements of the data matrix. Similar approaches are taken in Wold (1978) and Owen and Perry (2009).

Though *c*_{1} and *c*_{2} can always be chosen as described above, for certain applications cross-validation may not be necessary. If the PMD is applied to a data set as a descriptive method, in order to obtain an intuitive understanding of the data, then one might simply fix *c*_{1} and *c*_{2} based on some other criterion. For instance, one could select small values of *c*_{1} and *c*_{2} in order to obtain factors that have a desirable level of sparsity.

In the statistical and machine learning literature, a number of matrix decompositions have been developed. We present some of these decompositions here, as they are related to the PMD. The best known of these decompositions is the SVD, which takes the form of (2.1). The SVD has a number of interesting properties, but the vectors **u**_{k} and **v**_{k} of the SVD have (in general) no nonzero elements, and the elements may be positive or negative. These qualities result in vectors **u**_{k} and **v**_{k} that are often not interpretable.

Lee and Seung (1999), Lee and Seung (2001) developed the nonnegative matrix factorization (NNMF) in order to improve upon the interpretability of the SVD. The matrix **X** is approximated as

(2.16)

where the elements of **u**_{k} and **v**_{k} are constrained to be nonnegative. The factors **u**_{k} and **v**_{k} can be interpretable: the authors apply the NNMF to a database of faces and show that the resulting factors represent facial features. The SVD does not result in interpretable facial features.

Hoyer (2002, 2004) presents the nonnegative sparse coding (NNSC), an extension of the NNMF that results in nonnegative vectors **v**_{k} and **u**_{k}, one or both of which may be sparse. Sparsity is achieved using an *L*_{1}-penalty. Since NNSC enforces a nonnegativity constraint, the resulting vectors can be quite different from those obtained via PMD; moreover, the iterative algorithm for finding the NNSC vectors is not guaranteed to decrease the objective at each step.

Lazzeroni and Owen (2002) present the plaid model, which (in the simplest case) takes the form of (1.1). They seek *d*_{k}, **u**_{k}, **v**_{k} that minimize

(2.17)

Though the plaid model provides interpretable layers, it has the drawback that the criterion cannot be optimized exactly due to the nonconvex form of the constraints on **u**_{k} and **v**_{k}.

We now consider a simple example involving comparative genomic hybridization (CGH) data, which measures DNA copy number changes along a chromosome in cancer samples. It is known that some cancers are characterized by contiguous regions of chromosomal gain or loss. For this reason, the fused lasso criterion has been proposed as a way to denoise CGH data for a single sample (Tibshirani and Wang 2008):

(2.18)

In (2.18), **y** is a vector of length *p* corresponding to measured log copy number gain/loss, ordered along the chromosome, and is a smoothed estimate of the copy number. Note that *λ*_{1}, *λ*_{2} ≥ 0.

Now, suppose that multiple CGH samples are available. We expect some patterns of gain and loss to be shared between some of the samples, and we wish to identify those patterns and samples. Let **X** denote the data matrix; the *n* rows denote the samples and the *p* columns correspond to (ordered) CGH spots. In this case, the use of PMD(*L*_{1}, FL) is appropriate because we wish to encourage sparsity in **u** (corresponding to a subset of samples) and sparsity and smoothness in **v** (corresponding to chromosomal regions). The use of PMD(*L*_{1}, FL) in this context is related to ongoing work by Nowak *and others* (2009). One could apply PMD(*L*_{1}, FL) to all chromosomes together (making sure that smoothness in the fused lasso penalty is not required between chromosomes) or one could apply PMD(*L*_{1}, FL) to each chromosome separately.

We demonstrate this method on a simple simulated example. We simulate 12 samples, each of which consists of copy number measurements on 1000 spots on a single chromosome. Five of the 12 samples contain a region of gain from spots 100–500. In Figure 2, we compare the results of PMD(*L*_{1}, *L*_{1}) to PMD(*L*_{1}, FL). It is clear that the latter method precisely uncovers the region of gain and the set of samples in which that gained region is present. Simulation details are given in the Appendix (Section A.3).

Here, we begin with an *n* × *p* data matrix **X**, with centered columns. Several methods have been proposed for estimating sparse principal components, based on either the maximum variance property of principal components or the regression/reconstruction error property. In this section, we present 2 existing methods for sparse PCA from the literature, as well as a new method based on the PMD. We will then go on to show that these 3 methods are closely related to each other. We will use the connection between PMD and one of the other methods to develop a fast algorithm for what was previously a computationally difficult formulation for sparse PCA.

The 3 methods for sparse PCA are as follows:

- SPCA: Zou
*and others*(2006) exploit the regression/reconstruction error property of principal components in order to obtain sparse principal components. For a single component, their sparse principal components analysis (SPCA) technique solveswhere(3.1)*λ*_{1},*λ*_{2}≤ 0 and**v**and*θ*are*p*-vectors. The criterion can equivalently be written with an inequality*L*_{2}bound on*θ*, in which case it is biconvex in*θ*and**v**. - SCoTLASS: The SCoTLASS procedure of Jolliffe
*and others*(2003) uses the maximal variance characterization for principal components. The first sparse principal component solves the problemand subsequent components solve the same problem with the additional constraint that they must be orthogonal to the previous components. This problem is not convex, since a convex objective must be maximized, and the computations are difficult. Trendafilov and Jolliffe (2006) provide a projected gradient algorithm for optimizing (3.2). We will show that this criterion can be optimized much more simply by direct application of Algorithm 3 in Section 2.3.(3.2) - SPC: We propose a new method for sparse PCA. Consider the PMD criterion with
*P*_{2}(**v**) = ‖**v**‖_{1}, and no*P*_{1}-constraint on**u**. We call this criterion PMD(·,*L*_{1}), and it can be written as follows:The algorithm for PMD(·,(3.3)*L*_{1}) is obtained by replacing Step 2(a) of Algorithm 3 (the single-factor PMD(*L*_{1},*L*_{1}) algorithm) with the simpler update . We will refer to this as “sparse principal components,” or SPC.

Now, consider the SPC criterion in (3.3). It is easily shown that if **v** is fixed and we seek **u** to maximize (3.3), then the optimal **u** will be . Therefore, **v** that maximizes (3.3) also maximizes

(3.4)

We recognize (3.4) as the SCoTLASS criterion (3.2). Now, since we have a fast iterative algorithm for solving (3.3), this means that we have also developed a fast method to maximize the SCoTLASS criterion. We can extend SPC to find the first *K* sparse principal components, as in Algorithm 2. Note, however, that only the first component is the solution to the SCoTLASS criterion (since we are not enforcing the constraint that component **v**_{k} be orthogonal to components **v**_{1}, …, **v**_{k − 1}).

It is also not hard to show that PMD applied to a covariance matrix with symmetric *L*_{1}-penalties on the rows and columns, as follows,

(3.5)

results in solutions **u** = **v**. (This follows from the Cauchy–Schwarz inequality applied to vectors **X****v** and **X****u**.) As a result, these solutions solve the SCoTLASS criterion as well. This also means that SPC can be performed using the covariance matrix instead of the raw data in cases where this is more convenient (e.g. if *n* *p* or if the raw data are unavailable).

We have shown that the SPC criterion is equivalent to the SCoTLASS criterion for one component and that the fast algorithm for the former can be used to maximize the latter. It turns out that there also is a connection between the SPCA criterion and the SPC criterion. Consider a modified version of the SPCA criterion (3.1) that uses the bound form, rather than the Lagrange form, of the constraints on **v**:

(3.6)

With ‖*θ*‖_{2}^{2} = 1, we have

(3.7)

So solving (3.6) is equivalent to

(3.8)

or equivalently

(3.9)

Now, suppose we add an additional constraint to (3.6): that is, let us require also that ‖*θ*‖_{1} ≤ *c*. We maximize

(3.10)

with respect to *θ* and **v**. Note that for any vectors **w** and **z**, ‖**z** − **w**‖_{2}^{2} ≥ 0. This means that **w**^{T}**w** ≥ 2**w**^{T}**z** − **z**^{T}**z**. Let **w** = **X****v** and **z** = **X***θ*; it follows that *θ*^{T}**X**^{T}**X***θ* ≥ 2**v**^{T}**X**^{T}**X***θ* − **v**^{T}**X**^{T}**X****v**. So (3.10) is maximized when **v** = *θ*; that is, **v** that maximizes (3.10) also maximizes

(3.11)

which of course is simply the SCoTLASS criterion (3.2) again. Therefore, we have shown that if a symmetric *L*_{1}-constraint on *θ* is added to the bound form of the SPCA criterion, then the SCoTLASS criterion results. From this argument, it is also clear that the solution to the bound form of SPCA will give lower reconstruction error (defined as ‖**X** − **X****v***θ*^{T}‖_{F}^{2}) than the solution to the SCoTLASS criterion.

We compare the proportion of variance explained by SPC and SPCA on a publicly available gene expression data set from http://icbp.lbl.gov/breastcancer/, and described in Chin *and others* (2006), consisting of 19 672 gene expression measurements on 89 samples. (For consistency with Section 4.3, we use the subset of the samples for which both gene expression and CGH measurements are available.) For computational reasons, we use only the subset of the data consisting of the 5% of genes with highest variance. We compute the first 25 sparse principal components for SPC using the constraint on **v** that results in an average of 195 genes with nonzero elements per sparse component. We then perform SPCA on the same data, using tuning parameters such that each loading has the same number of nonzero elements obtained using the SPC method. Figure 3 shows the proportion of variance explained by the first *k* sparse principal components, defined as tr(**X**_{k}^{T}**X**_{k}), where **X**_{k} = **X****V**_{k}(**V**_{k}^{T}**V**_{k})^{ − 1}**V**_{k}^{T} and where **V**_{k} is the matrix that has the first *k* sparse principal components as its columns. (This definition is proposed in Shen and Huang, 2008.) SPC results in a substantially greater proportion of variance explained, as expected.

Breast cancer gene expression data: a greater proportion of variance is explained when SPC is used to obtain the sparse principal components, rather than SPCA. Multiple SPC components were obtained as described in Algorithm 2.

Our extension of PMD to the problem of identifying sparse principal components is closely related to the SPCA method of Shen and Huang (2008). They present a method for identifying sparse principal components via a regularized low-rank matrix approximation, as follows:

(3.12)

Then, is the first sparse principal component of their method. They present a number of forms for *P*_{λ}(**v**), including *P*_{λ}(**v**) = ‖**v**‖_{1}. This is very close in spirit to PMD(·, *L*_{1}), and in fact the algorithm is almost the same. Since Shen and Huang (2008) use the Lagrange form of the constraint on **v**, their formulation does not solve the SCoTLASS criterion. Our method unifies the regularized low-rank matrix approximation approach of Shen and Huang (2008) with the maximum variance criterion of Jolliffe *and others* (2003) and the SPCA method of Zou *and others* (2006).

To summarize, in our view, the SCoTLASS criterion (3.2) is the simplest, most natural way to define the notion of sparse principal components. Unfortunately, the criterion is difficult to optimize. Our SPC criterion (3.3) recasts this problem as a biconvex one, leading to an extremely simple algorithm for the solution of the first SCoTLASS component. Furthermore, the SPCA criterion (3.1) is somewhat complex. But we have shown that when a natural symmetric constraint is added to the SPCA criterion (3.1), it is also equivalent to (3.2) and (3.3). Taken as a whole, these arguments point to the SPC criterion (3.3) as the criterion of choice for this problem, at least for a single component.

As mentioned in Section 3.1, the first sparse principal component of our SPC method optimizes the SCoTLASS criterion. But subsequent sparse principal components obtained using SPC do not, since we do not enforce that **v**_{k} be orthogonal to **v**_{1},…,**v**_{k − 1}. It is not obvious that SPC can be extended to achieve orthogonality among subsequent **v**_{i}s, or even that orthogonality is desirable. However, SPC can be easily extended to give something similar to orthogonality.

Consider the criterion for the first factor of SPC, given in (3.3). One could extend to multiple factors as proposed in Algorithm 2. (This was done in Figure 3.) Alternatively, one could obtain multiple factors **u**_{k},**v**_{k} by optimizing the following criterion, for *k*>1:

(3.13)

With **u**_{k} fixed, one can solve (3.13) for **v**_{k} easily, as has been done throughout this paper (e.g. Step 2(b) of Algorithm 3). With **v**_{k} fixed, the problem is as follows: we must find **u**_{k} that maximizes

(3.14)

Let **U**_{k − 1}^{} denote an orthogonal basis that is orthogonal to **U**_{k − 1}, the matrix with columns **u**_{1},…,**u**_{k − 1}. It follows that **u**_{k} is in the column space of **U**_{k − 1}^{}, and so can be written as **u**_{k} = **U**_{k − 1}^{}*θ*. Note also that ‖**u**_{k}‖_{2} = ‖*θ*‖_{2}. So (3.14) is equivalent to solving

(3.15)

and so we find that the optimal *θ* is

(3.16)

Therefore, the value of **u**_{k} that solves (3.14) is

(3.17)

So we can use this update step for **u**_{k} to develop an iterative algorithm to find multiple factors for (3.3), the single-factor SPC criterion, that yields orthogonal **u**_{k}s. Though we have not guaranteed that the **v**_{k}s will be exactly orthogonal, they are unlikely to be very correlated since the different **v**_{k}s each are associated with orthogonal **u**_{k}s. Based on our initial investigation, this appears to be a promising path for obtaining multiple sparse principal components.

Suppose that we have *n* observations on *p* + *q* variables and that the variables are naturally partitioned into 2 sets of size *p* and *q*. Let **X** denote the *n* × *p* matrix that is comprised of the first set of variables, and let **Y** denote the *n* × *p* matrix that is comprised of the remaining variables; assume that the columns of **X** and **Y** have been centered and scaled. CCA, developed by Hotelling (1936), involves finding **u**, **v** that maximize cor(**X****u**,**Y****v**)—that is, that solve

(4.1)

There is a closed-form solution for **u** and **v** that involves the eigenvectors of some function of the covariance matrices of **X** and **Y**. We call **u** and **v** the canonical variates.

CCA results in vectors **u**, **v** that are not sparse, and these vectors are not unique if *p* or *q* exceeds *n*. In certain applications, especially if *p* or *q* is large, one might be interested in finding a linear combination of the variables in **X** and **Y** that has large correlation but is also sparse in the variables used.

One way to obtain penalized canonical variates would simply be to include penalties in (4.1):

(4.2)

It has been shown that in other high-dimensional problems, treating the covariance matrix as diagonal can yield good results (see, e.g. Dudoit *and others*, 2001; Tibshirani *and others*, 2003). For this reason, rather than using (4.2) as our penalized CCA criterion, we substitute in the identity matrix **I** for **X**^{T}**X** and **Y**^{T}**Y**; this gives what could be called “diagonal penalized CCA”:

(4.3)

Of course, this criterion is simply (2.7) with **X** replaced with **X**^{T}**Y**; it can be solved with Algorithm 1. But in practice, it can be solved more efficiently, without computation of **X**^{T}**Y**. To compute multiple canonical variates, we use Algorithm 2. Following the notation of Section 2.3, we refer to this method as PMD(*A*, *B*) if *A* is the penalty on **u** and *B* is the penalty on **v**.

PMD(*L*_{1}, *L*_{1}) yields sparse vectors **u** and **v** for *c*_{1} and *c*_{2} sufficiently small. Waaijenborg *and others* (2008) also present a sparse CCA method, and their algorithm for finding the sparse canonical variates is quite similar to PMD(*L*_{1}, *L*_{1}). However, they arrive at their algorithm in a circuitous way, as an approximation to the elastic net, and they do not state the exact criterion that they are solving. Moreover, their method involves the Lagrange form rather than the bound form of the *L*_{1}-constraints on **u** and **v**; as a result, it yields a different solution. Parkhomenko *and others* (2007) and Wiesel *and others* (2008) present sparse CCA algorithms that lack exact criteria and biconvexity, respectively. Our method for sparse CCA is very closely related to the method of Parkhomenko *and others* (2009), which we encountered after our paper was submitted.

When *L*_{1}-penalties are used for both *P*_{1} and *P*_{2}, then the values of *c*_{1} and *c*_{2} can be chosen by cross-validation, where *c*_{1} and *c*_{2} are chosen using a grid search to maximize (across the cross-validation folds) the quantity cor(**X****u**,**Y****v**), where **u** and **v** are computed on a training set and **X** and **Y** constitute an independent test set. Alternatively, values of *c*_{1} and *c*_{2} can simply be chosen to result in desired amounts of sparsity of **u** and **v**.

We demonstrate the PMD(*L*_{1}, *L*_{1}) method on a simple simulated example. In this simulation, *p*,*q*>*n*, so classical CCA cannot be used. There are 2 sparse latent factors that generate **X** and **Y**; the factors are orthogonal to each other. Results can be seen in Figure 4. For comparison, we also computed the SVD of **X**^{T}**Y**. Compared to the SVD, PMD(*L*_{1}, *L*_{1}) does fairly well at identifying linear combinations of the underlying factors. Details of the simulation are given in Section A. of the Appendix.

In genomic research, it is becoming increasingly common for researchers to use multiple assays in order to characterize a single set of samples. For instance, gene expression measurements and DNA copy number data might be available on the same set of tissue samples. The patients might also be genotyped. Examples of studies that combine gene expression and copy number and/or single-nucleotide polymorphism (SNP) data include Hyman *and others* (2002), Pollack *and others* (2002), Morley *and others* (2004), and Stranger *and others* (2005), (Stranger07). While much research has gone into developing methods for the identification of genes and SNPs that are associated with an outcome based on a single gene expression or SNP data set, the question of how to combine the results of multiple assays in order to perform inference across the data sets has not been thoroughly investigated.

If both gene expression data and genotype data are available on the same set of samples, then a natural question is to identify sets of genes that are correlated with sets of SNPs. Both Parkhomenko *and others* (2007), Parkhomenko *and others* (2009) and Waaijenborg *and others* (2008) demonstrate the use of sparse CCA for this purpose. Similarly, if CGH and gene expression data both are available for a set of cancer samples, then one might wish to identify a set of genes that have expression that is correlated with a set of chromosomal gains or losses.

In the case of gene expression and CGH data, it makes sense to perform penalized CCA with an *L*_{1}-penalty on the canonical variate corresponding to genes and a fused lasso penalty on the canonical variate corresponding to copy number—in other words, we might use PMD(*L*_{1}, FL). There are 2 ways that this could be done. Penalized CCA could be applied using all available gene expression data and copy number data on all chromosomes (being sure that the fused lasso smoothness penalty is not applied between chromosomes). Alternatively, one could perform penalized CCA once per chromosome, each time using copy number data on that chromosome and all the available gene expression data. (The expression data are not restricted to a particular chromosome.) We choose to pursue this latter approach.

We examine the performance of PMD(*L*_{1}, FL) on the breast cancer data set publicly available at http://icbp.lbl.gov/breastcancer/ and described in Chin *and others* (2006). In addition to the gene expression data described earlier, CGH measurements are also available on the same set of 89 samples. There are *p* = 19672 gene expression measurements and *q* = 2149 CGH measurements. For convenience and interpretability of the results, we ran PMD(*L*_{1}, FL) using values of the tuning parameters that resulted in a list of approximately 25 nonzero genes per chromosome (i.e. 25 nonzero elements of **u**) and very smooth and somewhat sparse **v**. Since we performed PMD(*L*_{1}, FL) once for each of the 23 chromosomes, we obtained 23 **v** vectors. The 23 **v**s are shown in the left panel of Figure 5. Nonzero **u**s and **v**s were found for all chromosomes except for chromosome 2. It is clear that PMD(*L*_{1}, FL) resulted in both sparsity and smoothness of the **v** vectors.

PMD(*L*_{1}, FL) was performed for the breast cancer data set. Left: for each chromosome, the weights of **v**obtained using PMD(*L*_{1}, FL) are shown. All the **v**weights shown are positive, but the results would not be affected by flipping the signs of both **v**and **u** **...**

The genes corresponding to nonzero weights in each **u** vector can also be examined. Consider Table 1, which shows the genes that had nonzero weights when sparse CCA was run using the CGH spots on chromosome 1 and all the available genes. Notably, only genes located on chromosome 1 were given nonzero weights. This is intuitive: a copy number change on chromosome 1 should be correlated with expression changes in the genes that were amplified or deleted. Similar results were seen when PMD(*L*_{1}, FL) was run using the CGH spots on other chromosomes. If one is interested only in discovering genes not on chromosome *k* that are correlated with copy number change on chromosome *k*, then one could perform PMD(*L*_{1}, FL) using the CGH spots on chromosome *k* and only genes that are not on chromosome *k*.

PMD(*L*_{1}, FL) was performed using the CGH spots on chromosome 1 and gene expression measurements on all chromosomes. The genes corresponding to nonzero elements of **u** are shown. Notably, all the genes with nonzero *u _{i}*s are located on chromosome 1. Similar

In order to assess whether PMD(*L*_{1}, FL) is capturing real structure in the breast cancer data, we computed *p*-values for the penalized canonical variates. For each chromosome, a *p*-value for the penalized canonical variates was computed as follows:

- Let
**u**,**v**denote the penalized canonical variates found for this chromosome, and record*c*= cor(**X****u**,**Y****v**). - For
*i*1, …,*B*, permute the samples in**X**to obtain**X**^{*}; then, compute**u**^{*}and**v**^{*}, the penalized canonical variates based on data (**X**^{*},**Y**). Record*c*_{i}= cor(**X**^{*}**u**^{*},**Y****v**^{*}). - The
*p*-value is given by .

With the exception of chromosome 2, all *p*-values were significant.

In order to further assess the penalized canonical variates that we obtained, we used a training set/test set approach, as follows:

- We repeatedly divided the 89 samples into a training set (
**X**_{tr},**Y**_{tr}) containing 3/4 of the samples, and a test set (**X**_{te},**Y**_{te}) containing the remaining samples. - Penalized CCA was performed on the training set to obtain the vectors
**u**_{tr}and**v**_{tr}. - cor(
**X**_{tr}**u**_{tr},**Y**_{tr}**v**_{tr}) and cor(**X**_{te}**u**_{tr},**Y**_{te}**v**_{tr}) were computed.

The right panel of Figure 5 shows the average values of cor(**X**_{tr}**u**_{tr},**Y**_{tr}**v**_{tr}) and cor(**X**_{te}**u**_{tr},**Y**_{te}**v**_{tr}) for each chromosome. Even on the test set, quite high correlations result for most chromosomes. In the absence of signal, the average value of cor(**X**_{te}**u**_{tr},**Y**_{te}**v**_{tr}) would be 0.

We have developed a method for finding a PMD in an efficient manner. This decomposition builds upon a variety of existing matrix decompositions, such as the SVD, the NNMF (Lee and Seung, 1999), (Lee and Seung, 2001), and the plaid model (Lazzeroni and Owen, 2002). We are most interested in obtaining a decomposition made up of sparse vectors. To do this, we use an *L*_{1}-penalty on the rows and columns of our decomposition. We also explore the use of an *L*_{1}-penalty on the rows and a fused lasso penalty on the columns; this is appropriate if the samples correspond to DNA copy number, ordered by chromosomal location. We exploit the biconvex nature of the PMD criterion in order to minimize it via an alternating algorithm.

We have applied the PMD to give attractive solutions to 2 additional problems: sparse PCA and sparse CCA. We used the resulting sparse CCA method to identify the sets of genes that are correlated with regions of DNA copy number change using a data set consisting of DNA copy number change and gene expression measurements on the same set of samples. In addition, we have established connections between 3 different methods for obtaining sparse principal components.

An R package implementing these methods, called PMA (for penalized multivariate analysis) is available on CRAN.

National Defense Science and Engineering Graduate Fellowship to D.W.; National Science Foundation (DMS-9971405 to R.T., DMS-0505676 to T.H.); National Institutes of Health (N01-HV-28183 to R.T., 2R01 CA 72028-07 to T.H.).

We thank an associate editor and 2 referees for helpful comments, Stephen Boyd for a useful discussion of biconvexity, and Holger Hoefling for the use of his R code for fused lasso regression. *Conflict of Interest:* None declared.

Let **u**_{k} and **v**_{k} denote column *k* of **U** and **V**, respectively. We prove the theorem by expanding out the squared Frobenius norm and rearranging terms:

(A.1)

We seek **u** that minimizes

(A.2)

First, we rewrite the criterion using Lagrange multipliers:

(A.3)

and we differentiate, set the derivative to 0, and solve for **u**:

(A.4)

where Γ_{i} = sgn(*u*_{i}) if *u*_{i} ≠ 0; otherwise, Γ_{i} [ − 1,1]. The Karush–Kuhn–Tucker conditions for optimality consist of (A.4), along with *λ*(‖**u**‖_{2}^{2} − 1) = 0 and Δ(‖**u**‖_{1} − *c*_{1}) = 0. Now if *λ*>0, we have

(A.5)

In general, we have either *λ* = 0 (if this results in a feasible solution) or *λ* must be chosen such that ‖**u**‖_{2} = 1. So we see that

(A.6)

Again by the Karush–Kuhn–Tucker conditions, either Δ = 0 (if this results in a feasible solution) or Δ must be chosen such that ‖**u**‖_{1} = *c*_{1}. So, Δ = 0 if this results in ‖**u**‖_{1} ≤ *c*_{1}; otherwise, we choose Δ such that ‖**u**‖_{1} = *c*_{1}. This completes the proof of the Lemma.

Let **X** be a 12 × 1000 matrix of data. The elements of **X** are generated as follows:

- For
*i*1, …, 5 and*j*100, …, 500,*X*_{ij}~*N*(1,1). - Otherwise,
*X*_{ij}~*N*(0,1).

In other words, the first 5 patients have a region of gain between positions 100 and 500.

We generate matrices **X** and **Y**, with *n* = 50 and *p* = 100.

- Let
**u**_{1}be a vector of length*p*, with 20 1s, 20 − 1s, and 60 0s. - Let
**u**_{2}be a vector of length*p*, with 10 − 1s, 10 1s, 10 − 1s, 10 1s, and 60 0s. - Let
**v**_{1}be a vector of length*p*, with 60 0s, 20 − 1s, and 20 1s. - Let
**v**_{2}be a vector of length*p*, with 60 0s, 10 1s, 10 − 1s, 10 1s, and 10 − 1s. - Let
**w**_{1}and**w**_{2}be orthonormal vectors of length*n*. - Generate the data matrices as follows:
*X*_{ij}~*N*(*w*_{1i}*u*_{1j}+*w*_{2i}*u*_{2j},0.3^{2}) and*Y*_{ij}~*N*(*w*_{1i}*v*_{1j}+*w*_{2i}*v*_{2j},0.3^{2}).

- Boyd S, Vandenberghe L. Convex Optimization. New York: Cambridge University Press; 2004.
- Chin K, DeVries S, Fridlyand J, Spellman P, Roydasgupta R, Kuo W-L, Lapuk A, Neve R, Qian Z, Ryder T. Genomic and transcriptional aberrations linked to breast cancer pathophysiologies. Cancer Cell. 2006;10:529–541.
*and others*. [PubMed] - Dudoit S, Fridlyand J, Speed T. Comparison of discrimination methods for the classification of tumors using gene expression data. Journal of the American Statistical Association. 2001;96:1151–1160.
- Eckart C, Young G. The approximation of one matrix by another of low rank. Psychometrika. 1936;1:211.
- Friedman J, Hastie T, Hoefling H, Tibshirani R. Pathwise coordinate optimization. Annals of Applied Statistics. 2007;1:302–332.
- Hoefling H. A path algorithm for the fused lasso (in preparation) 2009
- Hotelling H. Relations between two sets of variates. Biometrika. 1936;28:321–377.
- Hoyer P. Non-negative sparse coding. Proceedings of the IEEE Workshop on Neural Networks for Signal Processing. 2002:557–565.
- Hoyer P. Non-negative matrix factorization with sparseness constraints. Journal of Machine Learning Research. 2004;5:1457–1469.
- Hyman E, Kauraniemi P, Hautaniemi S, Wolf M, Mousses S, Rozenblum E, Ringner M, Sauter G, Monni O, Elkahloun A. Impact of DNA amplification on gene expression patterns in breast cancer. Cancer Research. 2002;62:6240–6245.
*and others*. [PubMed] - Jolliffe I, Trendafilov N, Uddin M. A modified principal component technique based on the lasso. Journal of Computational and Graphical Statistics. 2003;12:531–547.
- Lazzeroni L, Owen A. Plaid models for gene expression data. Statistica Sinica. 2002;12:61–86.
- Lee DD, Seung HS. Learning the parts of objects by non-negative matrix factorization. Nature. 1999;401:788. [PubMed]
- Lee DD, Seung HS. Algorithms for non-negative matrix factorization. In: Leen TK, Diatterich TG, Tresp V, editors. Advances in Neural Information Processing Systems. Cambridge, MA: MIT Press; 2001. pp. 556–562.
- Morley M, Molony C, Weber T, Devlin J, Ewens K, Spielman R, Cheung V. Genetic analysis of genome-wide variation in human gene expression. Nature. 2004;430:743–747. [PMC free article] [PubMed]
- Nowak G, Hastie T, Pollack J, Tibshirani R. Identifying copy number changes in CGH data for multiple samples (in preparation) 2009
- Owen AB, Perry PO. Bi-cross-validation of the SVD and the non-negative matrix factorization. Annals of Applied Statistics. 2009
- Parkhomenko E, Tritchler D, Beyene J. Genome-wide sparse canonical correlation of gene expression with genotypes. BMC Proceedings. 2007;1:S119. [PMC free article] [PubMed]
- Parkhomenko E, Tritchler D, Beyene J. Sparse canonical correlation analysis with application to genomic data integration. Statistical Applications in Genetics and Molecular Biology. 2009;8:1–34. [PubMed]
- Pollack J, Sorlie T, Perou C, Rees C, Jeffrey S, Lonning P, Tibshirani R, Botstein D, Borresen-Dale A, Brown P. Microarray analysis reveals a major direct role of DNA copy number alteration in the transcriptional program of human breast tumors. Proceedings of the National Academy of Sciences of the United States of America. 2002;99:12963–12968. [PubMed]
- Shen H, Huang J. Sparse principal component analysis via regularized low rank matrix approximation. Journal of Multivariate Analysis. 2008;101:1015–1034.
- Stranger B, Forrest M, Clark A, Minichiello M, Deutsch S, Lyle R, Hunt S, Kahl B, Antonarakis S, Tavare S. Genome-wide associations of gene expression variation in humans. PLoS Genetics. 2005:1–e78.
*and others*. [PubMed] - Stranger B, Forrest M, Dunning M, Ingle C, Beazley C, Thorne N, Redon R, Bird C, de Grassi A, Lee C. Relative impact of nucleotide and copy number variation on gene expression phenotypes. Science. 2007;315:848–853.
*and others*. [PMC free article] [PubMed] - Tibshirani R, Hastie T, Narasimhan B, Chu G. Class prediction by nearest shrunken centroids, with applications to DNA microarrays. Statistical Science. 2003;18:104–117.
- Tibshirani R, Saunders M, Rosset S, Zhu J, Knight K. Sparsity and smoothness via the fused lasso. Journal of the Royal Statistical Society Series B. 2005;67:91–108.
- Tibshirani R, Wang P. Spatial smoothing and hotspot detection for CGH data using the fused lasso. Biostatistics. 2008;9:18–29. [PubMed]
- Trendafilov N, Jolliffe I. Projected gradient approach to the numerical solution of the scotlass. Computational Statistics & Data Analysis. 2006;50:242–253.
- Troyanskaya O, Cantor M, Sherlock G, Brown P, Hastie T, Tibshirani R, Botstein D, Altman R. Missing value estimation methods for DNA microarrays. Bioinformatics. 2001;16:520–525. [PubMed]
- Waaijenborg S, Verselewel de Witt Hamer P, Zwinderman A. Quantifying the association between gene expressions and DNA-markers by penalized canonical correlation analysis. Statistical Applications in Genetics and Molecular Biology. 2008 7. Issue 1, Article 3. [PubMed]
- Wiesel A, Kliger M, Hero A. A greedy approach to sparse canonical correlation analysis (In preparation) Available at http://arxiv.org/abs/0801.2748. 2008
- Wold S. Cross-validatory estimation of the number of components in factor and principal components models. Technometrics. 1978;20:397–405.
- Zou H, Hastie T, Tibshirani R. Sparse principal component analysis. Journal of Computational and Graphical Statistics. 2006;15:265–286.

Articles from Biostatistics (Oxford, England) are provided here courtesy of **Oxford University Press**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |