A Note on Sparse Generalized Eigenvalue Problem

Abstract:
The sparse generalized eigenvalue problem (SGEP) aims to find the leading eigenvector with sparsity structure. SGEP plays an important role in statistical learning and has wide applications including, but not limited to, sparse principal component analysis, sparse canonical correlation analysis and sparse Fisher discriminant analysis, etc. Due to the sparsity constraint, the solution of SGEP entails interesting properties from both numerical and statistical perspectives. In this paper, we provide a detailed sensitivity analysis for SGEP and establish the rate-optimal perturbation bound under the sparse setting. Specifically, we show that the bound is related to the perturbation/noise level and the recovery of the true support of the leading eigenvector as well. We also investigate the estimator of SGEP via imposing a non-convex regularization. Such estimator can achieve the optimal error rate and can recover the sparsity structure as well. Extensive numerical experiments corroborate our theoretical findings via using alternating direction method of multipliers (ADMM)-based computational method.

heading: 1 Introduction
text: The sparse generalized eigenvalue problem (SGEP) is to solve the following constrained optimization problem:
max x∈Rp
xTÃx xTB̃x , subject to ‖x‖0 ≤ k, (1)
where Ã, B̃ are both p-by-p symmetric matrices and B̃ is (semi) positive definite; ‖x‖0 denotes the `0-norm of x, i.e., the number of nonzero entries of x; k is an integer which may be much smaller than p. For a fixed k, we let x∗ be the solution to (1). Moreover, matrices Ã and B̃ usually have the following decomposition:
Ã = A+ E, B̃ = B + F, (2)
where A, B are underlying symmetric matrices and B is positive definite. However, A and B are unobserved in practice and they can be contaminated by some noise/perturbation matrices E and F , respectively. Let λ1 be the largest generalized eigenvalue of matrix pair (A,B) and u1 be its corresponding eigenvector which is also known as the leading eigenvector. In addition, u1 is always assumed to be sparse, i.e., ‖u1‖0 p. Therefore, the SGEP is essentially to find an approximation of u1 based on observation (Ã, B̃).
The SGEP is challenging due to the following aspects. First, the sparsity constraint optimization problem (1) is NP-hard, since it is essentially a subset selection problem [19, 20]. To recover the underlying sparsity structure, we need to try different k which is computationally expensive.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Second, due to insufficient number of samples, Ã and B̃ can be both ill-conditioned or even share an approximate common null space, and perturbations E and F can be large. As a result, the leading eigenvector of pair (Ã, B̃) may not be a good approximation of u1.
Though challenging, especially when there is a limited number of data samples, SGEP and related problems have received more and more attentions in the past decade and various numerical methods are proposed to solve the SGEP. For example, in [24], the SGEP is framed as a d.c. (difference of convex functions) program and is solved as a sequence of convex programs by invoking the majorization-minimization method. In [23], the `0 norm constraint is approximated by a continuous surrogate function, an algorithm inspired by the majorization-minimization method is developed via iteratively majorizing the surrogate function by a quadratic separable function, and a systematic way based on smoothing is proposed to deal with the singularity issue. [9] proposes a semidefinite programming method for sparse canonical correlation analysis which is one of important SGEP applications. In [30], a two-stage computational framework is proposed to solve the SGEP, in which the first stage computes an initial guess via convex relaxation, and the second stage solves a nonconvex optimization via a fixed step size steepest ascent method, and followed by the simple truncation. In [22], a general framework called sparse estimation with linear programming is proposed for sparse canonical correlation analysis. Combining the idea of `1-penalized adaptive normalized quasi-Newton algorithm with nested orthogonal complement structure, an algorithm is proposed in [32] to solve the SGEP. [13] solves the SGEP by modifying the standard generalized orthogonal iteration with a sparsity-inducing penalty for the eigenvectors. For more comprehensive review of computations in the SGEP, please see [3, 42] and references therein.
Despite the efforts mentioned above, there is very little work investigating the optimal estimation error and no work provides theoretical guarantees whether the underlying sparsity structure of the leading eigenvector of SGEP can be recovered. In this paper, we bridge the theoretical gap in the literature. First, we provide both upper and lower bounds of the approximation error given noisy observations in the sparse generalized eigenvalue problem setting. To be specific, we prove that the difference between u1 and x∗ is bounded by two parts, the perturbation error and error of not recovering the true support of u1. Furthermore, we show that the error bound we obtained is rate optimal. That is, the lower bound matches the upper bound up to a multiplicative constant. Second, we consider a family of non-convex penalty functions and reformulate (1) into a regularization problem. We show that such regularized estimator enjoys the merits of sparsity and nearly unbiasedness. With these nice properties, the estimator can be shown to achieve the optimal estimation error rate and can recover the true support of leading eigenvector as well. In addition, we also present a new computational algorithm. The estimation procedure adopts the Alternating Direction Method of Multipliers (ADMM, [37, 21]) and can recover the sparsity structure well under mild conditions. Local convergence theory is established for the proposed method.
The rest of paper is organized as follows. In Section 2, we introduce the sparse generalized eigenvalue problem with its applications and technical tools. In Section 3, we provide a sensitivity analysis for SGEP and establish its upper and lower perturbation bounds. In Section 4, we propose a regularized estimator via using non-convex penalization method and established the estimation bound. We also propose a computational method, non-convex SGEP algorithm (NC-SGEP), to tackle with the estimation issues. Multiple numerical results are given in Section 5 and corroborate our theory. The concluding remark is given in Section 6.
Notations: The calligraphic letters I , J , K, and S are usually used to denote index sets. |I| denotes the cardinality of I, e.g., I = {i1, i2, . . . , is}, where i1, i2, . . . , is are distinct integers, then |I| = s. We use x to denote a generic vector in Rp and xT to denote its transpose. x[j] is the jth element in x. x[I] stands for the subvector of x with indices in set I. We use A to denote a matrix in Rp×p. A[i, j] is the entry in the ith row and the jth column. A[I,J ] stands for the submatrix of A with row indices in set I and column indices in set J . If I = J , A[I, I] is also denoted by AI . We also let (I,J ) be the union of sets I and J . For x ∈ Rp, supp(x) denotes the index set of all nonzero entries of x. For matrix A, supp(A) denotes the index set of all nonzero rows of A. Ip is the p× p identity matrix. diag(d1, . . . , dp) represents a diagonal matrix with d1, . . . , dp being its diagonal
elements. ‖x‖0 represents the number of non-zero entries in x; ‖x‖2 := √∑ j x 2[j]; ‖A‖2 denotes
the largest singular value of A; ‖A‖F and ‖A‖∗ is the Frobenius norm and nuclear norm of matrix A respectively. ‖A‖∞ := maxi,j |A[i, j]|. We write a b (a b) if a ≥ Kb (a ≤ b/K) for some sufficiently large number K.

heading: 2 Preliminary
text: Applications Sparse generalized eigenvalue problem has a wide application. Many high-dimensional multivariate statistical problems can be formulated as special instances of (1).
Example 1 [Sparse Principal Component Analysis [44, 5, 1, 35, 36]] Given n observations with p features, sparse principal component analysis (SPCA) seeks the best low dimensional projection of the observed data for increasing interpretability, minimizing information loss and achieving sparsity structure. Let Σ̃ be the empirical covariance matrix. SPCA aims to solve
max x∈Rp
xTΣ̃x, subject to ‖x‖0 ≤ k, ‖x‖2 = 1. (3)
Such problem is a special case of SGEP where Ã = Σ̃ and B̃ = I .
Example 2 [Sparse Fisher’s Discriminant Analysis [31, 12, 17, 14, 16]] Given n observations with p features from K different classes, Fisher’s discriminant analysis seeks a projection for mapping observations to the low dimensional space on which the between-class variance Σb is large while the within-class variance Σw is small. Let Σ̃b and Σ̃w be the sample estimates of Σb and Σw, respectively. To obtain the sparse leading discriminant vector, one can solve
max x∈Rp
xTΣ̃bx, subject to xTΣ̃wx = 1, ‖x‖0 ≤ k. (4)
Such a problem can be formulated as an SGEP with Ã = Σ̃b and B̃ = Σ̃w.
Example 3 [Sparse Canonical Correlation Analysis [39, 9, 10, 41]] Given two random vectors X ∈ Rp, Y ∈ Rp, let Σxx, Σyy , Σxy be the covariance matrices for X , Y , and the cross-covariance matrix between X and Y , respectively. Let Σ̃xx, Σ̃yy , Σ̃xy be estimators (constructed from samples) for Σxx, Σyy, Σxy, respectively. Sparse canonical correlation analysis (SCCA) aims to solve the constrained optimization problem:
max ux,uy
uTx Σ̃xyuy, subject to u T x Σ̃xxux = 1, u T y Σ̃yyuy = 1,
‖ux‖0 ≤ kx, ‖uy‖0 ≤ ky, where kx and ky are two small integers. Such a problem can be recast as an SGEP with
Ã =
[ 0 Σ̃xy
Σ̃Txy 0
] , B̃ = [ Σ̃xx 0
0 Σ̃yy
] , u1 = [ ux uy ] . 1
Strictly speaking, two problems are not exactly equivalent, because ‖u‖0 ≤ kx + ky does not necessarily imply ‖ux‖0 ≤ kx and ‖uy‖0 ≤ ky .
Technical preparation We first introduce some useful terminologies for describing the generalized eigenvalue problem (GEP). Let A ∈ Rp×p be a symmetric matrix and B ∈ Rp×p be a symmetric and positive definite matrix. Then GEP for the matrix pair (A,B) is defined as
Aui = λiBui; i = 1, . . . , p.
Without loss of generality, we can always assume that λ1 ≥ · · · ≥ λp. Here ui is ith eigenvector and λi is ith eigenvalue; (λi,ui) is called ith eigenpair of (A,B). Notice that ui is only determined up to a scale. In this paper, we always assume that ‖ui‖2 = 1. The perturbation analysis of generalized eigenvalue problem has been extensively studied since 1970s and many perturbation bounds have been developed [25, 27, 28, 33]. For reader convenience, we present several existing fundamental results before moving to our main theory immediately.
Definition 1 The angle between x,y( 6= 0) ∈ Rp is defined as θ(x,y) := arccos |x Ty|
‖x‖2‖y‖2 .
We assume x and y ∈ Rp have unit l2 norms and let [x, X2] be an orthogonal matrix. Then it holds | sin θ(x,y)| = ‖XT2 y‖2. (5) Furthermore, it can be checked that | sin θ(x,y)| ≤ ‖x−y‖2 ≤ √
2| sin θ(x,y)|. Thus | sin θ(x,y)| is a measure to quantify the distance between two vectors. In generalized eigenvalue problem, Crawford number is an important quantity for characterizing the perturbation bound.
Definition 2 The Crawford number of a symmetric matrix pair (A,B) is defined as
c(A,B) := min ‖x‖2=1
√ (xTAx)2 + (xTBx)2.
A symmetric matrix pair (A,B) is referred to as “definite” if c(A,B) > 0.
Obviously, the symmetric matrix pair (A,B) is definite for any positive definite B. By [25], it is known that the Crawford number is continuous with respect to the matrix pair.

heading: 3 Understanding of Perturbation Bound

heading: 3.1 On the General Setting
text: In order to understand the difference between the solution of (1) and the true leading eigenvector, we first describe the perturbation results without considering the sparsity. The following lemma is from [29] and is a generalization of the Davis-Kahn’s sin Θ theorem in [6] (also see [26, 15]). It gives an upper bound for | sin θ(u1, ũ1)|, where ũ1 is the leading eigenvector of (Ã, B̃). Lemma 1 Suppose (A,B), (Ã, B̃) = (A+E,B + F ) are both symmetric-definite pairs, and their eigenvalues satisfy λ1 ≥ · · · ≥ λp ≥ 0 and λ̃1 ≥ · · · ≥ λ̃p ≥ 0, respectively. Let φ1 = arctanλ1 > φ̃2 = arctan λ̃2, u1, ũ1 be the eigenvectors corresponding to λ1, λ̃1, respectively. Denote
Cu =
√ 2(‖A‖22 + ‖B‖22)
c(A,B) , ξ =
√ ‖Eu1‖22 + ‖Fu1‖22
c(Ã, B̃) ,
then
| sin θ(u1, ũ1)| ≤ Cu ξ
sin(φ1 − φ̃2) .
Therefore, it holds asymptotically that
| sin θ(u1, ũ1)| . Cu
c(A,B) sin(φ1 − φ2) . (6)
Although perturbation upper bound has been studied extensively, there are few results on lower bound especially in algebra literature. Next we establish the perturbation lower bound. That is, given any symmetric definite matrix pair (A,B) with λ1 > λ2 and any sufficiently small constant > 0, we can always find a matrix pair (E, F ) satisfying √ ‖E‖22 + ‖F‖22 < such that the distance between u1 and ũ1 is lower bounded by a quantity in the same order of the bound in Lemma 1. The result is stated in the next lemma. Lemma 2 Follow the notations in Lemma 1. For any small positive constant , it holds that
sup (E,F )∈F
| sin θ(u1, ũ1)| & Cl
c(A,B) sin(φ1 − φ2) , (7)
where Cl = c(A,B)√
2(‖A‖22+‖B‖22)κ 3 2
with κ := ‖B‖2‖B−1‖2, F := {(E,F )| √ ‖E‖22 + ‖F‖22 ≤ }.
Comparing (6) and (7) and noticing that Cu, Cl are two constants, we can see that upper and lower bounds of approximation error only differ up to a multiplicative constant when condition number κ is assumed to be bounded. Quantity /c(A,B) can be seen as the relative perturbation and 1sin(φ1−φ2) is a monotonically decreasing function of the gap φ1 − φ2. The upper and lower bounds are both proportional to c(A,B) sin(φ1−φ2) . Thus, we declare that the perturbation bound for the leading eigenvector u1 is rate-optimal.

heading: 3.2 On the Sparse Setting
text: Now we are ready to present the upper and lower perturbation bounds under the sparse setting. Throughout the rest of this section, the following assumptions are assumed.
A1 For any K ⊃ supp(u1) with |K| ≤ s+ k, it holds√ ‖ÃK −AK‖22 + ‖B̃K −BK‖22
c(AK, BK) < 1,
where c(AK, BK) is the Crawford number of (AK, BK) and s = | supp(u1)|.
A2 For any K ⊃ supp(u1) with |K| ≤ s+ k, ÃK and B̃K are positive definite. Assumption A1 requires that the perturbation within a small superset of u1 is tiny. It says that one can get a good approximation for u1 (according to Lemma 1 and 2) when a small superset of u1 is available. It is, in fact, a necessary condition for computing a good approximation for u1. Assumption A2 is a technical requirement to ensure positive definiteness for submatrices of B̃.
We further adopt the following notations. We let S := supp(u1) and define the perturbation level
:= max K:K⊃S;|K|≤s+k
√ ‖ÃK −AK‖22 + ‖B̃K −BK‖22,
the perturbation set F ,l := {(E,F )| √ ‖EK‖22 + ‖FK‖22 ≤ for any subset K with size less than l}.
Upper Bound of | sin θ(u1,x∗)|
Theorem 1 Let J = supp(x∗), K = S ∪ J and = K \ J . Denote ρ∗ = x T ∗ Ãx∗
xT∗ B̃x∗ , cK = c(AK, BK)
and c̃K = c(ÃK, B̃K). Let δ = ‖Ã(,J )[x∗]J−ρ∗B̃(,J )[x∗]J ‖2
c̃K . If δ < √ 1 + ρ2∗ and
arctan ρ∗ > arctanµ2 + arctan + arctan δ√
1 + ρ2∗ , (8)
then | sin θ(u1,x∗)| ≤ 1
c̃K ( C c(AK, BK) sin(φ1 − φ̃2) + δ Cδ sin(φ∗ − φ̃2) ) ,
where φ∗ = arctan ρ∗, µ2 is the second largest eigenvalue for (AK, BK), C and Cδ are some constants.
Remark: The upper bound has two terms. The first term is due to the perturbation, which is approximately proportional to the perturbation level . The second term is due to the failure in finding the true support set of u1. If supp(u1) ⊂ supp(x∗), then δ = 0, consequently, the second term vanishes. Therefore, it holds asymptotically that
| sin θ(u1,x∗)| ≤ Cu,K
c(AK, BK) sin(φ1 − φ2) (9)
with Cu,K = √ 2(‖AK‖22+‖BK‖22) c(AK,BK) , when supp(x∗) ⊃ supp(u1) and is sufficiently small.
Lower Bound of | sin θ(u1,x∗)| Next we present the lower bound for | sin θ(u1,x∗)|. It follows immediately from the proofs of Theorem 1 and Lemma 2.
Theorem 2 Follow the notations and assumptions in Theorem 1. Then the following result holds.
Case (a) If δ = 0, then
max (E,F )∈F ,s+k
| sin θ(u1,x∗)| & Cl,K
c(AK, BK) sin(φ1 − φ2) ,
where Cl,K = c(AK,BK)√
2(‖AK‖22+‖BK‖22)κ 3 2 2 (BK)
.
Case (b) If = 0, then
max (E,F )∈F ,s+k
| sin θ(u1,x∗)| & 1
√ 2κ 3 2 2 (BK)
ξ̂K sin(φ∗ − φ2) ,
where ξ̂K =
√ ‖Ê[u1]K‖22+‖F̂ [u1]K‖22
‖A‖22+‖B‖22 , and Ê, F̂ ∈ R|K|×|K| are some matrices controlled by δ.
Note that Crawford and condition number is well bounded on submatrix pair (AK,BK). Therefore, in Case (a), we can see that the lower bound is at the same order of the first term of the upper bound in Theorem 1; In Case (b), ξ̂K = O(δ), then the lower bound is at the same order of the second term of the upper bound in Theorem 1. We may declare that our perturbation bounds are rate-optimal.
Connection to the literature Although we do not make any statistical assumptions on E and F , we still make connections to the statistical literature here. First, it can be computed that
sin(φ1 − φ2) = λ1 − λ2√ λ21 + 1 √ λ22 + 1 .
When B = I is an identity matrix, Cu,K = O( √ λ21 + 1) and Cl,K = Ω(
1√ λ21+1 ). Thus our lower
bound can be simplified to √ λ22+1
λ1−λ2 . In sparse PCA [36], it has been established that the minimax
lower bound inf ũ1 supM(E‖u1 − ũ1‖22)1/2 is Θ( √ λ1λ2 λ1−λ2 ) with = √ s log p n andM is the model space where covariance matrix A satisfies that its principle singular vector is s-sparse and λ1 is larger than λ2 by certain margin. When λ1,λ2 and λ1/λ2 are bounded constants, this minimax rate matches ours. In [2], they provide the minimax perturbation bound for singular subspace of non-symmetric matrices. The minimax rate of ‖ũ1 − u1‖2 is
Θ( αz21 + βz12
α2 + β2 −min{z212, z221} ),
where the detailed definitions of α, β, z12, z21 can be found accordingly in [2]. Especially, for symmetric A and Ã, the above bound can be simplified to λ1−λ2 when λ1 − λ2. Again, when λ1 and λ2 are bounded constants, the rate is consistent with ours.

heading: 4 Understanding of Estimation Quality
text: In order to estimate the sparse leading eigenvector, it is straightforward to solve (1). It can be seen that the sparsity of solution to (1) depends on the choice of k, which may not be easily tuned. It is also computationally expensive to find the optimal solution to (1). Alternatively, such problem can be solved via regularization method. For example, [23] uses a smooth surrogate function to replace ‖ · ‖0 penalty for solving SGEP. A good estimator should have the following properties, 1. sparsity, 2. nearly unbiasedness, and 3. stability. In high dimensional statistical problem, B̃ is always singular when the sample size is smaller than the number of features. In order to obtain a good approximation of u1, we consider the following restricted problem
(P1′) min x∈Rp,‖x‖0≤sn −xTÃx + pλ(x) s.t. xTB̃x ≤ 1, (10)
where pλ(x) := ∑p j=1 pλ(x[j]) and pλ(x) is some univariate non-convex function, and sn is the restricted dimension. In practice, sn need not be a small number and it can grow with sample size n. The restriction ‖x‖0 ≤ sn is imposed for enforcing the solution to be nearly low-dimensional. Along with pλ, the estimator can recover the sparsity structure. This reformulation can be viewed as the counterpart of two stage methods [30, 9], where they need to find a good approximation of u1 in the first stage.

heading: 4.1 Penalization Function
text: For the choice of pλ, we consider a family of special non-convex penalties:
Pλ = {pλ(x) : pλ(x) satisfies (a1) - (a3)}, (11) where (a1) Function pλ(x) is an even function, i.e., pλ(x) = pλ(−x); (a2) The derivative of pλ(x), p ′
λ(x), exists in (0,∞); lim x↓0
p ′ λ(x) = λ and p ′ λ(x) ≡ 0 if x ≥ γλ for some constant γ; (a3) On
(0,∞), p′λ(x) is monotone decreasing and Lipschitz continuous, i.e., there exists a constant κ such that 0 ≤ p ′ λ(x1)−p ′ λ(x2)
x2−x1 ≤ κ for any 0 < x1 < x2.
Here (a1) requires pλ(x) to be symmetric; (a2) specifies the local property of derivatives of pλ(x) around 0 and assumes the flatness of pλ(x) for larger x; (a3) puts continuity constraints on p ′
λ(x). Many popular penalty functions are included in Pλ, for example, smoothly clipped absolute deviation (SCAD, [7]), minimax concave penalty (MCP, [43]) , etc. Based on the definitions, we can see that any pλ(x) ∈ Pλ is similar to `1-norm locally around the origin. On the other hand, pλ(x) puts a smaller penalization on the signal compared with `1 penalty. They do not give the penalization to those large values. Intuitively speaking, the regularized estimator with non-convex penalty may outperform the `1-norm-based estimator since it is nearly unbiased.

heading: 4.2 Support Recovery
text: In this section, we provide the error bound and oracle properties of the non-convex estimator under suitable conditions. Some additional notations are introduced as follows. We define û1 as arg maxx:‖x‖0≤sn,xTB̃x≤1 x
TÃx and denote ρ̃ = ûT1 Ãû1. Vector û1 can be viewed as the best approximation of the leading eigenvector in the restricted space. We define c(A,B, sn) := max|K|≤sn c(AK, BK) and c(Ã, B̃, sn) := max|K|≤sn c(ÃK, B̃K). We also define N(A,B, sn) := max|K|≤sn √ ‖A2K +B2K‖2. We define s := maxK:‖K‖0≤sn{‖EK‖2, ‖FK‖2}. A different set of assumptions is stated as follows.
B0 (Regularity) λ1 − λ2 is positive; ‖B‖, ‖B−1‖ are bounded by some constant. B1 (Signal) min{|u1[j]| : j ∈ supp(u1)} √ Sλ.
B2 (Penalization) λ s · CF , where CF is a constant depending on c(A,B, sn), c(Ã, B̃, sn) and N(A,B, sn).
B3 (Support size) |S|λ2 1. B4 (Identifiability) maxK:‖K‖0≤sn,S6⊂K λ1(AK, BK) < λ1 − 2 .
Condition B0 ensures that the leading eigenvector is unique and underlying matrix B is well-behaved. Condition B1 requires that the absolute values of entries in the true support is not too small so that S can be identified. Condition B2 specifies the relationship between penalty level λ and noise level . That is, penalty level λ should be at least larger than the noise level up to a multiplicative constant. Condition B3 makes sure that the support size is not too large, i.e., u1 should be sparse. The following theorem guarantees that the estimator is not far away from the re-scaled leading eigenvector. Condition B4 is for the identifiability of |S| to ensure that the support of u1 could be still identified after perturbation. Without loss of generality, we always assume sn ≥ |S|. Then we have the following results.
Theorem 3 Let x̂ be the optimizer of Problem (P1’). Under Conditions B0 - B4, it holds that min
sgn∈{−1,1} ‖sgn · x̂− u1s‖2 ≤ C1 sin(θ(û1,u1)) + C2 s,
for some constants C1 and C2.
Here u1s = u1(uT1 Bu1) −1/2 is the re-scaled leading eigenvector. Furthermore, we can recover the support of u1 under this restricted problem.
Theorem 4 Under the same conditions in Theorem 3, it holds that supp(x̂) = supp(u1). (12)
Based on Theorem 4, we actually have even stronger results. When all required conditions are met, the proposed estimator is equal to the oracle estimator which is the one estimated when the true support S is known to us. Therefore, the estimator achieves the optimal error bound which matches the one obtained in Theorem 1. Specifically, the proposed estimator achieves the optimal error rate ( |S| log pn ) 1/2 in sparse PCA and sparse CCA problems.
Theorem 5 Under the same set of conditions in Theorem 4, we have that
| sin θ(x̂,u1)| ≤ Cu,S c(AS , BS) sin(φ1 − φ2) , (13)
with Cu,S = C √ 2(‖AS‖22+‖BS‖22) c(AS ,BS) (C is a universal constant).
Remark 1 To achieve the optimal statistical error rate ( |S| log pn ) 1/2, most existing methods are two-stage based. For example, [38] studies sparse PCA problem and propose a “sparse orthogonal iteration pursuit" (SOAP) algorithm which involves “relax" and “tighten" stages; [9] considers sparse CCA problem and adopt the two-stage regularization approach with L1 penalty for the first stage and group LASSO penalty for the second stage. By comparison, our current estimator does not require convex relaxation in the first stage.

heading: 4.3 Computation
text: For computational purpose, we consider to solve (P1’) by using the alternating direction method of multipliers (ADMM, [37, 21]). Specifically, we relax the problem (P1’) by reformulating it to
min x,y,z
L(x, z,y), s.t. xTB̃x = 1, ‖x‖0 ≤ sn, (14)
where L(x, z,y) = −xTÃx + pλ(z) + yT(x− z) + η
2 ‖x− z‖22.
In (14), we introduce several auxiliary variables for the following reasons. We construct a new vector z which is a copy of x. It can help us to split the original problem to two simple separate sub-problems. y is the dual variable for the constraint x = z.
By formulation (14), we can optimize the objective function with respect to each variable iteratively. The steps for updating x,y, z are described as follows.
i Update x: At (t+ 1)-th iteration, we aim to find x(t+1) which is
arg min x∈D
η 2 ‖x− z(t)‖22 + (y(t))T(x− z(t))− xTÃx, (15)
where D = {x : xTB̃x = 1, ‖x‖ ≤ sn}. Solve sub-problem (15) to get x(t+1). ii Update z: We know that z(t+1) = arg minz η2‖z− x
(t)‖22 + (y(t))T(x(t+1) − z) + pλ(z). Each entry of z can be optimized separately. Specifically, if we take pλ as the MCP penalty, then z(t+1) has the following analytical form,
z(t+1)[j] = ž (t) if |ž(t)| > γλ,
sgn(ž(t)[j])(|ž(t)[j]|−λη )+ 1− 1ηγ if |ž(t)| ≤ γλ,
where ž(t) := x(t+1) + y (t)
η .
iii Update y: By dual variable update in [21], we have y(t+1) = y(t) + η(x(t+1) − z(t+1)).
We call the above procedure as non convex-SGEP (NC-SGEP) algorithm. Such proposed algorithm works on matrix-vector product and eigen-decomposition for submatrices of B̃, which is computationally efficient. Compared with other truncation methods, we do not need to make efforts to choose best sn due to the existence of non-convex regularization term. Thanks to the penalization term, our method can give a sparse estimator with a very wide range to choose the restricted dimension sn.
Here we propose a projection-based method for solving sub-problem (15). At (t+ 1)-th iteration, we aim to find x(t+1) which is the minimizer of (15). We construct the active set St and consider the following recursive formula,
b(t+1)m = z (t) − y(t) − Ãx(t+1)m−1 η , (16)
x(t+1)m [St] = (β(t+1)m B̃St + I)−1b(t+1)m [St],
where x(t+1)0 = x (t) and index m ∈ {1, 2, . . .}. Scalar β(t+1)m satisfies 1 = ∑ j dj(b̃ (t+1) m [j]) 2
(β (t+1) m dj+1)2
with
b̃ (t+1) m = UTb (t+1) m [St]; B̃St = UDUT and D = diag(d1, . . . , dsn). Such recursive formula is valid due to the following two observations (Propositions 1 - 2).
Proposition 1 Let y̌ be the projection of y on to the ellipsoid {x | xTBx = 1}. Then y̌ has the following form
y̌ = (βB + I)−1y, where β is a scalar which is the solution to the equation 1 = ∑ j dj(ỹ[j]) 2 (βdj+1)2 , where ỹ = UTy, B = UDUT and D = diag(d1, . . . , dp).
Proposition 2 The limiting point returned by (16) is the stationary point of (15).
Proposition 1 gives the explicit formula to project an arbitrary vector y to the convex body {x|xTBx = 1}. By contrast, it may lead to worse performance, if we just do the naive rescaling method, i.e., y̌ = y(yTBy)−1/2. Since (15) is highly non-convex, Proposition 2 only guarantees a way to find a stationary solution but not necessarily a optimal solution.
In this paper, we consider the following two possible constructions of active set St.
(C1) St is the set of indices corresponding to first sn largest absolute values of entries in b(t+1)1 (see (16)).
(C2) St is the set of indices corresponding to first sn largest absolute values of entries in x(t).
The following theorem gives the local convergence of NC-SGEP algorithm under additional assumptions. In general, ADMM-based method is extremely hard to analyze especially in highly non-convex problem. It remains an open question whether the global convergence result could be established.
Theorem 6 For the active set construction (C2), we take a large η value, set the initial support of x(0) includes supp(u1) and let z(0) = x(0),y(0) = 0. Then it holds
T ( ) ≤ C(L(x (0), z(0),y(0))− f̄) ,
where T ( ) = min{t : ‖x(t)− z(t)‖ ≤ } and f̄ = minx L(x,x,y) for some constant C which may depend on η and maxS:|S|≤sn ‖(BS)−1/2AS(BS)−1/2‖2. Moreover, supp(x(T ( ))) = supp(u1) if = o( s).
For construction method (C1), the same local convergence result also holds with an extra assumption that B̃ is diagonal.

heading: 4.4 Remarks
text: In practice, we find that construction (C1) has higher probabilities to find the global optimum compared with construction (C2). Especially in the application of sparse principle component analysis, construction (C1) is very efficient to recover the sparsity structure of leading component. Theorem 6 gives a local convergence result, which requires the initial value x(0) contains the support of true leading eigenvector u1. Such a good initial candidate of leading eigenvector could be obtained via using semidefinite programming (SDP)-based methods [35, 36, 9]. In other words, our proposed method can be easily merged to a two-stage-type method. The proposed algorithm only requires O(snp+ s3n) operations per iteration, while SDP-based methods have O(p
3) computational complexity. More discussions about SDP-based methods can be found in the supplemental material.

heading: 5 Numerical Experiments

heading: 5.1 Validation of Perturbation Bounds
text: We conduct perturbation analyses of proposed estimator under different settings. The matrix pair is set as A = 4u1uT1 + I − Pu1 , B = I . We sample n data which follows N(0, A) and sample another n data which follows N(0, B). The Ã and B̃ are constructed based on the sample covariance correspondingly. The leading eigenvector u1 has unit norm and has non-zero entries in first |S| positions. We fixed dimension p ≡ 100 and let number of samples (n) and support size (|S|) vary. Each setting is repeated for 100 times with fixed choice of λ = 0.3, η = 1, sn = 25. The mean and
standard deviation of ‖x̂ − u1‖2 are reported. We also compute the oracle estimator which is the best x̂ when true support S is known. From Table 1, when noise level √ log p/n is small, we can see that the estimation error is quite close to the optimal (oracle) error. In addition, we can see that estimation error is proportional to √ |S| and √ log p/n. This indicates that our method can achieve
the theoretical optimal error rate, i.e., O( √ |S| log p
n ).

heading: 5.2 Validation of Sparsity Recovery
text: The underlying matrix pair is set as A = 3u1uT1 + I − Pu1 , B = I . Thus, we naturally take B̃ = B = I and Ã as the sample covariance of data. We let dimension p grow from 16 to 256, fix the sample size n = 100 and set λ = 0.5, η = 1 and sn = 50. We compare the proposed method with the semidefinite programming method [34, 18, 40, SDP] with `1 (SDP_L1) and MCP (SDP_MCP) penalty. The estimation error and percentage of support recovery are reported in Table 2. We can see that the proposed method can have a slightly better performance. This is because that NC-SGEP optimizes objective within space Rp unlike those SDP methods work on space Rp×p. In addition, the proposed method can recover the sparsity structure pretty well, while SDP methods can never recover the true support set.

heading: 6 Conclusion
text: In this paper, we establish the upper and lower bounds for perturbation analysis of sparse generalized eigenvalue problem. We also consider a new statistical estimation method. The proposed method gives a sparse, nearly unbiased and stable solution to SGEP. We show that the proposed estimator can achieve the optimal estimation error rate. We further present a non-convex SGEP (NC-SGEP) algorithm to solve a non-convex regularization problem with guarantee of local convergence. Multiple numerical results validate our theories and also show the superior performance of the proposed method. In the future work, on the theoretical side, we may focus on extending the current results to the problem of finding multiple leading sparse eigenvectors and establishing the corresponding new lower bound theory. On the computational side, even though computing leading vectors one by one seems straightforward, it is not an easy problem since it requires the orthogonalization procedure. The computational complexity may require a special care.

Review summaries:
The paper studies a class of sparse generalized eigenproblems which encompasses several common problems in machine learning. For this class of problems it gives upper and lower perturbation bounds which are shown to be rate-optimal. The authors then present an optimization problem with a non-convex penalty function and show that it can be used to estimate the sparse leading eigenvector. The theoretical results are validated numerically using synthetic data in the experimental section.   

The paper is dedicated to the sparse generalised eigenvalue problem (SGEP), which consists in recovering the leading generalized eigenvector corresponding to the largest generalized eigenvalue of the matrix pair (A,B). The authors consider the noisy version of A and B: \tilde{A} and \tilde{B}. The main challenges are: problem NP hard because of the subset selection problem (the sparse support size is unknown); the noisy matrices may be ill-conditioned to that algorithms requiring their inverse may not work. The contributions of the authors are as follows: derivation of an error bound between the true sparse eigenvector and the solution of the SGEP; derivation of the conditions for support recovery for potentially non-convex penalty functions; an ADMM based algorithm is proposed to solve the SGEP.