Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
Authors: Jian-Feng Cai
Zhuozhi XIAN
Jiaxi Ying
Venue: ICML 2025 poster (2025)
Keywords: Sparse PCA
principal component analysis
single-spiked covariance model
sample complexity
Primary area: probabilistic_methods

Abstract: We explore the single-spiked covariance model within the context of sparse principal component analysis (PCA), which aims to recover a sparse unit vector from noisy samples. From an information-theoretic perspective, $O(k \log p)$ observations are sufficient to recover a $k$-sparse $p$-dimensional vector $\mathbf{v}$. However, existing polynomial-time methods require at least $O(k^2)$ samples for successful recovery, highlighting a significant gap in sample efficiency. To bridge this gap, we introduce a novel thresholding-based algorithm that requires only $\Omega(k \log p)$ samples, provided the signal strength $\lambda = \Omega(||\mathbf{v}||_\infty^{-1})$. We also propose a two-stage nonconvex algorithm that further enhances estimation performance. This approach integrates our thresholding algorithm with truncated power iteration, achieving the minimax optimal rate of statistical error under the desired sample complexity. Numerical experiments validate the superior performance of our algorithms in terms of estimation accuracy and computational efficiency.

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
Jian-Feng Cai Zhuozhi Xian
Jiaxi Ying !
Abstract
et al., 2006; Shen & Huang, 2008), which offers robust alternative to traditional PCA.
We explore the single-spiked covariance model within the context of sparse principal component analysis (PCA), which aims to recover a sparse unit vector from noisy samples. From informition-thhorrtiic perspective, ()) observations sufficient recover k-sparse p-dimensional vector v. However, existing polynomial-ting time methods require at least O(k²) samples for successful recovery, highlighting a significant gap in sample efficiency. To bridge this gap, we introduce novel thresholding-base based algorithm that requires only 12k log p) samples, provided the signal strength À = (. We also propose two-stage nonconvex algorithm that further enhances estimation performance. This approach integrates our thresholding algo- rithm with truncated power iteration, achieving the minimax optimal rate of statistical error under the desired sample complexity, Numerical experiments validate the superior performance of our algorithms in terms of estimation accuracy and computational efficiency.
In this paper, we investigate the single-spiked covariance model (Johnstone, 2001), which involves recovering a sparse vector from n noisy samples:
33 Vxg. + Eis i i=
where æ, E Rp denotes the i-th observation and v e R' is the unknown vector to be estimated, characterized by |/oloo S k and ||o||2 1. The coefficients 9i are sampled independently from standard normal distribution. (0, 1), and the noise vectors are likewise independently drawn from multivariate normal distribution A(,, I), with each Ii and &, being mutually independent. The parameter À > 0 quantifies the strength of the signal.
From an information -heeoorttcc perspective, it has been shown that n Q logp) samples are sufficient to estimate k-sparse p-dimensional vector up constant error (Vu & Lei, 2013; Berthet & Rigollet, 2013). However, achieving such an estimate are nonpolynomial time algorithms such as exhaustive search. In contrast, all known polynomial- time algorithms, including diagonal thresholding (Johnstone & Lu, 2009), covariance thresholding (Krauthgamer et al., 2015; Deshpande Montanari, 2016), and convex relaxation techniques (d' Aspremont al., 2004; Amini & Wainwright, 2009; Ma & Wigderson, 2015), need at least ?(A22) observations to estimate v. Moreover, the known lower bounds for these methods are considered tight (Krauthgamer et al., 2015; Ma & Wigderson, 2015; Deshpande & Montanari, 2016; Gao et al., 2017) when the signal strength is treated as a constant.
1. Introduction
Principal component analysis (PCA) is cornerstone technique for dimension reduction, data preprocessing, clus- tering and feature extraction. However, in dingh-ddmmsiinal regimes, its dependence on all variables often leads to overfitting and poor interpretability, obscuring the truly important features (Johnstone Lu, 2009). By introducing sparsity into principal components, sparse PCA can mitigate issues of overfitting and loss of interpretability that often arise in high- digeesiinni settings, particularly when the number of variables surpasses the number of samples (Zou
Thus, for A ?(log p). there exists notable gap between the number of samples needed from an cetrmaation-thhored perspective, 52(k: log p). and the number required by existing polynomial-time algorithms, 2(k22). Moreover, reductions from the planted-clique conjecture imply that, without further assumptions, no polynomial-time algorithm can attain the cotommtiion-hhorettii rate (Berthet & Rigollet, 2013; Krauthgamer et al., 2015; Wang et al., 2016; Gao et al., 2017; Brennan et al., 2018). This motivates the key question of our work:
'Department of Mathematics, Hong Kong University of Science and Technology, Hong Kong. Correspondence to: Jiaxi Ying <jx.ying www.yinn@@oomnnctuus......
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
Can we design a polynomial-time algorithm to bridge this sample-complexity gap, under certain mild assumptions on the model?
A notable gap exists between the information-theortiic sample limit of S2(k logp) achieve a constant estimation error and the sample complexity of (2(k2) required by existing polynomial-tin timee algorithms. We bridge this gap by introducing novel polynomial-time algorithm based on the thresholding technique, which requires only ( log p) samples to achieve constant error, provided the signal strength
We answer this question affirmatively by designing a novel thresholding- based algorithm. We demonstrate that Okk lo is sufficient for our algorithm to estimate o up to constant error (in euclidean norm), given that the signal strength
). By imposing a condition on A and the maximum absolute value of v, we effectively reduce the parameter space of the model. Consequently, our algorithm achieves the desired sample complexity within this constrained parameter space. While the theoretical analysis of sparse PCA conditioned on the signal strength has been studied the literature (Johnstone & Lu, 2009; Deshpande & Montanari, 2016; Novikov, 2023; Ding al., 2024), these studies predominantly focus scenarios with relatively large sparsity k and correspondingly high sample complexities, such as n s(p). In contrast, our work aims to achieve a reduced sample complexity of n = Q(k p), operating under a specific condition on A.
In existing work, the signal strength often treated as constant when deriving both the information- theoretic lower bound of S2 logp) and the sample complexities of standard algorithmms (c.g., Q(k² log for diagonal thresholding and S2(k2) for covariance thresholding). Our algorithm achieves the optimal sample complexity of 22 log provided that À = ([] Under this condition, remains constant provided that ||0/|>> = 2(1). This scenario is naturally satisfied by signals whose nonzero entries exhibit a power-law decay (Jagatap & Hegde, 2019), This property is known as compressibility in the compressive sensing literature (Donoho, 2006; Candes et al., 2006).
Although the methods previously mentioned, such as diagonal thresholding (Johnstone & Lu, 2009), covariance thresholding (Krauthgamer et al., 2015; Deshpande &
To enhance estimation performance, we propose two-stage nonconvex algorithm. The first stage uses our thresholding- based approach for initialization, followed by a refinement through truncated power iteration. We show that the total estimation error of our algorithm can be separated into two components: an optimization error, which decays at a linear rate, and a statistical error, which achieves the minimax optimal rate when Q(k log and X ((Ibol))))
Montanari, 2016), and semidefinite programming (SDP) (d' Aspremont et al., 2004; Amini & Wainwright, 2009; Vu et al., 2013; Ma & Wigderson, 2015), are all computable in polynomial time, the SDP approach entails a significantly higher computational burden compared to thresholding- based methods. For instance, the SDP implementation by (d' Aspremont et al., 2004) incurs computational complexity exceeding O(p"), rendering it intractable in the big data scenarios. Among the thresholding- based algorithms, covariance thresholding (Krauthgamer et al., 2015; Deshpande& & Montanari, 2016) generally offers better estimation performance than diagonal thresholding (Johnstone & Lu, 2009). However, it requires eigendecomposition of typically p X p matrix, costing Onp" + +y")) totally, whereas diagonal thresholding involves eigendecomposition of smaller l x k matrix, leading to overall computational cost of (np + ). In this paper, we aim to develop an algorithm that achieves both high estimation accuracy and computational efficiency.
Numerical experiments demonstrate that our proposed algorithm offers significant advantages in terms of both estimation accuracy and computational efficiency. In contrast, the diagonal thresholding method results in substantial estimation errors, while the two covariance thresholding approaches incur considerable computational costs.
Notations: ||æ/|oo counts the nonzero entries of vector æ, while ||>|| and || denote its lu and foo norms respectively. 27 represents the vector of the same length as æ. where elements indexed by T are retained, and elements indexed by T are set zero. [xi denotes the sub-vector of x containing only the elements indexed by T. AT denotes the matrix of the same size a square matrix A, where columns and rows indexed by T are retained, and those indexed by T are set to zero. [AlT refers to the sub- matrix of A composed of the columns and rows indexed by T. We write f(n) = O(g(n) means fn) S cg(n) and f((m) = Mgg(n) means f(n) 2 d'g(n) for positive constants and d' and use [p] to denote the set I.....p.
1.1. Contributions
We introduce two complementary algorithms to address both sample complexity and estimation performance the single- spiked covariance model. The first is a thresholding-base polynomial-time algorithm that bridges the existing gap in sample complexity. The second is two-stage nonconvex algorithm, which leverages our thresholding algorithm combined with truncated power iteration to further refine estimation. Our main contributions are:
2

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
2. Preliminaries and Related Work
S2 log (Amini & Wainwright, 2009), which is worse than the information theoretical bound S?(k log To attain lower sample complexity, covariance thresholding, in contrast, involves a more intricate process. It applies either a soft thresholding (Deshpande & Montanari, 2016) or hard thresholding (Krauthgamer et al., 2015) to the empirical covariance matrix Σ, followed by using the indices of the top k entries of its principal eigenvector to estimate the support of v. It has been proved that the sample complexity of covariance thresholding proposed by (Deshpande & Montanari, 2016) is (2(22, but its computational complexity is O(pp² + +p") as analysed in Section 1.
In this section, we begin by introducing the preliminaries of sparse PCA, then present related work.
2.1. Preliminaries
The single-spiked covariance model, initially introduced by (Johnstone, 2001), has become a foundational framework for sparse PCA. It outlines a methodology to estimate a single sparse vector e from n noisy observations, denoted as ....... Zn E IRP. Each noisy sample æ; is independently drawn as in (1) and is therefore zero-mean.
We compute the empirical covariance matrix as Σ 1 E a. and its expectation is given by (Ma & Wigderson, 2015; Krauthgamer et al., 2015):
Convex relaxation methods, such as SDP relaxation (d' (d'Aremmonn et al., 2004) and degree-4 Sum-of-Squares (SoS) algorithm (Ma & Wigderson, 2015), offer an alternative approach. They reframe the non-convex sparse PCA optimization problem as a convex optimization problem, so it can be solved in polynomial time. SDP relaxation rewrites Problem (3) by introducing a new variable W = ww and then relaxes the fo-norm constraint by a convex constraint. Its sample complexity is S2(k² log (Amini & Wainwright, 2009) and its computational complexity exceeds O(p") (d' Aspremont et al., 2004). The degree-4 SoS algorithm applies the well- established SoS method (Lasserre, 2015) to the sparse PCA optimization problem. which requires larger computational costs than SDP but proves to be unable to improve its sample complexity (Ma & Wigderson, 2015).
E Σ = \vo +I.
(2)
Given an empirical covariance matrix E Rexp, the objective of sparse PCA is identify sparse unit vector
w E Iº that maximizes the quadratic form w" Ew.
Formally, this can be expressed as:
max w Ew, subject to ||l||| 1, [/w]lo llo S k. (3) W
Although the solution to (3) yields an estimator for the spike vector v, its combinatorial fy-norm constraint renders the problem NP-hard. In this paper, we address these challenges by proposing polynomial-time algorithm that bridges the existing sample-complexity gap and computationally efficient nonconvex algorithm that achieves the minimax optimal rate under the desired sample complexity.
3. Main Results
In this section, we begin by introducing two algorithms in Section 3.1. Following this, we examine their computational complexity in Section 3.2. Finally, in Section 3.3, we present a detailed theoretical analysis of both algorithms. This analysis includes establishing the sample complexity required for the thresholding algorithm to achieve constant error, thereby bridging the existing gap regarding the sample complexity. We also provide bounds for the estimation error of the two-stage algorithm, detailing how it enhances performance through refined estimation process.
2.2. Related Work
We review related work, specifically thresholding based algorithms (Johnstone & Lu, 2009; Krauthgamer et al., 2015; Deshpande & Montanari, 2016) and convex relaxation techniques (d' Aspremont et al., 2004; Ma & Wigderson, 2015), which have made significant strides in solving the sparse PCA optimization problem.
Thresholding- based algorithms for sparse PCA have been widely adopted due to their simplicity and effectiveness. These algorithms operate by initially estimating the support of v and subsequently computing the principal eigenvector of the corresponding submatrix of Σ approximate v. Notable among these are diagonal thresholding (Johnstone & Lu, 2009) and covariance thresholding algorithms (Krauthgamer et al., 2015; Deshpande & Montanari, 2016). Diagonal thresholding leverages the indices of the top diagonal entries of the empirical covariance matrix S as an estimate of the support of v. With Ompp + m)?) computational complexity, its sample complexity is
3.1. Proposed Algorithms
The first algorithm is a thresholding-basee algorithm. an advanced version of diagonal thresholding. The second algorithm is a nonconvex two-stage algorithm for enhancing the estimation performance.
3.1.1. THRESHOLDING- BASED ALGORITHM
Our proposed algorithm extends the well-known diagonal thresholding method (Johnstone & Lu, 2009), which, while being one of the most computationally efficient approaches,
3

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
suffers from a high sample complexity of Q(k2 log p) (Amini & Wainwright, 2009). This complexity exceeds the information-theeortie lower bound of (k log To address this discrepancy, we first delve into the operational principles of diagonal thresholding.
Our method is motivated by the fact that, when / S, the expected values of the entries in the 1-th column of y can be expressed as:
dep +1, j=l, jES 14, 0, jES,
E Eeill
(6)
Diagonal thresholding traditionally estimates the true spike by identifying the support and magnitude of the non non-zero entries separately. This method estimates the support of v by using the indices of the top k diagonal entries of Σ in absolute value. This approach is based on the observation that the expectation values of the diagonal entries of Σ are given by:
where {e denotes the set of standard basis vectors in R Similar to (5) for diagonal thresholding, we now consider the following statistical gap:
9c = min JES
[ E [ei]]
max JES
E E
[ei]l,
Au +1, JES, jES,
E [
(4)
2 A lerl min |uyl
1,
JES
where S denotes the support of v. There exists gap between diagonal entries of the expectation in the sets S and S This statistical gap is defined by:
Upon comparison, we observe that the lower bound of the statistical gap in (7) consistently exceeds that in (5) when | E S. This indicates that our method more effectively distinguishes entries in S from those in S than diagonal thresholding, thereby suggesting a reduced sample complexity requirement (see Section 3.3 for details).
9a := min JES
-
max jES=
[E[ 5s| = A min ug.
[
jES
(5)
A simple analysis for diagonal thresholding is based on the following proposition about the statistical gap Sa-
To maximize this lower bound, I should ideally be selected as the index of the largest absolute element of v, which aligns with the index of the largest diagonal entry of T Σ as shown in (2). However, since we only have the sample covariance matrix Σ rather than the population covariance E[E]. we designate / as the index of the largest diagonal entry of E, denoted by ju-
Proposition 3.1. 1f| / E. - [R[ aal S 1g holds for every i E [p], then, for any j E S and any j' E S", we have:
أッ 2 y.
The proof of Proposition 3.1 is detailed in Appendix A.2. Given sufficient samples, the concentration error E.
Our thresholding-basel algorithm is summarized in Algorithm and consists of two main steps. First, we estimate the support of the spike vector v by selecting the indices corresponding to the k largest absolute values in Σej,, which we denote as S. Second, we estimate via a procedure analogous to diagonal thresholding; we perform an eigendecomposition on the k x k submatrix [É]g formed by restricting to rows and columns indexed by S.
[E [E]],,] is small enough and the condition in Proposi- tion 3.1 holds, hence the diagonal entries Edg for indices j E S are larger than those for indices E S in absolute values. Consequently, diagonal thresholding can estimate the support of by selecting the indices of the top k diagonal entries of Σ in absolute value. Additionally, the sample complexity of diagonal thresholding is governed by the statistical gap Ja- larger gap Iu permits greater tolerance for the concentration error [ E [S]] m] thereby reducing the number of samples required.
The support-estimation step is critical for achieving the optimal Q log p) sample complexity. By amplifying the separation between in-support and out-of-support entries relative to diagonal thresholding. our method recovers the true support with fewer samples. This enhanced separation constitutes the key innovation of our approach. As an alternative spike-estimation step, one may simply normalize the jo-th column of E f]g- This variant avoids the eigendecomposstion for faster runtime, at the expense of a modest loss in estimation accuracy, yet it still achieves the optimal ((k logr) sample complexity.
Based on the above observations, we aim to distinguish between two groups of entries corresponding to the sets S and S", aiming for a larger gap between them. Our proposed thresholding algorithm draws inspiration from spectral initialization methods designed for sparse phase retrieval (Wu & Rebeschini, 2021; Cai et al., 2023). Instead of focusing the diagonal entries, examine the magnitudes of the entries in a specific column of E. More precisely, we identify the support of v by selecting the indices of the k largest entries in the jo-th column of E in absolute value, where jo is the index of the maximum element among the
3.1.2. TWO-STAGE ALGORITIIM
To enhance the estimation performance of our thresholding algorithm presented in Section 3.1.1, we propose a two-stage
diagonal entries of E
4

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
Algorithm 1 'Thresholding- based algorithm for bridging the sample-complexity gap
The convergence of truncated power iteration was established in (Yuan & Zhang, 2013), where was shown to achieve the optimal statistical rate. However, it hinges on the assumption that the initial estimate va satisfies ( w)] < 1 - C with C' E (0, 1) a constant. Such an initial point was not provided in (Yuan & Zhang, 2013).
Input: Samples (ax))=1- the sparsity k Compute ( و and set jo = arg maxisisn Ed
Compute Eejo and set S as the indices of the top k elements of Σejo in absolute value.
Our proposed thresholding algorithm introduced in Section 3.1.1 tackles this problem perfectly, which achieves the desired sample complexity. Building this foundation, we propose two-stage algorithm to solve Problem (3). The full algorithm is summarized in Algorithm 2.
Compute [Σك g. set [00) as the unit leading eigenvector of [ and set [b0]s = O.
Output: v
Algorithm Two-stage algorithm for enhancing estimation
Input: Samples {.; Mi=1 the sparsity k, parameter k' // Initialization stage:
algorithm comprising both initialization and refinement stages. In the initialization stage, we utilize thresholding algorithm to obtain an initial estimate. During the refinement stage, the estimate is refined by employing a well-known method called truncated power iteration (Yuan & Zhang. 2013), which is widely used for finding a sparse eigenvector of a matrix. By combining these two stages, our algorithm aims achieve superior performance in addressing Problem (3).
Compute an initial estimate v by Algorithm 1.
// Refinement stage:
for t = 1,2, do
J' = Tx ( Ei= (x)
0r )z;)
v' 8' 1/0012 ll2
end for Output: o²
The truncated power method (Yuan & Zhang, 2013) extends the classic power commonyn used to extract a matrix's leading eigenvector-by by inserting a sparsity- enforcing truncation at each step. For reference, the standard power updates are
3.2. Computational Costs
The computational expenses associated with our thresh- olding algorithm, as proposed in Section 3.1.1, include: Computing the diagonal of E requires ((n)) operations: Identifying the largest diagonal entry involves O(p) operations; Computing Σe,, requires O(np) operations; Finding the indices of the top k elements of Eejo in absolute value requires Opp logp) operations: Computing the submatrix [ requires O(k²²) operations; Computing the leading eigenvector of [ demands ((1a)) operations. Combining these, the overall computational complexity of our thresholding algorithm totals O(np + nkk)..
-'=
Ev
(8)
where of is the intermediate update at t-th iteration and o' its unit-norm normalization. The iteration (8) typically converges to dense vector under suitable conditions. However, since the sparse PCA problem (3) seeks sparse principal cigenvector, we employ the truncated power method instead of the traditional power method. The truncated power method modifies each iteration by truncating to zero all but the top k' entries in terms of absolute value:
Remarkably, this computational complexity is equivalent to that of diagonal thresholding (Johnstone & Lu, 2009), yet significantly less than that of covariance thresholding (Krauthgamer et al., 2015; Deshpande & Montanari, 20, which scales as Onpp² + p²). Moreover, we will show in Section 3.3 that our approach offers improved sample complexity compared to both diagonal thresholding and covariance thresholding.
3' = -Te v²
where (-) represents the truncation operator, which retains only the top k'' entries of a vector absolute values and sets all others to zero. This modification makes the truncated power method a natural choice for Problem (3), as it inherently promotes sparsity in the resulting vector. It is common to set the truncation parameter k' in the same order as the sparsity k of the true spike. Moreover, since
The two-stage algorithm employs our thresholding algo- rithm, costing O(np + nık²), and uses the truncated power iterations with O(np) in each iteration. Given the limited number of iterations typically required for convergence, the overall computational expense of the refinement stage does not significantly increase the total cost of the algorithm.
S = 1 =1 x Ex in (9) can be simplified as
" ( - -))
5

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
Under the assumption À = [[ from Theorem 3.2, diagonal thresholding (Johnstone & Lu, 2009) cannot attain the optimal ( logp) sample complexity. The reason is that diagonal thresholding fails to exploit the decay effectiveness is quantifying the estimation error. To this end, we define a loss function that captures the error between the true vector e and its estimate 6:
3.3. Theoretical Results
In this section, we delve into the theoretical underpinnings of both the thresholding algorithm and the two-stage algorithm. A critical aspect of assessing these algorithms'
structure of v in estimation, whereas our method does. By Proposition 3.1, a larger statistical gap implies lower sample complexity. When |/ has constant order, our method's gup in (7) depends only linearly on minjes |tyl. yielding an order-wise improvement over diagonal thresholding' gap in (5); diagonal thresholding cannot exploit the constant-order |/ollso and thus remains fundamentally limited.
dist(v, 6) == min { 6112, 110 + ö1|2).
(10)
This loss function chosen for the sign ambiguity of 6. Specifically, if é is a solution of the sparse PCA problem only be estimated
Theorem 3.3 (Two-stage algorithm). Let Rp be k- sparse unit vector Let I; = vdg.u + E,, i= 1, 1...., n. where X > 0, 9i
(3), - is also solution of (3), up to sign.
~ NN(,, and & idd N() are independent. There exist universal constants C3, Ca. Cs > O such that if Z C3||ee||a and n 2 Cak logp. then with probability exceeding 1 5y- the output o' of Algorittm 2 with parameter k' Csk and an initial estimate vd generated by Algorithmm 1 satisfies
iid
Theorem 3.2 (Thresholding algorithm). Let v E Rp be a l-- sparse unit vector. Let x, = Vdg.e + E., i n,, where X > 0, 9i idd 10. and & ~ NN I,) are
independent. For any 7 E (0, 1]. there exists universal constants C1, C2 > 0 such that if x 2 Cillv|laa and n 2 -2y->2loog log p. then with probability exceeding 1 - 5p the output Dd of Algorithm satisfies dist(v, w²) S v.
dist((" v) S (D, v0) + +d' dd'kllog/y (1))
Optimization error
Statistical
The proof of Theorem 3.2 is detailed in Appendix A.3. Theorem 3.2 demonstrates that our proposed algorithm requires only n = S p) samples to estimate the true spike v up constant error, given that the signal strength = [[ ( This result is significant as it effectively bridges the gap between the enformation-theoreet lower bound of Q(k lo necessary to achieve a constant estimation error, and the considerably higher sample complexity of 12(k2) demanded by existing polynomial- time algorithms. Thus, our proposed method not only aligns with the optimal information-theoree sample limit but also significantly reduces the required sample size compared to other contemporary polynomial-tin techniques.
where 0 < d < 1 and d' >> 0 are constants.
The proof of Theorem 3.3 is detailed in Appendix A.4. Theorem 3.3 shows that the estimation error between the estimated and true spikes is bounded by two terms, i.e., the optimization error and statistical error. The optimization error, d' dist(v, vo), decays to zero at linear rate with respect to the iteration number t. The statistical error
of our algorithm is on the order of v /Aloke, achieving the minimax optimal rate (Vu & Lei, 2013; Wang et al., 2014; 2016). Notably, the algorithm attains this optimal rate under the desired sample complexity of 22kk: log p). provided X = ((llll)) We can see the statistical error is independent of t, implying that it will not decrease during iterations in the algorithm.
n
Attaining the optimal sample complexity Ω(k log in polynomial time requires extra structural assumptions. Indeed, assuming the planted-clique conjecture holds, reductions show that no polynomial-time algorithm can recover the spike with only 9(k log samples without additional conditions (Berthet & Rigollet, 2013; Krauthgamer et al., 2015; Wang et al., 2016; Gao et al., 2017; Brennan et al., 2018). We assume À = [2([v/]==') ). In much of the SCArse-PCA literature, À is treated as a constant when deriving both information-theoretie lower bounds and algorithmic cample-coppexxity guarantees. In settings where the nonzeros of follow a pawer-law decay- common model in compressive sensing (Donoho, 2006; Candes et al., 2006)-> naturally remains of constant order under this assumption. An analogous phenomenon appears in sparse phase retrieval (Wang et al., 2017; Cai et al., 2022; 2024), where power-law decay likewise enables optimal (kk l) recovery (Jagatap & Hegde, 2019).
4. Experimental Results
We validate our two-stage algorithm in MATLAB on 2.10 GHz Xeon Gold 6152 machines, confirming the theoretical results and demonstrating its efficiency. Estimation error is measured as
Error = dist(v, 6),
and support recovery by
2tp
F-score- 2tp + fp + fn'
where tp. fp, and fn denote the counts of true positives, false positives and false negatives; The F-score ranges from 0 to 1, with indicating perfect support recovery.
6

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
1000
Hard
0.6
Estimsion ertr
Estiaaiin
0.3
0
0.5
1.5
2.5
35
Dimension
Elogp
(a) Comparisons of estimation error
(a) /P
DT
505
SotCCT Hard-CT
1000
3
5500
p-10000
0.4 Esiialion
2.0
0
3.5
010
2000
5000
Dimension
Elog
(b) Comparisons of computational time
(b) k 0.02p
Figure 2. Estimation and computational time dimen- sion for different methods. Our algorithm demonstrates both estimation accuracy and computational efficiency.
Figure 1. Estimation error versus scaled sample size for our algorithm. We consider two sparsity settings: (a) 0.5 and (b) & 0.02p. Curves are averaged over 200 executions.
Figure shows that the estimation errors are nearly identical for different values of p. indicating that Ct.p does not depend on p. Moreover, since is function of y in our settings, Ck.p does depend on either. Therefore, Ck.pp is independent of both p and k. Consequently, empirical evidence shows that the estimation error of our algorithm is O( p/g)/n), provided the sample size 2 2k log
4.1. Validating Theoretical Results
We provide empirical evidence that our algorithm achieves an estimation eиror of the order /k logp/n, provided the sample size čk log p for some constant č.
In experiments, we set the dimension to range from 500 to 10000, with each curve corresponding to specific set to 0 in (a) and 0.02p
4.2. Demonstrating Estimation Performance
dimension. The sparsity in (b). We generate the true spike e with nonzero entries, where the support of the entries is chosen randomly, and the value of each entry either 1/ l -1/Vk. We then generate i..dd samples æ; from the single-spike model. The signal strength À is set to 10.
We demonstrate that algorithm surpasses state-of-the art methods in both estimation accuracy and computational efficiency. Diagonal thresholding (DT) (Johnstone & Lu, 2009) tends to produce large estimation errors, whereas covariance thresholding (CT) methods, including soft- thresholding (Soft-CT) (Deshpande Montanari, 2016) and hard-thresholding (Hard-CT) (Krauthgamer et al., 2015). require substantial computational resources. In contrast, our proposed algorithm consistently achieves superior performance in terms of both lower estimation error and reduced computational time.
Figure plots the estimation error of our algorithm as a x
n function of the scaled sample size &logp' When Flogp 2 2,
the curves corresponding different dimensions align well with each other, indicating they achieve nearly the estimation error. For instance, when Flng 2, we have
2k logpp Without loss of generality, we consider the estimation eттоr to be / lo. p/n, where Ck,p can be a function of and By substituting n 2 log p, the estimation error simplifies to p/.tp/2..
Figure 2 compares the estimation error and computational time different methods across various dimensions. We set the sparsity & 0.02p and the sample size 3k log p.
/c,,p/

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
1.2
4.3. Examining Computational Performance
SONCT
Haat-CT Proposed
We assess the scalability of our proposed algorithm by measuring its running time and iteration count functions of dimension p, sample size n, and sparsity k. In each experiment, we vary one of these parameters while holding the other two fixed to isolate its impact on per-iteration cost and convergence behavior.
0.8
Estmatin
0.6
0.4
Table reveals the algorithm's response to increasing dimension p. The results demonstrate that while the per- iteration computational time exhibits near-linear growth with dimension, the iteration count remains remarkably stable. This stability in iteration count is particularly noteworthy, as it results in a total running time that scales approximately linearly with p. Such behavior strongly indicates that our method maintains its computational efficiency even in high-dimensional regimes, a crucial advantage for modern large-scale applications.
0.2
300
300
500
700
5 Sample size
(a) À = 5
0.7
DT
SobCT Hart-CT Procosed
0.8
Ezimoton cer
Dimension
5000
7500
10000
12500
15000
Number of iterations
9.05
9.06
9.14
9.13
9.13
Time for an iteration(s)
0.018
0.029
0.040
0.048
0.055
Total running time(s)
0.165
0,264
0.367
0.440
0.506
0.1
Table I. Impact of dimension p on computational time in the
Sample size
refinement stage of Algorithm 2, with 2500 and k: 100.
(b) x = 10
presents an interesting interplay between sample size n and computational efficiency. While larger sample sizes induce linear increase in per-iteration computation time, they simultaneously yield a compensatory decrease in the required number of iterations. This trade-off results in only modest increase in total running time as grows. The reduction in iteration count can be attributed to the enhanced statistical accuracy of initial estimates with larger sample sizes- more samples provide a more precise initialization closer to the optimal solution, thereby accelerating convergence.
Table 2
Figure 3. Estimation error versus sample size for different methods under the settings of (a) À = 5 and (b) À= 10. The dimension is set to p- = 1000. Error bars indicate 95% confidence intervals for the mean estimation error, calculated from 200 independent runs.
As the dimension p increases from 1000 to 10000, we observe that the computational time for both Soft-CT and Hard-CT increases dramatically, while DT and our proposed method show steady increase and take significantly less time. Furthermore, while Soft-CT, Hard-CT, and our algorithm achieve substantially smaller estimation errors compared to DT, our algorithm consistently obtains the smallest estimation error among all methods.
Sample size
2000
2500
3000
3500
4000
Number of iterations
9.38
8.92
8.47
7.53
5.75
In summary, Figure shows our algorithm excels in both estimation accuracy and computational efficiency. In contrast, the diagonal thresholding method incurs large estimation errors, while the covariance thresholding methods are computationally intensive.
Time for iteration(s)
0.014
0.019
0,022
0.026
0.030
Total running time(s)
0.133
0.167
0.191
0.195
0.175
Table 2. Impact of sample size n on computational time the
refinement stage of Algorithm 2, with p = 5000 and k = 100.
Figure 3 compares the estimation errors of different algorithms as a function of sample size, with the signal strength set to 5 in (a) and 10 in (b). Our proposed algorithm consistently achieves smaller estimation errors than the other algorithms,, except when = 100, where the limited number of samples leads to large estimation errors across all methods. When the signal strength À is increased to 10 in Figure 3(b), similar conclusions can be drawn.
Table demonstrates the algorithm's robustness to varying sparsity levels k. The results show that increased sparsity leads to mild increase in iteration count while maintaining nearly constant per-iteration computation times. This translates to only slight growth in total running time with increasing k. indicating that our method efficiently handles varying degrees of problem sparsity without significant computational overhead.
8

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
94
8375
6175
0.335
0.5
0.35
9.35
0:15
01
01
011
0.755
2605
2.00
0425
0.23
0.025
3.005
.0025
(a) DT: Error 0.6283,
(b) Soft-CT: Error 0.2248,
(c) Hard-CT: Error = 0.3395,
(d) Proposed: Error 0.1907, F-score = 0.8325
F-score = 0.4804
F-score= 0.8038
F-score 0.7778
Figure Comparison of different algorithms in the "three-peak" experiment. The black curve represents the spike, while the red points indicate the estimated spike by each algorithm. Our algorithm achieves a smaller estimation error and a higher F-score compared to ofth-of-the-art algorithms, demonstrating superior performance in both estimation accuracy and support recovery.
5. Conclusions
Sparsity
100
125
150
175
200
Number iterations
8.55
9.24
9.85
10.14
10.39
In this paper, we have investigated the single-spiked covariance model for sparse PCA, developing two complementary algorithms. Our thresholding- algorithm operates polynomial time and successfully bridges the gap between information theoorttitt limits and computational efficiency, requiring only Ω2(k: log p) samples when À = [[[||| Building upon this foundation, our two- stage nonconvex algorithm combines the thresholding technique with truncated power iteration to achieve the minimax optimal rate of statistical error under the same sample complexity. Extensive experiments have validated the superior performance of our algorithms terms of estimation accuracy and computational efficiency.
Time for an iteration(s)
0.022
0.024
0.023
0.023
0.022
Total running time(s)
0.189
0.219
0.224
0.234
0.231
Table Impact of sparsity k computational time in the refinement stage of Algorithm 2, with p 5000 and n 3000.
The empirical results across Tables 1, 2, and 3 provide strong validation of our theoretical complexity analysis. The observed near-linear scaling of per-iteration time with both p and n aligns precisely with the thcoretical (pn) complexity bound. This correspondence between theoretical predictions and empirical measurements, coupled with the favorable iteration scaling properties, establishes our method as both theoretically sound and practically efficient for scale sparse principal component analysis.
While this work makes significant progress in addressing the estimation problem for sparse PCA, several promising directions merit further investigation. First, we have not addressed the support recovery establishing support recovery consistency the spike with optimal sample complexity-which poses distinct statistical and computational challenges. Second, evaluating the optimality of our signal strength condition and exploring whether this requirement can be relaxed would deepen our understanding of the fundamental limits in sparse PCA.
4.4. Experiments on "Three-Peake Peakk Data
The "three-peat peakkeeppaa"aa experiment (Johnstone & Lu, 2009) is widely recognized benchmark for assessing the performance of sparse PCA algorithms. This experiment aims to assess the effectiveness of an algorithm recovering a structured sparse signal from noisy data. In this setup, the true spike vector v is constructed as a mixture of three Beta densities on [0, 1). producing characteristic shape with three distinct peaks, as illustrated by the black curve in Figure 4. For our experiments, we set the dimension p 2000, the sample size to n = 1000, and the signal strength to X = 8. The choice of these parameters reflect a high-dimensional setting where the number of variables significantly exceeds the number of observations, posing substantial challenge for accurate signal recovery.
A Acknowledgement
This work was supported by the Hong Kong Research Grant Council GRFs 16306124, and 16307023, and Hong Kong RGC Postdoctoral Fellowship Scheme of Project No. PDFS2425-6S05. We would also like thank the anonymous reviewers for their valuable feedback on the manuscript.
Our proposed algorithm demonstrates superior performance compared to state-of-the-art methods. Specifically, it achieves smaller estimation error and higher F-score, reflecting superior accuracy in both estimating the true spike and correctly identifying the support of v. These results highlight the advantages of our algorithm terms of estimation accuracy and support recovery.
Impact Statement
This paper aims advance the theoretical and compu- tational aspects of sparse PCA, enabling more efficient high-dimensional data analysis across diverse application domains. We do not foresee any specific ethical concerns arising from this work.

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
References
Horn, R. A. and Johnson, C. R. Matrix analysis. Cambridge University Press, 2012.
Amini, A. A. and Wainwright, M. J. High-dimensional analysis of semidefiitte relaxations for sparse principal components. The Annals of Statistics, 37752)12877-2921, 2009.
Jagatap, G. and Hegde, C. Sample-efficient algorithms for recovering structured signals from magnitude-only measurements. IEEE Transactions on Information Theory, 65(7):4434 4456, 2019.
Baraniuk, R., Davenport, M., Vore, R., and Wakin, M. A simple proof the restricted isometry property for random matrices. Constructive Approximation, 28:253- 263, 2008.
Johnstone, I. M. On the distribution the largest eigenvalue in principal components analysis. The Annals of Statistics, 29(2):295-327, 2001.
Berthet, Q. and Rigollet, P. Optimal detection of sparse principal components in high dimension. The Annals of Statistics, 41(4):1780 1815, 2013.
Johnstone, I. M. and Lu, A. Y. On consistency and sparsity for principal components analysis in high dimensions. Journal of the American Statistical Association, 104(486): 682-693, 2009.
Brennan, M., Bresler, G., and Huleihel, W. Reducibility and computational lower bounds for problems with planted sparse structure. In Conference On Learning Theory, pp- 48-166. PMLR, 2018.
Krauthgamer, R., Nadler, B., and Vilenchik, D. Do semidefinite relaxations solve sparse PCA up to the information limit? The Annals of Statistics, 43(3):1300- 1322, 2015.
Cai, J.-F., Li, Lu, X., and You, J. Sparse signal recovery from phaseless measurements via hard thresholding pursuit. Applied and Computational Harmonic Analysis, 56:367-390, 2022.
Lasserre, J. B. An introduction to polynomial and semi- algebraic optimization, volume 52. Cambridge University Press, 2015.
Cai, J.-F., Li, J., and You, Provable sample-efficient sparse phase retrieval initialized by truncated power method. Inverse Problems, 39(7):075008, 2023.
Laurent, B. and Massart, P. Adaptive estimation of a quadratic functional by model selection. The Annals of Statistics, 28(5): 2000.
Cai. J.-F., Long, Y.. Wen, R., and Ying. J. fast and provable algorithm for sparse phase retrieval. In The Twelfth International Conference on Learning Representations, 2024.
Ma, T. and Wigderson, A. Sum-of-squares lower bounds for sparse PCA. Advances in Neural Information Processing Systems, 28, 2015.
Candes, E., Romberg, J., and Tao, T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on Information Theory, 52(2)489-550, 2006.
Novikov, G. Sparse PCA beyond covariance thresholding. In Proceedings of Thirty Sixth Conference on Learning Theory, volume 195, pp. 4737-4776, 2023.
Shen, H. and Huang, Z. Sparse principal component analysis via regularized low rank matrix approximation. Journal of Multivariate Analysis, 99(6)1005-11044 2008.
d' Aspremont, A., Ghaoui, L., Jordan, M., and Lanckriet, G. A direct formulation for sparse PCA using semidefinite programming. Advances in neural information processing systems, 17, 2004.
Vershynin, R.
High-dimensional probability: An introduction with applications in data science, volume 47.
Deshpande, Y. and Montanari, A. Sparse PCA via covariance thresholding. Journal of Machine Learning Research, 17(141):144,,211
Cambridge University Press, 2018
Vu, V. Q. and Lei, J. Minimax sparse principal subspace estimation in high dimensions. The Annals of Statistics, 4 (6)) 2013.
Ding, Y., Kumisky, D., Wein, A. S., and Bandeira, A. S. Subexponential-time algorithms for sparse PCA. Foundations of Computational Mathematics, 24(3):865- 914, 2024.
Vu, V. Q. Cho, J., Lei, J., and Rohe, K. Fantope projection and selection: A near-optimal convex relaxation of sparse PCA. In Advances Neural Information Processing Systems, volume 26, 2013.
Donoho, D. L. Compressed sensing. IEEE Transactions on Information Theory, 52(4):1289-133 1306, 2006.
Wang. G., Zhang, L., Giannakis, G. B., Akçakaya, M., and Chen, J. Sparse phase retrieval via truncated amplitude flow. IEEE Transactions on Signal Processing, 66(2): 479-491. 2017.
Gao, C., Ma, Z., and Zhou, H. H. Sparse CCA: Adaptive estimation and computational barriers. The Annals of Statistics, 45(5):2074-210) 2017.
10

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
Wang, T., Berthet, Q., and Samworth, R. J. Statistical and computational trade-offs in estimation of sparse principal components. The Annals of Statistics, 1465);189-1130,, 2016.
Wang, Z., Lu, H., and Liu, H. Tighten after relax: Minimax- optimal sparse PCA in polynomial time. In Advances in Neural Information Processing Systems, volume 27, 2014.
Wu, F. and Rebeschini, P. Hadamard wirtinger flow for sparse phase retrieval. In Proceedings of The 24th International Conference Artificial Intelligence and Statistics, volume 130, pp. 982-990, 2021.
Yuan, X.-I. and Zhang, T. Truncated power method for sparse eigenvalue problems. Journal of Machine Learning Research, 14(4), 2013.
Zou, H., Hastie, T., and Tibshirani, R. Sparse principal component analysis. Journal of Computational and Graphical Statistics, 15(2):265-266, 2006.

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
A. Proofs
We present the proofs for Proposition 3.1, Theorem 3.2 for our thresholding algorithm, and Theorem 3.3 for our two-stage algorithm. First, we introduce some technical lemmas Appendix A.I. Subsequently. we present the proofs for Proposition 3.1, Theorem 3.2, and Theorem 3.3 in Appendices A.2, A.3, and A.4, respectively.
Throughout Appendix A, we define matrix E : Σ - E Σ and smallest 8-sparse eigenvalue by
For any symmetric matrix A E R=xn, we define the largest
max mmax (A,s) wER" 1|t0 1.]/wllo=s
w T
Au,
min w T Aw, Amin (A, weR.. _|/0/|22 1//==l//ooll=== =s
respectively. Also, we define the maximum spectral norm of all x submatrices of A by
P(A.s) s) max Amaax (A 8)l Dminn (A, s)(}
(12)
A.1. Technical Lemmas
In this section, we introduce and prove some technical lemmas that will be used for the proofs of Theorem 3.2 and Theorem 3.3.
The first lemma bounds the quantity p(E,s) defined in (12).
Lemma A.1. For any E (0,
"ep"
22
((p(E,s) & 3tA] 2 1 1
exp
(A+ (>+
(13)
for some absolute constant C₁ > O.
Proof. Denote the set of s-sparse vectors in R2 by Tg == ( w1l|l||m2 = 1, w|l I s}. For any 8 E (0, 1), there exists a
set set No8 C Tp such that for any w Tg, there exists wy E Ns such that supp(w) supp(ws) and w - -wslp ||2 S 8 and [Ws] < (*) < (22)) (Baraniuk et al., 2008).
From the definition,
p(E,s)
max M.ETY
y/ Ez =: =. Ez..
supp())=<upp(s)
From the definition of No. there exists ys: 28 Nss such that supp(ys) supp(y-) = supp(z.) supp(zs), 11 ys la I 8 and ||=- zal12 < 5. Then we have
y.= Ez. = y. E(z. - 28) + (y- - y0)² Ezs + y Ezs
S 225y Ez. + 1/5 Ezá,
which implies that
p(E, S 25)
25) 22))-1y11 -y Ezs $ (1 - (5)
max
y" Ez, E'Ez.
(14)
g.sENs
supp(y) supp(=)
where the inequalities hold when - 28 > 0.
Now for any (y. 22) E Ng. we bound 1y Ez| as follows. Recall that
n Ez Ez
3 [ T - E [x
=:2.
where x; = VAg.u + &.. With the assumption of 9; and &,. we have y x; id N(0, + 1) and कद ud N(0,A(=T x(=T))2++) hence y x z is sub-exponential (Vershynin, 2018). Using the Centering Lemma for sub-exponential random variables (Vershynin, 2018), we have
1 v,, S Say Myy o)2 +1 +c
12

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
where 61 0 is some absolute constant and the last inequality holds since Bernstein's ss inequality (Vershynin, 2018), for any t (0, 1), it holds that
l
|
||0||2 1. Therefore, by
x²
(1 Ez 2 s 2expl-Cin
(15)
(A + 1)2
where Cj some absolute constant. By taking union bounds for all (y, 2) E No, we obtain that, with probability 12x2
exceeding 1 - 2( ) expp-CCin(4tx3)
yy Ez s 1X.
M.N.N supp()))supp(s
Setting 5 together with (14) leads to (13).
The next lemma gives the concentration inequality for x2 -qqure variables.
Lemma A.2 (LLaurent & Massart, 2000)). x(m) denotes central chi-squared variable with m degrees freedom. For all
t 12 0,
Plx(m m 2 2/mm + 26] S exp exp(-t). (-t), P(x(m) S -2 mt] exp(-t))
(16)
(17)
The next lemma bounds the error between and s-sparse largest eigenvector of Σ in unit length.
Lemma A.3. Let A [p] be such that A AUS+ + and A| = s. Let be an eigenvector of unit length corresponding the largest eigenvalue of Ex. pl E, 8) 1 then have
UU 1|
dist(u, 22
|ual 1 +
(²(Ess)
(Allea -12
Proof. Recall that E = E -ES. (Horn & Johnson, 2012), we have
Denote the largest eigenvalue of Es, i.c. = =>x(Ex)) From Weyl's inequality
2
A E ΣA
(A))
2 2l +1 -((Es),
(18)
where the last incquality holds since for all j 2 2,
E
Σ uook and A, (EA) 2 Es) from the definition. Similarly, we have
(F(EE()) My(EN) s /// -
+\\((E))
=1 max (A,) $14 p(E,)). (19)
Notice ||u||2 = but 1|a I 1. We divide u as
U
u u=
+ a2
|uA
with uTy= = 0. 11 and af + a² = 1. Then we have supp(y) C A, and
Σ
Aaj |u lz
+ Jazy = du = Σx = ai |ull2z
+ az Eyy.
By taking the inner product with y, we obtain
8" S
Xag
a2y" TEA.
| 2
13

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
Since VΛ is the eigenvector of EE [A] and o 0. we have y E
0. This leads
y"(E), + E (
uyy y" E
a2l = |a
PAlla
= =Iaıl
PAlla 12
i y- EAy
Ā y Σ Ayy
Since supp(y) C A, we have
y² EA Tual2
S p( (). Moreover, since y is perpendicular to dy. from (19) we have
y" Ei
I max {5(EN) ((E,).
j22
So from (18) and s)E < 2 1 we have
Wi y E 1
[2
p( E.ss)
a1
X 4/ 3T EN
u/ 12 -2p(E,
Then, since a$ + a² 1, we have
af 2 14
p² R.s)
(Allea (.)().)))
which implies that
dist(u, o^)² min {llu^ u1l2 |/ + 2 22
=I 2 + 2 2/a1l: | l
uu l|2
< 11al12 + 1 2
(R,5) (Allex 2p((
+
The following two lemmas are used to prove the convergence of truncated power method in our proposed two-stage algorithm.
Lemma A.4 (CYuan & Zhang, 2013)). Let z be the eigenvector with the largest eigenvalue (in absolute value) of a symmetric matrix A, and let < 1 be the ratio of the second to the largest eigenvalue in absolute values. Given any y such that ||yl| 1, let yy Ay/ Ayilz, then
T y
2
y/
(1+2(1- n²) (1-1
Lemma A.5 (CYuan & Zhang, 2013)). Consider y with 11y/1o s. Consider z and let J supp(z,s') s') be the s' indices with the largest absolute values in If ||y||2 || 1, then
1 T z - Vs/s' min {vi {
1 プイ
A.2. Proof of Proposition 3.1
Proof of Proposition 3.1. For any E S and any j 8" using (4) and (5) obtains:
1
Σ
وز - | و1 زز
2 Au3 +
2 Dd.
1
Egyl
[
+
[B{ j/gt
Eyye
$ $1
2 7d
Following from the fact that ga minjes v7, one has Σ
2
E.yrr.
14

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
A.3. Proof of Theorem 3.2
We organize the proof of Theorem 3.2 in three steps. First, we prove that the index jo chosen in Algorithm satisfies |%% 2 I| les with high probability. Second, we show that & chosen in Algorithm contains the indices of most of larger non-zero entries of v. Finally, we put everything together.
2
Step 1: Estimating /.ol Recall that jo arg max Σz 15)5n
Lemma A.6. Assume that 2 C 1
and 2
Cyk logp for some absolute constants C22 Cj3 > 0. Then, with
probability exceeding 1 - p- 1300 2 2
Proof. Recall that
1
Viguu + (V)))()) + E.,).
Then we have B [[E}a] Firstly, we consider Ej.... where j. satisfies |/,_1 = 1101loc. Then (17) in Lemma A.2 implies that for any > 0,
Xei + 1. Since Vāgeu + Eii lid (O A03 + 1), obtain n nLjs/ (Aoī +1) lid
x2my-
=X1lu &2e P S. -1- +1} F F x(m) n -1 21 13e +1
(20)
s expl- n 221210011
4 (Allo) + 1)
Secondly, we consider S = {i (0 <
|| }. By using (16) in Lemma A.2 and taking union bound, we have,
2
Enj jjj 2 (1 - e)ll +1}
P
max jEE
(1 - =| +
< (p - 1)IP xim) 2
1 +1
) |
72
7 All &l +1)²
(21)
I (p - exxp
1 2c ) 1il
(1++) +
x +1
( 2) =)>
s (p - 1) exp!
x + 1)(1 + (1 (1-E) (1-E)\||e|33
where the second inequality holds from the definition of S..
Now we combine (20)(21) and set = 1. The complementary events in (20)(11) lead to maxjes, EnS < fN1lo113 + +1< 1< S S jo.joo which implies that jo & S1, i.e. |0jol N || Therefore, we obtain
|10|
P
2
2 p pexp
(22)
2
16 v 2 1)2
Together with the assumptions of À and n, (22) leads to the desired result.
Remark A.7. From (22), a sufficient condition for constants C₁2 and C22 in Lemma A.6 is
C22
4
/lllalle +1
which can be simplified as
C2
C Ch +1
2 4V2.
15

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
Step 2: Estimating |1ll2. For any E (0, 1], we define S := {i { layl
2vb and S = S Sy. Then we have
|1s; 112 1k k and | 11 2 1
7- Since |/0 2 vie Lemma A.6 implies that 1//0 2 2vT 2
1 2 7 with
high probability, and thus jo E S. The following lemma shows that S C S with high probability, where is chosen in Algorithm 1.
Lemma A.8. For any E (0. 13 if 2 Gillellad and n 2 Cgz log p for some absolute constants C₁. C2 > 0, then, with probability exceeding - 3p 1 S C S.
Proof. It suffices to show that with high probability,
min
$.10
max jeS
JEST
To prove this, first, show that for any E S2, where S2 = {i E 81108l1 (10s/2 2
min jES.
E
max jeS
/
which needs bound E and E51 E
[ for all j E 8 and E S
For any E S2, we consider maxjes"
From the definition, for all j E 8
ΣS Σ
Ens (V Ew.).
hence we have EE [ 0. Since Eid iid N( 1) and VAgun + E.p ind N(O, Xuf 1). E.,(Vggm Ex,) is sub- exponential (Vershynin, 2018) with
(VAg.vn + &.))lv, S <caAAle +1.
where C3 >> 0 some absolute constant. Then, by taking union bound and using Bernstein's inequality (Vershynin. 2018).
for any E1 0, have
max PP jES=
Eye 2 E1
}
s 2(p k) exp[-cammi min
E1
(23)
Rl|l|k & +1'
// +
where 0₁ > 0 is some absolute constant.
For any / E S2, we consider max
Li.jo - E EE[ssoo]
From the definition, for all j E St.
Σn
1 Σ
(Vgee + (vg) Ec.).
1=1
which implies that E [[,,]] auyer for all j E S.. and thus (0, Au3 + 1), (Vigez + (V..)(Cggg + E.,) sub-exponential (Vershynin, 2018). Using the Centering Lemma for sub-exponential random variables (Vershynin, 2018), have
E Σ
On the other hand, since Vāg.uy+ iid
4yk
[I(VAgeu (V.))))) E..) - -E <Allolll 1),
where cs 0 an absolute constant. Then, from Bernstein's inequality (Vershynin, 2018), have
2
P max
Ez1r E ES
2 E2 < 2k explplccnmiin
(24)
jEST
[A
112 All
12
16

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
where Cg > 0 is an absolute constant.
Now we combine (23)(24) and set E1 =2
>1//e]sss, 8yk
The complementary event in (23) is
max JCS
Σ 3
8vk
Moreover, the complementary event in (24) leads to
/
>111o11o
--ppl>==>a|llll >>llu1lo
|
Σ E5
E
Σ
$50
vij E est.
4vh
8vk
8v
These two inequalities implies that min jEST
max jES
for any E S2.
Finally, by taking union bound and using (22)
P(s c
c8] =]=
=
P
min jEST
> max jES
Eg
:00 = 1 1
SESg
2 E
P
min
Σ 31
max jE8c
|
P [jo + 17)
jES
1652
|5,51 }
2 Σ
(P [jo- = 1]
P
min jEST
Σz11
S
max jes
(25)
SES
nd² 11
2 1 pexp!
16(A|| & + 1)
71/11111
2p expl-c7nm min
644 + 1)²' 8vl((1|001 +1)
where e7 = min (.), Together with the assumption of À and and since E (0, 1]. (25) leads to the desired result.
ca).
Remark A.9. From (25), sufficient condition for constants C and C33 in Lemma A.6 is
C2
4V2
Ca
+1
VkI
1 00
3
Cl bo
+1
8
V
C₁
24
C
Clllu
loc +1
v
can be simplified as
C₁
v6
1
CCH
2 4./2 max
C7 Cl |oc +1
V07
or
C₁ VCH
2 max max
C +1
Step 3: Putting everything together. Now we estimate dist (v",, v) and prove Theorem 3.2.
Proof of Theorem 3.2. For any 7 E (0, 1J. we assume À > C |||| and n 2 -2 loo p with some absolute constants
C1, C2 > 0. For simplicity, we denote p = p(E, k).. Under the above assumptions, by applying Lemma A.1 with
t 16 7, 8 = 1: and Lemma A.8, we have
3 7A, ju E E cs) 21 1 2p- -3p 3p--1 2 1 5p
PP p=
(26)
16
17

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
Under the event in (26), we estimate dist(v", v). Since supp(b") = S, we have
dist(v', v)² dist(v", vz)² + 1
(27)
Firstly, we estimate | Since S C S,) STUSU US., we have
2 < 4
1
2 11ug12 >
3
1 11 - ||s. 1 + ||ss
4
Secondly, we estimate dist vz)². Applying Lemma A.3 with A S and s k; we obtain
110s 112
dist(v", .u))² & |1 +
S 110 2
p²
+ (2)-20)2
0²
Allva113-20
20)2
where the last inequality holds since ||5ll 4. Therefore, by using Lemma A.3 and || l1 S T. we have
7²
2
2
dist(v", ,0g S max
2 -
2
4
V1
1
2p)2
(
2p)
1
1 1-
T
S max
2
2
4
1+
p²
1 (12 2p)2
P
(1)-2p)2
2
(2
2p)² + (2 1 ))
2p²
4
max
(1) 2p) + p²
- - 2p)² + p2
22
2p²
+
4
(+) 2p)² + p²
follows from (27) and p S 7 that
dist(o"; v) < 22 72 +
72
2 2
completing the proof.
Remark A.10. From (13)(55,, similar to Remark A.9, a sufficient condition for constants C₁ and C.2 Theorem 3.2 is
C
2 + 1ogg(9e)
VC2
2 /
Gi ||o|lo
Ci
C
VC2 1illoloo
2 4. max
+1
Vc
Vk|lullo
which can be simplified as
C
2 + (0g2(9e)
vb
VC2 C +1
max
16,
4v max
1
C/
V07
18

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
A.4. Proof of Theorem 3.3
We organize the proof of Theorem 3.3 into two parts. First, we show that v" falls into small constant neighborhood of Subsequently, we prove the convergence of the truncated power method.
Proof of Theorem 3.3. Throughout the proof, we assume À 2 Cs|lellad and 2 C.k log p with some absolute constants
C3, C4 > 0. We denote k- = k 2k', p = p E 1)) and F, = supp(e"), where k'' = Cşk for some absolute constant Cs 2 1. Setting t = V klog p and 8 = 1 in Lemma A. and 7 = 1 in Lemma A.8 for some constant 0 < Ç < 0.001 /Ca, together
with the above assumptions, similar to Theorem 3.2, with the probability exceeding - 2p 3p 2 1 - 5p 1 the following event holds:
k klogp X, dist(v v)) - 1 1
E = p ps
We will continue the proof under this event.
Step 1: Estimating |o"'o) Since 1 2 dist(v, 02)2 = - 2 -2| 0%/. we have /0"00'
D.|.
00 2 0.5.
Step 2: Convergence of truncated power method. To prove (11), we will first show that dist(v, v²) S 1 by induction.
We denote A, =
U.S then |Ael S k + 2k' = ki Also, we define
u E x u' /1 0' 1|2,
(28)
hence we have u²= = u/z// /125/ 122. Let be the ratio of the second largest (in absolute value) to the largest eigenvalue of EA,- Then, since S C Ar, similar to (18)(19), we obtain
{g(SKrn)
maxj*1
1 1+1 p
1 + 0.001A
||0|| 0o + 0.001C; B1,
A1 (Ex,)
> Al|UA, 1 +
+ 0.999X
| 8o + 0.999C3
where in the second inequality we use p S SV Flgg X < 0.001 and in the last inequality we use À 2 C31|v11=1 Obviously
have < B1 <! 1.
Let T be unit eigenvector corresponding to the largest eigenvalue of En, and satisfying vZ 2 0. hence we have dist(v, ū) = |/ ج||ا. Then, using (28) and Lemma A.4, we have
oTu'l o 0 (1+ -1- (1x3)) (1-1010000
which implies that
1 (1-1-
B²
1 - 0² un S (1 (1- 100 0--0-111
11"0 -120 | 101)))111))
(29)
2
Since S C Ar, Lemma A.3 gives
p²
log 2
| 0112 = dist(v, (0) 22
(30)
v 1 + (A 20)
(A - 2p) 0.9982
n
9982
where in the second inequality use 1 vI+a for 0, and in the last two inequalities we use p S CV k X S 0.001A and the conditions of and Ç. Note that the induction assumption dist(v, , w111) s implies that |0" o which with (30) further leads to
n
> 0.5.
nJ n 1
o
w
(v o 2 |00 o - 1/00 01|2||0* 1112 2 0.4989.
(31)
Plugging (31) into (29), we have
ग. uz (0.6261 0.374Bf) (1- o 1
19

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
which is equivalent
di u²) S V 00.6261 + 0.374B3 dist(o, v²=1), where we use ||ō||2 ||u2|22 | 1l2 = 1. For unit vectors v, v' U, obtain
(32)
dist(ü, 1) I dist(U, v) + dist(v"" v).
(33)
This is because
dist(ö, v) + dist(v' ,0) | 2112 + 11e r20 112
| + 1200~1122 2 dist(ö, o² )
where 71,72 E {#1} and we use (10). Similarly, for unit vectors v, u 6, it holds that dist(v, u') I dist(v, 8) + dist(u', ö).
(34)
Using (30), (32), (33), and (34), we have dist(v, u')
C (1+ + 0.6261 + 0.374B3)
klogg
(35)
$ 0.6261 + 0.374B} dist(v,
dist(v,v'11)) 1(=1)+
0.998
Since k Cyki Lemma A.5 generates
v T J |A
2 | u''l - C5 1/2 min {V1 1//1-1172 10T (1+ 1
1/21 (1
(1-11-
2 v='u" u²! - C5 1/2 (1 + 07 (1
which implies that
1 100 [u')) Fr S 1 - u' + C,5 1/2 1/21/1111 + 17.1/3 (1- 11- B (o) vF'2' u u
where B2 V 1 + 20, 1/2 12(11 1/3). Recall that rt u /t///luzz 112. Then have
dist(v, (2,0') v²)
2 2 2 0'1 /2 2- 2--2/0 [=])= 1
2 2 2 [u')se| v B2 2 - (1)
B₂ dist(v, v. u²)
(36)
CB₂
klog
< B2 / 0.6261 + 0.374 dist(v, )
1 + V 0.6261 + 0.37433
0.998
B₂
< 0.660.6261 + 0.3743 di v'-1) + 998
+ 0.6261 + 0.374B²
where in the last second inequality use (35), and the last incquality holds from the conditions of n and Ç. Since dist( -'-1)= < 1, the above inequality also implies that dist(e, w²) < 1 with suitable constants C3. Cs (constants B₁. Bz).
Therefore, we complete the induction, which proves that dist(o, v²) < 1 for all t. As a result, the above inequality holds for all t, which leads to
k klog
dist(u, 0²) = d. dist(v, w--1)= + Bay
klogg
k loog
< d2 dist(v, v 22 + d. BBy
+ Bs V
n
B3
k klogp
< d' .diistsv, 00) + 1- d
20

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
where := B2 0.6261 + 0.3744f and B3 == inequality is the desired result.
CB2 (1
V..6261
Bs Denoting d' == I-d
the above
0.998
Remark A.II. From (13), similar to Remark A.I0, sufficient condition for constants C3. C4, Co in Theorem 3.3
0.001
<
C3
4 (3e)
2
VCB C3
+
+||u|l00
3
30
C3
v6
V0+
4v max
C3 C3llu
Do +1
Ve7
V
999B, V 0.6261 + 0.374B3 + Bg S 998,
where B,
| +0 Tol. +0.9990a and B2
C 1/2). (1 / 1/ +
can be simplified as
<
0.0001 Va,
( C3
1
4 + 3 1ogz
VOB C3 +1
3/
3C"
C3
VG
VC
4v2 max
C3 +1
VoT
999B2 V 0.6261 + 0.374B
+ Bg S 998,
where B,
1+1,0001Ca
1+0.99003
21

---

Fast and Provable Algorithms for Sparse PCA with Improved Sample Complexity
B. Additional Experiments
To further validate and analyze the behavior of our proposed Algorithm 2, we conduct two complementary studies. First, we assess the robustness of the algorithm under suboptimal signal-strength conditions by varying the parameter using the same data-generation procedure as Section 4.2. Second, isolate the contribution of the refinement stage by comparing the estimation accuracy after the initialization step alone with that of the full two-stage pipeline.
B.1. Robustness of Algorithm
To assess the robustness of Algorithm when theoretical assumptions are relaxed, we conduct experiments across range of signal strength values A. Following the experimental protocol from Section 4.2, we generate sparse eigenvectors with dimension 1000 and sparsity level k = 20. All reported results represent averages over 200 independent trials.
Table 4. Estimation error verus sample size n under different À
Sample size n
100
200
300
400
500
600
700
800
900
1000
À = 0.5
1.4003
1.3971
1.3973
1.3996
1.3976
1,3903
1.3939
1.3943
1.3927
1,3897
À = 2.5
1.3174
1.1650
0.7896
0.5427
0.3980
0.2459
0.1736
0.1535
0.1331
0,1280
À
0.7188
0,2635
0.1347
0.1062
0.0943
0.0863
0.0789
0.0751
0.0701
0.0664
À 7.5
0.2762
0.1246
0.0961
0.0843
0.0749
0.0683
0.0626
0.0596
0.0556
0.0526
Table 4 presents the estimation error function of sample size for varying signal strengths. The results show that larger values yield substantially lower estimation errors, and that the advantage of strong-signal settings over weak-signal settings becomes increasingly pronounced as n grows. Crucially, even under weak signals (small X), the algorithm degrades gracefuly-mmainaaiiing robust accuracy rather than failing catastrophically. These findings underscore the practical resilience of our method across a broader range of scenarios than our theoretical guarantees alone would cover.
B.2. Improvement of Refinement Stage in Algorithm 2
To quantify the contribution of the refinement stage in Algorithm 2, we compare the estimation accuracy after the initialization stage alone versus the complete two-stage procedure. We maintain the experimental protocol from Section 4.2 for generating v, setting dimension = 1000 and sparsity level & = 20. Results are averaged over 200 independent trials.
Table Estimation error verus sample size after different stages Algorithm 2
Sample size
100
200
300
400
500
600
700
800
900
1000
Initialization
1.0553
0.6733
0,4183
0.2880
0.1981
0.1401
0.1028
0.0848
0.0750
0.0697
Refinement 0.7188
0.2635
0.1.47
0.1062
0.0943
0.0863
0.0789
0.07.5
0.0701
0.0664
Table 5 reveals the substantial improvement achieved by the refinement stage. The truncated power iteration consistently reduces estimation error across all sample sizes, with the most dramatic gains occurring at smaller sample sizes for instance, reducing error by over 60% when 200. As the sample size increases, both stages converge toward similar performance, yet the refinement stage maintains a consistent advantage. These empirical findings validate our two-stage design and demonstrate that combining careful initialization with iterative refinement is essential for achieving optimal statistical performance.
22


Reviews:
[{"review_id": "tvKMtgZ3Jv", "reviewer": "Reviewer_LyQJ", "summary": "This paper looks at sparse PCA with few samples in the spiked Gaussian\nmodel, under the assumption that the largest single coordinate of the\nspike has pretty high variance.\n\n## update after rebuttal\n\nI remain unenthusiastic about this paper, but I don't strongly object to it.  The main contribution is the proposed parameterization, in terms of the largest coordinate in the spike; given this formulation, getting an algorithm is fairly straightforward.  I do agree that identifying that this is a fairly simple formulation with a potential for significantly improved sample complexity is a real contribution.  But I'm not excited by *further* restrictions on the spiked covariance model, which is already toy.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 1. Under the assumption given, the result seems like it should be pretty straightforward to achieve. Why not just take $(\\hat{\\Sigma} e _{j _0}) _{\\hat{S}}$ as your estimate $v^0$? The scaling $(\\lambda v_j) > 1$ only helps; ignoring it, we have $v_j < 1$ so the error is $|v^0 _j - v_j| \\leq c \\sqrt{\\log p / m}$.\n\n**Reply:** We thank the reviewer for the interesting suggestion. Our response is organized in three parts.\n\n**1. Slight Modification Required for the Reviewer’s Estimator**\n\nThe reviewer suggested using $ (\\hat{\\Sigma} e _{j _0}) _{\\hat{S}} $ as an estimator for $v$, where $e _{j _0}$ is the canonical basis vector. However, under the spiked covariance model, we have $E[\\hat{\\Sigma}] = \\lambda v v^T + I$, so that $E[\\hat{\\Sigma}] e _{j _0} = \\lambda  v _{j _0} v + e _{j _0}$. Thus, while for $j \\neq j_0$ the entry is $\\lambda v _{j_0} v_j$, the $j_0$-th entry becomes $\\lambda v _{j_0}^2 + 1$, introducing a bias of $1$. To remove this bias, we subtract $e _{j_0}$ and define the modified estimator\n$$\n\\hat{v} _{\\text{new}} = \\frac{(\\hat{\\Sigma} e _{j _0} - e _{j _0}) _{\\hat{S}}}{ ||(\\hat{\\Sigma} e _{j _0} - e _{j _0}) _{\\hat{S}} ||_2},\n$$ \nthereby properly centering the estimator.\n\n**2. Technical Challenges**\n\nBoth the estimator $\\hat{v} _{\\text{new}}$ and our proposed estimator rely on the random index $j_0$ and the estimated support $\\hat{S}$. Because of this, establishing the sample complexity of $\\hat{v} _{\\text{new}}$ is nontrivial and requires extra technical development. \n\nFirst, since $\\hat{v} _{\\text{new}}$ involves $j_0$ and $\\hat{S}$, its analysis requires key results (e.g. Lemmas A.6 and A.8) to control their randomness. In particular, Lemma A.8 (see Lines 880--888) directly leads to the assumption $\\lambda = \\Omega ( ||v|| _{\\infty}^{-1} )$, which ensures the necessary probability bounds in concentration inequalities (as in Eq. (23)).\n\nSecond, since $j_0$ is random, one cannot treat $\\hat{\\Sigma} _{j,j_0}$ as if $j_0$ were fixed; similar techniques as in Lemma A.8 are needed to handle its variability.\n\nThird, since the true vector $v$ is unknown, $\\lambda v _{j _0}$ cannot directly serve as the scaling between $(\\hat{\\Sigma} e _{j _0} - e _{j _0}) _{\\hat{S}}$ and $v$. Instead, one must first accurately estimate the $L_2$ norm $|| (\\hat{\\Sigma} e _{j _0} - e _{j _0}) _{\\hat{S}} ||_2$ \nto properly normalize the estimator.\n\n**3. Comparison: Computational Load and Empirical Performance**\n\nThe two estimators differ mainly in their final computational steps. Our proposed estimator (via Algorithm 1) computes the leading eigenvector of the $k \\times k$ submatrix $\\hat{\\Sigma} _{\\hat{S}}$, which costs $O(nk^2)$, while the modified estimator $\\hat{v} _{\\text{new}}$ only normalizes $(\\hat{\\Sigma} e _{j_0} - e _{j_0}) _{\\hat{S}}$, costing $O(p)$. Although $\\hat{v} _{\\text{new}}$ is more efficient, our experiments show that the estimator via Algorithm 1 achieves lower estimation error. For example, Table 1 displays estimation error versus sample size (with dimension $p = 4000$, sparsity $k = 20$ and signal strength $\\lambda = 10$), and Table 2 reports the corresponding computational times.\n\n**Table 1. Estimation error versus sample size for two estimators, with $p = 4000$, $k = 100$, $\\lambda = 10$.**\n\n| Sample size             | 500   | 1000    | 1500    | 2000    | 2500    | 3000    | 3500    | 4000    |\n|-------------------------------|--------|--------|--------|--------|--------|--------|--------|--------|\n| Estimator via Algorithm 1     | 0.9593 | 0.6111 | 0.4309 | 0.3104 | 0.2274 | 0.1652 | 0.1272 | 0.0895 | \n| Estimator $\\hat{v}_{\\text{new}}$     | 1.0920 | 0.7596 | 0.5534 | 0.4136 | 0.3215 | 0.2578 | 0.2215 | 0.1867 | \n\n**Table 2. Computational time(s) versus sample size for two estimators, with $p = 4000$, $k = 100$, $\\lambda = 10$.**\n\n| Sample size             | 500   | 1000    | 1500    | 2000   | 2500    | 3000    | 3500    | 4000    |\n|-------------------------------|--------|--------|--------|--------|--------|--------|--------|--------|\n| Estimator via Algorithm 1     | 0.0887 | 0.1192 | 0.1479 | 0.1737 | 0.1907 | 0.2145 | 0.2416 | 0.2648 |\n| Estimator $\\hat{v}_{\\text{new}}$     | 0.0854 | 0.1168 | 0.1454 | 0.1707 | 0.1875 | 0.2113 | 0.2383 | 0.2622 |\n\nWe will add further discussion and experimental evaluations of this estimate in the revised manuscript.\n\n> 2. The single-spiked covariance model is already a toy model, so adding another assumption beyond the basics isn't that exciting.\n\n**Reply:** Although the single-spiked covariance model is a toy model, a significant gap still exists between the information-theoretic sample complexity and what polynomial-time algorithms achieve. Under the planted clique conjecture, it is known that no polynomial-time algorithm can recover the spike at the optimal sample complexity without additional assumptions. **Please refer to the response to the first question for Reviewer gGz9 for a detailed discussion on the necessity of this assumption.**", "claims_and_evidence": "The paper shows how to achieve O(k log p) sample complexity, which is in general optimal, under an assumption that the largest single coordinate has fairly large variance.", "methods_and_evaluation_criteria": "Under the assumption given, the result seems like it should be pretty\nstraightforward to achieve.  Why not just take (Sigma^ e_{j0})_{S^} as\nyour estimate v^0?\n\nFor any j != j0, (Sigma^ e_{j0})_{j} is going to be distributed as\n\n  (1/m) sum_i lambda v_j0 v_j g_i^2 + (N(0, 1) * N(0,1))\n\nThe noise term is independent and constant variance, so it sums to\nbasically N(0, O(m)), and will be O(sqrt(m log p)) for all\ncoordinates; the leading term has (sum g_i^2) being (m +/- sqrt(m log\np)) for all coordinates; and by assumption, lambda v_j0 > 1.  So we\nget the value\n\n  v^0_j = (lambda v_j) v_j +/- v_j sqrt(log p / m) +/- sqrt(log p / m)\n\nThe scaling (lambda v_j) > 1 only helps; ignoring it, we have v_j < 1 so the error is\n\n |v^0_j - v_j| <~ sqrt(log p / m)\n\nFor the given m = (1/gamma^2) k log p, this is coordinatewise error at\nmost gamma/sqrt(k).  That Linf error of a sparse vector, after\nrestricting to the largest k entries, gives gamma L2 error.\n\nAm I missing something?", "theoretical_claims": "I didn't check the proofs carefully, but the result seems pretty straightforward so I believe them.", "experimental_designs_or_analyses": "The experiments are just synthetic, presumably because the spiked covariance model is a bit of a toy.", "relation_to_broader_scientific_literature": "There's an annoying gap between what we can achieve computationally and information theoretically for sparse PCA; this paper looks at how to work around that gap by parameterizing in terms of a different value.", "essential_references_not_discussed": "Not that I know of.", "other_strengths_and_weaknesses": "The single-spiked covariance model is already a toy model, so adding another assumption beyond the basics isn't that exciting.", "other_comments_or_suggestions": "see above", "questions_for_authors": "Please address my question about a simpler algorithm above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "RyChReUUU1", "reviewer": "Reviewer_NHCS", "summary": "The paper presents efficient algorithms for solving the sparse Principal Component Analysis (PCA) problem, it is one of fundamental problems in machine learning. The proposed algorithm significantly reduces the required sample complexity compared to previous polynomial-time methods. Under typical sparse conditions, prior polynomial-time algorithms required sample complexity on the order of  $O(k^2)$ , whereas this paper achieves near-optimal complexity of $k\\log p$, aligning closely with theoretical lower bounds. The authors propose a thresholding-based algorithm along with a two-stage nonconvex approach, combining thresholding initialization with truncated power iteration to ensure both theoretical rigor and computational efficiency. Rigorous proofs demonstrate that the proposed methods significantly improve upon existing polynomial-time methods in terms of both runtime and required sample size. Experiments on synthetic data illustrate practical effectiveness, validating the theoretical claims regarding improved scalability and accuracy.\n\n## update after rebuttal\nThank you for the authors’ responses. I maintain my rating.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 1. Was Algorithm 1 used in Section 4?\n\n**Reply:** Thank you for your comment. In Section 4, all experiments are conducted using Algorithm 2. This is because Algorithm 2 is our full two-stage procedure, which first employs Algorithm 1 for initialization and then refines the estimate via truncated power iteration. Since the refinement stage in Algorithm 2 consistently improves upon the initialization provided by Algorithm 1, we report only the results from Algorithm 2 in the main experimental section.\n\n> 2. In Algorithm 1, when selecting the top $k$ elements, how do you handle cases where some of the last elements are equal? Does this affect the theoretical analysis?\n\n**Reply:** Thank you for your question. We address the issue in two parts:\n\n**1. Potential Ties Do Not Affect Our Theoretical Analysis:** \n\nIn our theoretical analysis (e.g., see Lemma A.8), we establish that with high probability there is a strict separation between the values corresponding to the indices in the true support and those outside of it. Specifically, the minimum value among the selected entries is shown to be strictly greater than the maximum value among the unselected ones. This strict inequality guarantees correct support recovery and renders any ambiguities in the ordering of lower-ranked elements irrelevant. In other words, even if ties were to occur among the lower-ranked entries, they would not compromise the validity of our theoretical guarantees.\n\n**2. The Probability That Any Two Entries Are Exactly Equal Is Zero:** \n\nThe values used in our selection process are computed from continuous random variables (i.e., the entries in the empirical covariance matrix derived from i.i.d. samples). In any continuous distribution, the probability that any two independently drawn real numbers are exactly equal is zero, so the occurrence of two identical values is of zero probability.\n\nIn summary, our support recovery analysis is robust: it relies on a high-probability strict separation that ensures correct selection, and the continuity of the underlying distributions guarantees that exact ties effectively never occur. We will add this discussion to the revised paper.\n\n> 3. What is the stopping criterion for Algorithm 2 in Tables 1-3 in Appendix B? Additionally, how many samples were used for warm-starting (i.e., running the initialization stage)?\n\n**Reply:** Thank you for your comment. First, the stopping criterion for Algorithm 2 in Tables 1-3 in Appendix B is that the relative error between $v^{t-1}$ and $v^t$ is less than $10^{-8}$. Second, Algorithm 2 employs the same set of samples for both the initialization and refinement stages. Hence, the samples used in the warm-start (initialization stage) are exactly the same as those used in the entire experiment. For example, in the experiments reported in Table 1 of Appendix B, the total sample size is 2500; therefore, 2500 samples are used to warm-start the algorithm.\n\n>  4. Conducting empirical studies or discussing how robust the method is when this condition is not fully met could significantly enrich the paper.\n\n**Reply:** Thank you for your comment. To assess the robustness of our method when the ideal condition is not fully met, we conducted experiments with various values of $ \\lambda $. The results below (with $p = 1000$ and $k = 20$) report the estimation error as a function of the sample size:\n\n| Sample size $n$  | 100    | 200    | 300    | 400    | 500    | 600    | 700    | 800    | 900    | 1000   |\n|------------------|--------|--------|--------|--------|--------|--------|--------|--------|--------|--------|\n| $\\lambda = 0.5$  | 1.4003 | 1.3971 | 1.3973 | 1.3996 | 1.3976 | 1.3903 | 1.3939 | 1.3943 | 1.3927 | 1.3897 |\n| $\\lambda = 2.5$  | 1.3174 | 1.1650 | 0.7896 | 0.5427 | 0.3980 | 0.2459 | 0.1736 | 0.1535 | 0.1331 | 0.1280 |\n| $\\lambda = 5$    | 0.7188 | 0.2635 | 0.1347 | 0.1062 | 0.0943 | 0.0863 | 0.0789 | 0.0751 | 0.0701 | 0.0664 |\n| $\\lambda = 7.5$  | 0.2762 | 0.1246 | 0.0961 | 0.0843 | 0.0749 | 0.0683 | 0.0626 | 0.0596 | 0.0556 | 0.0526 |\n\nThese results show that the estimation error decreases with increasing sample size and that the method remains robust across different $ \\lambda $ values, even under less than ideal conditions. We will add this discussion to the revised paper.", "claims_and_evidence": "The claims presented in this submission are generally supported by the provided evidence.", "methods_and_evaluation_criteria": "The proposed methods and evaluation criteria are suitable and aligned with the problem addressed in the paper.", "theoretical_claims": "I took a glimpse at the proofs; they seem to be correct.", "experimental_designs_or_analyses": "The authors utilized synthetic data to examine the trade-offs between solution quality (near-optimality) and computational effort, highlighting the advantages of their proposed algorithms over state-of-the-art methods in achieving a balance between these competing factors.", "relation_to_broader_scientific_literature": "The paper presents a significant advancement in the design and analysis of sparse PCA algorithms by reducing sample complexity while maintaining computational efficiency.", "essential_references_not_discussed": "NA", "other_strengths_and_weaknesses": "Strengths:\n\n1) Handling $\\ell_0$-constrained sparsity PCA is a challenging task. The paper makes a significant contribution by bridging the gap between the optimal information-theoretic sample complexity of Ω(k log p) and the higher complexities, typically Ω(k²), required by earlier polynomial-time approaches. The authors rigorously support their claims with novel theoretical results and numerical experiments—specifically Theorems 3.2 and 3.3—that establish robust performance guarantees under reasonable assumptions.\n\n2) The innovative integration of thresholding-based initialization with refined truncated power iteration represents a creative approach. This two-stage framework not only enhances statistical performance but also preserves computational efficiency, a crucial factor in high-dimensional settings. A theoretical analysis for the combined algorithm is provided. \n\nWeaknesses:\n\n1) The condition on signal strength that $\\lambda = \\Omega(\\||v\\||^{-1}_\\infty)$ is crucial for the proposed improvements in sample complexity. Nonetheless, the paper could be enhanced by a more thorough examination of how practical this assumption is. It would be particularly useful to investigate if this condition is frequently met in real-world scenarios or datasets. Furthermore, conducting empirical studies or discussing how robust the method is when this condition is not fully met could significantly enrich the paper.\n\n2) All experiments for the two proposed algorithms are performed using synthetic data. Although these synthetic experiments are crucial for validating theoretical guarantees, demonstrating the algorithms’ performance on real-world, high-dimensional datasets would greatly enhance the practical relevance of the research.", "other_comments_or_suggestions": "Some numerical experimental results for Algorithm 2 should be moved to the main paper.", "questions_for_authors": "1) In Lines 276-279, the authors wrote, “We present a series of numerical experiments designed to verify the theoretical results and validate the efficiency and effectiveness of our proposed two-stage algorithm.” Was Algorithm 1 used in Section 4?\n\n2) In Algorithm 1, when selecting the top $k$ elements, how do you handle cases where some of the last elements are equal? Does this affect the theoretical analysis?\n\n3) What is the stopping criterion for Algorithm 2 in Tables 1-3 in Appendix B? Additionally, how many samples were used for warm-starting (i.e., running the initialization stage)?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "zGNDGMEvH2", "reviewer": "Reviewer_gGz9", "summary": "This paper presents their algorithms for Sparse Principal Component Analysis (Sparse PCA) under the setting of the Single-Spiked Covariance Model. With the assumption of signal strength, the authors introduce a thresholding-based algorithm with better (big-Omega) sample complexity and show it merges the gap of existing polynomial-time algorithms and information-theoretic lower bounds. Additionally, they propose a two-stage nonconvex optimization algorithm that refines the estimation using truncated power iteration, achieving minimax-optimal statistical error rates. The author provided theoretical analysis for the proposed algorithms and showed the numerical experiments that verified their performance claim in terms of estimation accuracy, sample complexity, and computational cost.\n\n## update after rebuttal\n\nI have read the author's response and would like to thank them for their detailed rebuttal. After reading other reviewers' opinions and the rebuttal, I acknowledge that the proposed method demonstrates a reasonable degree of practical robustness, even though the underlying assumptions may be somewhat idealized. Hence, I have decided to maintain my original score of Weak Accept.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> **1. The strength and necessity of the additional assumption on $\\lambda$**\n\n**Reply:** The single-spiked covariance model is a well-studied framework in high-dimensional statistics. Nevertheless, there remains a substantial gap between the information-theoretic sample complexity, $\\Omega(k \\log p)$, and the best known sample complexity of existing polynomial-time algorithms, $\\Omega(k^2)$. \n\n**To bridge this gap under the planted clique conjecture, it is necessary to introduce additional assumptions**. In fact, prior work—based on reductions from the planted clique conjecture—provides strong evidence that without extra conditions no polynomial-time algorithm can recover the spike with sample complexity $\\Omega(k \\log p)$ [R1-R5].\n\nOur assumption $\\lambda = \\Omega(||v|| _\\infty^{-1})$ ensures that the nonzero entries of the spike $v$ decay sufficiently fast. In many theoretical analyses, $\\lambda$ is taken to be a constant when deriving both the information-theoretic lower bounds and the sample complexities of various algorithms. Moreover, **in models where the nonzero entries of $v$ follow a power-law decay—a scenario well-known in compressed sensing—$\\lambda$ naturally remains of constant order**. \n\nA similar phenomenon was observed in the sparse phase retrieval problem [R6], where signals with power-law decay could be recovered with optimal sample complexity.\n\nWe will expand upon this discussion in the revised manuscript.\n\nReference:\n\n[R1] Statistical and computational trade-offs in estimation of sparse principal components, Annals of Statistics, 2016.\n\n[R2] Do semidefinite relaxations solve sparse PCA up to the information limit, Annals of Statistics, 2015.\n\n[R3] Optimal detection of sparse principal components in high dimension, Annals of Statistics, 2013.\n\n[R4] Sparse CCA: Adaptive estimation and computational barriers, Annals of Statistics, 2017.\n\n[R5] Reducibility and computational lower bounds for problems with planted sparse structure, COLT, 2018.\n\n[R6] Sample-efficient algorithms for recovering structured signals from magnitude-only measurements, IEEE Transactions on Information Theory, 2019.\n\n\n> 2. It assumes i.i.d. Gaussian Noise, is there other cases that this is not satisfied and then your analysis fails? How does it work in more real applications?\n\n**Reply:** Thank you for your question. The single-spiked covariance model is a well-known model in high-dimensional statistic analysis. Its research is mainly from a theoretical perspective. In the single-spiked model, for the sample $x_i = \\sqrt{\\lambda} g_i v + \\xi_i$, the noise $\\xi_i$ is generally assumed to follow a standard Gaussian distribution $N(0,I)$, or sometimes it is assumed to follow a sub-Gaussian distribution with variance proxy $O(1)$ [R7]. Therefore, in our paper, we also assume $\\xi_i \\sim N(0,1)$, which is the same as most of the previous work. \n\nIn our theoretical analysis, the concentration inequalities (Line 607-610 in Lemma A.1, Eq.(19) in Lemma A.6, Eq.(22)(23) in Lemma A.8) hold based on this assumption of $\\xi_i$. If we assume $\\xi_i$ follows a sub-Gaussian distribution with variance proxy $O(1)$ (a generalized Gaussian), we can still derive similar concentration inequalities to complete our analysis. However, if $\\xi_i$ is assumed to follow other kinds of distributions, the same or similar concentration inequalities may not hold, and thus the analysis may fail.\n\nWe will add this discussion in the revised manuscript.\n\nReference: \n\n[R7] Sum-of-squares lower bounds for sparse PCA, Advances in Neural Information Processing Systems, 2015.\n\n\n> 3. What is the robustness of different cases of SNR?\n\n**Reply:** **Please refer to the response to the fourth question for Reviewer NHCS for a detailed discussion on the robustness of our method under various SNR regimes.**\n\n> 4. The truncated power method for sparse eigenvalue problems has already been published. Could the author clarify what is their unique contribution in this part?\n\n**Reply:** While the truncated power method was introduced in [R8] and its convergence shown under a good initialization, that work does not provide a practical procedure for obtaining such an initialization nor does it establish the optimal sample complexity. Our unique contribution is the development of a novel thresholding algorithm that produces an initialization satisfying the necessary conditions for the truncated power iterations. This practical initialization enables our two-stage algorithm to achieve the optimal sample complexity of $\\Omega(k \\log p)$, a significant improvement over the $\\Omega(k^2)$ samples typically required by existing polynomial-time methods. In short, our work bridges the gap between theoretical conditions and practical implementation, thereby enhancing both the statistical and computational efficiency of sparse PCA.\n\nReference:\n\n[R8] Truncated power method for sparse eigenvalue problems, Journal of Machine Learning Research, 2013.", "claims_and_evidence": "The reviewer is not from the community of sparse PCA. The reviewer is convinced by the claims in general but some concerns need to be justified.\n\nFrom intuition, the ‘Single-Spiked Covariance’ may be not suitable for some cases with multiple principle components, which the analysis will not work for sure. It seems quite strong combining the signal strength assumption. Could the author have explained and verified this is a common setting in sparse PCA? \n\nBesides, it assumes i.i.d. Gaussian Noise, is there other cases that this is not satisfied and then your analysis fails? How does it work in more real applications? \n\nWhat is the robustness of different cases of SNR?", "methods_and_evaluation_criteria": "The synthetic experiments are based on their tailored designed setting, this is fine since it's hard to find ground truth for sparse PCA. However, it may be not practical since the assumption is quite strong.\n\nAlso, the reviewer would suggest verifying the proposed method for some real applications that heavily rely on the sample complexity, computational cost, and accuracy.", "theoretical_claims": "Under their assumption, the result looks correct and makes sense to the reviewer.", "experimental_designs_or_analyses": "The synthetic experiments are based on their tailored designed setting, this is fine since it's hard to find ground truth for sparse PCA. However, it may be not practical since the assumption is quite strong.\n\nAlso, the reviewer would suggest verifying the proposed method for some real applications that heavily rely on the sample complexity, computational cost, and accuracy.", "relation_to_broader_scientific_literature": "The reviewer believes this is an interesting work with significant improvement. However, the author needs to justify how general is their algorithm in more practice scenarios in the broader areas.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "The truncated power method for sparse eigenvalue problems has already been published. Could the author clarify what is their unique contribution in this part?", "other_comments_or_suggestions": "N/A", "questions_for_authors": "See my previous discussion.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "1yWyUn0Ax3", "reviewer": "Reviewer_z5wB", "summary": "The paper proposes a two-stage algorithm to obtain the principal component of the single-spiked covariance model (Sparse PCA). The first stage, called the thresholding-based algorithm, obtains a first estimation of the principal component and, most importantly, accomplishes reduced computational cost compared to competing algorithms such as diagonal thresholding. In particular, it reduces the number of data samples required for the computation of the principal component – assuming a restricted strength of the signal of the single-spiked covariance model. To further enhance this estimation, the second stage utilizes the truncated power iteration to refine the solution, accomplishing minmax optimal rate.\n\nUPDATE AFTER REBUTTAL\nThanks for replying to all the questions. Most of them were properly addressed, and I appreciate the effort to clarify both theoretical and empirical aspects. However, I would prefer that the answers which mention future inclusion of additional experiments or material in the final version already present that content within the rebuttal. Promising to include something later is not entirely sufficient, as it leaves the possibility that the addition may not be implemented after acceptance. Including this material now would allow reviewers to properly assess its quality and relevance during the review process. In any case, after reading the rest of the comments and the corresponding response, I would like to update my final score.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 1. You mention that \"operating under a specific condition on $\\lambda$\", which is the implication of this restriction? It would be interesting to include it this discussion in the paper, either theoretically or experimentally.\n\n**Reply:** **Please refer to the response to the first question for Reviewer gGz9 for a detailed discussion on the strength and necessity of this condition.**\n\n\n> 2. Include an experiment showing before and after the refinement using your error metric.\n\n**Reply:** Thank you for your question. Our experimental results clearly demonstrate that the refinement stage significantly enhances the performance of the initial estimator. In the table below, you can see that the estimation error is consistently reduced after applying the refinement across various sample sizes:\n\n| Sample Size **n** | 100    | 200    | 300    | 400    | 500    | 600    | 700    | 800    | 900    | 1000   |\n|-------------------|--------|--------|--------|--------|--------|--------|--------|--------|--------|--------|\n| **Initialization**   | 1.0553 | 0.6733 | 0.4183 | 0.2880 | 0.1981 | 0.1401 | 0.1028 | 0.0848 | 0.0750 | 0.0697 |\n| **Refinement**   | 0.7188 | 0.2635 | 0.1347 | 0.1062 | 0.0943 | 0.0863 | 0.0789 | 0.0751 | 0.0701 | 0.0664 |\n\nThese results indicate that while the initial (thresholding) algorithm provides a solid baseline, the refinement stage significantly reduces the estimation error. In the final version, we will include additional experiments to compare the performance before and after refinement.\n\n> 3. I understood that Diagonal Thresholding (DT) is computationally slower than your proposed algorithm. However, in Figure 2(b), we see that it is as fast as yours. Then, if we include the refinement using the Truncated power method (3.1.2), then we could accomplish the same performance with the same computationally efficiency?\n\n**Reply:** Thank you for your question. As explained in Section 3.2, both Diagonal Thresholding (DT) and our initialization algorithm incur the same order of computational cost of $O(np + nk^2)$, since they involve computing the entries of the sample covariance matrix and performing an eigenvalue decomposition.\n\nOur refinement stage utilizes truncated power iterations, each costing $O(np)$. In our experiments (see Tables 1-3), the number of iterations required for convergence is very small---typically no more than 10 iterations. Consequently, because the initialization stage dominates the overall cost, the additional refinement improves accuracy without compromising computational efficiency. This is further supported by Figure 2(b), where our two-stage method (red line) runs as fast as DT (green line).\n\n> 4.  Start with an introductory text explaining the single-spiked covariance model. It is missing the basic differences between PCA and sparse PCA.\n\n**Reply:** Thank you for your valuable suggestion. In the final version of our manuscript, we plan to revise the introductory text in Section 1 to improve clarity and provide a more comprehensive overview of the single-spiked covariance model. Specifically, we will:\n\n- Introduce principal component analysis (PCA) by emphasizing its importance, typical applications, and inherent limitations (e.g., the production of dense principal components).\n\n- Present sparse PCA as a remedy to these limitations by highlighting its role in generating more interpretable components.\n\n- Provide a detailed explanation of the single-spiked covariance model and add the citation (Johnstone, 2001) when presenting the mathematical formulation.\n\n- Include one or two examples of applications or problems where this model proves particularly useful.\n\nWe believe that these revisions will enhance the clarity and accessibility of the introductory section.\n\n>5.  Why diagonal thresholding does not fullfill Theorem 3.2.\n\n**Reply:** Thank you for your comment. With assumption $\\lambda = \\Omega({||v||} _{\\infty}^{-1})$, diagonal thresholding (DT) cannot achieve sample complexity $\\Omega(k \\log p)$. This is based on the following observations and analysis. \n\nThe sample complexity of DT is governed by the statistical gap defined in Eq.(5) of our paper. In contrast, the sample complexity of Algorithm 1 is governed by the statistical gap as Eq.(7) in our paper. A larger gap leads to a smaller sample complexity. If the infinity norm of $v$ is of constant order, then our gap depends only linearly on $\\min_{j \\in S} |v_j|$, improving order-wisely the gap of DT and consequently reducing the sample complexity $\\Omega(k^2 \\log p)$ to $\\Omega(k \\log p)$. However, for DT, the constant order of ${||v||} _{\\infty}$ does not improve its gap, suggesting that its sample complexity cannot be improved. \n\nIn conclusion, Algorithm 1 provides a way to incorporate the property of $v$ into estimation, while DT cannot do it according to our analysis. We will add this discussion in the final version.\n\n> 6. Other Comments:\n\n**Reply:** Thank you—we'll fix the typos.", "claims_and_evidence": "Overall, the paper mainly provides proper theoretical proves and derivations to theoretically demonstrate the computational complexity of the algorithm and its theoretical error. Additionally, it evaluates with synthetic data the main theoretical claims; one is missing regarding diagonal thresholding and I also miss experiments regarding the two stage refinement (see questions for authors).", "methods_and_evaluation_criteria": "Overall, the paper mainly provides proper theoretical proves and derivations to theoretically demonstrate the computational complexity of the algorithm and its theoretical error. Additionally, it evaluates with synthetic data the main theoretical claims; one is missing regarding diagonal thresholding and I also miss experiments regarding the two stage refinement (see questions for authors).", "theoretical_claims": "Theoretical proves are included in the supplementary, which I have not reviewed.", "experimental_designs_or_analyses": "I have checked the experimental analysis. From my perspective, experiments are properly aligned with the paper, but two main questions arises from two experiments that lack: one is missing regarding diagonal thresholding and I also miss experiments regarding the two stage refinement (see questions for authors).", "relation_to_broader_scientific_literature": "The key contributions of this paper build upon previous work in sparse PCA, particularly in improving sample complexity and computational efficiency. By introducing a novel thresholding algorithm for principal component initialization and truncated power iteration, the paper provides an efficient solution that achieves optimal statistical guarantees with significantly reduced sample complexity.", "essential_references_not_discussed": "Main references are included.", "other_strengths_and_weaknesses": "Overall, I think that the paper is good, but it requires addressing a few things to be a clearly ICML publcication. Please address the two comments in “Other Strengths And Weaknesses” and the three questions in “Questions For Authors”. If they are addressed, I would modify my recommendation to accept:\n1) I would rather start with an introductory text explaining the single-spiked covariance model, instead of just starting with the mathematical formulation; also missing including the citation (Johnstone, 2001). Moreover, for me, it is missing in this introductory text part the basic differences between PCA and sparse PCA/single-spiked covariance model: we assume that the data can be explained with a unique principal component, coefficients follow a standard normal (instead of are \"sampled\"), epsilon_i gathers the inherent noise or the information not explained by the single component, an example or two of applications/problems where this model is useful, etc.\n2) Sorry if I miss it in the paper, but I think it would be interesting to either include the prove that diagonal thresholding (Johnstone & Lu, 2009) requires n=Ω (k 2 logp) or either to include a discussion/counterexample illustrating that diagonal thresholding does not fulfill Theorem 3.2 (maybe this will be ever more interesting, as the other one is included in the original paper). This will showcase that your algorithm is effectively better than diagonal thresholding. It would be ok to prove it in the appendix and just reference it with one sentence.", "other_comments_or_suggestions": "3) \"reduce the parameter space of the model\" --> include here the notation that you use\n4) second stage incurs O(np) costs.. --> typo, two dots\n5) 2.1. Preliminaries --> explain you work with zero-mean observations\n6) to estimate a single sparse vector v, to identify a sparse unit vector w --> explain better the differences between v and w\n7) This approach is based on the observation that the expectation values of the diagonal entries of --> missing a reference to Eq. 2\n8) Is v^~ notation introduced in Equation (8)?\n9) When you define the Error in Section 4, you are repeating the same equation than Eq. 10, just reference it. Moreveor, this may be just a styling thing, but Error and F-score appear in another font and font size, which does not look really formal for a paper.", "questions_for_authors": "Overall, I think that the paper is good, but it requires addressing a few things to be a clearly ICML publcication. Please address the two comments in “Other Strengths And Weaknesses” and the three questions in “Questions For Authors”. If they are addressed, I would modify my recommendation to accept:\n10) You mention that \"operating under a specific condition on λ\", which is the implication of this restriction? It would be interesting to include it this discussion in the paper, either theoretically or experimentally.\n11) I understand that the two-stage algorithms theoretically guarantees the minimax optimal rate. However, until which level is this refinement necessary? In other words, until which level refines the your proposed thresholding algorithm? This is missing. For instance, include an experiment showing before and after the refinement using your Error metric.\n12) I understood that Diagonal Thresholding (DT) is computationally slower than your proposed algorithm. However, in Figure 2 b, we see that it is as fast as yours. Then, if we include the refinement using the Truncated power method (3.1.2), then we could accomplish the same performance with the same computationally efficiency?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper proposes a two-stage algorithm for sparse PCA under the single-spiked covariance model, achieving near-optimal 
$O(k \log p)$ sample complexity by assuming a large leading coordinate in the spike. The method combines a thresholding-based support recovery step with truncated power iteration for refinement.

The contribution is technically solid and closes a known computational-statistical gap under a plausible structural assumption. While the estimator is conceptually simple, the authors justify its analysis and show improved empirical accuracy compared to baseline variants. The rebuttal clarified several points, including robustness to the assumption, diagonal thresholding comparisons, and the role of the refinement step. The main limitation is that the setting is still a simplified model and results are only shown on synthetic data. The extra assumption helps the theory but may not always hold in real problems. However, given the strong theoretical contribution, I recommend an accept.
Decision: Accept (poster)