Efficient Sparse PCA via Block-Diagonalization
Authors: Alberto Del Pia
Dekun Zhou
Yinglun Zhu
Venue: ICLR 2025 Poster (2025)
Keywords: Sparse PCA
Block Diagonalization
Compurational Efficiency
Approximation Algorithms
Primary area: optimization

Abstract: Sparse Principal Component Analysis (Sparse PCA) is a pivotal tool in data analysis and dimensionality reduction. However, Sparse PCA is a challenging problem in both theory and practice: it is known to be NP-hard and current exact methods generally require exponential runtime. In this paper, we propose a novel framework to efficiently approximate Sparse PCA by (i) approximating the  general input covariance matrix with a re-sorted block-diagonal matrix, (ii) solving the Sparse PCA sub-problem in each block, and (iii) reconstructing the solution to the original problem. Our framework is simple and powerful: it can leverage any off-the-shelf Sparse PCA algorithm and achieve significant computational speedups, with a minor additive error that is linear in the approximation error of the block-diagonal matrix. Suppose $g(k, d)$ is the runtime of an algorithm (approximately) solving Sparse PCA in dimension $d$ and with sparsity constant $k$. Our framework, when integrated with this algorithm, reduces the runtime to $\mathcal{O}\left(\frac{d}{d^\star} \cdot g(k, d^\star) + d^2\right)$, where $d^\star \leq d$ is the largest block size of the block-diagonal matrix. For instance, integrating our framework with the Branch-and-Bound algorithm reduces the complexity from $g(k, d) = \mathcal{O}(k^3\cdot d^k)$ to $\mathcal{O}(k^3\cdot d \cdot (d^\star)^{k-1})$, demonstrating exponential speedups if $d^\star$ is small. We perform large-scale evaluations on many real-world datasets: for exact Sparse PCA algorithm, our method achieves an average speedup factor of 100.50, while maintaining an average approximation error of 0.61%; for approximate Sparse PCA algorithm, our method achieves an average speedup factor of 6.00 and an average approximation error of -0.91%, meaning that our method oftentimes finds better solutions.

Published as a conference paper at ICLR 2025
EFFICIENT SPARSE PCA VIA BLOCK- DIAGONALIZATION
Alberto Del Pia & Dekun Zhou University of Lisconsin-Madison (delpia,
Yinglun Zhu
University of California, Riverside yzhu@ucr.edu
dzhou4 4 4 @ sc. edu
ABSTRACT
Sparse Principal Component Analysis (Sparse PCA) is a pivotal tool in data analy- sis and dimensionality reduction. However, Sparse PCA is a challenging problem in both theory and practice: it is known to be NP-hard and current exact methods generally require exponential runtime. this paper, we propose a novel framework to efficiently approximate Sparse PCA by (i) approximating the general input covariance matrix with re-sorted block-diagonal matrix, (ii) solving the Sparse PCA sub-problem in each block, and (iii) reconstructing the solution to the original problem. Our framework is simple and powerful: it can leverage any off-the-shelf Sparse PCA algorithm and achieve significant computational speedups, with minor additive error that is linear in the approximation error of the block-diagonal matrix. Suppose g(k, d) is the runtime of an algorithm (approximately) solving Sparse PCA in dimension d and with sparsity constant k. Our framework, when in- tegrated with this algorithm, reduces the runtime to O (4 g(k, d²) + dd), where d² S d is the largest block size of the block -iik--iaaonal matrix. For instance, integrat- ing our framework with the Branch-and and--oundd algorithm reduces the complexity from s(k, d) O(k3 d to O d. (d* y-)), demonstrating exponential speedups if d is small. We perform large-scale evaluations on many real-world datasets: for exact Sparse PCA algorithm, method achieves average speedup factor of 100.50, while maintaining an average approximation error of 0.61 %; for approximate Sparse PCA algorithm, method achieves average speedup factor of 6.00 and an average approximation error of -0.91 %, meaning that our method oftentimes finds better solutions.
INTRODUCTION
In this paper, we study the Sparse Principal Component Analysis (Sparse PCA) problem, a variant of the well-known Principal Component Analysis (PCA) problem. Similar to PCA, Sparse PCA involves finding a linear combination of d features that explains most variance. However, Sparse PCA distinguishes itself by requiring the use of only < d many features, thus integrating sparsity constraint. This constraint significantly enhances interpretability, which is essential in data analysis when dealing with large number of features. Sparse PCA has found widespread application across various domains, including text data analysis (Zhang & Ghaoui| 2011). cancer research (Hsu et al.) 2014), bioinformatics (Ma & Dai (20TT), and neuroscience (Zhuang al. 2020). For further reading, we refer interested readers to comprehensive survey on Sparse PCA (Zou & Xue 2018).
While being important and useful, Sparse PCA is a challenging problem- it is known be NP- hard (Magdon-Ismall (007)). Solving Sparse PCA exactly, such as through the Branch-and- -ound algorithm (Berk& & Bertsimas, 2019), requires worst-case exponential runtime. Moreover, while there are near-optimal approximation algorithms that exhibit faster performance empirically, they still require worst-case exponential runtime (Bertsimas et al.| 2022 Li & Xie (2000) (0000/Cory-oaa Kory-Wright & Pauphilet 2022 Dey et al. 2022b Zou et 2006 Dey et al. 2022a). This complexity impedes their applicability to large-scat scale datasets. the other hand, there are numerous polynomial time
approximation algorithms for Sparse PCA (Chowddhry al. 2020 Li & Xie 2020 Papailiopoulos
cet al. 2013 Chan et al. 2013) Del Pial [s22]] et al. 2015). These algorithms typically face
trade-ofis between efficiency and solution quality. Some algorithms are fast but yield sub-optimal

---

Published as a conference paper at ICLR 2025
solutions, others provide high-quality solutions at the cost of greater computational complexity, and few achieve both efficiency and accuracy but only under specific statistical assumptions.
In this paper, introduce novel framework designed to efficiently approximate Sparse PCA through matrix dinnokkiinonliizationn Our framework facilitates the use of off-the-shelf algorithms in plug-and- aad-play manner. Specifically, the framework comprises three principal steps: (i) Matrix Preparation: Given a general input covariance matrix to Sparse PCA, we generate a re-sorted block-diagonal matrix approximation. This includes finding a denoised input matrix by thresholding, grouping non-zero entries into blocks, and slicing corresponding blocks in original input matrix. (ii) Sub-Problem Solution: Solve Sparse PCA sub-problems using any known algorithm in each block with lower dimension, and find out the best solution among sub-problems; (iii) Solution Reconstruction:
Reconstruct the best solution to the original Sparse PCA problem. We illustrate these steps in Fig. [ At a high level, our approach involves solving several sub-problems in smaller matrices instead of directly addressing Sparse PCA in a large input matrix (Step (ii)), which significantly reduces computational runtime. This is made possible by carefully constructing sub-problems (Step (i)) and efficiently reconstructing an approximate solution to the original problem (Step (ii)).
@ Matrix Preparation
Sel to sub- probioms
Sal original problem
Algonthm
Support resovery using inverse (ar-l)
iBeat obj valuel
denolsing sorting
Figure 1: Illustration of our proposed approach, given a 9 x 9 covariance input matrix A. (i) Entries away from zero are highlighted the upper matrix (original input A). Then group those entries in blocks, zero out outside entries, sort the matrix, and obtain the lower block-diagonal approximation. Heatmaps present values of matrix entries. The axes are indices of A; (ii) Extract sub-matrices from the block-diate diagonal approximation, and solve sub- sub-problems via suitable Sparse PCA algorithm; (iii) Select the solution with the highest objective value obtained from the sub-problems. Construct solution for the original Sparse PCA problem by mapping its non-zero entries to their original locations using the inverse mapping of the sorting process, and setting all other entries to zero.
Our framework significantly speeds up the computation of Sparse PCA, with negligible approximation of our framework
error. We theoretically quantify the computational speedup and approximation in Sections 3 and and conduct large scale empirical evaluations in Section 5 Next, we illustrate the performance of method via concrete example. Suppose matrix being noisy sorted block-diagonal matrix with 10 blocks, each of size 100, and is asked to solve Sparse PCA with a sparsity constant The Branch-and -ouund algorithm (Berk & Bertsimas 2019) takes over 3600 seconds and obtains solution with objective value 13.649. Our framework integrated with Branch-ant Bound algorithm reduces runtime 17 seconds, speedup by a factor of 211, and obtains solution with objective value 13.637, with an approximation enror of 0.09%.
is given an input covariance
Our contributions. We summarize our contributions in the following:
We propose novel framework for approximately solving Sparse PCA via matrix block- diagonalization. This framework allows users to reuse any known algorithm designed for Sparse PCA in the literature in very easy way, i..., simply applying the known algorithm to the blocks of the approximate matrix, and accelerates the computation. To the best of our knowledge, our work is the first to apply block-diagonaliizatin to solving Sparse PCA problem.
2

---

Published as a conference paper at ICLR 2025
We provide approximation guarantees for our framework. We show that when integrated with an approximation algorithm, our method can generally obtain a better multiplicative factor, with an additional cost of an additive error that is linear in the approximation error of the block-diagonal approximate matrix. Under certain assumptions, the additíve error can also be improved.
We also perform time complexity analysis showing that for an (approximation) algorithm with runtime g(,, d), our framework could reduce the runtime to ( 9(k; d²) + d²), where d is the dimension of Sparse PCA input covariance matrix, k is the sparsity constant, and d' is the largest
block size in the block-diagonal approximate matrix.
We show that for statistical model that is widely studied in robust statistics, there exists an efficient way to find a block-diagonal approximation matrix. Moreover, we provide a very simple way extend to the setting where the statistical assumptions are dropped and where computational resources are taken into consideration.
We conduct extensive experiments to evaluate our framework, and our empirical study shows that (i) when integrated with Branch-and-Bound algorithm, our framework achieves an average speedup factor of 100.50, while maintaining an average approximation error of 0.61 %, and (ii) when integrated with Chan's 's algorithm (Chan al.) 2015), our framework achieves average speedup factor of 6.00 and an average approximation error of -0.91 %, meaning that our method oftentimes finds better solutions.
Our techniques. In this paper, we also introduce several novel techniques that may be of independent interest. We characterize structural property of an optimal solution Sparse PCA, when given block-diagonal input matrix. We establish that solving Sparse PCA for block-diagonal matrix input equivalent to addressing individual Sparse PCA sub-problems within blocks, which consequently provides significant computational speedups. This result is non-trivial: while the support of optimal solutions might span multiple blocks, we prove that there always exists an optimal solution whose support is contained within single block, ensuring the efficiency of our framework. Additionally, leveraging the defined structure, we propose novel binary- search based method to identify the best block-diagonal upproximation, particularly when computational resources limited.
Paper organization. This paper is organized as follows. In Section 2 we introduce Sparse PCA problem formally, and provide motivation of our proposed framework. In Section we provide details of operational procedures of our proposed framework, and develop guarantees of approximation errors well as runtime complexity. In Section 4 we further extend our algorithms to learning problems either under (i) a statistical model, or under (ii) model-free setting, and provide theoretical guarantees respectively. In Section 5) run large-scale empirical evaluations and demonstrate the efficacy of our proposed framework when integrated with both exact and approximate Sparse PCA algorithms. We defer discussion of additional related work, including existing block- binok-iignnalizztioon methods in the literature, detailed proofs, and additional empirical results to the appendix.
2 PROBLEM SETTING AND MOTIVATION
In this section, we define formally the problem of interest in this paper, and then we explain the motivation behind our framework. We first define the Sparse PCA problem as follows:
OPT == max x Ax s.t. ||0||2 = 1, || < k.
(SPCA)
Here, a symmetric and positive semidefinite input matrix p and a positive integer S d are given. & is known as the sparsity constant, the upper bound on the number of nonzero entries in solution x.
Motivating idea in this paper. The motivation for our method is illustrated through the following example, which demonstrates the efficiency of solving sub-problems: Suppose that user is given an input block-diagonal covariance matrix A of size d, comprising d/d' many d' x d' blocks (assume that d divisible by d²), and is asked to solve [SCCA] with sparsity constant k. The user solve
can simply apply exact algorithm (c.g., Branch and-Bound algorithm) to solve the large problem in time Okk (t) = O d²).. The user could also solve d/d' many sub- broblemssems in each block, potentially obtaining an optimal solution (a fact we rigidly prove in Section 3.2.1) in time ((d/d")) d 3.3 (d*)*). This segmented approach yields significant speedup factor of O((d/ ((dd/=*)* d=*)-
However, not all input matrices are of the diagonla structure, even after permuting its rows and columns. One goal of our work is to develop a reasonable way construct (sorted)
3

---

Published as a conference paper at ICLR 2025
block-diagonal matrix approximation À to A such that they close enough. In this paper, interested in finding the following E-matrix approximation to A:
Definition (e-matrix approximation). Given matrix E Rdxd, denote by Aig the (1, entry of A. An e-matrix approximation of A is defined as a d x d matrix A such that
4- A = max i.jE(d
Aij
Ass Sc.
Moreover, we define the set of all c-matrix approximations of A to be B(A,e).
In other words, A E ((A, s) if and only if Å belongs to an -bal centered at A with radius e. Our overarching strategy that if the input matrix not block-diagonal, we find a matrix approximation in the ball so that it could be resorted to be block block-diagonal. We then solve the Sparse PCA problem using this resorted matrix and reconstruct the approximate solution to the original problem, taking advantage of the computational speedups via addressing smaller sub-problems.
Additional notation. We adopt con-asymptotic big-oh oh notation: For functions f.g: + R..
we write f (g) (resp. ʃ R(g)) if there exists constant C >> 0 such that [(2) < Cg(z)
(resp. f(z) 2 Cg(z)) for all 2. For an integer E N, we let [x] denote the set [1. (1.2....,1) ...
For vector E Fd
2; to denote its i-th entry. We denote supp(z) == {i [d] 1 & 0}
the support of z. For S q S co, we denote |z]lg the g-norm of 2, i... ||~||g : E 1//
and || == [ | For a matrix A E Ruxm, we use Ag to denote its (i.j)-th entry. We
detine ||A|| == max,j |Ayl- For index sets S [ [n.] T C [m]. we denote As, the |S| x [2] sub-matrix of A with row index S and column index T. For matrices A; E Rd.xd for i = 1.2..... we denote diag(A, : Ap) the (Ei= d,) (Ei= d,) block- diigonaa matrix with blocks A1, .2..A,,..4444........ Ap. For square matrix B E Raxd we define the size of B be
OUR APPROACH
In this section, we introduce our framework for solving SPCA diagonalization technique. We explain in detail the operational procedures of our framework in Section 3.1 In 3.22 we characterize the approximation guarantees and computational speedups of our framework. Note that, in this section, our proposed procedures require predefined threshold we extend our algorithms to learning the threshold in statistical model, as well
approximately via block-
E > 0. In Section 4 as in a model-free setting.
3.1 ALGGORITHMS
In this section, we outline the operational procedure of our proposed framework for solving SPCA approximately. The procedure commences with the input matrix A E Rgdxd and predefined threshold value > 0. Initially, a denoised matrix A is derived through thresholding as described in Algorithm| D Subsequently, several sub-matrices and their corresponding sets of indices are extracted from A and A by employing Algorithm 2 which groups index pair (i,j) )) if A,, * O. This grouping mechanism can be efficiently implemented using Depth-First Search approach, which iterates over each unvisited index i and examines all non-zero AC,, and then visit each corresponding index
that is unvisited. The computational complexity of Algorithm 2 is bounded by (dd). Finally, an approximate solution to iPCAAs obtained by solving several Sparse PCA sub-problems in each block using certain (approximation) algorithm A, mapping the solution back to the original problem, and finding out the best solution. The operational details are further elaborated in Algorithm
3.2 APPROXIMATION GUARANTEES AND COMPUTATIONAL SPEEDUPS
In this section, we provide analysis of approximation guarantees and speedups of our framework Al- Before showing approximation guarantces, we first define the e-intrinsic dimension of A:
gorithm 3
We that line 6 in Algorithm if k is larger than the size of A., a PCA problem is solved instead.
4

---

Published as a conference paper at ICLR 2025
Algorithm Denoising procedure via thresholding
Input: Matrix RRddd, threshold e
2: Output: Denoised matrix A E Rdxd
3: for (i,j) (1, to (d, do
4: If [Ayy > set A f Aij. otherwise set Aig +- 0
5: Return A
Algorithm 2 Matrix block boack-diagonlization
1: Input: Matrices A, .A Rdxd
2: Output: Sub-matrices and their corresponding lists of indices
3: for each pair of indices () [d] [d] do
4: Group and j together if A
5: Obtain sets of indices {S.) that are grouped together
6: Return sub-matrices As,5s, and the corresponding S;, for E [p)
Definition 2 (s-intrinsic dimension). Let A E Rddxd The e-intrinsic dimension of A, denoted by int dim(A,)), is defined as:
int (A,(), min lbs(B), BEB(A.e)
where Ibs(B) denotes the largest block size of B, ie., the largest among {[St|eett after rumning B.
lines 3-5 Algorithm @ with A
As we will see later this section, this concept plays erucial role in characterizing both the quality of the approximate solution obtained through our framework, Algorithm and the computational speedup achieved by our framework. Moreover, matrix A can be efficiently identified such that (Ā) int dim (A,e) E) using Algorithm shown in Lemma
Lemma 1. Let Raxdd and 0. Given inpul (A, e), Algorithm oudputs an e-approximation of A, denoted as A, such that lbs A) = int dim(A, 8) in time (dd)).
3.2.1 APPROXIMATION GUARANTEES
In this section, we provide approximation guarantees for our framework Algorithm 3 that Theorem 1 that, an optimal solution to SSCAA with input (A, and an optimal solution to (A k) are close to each other in terms of objective value:
We first show
Theorem 1. Let A E Rdxd be symmetric, Is S d, and e >> O. Denote by A the output of Algorithm / with imput (A, T an optimal solution to SPCS with input (A k), and x² an optimal solution to SPCA with input (A, k). Then, follows that
((x' Ax" Txx)
AT
2k E.
Algorithm 3 Efficient Sparse PCA via block-diagonalizztion
1: Input: Matrix A dxd positive integer d, threshold 0, algorithm A for (approxi- mately) solving SPCCA
SSCAA
2: Output: Approximate solution to SPCA
3: Obtain denoised matrix A using Algorithm with input (A,s)
4: Obtain sub-matrices
(A, ) and corresponding index sets {S.). via Algorithm with AA)
5: for i = to do
6: Solve SPCA approximately with input (Ā,, k) using algorithm and obtain solution x;
7: Construct y; Rd by placing (x:)y at (yr) (.)) for j in block indices St and zero elsewhere
8: Return the best solution among {y:P_ that gives the largest v.T Ay

---

Published as a conference paper at ICLR 2025
Then, in Theorem 2 we show that solving SPCA exactly with block-diagonal input matrix Á is equivalent to addressing sub-problems in blocks within A:
Theorem 2. Let À = diag( A1, A Ap) be a symmetric block-diagone matrix. Denote OPT to be the optimal value to SPCA with input pair (Ā, k)). Let OPT; to be the optimal value to SPCA with imput pair (Ai, k). for i E [p]. Then, one has OPT = maxie(p OPT;
Remark 1. Theorem? 2 is particularly non-trivial: without this result, one might think the support of any optimal solution to SSCA may span different blocks. Solving SScpssSCAWWAA with a block- block-diagonal input would consequently involve iterating over all possible sparsity constants (that sum up to k) for each
blockddiggonal sub-matrix, making our framework proposed in Section 3.7 nearly impractical.
We now extend Theorem to incorporate the use of various approximation algorithms. We note
that our extension improves the general bound proposed in Theorem [ since Algorithm 3 finds specific block-diagonal approximation of A. Furthermore, Theorem 3 suggests that our frame- work Algorithm [ can also be highly effective when combined with other approximation algorithms. Specifically, it has the potential to find better solution, particularly when the value of e is small.
Theorem 3. Let A E Rdxd be symmetric and positive semidefinite, let k be positive integer, and denote by x* E Hd an optimal solution to SPCA with input (A, k). Suppose that an algorithm for SPCA with input (A, k) finds an approximate solution x E R to SPCA with multiplicative factor m(k, d) 2 1 and additive error a(k.d) d) Z i.e., one has TT Ax 2 x Ax" /m[(k, d) a(k, d). Furthermore, we assume that m(k, d) and a(k, d) is non-decre with respect to d. For E > 0, denote by y E I&d be the outpul of Algorithm with imput tuple (A, k, A). Then, one has
(xx)T Axx
y' Ay 2 m. (ke, int dim (A))) a (Ic, int din((,,
m (k, int ()) KE-
Remark 2. Denote by d' := int dim((,, s) S d. Theorem 3 iimplies that our framework, Algorithm 3 when imputted with (A, e, A), could find approximate solution with better multiplicative factor m (k, ). and additive error (Ic, d²)) ) m(k.as).
ke
Note that the additive error is also improved if
a(k, d) 2 a (ki, int dim(), +
m. (ki int dim A, =)) kee
In addition, we have found instances where our framework can indeed obtain better solutions when integrated with some approximation algorithm, as shown in Table | in Appendix C Finally, we summarized in Table 7 the change of multiplicative factor and additive error when integrated our framework to some existing (approximation) algorithms for SSPAA
Table 1: Comparison of multiplicative factors (MF) and additive errors (AE) when integrating our framework with different algorithms. The MISDP and Greedy algorithms are both studied in Li &
Xie (2020). Note that our MF is always no larger than the original MF.
Algorithm
MF
AE
Our MF
Our AE
Chan et al. 2015)
min( Vk, dd/3) min (k, d/kb) k
min( VR, int dim (a,,(A, A,,)//3) min (kc int dim A,k)/k)
k MF
MISDP
ke/Our MF
Greedy
k
Exact algorithms
ke ke + 5
Asteris et al. (2015)
1
8
3.2.2 COMPUTATIONAL SPEEDUPS
In this section, we characterize the runtime complexity of Algorithm 3 We will see that the e-intrinsic dimension also plays an important role in this section.
Proposition 1. Let A Raddd and let k be a positive integer. Suppose that for specific algorithm A. the time complexity of using for (approximately) solving SPCA with input (A, k)) is a function Then, for
g(k, d) which is convex and non-decreasing with respect to d. We assume that 9(k, 1) a given threshold e > O, the runtime of Algorithm 3 with input tuples (A, K, E, A, is at most
0 ( int dim(A,e)
g (ks int dim(A, (A,e)) + d²

---

Published as a conference paper at ICLR 2025
We remark that the assumption g(k, d) is convex and non-decreasing with respect d in Proposition iss a very weak assumption. Common examples include d' and do log d with a 2 1.
Remark 3.
Remark 4. this remark, we explore the computational speedups provided by Algorithm 3 by modeling int dim(A, E) as a function of A and e. We assume without loss of generality that ||
1, thereby setting int dim (A, 1) = 0. It is important to note that while int dim(A, E) is discontinuous, for the purpose of illustration in Table 2 it treated continuous function of e.
Table 2: Comparison of original runtimes and current runtimes in our framework. We assume that the algorithm has runtime (dd)), ) for some a 2 2.
int dim dim(A,s)
Runtime
Our Runtime
Speedup factor
d(1 e ((1--))
(dd))
Odd - - -001)=))a- a-11)
O((1 e -[1-=))] 1-a)
d(1 E)
(dd))
Old (1 =1
O((1 - 8)1
d(1 -)
(dd))
Old 10 O(11 (m)m
4 EXTENSIONS
Our algorithm presented in Section 3.1 takes a known threshold as an input. In this section, remove this constraint by extending our algorithm to various settings. Specifically, in Section we consider learning under statistical model, where one can do the denoising of a matrix easily. In we propose a practical learning approach without any assumption.
Section 4.2
4.1 LEARNING UNDER A STATISTICAL MODEL
A crucial aspect of our proposed framework involves determining an appropriate threshold, s, as required by Algorithm While specifying such threshold is generally necessary, there are many instances where it can be effectively estimated from the data. In this section, we discuss model
commonly employed in robust statistics (Comminges et al. [221] Kotekal & Gao [023 2024). This
model finds broad applications across various fields, including clustering (Chen & Witten 2022 and sparse linear regression (Minsker et al.] 2022).
Before presenting the formal definition, we first establish some foundational intuition about the model. Conceptually, if the input data matrix A is viewed as an empirical covariance matrix, and based on the intuition that, a feature is expected to exhibit strong correlations with only few other features, it follows that should be, or at least be close to, block-diagonal matrix. This observation motivates us for considering Model[]
Model 1. The input matrix A E [Rdxd is close to a block-diagonal symmetric matrix À dRd i.e., A = À E for some symmetric matrix E E RRddd such that A is positive semidefinite and
(i) for S j, each entry Fag is drawn from i.i.d. .-sub-Gaussinn distribution, with mean zero. and variance o2, where @ >> O and o² being unknown, but an upper bound u 2 e² /o2 is given; (i) write Ā = diag (A1, A2,,,,,.,,4 A,). the size of each diagonal block Ag is upper bounded by d', where d' < Cd for some known constant C > 0 and a E (0, 1), but d' is unknown.
In assumption (i) of Model [ the sub-Gaussian distribution is particularly favored as this class of distributions is sufficiently broad to encompass not only the Gaussian distribution but also distributions of bounded variables, which are commonly observed in noise components of real-world datasets. Assumption (ii) of Model uses the constant quantify the extent which features are clustered.
The procedure for estimating ||] >c in Model involves an estimation of the variance a2, which is adapted from Comminges et al. (022)), and then multiply a constant multiple of vlog d. We present the details in Algorithm|4]
In the next theorem, we provide a high probability approximation bound, as well as a runtime analysis, for Algorithm 3 juiig the threshold found via Algorithm 4 The results imply that in Model one can efficiently solve SSCA) by exploiting the hidden sparse structure.
Proposition 2. Consider Model | and denote by the output of Algorithm with input tuple (A, C, w). Let k S d. be positive integer, and assume that d satisfies that dl > Co (C + 1)
7

---

Published as a conference paper at ICLR 2025
Algorithm 4 Estimation of || || in Model
1: Input: Matrix A E Rdxd parameters C 0, E (0, 1), and in Model
2: Output: A number
as an estimate of F
Divide the set of indices {(i.j) E [d] [d] : into m [2CC 2)d disjoint subsets
B1, B2........ Bm randomly, ensuring each subset B; has cardinality lower bounded by
d'-d
2m
4: Initialize an array S to store subset variances
5: for 1 to do
6: S: t [B (L,)EB A
7: x² (S)
8: Return t 2uō v logd
u logg(CC 8) for some large enough absolute constant Co >> 0. Then, the following holds with probability least 2d exp dl+a dl+a/(4C 4)}:
(i) Denote by X= the optimal solution to SPCA with with input (A, k). For an (approximation) algorithm A with multiplicative factor m (k, d) and additive error a(k, d) 2 O, where the functions
m and a is non non-decean den-aiinnnng with respect to d, the output of Algorithm with input tuple
(A, k, E, satisfies that
T A
(1
2 m (k, dz)
V Ay in (ki dk)
(k,
ke.
(i) If, in addition, the time complexity of using A for (approximately) solving SSCA with input (A, k) is function g(k, d) which is convex and non-decreusing with respect to d and satisfies g(k, Then, the runtime of Algorithm with input tuples (A, k, 8, A) most
o
(I d d (ke dk) +
Finally, we remark that, the assumption d' 21-o > Co (C' + 1) log(8C'+ + +8) is true when d is large enough, as Co. C, and u are usually some constants not changing with respect to the dimension d.
4.2 PRACTICAL LEARNING APPROACH
In this section, we describe practical algorithm to estimate without prior statistical assumptions while simultaneously providing an approximate solution to (SCCA)
thresholded matrix A is generated via Algorithm| from the original matrix A. This matrix A is then utilized to approximate (SCCA using specified algorithm A. Our proposed approach in this section involves initially selecting an value within a predefined range, executing Algorithm with this E, and subsequently adjusting E in the following iterations through binary search process. Since efficiency is of our major interest; thus, ask users to provide parameter do that indicates the maximum dimension they are willing to compute due computational resource limit. We defer the detailed description of this binary-search based algorithm to Algorithm 5 For clarity, let us denote by OPT(E) the objective value of (x*)7 Ar where rf is the output of Algorithm 3 with inputs (A, k. E, A).
Recall from Algorithm 3 that for given threshold
We have the characterization of approximation eπror and runtime complexity of Algorithm 5
Theorem 4. Let A d be symmetric and positive semidefinite, and let be a positive integer.
Suppose that an algorithm A for SSCAA with input (A,k) finds an approximate solution x E IRd to SPCA such that has multiplicative factor m(k, d) Z and additive error a(k, d) 2 0, with m(k, d) and alk, d) being non-de non-decreasing with respect to d. Suppose that Algorithm with input (A, k. A, do) terminates with OPT(=*) 2 Then one has
OPT(E*) Z m(k: do) OPT(O) - a(k, do) m(k, do)
ke*
Suppose that for A, the time complexity of using for (approximately) solving SPCA with inpur (A, k) is a function g(k, which convex and non-decreasing with respect d, g(k, 6, = and
8

---

Published as a conference paper at ICLR 2025
Algorithm 5 A practical version of Algorithm without given threshold
1: Input: Matrix A E [Rdxd, positive integer S d, tolerance parameter 6, (approximate) algorithm
A for SPPA and integer do
2: Result: Approximate solution to Problem SPCA
3: L+ 0, U <- ||AA]] find OPT(U) via A
4: while U -
do
5:
(1/
6: Find A via Algorithm] 1 and the largest block size d' := lbs(A*)
7: if problems with largest block size 2 d² has been solved before then U E; continue
8: if d' > dg then I E; continue
9: Find (PT(s) via Algorithm with algorithm
10: if (e) does not improve then U +
return best OPT(E) found and the corresponding x²
g(k, 0) 0. The runtime for Algorithm
5
is most
log
(14//2
||Al| 6
(1 do
(ke do) + de) )
Remark 5. In this remark, we discuss the practical advantages of Algorithm
The first direct advantage is that Algorithm d not rely on any assumptions, enabling users to apply the frumework to any dataset. Moreover, the framework can be applied to scenarios where an improved estimation of the binary search range is known. For instance, if users possess prior knowledge about an appropriate threshold E for binary search, they can refine the search interval in Algorithm by setting the lower and upper bounds to L := aē and U := bĚ for some proper 0 I a < b to improve the efficiency by a speedup factor of (ll (oo)//oo((b-- Secondly, detailed in Theorem 4 is that Algorithun achieves similar trade-off between computational efficiency and approximation error as Algorithm which are proposed in Theorem 3 and Proposition] 7 The time complexity of Algorithm compared Algorithm increases only by a logarithmic Jactor due to the binary search step, thus preserving its efficiency.
a)E/8)).
The third advantage is that this algorithm considers computational resources, ensuring that each invocation of Algorithm remains efficient.
5 EMPIRICAL RESULTS
In this section, we report the numerical performance of our proposed framework.
Datasets. We perform numerical tests in datasets widely studied in feature selection. including datasets CovColon from Alon et al. (1999). LymphomaCovl from Alizadeh et al. (1998), Red- 500 and Reddit2000 from Dey al. (222b)) LeukemiaCov, LymphomaCov2, and ProstateCov from Dettling (2004), ArceneCov and DorotheaCov from Guyon et al. (2007), and GLI85Cov and GLABRAI80Cov from Li (2020). We take the sample covariance matrices of features in these datasets as our input matrices A. The dimensions of A range from 500 to 100000. Most of the matrices are dense and hence are not block-diagonal matrices, as discussed in Appendix C.T
We use the state-of-the-an exact algorithm Branch-and-Boot (Berk & Bertsimas 2019) and the state-of-the-art polynomial-time approximation algorithm (in terms multiplicative factor) Chan's algorithm (Chan et al. 2015) as our baselines.
Baselines.
Notation. We calculate the appreximation error between two methods using the formula (obj 100%, where objnase is the objective value obtained by a bascline algorithm, while objours is the objective value obtained by our framework integrated with the baseline algorithm. We calculate the speedup factor by dividing the runtime of a baseline algorithm by the runtime of Algorithm integrated with the baseline algorithm.
Results. In Table 3 report the performance of Algorithm 5 when integrated with Branch-and Bound algorithm. We summarized the performance results for = 2, 5, 10, 15 across various datasets. For each dataset, we report the average approximation error, average speedup factors, and

---

Published as a conference paper at ICLR 2025
their corresponding standard deviations. Notably, the average speedup factor across all datasets is 100.50, demonstrating significant computational efficiency. Our theoretical analysis anticipates a trade-off between speedups and approximation enrors; however, the empirical performance of our framework in terms of approximation accuracy is remarkably strong, with an average error of just 0.61% across all datasets.
In Table we report the performance of Algorithm 5 larger datasets with larger k's, when integrated with Chan's algorithm (Chan al. 2015). We summarized the performance results for 200, 500, 1000, 2000 across various datasets. Chan's algorithm terminates in time O(d²) and
returns approximate solution to SPCA with a multiplicative factor min ( Vk, d'/3). The average speedup factor 6.00, and grows up to 14.51 as the dimension of the dataset grows, while the average emror across all four datasets remains remarkably below 0%, meaning that our framework finds better solutions than the standalone algorithm.
These findings affirm the ability of our framework to handle diverse datasets and various problem sizes both efficiently and accurately. For a more detailed discussion, additional information on settings, and further empirical results, refer readers to Appendix [
Table 3: Summary of average approximation errors and average speedup factors for each dataset, compared to Branch-and-Boud Bound. Standard deviations are reported the parentheses.
Dataset
Avg. Error (%) (Std. Dev)
Avg. Speedup (Std. Dev)
CovColonCov
0.00 (0.00)
25.84 (38.69)
LymphomaCovi
0.00 (0.00)
22.64 (15.87)
Reddit1500
0.00 (0.00)
365.12 12 (238,24)
Reddit2000
0.06 (0.13)
278.69 (164.93)
LeukemiaCov
0.00 (0.00)
71.39 (83.17)
LymphomaCov2
4.82 (5.78)
6.95 (4.35)
ProstateCov
0,00 (0.00)
19.90 (30,82)
ArceneCov
0.00 (0.00)
13.50 (14.17)
Overall
0.61 (2.42)
100.50 (163.60)
Table 4: Summary of average approximation errors and average speedup factors for each dataset, compared to Chan's is algorithm. Standard deviations are reported in the parentheses. Datasets are sorted in ascending order according their dimensions.
Dataset
Avg. Error (%) (Std. Dev)
Avg. Speedup (Std. Dev)
ArceneCov
-0.08 (0.10)
0.86 (0.28)
GLI85Cov
-1.06 (1.22)
1.58 (0.21)
ELABBRAI8COO
-0.19 (0.26)
7,07 (0.62)
DorotheaCov
-2.30 (0.62)
14.51 (0.50)
Overall
-0.91 (1.11)
6.00 (5.66)
6 DISCUSSION
In this paper, we address the Sparse PCA problem by introducing an efficient framework that accommodates any off-the-shelf algorithm designed for Sparse PCA in plug-and- fashion. Our approach leverages matrix block- diagonnalizzztinn, thereby facilitating significant speedups when integrated with existing algorithms. The implementation of our framework is very easy, enhancing its practical applicability. We have conducted a thorough theoretical analysis of both the approximation error and the time complexity associated with our method. Moreover, extensive empirical evaluations on large-scale datasets demonstrate the efficacy our framework, which achieves considerable improvements in runtime with only minor trade-offs in approximation error.
Looking ahead, an important direction for future research involves extending our framework to handle the Sparse PCA problem with multiple principal components, while maintaining similar approximation and runtime guarantees. This extension could further broaden the applicability of our approach, addressing more complex scenarios in real-world applications.
10

---

Published as a conference paper at ICLR 2025
7 ACKNOWLEDGEMENTS
A. Del Pia and D. Zhou are partially funded by AFOSR grant FA9550-23-1-0433.. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the Air Force Office of Scientific Research.
REFERENCES
Farid Alizadeh, Saan-Pierre A Haeberly, and Michael L Overton. Primal-dual interior prriioor-ponnttte methods for semidefinite programming: convergence rates, stability and numerical results. SIAM Journal Optimization, 8())746-766,, 1998.
Uri Alon, Naama Barkai, Daniel A Notterman, Kurt Gish, Suzanne Ybarra, Daniel Mack, and Amold J Levine. Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays. Proceedings of the National Academy of Sciences, 96(1 96(12))6745-6750 1999.
Arash Amini and Martin J Wainwright. high-dimensional analysis of semidefinite relaxations for sparse principal components, In IEEE International Symposium on Information Theory, pp. 2454-2458, 2008.
Megasthenis Asteris, Dimitris Papailiopoulos, Anastasios Kyrillidis, and Alexandros G Dimakis. Sparse pea via bipartite matchings. Advances in Neural Information Processing Systems, 28, 2015.
Lauren Berk and Dimitris Bertsimas. Certifiably optimal sparse principal component analysis. Mathematical Programming Computation, 11:381-420) 2019.
Dimitris Bertsimas, Ryan Cory-Wright, and Jean Pauphilet. Solving large-scale sparse pea to certifiable (near) optimality. Journal of Machine Learning Research, 23(13):1-35, 2022.
Siu On Chan, Dimitris Papailiopoulos, and Aviad Rubinstein. On the worst-case approximability of sparse pea, arXiv preprini urXiv:1507.05530 2015.
Yiqun Chen and Daniela M Witten. Selective inference for k-means clustering. arXiv preprint arXiv:2203.15267, 2022.
Agniva Chowdhury, Petros Drineas, David Woodruff, and Samson Zhou. Approximation algorithms for sparse principal component analysis. arXiv preprint (xXi::2006.12748, 2020.
Laëtitia Comminges, Olivier Collier, Mohamed Ndaoud, and Alexandre Tsybakov. Adaptive robust estimation in sparse vector model. The Annals of Statistics, 4943):1347-1377, 2021.
Ryan Wtight and Jean Pauphilet. Sparse pea with multiple components. arXiv preprint
arXiv:2209.14790) 2022.
Alberto Del Pia. Sparse pea fixed-rank matrices. Mathematical Programming, pp. 1-19, 2022.
Yash Deshp, Andrea Montanari, et al. Sparse pea via covariance thresholding. Journal of Machine Learning Research, 17(141):1-4,, 2016..
Marcel Dettling. Bagboosting for tumor classification with gene expression data. Bioinformatics, 20 (18):3583-3597 2004.
Emilie Devijver and Mélina Gallopin. Block-diagonal covariance selection for high-dimensional gaussian graphical models. Journal of the American Statistical Association. 1136521):306314,, 2018.
Santanu Dey, Rahul Mazumder, and Guanyi Wang. Using relaxxatio and integer programming to obtain dual bounds for sparse pea. Operations Research, 1013):114-1192,, 2022a.
Santanu S Dey, Marco Molinaro, and Guanyi Wang. Solving sparse principal component analysis with global support. Mathematical Programming, pp. 1-39, 2022b.

---

Published as a conference paper at ICLR 2025
Yunzi Ding, Dmitriy Kunisky, Alexander S Wein, and Afonso S Bandeira. Subexponential time algorithms for sparse pea. Foundations of Computational Mathematics, pp. 1-50, 2023.
Tommaso d'Orsi, Pravesh K Kothari, Gleb Novikov, and David Steurer. Sparse pea: algorithms, adversarial perturbations and certificates. In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pp. 553-564. IEEE, 2020.
Salar Fattahi and Somayeh Sojoudi. Graphical lasso and thresholding: Equivalence and closed-form solutions. Journal of machine learning research, 20(10):1-44, 2019.
Jiashi Feng. Zhouchen Lin, Huan Xu, and Shuicheng Yan. Robust subspace segmentation with block-diagonal prior In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 38 3818-3825, 2014.
Isabelle Guyon, Jiwen Li, Theodor Mader, Patrick A Pletscher, Georg Schneider, and Markus Uhr. Competitive baseline methods set new standards for the nips 2003 feature selection benchmark. Puttern recognition letters, (8(12):1438-14 2007.
William W Hager, Dzung T Phan, and Jiajie Zhu. Projection algorithms for nonconvex minimization with application to sparse principal component analysis. Journal of Global Oplimization, 65: 657-676, 2016.
Insu Han, Rajesh Jayaram, Amin Karbasi, Vahab Mirrokni, David Woodruff, and Amir Zandich. Hyperattention: Long-context attention in near-linear time. In The Twelfth International Conference on Learning Representations, 2023.
Ying-Lin Hsu, Po- Yu Huang, and Dung-Tsa Chen. Sparse principal component analysis in cancer research. Translational cancer research, 3(3):182, 2014.
Subhodh Kotekal and Chao Gao. Sparsity meets correlation in gaussian sequence model, 2023.
Subhodh Kotekal and Chao Gao. Optimal estimation of the null distribution in large-scale inference, 2024.
Jundong Li. Feature selection datasets at asu. https:///undorgl..iiiee ío/
sci kit meittffaaee --e.turreeeLLDD httpa//ddaaattt_odddddoooo www. html, 2020. Accessed: 2024-09-27.
Yongchun Li and Weijun Xie. Exact and approximation algorithms for sparse pea. arXiv preprint arXiv:2008. /2438, 2020.
Shuangge Ma and Ying Dai. Principal component analysis based methods in bioinformatics studies. Briefings in bioinformatics, 12(6):714-722, 2011.
Malik Magdon-Ismail. Np-hardness and inapproximability of sparse PCA. Information Processing Letters, 126:35-33 38, 2017.
Rahul Mazumder and Trevor Hastie. Exact covariance thresholding into connected components for large-scale graphical lasso. The Journal of Machine Learning Research, 13(1):781-794, 2012.
Stanislav Minsker, Mohamed Ndaoud, and Lang Wang. Robust and tuning- sparse linear regression via square-root slope. arXiv preprint arXiv:22 16808, 2022.
Dimitris Papailiopoulos, Alexandros Dimakis, and Stavros Korokythakis. Sparse pea through low- rank approximations. In International Conference on Muchine Learning, pp. 747-755. PMLR, 2013.
Roman Vershynin. High-dimensional probability: An introduction with applications in data science, volume 47, Cambridge university press, 2018.
Xiwen Wang, Jiaxi Ying, and Daniel Palomar. Learning large-scale mtp .2 gaussian graphical models via bridge-block decomposition. Advances in Neural Information Processing Systems, 36: 73211-73231, 2023.
Xiao- Xiao-Tong Yuan and Tong Zhang. Truncated power method for sparse eigenvalue problems. Journal of Machine Learning Research, 14(4), 2013.
12

---

Published as a conference paper at ICLR 2025
Youwei Zhang and Laurent Ghaoui. Large-scale sparse principal component analysis with application to text data. Advances in Neural Information Processing Systems, 24, 2011.
Xiaowei Zhuang, Zhengshi Yang, and Dietmar Cordes, A technical review of canonical correlation analysis for neuroscience applications. Human brain mapping, 41(13):3807-3833, 2020.
Hui Zou and Lingzhou Xue. selective overview of sparse principal component analysis. Proceedings of the IEEE, 106(8):1111-3300 2018.
Hui Zou, Trevor Hastie, and Robert Tibshirani. Sparse principal component analysis. Journal of computational and graphical statistics, 15(2):265-285 286. 2006.
13

---

Published as a conference paper at ICLR 2025
A ADDITIONAL RELATED WORK
In this section, we discuss additional related work. There exist many attempts in solving Sparse PCA in the literature, and we outlined only selective overview of mainstream methods. Algorithms require worst-case exponential runtimes, including mixed-integer semidefinite programs (Bertsimas et al. 2022 Li & Xie 2020; Cory-Wright & Pauphilet 2022). mixed-integer quadratic programs (Dey ct al. 2022b), non non-convex programs (Zu et al. 2006 Dey et al.] 2022a), and sub-cxponential algorithms (Ding al. 2023). On the other hand, there also exist number of algorithms that take polynomial runtime to approximately solving Sparse PCA. To the best of our knowledge, the state-of-the-ar approximation algorithm is proposed by/Chan al. (2015), in which an algorithm with multiplicative factor min (vk, d-1/33 is proposed. There also exists tons of other (approximation) algorithms that have polynomial runtime. including semidefinite relaxations (Chowdhury et al. 2020), greedy algorithms (Li & Xie, 2020). low-rank approximations (Papailiopoulos al. (013)). and fixed-rank Sparse PCA (Del Pia, 2202). Apart from the work we have mentioned, there also exists a streamline of work focusing on efficient estimation of sparse vectors in some specific statistical models. For instance, for a statistical model known as Spiked Covariance model, methods such as semidefinite relaxations (Amini & Wainwright, 2008: d'Orsi et al., 2020) and covariance thresholding (Deshp et al. 2016). Finally, there are also local approximation algorithms designed for Sparse PCA, and iterative algorithms are proposed (Yuan & Zhang 2013 Hager et al. 2016).
The idea of decomposing large-scale optimization problem into smaller sub-problems is widely adopted approach in the optimization literature. For instance, Mazumder & Hastie (2012) leverages the connection between the thresholded input matrix and the support graph of the optimal solution in graphical lasso problem (GL), enabling significant computational speedups by solving sub-problems. Building on this, Fattahi & Sojoudi (2019) improves the computational efficiency of these subroutines by providing closed-form solutions under an acyclic assumption and analyzing the graph's edge structure. Wang al. (203)) addresses GL in Gaussian graphical models through bridge-block decomposition and derives closed-fore solutions for these blocks. In addition, several methods in the literature involve or are similar to matrix bloek-diagonllization For instance, (Fng et al. (2014)
employ block-diagonalizze
for subspace clustering and segmentation problems by mapping data matrices into block diigonna matrices through linear transformations. More recently, |Han al. (2023)
introduce the sortedLSH algorithm to sparsify the softmax matrix, thus reducing the computational time of the attention matrix through random projection-bate remark that in all studies above, the underlying problems of interests, well as the techniques, different from ours, and the techniques developed therein cannot be directly applied (o/SPCA] A more related work (Devijver & Gallopin 2018) explores block-diagonal covariance matrix estimation using hard thresholding techniques akin to ours. However, our framework diverges in two critical aspects: (i) They assume data is drawn from Gaussian distribution, while we do not we show in Section 4.2 that our method can by applied to any data inputs; and (ii) they use maximum likelihood estimation, whereas our Algorithm 5 involves iterative binary search determine the optimal threshold
bucketing of "block indices". We
B DEFERRED PROOFS
In this section, we provide proofs that we defer in the paper.
B.1 PROOFS IN SECTION
Lemma 1. Let A E Rdxd and e > 0. Given inpul (A, e), Algorithm outputs an e-approximation of A, denoted as A, such that lbs( A) int dim(A,e) in time (d2))
Proof of Lemma 7 It clear that Algorithm [ has time complexity Odd)) and it suffices show that the output À of Algorithm with input (A,e) satisfies lbs (A) int dim (A,
Suppose that A is an e-approximation of A, and satisfies lbs (A*) = int dim(A, We intend to
show lbs A < lbs (A*). WLOG we assume that A E B(A, e) is a block-diagonal matrix with
14

---

Published as a conference paper at ICLR 2025
An
A12
A1p
A21
A22
App
blocks A1. A2, A2,..., Ag. i.e., A = =diag A2......... ,,)). Write A =
Ap1
Ap2
App
A1
A12
A1p
A21
Ag2
App
and =
where Asi Asc, and A have the same number of rows and
Appl columns, by the definition of (3(A, e), it is clear that |Aill I for any 1 S i S This further implies that for any s i * S p. Aig is zero matrix by the procedure described in
Apl
Ap2
Algorithm| and therefore clear that lbs A I lbs (A*). By minimality of ^*, we obtain that
lbs lbs (A*).
B.1.1 PROOFS SECTION UPA
Theorem Let A Raxd be symmetric, d, and E > O. Denote by A the output of Algorithm with input EA, E an optimal solution to ISPAA input (A k), and x² optimal solution SPCA with inpul (A, Then, it follows that
1(x') Ax" 5. AT 28. E.
Proof of Theorem
7
Note that
(x*)T Ar
I Az S (x*) Ax 2 A T (At A)T
TA
(x*) Ax" - F.. A
+ k. ||A A floc
where the last inequality follows from Holder's inequality and the fact that Therefore, suffices to show that
|| S V || = Vh.
T Ar
=T S k ||A A Allloo
Observe that
[(=') Axx - -xT
max T Ar
max x -Af
||
| -=
||r sk
< max
x (AA A)xz
|
| Sk
< k: ||A4 A Allooo
where the first inequality follows from the following facts:
max a 111,
x' Ar max =|
T A S max x. T(AA Af)x - Inllosk
| 1
max x T x L 11
max x Ar < max x TTA -Aex Ix |
|
| Sk
|allosk
Finally, combining the fact that ||A A E, we are done.
Theorem 2. Let À diag A Àz.... Ap) be a symmetric block diakonaa matrix. Denote OPT to be the optimal value to SPCA with input pair (Ā, k). Let OPT; ; to be the optimal value to SPCA with input pair (As, k), for i [p]. Then, one has OPT max OPT;
15

---

Published as a conference paper at ICLR 2025
Proof of Theorem OPT & maxie(p] O
It is clear that OPT 2 ma OPT, It remains to show the reverse, i.e.,
Suppose not, have OPT max OPT, Then, there must exist principal submatrices A Aly,..., A' contained in A1, A Ap such that diag( A1, A2....,, A²p) is a & x matrix, and
OPT max diag(A Ag. A/)xr max OPT; |
sE [p
This gives a contradiction,
max OPT;
max max x. T A,, :Elp] |
max x TT diag( A1, A22 ......... Apya, =1
P
where the inequality follows from the fact that
max TT diag(A', A2..... A²)== max max x T 1/== max max x A/. = max OPT,
|=||==1
#Elpl
1,1
selpl |
#Elpl
Theorem 3. Let A gdxd be symmetric and positive semidefinite, let k be positive integer, and denote by x² E Rd an optimal solution to SPCA with input (A, Suppose that an algorithm for SPCA with input (A. k)) finds an approximate solution x E R2d to SPCA with multiplicative factor mik, d) 2 1 and additive error alk, d) 2 O, i.e. one has x Ar 2 (x*) Axx /m(k, d) - a(k, d). Furthermore, we ussume that mkk, dl) and a(k, d) is non-decreusing with respect to For e > denote by E T& be the ouiput of Algorithm with input tuple (A, e, A). Then, one has
( Ax y' TTAy Ay 2 m. (hi, int dim (A,e)) (k, int dimndi,nn,.. (ii)(A m (k, int (A,)) ke.
Proof of Theorem 3 In the proof, we use the same notation as in Algorithm 3 We assume WLOG that the thresholded matrix A is a block -iooo-da matrix, i.e., A diag(A), A5. A5). where
Aç = A5 We also define À == diag (Āi, Āg..,, Ap) with A, = As..s. It is clear that
A 7 < Denote by F- R an optimal solution to SPCA with input (Ā, k). By Theorem 2
one can assume WLOG that E* has zero entries for indices greater than d1 == |S1l. In other words, K is "found" in the block A₁. Recall that we denote x, the solution found by A with input (Ā,, k) Then, is clear that
max a. A,z, 2 2I A1
sE(p!
(7*) Ax
Ar
a(k,dd)
m(k d1)
(7*) AF
a(k, int dim(A, (A))))
mkk int dim (A,
Recall that y, the reconstructed solution. By the proof of Theorem [ one obtains that
max V.
Ay:
AE|p
max sE(p]
2. ApT:
(1)
(3*) (T*)TA"" Ar -a(k, int dim dim(A,
m(k, int dim dim(A, A,
1
((x*) Axt xa|1A al_ a a(k, int dim(A, ())
mk, int dim (A, =))
16

---

Published as a conference paper at ICLR 2025
B.1.2 PROOFS IN SECTION 3.2.2
Proposition Let E Raxdd and let i. be a positive integer. Suppose that for specific algorithm
A. the time complexity of using for (approximately) solving SonegsSCAAWWWr with input (A, k) is function g(k, d) which is convex and non-decreasing with respect to We assume that g(k, 1) = Then, for given threshold & > O, the runlime of Algorithm with input tuples (A, k, A) is at most
o ( dim((, ([intcim(A)))
9 (ks int (im((A de
Proof of Proposition Suppose that the largest block size of Á B(A,E) is equal to int dim (A, e). and A can be sorted a block-diagonal matrix with p blocks A1, A2......, Ap and block sizes d₁, ... dp. respectively. WLOG we assume that p < d, otherwise the statement holds trivially. It clear that the time complexity of running line 5 in Algorithm 3is upper bounded by a constant multiple of
Σ g(k, dx).
15d, <int eim)),e) Em dend
Since & is fixed, the optimization problem polytope
is maximizing convex objective function over
P P== P dd.. RP 15 d; S int dim,A,)), Vi E [pl,
Σ d: 1=1
Therefore, the optimum to 2 is obtained at extreme point of P. It is clear that extreme point (di, dy)....()) of P is active at a linearly independent system, i.c., among all di ss there are p of them must be equal to or int dim Lim(A,9 A) and the last one is taken such that the euqality D di holds. Since there could be most [d/ int dim(A, e) many blocks in À with block size equal to int dim(A, (A,)). E). and that gk: d) increasing with respect to d, we obtain that the optimum of the optimization problem 2 upper bounded by
(k, int dim(A,s)) + (p- (p-
d (
int dim A, 8)
int dimn (A,
d
9 (k, int dimmaii))))))
int dim A,
We are then done by noticing the fact that lines 3 in Algorithm have a runtime bounded by (dd).
B.2 PROOFS IN SECTION
To establish a connection with robust statisties, we begin by expressing each element of the input matrix A in Model aa follows:
Ag Aij Eij Aij o i
where Zij is an i.i.d. centered (e²/o2)-sub -/bo--suu-Gaaasiia Gaussian random variable with unit variance. Considering that Á is block-diagonal and, according to (ii) in Model has at most [d/d"] (d*)² +a) nonzero entries, the methodology proposed by Comminges et al. (2021) can be employed to estimate σ² accurately using Algorithm
Inspired by the methods proposed in Comminges et al. (2021), we show that the number ō found in Algorithm 4 could provide constant approximation ratio to the true variance σ² in Model with high probability:
Proposition 3 (adapted from Proposition in Deshp et al. (0166). Consider Model and denote Pz to be the distribution of the random variable Fry/e 0 Let 2 be the number found in
17

---

Published as a conference paper at ICLR 2025
Algorithmn with input tuple (A, C, a.). There exist absolute constant dimension of the input matrix A satisfying
0 such that for the
log (8C + +8)
+1
log d > log (4C +4
/11 /(-aa
20 +11
we have that for distribution Pz of Ziz. the variance o², one has uniform lower bound on the accuracy of estimation:
( ō²
2 1 exp {- da+1
inf inf
Pz a20
inf losalta
P Przaau
2 σ²
S 2
{ 40 +4
Proof of Proposition Recall that we write
A.j = Aij + E Aij + Zij, is
where Zy is i.i.d. centered (d%/o")-sub-Gaussin
random variable with unit variance. Since d' & Cd for some E (0, 1), WLOG we assume that d* = Cdd for some C > 0, and for simplicity we assume that to and C all integers to avoid countless discussions of integrality, By Algorithm 4 one divides the index set ((i,j) E [d] [d] i S j} into = (2C + 2) disjoint subsets B1, B2...., Bm, and each B, has cardinality at lease k == [dd + d)/ (2m)J. Denote the set J to be the set of indices such that the set B; contains only indices such that for any (k, () B..
Ast 0. It is clear that [J] is lower bounded by Cd (2C 1)d²+a
By Bernstein's inequality (see, e.g., Corollary 2.8.3 Vershynin 2018 ), is clear that there exists an absolute constant 0 such that
1
P
Ba
Rdl
0²
=
P
Ba
kt
2
R:
(A:1)E
0 0²
< exp
-cmin
BB
40 20²
0 0²
de + d
}
< 2exp
-emin
4g
20²
Then, we denote 1 == 2/ denote A; to be the random event and denote A, be the indicator function of Aç. One obtains that
Abs 1).
[1.]
P P(0²
1)
P
1
N
m
2
1=1
=
A.
> m
Σ A
2
iEJ
it
< P
1 2
m
Cd
2
VEJ
Denote
00 o²
d2 +
f =
-Cd++1 -2|J exp
-emin
40 20"
2m
18

---

Published as a conference paper at ICLR 2025
and by Hoeffding's inequality, obtains that
P
Σ
2 m
Cd
2
VEJ
( 0 o²
)
de -
}
P
Σ A, BIA.) 2 2
Cd -// exp -emin
40 20²
2m
(EJ
Σ (LA, E1A)22
VEJ
$ exp 26
Since (2C + 2) it clear that | (2C' + 1)d and
exp exp { {-omin
o²
411 2d (1100++))2001111 |d
t t d
11- 2(2C +
ld' [d'
4
20²
where the last inequality follows from the fact that
log (SC +
log > log (4C 4)
((1-a).
og2 20
0²
+ +1]
20
Hence, the probability of x² & is upper bounded by exp( -d+t 1/C= + 4)}.
Proposition 2. Consider Model 7 and denote by the output of Algorithm with input tuple (A, (A,C,a, u.. Let k S d be positive integer, and assume that d satisfies that d' >> Co (C-(C+1) log(8C 8) for some large enough absolute constant Co 0. Then, the following holds with probubility at least 1 2d-1 exp 22gg2log - di+a//(a /(AC 4)}:
(i) Denote by
X= the optimal solution to SPCA with impul (Ā, k). For (approximation) algorithm A with multiplicatiie factor m(k, 6,, d) 2 and additive error a(k, d) 2 0, where the functions and a is non non-deceasing with respect to the output y of Algorithm 3 with input tuple (A, k, A) satisfies that
T Ar
Ay 2 m (k, (k,d') a (k, d")
m (k, (k, Kg
(ii) in addition, the time complexity of using for (approximately) solving SPCA with input (A, k) is a function g(k, d) which is convex and non-decreasing with respect to d, and satisfies g 1) 1. Then, the runtime of Algorithm with input tuples (A, k, E, A) most
o
((
d
g (ki d') + d2
Proof of Proposition 2 We first show that, the probability that |bll S E is at least - 2d 2d-1 exp 2 log d I0+1 This comes directly from Proposition 3 that 8² satisfies that 1/2 - aa with probability at least - -x +4), the definition of Ext being i.i.d. g-sub-Gaussian variables, and union bound argument, Indeed, define the random event A == & a ) one
19

---

Published as a conference paper at ICLR 2025
see that
PP ([ E) PP ([EM)] > E A) P ([ A())(
(IEMl
2 2 3uō log d A) P(A) + P ( > E P9))(AA)
[Eaal
VBua log
P(A) + PPA4)
+ +exp {
d'
S (EEul
VBuo Vlog d)
4C+
6u2o²
+ texp {
d
< 2exp
log d
{
20²
40 +4
d=ta
2d + exp
40 +4
By union bound, one has
d
P[ > &f > >E 2d + exp log
4C +
This implies that, with high probability, int dim (A, S dt In the remainder of the proof, we will assume that the random event |ill I E holds true.
We are now ready to show the desired approximation error. Denote by x² the optimal solution to SPCA with input (A, k), by the proof of Theorem is clear that
Ax
(T kē
due to the fact that |
A A
5. Combining with Theorem we have
y Ay 2 v Ay E
(x*) Ar
(k, int dim(A,))) - -(1+ (k, int dim(m,))))))))))))))))))))))))))
Ke
(k, int (im((,))))
(7*)* AF
(1+
m (k, int timm(A,)) a (k, int dim (A,,))
KE.
m (k, int dim (A,
Finally, we obtain the desired runtime complexity via Proposition
B.3 PROOFS IN SECTION 4.2
Theorem 4. Let Raxd be symmetric and positive semidefinite, and let le be positive integer.
Suppose that an algorithm A for SSCCA with input (A, k) finds an approximate solution x E RRd
SPCA such that has maltiplicative factor m(k,d) d) 2 and additive error a(k, > 0, with mik, d) and a(k, d) being non non-decreasing with respect to d. Suppose that Algorithm with input (A, 8, A, do) terminates with OPT(=*) 2 Then one has
OPT(E*) 2 m(k, do) OPT(O) - a(k, do) m(k, do) k=*
Suppose that for A, the time complexity of using for (approximately) solving SPCA with input (A, k) is function g(k, d' which is convex and non-decreasing with respect d, g 1, and g(k 0) = 0. The runtime for Algorithm is at most
|A 8
(T
[ do
g (k, do) de )
o log
Proof of Theorem
We will use Theorem 3 and Proposition
to prove this theorem.
20

---

Published as a conference paper at ICLR 2025
We first prove the approximation bound. First note that | - A² 11 = By Theorem
OPT(E*) 2 (x)) Ar
m (k, int dim (Am((A a (k. int (elddimA,,, (A,,*)) (k, int (im(A,
ke*
OPT(O)
171 (k, int dim(A,,*)) a (k: int dim(A,
m (k, int dim(A, s*)) ks".
Note that in Algorithm 5 the algorithm is guaranteed to terminate with an et such that int dim(A, S do. Since the functions m(k, d) and a(k, d) are all non-decreasing with respect to d, we have that m (k, int dim (A,*)) & m(k, do) and a (k, int dim(A, S a(k, do), and therefore we obtain our desired approximation bound
OPT(0) a (ki do) m (k: do)
OPT(E*)
ke*
m (ki do)
Here, the inequality holds due to the fact that (i) OPT I ket, then since OPT(e*) is non- negative, the inequality trivially holds; (ii) if (P)(0) > ke* then the inequality also holds since a(k, d) and m(k, d) are non-decreasing with respect to d.
Then, we show the time complexity. By Proposition T it is clear that the runtime for a single iteration in the while loop in Algorithm S is upper bounded by
o for some 2 0 such that int dim(A, < do. Let E 2 E2 be two positive numbers, and it is clear that int dim(A, E1) I int dim E2). Since
d ( int (A,E)
9 (ki, int dim(A,E)) + d²)
(3)
int dim E1)
int dim(A, 0,
int dindiii,,,a s1) = int dim limm(,s5) 82)
int dim 52) + 1
int dim (, 92)
and by convexity of g(k; d), along with the fact that g(:, 0) = 0, we obtain that g(k, int dim (iim(A,s)))
< int dim (1) E1)
int dim (A, g(k, k, int dim A, E2)) + O.
This implies that for any 2 0 such that int dim(A, < dg. we have that
d
g(,, int dim(A, < d
g(k, do),
int (A,
dg
and thus
d int dim A, $)
g(k, int dim(A.=)
int dim (A, g(k: int (imm(A,)) E)) + g(k, int dim(A, (A,s))
dg (k:, do) + g(k, int (Aim(A,))))
d
d
g(k: do) + g(k, do)
g(k, do).
do
do
Combining with (3). and the fact that at most (log(( (1ll./8)) iterations would be executed in the while loop in Algorithm 5) we are done.
C ADDITIONAL EMPIRICAL SETTINGS AND RESULTS
In this section, report additional empirical settings in Appendix and detailed empirical
results in Appendices C.2 and C.3 We report empirical results under M Model in Appendix CC.] In
Appendix C.5 we discuss impacts of parameters do and 5 in Algorithm
21

---

Published as a conference paper at ICLR 2025
C.I EMPIRICAL SETTINGS
Datasets, We summarize the dimensions and percentage of non-zero entries in the input matrices in T it important note that more than 70% of these input covariance matrices are dense. We utilize these matrices to evaluate the practicality of our approach, demonstrating its effectiveness irrespective of whether the input matrix has a block-diagonal structure.
Table 5: Dimensions for each dataset and percentage of non-zero entries.
Dataset
Dimension (d)
Percentage of non-zero entries
CovColonCov
500
100%
LymphomaCovl
500
100%
Reddit1500
1500
5.30%
Reddit2000
2000
6,20%
LeukemiaCov
3571
100%
yymphmmCov2
4026
100%
ProstateCov
6033
100%
ArceneCov
10000
100%
GLI85Cov
22283
100%
GLABRAI8OCAS
49151
100%
DorotheaCov
100000
77.65%
Parameters. (1) When integrated Algorithm with
with Branch-and- and-anudddBdddddd- algorithm: We choose the == || loc+ dg =30, and 0.01 ||
parameter ==
(ii) When integrated
Algorithm 5 with Chan's 's algorithm: We choose the parameter E == 1, a = 0, b : |||| do == 2k
(twice the sparsity constant), and 0.01 A|
Early Stopping. For experiments conducted with positive semidefinite input matrix A, we add very simple step Algorithm [ we set an extra stopping criteria on Algorithm [ which breaks the while loop as soon as problem with largest block size do has been solved. This additional step helps further reduce computational time Algorithm
Compute resources. We conducted all tests computing cluster equipped with 36 Cores (2x 3.1G Xeon Gold 6254 CPUs) and 768 GB of memory.
C.2 EMPIRICAL RESULTS WHEN INTEGRATED WITH BRANCH-AND-BOURT ALGORITHM
In this section, we provide the detailed empirical results of our framework when integrated with Branch-and -Bund algorithm (previously summarized in Table 3, in Table We provide statistics related to the best solution in Table 7 including the ratio of corresponding =* (the best threshold found) and ||Al| the percentage zeros in the matrix A= and the Jaccard index between the solution obtained by baseline algorithm and that obtained by Algorithm 5 which is defined as
Jaccard index := suupp(znase O supp(rouus)
(upp(eRRase supp(rooar)
Here, .Base denotes the solution obtained by the baseline algorithm, which refers to Branch-and Bound algorithm in this section and Chan's 's algorithm in Appendix [C.3 and ZOurs denotes the solution obtained by Algorithm when integrated with the baseline algorithm. We also summarize the average approximation errors and average speedup factors for with respect to in Table
From the data presented in T it is evident that framework achieves exceptional performance metrics. The median approximation error is O, with the 90th percentile approximation error not exceeding 19%. Additionally, the median speedup factor is recorded 20.87. These statistics further underscore the efficiency and effectiveness of our framework, demonstrating its capability to provide substantial computational speedups while maintaining small approximation errors.
From the data in Table we conclude that the Jaccard index is oftentimes one, and the change of Jaccard index aligns with the change of the approximation error. Additionally, for each dataset, the relative threshold value e*/ ||44] and the percentage of zeros in A remain relatively stable across
22

---

Published as a conference paper at ICLR 2025
Table Comparison of runtime and objective values on real-world datasets, between Branch-and- Bound and Algorithm integrated with Branch-and-Boumd problem, and k is the sparsity constant. We set the time limit to 3600 seconds for each method. Approximation emrors between two methods, and speedup factors are reported. We note that the objective value ArceneCov scaled by 10
is the dimension of Sparse PCA
Dataset
d
Branch-anc and Bod-BBundd
Ours
Error (%)
Speedup
Time
Objective
Time(s)
Objective
CovColonCov
500
8
1059.49
1.20
1059.49
0.00
6.67
CovColonCov
500
8
1646.45
1.08
1646.45
0.00
7.41
CovColonCov
500
10
3600
2641.23
42.93
2641.23
0.00
83.86
CovColonCov
500
15
3600
3496.6
666.95
3496.6
0.00
5.40
LymphomaCov
500
11
2701.61
1.32
2701.61
0.00
8.33
LymphomaCov
500
44
4300.5
1.51
4300.5
0.00
29.14
LymphomaCov
500
10
3600
6008.74
85.72
6008.74
0.00
42.00
LymphomaCovl
500
3600
7628.66
324.6
7628.66
0.00
11.09
Reddit1550
1500
3
75
946.12
3.68
946.12
0.00
20.38
Reddit1
1500
3600
980.97
6.86
980.97
0.00
524.78
Reddit1500
1500
10
3600
1045.74
9.20
1045.74
0.00
391.30
Reddit1500
1500
15
3600
1082.12
6.87
1082.12
0.00
524.02
Reddit2000
2000
412
1311.36
8.70
1311.36
0.00
47.36
Reddit2000
2000
3600
1397.36
9.03
1397.36
0.00
398.67
Reddit2000
2000
10
3600
1523.82
13.18
1523.82
0.00
273.14
Reddit2000
2000
15
3600
1605.48
9.10
1601.32
0.26
395.60
LeukemiaCov
3571
280
7.13
18.89
7.13
0.00
14.82
LeukemiaCov
3571
5
3600
10.29
18.64
10.29
0.00
193.13
LeukemiaCov
3571
3600
17.22
64.03
17.22
0.00
56.22
LeukemiaCov
3571
15
3600
21.87
168.47
21.87
0.00
21.37
LymphomaCov2
4026
210
40.62
25.73
40.62
0.00
8.16
LymphomaCov2
4026
5
144
63.66
26.75
63.66
0.00
5.38
LymphomaCov2
4026
10
3600
78.29
293.42
69.27
11.52
12.27
LymphomaCov2
4026
15
3600
93.46
1810.19
86.20
7.77
1.99
ProstateCov
6033
102
8.19
57.35
8.19
0.00
1.78
ProstateCov
6033
5
3600
12.92
54.63
12.92
0.00
65.90
ProstateCov
6033
10
3600
24.38
407.14
24.38
0.00
8.84
ProstateCov
6033
15
3600
34.98
1175.40
34.98
0.00
3.06
ArceneCov
10000
88
3.36
138.07
3.36
0.00
0.64
ArceneCov
10000
248
5,56
135.47
5.56
0.00
1.83
ArceneCov
10000
10
3600
10.81
138.56
10.81
0.00
25.98
ArceneCov
10000
3600
15.30
140.92
15.30
0.00
25.55
k. as the largest block size dg 30, and the best solutions are typically found in blocks near this size.
From the data in Table 3 it clear that our framework consistently obtains significant computational speedups across a range of settings. Notably, achieves an average speedup factor of 13.52 for & = 3. with this factor increasing dramatically to over 153 for k 2 the relationship between & and performance metrics: although the speedup factor increases with larger k there is corresponding rise the average approximation errors. Specifically, for and & 5, the average errors remain exceptionally low, which are all zero, but they escalate 1.47% for & = 10 and 1.05% for I = 15. This indicates a consistent pattern of increased approximation errors as k grows, which is aligned with Theorem despite the significant gains in speed.
observation from the analysis
C.3 EMPIRICAL RESULTS WHEN INTEGRATED WITH CHAN'S ALGORITHM
In this section, we provide the detailed numerical test results of our framework when integrated with Chan's algorithm (previously summarized in Table in Table We provide statisties related the
23

---

Published as a conference paper at ICLR 2025
Table Additional statistics regarding the best solution obtained by Algorithm in Table 6 Denote by e* the threshold found that yields the best solution in Algorithm S A the input covariance matrix, A is thresholded matrix obtained by Algorithm with input (A,**).
Dataset
Error (%)
E*/ //||||l
Zeros in A (%) Jaccard Index
CovColonCov
0
0.75
99.988
1.00
CovColonCov
0
0,59
99.94
1.00
CovColonCov
10
0
0.58
99,926
1.00
CovColonCov
15
0
0,59
99.941
1.00
LymphomaCovl
0
0.50
99,982
1.00
LymphomaCov
0
0.50
99.982
1.00
LymphomaCov
0
0.38
99,951
1.00
LymphomaCov
15
0
0.35
99.940
1.00
Reddit1500
0
0.13
99.999
1.00
Reddit1500
0
0.06
99.997
1.00
Reddit 150
10
0
0.06
99.997
1.00
Reddit1500
15
0
0.06
99.997
1.00
Reddir2000
0
0.13
99.999
1.00
Reddit2000
0
0.09
99.998
1.00
Reddic2000
10
0
0.09
99.998
1.00
Reddit2000
0.26
0.09
99.998
0.88
LeukemiaCov
3
0
0.41
99.998
1.00
LeukemiaCov
0.50
99,999
1.00
LeukemiaCov
10
0
0.50
99.999
1.00
CoukemiaCov
0
0.44
99.999
1.00
LymphomaCov2
3
0
0.50
99.999
1.00
LymphomaCovv
0.50
99.999
1.00
LymphomaCov2
10
11.52
0.42
99.999
0.42
LymphomaCov?
15
7.77
0.44
99.999
0.30
ProstateCov
0
0.75
99.999
1.00
ProstateCov
0
0.68
99.999
1.00
ProstateCov
10
0
0.69
99,999
1.00
ProstateCov
15
0
0.69
99.999
1.00
ArceneCov
0.75
99.999
1.00
ArceneCov
0
0.75
99.999
1.00
ArceneCov
0.75
99.999
1.00
ArceneCov 15
0
0.75
99,999
1.00
Table 8: Summary of average approximation errors and average speedup factors for each ks when integrated our framework with Branch-and- Bounnd algorithm. The standard deviations are reported in the parentheses.
k
Avg. Error (Std. Dev)
Avg. Speedup (Std. Dev)
3
0.00 (0.00)
13.52 (15.12)
5
0.00 (0.00)
153.28 (203.17)
10
1.47 (4.06)
111.70 (141.81)
15
1.05 (2.72)
123.51 (210.55)
best solution in Table 10 including the ratio of corresponding =* and |Allsc the percentage of zeros the matrix A=", and the Jaccard index. We also summarized the average approximation errors and average speedup factors for with respect to l: in Table 11
In Table 9 we once again demonstrate the substantial speedup improvements our framework achieves when integrated with Chan's 's algorithm. The results show median approximation error of -0.38% with the maximum approximation enror not exceeding 0.28%, and a median speedup factor of 3.96. We observe that the speedup factor grows (linearly) as the dimension of the dataset goes up,
24

---

Published as a conference paper at ICLR 2025
highlighting the efficacy of our framework, especially in large datasets. Notably, in all instances within the DorotheaCov dataset, the instances = 500, 1000, 2000 in GLI85Cov and ELABRAI80COM datasets, the instances & 1000, 2000 in ArceneCov dataset, our framework attains solutions with negative approximation errors, indicating that it consistently finds better solutions than Chan's algorithm alone. These instances provide tangible proof of our framework's potential to find better solutions, as highlighted in Remark
In Table IO, the Jaccard index does not vary significantly with k, nor does its change align with the change of the approximation error. A possible reason is that Chan's algorithm an approximation algorithm, and the Jaceard index may not fully capture the similarity between the solutions obtained and the true optimal solutions. As expected, since do 2k, the relative threshold value =*/ /| || and the percentage of zeros decrease k increases.
In Table 11] we do not observe an increase in speedup factors as k increases. This phenomenon can be attributed the runtime of Chan's 's algorithm, which is (dd) and independent k. Consequently, as predicted by Proposition] [ the speedup factors do not escalate with larger values of k. Additionally, unlike the trends noted in T the average errors do not increase with k. This deviation likely due to Chan's algorithm only providing approximate solution with a certain multiplicative factor: As shown Theorem 3. our algorithm has the capability to improve this multiplicative factor at the expense of some additive error. When both & and d are large, the benefits from the improved multiplicative factor tend to balance out the impact of the additive error. Consequently, our framework consistently achieves high-quality solutions that comparable to those obtained solely through Chan's algorithm, maintaining an average error of less than 0.55% for each k.
Finally, note that Chan's algorithm is considerably more scalable than the Branch-and-Bonl Bound algorithm. The statistics underscore our framework's ability not only to enhance the speed of scalable algorithms like Chan's ss but also maintain impressively low approximation errors in the process.
Table 9: Comparison of runtime and objective values on real-world datasets, between Chan's algorithm and Algorithm integrate with Chan's algorithm. d is the dimension of Sparse PCA problem, and & is the sparsity constant. We note that the objective value for GLI85Cov is scaled by 10 and that for GLABRA B.ABBRA1OCOv is scaled by 10 The error being negative means that our framework finds a better solution than the standalone algorithm.
Dataset
d
k
Chan's algorithm
Ours
Error
Speedup
Time
Objective
Time
Objective
ArceneCov
10000
200
228
91.7
288
91.7
0.00
0.79
ArceneCov
10000
500
239
140.4
188
140.4
0.00
1.27
ArceneCov
10000
1000
249
190.4
349
190.8
-0.21
0.71
ArceneCov
10000
2000
256
233.3
377
233.5
-0.09
0.68
GLI85Cov
22283
200
2836
35.8
1660
35.7
0.28
1.71
GLI85Cov
22283
500
2728
38.8
1528
39.4
-1.55
1.79
GLI85Cov
22283
1000
2869
39.9
1942
40.9
-2.51
1.48
GLI85Cov
22283
2000
2680
42.3
2034
42.5
-0.47
1.32
G00B
49151
200
35564
44.7
4890
44.6
0.11
7.27
ILABRA180Com
49151
500
33225
60.4
4475
60.7
-0.50
7.42
GLABRAI 180Cov
49151
1000
30864
70.5
4142
70.7
-0.28
7.45
GLABRAI80CA
49151
2000
29675
78.1
4830
78.2
-0.10
6.14
Dorothea
100000
200
246348
3.25
16194
3.34
-2.77
15.21
Dorothea
100000
500
246539
4.94
17554
5.07
-2.63
14.04
Dorothea
100000
1000
244196
6,23
16934
6.38
-2.41
14.42
Dorothea 100000
2000
255281
7.18
17774
7.28
-1.39
14.36
C.4 EMPIRICAL RESULTS UNDER MODEL
In this section, report the numerical results in Model Algorithm with threshold obtained by Algorithm 4 In Model I we set to have i.i.d. centered Gaussian variables with standard
deviation 0.1 in its lower triangle entries, i.e., Fij for 1 < S and set F.j- Fgs for to
25

---

Dataset Published as a conference paper at ICLR 2025
Table 10: Additional statistics regarding the best solution obtained by Algorithm 5 in Table Denote by e* the threshold found that yields the best solution in Algorithm S A the input covariance matrix, A= is thresholded matrix obtained by Algorithm with input (A, e
k
Error (%)
<*/ /||||||
Zeros in A (%)
Jaccard Index
ArceneCov
200
0.00
0.38
99.974
0.82
ArceneCov
500
0.00
0.26
99.868
0.82
ArceneCov
1000
-0.21
0.16
99.469
0.79
ArceneCov
2000
-0.09
0.05
95.891
0,94
GLI85Cov
200
0.28
0.09
99.997
0.67
GLI85Cov
500
-1.55
0.04
99.985
0.49
GLI85Cov
1000
-2.51
0.02
99.962
0.62
GLI85Cov
2000
-0.47
0.01
99.775
0.73
GLABRAI 80Cov
200
0.11
0.18
99.999
0.82
GLABRAI 80Cov
500
-0.50
0.10
99.997
0.81
GLABRA 10CO
1000
-0.28
0.05
99.990
0.82
GLABRAI 1000 180Cov
2000
-0.10
0.04
99.980
0.80
DorotheaCov
200
-2.77
0.12
99.999
0.78
DorotheaCov
500
-2.63
0.09
99.999
0.79
DorotheaCov
1000
-2.41
0,06
99,999
0.76
DorotheaCov
2000
-1.39
0.05
99.998
0.67
Table Summary of average approximation errors and average speedup factors for each k, when integrated our framwork with Chan's algorithm. The standard deviations are reported in the parenthe- ses.
k
Avg. Error (Std. Dev)
Avg. Speedup (Std. Dev)
200
-0.60 (1.45)
6.25 (6.63)
500
-1.17 (1.17)
6.13 (5.96)
1000
-1.35 (1.28)
6.02 (6.36)
2000
-0.51 (0.61)
5.63 (6.31)
make it symmetric. We generate 30 independent random blocks A, each of size 20, with each block defined as M,T M/100, where M E TR 100x20 has i.i.d. standard Gaussian entries (which implies that u in MT
We Algorithm with C 1, 0.7, and 1, obtaining the threshold È the output. Subsequently, we execute Algorithm 3 with A, k E 3,, 10), the Baunchandd-oumt algorithm. and ह We denote the solution output by Algorithm 3 as .rOurs
In Table 12 we compare Algorithm 3 integrated with Branch-and unndBound algorithm, with vanilla Branch-and Maacchandd-Bl Bondd algorithm. We report the optimality gap, speed up factor, and the value of g outputted by Algorithm 4 The optimality gap defined as:
Gap = ObjEs Objours:
where Obj := - AugB, and ΣBB is the output of the Baaechhand-Bound algorithm with input (Ā, k),, and where Objoun is the objective value y Ay in Proposition 2
To ensure reproducibility, we set the random seed 42 and run the experiments ten times for each
observe that the average optimality gap increases & grows. Additionally, does not depend on at all. The optimality gap is much smaller than the predicted bound 3k E, providing computational verification of the bound proposed in Proposition] 2 The speedup factor is exceptionally high, often exceeding thousand when 5.
From Table I2
remains relatively stable across different values of k as the calculation
26

---

Published as a conference paper at ICLR 2025
Table 12: Summary of average optimality gaps, average speedup factors, and average threshold for each kt when integrated our framework with Branch-and-Boume algorithm. The standard deviations are reported in the parentheses. The time limit for Branch-and aaaechand-Boundd is set 600 seconds.
k
Avg. Gap (Std. Dev)
Avg. Spdup (Std. Dev)
Avg. (Std. Dev)
2
0.29 (0.08)
16.52 12)
0.96 (0.0016)
3
0.35 (0.13)
94.78 (78.52)
0,96 (0.0030)
5
0.48 (0.12)
2010.02 (625.55)
0,96 (0.0028)
7
0.59 (0.14)
1629.38 (869.91)
0.96 (0.0040)
10
0.64 (0.14)
1900.05 14)
0.96 (0.0036)
C.5 IMPACTS PARAMETERS do AND å IN ALGORITHME
In this section, we discuss the impacts of parameters do and 6 in Algorithm We showease the impacts by comparing the approximation errors and speedups on certain datasets when integrating Algorithm with Branch-and and-Bound algorithm.
C.5.1
IMPACT OF do
we summarize the numerical results when increasing do from 30 to 45 and 55 in the LymphomaCov2 dataset. We observe significant improvement on the approximation error, going from an average of 4.82% to 1.94% and 0.03%, at the cost of lowering the average speedup factor from 6.95 to around 4.
In Table 13
Table 13: Comparison of errors and speedups for different do values across LymphomaCov2 dataset.
Dataset
d
k
dp
30
dg = 45
dg
55
Error
Speedup
Error
Speedup
Error
Speedup
LymphomaCov2
4026
0.0
8.16
0.0
8.29
0,0
8.62
LymphomaCov3
4026
0.0
5.38
0.0
5.35
0.0
5.95
LymphomaCov4
4026
10
11.52
12.27
0.0
1.84
0.0
1.00
LymphomaCov5
4026
15
7.77
1.99
7.77
0.92
0.13
1.00
Overall
4.82
6.95
1.94
4.10
0.03
4.14
C.5.2
IMPACT OF 8
In Table 14 we summarize the numerical results when increasing 8 from 0.01 100-|140 ||A4]
0.1 ||A
and 0.15 || in the ProstateCov dataset. We observe that, by increasing 8, the speedup factors dramatically improve, from an average of 19.90 38.41 and 99.49, at the potential cost of lowering the approximation error from 0 to 6.07%.
Table 14: Comparison of errors and speedups for different 8 values across ProstateCov dataset.
Dataset
d
8 = 0.01 ||A||
8
0.1 ||A|||
0.15
||]]
Error
Speedup
Error
Speedup
Error
Speedup
ProstateCov
6033
0.0
1.78
0.0
2.79
0.0
3.62
ProstateCov
6033
0.0
65.90
0.0
106.32
0.0
128.16
ProstateCov
6033
10
0.0
8.84
0.0
29.56
0.0
128.25
ProstateCov
6033
15
0.0
3,06
0.0
14.99
24.27
137.93
Overall
0.0
19.90
0.0
38.41
6.07
99.49
27

---

Published as a conference paper at ICLR 2025
D EXTENSION TO NON-POSITIVE SEMIDEFINITE INPUT MATRIX
In this paper, we focus on solving SPCA where the input matrix is symmetric and positive semidefinite, following conventions the Sparse PCA literature (Chan et al., 2015 Amini & Wainwright 2008 Berk & Bertsimas 2019 Dey et al. 2022a b). We remark that, our framework remains effective for symmetric and non-positive semidefiite imputs, provided that the algorithm A for (approximately) solving SPCA supports symmetric and non-positive semidefinite inputs [ In other words, if A compatible with symmetric and non-positive semidefinite inputs, the assumption of A being positive semidefinite can be removed in Theorems and 4 and in Proposition 3 (stated in Model T. and all theoretical guarantees still hold.
"Such algorithms do exist; for example, naive exact algorithm that finds eigenvector corresponding to the largest eigenvalue of submatrices Ag,s for all subsets with |S| < k.
28


Reviews:
[{"review_id": "bq7R2mzZLb", "reviewer": "Reviewer_rGCR", "summary": "The paper presents a novel framework for approximating Sparse PCA by decomposing the original large-scale problem into smaller subproblems through matrix block-diagonalization. The method involves three key steps: (1) transforming the input covariance matrix into a block-diagonal form by thresholding and grouping non-zero entries, (2) applying any existing Sparse PCA algorithm to each block independently, and (3) reconstructing the overall solution from the subproblem results. This approach aims to achieve significant computational speedups while maintaining high solution quality, supported by theoretical guarantees and extensive empirical evaluations on various large-scale datasets.", "questions": "Please address the concerns outlined in Weaknesses.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The framework significantly reduces runtime by decomposing the original problem into smaller, more manageable subproblems.\n\n2. The paper provides approximation guarantees and time complexity analyses, ensuring that the method preserves solution quality.\n\n3. Extensive experiments on diverse, large-scale datasets demonstrate the framework’s effectiveness in reducing computational time with minimal approximation errors.", "weaknesses": "1. The core idea of decomposing a large-scale optimization problem into smaller subproblems and then combining their solutions has been previously explored in the literature (e.g., [1,2,3]). The paper does not sufficiently acknowledge or discuss these existing approaches, missing an opportunity to position its contribution within the broader context of optimization techniques.\n\n2. While each subproblem is solved under a $k$-sparsity constraint, it is not explicitly clear how the framework ensures that the combined solution across all subproblems maintains the overall $k$-sparsity required by the Sparse PCA problem. This raises concerns about the validity of the final solution's sparsity, which is crucial for interpretability and efficiency.\n\n[1] \"Exact covariance thresholding into connected components for large-scale graphical lasso.\" Journal of Machine Learning Research, 13(27):781−794, 2012.\n\n[2] \"Graphical lasso and thresholding: Equivalence and closed-form solutions.\" Journal of Machine Learning Research, 20(10):1−44, 2019.\n\n[3] \"Learning large-scale MTP2 Gaussian graphical models via bridge-block decomposition.\" Advances in Neural Information Processing Systems 36:73211-73231, 2023.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "c9F3uxROy4", "reviewer": "Reviewer_ZrrB", "summary": "This paper solves approximate Sparse PCA with matrix block-diagonalization. The proposed algorithm has 3 steps:\napproximate a given matrix by a block diagonal, \nsolve separate problems and then, combine them. \n \nThe matrix preparation step re-arranges the rows and columns so that most of the energy of the entries lies in a block diagonal structure and zeros out the rest of the entries. This part is what I think is the hard part to compute and the hardness of this approximation ( possibly harder than Sparse PCA to start with) is one confusing thing. \n\nThe central theoretical innovation of the paper is the following structural property of an SPCA: The support of the optimal solution must always be contained within one block.  I think this is a cool observation (but not some deep theorem) and I have a few questions about it that I discuss later. \n\nThe combination of solutions, after this observation is realized, is clear: you only pick the solution from the correct sub-block, since there are no interaction terms to worry about across blocks. This is cool, and all subsequent running time results and combinations with other algorithms are clear.\n\nEdit after reading replies and discussion: I think the authors sort of addressed my concern, I give a small upgrade in my score", "questions": "q1) I was very confused about the matrix preparation part of the algorithm (i.e. best approximation by block -diagonal). I would think that finding the best block-diagonal matrix that approximates a given matrix A (using any reasonable metric of distance) is NP-hard. Are the authors claiming to achieve this in polynomial time for any matrix using their L_inf distance? \nI was a little lost on the notation and didn't understand where the search over all permutations of rows and columns is addressed. \nPlease clarify the computational complexity of this problem and how the algorithm gets around checking all row and column permutations. \n\n\nq2) I do not see why theorem 2 is so challenging? I think its a cool observation but I think of it simply as: There is no reason for the principal component to put any energy on positions where the off-diagonal entries are zero, since it gets no benefit from that. Therefore it will pick the block with the most energy and put all its l2 mass there. For example if a matrix is just diagonal, it will put all its energy on the largest diagonal entry in absolute value. \n\nq3) Theorem 2 doesn't really use sparsity anywhere right? Its really a result about PCA for block-diagonal matrices, or do i misunderstand something?\n\nRemark: Approximating matrices with block-diagonal is indeed a very interesting topic.See this interesting question and answer by Gowers on discovering near-block diagonal structure of a matrix by re-ordering and its connections to combinatorics: \nhttps://mathoverflow.net/questions/68041/showing-block-diagonal-structure-of-matrix-by-reordering\nAnother flavor is the \"Bandwidth minimization problem\". There instead of blocks one tries to re-arrange columns and rows so that all the non-zeros are on a few bands around the diagonal. This is known to be np-hard even for 3 bands, see: \nComplexity Results for Bandwidth Minimization\nM. R. Garey, R. L. Graham, D. S. Johnson, D. E. Knuth\n\nq4) On experimental evaluation. \nWhy don't the authors compare to other algorithms as done in Papailiopoulos et al. ICML 2013 and the algorithms compared from that paper? There is (FullPath) of (d’Aspremont et al., 2007b) and  \nthe truncated power method (TPower) of (Yuan & Zhang, 2011). I believe Tpower is usually the fastest in practice. \nI think the Chan paper may have the best theoretical performance but it's not clear how it performs in practice. I did not see any performance plots in that paper. \n\nI'm willing to improve my score if the authors manage to convince me on these concerns.", "limitations": "", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper has a cool observation and I think is well written. \nI have some concerns about the complexity of the approximation of a matrix by a block diagonal that I would really like clarified. \nThe experimental section is reasonable some it misses additional experiments that I would like to see.", "weaknesses": "I have issues with the complexity of the block-diagonal matrix approximation and also with the experimental section, as I clarify below. \n\nAlso some minor typos here:\nDefinition 1, I assume Aij and \\tilde{A}_ij are the entries of the matrix, would be good to mention. I know its mentioned later in additional notation but would be better to include in this definition. \n\ntypo: On the other hand, there also exist a number of algorithms that takes* polynomial runtime\n\ntypo: The axies* are indices of A; \n\ntypo: th3, positive integer, and donote*", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "8BBbsQTNZj", "reviewer": "Reviewer_vGQN", "summary": "The paper introduces a new framework for efficiently approximating Sparse Principal Component Analysis (Sparse PCA) by transforming the input covariance matrix into a block-diagonal structure, allowing for computationally manageable sub-problems.  The proposed method involves three key steps: creating a block-diagonal approximation of the matrix, solving Sparse PCA for each block, and then reconstructing an approximate solution for the entire matrix. By focusing on smaller blocks, this approach achieves substantial speedups with minimal loss in accuracy. The framework is adaptable, integrating well with existing Sparse PCA algorithms and reducing overall complexity theoretical and emperically.", "questions": "See above.", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The idea is well-motivated, and the problem is relevant to the community. Despite the NP-hardness of sparse PCA (SPCA), the authors propose addressing it through matrix block diagonalization. This framework demonstrates advantages in time complexity over traditional methods, both theoretically and empirically. Additionally, the authors discuss how to determine the appropriate SNR threshold, $\\epsilon$, within a statistical model $A = \\widetilde{A} + E $ using the proposed algorithm.", "weaknesses": "The authors investigate the recovery of the first individual eigenvector and ensure correctness by establishing an upper bound on the gap between the corresponding eigenvalues in Theorem 1. However, the situation changes when considering the principal subspaces of the covariance matrix $\\Sigma$, which are spanned by sparse leading eigenvectors. When leading eigenvalues are identical or close to each other, individual eigenvectors may become unidentifiable. Could the analysis in Theorem 1 be extended to handle the aformentioned case?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "kEHVxcubuX", "reviewer": "Reviewer_uKtp", "summary": "This paper presents a block-diagonal-based approach to speed up Sparse PCA, preserving an average approximation error. The authors provide mathematical proof for this algorithm, and extensive testing on real-world datasets demonstrates consistent improvements that align with theoretical predictions.", "questions": "1. In your algorithm, you noted that the input matrix could be solved using $p$ sub-matrices. Could you please give more details on determining $p$ or explaining the relationship of $d_i$ (besides stating that $\\sum d_i = d$)?\n2. Could you explain why $|| E ||_\\infty$ is the estimation in Model 1?\n3. In algorithm 4, why $m=\\lfloor (2C+2)d^{1+\\alpha}  \\rfloor$?\n4. I'm curious about the inequality in Line 71. Can you prove it?\n5. The proof of Theorem 1 is incomplete. Could you complete it?\n6. In the proof of Theorem 3, how do you come up with $k || A-A^\\epsilon ||_\\infty \\leq \\epsilon $?\n7. Could you prove why the g() function in Eq.(2) is convex?\n8. Do you think the inequality in Line 41-45 is too relaxed?\n9. Some descriptions are not rigorous. For example, \"an absolute constant $c > 0$\". Apparently, an absolute value is greater than 0 in this case.\n10. Including additional details in the proof of Theorem 4 would improve clarity. For instance, explaining how the running time is derived and whether the relaxation is reasonable would be helpful.\n11. I noticed that you frequently use expressions like $ \\max \\max $ or $ \\inf \\inf \\inf $ in your proof, which appears unprofessional. Could you consider rephrasing these using constraints for clarity?\n12. Some sentences are too long for readers to understand. For example, Line 51-53. Could you rephrase them?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1. This paper is the first to apply a block-diagonal-based method to the Sparse PCA problem, enhancing the balance between runtime efficiency and average approximation error.  \n2. The authors present theoretical analysis, providing guarantees for both approximation error and runtime complexity.  \n3. This approach is evaluated across many real-world datasets, demonstrating consistent improvements in runtime and approximation error.", "weaknesses": "1. In your algorithm, you noted that the input matrix could be solved using $p$ sub-matrices, but you did not provide a method to determine $p$ or explain the relationship of $d_i$ (besides stating that $\\sum d_i = d$).\n2. There is not much discussion about the impact of hyperparameters.\n3. In your notation, $p$ in $p$-norm should not represent the same value as $p$ in diag($A_1, \\cdots, A_p$).\n4. Some descriptions in your Theorem, Proposition, and Proof are not rigorous.\n5. Some proofs are not completely correct.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "bzWQsD4h4b", "reviewer": "Reviewer_fTk5", "summary": "This paper proposes to address the one-component sparse PCA (SPCA) problem by first approximating the target covariance matrix by a block-diagonal matrix (by thresholding its entries and permuting rows and columns), and then looking for solutions within each one of the blocks. These smaller SPCA problems associated with each block can be (exactly or approximately) solved by means of existing SPCA algorithms, and are expected to require (potentially much) less compute in comparison with the standard approach. The authors state and prove some results bounding the suboptimality of this scheme (in terms of objective function values) in terms of the infinity norm of the matrix approximation error, the sparsity level and the largest size among the obtained blocks, including in their analysis the case in which an approximate algorithm with additive and multiplicative errors is used. Simulation results are given, comparing the proposed scheme with two other approaches: exact solution by branch-and-bound, and approximate solution by Chan's algorithm (reference Chan et al., 2015 cited by the authors).", "questions": "I suggest that the authors take the above stated weaknesses into account to improve their manuscript. Apart from this suggestion:\n\nQ1) How far from the bound implied by Theorem 1 are the measured suboptimality gaps in practice?\n\nQ2) How realistic and how restrictive is the bound on $d^{1-\\alpha}$ of Proposition 1? How can it be interpreted? A discussion should be presented on this point, as it is not obvious.\n\nQ3) How hard is it to extend your results to the multi-component case?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "S1) The proposed scheme is simple and intuitively well-motivated.\n\nS2) The theoretical results are also fairly simple and support the intuition.\n\nS3) The analysis includes approximate algorithms, which have great practical relevance since the problem is computationally hard.\n\nS4) The simulation results demonstrate a significant speedup in large problem instances, most often with only a small relative suboptimality gap. Hence, the proposed algorithm is clearly useful for practitioners.", "weaknesses": "W1) Some of the technical contents, including the problem formulation, are not properly formalized, with some important aspects completely omitted. Notably, it should be stated that matrix $A$ in (SPCA) is necessarily symmetric positive semidefinite, implying the same for the blocks in the block-diagonal approximation. The absence of this information and the fact that throughout the paper $A$ is simply referred to as an \"input matrix\" rather than a covariance matrix may mislead the reader into thinking that the problem is more general than it actually is.\n\nW2) The presentation of the simulation results is somewhat superficial, focusing only on presenting and briefly commenting the two quantitative criteria used for comparison, without much discussion or insight into what is going on. Specifically:\n- Separate values for the different used $k$ should be reported (see point W2 below).\n- It would be interesting to report the threshold value $\\varepsilon$ effectively chosen by the algorithm (relatively to the infinity norm of the input matrix), as well as the proportion of zero entries after the thresholding. \n- It would also be interesting to compare the support of the solution obtained by the proposed scheme with that obtained by the baseline methods (e.g., using a Jaccard index).\n\nW3) Reporting average results with respect to $k$ is not a reasonable choice in my opinion, as the statistics of the chosen metrics probably behave very differently for different values of $k$. \n\nW4) As it is, I don't see the utility of Section 4.1. First, this model is not applied to any of the datasets used in the experiments. This leads one to suspect that in practice it is quite hard to come up with estimates of the parameters required by Algorithm 4. Second, the authors do not even generate synthetic data following such a model (which is one typical use of a model) in order to illustrate the obtained results. In my view results of this kind should be added, or else the contents of 4.1 should be moved to the appendices as they're not really important (or both).\n\nW5) Though generally well-written, the paper lacks some clarity at times, and the notation is not always consistent/clear. In particular:\n- The sentence \"This result is non-trivial: while the support of an optimal solution could span multiple blocks, we theoretically show that there must exist a block that contains the support of an optimal solution, which guarantees the efficiency of our framework.\" seems to be contradicatory. I believe that the authors mean the following: *one could think* that the solution could span multiple blocks, but they show this is not true. The same goes for Remark 1.\n- What is the difference between $A^\\varepsilon$ and $\\tilde{A}$ in Theorem 1? It seems that two different symbols are used to denote the same object.\n- The constraint $\\|x\\|_0 \\le k$ in the last problem that appears in the proof of Theorem 2 is inocuous, since the size of each block $\\tilde{A}_i'$ is at most $k$ anyway. This should be commented.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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:

Decision: Accept (Poster)