In this paper, we consider two important problems arising from high-dimensional statistics and machine learning: sparse principal component analysis (PCA) and sparse canonical correlation analysis (CCA). Both problems are special cases of the following more general manifold optimization problem:
where H(A, B) is a smooth function of A, B with a Lipschitz continuous gradient, f (•) and g(•) are lower semicontinuous (possibly nonsmooth) convex functions, and } 1 , } 2 are two embedded submanifolds in the Euclidean space. Problem (1) is numerically challenging because of the combination of nonsmoothness and nonconvexity.
The contributions of this paper lie in several folds. First, we propose a new alternating manifold proximal gradient method for solving (1) and establish its convergence guarantee. Second, we show how to apply the proposed algorithm to solve sparse PCA and sparse CCA. Third, we conduct extensive numerical experiments to compare the performance of the proposed algorithm and existing methods for solving these two problems. Fourth, we analyze the convergence guarantee of an inexact (and more practical) version of the proposed method. PCA was invented by Pearson (1901) and is widely used in dimension reduction. Let X [X 1 , . . . , X p ] ∈ R n×p be a given data matrix whose column means are all zero. Assume that the singular-value decomposition (SVD) of X is X UDV . Then it is known that Z UD are the principal components (PCs), and the columns of V are the corresponding loadings of the PCs. In other words, the first PC can be defined as Z 1 ∑p j 1 α 1j X j with α 1 (α 11 , . . . , α 1p ) maximizing the variance of Z 1 ; that is,
where Σ (X X)/(n -1) is the sample covariance matrix. The remaining PCs are defined as α k+1 arg max α α Σα s.t. α 2 1, α α l 0, ∀1 ≤ l ≤ k.
CCA, introduced by Hotelling (1936), is another widely used tool that explores the relation between two sets of variables. For random variables x ∈ R p and y ∈ R q , CCA seeks linear combinations of x and y such that the resulting values are mostly correlated. That is, it targets to solve the following optimization problem:
where Σ x and Σ y are variances of x and y, Σ xy is their covariance matrix, and u ∈ R p and v ∈ R q are the first canonical vectors. It can be shown that solving (2) corresponds to computing the SVD of Σ -1/2
x Σ xy Σ -1/2 y . In practice, given two centered data sets X ∈ R n×p , Y ∈ R n×q with a joint covariance matrix [Σ x , Σ xy ; Σ yx , Σ y ], CCA seeks the coefficients u, v such that the correlation of Xu and Yv is maximized. The classical CCA (Hotelling 1936) can be formulated as max
where X Y, X X, Y Y are used to estimate the true parameters Σ xy , Σ x , Σ y after scaling. However, PCA and CCA perform poorly and often lead to wrong findings when modeling with highdimensional data. For example, when the dimension is proportional to the sample size such that lim n→∞ p/n γ ∈ (0, 1) and the largest eigenvalue λ 1 ≤ ̅̅ γ √ , the leading sample principal eigenvector could be asymptotically orthogonal to the leading population principal eigenvector under a multivariate Gaussian model (Paul 2007) or a spiked covariance model with a single component and finite fourth moment (Nadler 2008). The orthogonality result is built on the almost sure limits of sample eigenvalues (Baik and Silverstein 2006). Sparse PCA and sparse CCA are proposed as the more interpretable and reliable dimension-reduction and featureextraction techniques for high-dimensional data. In what follows, we provide a brief overview of their methodological developments.
Sparse PCA seeks a sparse basis (loadings) of the subspace spanned by the data so that the leading PCs obtained are easier to interpret. Jolliffe et al. (2003) proposed the SCoTLASS procedure (Simplified Component Technique-LASSO) by imposing an 1 norm on the loading vectors, which can be formulated as the following optimization problem for given data X ∈ R n×p :
where Tr(Z) denotes the trace of matrix Z, μ > 0 is a weighting parameter, A 1 ∑ ij |A ij |, and St(p, r) : {A ∈ R p×r : A A I r } is the Stiefel manifold, where I r is the r × r identity matrix. Note that the original SCoTLASS model in Jolliffe et al. (2003) uses an 1 constraint A 1 ≤ t instead of penalizing A 1 in the objective. The SCoTLASS model ( 4) is numerically challenging. Algorithms for solving it have been limited. As a result, a new formulation of sparse PCA has been proposed by Zou et al. (2006), and it has been the main focus in the literature on this topic. Zou et al. (2006) formulate the sparse PCA problem as the following ridge regression problem plus a LASSO penalty:
where x i denotes the transpose of the ith row vector of X, B j is the jth column vector of B, and μ > 0 and μ 1,j > 0 are weighting parameters. Given A ∈ St(p, r), we have the following:
Therefore, (5) can be equivalently written as follows:
min
where
The sparse PCA ( 6) is in the form of (1) with r), and } 2 R p×r . However, it should be noted that ( 6) is indeed still numerically challenging. Zou et al. (2006) proposed to solve it using an alternating minimization algorithm (AMA) that updates A and B alternatingly with the other variable fixed as the current iterate. A typical iteration of AMA is as follows:
The A-subproblem in ( 7) is known as a Procrustes rotation problem and has a closed-form solution given by a SVD. The B-subproblem in ( 7) is a linear regression problem with an elastic net regularizer, and it can be solved by many existing solvers such as elastic net (Zou and Hastie 2005), coordinate descent (Friedman et al. 2010), and the fast iterative shrinkage-thresholding algorithm (FISTA; Beck and Teboulle 2009). However, there is no convergence guarantee of AMA ( 7). Recently, some new algorithms are proposed in the literature that can solve (6) with guarantees of convergence to a stationary point. We provide a summary of some representative algorithms in the next section. We need to point out that there are other ways to formulate sparse PCA such as the ones in d 'Aspremont et al. (2007'Aspremont et al. ( , 2008)), Ma (2013), Lu and Zhang (2012), Vu et al. (2013), d'Aspremont (2011), Shen and Huang (2008), Witten et al. (2009), Journee et al. (2010), Yuan andZhang (2013), andMoghaddam et al. (2005). We refer interested readers to the recent survey paper by Zou and Xue (2018) for more details on these works on sparse PCA. In this paper, we focus on the formulation of (6) to estimate multiple PCs, which is a manifold optimization problem with nonsmooth objective function.
Sparse CCA (Wiesel et al. 2008, Parkhomenko et al. 2009, Witten et al. 2009, Hardoon and Shawe-Taylor 2011) is proposed to improve the interpretability of CCA, which can be formulated as follows:
where X ∈ R n×p , Y ∈ R n×q , f , and g are regularization terms promoting the sparsity of u and v, and common choices for them include the 1 norm for sparsity and the 2,1 norm for group sparsity. When multiple canonical vectors are needed, one can consider the matrix counterpart of (8) that can be formulated as follows:
min
where r is the number of canonical vectors needed. From now on, we call (8) the single sparse CCA model and ( 9) the multiple sparse CCA model. Moreover, motivated by Gao et al. (2017), in this paper we choose f and g to be the 2,1 norm to promote the group sparsity of A and B in (9). Specifically, we choose f (A) τ 1 A 2,1 and g(B) τ 2 B 2,1 , where the 2,1 norm is defined as A 2,1 ∑ p j 1 A j• 2 , and A j• denotes the jth row vector of matrix A, and τ 1 > 0 and τ 2 > 0 are weighting parameters. In this case, the multiple sparse CCA (9) reduces to min A∈R p×r ,B∈R q×r -Tr A X YB
which is in the form of (1) with
Note that in ( 8) and ( 9) we assumed that both X X and Y Y are positive definite to guarantee that } 1 and } 2 are submanifolds. If they are not positive definite, we can always add a small perturbation to make them so. These manifolds used in ( 8) and ( 9) are generalized Stiefel manifolds. Note that when r 1, the 2,1 norm becomes the 1 norm of the vector because u 2,1 u 1 for any u ∈ R p . In this case, (8) reduces to min
Manifold optimization has recently drawn a lot of research attention because of its success in a variety of important applications, including low-rank matrix completion (Boumal andAbsil 2011, Vandereycken 2013), phase retrieval (Bendory et al. 2018, Sun et al. 2018), phase synchronization (Boumal 2016, Liu et al. 2017), blind deconvolution (Huang and Hand 2018), and dictionary learning (Cherian andSra 2017, Sun et al. 2017).
Most existing algorithms for solving manifold optimization problems rely on the smoothness of the objective (see the recent monograph by Absil et al. (2009)). Studies on manifold optimization problems with nonsmooth objectives such as ( 6), (10), and (11) have been limited. This urges us to study efficient algorithms that solve manifold optimization problems with nonsmooth objectives, and this is the main focus of this paper. The rest of this paper is organized as follows. We review existing methods for sparse PCA and sparse CCA in Section 2. We propose a unified alternating manifold proximal gradient method with provable convergence guarantees for solving both sparse PCA and sparse CCA in Section 3. The numerical performance is demonstrated in Section 4. We provide preliminaries on manifold optimization and details of the global convergence analysis of our proposed method in the online appendix. We also discuss the convergence result of an inexact version of our algorithm in the online appendix.
Before proceeding, we review existing methods for solving sparse PCA (6) and sparse CCA (8) and (9).
2.1. Existing Methods for Sparse PCA For sparse PCA (6), other than the AMA algorithm suggested in Zou et al. (2006), there exist some other efficient algorithms for solving this problem. We now give a brief review of these works. We first introduce two powerful optimization algorithms for solving nonconvex problems: the proximal alternating minimization (PAM) algorithm (Attouch et al. 2010) and the proximal alternating linearization method (PALM) (Bolte et al. 2014). Surprisingly, it seems that these two methods have not been used to solve (6) yet. We now briefly describe how these two methods can be used to solve (6). PAM for (6) solves the following two subproblems in each iteration:
where t 1 > 0, t 2 > 0 are step sizes. Note that in each subproblem, PAM minimizes the objective function with respect to one variable by fixing the other, and a proximal term is added for the purpose of convergence guarantee. It is shown in Attouch et al. (2010) that the sequence of PAM converges to a critical point of (6) under the assumption that the objective function satisfies the Kurdyka-Łojasiewicz (KL) inequality. 1 We need to point out that the only difference between PAM ((12) and ( 13)) and the AMA (7) is the proximal terms, which together with the KL inequality helps establish the convergence result. Note that the A-subproblem in (12) corresponds to the reduced-rank Procrustes rotation and can be solved by a SVD. The B-subproblem in ( 13) is a LASSO-type problem and can be solved efficiently by first-order methods such as FISTA or block coordinate descent. A better algorithm that avoids an iterative solver for the subproblem is PALM, which linearizes the quadratic functions in the subproblems of ( 12) and (13). A typical iteration of PALM is the following:
where ∇ A H and ∇ B H denote the gradient of H with respect to A and B, respectively. The two subproblems in ( 14) and ( 15) are easier to solve than the ones in ( 12) and ( 13) because they both admit closed-form solutions. In particular, the solution of the A-subproblem in ( 14) corresponds to the projection onto the orthogonality constraint, which is given by a SVD; the solution of the B-subproblem in ( 15) is given by the 1 softthresholding operation. It is shown in Bolte et al. (2014) that the sequence of PALM converges to a critical point of ( 6) under the assumption that the objective function satisfies the KL inequality. Recently, Erichson et al. (2020) proposed a projected gradient method based on variable projection (VP) for solving (6). Though the motivation of this algorithm is different, it can be viewed as a variant of PAM and PALM. Roughly speaking, the VP algorithm combines the A-subproblem (without the proximal term) in ( 12) and the B-subproblem in (15). That is, it updates the iterates as follows:
Note that the difference of PALM (( 14) and ( 15)) and VP (( 16) and ( 17)) lies in the A-subproblem. The A-subproblem linearizes the quadratic function H(A, B k ) in ( 14) but not in ( 16). This does not affect the performance of the algorithms much because in this specific problem, the A-subproblems correspond to a SVD in both algorithms. It is shown in Erichson et al. (2020) that the VP (( 16) and ( 17)) converges to a stationary point of ( 6). Another recent work that can solve ( 6) is the manifold proximal gradient (ManPG) algorithm proposed by Chen et al. (2020). We will discuss it in more details later because it is closely related to the algorithm we propose in this paper. For other algorithms for solving sparse PCA, we refer interested readers to a recent survey paper (Zou and Xue 2018) for more details.
Chen et al. ( 2013) studied a canonical correlation analysis via the precision-adjusted iterative thresholding (CAPIT) algorithm for solving the single sparse CCA (11). CAPIT alternates between an iterative thresholding step and a power method step to deal with the sparsity regularization and orthogonality constraints, respectively. The convex program with group LASSO refinement (CoLaR) method proposed by Gao et al. (2017) targets the multiple sparse CCA (9). CoLaR is a two-stage algorithm. In the first stage, a convex relaxation of (9) based on the matrix lifting technique is solved. In the second stage, the solution obtained from the first stage is refined by solving a group LASSO-type problem. Wiesel et al. (2008) suggested a greedy approach for solving (2) with cardinality constraints on u and v. Recently, Suo et al. (2017) developed an AMA for solving the single sparse CCA (11) that solves two subproblems in each iteration by solving (11) with respect to u (respectively, v) with v (respectively, u) fixed as v k (respectively, u k ). The subproblems were then solved by a linearized alternating direction method of multipliers (ADMM) algorithm. We need to point out that none of these algorithms for sparse CCA has a convergence guarantee. There exist some other methods for sparse CCA (see, e.g., Witten et al. 2009, Hardoon andShawe-Taylor 2011), but we omit their details here because they are not directly related to (8) and ( 9). We also point out that the PAM, PALM, and VP algorithms discussed in Section 2.1 are not good choices for sparse CCA (( 8) and ( 9)) because the resulting subproblems are not easy to solve and require expensive subroutines. For instance, to apply PALM to (8), one needs to compute the proximal mapping of f (u) + ι(u X Xu 1), which does not admit a closed-form solution and is thus computationally expensive, where ι(•) denotes the indicator function.
In this section, we introduce our alternating manifold proximal gradient (A-ManPG) algorithm for solving (1).
The following assumptions are made for problem (1).
Assumption 1. We assume that the following are true: a. The functions f and g are convex and Lipschitz continuous with Lipschitz constants L f and L g , respectively, and their proximal mappings are semismooth. 2 b. The function ∇ A H(A, B) is Lipschitz continuous with respect to A when fixing B, and the Lipschitz constant is L A (B). Similarly, ∇ B H(A, B) is Lipschitz continuous with respect to B when fixing A, and the Lipschitz constant is L B (A). Moreover, there exist constants L A > 0 and L
c. The function F is lower bounded by a constant F * . d. The function F is coercive with respect to (A, B). An immediate consequence of this assumption is that the sublevel set
Remark 1. In Assumption 1(a), the Lipschitz continuity and convexity are all defined in the Euclidean space. In Assumption 1(b), it is easy to obtain that L A (B) ≡ 0 and L B (A) 2λ max (X X) for problem (6) and L A (B) L B (A) ≡ 0 for problems (8) and(10). Assumptions 1(c) and1(d) can be removed if } 1 and} 2 are compact submanifolds. The ManPG algorithm proposed by Chen et al. (2020) can be applied to solve (1). In each iteration, ManPG linearizes H(A, B) and solves the following convex subproblem:
where t < 1/L, and L is the Lipschitz constant of ∇H(A, B) on the tangent space
Here A, B Tr(A B) denotes the inner product. Note that ( 18) is actually separable for D A and D B and thus reduces to two subproblems for D A and D B , respectively. As a result, ManPG (18) can be viewed as a Jacobi-type algorithm in this case because it computes D A and D B in parallel. We found from our numerical experiments that the algorithm converges faster if D A and D B are computed in a Gauss-Seidel manner. This leads to the following updating scheme, which is the basis of our A-ManPG algorithm:
where A k+1 is obtained via a retraction operation (see Algorithm 1), t 1 ≤ 1/L A , and t 2 ≤ 1/L B . The Gauss-Seideltype algorithm A-ManPG usually performs much better than the Jacobi-type algorithm ManPG, because the Lipschitz constants are smaller, and thus, larger step sizes are allowed. The details of the A-ManPG algorithm are described in Algorithm 1.
Algorithm 1 (A-ManPG Method)
). 15: end for Remark 2. Note that the iterates A k and B k are kept on the manifolds through the retraction operations Retr A and Retr B . There exist many choices for the retraction operations, and in Algorithm 1, we did not specify which ones to use. We discuss common retractions for the Stiefel manifold and the generalized Stiefel manifold in the online appendix. In our numerical experiments in Section 4, we chose polar decomposition as the retraction. Lines 4-7 and 10-13 in Algorithm 1 are backtracking line search procedures. These are necessary to guarantee that the objective function has a sufficient decrease in each iteration, which is needed for the convergence analysis (see the online appendix).
From Lemma 3 (see the online appendix), we know that D A k 0 and D B k 0 imply that (A k , B k ) is a stationary point (see the online appendix, Definition 5) for problem (1). As a result, we can define an -stationary point of (1) as follows.
Definition 1. We call (A, B) an -stationary point of (1) if D A and D B returned by ( 19) and (20) with step sizes
We have the following convergence results for the A-ManPG algorithm (Algorithm 1), whose proof is given in the online appendix.
Theorem 1. Under Assumption 1, any limit point of the sequence {(A k , B k )} generated by Algorithm 1 is a stationary point of problem (1). Moreover, Algorithm 1 with step sizes t 1 1/L A , t 2 1/L B returns an -stationary point (A k , B k ) in at most 2(F(A 0 , B 0 ) -F * )/((γ ᾱ1 t 1 + γ ᾱ2 t 2 ) 2 ) iterations, where ᾱ1 > 0 and ᾱ2 > 0 are constants.
The main computational effort in each iteration of Algorithm 1 is to solve ( 19) and (20). For the Stiefel manifold and the generalized Stiefel manifold, ( 19) and ( 20) are both equality-constrained convex problems, given that both f and g are convex functions. Note that if f (respectively, g) vanishes, the A-subproblem (respectively, B-subproblem) becomes the projection onto the tangent space of } 1 (respectively, } 2 ), which reduces to a Riemannian gradient step and can be easily done. Here we discuss the general case where f and g do not vanish. The subproblem solver should be efficient and yields highly accurate solutions. We have tested several possible candidates of subproblem solvers including ADMM and the augmented Lagrangian method but found that they are either less efficient or do not yield a highly accurate solution. In the end, we found that an adaptive semismooth Newton (SSN) method (Xiao et al. 2018) is suitable for solving this kind of problem. The notion of semismoothness was originally introduced by Mifflin (1977) for real-valued functions and extended to vector-valued mappings by Qi and Sun (1993). In a pioneering work on the SSN method, Solodov and Svaiter (1998) proposed a globally convergent Newton's method by exploiting the structure of monotonicity and established a local superlinear rate under the conditions that the operator is semismooth and the generalized Jacobian is nonsingular at the global optimal point. The convergence rate is extended in Zhou and Toh (2005) to the setting where the generalized Jacobian is not necessarily nonsingular. Recently, the SSN method has received significant attention because of its success in solving structured convex problems to high accuracy (see Wang et al. 2010;Zhao et al. 2010;Qi and Sun 2011;Yang et al. 2013Yang et al. , 2015;;Li et al. 2018;Xiao et al. 2018).
We now describe how to apply the adaptive SSN (ASSN) method in Xiao et al. (2018) to solve ( 19) and ( 20). The derivation here closely follows the one in Chen et al. (2020). For brevity, we only focus on the A-subproblem with } 1 {A | A X XA I r } and f (A) τ 1 A 2,1 as used in (10). Note that for the Stiefel manifold } St(p, r), its tangent space is given by T A } {D | D A + A D 0}, and for the generalized Stiefel manifold } {A | A MA I r }, its tangent space is given by T A } {D | D MA + A MD 0}. For ease of notation, we denote t t 1 , D D A , M : X X, and h(A) : H(A, B k ). In this case, the A-subproblem in (19) reduces to
By associating a Lagrange multiplier Λ with the linear equality constraint, the Lagrangian function of ( 21) can be written as follows:
and the Karush-Kuhn-Tucker system of ( 21) is given by
The first condition in (23) implies that D can be computed by
where prox f (A) denotes the proximal mapping of function f at point A. By substituting (24) into the second condition in ( 22), we obtain that Λ satisfies
and thus the problem reduces to finding a root of function E. Because E is a monotone operator (see Chen et al. 2020) and the proximal mapping of the 2 norm is semismooth, 3 we can apply SSN to find the zero of E. The SSN method requires computation of the generalized Jacobian of E, and in the following, we show how to compute it. We first derive the vectorization of E(Λ):
,
where K rp and K rr denote the commutation matrices. We define the following matrix:
where ⊗ denotes the Kronecker product, and )(y) is the generalized Jacobian of prox tf (y). The matrix )(y) is defined as
, where the matrices Δ j , j 1, . . . , p are defined as follows: (Xiao et al. 2018). For the local convergence rate, one needs more conditions on ∂vec(E(vec(Λ))). We refer to Xiao et al. (2018) for more details. Note that because Λ is a symmetric matrix, we can work with the lower triangular part of Λ only and remove the duplicated entries in the upper triangular part. To do so, we use vec(Λ) to denote the 1 2 r(r + 1)-dimensional vector obtained from vec(Λ) by eliminating all superdiagonal elements of Λ. It is known that there exists a unique r 2 × 1 2 r(r + 1) matrix U r , which is called the duplication matrix (Magnus and Neudecker 2018, chapter 3.8), such that U r vec(Λ) vec(Λ). The Moore-Penrose inverse of U r is U + r (U r U r ) -1 U r , and it satisfies U + r vec(Λ) vec(Λ). The alternative of the generalized Jacobian of vec(E(U r vec(Λ))) is given by
where we use the identity K rr + I r 2 2U r U + r . Therefore, (26) can be simplified to G vec Λ ( ) ( )
Because both U r and U + r have only r 2 nonzero elements, the computation flop of matrix G(vec(Λ)) is 2pr 4 + 2r(r + 1)r 2 . The ASSN method in Xiao et al. (2018) first computes Newton's direction d k by solving
where η > 0 is a regularization parameter. Note that η is necessary here because G(vec(Λ)) could be singular if Δ j 0 for some j. We then update
Thus, the main computational cost of the SSN algorithm lies in solving (28), which is 2(r 6 ) in each iteration. Therefore, performing a SSN step is faster than computing the gradient of the objective function in (1), which is 2(p 2 r), if r is much smaller than p. We describe the detailed ASSN algorithm proposed in Xiao et al. (2018) in the online appendix.
4.1. Sparse PCA In this section, we apply our A-ManPG algorithm 5 to solve sparse PCA ( 6) and compare its performance with that of three existing methods: AMA (Zou et al. 2006), PALM (Bolte et al. 2014), and VP (Erichson et al. 2020).
The details of the parameter settings of these algorithms are as follows:
• AMA (7): Maximum iteration number of AMA is set to 1,000. We use FISTA (Beck and Teboulle 2009) to solve the B-subproblem, and we set the maximum iterations of FISTA to 1,000 and stop it if A j -A j-1 / A j-1 ≤ 10 -6 , where {A j } is the sequence in FISTA.
• PALM (( 14) and ( 15)): Here t 1 : 1, t 2 : 1/(2λ max (X X)). Maximum iteration number is set to 10,000.
• VP (( 16) and ( 17)): Here t 2 : 1/(2λ max (X X)). iteration number is set to 10,000.
• A-ManPG: Here t 1 100/p, t 2 : 1/(2λ max (X X)). Maximum iteration number is set to 10,000. The algorithms are terminated using the following criteria. First, we use the PALM algorithm as a baseline, and we denote the objective function value in (6) as F(A, B), that is., F(A, B) H(A, B) + μ ∑ r j 1 B j 2 + ∑ r j 1 μ 1,j B j 1 . We terminate PALM when we find that
We then terminate the AMA, A-ManPG, and VP algorithms when their objective function values are smaller than F PALM and the change of their objective values in two consecutive iterations is less than 10 -5 . We generate the data matrix X in the following manner. First, the entries of X are generated following the standard normal distribution 1(0, 1). The columns of X are then centered so that the columns have zero mean, and they are then scaled by dividing the largest 2 norm of the columns. We report the comparison results of the four algorithms in Tables 1 and2, where r 6 for all cases. In particular, Table 1 reports the results for n < p, and we tested μ 1 and μ 10 because it is suggested in Zou et al. (2006) that μ should be relatively large in this case. Table 2 reports the results for n > p, and we set μ 10 -6 because it is suggested in Zou et al. (2006) that μ should be sufficiently small in this case. In these tables, central processing unit (CPU) times are in seconds, and "sp" denotes the percentage of zero entries of matrix B. From Tables 1 and2, we see that the four algorithms generated solutions with similar objective function value F(A, B) and similar sparsity sp. In terms of CPU time, AMA is the slowest, and the other three are comparable and are all much faster than AMA. This is because AMA needs an iterative solver to solve the B-subproblem, which is time-consuming in practice. The comparison results are reported in Figures 1 and2. These results indicate that A-ManPG is at least two times faster than ManPG-Ada on the tested synthetic problems. More comparison results are presented in Figures 3 and4 in the online appendix.
In this section, we report the numerical results of the A-ManPG algorithm for solving the single sparse CCA (11) and compare its performance with a recent approach proposed by Suo et al. (2017): AMA + linearized alternating direction method of multipliers (LADMM). More specifically, the AMA + LADMM method aims to solve the relaxation of (11) as follows:
We generate the data in the same manner as Suo et al. (2017) and v are the true canonical vectors, and ρ is the true canonical correlation. In our numerical tests, û and v are generated randomly such that they both have five nonzero entries and the nonzero coordinates are set at the {1, 6, 11, 16, 21}th coordinates. The nonzero entries are obtained from normalizing (with respect to Σ x and Σ y ) random numbers drawn from the uniform distribution on the finite set {-2, -1, 0, 1, 2}. We set ρ 0.9 in all tests. We tested three different ways to generate the covariance matrices Σ x and Σ y that are exactly the same as presented in section 4.1 of Suo et al. (2017): identity matrices, Toeplitz matrices, and sparse inverse matrices.
For succinctness, we omit them here but include them in the online appendix. The matrices X and Y are both divided by ̅̅̅̅̅̅̅ n -1 √ such that X Y is the estimated covariance matrix. Note that if n < p or n < q, the covariance matrix X X or Y Y is not positive definite. In this case, we replace X X by (1 -α)X X + αI p and Y Y by (1 -α)Y Y + αI q in the constraints of ( 11), so we can still keep them as manifold constraints. In our experiments, we chose α 10 -4 . We adopted the same initialization procedure as suggested in Suo et al. (2017). We set 11), where b was set to b {1, 1.2, 1.4, 1.6}. We report the best result among all the candidates. For each b, we solved (11) by A-ManPG with δ 10 -4 , γ 0.5, and t 1 t 2 1. The A-ManPG algorithm was stopped if max{
, and the ASSN algorithm was stopped if E(Λ k ) F ≤ 10 -5 in (26). For AMA + LADMM, we set the stopping criteria of LADMM as u ju j-1 ≤ 10 -3 and v jv j-1 ≤ 10 -3 , where u j and v j are iterates in LADMM. We set the stopping criteria of AMA as u ku k-1 ≤ 10 -3 and v kv k-1 ≤ 10 -3 , where u k and v k are iterates in AMA. We report the numerical results in Table 3, where lossu and lossv are defined the same way as in Suo et al. (2017), nu and nv denote the number of nonzeros in u and v after setting their entries whose magnitudes are smaller than 10 -4 to zero, and ρ denotes the canonical correlation computed from the solution returned by the algorithms. All reported values in Table 3 are the medians from 20 repetitions. From Table 3, we see that A-ManPG and AMA + LADMM achieve similar loss function values lossu and lossv, but A-ManPG is usually faster than AMA+LADMM, and for some cases, it is even two to three times faster. More important, AMA + LADMM lacks convergence analysis, but A-ManPG is guaranteed to converge to a stationary point (see the Note. CPU times are given in seconds; lossu and lossv are defined the same way as in Suo et al. (2017); nu and nv denote the number of nonzeros in u and v after setting their entries whose magnitudes are smaller than 10-4 to zero; and ρ denotes the canonical correlation computed from the solution returned by the algorithms. The best ones among the four algorithms are highlighted using bold text. online appendix). Furthermore, AMA + LADMM is time-consuming for the multiple sparse CCA (10), but A-ManPG is suitable for (10), as we show in the next section.
In this section, we apply the A-ManPG algorithm to solve the multiple sparse CCA (10) and compare its performance with the Convex program with group-LASSO Refinement (CoLaR) method proposed by Gao et al. (2017). CoLaR is a two-stage method based on convex relaxations. In all tests, we chose r 2 and generated Σ xy Σ x UΛV Σ y , where Λ ∈ R r×r is a diagonal matrix with diagonal entries Λ 11 0.9 and Λ 22 0.8. The nonzero rows of both U and V are set at the {1, 6, 11, 16, 21}th rows. The values at the nonzero coordinates are obtained from normalizing (with respect to Σ x and Σ y ) random numbers drawn from the uniform distribution on the finite set {-2, -1, 0, 1, 2}. The two data sets X ∈ R n×p and Y ∈ R n×q are generated in the same way as in Section 4.2. The matrices X and Y are both divided by ̅̅̅̅̅̅̅ n -1 √ such that X Y is the estimated covariance matrix. The loss between the estimation A and the ground truth U is measured by the subspace distance lossu P U -P A 2 F , where P U denotes the projection matrix onto the column space of U. Similarly, the loss for B and V is measured as lossv P V -P B
F . The codes of CoLaR were downloaded from the authors' web page. 6 We used all default settings of their codes. In particular, ADMM is used to solve the first-stage problem ((E.3) in the online appendix), and it is terminated when it does not make much progress or it reaches the maximum iteration number 100. For our A-ManPG algorithm, we run only one iteration of ADMM for (E.3) and use the returned solution as the initial point of A-ManPG, because we found that this already generates a very good solution for A-ManPG. To be fair, we also compare the same case for CoLaR, where only one iteration of ADMM is used for (E.3). The parameter τ in (E.3) is set to τ 0.55
10) and τ b in (E.4) (in the online appendix), where b was set to b {0.8, 1, 1.2, 1.4, 1.6}. For each b, we solved (10) by A-ManPG with δ 10 -4 , γ 0.5, t 1 t 2 1. The A-ManPG algorithm was stopped if max{ D
and the ASSN algorithm was stopped if E(Λ k ) F ≤ 10 -5 in (25).
We report the numerical results in Tables 4567, where CPU times are in seconds, and nA and nB denote the number of nonzeros of A and B, respectively, after truncating the entries whose magnitudes are smaller than 10 -4 to zeros. The scalars ρ 1 and ρ 2 are the two canonical correlations, and they should be close to 0.9 and 0.8, respectively. All reported values in Tables 4-7 are the medians from 20 repetitions. More specifically, Table 4 reports the results obtained from the first stage where ADMM was used to solve (E.3). We use Init-1 to indicate that we only run one iteration of ADMM, and Init-100 to indicate the case where we run ADMM until it does not make much progress or the maximum iteration number 100 is reached. From Table 4, we see that solving the first-stage problem by running ADMM for 100 iterations indeed improves the two losses significantly. Tables 5-7 report the results for the three different types of covariance matrices. A-ManPG-1 and CoLaR-1 are the versions where we run only one iteration of ADMM for the first stage, and CoLaR-100 is the version where the first-stage problem ((E.3) in the online appendix) is solved more accurately by ADMM, as discussed earlier.
We observed that running more iterations of ADMM in the first stage does not help much for A-ManPG; we thus only report the results of A-ManPG-1. From Tables 5-7, we see that CoLaR-100 gives much better results than CoLaR-1 in terms of the two losses lossu and lossv, especially when the sample size is relatively small compared with the matrix sizes. Moreover, we see that A-ManPG-1 outperforms both CoLaR-1 and CoLaR-100 significantly. In particular, A-ManPG-1 generates comparable and often better solutions than CoLaR-1 and CoLaR-100 in terms of solution sparsity and losses lossu and lossv. Furthermore, A-ManPG-1 is usually faster than CoLaR-1 and much faster than CoLaR-100.
In this paper, we proposed an efficient algorithm for solving two important and numerically challenging optimization problems arising from statistics: sparse PCA and sparse CCA. These two problems are challenging to solve because they are manifold optimization problems with nonsmooth objectives, a topic that is still underdeveloped in optimization. We proposed an A-ManPG method to solve these two problems.
Convergence and convergence rate to a stationary point of the proposed algorithm are established. Numerical results on statistical data demonstrate that A-ManPG is comparable to existing algorithms for solving sparse PCA and is significantly better than existing algorithms for solving sparse CCA.
INFORMS Journal on Optimization, Articles inAdvance, pp. 1-17, © 2020 INFORMS