# A recursive divide-and-conquer approach for sparse principal component analysis

## Abstract

In this paper, a new method is proposed for sparse PCA based on the recursive divide-and-conquer methodology. The main idea is to separate the original sparse PCA problem into a series of much simpler sub-problems, each having a closed-form solution. By recursively solving these sub-problems in an analytical way, an efficient algorithm is constructed to solve the sparse PCA problem. The algorithm only involves simple computations and is thus easy to implement. The proposed method can also be very easily extended to other sparse PCA problems with certain constraints, such as the nonnegative sparse PCA problem. Furthermore, we have shown that the proposed algorithm converges to a stationary point of the problem, and its computational complexity is approximately linear in both data size and dimensionality. The effectiveness of the proposed method is substantiated by extensive experiments implemented on a series of synthetic and real data in both reconstruction-error-minimization and data-variance-maximization viewpoints.

## 1Introduction

Principal component analysis (PCA) is one of the most classical and popular tools for data analysis and dimensionality reduction, and has a wide range of successful applications throughout science and engineering [?]. By seeking the so-called principal components (PCs), along which the data variance is maximally preserved, PCA can always capture the intrinsic latent structure underlying data. Such information greatly facilitates many further data processing tasks, such as feature extraction and pattern recognition.

Despite its many advantages, the conventional PCA suffers from the fact that each component is generally a linear combination of all data variables, and all weights in the linear combination, also called loadings, are typically non-zeros. In many applications, however, the original variables have meaningful physical interpretations. In biology, for example, each variable of gene expression data corresponds to a certain gene. In these cases, the derived PC loadings are always expected to be sparse (i.e. contain fewer non-zeros) so as to facilitate their interpretability. Moreover, in certain applications, such as financial asset trading, the sparsity of the PC loadings is especially expected since fewer nonzero loadings imply fewer transaction costs.

Accordingly, sparse PCA has attracted much attention in the recent decade, and a variety of methods for this topic have been developed [?]. The first attempt for this topic is to make certain post-processing transformation, e.g. rotation [?] by Jolliffe and simple thresholding [?] by Cadima and Jolliffe, on the PC loadings obtained by the conventional PCA to enforce sparsity. Jolliffe and Uddin further advanced a SCoTLASS algorithm by simultaneously calculating sparse PCs on the PCA model with additional -norm penalty on loading vectors [?]. Better results have been achieved by the SPCA algorithm of Zou et al., which was developed based on iterative elastic net regression [?]. D’Aspremont et al. proposed a method, called DSPCA, for finding sparse PCs by solving a sequence of semidefinite programming (SDP) relaxations of sparse PCA [?]. Shen and Huang developed a series of methods called sPCA-rSVD (including sPCA-rSVD, sPCA-rSVD, sPCA-rSVD), computing sparse PCs by low-rank matrix factorization under multiple sparsity-including penalties [?]. Journée et al. designed four algorithms, denoted as GPower, GPower, GPower, and GPower, respectively, for sparse PCA by formulating the issue as non-concave maximization problems with - or -norm sparsity-inducing penalties and extracting single unit sparse PC sequentially or block units ones simultaneously [?]. Based on probabilistic generative model of PCA, some methods have also been attained [?], e.g. the EMPCA method derived by Sigg and Buhmann for sparse and/or nonnegative sparse PCA [?]. Sriperumbudur et al. provided an iterative algorithm called DCPCA, where each iteration consists of solving a quadratic programming (QP) problem [?]. Recently, Lu and Zhang developed an augmented Lagrangian method (ALSPCA briefly) for sparse PCA by solving a class of non-smooth constrained optimization problems [?]. Additionally, d’Aspremont derived a PathSPCA algorithm that computes a full set of solutions for all target numbers of nonzero coefficients [?].

There are mainly two methodologies utilized by the current research on sparse PCA problem. The first is the greedy approach, including DSPCA [?], sPCA-rSVD [?], EMPCA [?], PathSPCA [?], etc. These methods mainly focus on the solving of one-sparse-PC model, and more sparse PCs can be sequentially calculated on the deflated data matrix or data covariance [?]. Under this methodology, the first several sparse PCs underlying the data can generally be properly extracted, while the computation for more sparse PCs tends to be incrementally invalidated due to the cumulation of computational error. The second is the block approach. Typical methods include SCoTLASS [?], GPower, GPower [?], ALSPCA [?], etc. These methods aim to calculate multiple sparse PCs at once by utilizing certain block optimization techniques. The block approach for sparse PCA is expected to be more efficient than the greedy one to simultaneously attain multiple PCs, while is generally difficult to elaborately rectify each individual sparse PC based on some specific requirements in practice (e.g. the number of nonzero elements in each PC).

In this paper, a new methodology, called the recursive divide-and-conquer (ReDaC briefly), is employed for solving the sparse PCA problem. The main idea is to decompose the original large and complex problem of sparse PCA into a series of small and simple sub-problems, and then recursively solve them. Each of these sub-problems has a closed-form solution, which makes the new method simple and very easy to implement. On one hand, as compared with the greedy approach, the new method is expected to integratively achieve a collection of appropriate sparse PCs of the problem by iteratively rectifying each sparse PC in a recursive way. The group of sparse PCs attained by the proposed method is further proved being a stationary solution of the original sparse PCA problem. On the other hand, as compared with the block approach, the new method can easily handle the constraints superimposed on each individual sparse PC, such as certain sparsity and/or nonnegative constraints. Besides, the computational complexity of the proposed method is approximately linear in both data size and dimensionality, which makes it well-suited to handle large-scale problems of sparse PCA.

In what follows, the main idea and the implementation details of the proposed method are first introduced in Section 2. Its convergence and computational complexity are also analyzed in this section. The effectiveness of the proposed method is comprehensively substantiated based on a series of empirical studies in Section 3. Then the paper is concluded with a summary and outlook for future research. Throughout the paper, we denote matrices, vectors and scalars by the upper-case bold-faced letters, lower-case bold-faced letters, and lower-case letters, respectively.

## 2The recursive divide-and-conquer method for sparse PCA

In the following, we first introduce the fundamental models for the sparse PCA problem.

### 2.1Basic models of sparse PCA

Denote the input data matrix as , where and are the size and the dimensionality of the given data, respectively. After a location transformation, we can assume all to have zero mean. Let be the data covariance matrix.

The classical PCA can be solved through two types of optimization models [?]. The first is constructed by finding the -dimensional linear subspace where the variance of the input data is maximized [?]. On this data-variance-maximization viewpoint, the PCA is formulated as the following optimization model:

where denotes the trace of the matrix and denotes the array of PC loading vectors. The second is formulated by seeking the -dimensional linear subspace on which the projected data and the original ones are as close as possible [?]. On this reconstruction-error-minimization viewpoint, the PCA corresponds to the following model:

where is the Frobenius norm of , is the matrix of PC loading array and is the matrix of projected data. The two models are intrinsically equivalent and can attain the same PC loading vectors [?].

Corresponding to the PCA models (Equation 1) and (Equation 2), the sparse PCA problem has the following two mathematical formulations^{1}

and

where or and the corresponding denotes the - or the -norm of , respectively. Note that the involved or penalty in the above models (Equation 3) and (Equation 4) tends to enforce sparsity of the output PCs. Methods constructed on (Equation 3) include SCoTLASS [?], DSPCA [?], DCPCA [?], ALSPCA [?], etc., and those related to (Equation 4) include SPCA [?], sPCA-rSVD [?], SPC [?], GPower [?], etc. In this paper, we will construct our method on the reconstruction-error-minimization model (Equation 4), while our experiments will verify that the proposed method also performs well based on the data-variance-maximization criterion.

### 2.2Decompose original problem into small and simple sub-problems

The objective function of the sparse PCA model (Equation 4) can be equivalently formulated as follows:

where . It is then easy to separate the original large minimization problem, which is with respect to and , into a series of small minimization problems, which are each with respect to a column vector of and of for , respectively, as follows:

and

Through recursively optimizing these small sub-problems, the recursive divide-and-conquer (ReDaC) method for solving the sparse PCA model (Equation 4) can then be naturally constructed.

It is very fortunate that both the minimization problems in (Equation 5) and (Equation 6) have closed-form solutions. This implies that the to-be-constructed ReDaC method can be fast and efficient, as presented in the following sub-sections.

### 2.3The closed-form solutions of () and ()

For the convenience of denotation, we first rewrite (Equation 5) and (Equation 6) as the following forms:

and

where is -dimensional and is -dimensional. Since the objective function can be equivalently transformed as:

(Equation 7) and (Equation 8) are equivalent to the following optimization problems, respectively:

and

The closed-form solutions of (Equation 9) and (Equation 10), i.e. (Equation 7) and (Equation 8), can then be presented as follows.

We present the closed-form solution to (Equation 8) in the following theorem.

The theorem is very easy to prove by calculating where the gradient of is equal to zero. We thus omit the proof.

In the case, the closed-form solution to (Equation 9) is presented in the following theorem. Here, we denote , and the hard thresholding function, whose -th element corresponds to , where is the -th element of and (equals if is ture , and otherwise) is the indicator function. The proof of the theorem is provided in Appendix A.

In the above theorem, denotes the empty set, implying that when , the optimum of ( ?) does not exist.

In the case, (Equation 7) has the following closed-form solution. In the theorem, we denote , where represents the soft thresholding function , where represents the vector attained by projecting to its nonnegative orthant, and denotes the permutation of based on the ascending order of .

It should be noted that we have proved that and is a monotonically decreasing function with respect to in Lemma 1 of the appendix. This means that we can conduct the optimum of the optimization problem (Equation 7) for any based on the above theorem.

The ReDaC algorithm can then be easily constructed based on Theorems 1-3.

### 2.4The recursive divide-and-conquer algorithm for sparse PCA

The main idea of the new algorithm is to recursively optimize each column, of or of for , with other s and s () fixed. The process is summarized as follows:

Update each column of for by the closed-form solution of (Equation 5) attained from Theorem 2 (for ) or Theorem 3 (for ).

Update each column of for by the closed-form solution of (Equation 6) calculated from Theorem 1.

Through implementing the above procedures iteratively, and can be recursively updated until the stopping criterion is satisfied. We summarize the aforementioned ReDaC technique as Algorithm ?.

We then briefly discuss how to specify the stopping criterion of the algorithm. The objective function of the sparse PCA model (Equation 4) is monotonically decreasing in the iterative process of Algorithm 1 since each of the step 5 and step 6 in the iterations makes an exact optimization for a column vector of or of , with all of the others fixed. We can thus terminate the iterations of the algorithm when the updating rate of or is smaller than some preset threshold, or the maximum number of iterations is reached.

Now we briefly analyze the computational complexity of the proposed ReDaC algorithm. It is evident that the computational complexity of Algorithm ? is essentially determined by the iterations between step 5 and step 6, i.e. the calculation of the closed-form solutions of and of and , respectively. To compute , only simple operations are involved and the computation needs cost. To compute , a sorting for the elements of the -dimensional vector is required, and the total computational cost is around by applying the well-known heap sorting algorithm [?]. The whole process of the algorithm thus requires around computational cost in each iteration. That is, the computational complexity of the proposed algorithm is approximately linear in both the size and the dimensionality of input data.

### 2.5Convergence analysis

In this section we evaluate the convergence of the proposed algorithm.

The convergence of our algorithm can actually be implied by the monotonic decrease of the cost function of (Equation 4) during the iterations of the algorithm. In specific, in each iteration of the algorithm, step 5 and step 6 optimize the column vector of or of , with all of the others fixed, respectively. Since the objective function of (Equation 4) is evidently lower bounded (), the algorithm is guaranteed to be convergent.

We want to go a further step to evaluate where the algorithm converges. Based on the formulation of the optimization problem (Equation 4), we can construct a specific function as follows:

where

and for each of , is an indicator function defined as:

It is then easy to show that the constrained optimization problem (Equation 4) is equivalent to the unconstrained problem

The proposed ReDaC algorithm can then be viewed as a block coordinate descent (BCD) method for solving (Equation 11) [?], by alteratively optimizing , , respectively. Then the following theorem implies that our algorithm can converge to a stationary point of the problem.

In the above theorem, the assumption that the function , as defined in (Equation 11), is regular holds under the condition that is open and is Gateaux-differentiable on (Lemma 3.1 under Condition A1 in [?]). Based on Theorems 1-3, we can also easily see that has unique minimum in each and with others fixed. The above theorem can then be naturally followed by Theorem 4.1(c) in [?].

Another advantage of the proposed ReDaC methodology is that it can be easily extended to other sparse PCA applications when certain constraints are needed for output sparse PCs. In the following section we give one of the extensions of our methodology — nonnegative sparse PCA problem.

### 2.6The ReDaC method for nonnegative sparse PCA

The nonnegative sparse PCA [?] problem differs from the conventional sparse PCA in its nonnegativity constraint imposed on the output sparse PCs. The nonnegativity property of this problem is especially important in some applications such as microeconomics, environmental science, biology, etc. [?]. The corresponding optimization model is written as follows:

where means that each element of is greater than or equal to 0.

By utilizing the similar recursive divide-and-conquer strategy, this problem can be separated into a series of small minimization problems, each with respect to a column vector of and of for , respectively, as follows:

and

where or . Since (Equation 14) is of the same formulation as (Equation 6), we only need to discuss how to solve (Equation 13). For the convenience of denotation, we first rewrite (Equation 13) as:

The closed-form solution of (Equation 15) is given in the following theorem.

By virtue of the closed-form solution of (Equation 15) given by Theorem 5, we can now construct the ReDaC algorithm for solving nonnegative sparse PCA model (Equation 12). Since the algorithm differs from Algorithm 1 only in step 5 (i.e. updating of ), we only list this step in Algorithm 2.

We then substantiate the effectiveness of the proposed ReDaC algorithms for sparse PCA and nonnegative sparse PCA through experiments in the next section.

## 3Experiments

To evaluate the performance of the proposed ReDaC algorithm on the sparse PCA problem, we conduct experiments on a series of synthetic and real data sets. All the experiments are implemented on Matlab 7.11(R2010b) platform in a PC with AMD Athlon(TM) 64 X2 Dual 5000+@2.60 GHz (CPU), 2GB (memory), and Windows XP (OS). In all experiments, the SVD method is utilized for initialization. The proposed algorithm under both and was implemented in all experiments and mostly have a similar performance. We thus only list the better one throughout.

### 3.1Synthetic simulations

Two synthetic data sets are first utilized to evaluate the performance of the proposed algorithm on recovering the ground-truth sparse principal components underlying data.

#### Hastie data

Hastie data set was first proposed by Zou et al. [?] to illustrate the advantage of sparse PCA over conventional PCA on sparse PC extraction. So far this data set has become one of the most frequently utilized benchmark data for testing the effectiveness of sparse PCA methods. The data set is generated in the following way: first, three hidden factors , and are created as:

where , and , and are independent; afterwards, observable variables are generated as:

where and all s are independent. The data so generated are of intrinsic sparse PCs [?]: the first recovers the factor only using , and the second recovers only utilizing .

We generate sets of data, each contains data generated in the aforementioned way, and apply Algorithm 1 to them to extract the first two sparse PCs. The results show that our algorithm can perform well in all experiments. In specific, the proposed ReDaC algorithm faithfully delivers the ground-truth sparse PCs in all experiments. The effectiveness of the proposed algorithm is thus easily substantiated in this series of benchmark data.

#### Synthetic toy data

As [?] and [?], we adopt another interesting toy data, with intrinsic sparse PCs, to evaluate the performance of the proposed method. The data are generated from the Gaussian distribution with mean and covariance , which is calculated by

Here, , the eigenvalues of the covariance matrix , are pre-specified as , respectively, and are -dimensional orthogonal vectors, formulated by

and the rest being generated by applying Gram-Schmidt orthonormalization to randomly valued -dimensional vectors. It is easy to see that the data generated under this distribution are of first two sparse PC vectors and .

Four series of experiments, each involving sets of data generated from , are utilized, with sample sizes , , , , respectively. For each experiment, the first two PCs, and , are calculated by a sparse PCA method and then if both and are satisfied, the method is considered as a success. The proposed ReDaC method, together with the conventional PCA and current sparse PCA methods, including SPCA [?], DSPCA [?], PathSPCA [?], sPCA-rSVD, sPCA-rSVD, sPCA-rSVD [?], EMPCA [?], GPower, GPower, GPower, GPower [?] and ALSPCA [?], have been implemented, and the success times for four series of experiments have been recorded and summarized, respectively. The results are listed in Table 1.

PCA | ||||

SPCA | ||||

DSPCA | ||||

PathSPCA | ||||

sPCA-rSVD | ||||

sPCA-rSVD | ||||

sPCA-rSVD | ||||

EMPCA | ||||

GPower | 155 | 154 | 155 | 139 |

GPower | 122 | 127 | 126 | 126 |

GPower | 91 | 76 | 71 | 16 |

GPower | 90 | 92 | 88 | 82 |

ALSPCA | 669 | 749 |
826 | 927 |

ReDaC | ||||

The advantage of the proposed ReDaC algorithm can be easily observed from Table 1. In specific, our method always attains the highest or second highest success times (in the size case, less than ALSPCA) as compared with the other utilized methods in all of the four series of experiments. Considering that the ALSPCA method, which is the only comparable method in these experiments, utilizes strict constraints on the orthogonality of output PCs while the ReDaC method does not utilize any prior ground-truth information of data, the capability of the proposed method on sparse PCA calculation can be more prominently verified.

### 3.2Experiments on real data

In this section, we further evaluate the performance of the proposed ReDaC method on two real data sets, including the pitprops and colon data. Two quantitative criteria are employed for performance assessment. They are designed in the viewpoints of reconstruction-error-minimization and data-variance-maximization, respectively, just corresponding to the original formulations (Equation 4) and (Equation 3) for sparse PCA problem.

Reconstruction-error-minimization criterion: RRE.

Once sparse PC loading matrix is obtained by a method, the input data can then be reconstructed by , where , attained by the least square method. Then the relative reconstruction error (RRE) can be calculated by

to assess the performance of the utilized method in data reconstruction point of view.

Data-variance-maximization criterion: PEV.

After attaining the sparse PC loading matrix , the input data can then be reconstructed by , as aforementioned. And thus the variance of the reconstructed data can be computed by . The percentage of explained variance (PEV, [?]) of the reconstructed data from the original one can then be calculated by

to evaluate the performance of the utilized method in data variance point of view.

#### Pitprops data

The pitprops data set, consisting of observations and measured variables, was first introduced by Jeffers [?] to show the difficulty of interpreting PCs. This data set is one of the most commonly utilized examples for sparse PCA evaluation, and thus is also employed to testify the effectiveness of the proposed ReDaC method. The comparison methods include SPCA [?], DSPCA [?], PathSPCA [?], sPCA-rSVD, sPCA-rSVD, sPCA-rSVD [?], EMPCA [?], GPower, GPower, GPower, GPower [?] and ALSPCA [?]. For each utilized method, sparse PCs are extracted from the pitprops data, with different cardinality settings: 8-5-6-2-3-2 (altogether nonzero elements), 7-4-4-1-1-1 (altogether nonzero elements, as set in [?]) and 7-2-3-1-1-1 (altogether nonzero elements, as set in [?]), respectively. In each experiment, both the RRE and PEV values, as defined above, are calculated, and the results are summarized in Table 2. Figure ? further shows the the RRE and PEV curves attained by different sparse PCA methods in all experiments for more illumination. It should be noted that the GPower, GPower and ALSPCA methods employ the block methodology, as introduced in the introduction of the paper, and calculate all sparse PCs at once while cannot sequentially derive different numbers of sparse PCs with preset cardinality settings. Thus the results of these methods reported in Table 2 are calculated with the total sparse PC cardinalities being 26, 18 and 15, respectively, and are not included in Figure ?.

RRE | PEV | RRE | PEV | RRE | PEV | |

SPCA | 0.4162 | 82.68% | 0.4448 | 80.22% | 0.4459 | 80.11% |

DSPCA | 0.4303 | 81.48% | 0.4563 | 79.18% | 0.4771 | 77.23% |

PathSPCA | 0.4080 | 83.35% | 0.4660 | 80.11% | 0.4457 | 80.13% |

sPCA-rSVD | 0.4139 | 82.87% | 0.4376 | 80.85% | 0.4701 | 77.90% |

sPCA-rSVD | 0.4314 | 81.39% | 0.4427 | 80.40% | 0.4664 | 78.25% |

sPCA-rSVD | 0.4306 | 81.45% | 0.4453 | 80.17% | 0.4762 | 77.32% |

EMPCA | 0.4070 | 83.44% | 0.4376 | 80.85% | 0.4451 | 80.18% |

GPower | 0.4092 | 83.26% | 0.4400 | 80.64% | 0.4457 | 80.13% |

GPower | 0.4080 | 83.35% | 0.4460 | 80.11% | 0.4457 | 80.13% |

GPower | 0.4224 | 82.16% | 0.5089 | 74.10% | 0.4644 | 78.44% |

GPower | 0.4187 | 82.46% | 0.4711 | 77.81% | 0.4589 | 78.94% |

ALSPCA | 0.4168 | 82.63% | 0.4396 | 80.67% | 0.4537 | 79.42% |

ReDaC | 0.4005 |
83.50% |
0.4343 |
81.14% |
0.4420 |
80.46% |

It can be seen from Table 2 that under all cardinality settings of the first PCs, the proposed ReDaC method always achieves the lowest RRE and highest PEV values among all the competing methods. This means that the ReDaC method is advantageous in both reconstruction-error-minimization and data-variance-maximization viewpoints. Furthermore, from Figure ?, it is easy to see the superiority of the ReDaC method. In specific, for different number of extracted sparse PC components, the proposed ReDaC method can always get the smallest RRE values and the largest PEV values, as compared with the other utilized sparse PCA methods, in the experiments. This further substantiates the effectiveness of the proposed ReDaC method in both reconstruction-error-minimization and data-variance-maximization views.

#### Colon data

The colon data set [?] consists of tissue samples with the gene expression profiles of genes extracted from DNA micro-array data. This is a typical data set with high-dimension and low-sample-size property, and is always employed by sparse methods for extracting interpretable information from high-dimensional genes. We thus adopt this data set for evaluation. In specific, sparse PCs, each with nonzero loadings, are calculated by different sparse PCA methods, including SPCA [?], PathSPCA [?], sPCA-rSVD, sPCA-rSVD, sPCA-rSVD [?], EMPCA [?], GPower, GPower, GPower, GPower [?] and ALSPCA [?], respectively. Their performance is compared in Table 3 and Figure ? in terms of RRE and PEV, respectively. It should be noted that the DSPCA method has also been tried, while cannot be terminated in a reasonable time in this experiment, and thus we omit its result in the table. Besides, we have carefully tuned the parameters of the GPower methods (including GPower, GPower, GPower and GPower), and can get sparse PCs with total cardinality around , similar as the total nonzero elements number of the other utilized sparse PCA methods, while cannot get sparse PC loading sequences each with cardinality as expected. The results are thus not demonstrated in Figure ?.

SPCA | PathSPCA | sPCA-rSVD | sPCA-rSVD | |
---|---|---|---|---|

RRE. | 0.7892 | 0.5287 | 0.5236 | 0.5628 |

PEV. | 37.72% | 72.05% | 72.58% | 68.32% |

sPCA-rSVD | EMPCA | GPower | GPower | |

RRE. | 0.5723 | 0.5211 | 0.5042 | 0.5076 |

PEV. | 67.25% | 72.84% | 74.56% | 74.23% |

GPower | GPower | ALSPCA | ReDaC | |

RRE. | 0.4870 | 0.4904 | 0.5917 | 0.4737 |

PEV. | 76.29% | 75.95% | 64.99% | 77.56% |

From Table 3, it is easy to see that the proposed ReDaC method achieves the lowest RRE and highest PEV values, as compared with the other employed sparse PCA methods. Figure ? further demonstrates that as the number of extracted sparse PCs increases, the advantage of the ReDaC method tends to be more dominant than other methods, with respect to both the RRE and PEV criteria. This further substantiates the effectiveness of the proposed method and implies its potential usefulness in applications with various interpretable components.

### 3.3Nonnegative sparse PCA experiments

We further testify the performance of the proposed ReDaC method (Algorithm 2) in nonnegative sparse PC extraction. For comparison, two existing methods for nonnegative sparse PCA, NSPCA [?] and Nonnegative EMPCA (N-EMPCA, briefly) [?], are also employed.

#### Synthetic toy data

As the toy data utilized in Section 3.2, we also formulate a Gaussian distribution with mean and covariance matrix . Both the leading two eigenvectors of are specified as nonnegative and sparse vectors as:

and the rest are then generated by applying Gram-Schmidt orthonormalization to randomly valued -dimensional vectors. The 10 corresponding eigenvalues are preset as , respectively. Four series of experiments are designed, each with data sets generated from , with sample sizes , , and , respectively. For each experiment, the first two PCs are calculated by the conventional PCA, NSPCA, N-EMPCA and ReDaC methods, respectively. The success times, calculated in the similar way as introduced in Section ?, of each utilized method on each series of experiments are recorded, as listed in Table 4.

PCA | ||||

NSPCA | ||||

N-EMPCA | ||||

ReDaC | ||||

From Table 4, it is seen that the ReDaC method achieves the highest success rates in all experiments. The advantage of the proposed ReDaC method on nonnegative sparse PCA calculation, as compared with the other utilized methods, can thus been verified in these experiments.

#### Colon data

The colon data set is utilized again for nonnegative sparse PCA calculation. The NSPCA and N-EMPCA methods are adopted as the competing methods. Since the NSPCA method cannot directly pre-specify the cardinalities of the extracted sparse PCs, we thus first apply NSPCA on the colon data (with parameters and ) and then use the cardinalities of the nonnegative sparse PCs attained by this method to preset the N-EMPCA and ReDaC methods for fair comparison. sparse PCs are computed by the three methods, and the performance is compared in Table 5 and Figure ?, in terms of RRE and PEV, respectively.

NSPCA | N-EMPCA | ReDaC | |
---|---|---|---|

RRE | 0.3674 | 0.3399 | 0.2706 |

PEV | 86.50% | 88.45% | 92.68% |

Just as expected, it is evident that the proposed ReDaC method dominates in both RRE and PEV viewpoints. From Table 5, we can observe that our method achieves the lowest RRE and highest PEV on extracted nonnegative sparse PCs than the other two utilized methods. Furthermore, Figure ? shows that our method is advantageous, as compared with the other methods, for any preset number of extracted sparse PCs, and this advantage tends to be more significant as more sparse PCs are to be calculated. The effectiveness of the proposed method on nonnegative sparse PCA calculation can thus be further verified.

#### Application to face recognition

In this section, we introduce the performance of our method in face recognition problem [?]. The proposed ReDaC method, together with the conventional PCA, NSPCA and N-EMPCA methods, have been applied to this problem and their performance is compared in this application. The employed data set is the MIT CBCL Face Dataset #1, downloaded from “http://cbcl.mit.edu/software-datasets/FaceData2.html”. This data set consists of aligned face images and non-face images, each with resolution . For each of the four utilized methods, PC loading vectors are computed on face images, as shown in Figure ?, respectively. For easy comparison, we also list the RRE and PEV values of three nonnegative sparse PCA methods in Table 6.

NSPCA | N-EMPCA | ReDaC | |
---|---|---|---|

RRE | 0.6993 | 0.6912 | 0.6606 |

PEV | 51.10% | 52.22% | 56.36% |

As depicted in Figure ?, the nonnegative sparse PCs obtained by the ReDaC method more clearly exhibit the interpretable features underlying faces, as compared with the other utilized methods, e.g. the first five PCs calculated from our method clearly demonstrate the eyebrows, eyes, cheeks, mouth and chin of faces, respectively. The advantage of the proposed method can further be verified quantitatively by its smallest RRE and largest PEV values, among all employed methods, in the experiment, as shown in Table 6. The effectiveness of the ReDaC method can thus be substantiated.

To further show the usefulness of the proposed method, we apply it to face classification under this data set as follows. First we randomly choose face images and non-face images from MIT CBCL Face Dataset #1, and take them as the training data and the rest images as testing data. We then extract PCs by utilizing the PCA, NSPCA, N-EMPCA and ReDaC methods to the training set, respectively. By projecting the training data onto the corresponding PCs obtained by these four methods, respectively, and then fitting the linear Logistic Regression (LR) [?] model on these dimension-reduced data (-dimensional), we can get a classifier for testing. The classification accuracy of the classifier so obtained on the testing data is then computed, and the results are reported in Table 7. In the table, the classification accuracy attained by directly fitting the LR model on the original training data and testing on the original testing data is also listed for easy comparison.

Face (%) | Non-face (%) | Total (%) | |
---|---|---|---|

LR | 96.71 | 93.57 | 94.47 |

PCA + LR | 96.64 | 94.17 | 94.88 |

NSPCA + LR | 94.89 | 93.49 | 93.89 |

N-EMPCA + LR | 96.71 | 94.39 | 95.06 |

ReDaC + LR | 96.78 |
94.46 |
95.84 |

From Table 7, it is clear that the proposed ReDaC method attains the best performance among all implemented methods, most accurately recognizing both the face images and the non-face images from the testing data. This further implies the potential usefulness of the proposed method in real applications.

## 4Conclusion

In this paper we have proposed a novel recursive divide-and-conquer method (ReDaC) for sparse PCA problem. The main methodology of the proposed method is to decompose the original large sparse PCA problem into a series of small sub-problems. We have proved that each of these decomposed sub-problems has a closed-form global solution and can thus be easily solved. By recursively solving these small sub-problems, the original sparse PCA problem can always be very effectively resolved. We have also shown that the new method converges to a stationary point of the problem, and can be easily extended to other sparse PCA problems with certain constraints, such as nonnegative sparse PCA problem. The extensive experimental results have validated that our method outperforms current sparse PCA methods in both reconstruction-error-minimization and data-variance-maximization viewpoints.

There are many interesting investigations still worthy to be further explored. For example, when we reformat the square -norm error of the sparse PCA model as the -norm one, the robustness of the model can always be improved for heavy noise or outlier cases, while the model is correspondingly more difficult to solve. By adopting the similar ReDaC methodology, however, the problem can be decomposed into a series of much simpler sub-problems, which are expected to be much more easily solved than the original model. Besides, although we have proved the convergence of the ReDaC method, we do not know how far the result is from the global optimum of the problem. Stochastic global optimization techniques, such as simulated annealing and evolution computation methods, may be combined with the proposed method to further improve its performance. Also, more real applications of the proposed method are under our current research.

## Appendix A. Proof of Theorem 2

In the following, we denote , and the hard thresholding function, whose -th element corresponds to , where is the -th element of and (equals if is ture , and otherwise) is the indicator function

**Theorem 2.** The optimal solution of

is given by:

where denotes the -th largest element of .

## Appendix B. Proof of Theorem 3

We denote the permutation of based on the ascending order of , the soft thresholding function , and throughout the following.

**Theorem 3.** The optimal solution of

is given by:

where for ,

where , .

## Appendix C. Proof of Theorem 5

**Theorem 5.** The global optimal solution to (Equation 15) is (), where , and and are defined in Theorem 2 and Theorem 3, respectively.

It is easy to prove this theorem based on the following lemma.

### Footnotes

- It should be noted that the orthogonality constraints of PC loadings in (Equation 1) and (Equation 2) are not imposed in (Equation 3) and (Equation 4). This is because simultaneously enforcing sparsity and orthogonality is generally a very difficult (and perhaps unnecessary) task. Like most of the existing sparse PCA methods [?], we do not enforce orthogonal PCs in the models.