A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis

Abstract

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₁-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₁-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.

1. INTRODUCTION

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 ℝⁿ 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:

1. Applying PMD to a data matrix can yield interpretable factors that provide insight into the data.

2. Applying PMD to a data matrix with L₁-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).

3. 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^TY 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.

2. THE PMD

2.1 General form of PMD

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 kth 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-rn × p matrices and ‖·‖_F² 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₁ and P₂ are convex penalty functions, which can take on a variety of forms. Useful examples areOnly certain ranges of c₁ and c₂ 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.

• lasso: P₁(u) = ∑_{i = 1}ⁿ|u_i| and

• fused lasso: P₁(u) = ∑_{i = 1}ⁿ|u_i| + λ∑_{i = 2}ⁿ|u_i − u_{i − 1}|, where λ > 0

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^TXv. The objective function u^TXv 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₂-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‖₂² = 1, provided that c₁ is chosen so that (for fixed v) the u that maximizes 

(2.9)

has L₂-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₁ 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 modelIn Section 2.2, we present an algorithm for obtaining multiple-factor solutions for the PMD. When P₁ and P₂ both are L₁-penalties, maximizations in Steps 2(a) and 2(b) are simple. This is explained in Algorithm 3 in Section 2.3.

1. Initialize v to have L₂-norm 1. 2.

2. Iterate until convergence:

   • (a) u←argmax_uu^TXv subject to P₁(u) ≤ c₁ and ‖u‖₂² ≤ 1.

   • (b) v←argmax_vu^TXv subject to P₂(v) ≤ c₂ and ‖v‖₂² ≤ 1.

3. d←u^TXv.

It can be seen that without the P₁- and P₂-constraints, the algorithm above leads to the usual rank-1 SVD. Starting with v⁽⁰⁾, 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^TX, 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.

2.2 PMD for multiple factors

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 PMDWithout the P₁- and P₂-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₁ and/or P₂ 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.

1. Let X¹ ← X.

2. 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_ku_kv_k^T.

2.3 Forms of PMD of special interest

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

(2.11)

This method results in factors u and v that are sparse for c₁ and c₂ chosen appropriately. As shown in Figure 1, we restrict c₁ and c₂ to the ranges ⁠.

Fig. 1.

A graphical representation of the L₁- and L₂-constraints on u in the PMD(L₁, L₁) criterion. The constraints are as follows: ‖u‖₂² ≤ 1 and ‖u‖₁ ≤ c. Here, u is two-dimensional, and the grey lines indicate the coordinate axes u₁ and u₂. Left: the L₂-constraint is the solid circle. For both the L₁- and L₂-constraints to be active, c must be between 1 and ⁠. The constraints ‖u‖₁ = 1 and ‖u‖₁= are shown using dashed lines. Right: The L₂- and L₁-constraints on u are shown for some c between 1 and ⁠. Small circles indicate the points where both the L₁- and the L₂-constraints are active. The solid arcs indicate the solutions that occur when Δ₁ = 0 in Algorithm 3. The figure shows that in 2D, the points where both the L₁- and L₂-constraints are active do not have either u₁ or u₂ equal to 0.

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‖₁ ≤ c; otherwise, Δ is chosen so that ‖u‖₁ = 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₁, L₁) modelIf 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, Δ₁ and Δ₂ are chosen by a binary search.

1. Initialize v to have L₂-norm 1.

2. Iterate until convergence:

   • (a) ⁠, where Δ₁ = 0 if this results in ‖u‖₁ ≤ c₁; otherwise, Δ₁ is chosen to be a positive constant such that ‖u‖₁ = c₁.

   • (b) ⁠, where Δ₂ = 0 if this results in ‖v‖₁ ≤ c₂; otherwise, Δ₂ is chosen to be a positive constant such that ‖v‖₁ = c₂.

3. d ← u^TXv.

Figure 1 shows a graphical representation of the L₁- and L₂-constraints on u that are present in the PMD(L₁, L₁) criterion: namely, ‖u‖₂² ≤ 1 and ‖u‖₁ ≤ c₁. From the figure, it is clear that in two dimensions, the intersection of the L₁- and L₂-constraints results in both u₁ and u₂ 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₁, 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₁, FL) modelStep 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).

1. Initialize v to have L₂-norm 1.

2. Iterate until convergence:

   • (a) ⁠, where Δ = 0 if this results in ‖u‖₁ ≤ c; otherwise, Δ is chosen to be a positive constant such that ‖u‖₁ = c.

   • (b) ⁠.

3. d ← u^TXv.

2.4 PMD for missing data and choice of c₁ and c₂

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₁ and c₂ in the PMD criterion. We can treat c₁ and c₂ 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 PMDNote 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).

1. From the original data matrix X, construct 10 data matrices X₁,…,X₁₀, each of which is missing a nonoverlapping one-tenth of the elements of X, sampled at random from the rows and columns.

2. For each candidate value of c₁ and c₂:

   • (a) For i ∈ 1, …, 10:

     • i) Fit the PMD to X_i with tuning parameters c₁ and c₂ and calculate _i = duv^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₁,…,X₁₀ for tuning parameters c₁ and c₂.

3. The optimal values of c₁ and c₂ are those which correspond to the lowest mean squared error.

Though c₁ and c₂ 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₁ and c₂ based on some other criterion. For instance, one could select small values of c₁ and c₂ in order to obtain factors that have a desirable level of sparsity.

2.5 Relationship between PMD and other matrix decompositions

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₁-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.

2.6 An example: PMD for DNA copy number data

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 λ₁, λ₂ ≥ 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₁, 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₁, FL) in this context is related to ongoing work by Nowak and others (2009). One could apply PMD(L₁, 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₁, 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₁, L₁) to PMD(L₁, 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).

Fig. 2.

Simulated CGH data. Top: results of PMD(L₁, FL); middle: results of PMD(L₁, L₁); bottom: generative model. PMD(L₁, FL) successfully identifies both the region of gain and the subset of samples for which that region is present.

3. SPARSE PCA VIA PMD

3.1 Three methods for sparse PCA

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:

1. 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 solves 

   

   (3.1)

   where λ₁, λ₂ ≤ 0 and v and θ are p-vectors. The criterion can equivalently be written with an inequality L₂ bound on θ, in which case it is biconvex in θ and v.
2. SCoTLASS: The SCoTLASS procedure of Jolliffe and others (2003) uses the maximal variance characterization for principal components. The first sparse principal component solves the problem 

   

   (3.2)

   and 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. SPC: We propose a new method for sparse PCA. Consider the PMD criterion with P₂(v) = ‖v‖₁, and no P₁-constraint on u. We call this criterion PMD(·, L₁), and it can be written as follows: 

   

   (3.3)

   The algorithm for PMD(·, L₁) is obtained by replacing Step 2(a) of Algorithm 3 (the single-factor PMD(L₁, L₁) 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₁, …, v_{k − 1}).

It is also not hard to show that PMD applied to a covariance matrix with symmetric L₁-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 Xv and Xu.) 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 ‖θ‖₂² = 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 ‖θ‖₁ ≤ c. We maximize 

(3.10)

with respect to θ and v. Note that for any vectors w and z, ‖z − w‖₂² ≥ 0. This means that w^Tw ≥ 2w^Tz − z^Tz. Let w = Xv and z = Xθ; it follows that θ^TX^TXθ ≥ 2v^TX^TXθ − v^TX^TXv. 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₁-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 − Xvθ^T‖_F²) 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^TX_k), where X_k = XV_k(V_k^TV_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.

Fig. 3.

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‖₁. This is very close in spirit to PMD(·, L₁), 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.

3.2 Another option for SPC with multiple factors

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₁,…,v_{k − 1}. It is not obvious that SPC can be extended to achieve orthogonality among subsequent v_is, 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₁,…,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‖₂ = ‖θ‖₂. 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_ks. Though we have not guaranteed that the v_ks will be exactly orthogonal, they are unlikely to be very correlated since the different v_ks each are associated with orthogonal u_ks. Based on our initial investigation, this appears to be a promising path for obtaining multiple sparse principal components.

4. PENALIZED CCA VIA PMD

4.1 PMD and other approaches to penalized CCA

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(Xu,Yv)—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^TX and Y^TY; this gives what could be called “diagonal penalized CCA”: 

(4.3)

Of course, this criterion is simply (2.7) with X replaced with X^TY; it can be solved with Algorithm 1. But in practice, it can be solved more efficiently, without computation of X^TY. 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₁, L₁) yields sparse vectors u and v for c₁ and c₂ 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₁, L₁). 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₁-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₁-penalties are used for both P₁ and P₂, then the values of c₁ and c₂ can be chosen by cross-validation, where c₁ and c₂ are chosen using a grid search to maximize (across the cross-validation folds) the quantity cor(Xu,Yv), where u and v are computed on a training set and X and Y constitute an independent test set. Alternatively, values of c₁ and c₂ can simply be chosen to result in desired amounts of sparsity of u and v.

4.2 Sparse CCA applied to simulated data

We demonstrate the PMD(L₁, L₁) 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^TY. Compared to the SVD, PMD(L₁, L₁) does fairly well at identifying linear combinations of the underlying factors. Details of the simulation are given in Section A. of the Appendix.

Fig. 4.

The efficacy of PMD(L₁, L₁) is demonstrated using a simulation in which Xis generated from 2 sparse latent factors, called u₁and u₂, and Yis generated from 2 sparse latent factors, called v₁and v₂. The PMD(L₁, L₁) method identifies linear combinations of these sparse factors. Details of the simulation are given in Section A. of the  Appendix.

4.3 Application of penalized CCA to genomic data

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₁-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₁, 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₁, 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₁, 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₁, FL) once for each of the 23 chromosomes, we obtained 23 v vectors. The 23 vs are shown in the left panel of Figure 5. Nonzero us and vs were found for all chromosomes except for chromosome 2. It is clear that PMD(L₁, FL) resulted in both sparsity and smoothness of the v vectors.

Fig. 5.

PMD(L₁, FL) was performed for the breast cancer data set. Left: for each chromosome, the weights of vobtained using PMD(L₁, FL) are shown. All the vweights shown are positive, but the results would not be affected by flipping the signs of both vand u. On chromosome 2, vhas no nonzero elements. Right: for each chromosome, uand vwere computed on a training set consisting of 3/4 of the samples, and cor(Xu,Yv)is plotted, where Xand Yare the training (dashed) and test (solid) data.

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₁, 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₁, FL) using the CGH spots on chromosome k and only genes that are not on chromosome k.

Table 1.

PMD(L₁, 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_is are located on chromosome 1. Similar results hold for the other chromosomes

i	Gene	Chromosome	ui
1	Jumping translocation break point	1	0.039
2	Translocated promoter region (to activated MET oncogene)	1	0.153
3	Glyceronephosphate O-acyltransferase	1	0.255
4	NADH dehydrogenase (ubiquinone) Fe-S protein 2	1	0.265
5	Nucleoporin 133kD	1	0.007
6	Geranylgeranyl diphosphate synthase 1	1	0.131
7	Rab3 GTPase-activating protein, noncatalytic subunit (150 kD)	1	0.283
8	Peroxisomal biogenesis factor 11B	1	0.154
9	Phosphatidylinositol glycan, class C	1	0.124
10	Tubulin-specific chaperone e	1	0.069
11	Protoporphyrinogen oxidase	1	0.052
12	Tuftelin 1	1	0.037
13	Papillary renal cell carcinoma (translocation associated)	1	0.055
14	Splicing factor 3b, subunit 4, 49 kD	1	0.469
15	UDP-Gal:betaGlcNAc beta 1,4- galactosyltransferase	1	0.27
16	Hypothetical protein FLJ12671	1	0.229
17	Hypothetical protein HSPC155	1	0.168
18	Mitochondrial ribosomal protein L24	1	0.195
19	HSPC003 protein	1	0.391
20	Hypothetical protein FLJ10876	1	0.091
21	CGI-78 protein	1	0.154
22	Chromosome 1 open reading frame 27	1	0.133
23	Hypothetical protein My014	1	0.278

In order to assess whether PMD(L₁, 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:With the exception of chromosome 2, all p-values were significant.

1. Let u, v denote the penalized canonical variates found for this chromosome, and record c = cor(Xu,Yv).

2. 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^*,Yv^*).

3. The p-value is given by ⁠.

In order to further assess the penalized canonical variates that we obtained, we used a training set/test set approach, as follows:The right panel of Figure 5 shows the average values of cor(X_tru_tr,Y_trv_tr) and cor(X_teu_tr,Y_tev_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_teu_tr,Y_tev_tr) would be 0.

1. 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.

2. Penalized CCA was performed on the training set to obtain the vectors u_tr and v_tr.

3. cor(X_tru_tr,Y_trv_tr) and cor(X_teu_tr,Y_tev_tr) were computed.

5. DISCUSSION

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₁-penalty on the rows and columns of our decomposition. We also explore the use of an L₁-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.

FUNDING

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.

APPENDIX

A.1 Proof of Theorem 2.1

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)

A.2 Proof of Lemma 2.2

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‖₂² − 1) = 0 and Δ(‖u‖₁ − c₁) = 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‖₂ = 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‖₁ = c₁. So, Δ = 0 if this results in ‖u‖₁ ≤ c₁; otherwise, we choose Δ such that ‖u‖₁ = c₁. This completes the proof of the Lemma.

A.3 Simulation details for Figure 2

Let X be a 12 × 1000 matrix of data. The elements of X are generated as follows:In other words, the first 5 patients have a region of gain between positions 100 and 500.

1. For i ∈ 1, …, 5 and j ∈ 100, …, 500, X_ij ∼ N(1,1).

2. Otherwise, X_ij ∼ N(0,1).

A.4 Simulation details for Figure 4

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

1. Let u₁ be a vector of length p, with 20 1s, 20 − 1s, and 60 0s.

2. Let u₂ be a vector of length p, with 10 − 1s, 10 1s, 10 − 1s, 10 1s, and 60 0s.

3. Let v₁ be a vector of length p, with 60 0s, 20 − 1s, and 20 1s.

4. Let v₂ be a vector of length p, with 60 0s, 10 1s, 10 − 1s, 10 1s, and 10 − 1s.

5. Let w₁ and w₂ be orthonormal vectors of length n.

6. Generate the data matrices as follows: X_ij ∼ N(w_1iu_1j + w_2iu_2j,0.3²) and Y_ij ∼ N(w_1iv_1j + w_2iv_2j,0.3²).