Efficient Sparse PCA via Block-Diagonalization

## Abstract

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 ) 𝑔 𝑘 𝑑 g(k,d) italic_g ( italic_k , italic_d ) is the runtime of an algorithm (approximately) solving Sparse PCA in dimension d 𝑑 d italic_d and with sparsity value k 𝑘 k italic_k .
Our framework, when integrated with this algorithm, reduces the runtime to 𝒪 ⁢ ( d d ⋆ ⋅ g ⁢ ( k , d ⋆ ) + d 2 ) 𝒪 ⋅ 𝑑 superscript 𝑑 ⋆ 𝑔 𝑘 superscript 𝑑 ⋆ superscript 𝑑 2 \mathcal{O}\left(\frac{d}{d^{\star}}\cdot g(k,d^{\star})+d^{2}\right) caligraphic_O ( divide start_ARG italic_d end_ARG start_ARG italic_d start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT end_ARG ⋅ italic_g ( italic_k , italic_d start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT ) + italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) , where d ⋆ ≤ d superscript 𝑑 ⋆ 𝑑 d^{\star}\leq d italic_d start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT ≤ italic_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 ) = 𝒪 ⁢ ( k 3 ⋅ d k ) 𝑔 𝑘 𝑑 𝒪 ⋅ superscript 𝑘 3 superscript 𝑑 𝑘 g(k,d)=\mathcal{O}(k^{3}\cdot d^{k}) italic_g ( italic_k , italic_d ) = caligraphic_O ( italic_k start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ⋅ italic_d start_POSTSUPERSCRIPT italic_k end_POSTSUPERSCRIPT ) to 𝒪 ⁢ ( k 3 ⋅ d ⋅ ( d ⋆ ) k − 1 ) 𝒪 ⋅ superscript 𝑘 3 𝑑 superscript superscript 𝑑 ⋆ 𝑘 1 \mathcal{O}(k^{3}\cdot d\cdot(d^{\star})^{k-1}) caligraphic_O ( italic_k start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ⋅ italic_d ⋅ ( italic_d start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT italic_k - 1 end_POSTSUPERSCRIPT ) , demonstrating exponential speedups if d ⋆ superscript 𝑑 ⋆ d^{\star} italic_d start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT 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 93.77, while maintaining an average approximation error of 2.15%; for approximate Sparse PCA algorithm, our method achieves an average speedup factor of 6.77 and an average approximation error of merely 0.37%.

## 1 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 $k\ll d$ many features, thus integrating a sparsity constraint.
This constraint significantly enhances interpretability, which is essential in data analysis when dealing with a large number of features.
Sparse PCA has found widespread application across various domains, including text data analysis (Zhang and Ghaoui, ), cancer research (Hsu et al., ), bioinformatics (Ma and Dai, ), and neuroscience (Zhuang et al., ).
For further reading, we refer interested readers to a comprehensive survey on Sparse PCA (Zou and Xue, ).

While being important and useful, Sparse PCA is a challenging problem—it is known to be NP-hard (Magdon-Ismail, ).
Solving Sparse PCA exactly, such as through well-known Branch-and-Bound algorithm (Berk and Bertsimas, ), requires worst-case exponential runtime.
Moreover, while there are near-optimal approximation algorithms that exhibit faster performance empirically, they still require exponential runtime in the worst-case scenarios (Bertsimas et al., ; Li and Xie, ; Cory-Wright and Pauphilet, ; Dey et al., ; Zou et al., ; Dey et al., ).
This complexity impedes their applicability to large-scale datasets.
On the other hand, numerous polynomial-time approximation algorithms have been proposed for Sparse PCA (Chowdhury et al., ; Li and Xie, ; Papailiopoulos et al., ; Chan et al., ; Del Pia, ; Asteris et al., ).
These algorithms typically face trade-offs between efficiency and solution quality.
They either achieve sub-optimal solutions with high efficiency, deliver high-quality solutions at the cost of greater computational complexity, or produce high-quality solutions under specific assumptions.

In this paper, we introduce a novel framework designed to efficiently approximate Sparse PCA through matrix block-diagonalization.
Our framework
facilitates the use of off-the-shelf algorithms in a plug-and-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 process includes denoising by zeroing out entries below a specific threshold and then grouping indices that correspond to non-zero entries into blocks;
(ii) *Sub-Problem Solution:* Solve Sparse PCA sub-problems using any known algorithm in each block with a smaller dimension, and find out the best solution among these sub-problems;
(iii) *Solution Reconstruction:* Reconstruct the best solution to the original Sparse PCA problem.
We illustrate these steps in Fig. 1.
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 (iii)).

Figure: Figure 1: Illustration of our proposed Sparse PCA approach, given a $9\times 9$ covariance input matrix. (i) Entries away from zero are highlighted in the upper matrix (original input matrix $A$). Then zero out other entries, sort the matrix, and obtain the lower block-diagonal approximation matrix. Heatmaps are used to present the values of matrix entries. The axies are indices of $A$; (ii) Extract sub-matrices from the block-diagonal approximation, and solve sub-problems via a suitable Sparse PCA algorithm; (iii) Select the solution with the highest objective value obtained from the sub-problems. Construct a 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.
Refer to caption: extracted/5903005/illustration1.png

Our framework significantly speeds up the computation of Sparse PCA, with negligible approximation error.
We theoretically quantify the computational speedup and approximation error of our framework in Sections 3 and 4, and conduct large scale empirical evaluations in Section 5.
Next, we illustrate the performance of our method via a concrete example.
Suppose one is given an input covariance matrix being a noisy sorted block-diagonal matrix with 10 blocks, each of size 100, and is asked to solve Sparse PCA with a sparsity constant $4$.
The Branch-and-Bound algorithm (Berk and Bertsimas, ) takes over 3600 seconds and obtains a solution with objective value 13.649.
Our framework integrated with Branch-and-Bound algorithm reduces runtime to 17 seconds, a speedup by a factor of 211, and obtains a solution with objective value 13.637, with an approximation error of 0.09%.

Our contributions.
We summarize our contributions in the following:

- •
We propose a 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 a very easy way, i.e., 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-diagonalization to solving Sparse PCA problem.
- •
We provide approximation guarantees for our framework, showing that we can obtain a solution with an additive error that is linear in the approximation error of the block-diagonal approximate matrix.
We show that when integrated with an approximation algorithm, we can generally obtain a better multiplicative factor, with an additional cost of an additive error.
Under certain assumptions, the additive error can also be improved.
- •
We also perform a time complexity analysis showing that for an (approximation) algorithm with runtime $g(k,d)$, our framework could reduce the runtime to $\mathcal{O}(\frac{d}{d^{\star}}\cdot g(k,d^{\star})+d^{2})$, where $d$ is the dimension of Sparse PCA input covariance matrix, $k$ is the sparsity constant, and $d^{\star}$ is the largest block size in the block-diagonal approximate matrix.
- •
We show that for a 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 to 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 93.77, while maintaining an average approximation error of 2.15%, and (ii) when integrated with Chan’s algorithm (Chan et al., ), our framework achieves an average speedup factor of 6.77 and an average approximation error of 0.37%.

Our techniques.
In this paper, we also introduce several novel techniques that may be of independent interest.
We characterize a structural property of an optimal solution to Sparse PCA, when given a block-diagonal input matrix.
We establish that solving Sparse PCA for a block-diagonal matrix input is 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 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.
Additionally, leveraging the defined structure, we propose a novel binary-search based method to identify the best block-diagonal approximation, particularly when computational resources are 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 3, we provide details of operational procedures of our proposed framework, and develop guarantees of approximation errors as well as runtime complexity.
In Section 4 we further extend our algorithms to learning problems either under (i) a statistical model, or under (ii) a model-free setting, and provide theoretical guarantees respectively.
In Section 5, we 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-diagonalization methods in the literature, detailed proofs across the sections, 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:

$$ $\displaystyle\mathsf{OPT}\vcentcolon=\max\ x^{\top}Ax\quad\text{s.t. }\left \lVert x\right\rVert_{2}=1,\ \ \left\lVert x\right\rVert_{0}\leq k.$ (SPCA) $$

In the problem, one is given $A\in{\mathbb{R}}^{d\times d}$ as an input data matrix, and a positive integer $k\leq d$ denoting the sparsity constant, i.e., the upper bound on the number of non-zero entries in the solution vector $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 a user is given an input block-diagonal covariance matrix $A$ of size $d$, comprising $d/d^{\star}$ many $d^{\star}\times d^{\star}$ blocks (assume that $d$ is divisible by $d^{\star}$), and is asked to solve SPCA with sparsity constant $k$.
The user can simply apply an exact algorithm (e.g., Branch-and-Bound algorithm) to solve the large problem in time $\mathcal{O}(k^{3}\cdot\binom{d}{k})=\mathcal{O}(k^{3}\cdot d^{k})$.
The user could also solve $d/d^{\star}$ many sub-problems in each block, potentially obtaining an optimal solution (a fact we rigidly prove in Section 3.2.1), in time $\mathcal{O}((d/d^{\star})\cdot k^{3}\cdot(d^{\star})^{k})$.
This segmented approach yields a significant speedup factor of $\mathcal{O}((d/d^{\star})^{k-1})$.

However, not all input matrices $A$ are of the block-diagonal structure, even after permuting its rows and columns.
One goal of our work is to develop a reasonable way to construct a (sorted) block-diagonal matrix approximation $\widetilde{A}$ to $A$ such that they are close enough.
In this paper, we are interested in finding the following *$\varepsilon$-matrix approximation to $A$*:

###### Definition 1 ( $\varepsilon$ -matrix approximation) .

Given a matrix $A\in{\mathbb{R}}^{d\times d}$, an $\varepsilon$-matrix approximation of $A$ is defined as a $d\times d$ matrix $\widetilde{A}$ such that

$$ $\displaystyle\max_{i,j\in[d]}\left\lvert A_{ij}-\widetilde{A}_{ij}\right\rvert \leq\varepsilon.$ $$

Moreover, we define the set of all $\varepsilon$-matrix approximations of $A$ to be $\mathcal{B}(A,\varepsilon)$.

In other words, $\widetilde{A}\in\mathcal{B}(A,\varepsilon)$ if and only if $\widetilde{A}$ belongs to an $\ell_{\infty}$-ball centered at $A$ with radius $\varepsilon$.
Our overarching strategy is that if the input matrix is not block-diagonal, we find a matrix approximation in the ball so that it could be resorted to be 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 non-asymptotic big-oh notation: For functions
$f,g:\mathcal{Z}\to\mathbb{R}_{+}$, we write $f=\mathcal{O}(g)$ (resp. $f=\Omega(g)$) if there exists a constant
$C>0$ such that $f(z)\leq{}Cg(z)$ (resp. $f(z)\geq{}Cg(z)$)
for all $z\in\mathcal{Z}$.
For an integer $n\in\mathbb{N}$, we let $[n]$ denote the set $\{1,2,\dots,n\}$.
For a vector $z\in\mathbb{R}^{d}$, we use $z_{i}$ to denote its $i$-th entry.
We denote $\operatorname{supp}(z)\vcentcolon=\{i\in[d]:z_{i}\neq 0\}$ the *support* of $z$.
For $1\leq p\leq\infty$, we denote $\|x\|_{p}$ the $p$-norm of $z$, i.e., $\|x\|_{p}\vcentcolon=(\sum_{i=1}^{d}\left\lvert x_{i}\right\rvert^{p})^{1/p}$, and $\left\lVert x\right\rVert_{\infty}\vcentcolon=\max_{i\in[d]}\left\lvert x_{i}\right\rvert$.
For a matrix $A\in{\mathbb{R}}^{n\times m}$, we use $A_{ij}$ to denote its $(i,j)$-th entry.
We define $\left\lVert A\right\rVert_{\infty}\vcentcolon=\max_{i,j}|A_{ij}|$.
For index sets $S\subseteq[n]$, $T\subseteq[m]$, we denote $A_{S,T}$ the $|S|\times|T|$ sub-matrix of $A$ with row index $S$ and column index $T$.
For matrices $A_{i}\in{\mathbb{R}}^{d_{i}\times d_{i}}$ for $i=1,2,\ldots,p$, we denote $\mathsf{diag}(A_{1},A_{2},\ldots,A_{p})$ the $(\sum_{i=1}^{p}d_{i})\times(\sum_{i=1}^{p}d_{i})$ block-diagonal matrix with blocks $A_{1},A_{2},\ldots,A_{p}$.
For a square matrix $B\in{\mathbb{R}}^{d\times d}$, we define the *size* of $B$ to be $d$.

###### Definition 1 ( ε 𝜀 \varepsilon italic_ε -matrix approximation) .

## 3 Our approach

In this section, we introduce our framework for solving SPCA approximately via block-diagonalization technique.
We explain in detail the operational procedures of our framework in Section 3.1.
In Section 3.2, we characterize the approximation guarantees and computational speedups of our framework.
Note that, in this section, our proposed procedures require a predefined threshold $\varepsilon>0$.
In Section 4, we extend our algorithms to learning the threshold in a statistical model, as well as in a model-free setting.

### 3.1 Algorithms

In this section, we outline the operational procedure of our proposed framework for solving SPCA approximately.
The procedure commences with the input matrix $A\in\mathbb{R}^{d\times d}$ and a predefined threshold value $\varepsilon>0$.
Initially, a denoised matrix $A^{\varepsilon}$ is derived through thresholding as described in Algorithm 1.
Subsequently, several sub-matrices and their corresponding sets of indices are extracted from $A^{\varepsilon}$ by employing Algorithm 2, which groups index pair $(i,j)$ if $A^{\varepsilon}_{ij}\neq 0$.
This grouping mechanism can be efficiently implemented using a Depth-First Search approach, which iterates over each unvisited index $i$ and examines all non-zero $A_{ij}^{\varepsilon}$, and then visit each corresponding index $j$ that is unvisited.
The computational complexity of Algorithm 2 is bounded by $\mathcal{O}(d^{2})$.
Finally, an approximate solution to SPCA is obtained by solving several Sparse PCA sub-problems in each block using a certain (approximation) algorithm $\mathcal{A}$, mapping the solution back to the original problem, and finding out the best solution.
The operational details are further elaborated in Algorithm 3(^1^11We note that on line 6 in Algorithm 3, if $k$ is larger than the size of $\widetilde{A}_{i}$, a PCA problem is solved instead.).

Figure: Algorithm 1 Denoising procedure via thresholding

Figure: Algorithm 2 Matrix block-diagonlization

Figure: Algorithm 3 Efficient Sparse PCA via block-diagonalization

### 3.2 Approximation guarantees and computational speedups

In this section, we provide analysis of approximation guarantees and speedups of our framework Algorithm 3.
Before showing approximation guarantees, we first define the *$\varepsilon$-intrinsic dimension of $A$*:

###### Definition 2 ( $\varepsilon$ -intrinsic dimension) .

Let $A\in{\mathbb{R}}^{d\times d}$.
The $\varepsilon$-intrinsic dimension of $A$, denoted by $\operatorname*{int\,dim}(A,\varepsilon)$, is defined as:

$$ $\displaystyle\operatorname*{int\,dim}(A,\varepsilon)\vcentcolon=\min_{B\in \mathcal{B}(A,\varepsilon)}\textup{lbs}(B),$ $$

where $\textup{lbs}(B)$ denotes the largest size of sub-matrices outputted by Algorithm 2 given $B$ as input.

As we will see later in this section, this concept plays a crucial role in characterizing both the quality of the approximate solution obtained through our framework Algorithm 3, and the computational speedup achieved by our framework.
Moreover, a matrix $\widetilde{A}$ can be efficiently identified such that $\textup{lbs}(\widetilde{A})=\operatorname*{int\,dim}(A,\varepsilon)$ using Algorithm 1, as shown in Lemma 1.

###### Lemma 1 .

Let $A\in{\mathbb{R}}^{d\times d}$ and $\varepsilon>0$.
Given input $(A,\varepsilon)$, Algorithm 1 outputs an $\varepsilon$-approximation of $A$, denoted as $\widetilde{A}$, such that $\textup{lbs}(\widetilde{A})=\operatorname*{int\,dim}(A,\varepsilon)$ in time $\mathcal{O}(d^{2})$.

###### Definition 2 ( ε 𝜀 \varepsilon italic_ε -intrinsic dimension) .

###### Lemma 1 .

#### 3.2.1 Approximation guarantees

In this section, we provide approximation guarantees for our framework Algorithm 3.
We first show that in Theorem 1 that, an optimal solution to SPCA with input $(A,k)$ and an optimal solution to $(A^{\varepsilon},k)$ are close to each other in terms of objective value:

###### Theorem 1 .

Let $A\in{\mathbb{R}}^{d\times d}$, $k\leq d$, and $\varepsilon>0$.
Denote by $A^{\varepsilon}$ the output of Algorithm 1 with input $(A,\varepsilon)$, $\widetilde{x}$ the optimal solution to SPCA with input $(\widetilde{A},k)$, and $x^{\star}$ the optimal solution to SPCA with input $(A,k)$.
Then, it follows that

$$ $\displaystyle\left\lvert(x^{\star})^{\top}Ax^{\star}-\widetilde{x}^{\top}A \widetilde{x}\right\rvert\leq 2k\cdot\varepsilon$ $$

Then, in Theorem 2, we show that solving SPCA exactly with input $(\widetilde{A},k)$ is equivalent to addressing SPCA sub-problems in blocks within $\widetilde{A}$:

###### Theorem 2 .

Let $\widetilde{A}=\mathsf{diag}(\widetilde{A}_{1},\widetilde{A}_{2},\ldots,
\widetilde{A}_{p})$.
Denote $\mathsf{OPT}$ to be the optimal value to SPCA with input pair $(\widetilde{A},k)$.
Let $\mathsf{OPT}_{i}$ to be the optimal value to SPCA with input pair $(\widetilde{A}_{i},k)$, for $i\in[p]$.
Then, one has $\mathsf{OPT}=\max_{i\in[p]}\mathsf{OPT}_{i}.$

###### Remark 1 .

Theorem 2 is particularly non-trivial: without this result, the support of any optimal solution to SPCA may span different blocks.
Solving SPCA with a block-diagonal input would consequently involve iterating over all possible sparsity constants (that sum up to $k$) for each block-diagonal sub-matrix, making our framework proposed in Section 3.1 nearly impractical.

With these two theorems established, we now extend Theorem 1 to incorporate the use of various approximation algorithms.
We note that the established approximation guarantee follows from subroutines as described in Theorem 1.
Furthermore, Theorem 3 suggests that our framework Algorithm 3, can also be highly effective when combined with other approximation algorithms.
Specifically, it has the potential to find a better solution, particularly when the value of $\varepsilon$ is small.

###### Theorem 3 .

Let $A\in{\mathbb{R}}^{d\times d}$ and let $k$ be a positive integer, and donote by $x^{\star}\in{\mathbb{R}}^{d}$ an optimal solution to SPCA with input $(A,k)$.
Suppose that an algorithm $\mathcal{A}$ for SPCA with input $(A,k)$ finds an approximate solution $x\in{\mathbb{R}}^{d}$ to SPCA such that has multiplicative factor $m(k,d)\geq 1$ and additive error $a(k,d)\geq 0$, i.e., one has $x^{\top}Ax\geq(x^{\star})^{\top}Ax^{\star}/m(k,d)-a(k,d)$.
Furthermore, we assume that $m(k,d)$ and $a(k,d)$ is non-decreasing with respect to $d$.
For $\varepsilon>0$, denote by $y\in{\mathbb{R}}^{d}$ be the output of Algorithm 3 with input tuple $(A,k,\varepsilon,\mathcal{A})$.
Then, one has

$$ $\displaystyle y^{\top}Ay\geq\frac{(x^{\star})^{\top}Ax^{\star}}{m\left(k, \operatorname*{int\,dim}(A,\varepsilon)\right)}-a\left(k,\operatorname*{int\, dim}(A,\varepsilon)\right)-\left(1+\frac{1}{m\left(k,\operatorname*{int\,dim}( A,\varepsilon)\right)}\right)\cdot k\varepsilon.$ $$

###### Remark 2 .

Theorem 3 implies that our framework, Algorithm 3, when inputted with $(A,k,\varepsilon,\mathcal{A})$, could find an approximate solution with a better multiplicative factor of $m\left(k,\operatorname*{int\,dim}(A,\varepsilon)\right)$, and an additive error of $a\left(k,\operatorname*{int\,dim}(A,\varepsilon)\right)+k\varepsilon$.
Note that the additive error is also improved if

$$ $\displaystyle a(k,d)\geq a\left(k,\operatorname*{int\,dim}(A,\varepsilon) \right)+\left(1+\frac{1}{m\left(k,\operatorname*{int\,dim}(A,\varepsilon) \right)}\right)k\varepsilon.$ $$

In addition, we have found instances where our framework can indeed obtain better solutions when integrated with some approximation algorithm, as shown in Table 10 in Appendix C.
Finally, we summarized in Table 1 the change of multiplicative factor and additive error when integrated our framework to some existing (approximation) algorithms for SPCA.

**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 and Xie (). Note that our MF is always no larger than the original MF.**
| Algorithm | MF | AE | Our MF | Our AE |
| --- | --- | --- | --- | --- |
| Chan et al. () | $\min(\sqrt{k},d^{1/3})$ | 0 | $\min(\sqrt{k},\operatorname*{int\,dim}(A,\varepsilon)^{1/3})$ | $(1+1/\text{Our MF})\cdot k\varepsilon$ |
| MISDP | $\min(k,d/k)$ | 0 | $\min(k,\operatorname*{int\,dim}(A,\varepsilon)/k)$ | $(1+1/\text{Our MF})\cdot k\varepsilon$ |
| Greedy | $k$ | 0 | $k$ | $(1+1/k)\cdot k\varepsilon$ |
| Exact algorithms | 1 | 0 | 1 | $2k\varepsilon$ |
| Asteris et al. () | 1 | $\delta$ | 1 | $2k\varepsilon+\delta$ |

###### Theorem 1 .

###### Theorem 2 .

###### Remark 1 .

###### Theorem 3 .

###### Remark 2 .

#### 3.2.2 Computational speedups

In this section, we characterize the runtime complexity of Algorithm 3.
We will see that the $\varepsilon$-intrinsic dimension also plays an important role in this section.

###### Proposition 1 .

Let $A\in{\mathbb{R}}^{d\times d}$, and let $k$ be a positive integer.
Suppose that for a specific algorithm $\mathcal{A}$, the time complexity of using $\mathcal{A}$ for (approximately) solving SPCA with input $(A,k)$ is a function $g(k,d)$ which is convex and non-decreasing with respect to $d$.
We assume that $g(k,1)=1$.
Then, for a given threshold $\varepsilon>0$, the runtime of Algorithm 3 with input tuples $(A,k,\varepsilon,\mathcal{A})$ is at most

$$ $\displaystyle\mathcal{O}\left(\left\lceil\frac{d}{\operatorname*{int\,dim}(A, \varepsilon)}\right\rceil\cdot g\left(k,\operatorname*{int\,dim}(A,\varepsilon )\right)+d^{2}\right)$ $$

###### Remark 3 .

We remark that the assumption $g(k,d)$ is convex and non-decreasing with respect to $d$ in Proposition 1 is a very weak assumption.
Common examples include $d^{\alpha}$ and $d^{\alpha}\log{d}$ with $\alpha\geq 1$.

###### Remark 4 .

In this remark, we explore the computational speedups provided by Algorithm 3 by modeling $\operatorname*{int\,dim}(A,\varepsilon)$ as a function of $A$ and $\varepsilon$.
We assume without loss of generality that $\left\lVert A\right\rVert_{\infty}=1$, thereby setting $\operatorname*{int\,dim}(A,1)=0$.
It is important to note that while $\operatorname*{int\,dim}(A,\varepsilon)$ is discontinuous, for the purpose of illustration in Table 2, it is treated as a continuous function of $\varepsilon$.

**Table 2: Comparison of original runtimes and current runtimes in our framework. We assume that the algorithm has runtime $\mathcal{O}(d^{\alpha})$, for some $\alpha\geq 2$.**
| $\operatorname*{int\,dim}(A,\varepsilon)$ | Runtime | Our Runtime | Speedup factor |
| --- | --- | --- | --- |
| $d(1-e^{-c(1-\varepsilon)})$ | $\mathcal{O}(d^{\alpha})$ | $\mathcal{O}(d^{\alpha}(1-e^{-c(1-\varepsilon)})^{\alpha-1})$ | $\mathcal{O}((1-e^{-c(1-\varepsilon)})^{1-\alpha})$ |
| $d(1-\varepsilon)$ | $\mathcal{O}(d^{\alpha})$ | $\mathcal{O}(d^{\alpha}(1-\varepsilon)^{\alpha-1})$ | $\mathcal{O}((1-\varepsilon)^{1-\alpha})$ |
| $d(1-\varepsilon)^{m}$ | $\mathcal{O}(d^{\alpha})$ | $\mathcal{O}(d^{\alpha}(1-\varepsilon)^{(\alpha-1)m})$ | $\mathcal{O}((1-\varepsilon)^{(1-\alpha)m})$ |

###### Proposition 1 .

###### Remark 3 .

###### Remark 4 .

## 4 Extensions

Our algorithm presented in Section 3.1 takes a known threshold as an input. In this section, we remove this constraint by extending our algorithm to various settings.
Specifically, in Section 4.1, we consider learning under a statistical model, where one can do the denoising of a matrix easily.
In Section 4.2, we propose a practical learning approach without any assumption.

### 4.1 Learning under a statistical model

A crucial aspect of our proposed framework involves determining an appropriate threshold, $\varepsilon$, as required by Algorithm 3.
While specifying such a threshold is generally necessary, there are many instances where it can be effectively estimated from the data.
In this section, we discuss a model commonly employed in robust statistics (Comminges et al., ; Kotekal and Gao, , ).
This model finds broad applications across various fields, including clustering (Chen and Witten, ) and sparse linear regression (Minsker et al., ).

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 a few other features, it follows that $A$ should be, or at least be close to, a block-diagonal matrix.
This observation motivates us for considering 1:

###### Model 1 .

The input matrix $A\in{\mathbb{R}}^{d\times d}$ is close to a block-diagonal symmetric matrix $\widetilde{A}\in{\mathbb{R}}^{d\times d}$, i.e., $A=\widetilde{A}+E$ for some symmetric matrix $E\in{\mathbb{R}}^{d\times d}$ such that

- (i)
for $i\leq j$, each entry $E_{ij}$ is drawn from an i.i.d. $\varrho$-sub-Gaussian distribution, with mean zero, and variance $\sigma^{2}$, where $\varrho>0$ and $\sigma^{2}$ being unknown, but an upper bound for $\varrho^{2}/\sigma^{2}$ is given;
- (ii)
write $\widetilde{A}=\mathsf{diag}\left(\widetilde{A}_{1},\widetilde{A}_{2},\ldots,
\widetilde{A}_{p}\right)$, the size of each diagonal block $\widetilde{A}_{i}$ is upper bounded by $d^{\star}$, where $d^{\star}\leq Cd^{\alpha}$ for some known constant $C>0$ and $\alpha\in(0,1)$, but $d^{\star}$ is unknown.

In assumption (i) of 1, 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 1 uses the constant $\alpha$ to quantify the extent to which features are clustered.

The procedure for estimating $\left\lVert E\right\rVert_{\infty}$ in 1 involves an estimation of the variance $\sigma^{2}$, which is adapted from Comminges et al. (), and then multiply a constant multiple of $\sqrt{\log{d}}$.
We present the details in Algorithm 4:

Figure: Algorithm 4 Estimation of $\left\lVert E\right\rVert_{\infty}$ in 1

In the next theorem, we provide a high probability approximation bound, as well as runtime analysis, for Algorithm 3 using the threshold $\bar{\varepsilon}$ found via Algorithm 4, suggesting that in 1 one can efficiently solve SPCA by exploiting the hidden sparse structure.

###### Proposition 2 .

Consider 1, and denote by $\bar{\varepsilon}$ the output of Algorithm 4 with input tuple $(A,C,\alpha,u)$.
Let $k\leq d$ be a positive integer, and assume that $d$ satisfies that $d^{1-\alpha}>C_{0}\cdot(C+1)\log((C+1)/[u(u+1)]+1)$
for some large enough absolute constant $C_{0}>0$.
Then, the following holds with probability at least $1-2d^{-1}-\exp\{2\log{d}-{d^{1+\alpha}}/{(4C+4)}\}$:

- (i)
Denote by $\widetilde{x}^{\star}$ the optimal solution to SPCA with input $(\widetilde{A},k)$.
For an (approximation) algorithm $\mathcal{A}$ with multiplicative factor $m(k,d)\geq 1$ and additive error $a(k,d)\geq 0$, where the functions $m$ and $a$ is non-deceasing with respect to $d$, the output $y$ of Algorithm 3 with input tuple $(A,k,\bar{\varepsilon},\mathcal{A})$ satisfies that
$\displaystyle y^{\top}\widetilde{A}y\geq\frac{(\widetilde{x}^{\star})^{\top}
\widetilde{A}\widetilde{x}^{\star}}{m\left(k,d^{\star}\right)}-a\left(k,d^{
\star}\right)-\left(1+\frac{1}{m\left(k,d^{\star}\right)}\right)\cdot 2k\bar{
\varepsilon}.$
- (ii)
If, in addition, the time complexity of using $\mathcal{A}$ 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(k,1)=1$.
Then, the runtime of Algorithm 3 with input tuples $(A,k,\varepsilon,\mathcal{A})$ is at most
$\displaystyle\mathcal{O}\left(\left\lceil\frac{d}{d^{\star}}\right\rceil\cdot g
\left(k,d^{\star}\right)+d^{2}\right)$

###### Model 1 .

###### Proposition 2 .

### 4.2 A practical learning approach

In this section, we describe a practical algorithm to estimate $\varepsilon$ without prior statistical assumptions while simultaneously providing an approximate solution to SPCA.

Recall from Algorithm 3 that for a given threshold $\varepsilon$, a thresholded matrix $A^{\varepsilon}$ is generated via Algorithm 1 from the original matrix $A$.
This matrix $A^{\varepsilon}$ is then utilized to approximate SPCA using a specified algorithm $\mathcal{A}$.
Our proposed approach in this section involves initially selecting an $\varepsilon$ value within a predefined range, executing Algorithm 3 with this $\varepsilon$, and subsequently adjusting $\varepsilon$ in the following iterations through a binary search process.
Since efficiency is of our major interest; thus, we ask users to provide a parameter $d_{0}$ that indicates the maximum dimension they are willing to compute due to computational resource limit.
We defer the detailed description of this binary-search based algorithm to Algorithm 5.
For clarity, let us denote by $\mathsf{OPT}(\varepsilon)$ the objective value of $(x^{\varepsilon})^{\top}Ax^{\varepsilon}$, where $x^{\varepsilon}$ is the output of Algorithm 3 with inputs $(A,k,\varepsilon,\mathcal{A})$.

Figure: Algorithm 5 A practical version of Algorithm 3 without a given threshold

We have the characterization of approximation error and runtime complexity of Algorithm 5:

###### Theorem 4 .

Let $A\in{\mathbb{R}}^{d\times d}$ and let $k$ be a positive integer.
Suppose that an algorithm $\mathcal{A}$ for SPCA with input $(A,k)$ finds an approximate solution $x\in{\mathbb{R}}^{d}$ to SPCA such that has multiplicative factor $m(k,d)\geq 1$ and additive error $a(k,d)\geq 0$, with $m(k,d)$ and $a(k,d)$ being non-decreasing with respect to $d$.
Suppose that Algorithm 5 with input $(A,k,\delta,\mathcal{A},d_{0})$ terminates with $\mathsf{OPT}(\varepsilon^{\star})$.
Then one has

$$ $\displaystyle\mathsf{OPT}(\varepsilon^{\star})\geq\frac{1}{m(k,d_{0})}\cdot \mathsf{OPT}(0)-a(k,d_{0})-\left(1+\frac{1}{m(k,d_{0})}\right)\cdot k \varepsilon^{\star}.$ $$

Suppose that for $\mathcal{A}$, the time complexity of using $\mathcal{A}$ 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 $g(k,1)=1$.
The runtime for Algorithm 5 is at most

$$ $\displaystyle\mathcal{O}\left(\log\left(\frac{\left\lVert A\right\rVert_{ \infty}}{\delta}\right)\cdot\left(\left\lceil\frac{d}{d_{0}}\right\rceil\cdot g \left(k,d_{0}\right)+d^{2}\right)\right)$ $$

###### Remark 5 .

In this remark, we discuss the practical advantages of Algorithm 5:

- •
The first direct advantage is that Algorithm 5 does not rely on any assumptions, enabling users to apply the framework 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 $\bar{\varepsilon}$ for binary search, they can refine the search interval in Algorithm 5 by setting the lower and upper bounds to $L\vcentcolon=a\bar{\varepsilon}$ and $U\vcentcolon=b\bar{\varepsilon}$ for some proper $0\leq a<b$ to improve the efficiency by a speedup factor of
$\mathcal{O}(\log(\left\lVert A\right\rVert_{\infty}/\delta)/\log((b-a)\bar{
\varepsilon}/\delta))$.
- •
Secondly, as detailed in Theorem 4, is that Algorithm 5 achieves a similar trade-off between computational efficiency and approximation error as Algorithm 1, which are proposed in Theorems 3 and 1. The time complexity of Algorithm 5, compared to Algorithm 3, increases only by a logarithmic factor due to the binary search step, thus preserving its efficiency.
- •
The third advantage is that this algorithm considers computational resources, ensuring that each invocation of Algorithm 3 remains efficient.

###### Theorem 4 .

###### Remark 5 .

## 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. (), LymphomaCov1 from Alizadeh et al. (), Reddit1500 and Reddit2000 from Dey et al. (), LeukemiaCov, LymphomaCov2, and ProstateCov from Dettling (), ArceneCov and DorotheaCov from Guyon et al. (), and GLI85Cov and GLABRA180Cov from Li ().
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 $A$ are dense and hence are not block-diagonal matrices, as discussed in Section C.1.

Baselines. We use the state-of-the-art exact algorithm Branch-and-Bound (Berk and Bertsimas, ) and the state-of-the-art polynomial-time approximation algorithm (in terms of multiplicative factor) Chan’s algorithm (Chan et al., ) as our baselines.

Notation.
We calculate the *approximation error* between two methods using the formula $(\textup{obj}_{\textup{Base}}-\textup{obj}_{\textup{Ours}})/\textup{obj}_{
\textup{Base}}\times 100\%$, where $\textup{obj}_{\textup{Base}}$ is the objective value obtained by a baseline algorithm, while $\textup{obj}_{\textup{Ours}}$ 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 5 integrated with the baseline algorithm.

Results. In Table 3, we report the performance of Algorithm 5 when integrated with Branch-and-Bound algorithm.
We summarized the performance results for $k=2,5,10,15$ across various datasets.
For each dataset, we report the average approximation error, average speedup factors, and their corresponding standard deviations.
Notably, the average speedup factor across all datasets is 93.77, demonstrating significant computational efficiency.
Our theoretical analysis anticipates a trade-off between speedups and approximation errors; however, the empirical performance of our framework in terms of approximation accuracy is remarkably strong, with an average error of just 2.15% across all datasets.

In Table 4 we report the performance of Algorithm 5, in larger datasets with larger $k$’s, when integrated with Chan’s algorithm (Chan et al., ).
We summarized the performance results for $k=200,500,1000,2000$ across various datasets.
Chan’s algorithm terminates in time $\mathcal{O}({d^{3}})$ and returns an approximate solution to SPCA with a multiplicative factor $\min\{\sqrt{k},d^{1/3}\}$.
The average speedup factor is 6.77, and grows up to 14.70 as the dimension of the dataset grows, while the average error across all four datasets remains below 1.35%.
Remarkably, for the DorotheaCov dataset, our framework integrated with Chan’s algorithm achieves negative approximation errors, 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.
While generally successful, there are rare instances where our framework may underperform. However, we have found that fine-tuning Algorithm 5 readily enhances performance in these instances.
For a more detailed discussion, additional information on settings, and further empirical results, we refer readers to Appendix C.

**Table 3: Summary of average approximation errors and average speedup factors for each dataset, compared to Branch-and-Bound. Standard deviations are reported in the parentheses.**
| Dataset | Avg. Error (Std. Dev) | Avg. Speedup (Std. Dev) |
| --- | --- | --- |
| CovColonCov | 0.00 (0.00) | 158.15 (179.88) |
| LymphomaCov1 | 0.53 (0.77) | 239.14 (336.83) |
| Reddit1500 | 0.05 (0.05) | 59.95 (94.12) |
| Reddit2000 | 0.21 (0.28) | 90.37 (51.84) |
| LeukemiaCov | 1.50 (2.28) | 117.84 (70.87) |
| LymphomaCov2 | 11.99 (16.50) | 41.52 (46.43) |
| ProstateCov | 0.39 (0.78) | 38.41 (24.77) |
| ArceneCov | 2.55 (4.68) | 4.77 (5.02) |
| Overall | 2.15 (6.64) | 93.77 (145.74) |

**Table 4: Summary of average approximation errors and average speedup factors for each dataset, compared to Chan’s algorithm. Standard deviations are reported in the parentheses. Datasets are sorted in ascending order according to their dimensions.**
| Dataset | Avg. Error (Std. Dev) | Avg. Speedup (Std. Dev) |
| --- | --- | --- |
| ArceneCov | 1.35 (1.24) | 1.30 (0.08) |
| GLI85Cov | 0.25 (1.85) | 2.97 (0.28) |
| GLABRA180Cov | 0.76 (0.33) | 8.12 (0.87) |
| DorotheaCov | -0.87 (0.68) | 14.70 (0.66) |
| Overall | 0.37 (1.35) | 6.77 (5.42) |

## 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 a plug-and-play fashion.
Our approach leverages matrix block-diagonalization, 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 of 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.

## Appendix 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 a selective overview of mainstream methods.
Algorithms require worst-case exponential runtimes, including mixed-integer semidefinite programs [Bertsimas et al., , Li and Xie, , Cory-Wright and Pauphilet, ], mixed-integer quadratic programs [Dey et al., ], non-convex programs [Zou et al., , Dey et al., ], and sub-exponential algorithms [Ding et al., ].
On the other hand, there also exist a number of algorithms that takes polynomial runtime to approximately solving Sparse PCA.
To the best of our knowledge, the state-of-the-art approximation algorithm is proposed by Chan et al. [], in which an algorithm with multiplicative factor $\min\{\sqrt{k},d^{-1/3}\}$ is proposed.
There also exists tons of other (approximation) algorithms that have polynomial runtime, including semidefinite relaxations [Chowdhury et al., ], greedy algorithms [Li and Xie, ], low-rank approximations [Papailiopoulos et al., ], and fixed-rank Sparse PCA [Del Pia, ].
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 and Wainwright, , d’Orsi et al., ] and covariance thresholding [Deshp et al., ].
Finally, there are also local approximation algorithms designed for Sparse PCA, and iterative algorithms are proposed [Yuan and Zhang, , Hager et al., ].

We also notice that several methods involving or similar to matrix block-diagonalization are present in the literature.
For instance, Feng et al. [] employ block-diagonalization for subspace clustering and segmentation problems by mapping data matrices into block-diagonal matrices through linear transformations. More recently, Han et al. [] introduce the sortedLSH algorithm to sparsify the softmax matrix, thus reducing the computational time of the attention matrix through random projection-based bucketing of “block indices”.
We remark that in these studies, the underlying problems of interests, as well as the techniques, are different from ours, and the techniques developed therein cannot be directly applied to SPCA.
A more related work [Devijver and Gallopin, ] 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 a 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 to determine the optimal threshold $\varepsilon$.

## Appendix B Deferred proofs

In this section, we provide proofs that we defer in the paper.

### B.1 Proofs in Section 3.2

See 1

###### Proof of Lemma 1 .

It is clear that Algorithm 1 has time complexity $\mathcal{O}(d^{2})$, and it suffices the show that the output $\widetilde{A}$ of Algorithm 1 with input $(A,\varepsilon)$ satisfies $\textup{lbs}\left(\widetilde{A}\right)=\operatorname*{int\,dim}(A,\varepsilon)$.

Suppose that $A^{\star}$ is an $\varepsilon$-approximation of $A$, and satisfies $\textup{lbs}\left(A^{\star}\right)=\operatorname*{int\,dim}(A,\varepsilon)$.
We intend to show $\textup{lbs}\left(\widetilde{A}\right)\leq\textup{lbs}\left(A^{\star}\right)$.
WLOG we assume that $A^{\star}\in\mathcal{B}(A,\varepsilon)$ is a block-diagonal matrix with blocks $A_{1}^{\star},A_{2}^{\star},\ldots,A_{p}^{\star}$, i.e., $A^{\star}=\mathsf{diag}\left(A_{1}^{\star},A_{2}^{\star},\ldots,A_{p}^{\star}\right)$.
Write $A=\begin{pmatrix}A_{11}&A_{12}&\ldots&A_{1p}\\
A_{21}&A_{22}&\ldots&A_{2p}\\
\vdots&\vdots&\vdots&\vdots\\
A_{p1}&A_{p2}&\ldots&A_{pp}\end{pmatrix}$ and $\widetilde{A}=\begin{pmatrix}\widetilde{A}_{11}&\widetilde{A}_{12}&\ldots&
\widetilde{A}_{1p}\\
\widetilde{A}_{21}&\widetilde{A}_{22}&\ldots&\widetilde{A}_{2p}\\
\vdots&\vdots&\vdots&\vdots\\
\widetilde{A}_{p1}&\widetilde{A}_{p2}&\ldots&\widetilde{A}_{pp}\end{pmatrix}$, where $A_{ii}$, $\widetilde{A}_{ii}$, and $A^{\star}_{i}$ have the same number of rows and columns, by the definition of $\mathcal{B}(A,\varepsilon)$, it is clear that $\left\lVert A_{ij}\right\rVert_{\infty}\leq\varepsilon$, for any $1\leq i\neq j\leq p$.
This further implies that for any $1\leq i\neq j\leq p$, $\widetilde{A}_{ij}$ is a zero matrix by the procedure described in Algorithm 1, and therefore it is clear that $\textup{lbs}\left(\widetilde{A}\right)\leq\textup{lbs}\left(A^{\star}\right)$.
By minimality of $A^{\star}$, we obtain that $\textup{lbs}\left(\widetilde{A}\right)=\textup{lbs}\left(A^{\star}\right)$.
∎

###### Proof of Lemma 1 .

#### B.1.1 Proofs in Section 3.2.1

See 1

###### Proof of Theorem 1 .

Note that

$$ $\displaystyle\left\lvert(x^{\star})^{\top}Ax^{\star}-\widetilde{x}^{\top}A \widetilde{x}\right\rvert$ $\displaystyle\leq\left\lvert(x^{\star})^{\top}Ax^{\star}-\widetilde{x}^{\top}A ^{\varepsilon}\widetilde{x}\right\rvert+\left\lvert\widetilde{x}^{\top}(A^{ \varepsilon}-A)\widetilde{x}\right\rvert$ $\displaystyle\leq\left\lvert(x^{\star})^{\top}Ax^{\star}-\widetilde{x}^{\top}A ^{\varepsilon}\widetilde{x}\right\rvert+k\cdot\left\lVert A-A^{\varepsilon} \right\rVert_{\infty},$ $$

where the last inequality follows from Holder’s inequality and the fact that $\left\lVert\widetilde{x}\right\rVert_{1}\leq\sqrt{k}\cdot\left\lVert\widetilde
{x}\right\rVert_{2}=\sqrt{k}$.
Therefore, it suffices to show that

$$ $\displaystyle\left\lvert(x^{\star})^{\top}Ax^{\star}-\widetilde{x}^{\top}A^{ \varepsilon}\widetilde{x}\right\rvert\leq k\cdot\left\lVert A-A^{\varepsilon} \right\rVert_{\infty}.$ $$

Observe that

$$ $\displaystyle\left\lvert(x^{\star})^{\top}Ax^{\star}-\widetilde{x}^{\top}A^{ \varepsilon}\widetilde{x}\right\rvert$ $\displaystyle=\left\lvert\max_{\begin{subarray}{c}\left\lVert x\right\rVert_{2 }=1,\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}x^{\top}Ax-\max_{\begin{ subarray}{c}\left\lVert x\right\rVert_{2}=1,\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}x^{\top}A^{\varepsilon}x\right\rvert$ $\displaystyle\leq\max_{\begin{subarray}{c}\left\lVert x\right\rVert_{2}=1,\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}\left\lvert x^{\top}(A^{ \varepsilon}-A)x\right\rvert$ $\displaystyle\leq k\cdot\left\lVert A-A^{\varepsilon}\right\rVert_{\infty},$ $$

where the first inequality follows from the following facts:

$$ $\displaystyle\quad\max_{\begin{subarray}{c}\left\lVert x\right\rVert_{2}=1,\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}x^{\top}Ax-\max_{\begin{ subarray}{c}\left\lVert x\right\rVert_{2}=1,\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}x^{\top}A^{\varepsilon}x\leq \max_{\begin{subarray}{c}\left\lVert x\right\rVert_{2}=1,\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}x^{\top}(A-A^{\varepsilon})x$ $\displaystyle\quad\max_{\begin{subarray}{c}\left\lVert x\right\rVert_{2}=1,\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}x^{\top}A^{\varepsilon}x- \max_{\begin{subarray}{c}\left\lVert x\right\rVert_{2}=1,\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}x^{\top}Ax\leq\max_{\begin{ subarray}{c}\left\lVert x\right\rVert_{2}=1,\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}x^{\top}(A^{\varepsilon}-A)x$ $$

∎

See 2

###### Proof of Theorem 2 .

It is clear that $\mathsf{OPT}\geq\max_{i\in[p]}\mathsf{OPT}_{i}$.
It remains to show the reverse, i.e., $\mathsf{OPT}\leq\max_{i\in[p]}\mathsf{OPT}_{i}$.

Suppose not, we have $\mathsf{OPT}>\max_{i\in[p]}\mathsf{OPT}_{i}$.
Then, there must exist principal submatrices $\widetilde{A}_{1}^{\prime},\widetilde{A}_{2}^{\prime},\ldots,\widetilde{A}_{p}
^{\prime}$ contained in $\widetilde{A}_{1},\widetilde{A}_{2},\ldots,\widetilde{A}_{p}$ such that $\mathsf{diag}(\widetilde{A}_{1}^{\prime},\widetilde{A}_{2}^{\prime},\ldots,
\widetilde{A}_{p}^{\prime})$ is a $k\times k$ matrix, and

$$ $\displaystyle\mathsf{OPT}=\max_{\begin{subarray}{c}\left\lVert x\right\rVert_{ 2}=1\end{subarray}}x^{\top}\mathsf{diag}(\widetilde{A}_{1}^{\prime},\widetilde {A}_{2}^{\prime},\ldots,\widetilde{A}_{p}^{\prime})x>\max_{i\in[p]}\mathsf{OPT }_{i}.$ $$

This gives a contradiction, as

$$ $\displaystyle\max_{i\in[p]}\mathsf{OPT}_{i}=\max_{i\in[p]}\max_{\begin{ subarray}{c}\left\lVert x\right\rVert_{2}=1,\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}x^{\top}\widetilde{A}_{i}x \geq\max_{\begin{subarray}{c}\left\lVert x\right\rVert_{2}=1\end{subarray}}x^{ \top}\mathsf{diag}(\widetilde{A}_{1}^{\prime},\widetilde{A}_{2}^{\prime}, \ldots,\widetilde{A}_{p}^{\prime})x,$ $$

where the inequality follows from the fact that

$$ $\displaystyle\max_{\begin{subarray}{c}\left\lVert x\right\rVert_{2}=1\end{ subarray}}x^{\top}\mathsf{diag}(\widetilde{A}_{1}^{\prime},\widetilde{A}_{2}^{ \prime},\ldots,\widetilde{A}_{p}^{\prime})x=\max_{i\in[p]}\max_{\begin{ subarray}{c}\left\lVert x\right\rVert_{2}=1\end{subarray}}x^{\top}\widetilde{A }_{i}^{\prime}x=\max_{i\in[p]}\max_{\begin{subarray}{c}\left\lVert x\right \rVert_{2}=1\\ \left\lVert x\right\rVert_{0}\leq k\end{subarray}}x^{\top}\widetilde{A}_{i}^{ \prime}x.$ $$

∎

See 3

###### Proof of Theorem 3 .

In the proof, we use the same notation as in Algorithm 3.
We assume WLOG that the thresholded matrix $A^{\varepsilon}$ is a block-diagonal matrix, i.e., $A^{\varepsilon}=\mathsf{diag}\left(\widetilde{A}_{1},\widetilde{A}_{2},\dots,
\widetilde{A}_{p}\right)$.
Denote $\widetilde{x}^{\star}\in{\mathbb{R}}^{d}$ to be an optimal solution to SPCA with input $(A^{\varepsilon},k)$.
By Theorem 2, one can assume WLOG that $\widetilde{x}^{\star}$ has zero entries for indices greater than $d_{1}$.
In other words, $\widetilde{x}^{\star}$ is “found” in the block $A_{1}$.
Recall that we denote $x_{i}$ the solution found by $\mathcal{A}$ with input $\left(\widetilde{A}_{i},k\right)$.
Then, it is clear that

$$ $\displaystyle\max_{i\in[p]}x_{i}^{\top}\widetilde{A}_{i}x_{i}$ $\displaystyle=x_{1}^{\top}\widetilde{A}_{1}x_{1}$ $\displaystyle\geq\frac{1}{m(k,d_{1})}\cdot(\widetilde{x}^{\star})^{\top}A^{ \varepsilon}\widetilde{x}^{\star}-a(k,d_{1})$ $\displaystyle\geq\frac{1}{m(k,\operatorname*{int\,dim}(A,\varepsilon))}\cdot( \widetilde{x}^{\star})^{\top}A^{\varepsilon}\widetilde{x}^{\star}-a(k, \operatorname*{int\,dim}(A,\varepsilon)).$ $$

Recall that $y_{i}$ is the reconstructed solution.
By the proof of Theorem 1, one obtains that

$$ $\displaystyle\quad\max_{i\in[p]}y_{i}^{\top}A^{\varepsilon}y_{i}$ (1) $\displaystyle=\max_{i\in[p]}x_{i}^{\top}\widetilde{A}_{i}x_{i}$ $\displaystyle\geq\frac{1}{m(k,\operatorname*{int\,dim}(A,\varepsilon))}\cdot( \widetilde{x}^{\star})^{\top}A^{\varepsilon}\widetilde{x}^{\star}-a(k, \operatorname*{int\,dim}(A,\varepsilon))$ $\displaystyle\geq\frac{1}{m(k,\operatorname*{int\,dim}(A,\varepsilon))}\cdot \left((x^{\star})^{\top}Ax^{\star}-k\left\lVert A-A^{\varepsilon}\right\rVert_ {\infty}\right)-a(k,\operatorname*{int\,dim}(A,\varepsilon))$ $$

Finally, since

$$ $\displaystyle\left\lvert y_{i}^{\top}A^{\varepsilon}y_{i}-y_{i}^{\top}Ay_{i} \right\rvert\leq k\left\lVert A-A^{\varepsilon}\right\rVert_{\infty}\leq\varepsilon,$ $$

combining with (1), one obtains the desired result

$$ $\displaystyle\max_{i\in[p]}y_{i}^{\top}Ay_{i}\geq\frac{1}{m(k,\operatorname*{ int\,dim}(A,\varepsilon))}\cdot\left((x^{\star})^{\top}Ax^{\star}-k\varepsilon \right)-a(k,\operatorname*{int\,dim}(A,\varepsilon))-k\varepsilon.$ $$

∎

###### Proof of Theorem 1 .

###### Proof of Theorem 2 .

###### Proof of Theorem 3 .

#### B.1.2 Proofs in Section 3.2.2

See 1

###### Proof of Proposition 1 .

Suppose that the largest block size of $\widetilde{A}\in\mathcal{B}(A,\varepsilon)$ is equal to $\operatorname*{int\,dim}(A,\varepsilon)$, and $\widetilde{A}$ can be sorted to a block-diagonal matrix with $p$ blocks $\widetilde{A}_{1},\widetilde{A}_{2},\ldots,\widetilde{A}_{p}$ and block sizes $d_{1},d_{2},\ldots,d_{p}$, respectively.
WLOG we assume that $p<d$, otherwise the statement holds trivially.
It is clear that the time complexity of running line 5 - 8 in Algorithm 3 is upper bounded by a constant multiple of

$$ $\displaystyle\max_{\begin{subarray}{c}1\leq d_{i}\leq\operatorname*{int\,dim}( A,\varepsilon)\\ \sum_{i=1}^{m}d_{i}=d\end{subarray}}\sum_{i=1}^{p}g(k,d_{k}).$ (2) $$

Since $k$ is fixed, the optimization problem (2) is maximizing a convex objective function over a polytope

$$ $\displaystyle P\vcentcolon=\left\{(d_{1},d_{2},\ldots,d_{p})\in{\mathbb{R}}^{p }:1\leq d_{i}\leq\operatorname*{int\,dim}(A,\varepsilon),\ \forall i\in[p],\ \sum_{i=1}^{p}d_{i}=d\right\}.$ $$

Therefore, the optimum to (2) is obtained at an extreme point of $P$.
It is clear that an extreme point $(d_{1}^{\star},d_{2}^{\star},\ldots,d_{p}^{\star})$ of $P$ is active at a linearly independent system, i.e., among all $d_{i}^{\star}$’s, there are $p-1$ of them must be equal to $1$ or $\operatorname*{int\,dim}(A,\varepsilon)$, and the last one is taken such that the euqality $\sum_{i=1}^{p}d_{i}^{\star}=d$ holds.
Since there could be at most $\lfloor d/\operatorname*{int\,dim}(A,\varepsilon)\rfloor$ many blocks in $\widetilde{A}$ with block size equal to $\operatorname*{int\,dim}(A,\varepsilon)$, and that $g(k,d)$ is increasing with respect to $d$, we obtain that the optimum of the optimization problem (2) is upper bounded by

$$ $\displaystyle\quad\left\lceil\frac{d}{\operatorname*{int\,dim}(A,\varepsilon)} \right\rceil\cdot g\left(k,\operatorname*{int\,dim}(A,\varepsilon)\right)+ \left(m-\left\lfloor\frac{d}{\operatorname*{int\,dim}(A,\varepsilon)}\right \rfloor\right)$ $\displaystyle<\left\lceil\frac{d}{\operatorname*{int\,dim}(A,\varepsilon)} \right\rceil\cdot g\left(k,\operatorname*{int\,dim}(A,\varepsilon)\right)+d.$ $$

We are then done by noticing the fact that lines 3 - 4 in Algorithm 3 have a runtime bounded by $\mathcal{O}(d^{2})$.
∎

###### Proof of Proposition 1 .

### B.2 Proofs in Section 4.1

To establish a connection with robust statistics, we begin by expressing each element of the input matrix $A$ in 1 as follows:

$$ $\displaystyle A_{ij}=\widetilde{A}{ij}+E{ij}\vcentcolon=\widetilde{A}{ij}+ \sigma^{2}Z{ij},\quad i\leq j,$ $$

where $Z_{ij}$ is an i.i.d. centered $(\varrho^{2}/\sigma^{2})$-sub-Gaussian random variable with unit variance.
Considering that $\widetilde{A}$ is block-diagonal and, according to (ii) in 1, has at most $\lceil d/d^{\star}\rceil\cdot(d^{\star})^{2}=\mathcal{O}(d^{1+\alpha})$ nonzero entries, the methodology proposed by Comminges et al. [] can be employed to estimate $\sigma^{2}$ accurately using Algorithm 4.

Inspired by the methods proposed in Comminges et al. [], we show that the number $\bar{\sigma}$ found in Algorithm 4 could provide an constant approximation ratio to the true variance $\sigma^{2}$ in 1 with high probability:

###### Proposition 3 (adapted from Proposition 1 in Deshp et al. [ 2016 ] ) .

Consider 1, and denote $\mathcal{P}_{Z}$ to be the distribution of the random variable $E_{ij}/\sigma^{2}$.
Let $\bar{\sigma}^{2}$ be the same number found in Algorithm 4 with input tuple $(A,C,\alpha)$.
There exist an absolute constant $c>0$ such that for the dimension $d$ of the input matrix $A$ satisfying

$$ $\displaystyle\log{d}>\log\left[(2C+2)\cdot\left(\frac{\log{(16C+16)}}{\frac{c \sigma^{2}}{2\varrho^{2}}\cdot\left(\frac{\sigma^{2}}{2\varrho^{2}}+1\right)}+ 1\right)\right]/(1-\alpha)$ $$

we have that for distribution $\mathcal{P}_{Z}$ of $Z_{ij}$, the variance $\sigma^{2}$, one has a uniform lower bound on the accuracy of estimation:

$$ $\displaystyle\inf_{\mathcal{P}_{Z}}\inf_{\sigma>0}\inf_{\|\widetilde{A}\|_{0} \leq Cd^{1+\alpha}}\mathbb{P}_{\mathcal{P}_{Z},\sigma,\widetilde{A}}\left( \frac{1}{2}\leq\frac{\bar{\sigma}^{2}}{\sigma^{2}}\leq\frac{3}{2}\right)$ $\displaystyle\geq 1-\exp\left\{-\frac{d^{\alpha+1}}{4C+4}\right\}.$ $$

###### Proof of Proposition 3 .

Recall that we write

$$ $\displaystyle A_{ij}=\widetilde{A}{ij}+E{ij}\vcentcolon=\widetilde{A}{ij}+ \sigma^{2}Z{ij},\ i\leq j,$ $$

where $Z_{ij}$ is an i.i.d. centered $(\varrho^{2}/\sigma^{2})$-sub-Gaussian random variable with unit variance.
Since $d^{\star}\leq Cd^{\alpha}$ for some $\alpha\in(0,1)$, WLOG we assume that $d^{\star}=Cd^{\alpha}$ for some $C>0$, and for simplicity we assume that $d^{1+\alpha}$ and $C$ are all integers to avoid countless discussions of integrality.
By Algorithm 4, one divides the index set $\{(i,j)\in[d]\times[d]:i\leq j\}$ into $m=(2C+2)d^{1+\alpha}$ disjoint subsets $B_{1},B_{2},\ldots,B_{m}$, and each $B_{i}$ has cardinality at lease $k:=\lfloor(d^{2}+d)/(2m)\rfloor$.
Denote the set $J$ to be the set of indices $i$ such that the set $B_{i}$ contains only indices such that for any $(k,l)\in B_{i}$, $\widetilde{A}_{kl}=0$.
It is clear that $|J|$ is lower bounded by $m-Cd^{\alpha+1}=(2C+1)d^{1+\alpha}$.

By Bernstein’s inequality (see, e.g., Corollary 2.8.3 in Vershynin []), it is clear that there exists an absolute constant $c>0$ such that

$$ $\displaystyle\mathbb{P}\left(\left\lvert\frac{1}{|B_{i}|}\sum_{(k,l)\in B_{i}} E_{kl}^{2}-\sigma^{2}\right\rvert>\frac{\sigma^{2}}{2}\right)$ $\displaystyle=\mathbb{P}\left(\left\lvert\frac{1}{|B_{i}|}\sum_{(k,l)\in B_{i} }Z_{kl}^{2}-1\right\rvert>\frac{1}{2}\right)$ $\displaystyle\leq 2\exp\left\{-c\min\left(\frac{\sigma^{4}}{4\varrho^{4}}, \frac{\sigma^{2}}{2\varrho^{2}}\right)|B_{i}|\right\}$ $\displaystyle\leq 2\exp\left\{-c\min\left(\frac{\sigma^{4}}{4\varrho^{4}}, \frac{\sigma^{2}}{2\varrho^{2}}\right)\left\lfloor\frac{d^{2}+d}{2m}\right \rfloor\right\}.$ $$

Then, we denote $I\vcentcolon=[\sigma^{2}/2,3\sigma^{2}/2]$, denote $A_{i}$ to be the random event $\left\{\frac{1}{|B_{i}|}\sum_{(k,l)\in B_{i}}A_{kl}^{2}\not\in I\right\}$, and denote $\mathds{1}_{A_{i}}$ to be the indicator function of $A_{i}$.
One obtains that

$$ $\displaystyle\mathbb{P}\left(\bar{\sigma}^{2}\not\in I\right)$ $\displaystyle\leq\mathbb{P}\left(\sum_{i=1}^{m}\mathds{1}_{A_{i}}\geq\frac{m}{ 2}\right)$ $\displaystyle=\mathbb{P}\left(\sum_{i\in J}\mathds{1}_{A_{i}}\geq\frac{m}{2}- \sum_{i\not\in J}\mathds{1}_{A_{i}}\right)$ $\displaystyle\leq\mathbb{P}\left(\sum_{i\in J}\mathds{1}_{A_{i}}\geq\frac{m}{2 }-Cd^{\alpha+1}\right)$ $$

Denote

$$ $\displaystyle t\vcentcolon=\frac{m}{2}-Cd^{\alpha+1}-2|J|\exp\left\{-c\min \left(\frac{\sigma^{4}}{4\varrho^{4}},\frac{\sigma^{2}}{2\varrho^{2}}\right) \left\lfloor\frac{d^{2}+d}{2m}\right\rfloor\right\},$ $$

and by Hoeffding’s inequality, one obtains that

$$ $\displaystyle\quad\mathbb{P}\left(\sum_{i\in J}\mathds{1}_{A_{i}}\geq\frac{m}{ 2}-Cd^{\alpha+1}\right)$ $\displaystyle\leq\mathbb{P}\left(\sum_{i\in J}(\mathds{1}_{A_{i}}-\mathbb{E} \mathds{1}_{A_{i}})\geq\frac{m}{2}-Cd^{\alpha+1}-2|J|\exp\left\{-c\min\left( \frac{\sigma^{4}}{4\varrho^{4}},\frac{\sigma^{2}}{2\varrho^{2}}\right)\left \lfloor\frac{d^{2}+d}{2m}\right\rfloor\right\}\right)$ $\displaystyle=\mathbb{P}\left(\sum_{i\in J}(\mathds{1}_{A_{i}}-\mathbb{E} \mathds{1}_{A_{i}})\geq t\right)$ $\displaystyle\leq\exp\left\{-\frac{2t^{2}}{|J|}\right\}$ $$

Since $m=(2C+2)d^{\alpha+1}$, it is clear that $m\geq|J|\geq(2C+1)d^{\alpha+1}$, and

$$ $\displaystyle t$ $\displaystyle\geq d^{\alpha+1}\left(1-2(2C+2)\exp\left\{-c\min\left(\frac{ \sigma^{4}}{4\varrho^{4}},\frac{\sigma^{2}}{2\varrho^{2}}\right)\lfloor d^{1- \alpha}/(4C+4)\rfloor\right\}\right)\geq d^{\alpha+1}\cdot\frac{1}{2},$ $$

where the last inequality follows from the fact that

$$ $\displaystyle\log{d}>\log\left[(2C+2)\cdot\left(\frac{\log{(16C+16)}}{\frac{c \sigma^{2}}{2\varrho^{2}}\cdot\left(\frac{\sigma^{2}}{2\varrho^{2}}+1\right)}+ 1\right)\right]/(1-\alpha).$ $$

Hence, the probability of $\bar{\sigma}^{2}\not\in I$ is upper bounded by $\exp\{-d^{\alpha+1}/(4C+4)\}$.
∎

See 2

###### Proof of Proposition 2 .

We first show that, the probability that $\left\lVert E\right\rVert_{\infty}\leq\bar{\varepsilon}$ is at least $1-2d^{-1}-\exp\left\{2\log{d}-\frac{d^{1+\alpha}}{4C+4}\right\}$.
This comes directly from Proposition 3 that $\bar{\sigma}^{2}$ satisfies that $\frac{1}{2}\leq\frac{\bar{\sigma}^{2}}{\sigma^{2}}\leq\frac{3}{2}$ with probability at least $1-\exp\{-d^{1+\alpha}/(4C+4)\}$, the definition of $E_{kl}$ being i.i.d. $\varrho$-sub-Gaussian variables, and a union bound argument.
Indeed, define the random event $A\vcentcolon=\{\frac{1}{2}\leq\frac{\bar{\sigma}^{2}}{\sigma^{2}}\leq\frac{3}{
2}\}$, one can see that

$$ $\displaystyle\mathbb{P}\left(|E_{kl}|>\bar{\varepsilon}\right)$ $\displaystyle=\mathbb{P}\left(|E_{kl}|>\bar{\varepsilon}\mid A\right)\mathbb{P }(A)+\mathbb{P}\left(|E_{kl}|>\bar{\varepsilon}\mid A^{c}\right)\mathbb{P}(A^{ c})$ $\displaystyle=\mathbb{P}\left(|E_{kl}|>2\sqrt{3}u\bar{\sigma}\sqrt{\log{d}} \mid A\right)\mathbb{P}(A)+\mathbb{P}\left(|E_{kl}|>\bar{\varepsilon}\mid A^{c }\right)\mathbb{P}(A^{c})$ $\displaystyle\leq\mathbb{P}\left(|E_{kl}|>\sqrt{6}u\sigma\sqrt{\log{d}}\mid A \right)\mathbb{P}(A)+\mathbb{P}(A^{c})$ $\displaystyle\leq\mathbb{P}\left(|E_{kl}|>\sqrt{6}u\sigma\sqrt{\log{d}}\right) +\exp\left\{-\frac{d^{1+\alpha}}{4C+4}\right\}$ $\displaystyle\leq 2\exp\left\{-\frac{6u^{2}\sigma^{2}}{2\varrho^{2}}\log{d} \right\}+\exp\left\{-\frac{d^{1+\alpha}}{4C+4}\right\}$ $\displaystyle=2d^{-3}+\exp\left\{-\frac{d^{1+\alpha}}{4C+4}\right\}.$ $$

By union bound, one has

$$ $\displaystyle\mathbb{P}\left(\left\lVert E\right\rVert_{\infty}>\bar{ \varepsilon}\right)\leq d^{2}\cdot\mathbb{P}\left(|E_{kl}|>\bar{\varepsilon} \right)\leq 2d^{-1}+\exp\left\{2\log{d}-\frac{d^{1+\alpha}}{4C+4}\right\}.$ $$

This implies that, with high probability, $\operatorname*{int\,dim}\left(A,\bar{\varepsilon}\right)\leq d^{\star}$.
In the remainder of the proof, we will assume that the random event $\left\lVert E\right\rVert_{\infty}\leq\bar{\varepsilon}$ holds true.

We are now ready to show the desired approximation error.
We first notice that

$$ $\displaystyle\left\lVert\widetilde{A}-A^{\bar{\varepsilon}}\right\rVert_{ \infty}\leq\left\lVert\widetilde{A}-A\right\rVert_{\infty}+\left\lVert A-A^{ \bar{\varepsilon}}\right\rVert_{\infty}\leq 2\bar{\varepsilon}.$ $$

Follow a similar argument in the proof of Theorem 3, one can obtain the desired approximation error.

Finally, we obtain the desired runtime complexity via Proposition 1.
∎

###### Proposition 3 (adapted from Proposition 1 in Deshp et al. [ 2016 ] ) .

###### Proof of Proposition 3 .

###### Proof of Proposition 2 .

### B.3 Proofs in Section 4.2

See 4

###### Proof of Theorem 4 .

We only need to note that for such $\varepsilon^{\star}$, one has $\operatorname*{int\,dim}(A,\varepsilon^{\star})\leq d_{0}$, and $\left\lVert A-A^{\varepsilon^{\star}}\right\rVert_{\infty}\leq\varepsilon^{\star}$.
Also, for while loop in Algorithm 5, at most $\mathcal{O}(\log(\left\lVert A\right\rVert_{\infty}/\delta))$ iterations would be executed, accounting for the extra logarithmic factor.

The rest of the proof is identical to the proofs of Theorems 3 and 1.
∎

###### Proof of Theorem 4 .

## Appendix C Additional empirical settings and results

In this section, we report additional empirical settings in Section C.1 and detailed empirical results in Sections C.2 and C.3.

### C.1 Empirical Settings

Datasets. We summarize the dimensions and percentage of non-zero entries in the input matrices in Table 5.
It is important to 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% |
| LymphomaCov1 | 500 | 100% |
| Reddit1500 | 1500 | 5.30% |
| Reddit2000 | 2000 | 6.20% |
| LeukemiaCov | 3571 | 100% |
| LymphomaCov2 | 4026 | 100% |
| ProstateCov | 6033 | 100% |
| ArceneCov | 10000 | 100% |
| GLI85Cov | 22283 | 100% |
| GLABRA180Cov | 49151 | 100% |
| DorotheaCov | 100000 | 77.65% |

Parameters.
(i) When integrated Algorithm 5 with Branch-and-Bound algorithm: For experiments involving datasets CovColon [Alon et al., ], LymphomaCov1 [Alizadeh et al., ], Reddit1500 and Reddit2000 [Dey et al., ], we use the method described in Remark 5: first set $\bar{\varepsilon}\vcentcolon=\bar{\sigma}\sqrt{\log{d}}$, where $\bar{\sigma}$ is estimated using Algorithm 4 with $C\vcentcolon=1$, $\alpha\vcentcolon=0.5$, $u\vcentcolon=1$.
Additionally, we set parameters $a\vcentcolon=0.01$, $b\vcentcolon=1000$, $d_{0}\vcentcolon=40$, and $\delta\vcentcolon=0.1$.
For experiments with datasets LeukemiaCov, LymphomaCov2, and ProstateCov [Dettling, ], and ArceneCov [Guyon et al., ], we adjust the parameters to $\bar{\varepsilon}\vcentcolon=1$, $a\vcentcolon=0$, $b\vcentcolon=\left\lVert A\right\rVert_{\infty}$, $d_{0}\vcentcolon=40$, and $\delta\vcentcolon=0.01$.
This adjustment was made after observing that the $\bar{\varepsilon}$ estimated from Algorithm 5 resulted in significantly longer runtimes for these large-scale datasets, whereas a simpler, fixed value of $\bar{\varepsilon}$ proved to be more effective.
(ii) When integrated Algorithm 5 with Chan’s algorithm: We choose the parameter $\bar{\varepsilon}\vcentcolon=1$, $a\vcentcolon=0$, $b\vcentcolon=\left\lVert A\right\rVert_{\infty}$, $d_{0}\vcentcolon=2k$ (twice the sparsity constant), and $\delta\vcentcolon=0.01\cdot\left\lVert A\right\rVert_{\infty}$.

Solution improvement.
For experiments conducted, we add a very simple step in Algorithm 5 - solving PCA exactly after each $x^{\varepsilon}$ is obtained, for matrix $A_{S,S}$, where $S\vcentcolon=\operatorname{supp}(x^{\varepsilon})$.
This additional step is efficient for small $k$, and could further decrease the approximation error.

Compute resources.
We conducted all tests on a 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-Bound algorithm

In this section, we provide the detailed empirical results of our framework when integrated with Branch-and-Bound algorithm (previously summarized in Table 3), in Table 6.
We also summarized the average approximation errors and average speedup factors for with respect to $k$ in Table 7.
In Sections C.2.1 and C.2.2, we discuss improvements to Algorithm 5 to obtain smaller approximation errors and higher speedup factors.

**Table 6: Comparison of runtime and objective values on real-world datasets, between Branch-and-Bound and Algorithm 5 integrated with Branch-and-Bound. $d$ is the dimension of Sparse PCA problem, and $k$ is the sparsity constant. We set the time limit to 3600 seconds for each method. Approximation errors between two methods, and speedup factors are reported. We note that the objective value for ArceneCov is scaled by $10^{5}$.**
| Dataset | $d$ | $k$ | Branch-and-Bound | Ours | Error (%) | Speedup |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  |  | Time | Objective | Time(s) | Objective |  |  |
| CovColonCov | 500 | 3 | 8 | 1059 | 3 | 1059 | 0 | 2.67 |
| CovColonCov | 500 | 5 | 8 | 1646 | 3 | 1646 | 0 | 2.67 |
| CovColonCov | 500 | 10 | 3600 | 2641 | 12 | 2641 | 0 | 300 |
| CovColonCov | 500 | 15 | 3600 | 3496 | 11 | 3496 | 0 | 327.27 |
| LymphomaCov1 | 500 | 3 | 11 | 2701 | 4 | 2701 | 0 | 2.75 |
| LymphomaCov1 | 500 | 5 | 44 | 4300 | 5 | 4300 | 0 | 8.8 |
| LymphomaCov1 | 500 | 10 | 3600 | 6008 | 5 | 5911 | 1.64 | 720 |
| LymphomaCov1 | 500 | 15 | 3600 | 7628 | 16 | 7593 | 0.46 | 225 |
| Reddit1500 | 1500 | 3 | 75 | 946 | 13 | 946 | 0 | 5.77 |
| Reddit1500 | 1500 | 5 | 3600 | 980 | 18 | 980 | 0 | 200 |
| Reddit1500 | 1500 | 10 | 3600 | 1045 | 120 | 1044 | 0.10 | 30 |
| Reddit1500 | 1500 | 15 | 3600 | 1082 | 896 | 1081 | 0.09 | 4.02 |
| Reddit2000 | 2000 | 3 | 412 | 1311 | 25 | 1311 | 0 | 16.48 |
| Reddit2000 | 2000 | 5 | 3600 | 1397 | 29 | 1396 | 0.07 | 124.14 |
| Reddit2000 | 2000 | 10 | 3600 | 1523 | 28 | 1521 | 0.13 | 128.57 |
| Reddit2000 | 2000 | 15 | 3600 | 1605 | 39 | 1595 | 0.63 | 92.31 |
| LeukemiaCov | 3571 | 3 | 280 | 7.1 | 22 | 7.1 | 0 | 12.73 |
| LeukemiaCov | 3571 | 5 | 3600 | 10.2 | 22 | 10.2 | 0 | 163.64 |
| LeukemiaCov | 3571 | 10 | 3600 | 17.2 | 23 | 17.0 | 1.18 | 156.52 |
| LeukemiaCov | 3571 | 15 | 3600 | 21.8 | 26 | 20.8 | 4.81 | 138.46 |
| LymphomaCov2 | 4026 | 3 | 210 | 40.6 | 32 | 40.6 | 0 | 6.56 |
| LymphomaCov2 | 4026 | 5 | 144 | 63.6 | 32 | 63.6 | 0 | 4.5 |
| LymphomaCov2 | 4026 | 10 | 3600 | 78.2 | 35 | 69.2 | 13.01 | 102.86 |
| LymphomaCov2 | 4026 | 15 | 3600 | 93.4 | 69 | 69.2 | 34.97 | 52.17 |
| ProstateCov | 6033 | 3 | 102 | 8.1 | 64 | 8.1 | 0 | 1.59 |
| ProstateCov | 6033 | 5 | 3600 | 12.9 | 65 | 12.9 | 0 | 55.38 |
| ProstateCov | 6033 | 10 | 3600 | 24.3 | 75 | 24.3 | 0 | 48 |
| ProstateCov | 6033 | 15 | 3600 | 35.0 | 74 | 34.4 | 1.57 | 48.65 |
| ArceneCov | 10000 | 3 | 88 | 3.3 | 376 | 3.3 | 0 | 0.23 |
| ArceneCov | 10000 | 5 | 248 | 5.5 | 388 | 5.5 | 0 | 0.64 |
| ArceneCov | 10000 | 10 | 3600 | 10.8 | 377 | 10.8 | 0 | 9.55 |
| ArceneCov | 10000 | 15 | 3600 | 15.3 | 416 | 15.2 | 0.66 | 8.65 |

**Table 7: Summary of average approximation errors and average speedup factors for each $k$, when integrated our framework with Branch-and-Bound algorithm. The standard deviations are reported in the parentheses.**
| k | Avg. Error (Std. Dev) | Avg. Speedup (Std. Dev) |
| --- | --- | --- |
| 3 | 0.00 (0.00) | 6.10 (5.73) |
| 5 | 0.01 (0.03) | 69.97 (81.23) |
| 10 | 3.20 (5.11) | 186.94 (234.14) |
| 15 | 5.40 (12.05) | 112.07 (113.47) |

From the data presented in Table 6, it is evident that our framework achieves exceptional performance metrics.
The median approximation error is 0, with the 90th percentile approximation error not exceeding 3.85%.
Additionally, the median speedup factor is recorded at 39.
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 7, it is clear that our framework consistently obtains significant computational speedups across a range of settings. Notably, it achieves an average speedup factor of 6.1 for $k=3$, with this factor increasing dramatically to over 69 for $k\geq 5$.
An observation from the analysis is the relationship between $k$ and performance metrics: although the speedup factor increases with larger $k$, there is a corresponding rise in the average approximation errors.
Specifically, for $k=3$ and $k=5$, the average errors remain exceptionally low, under 0.01%, but they escalate to 3.2% for $k=10$ and 5.4% for $k=15$. Even after adjusting for outlier impacts from the LymphomaCov2 datasets, the trend persists, with average errors for $k=10$ and $k=15$ at 0.43% and 1.17%, respectively.
This indicates a consistent pattern of increased approximation errors as $k$ grows, which is aligned with Theorem 3, despite the significant gains in speed.

Although the overall performance of our framework looks promising, we observe that the performance of our framework in some instances presented in Table 6 is not ideal.
For example, for the instances $k=10$ and $k=15$ in the LymphomaCov2 dataset, our framework yields notably large approximation errors.
Additionally, we find that on ArceneCov dataset, specifically for $k=3$ and $k=5$, Algorithm 5 operates slower than the Branch-and-Bound algorithm.
In the following sections, we show that performance of our framework on these instances can be greatly improved by fine-tuning Algorithm 5.

#### C.2.1 Improving approximation errors

In the LymphomaCov2 dataset, instances with $k=10$ and $k=15$ exhibit notably large approximation errors when using Algorithm 5.
This issue primarily arises because the maximum block size constraint of $d_{0}=40$ is insufficient to include the support of the optimal solution for SPCA for these two instances, leading to sub-optimal solutions.
In Table 8, we summarize the numerical results when increasing $d_{0}$ to 50 and 60.
We observe a significant improvement on the approximation error, going from an average of 12.0% to less than 1%, at the cost of slightly lowering the speedups.

**Table 8: Comparison of errors and speedups for different $d_{0}$ values across LymphomaCov2 dataset.**
| Dataset | $d$ | $k$ | $d_{0}=40$ | $d_{0}=50$ | $d_{0}=60$ |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  |  | Error | Speedup | Error | Speedup | Error | Speedup |
| LymphomaCov2 | 4026 | 3 | 0.0 | 6.6 | 0.0 | 5.8 | 0.0 | 5.7 |
| LymphomaCov2 | 4026 | 5 | 0.0 | 4.5 | 0.0 | 4.4 | 0.0 | 4.2 |
| LymphomaCov2 | 4026 | 10 | 13.0 | 102.9 | 1.2 | 83.7 | 0.3 | 90.0 |
| LymphomaCov2 | 4026 | 15 | 35.0 | 52.2 | 3.0 | 43.9 | 2.4 | 19.5 |
| Overall |  |  | 12.0 | 41.5 | 1.0 | 34.4 | 0.7 | 29.8 |

#### C.2.2 Improving speedup factors

In ArceneCov dataset, Algorithm 5 exhibits relatively low speedup factors compared to other datasets.
This slower performance is due to the fact that Algorithm 5 necessitating numerous calls to Algorithm 3 compared to its solving process with other datasets, thereby significantly extending the runtime.
We managed to improve the performance of Algorithm 5 on the ArceneCov dataset by the following two simple improvements on Algorithm 5: (i) We set an extra stopping criteria on Algorithm 5, which breaks the while loop as soon as a problem with largest block size $d_{0}$ has been solved, and (ii) We set the tolerance parameter $\delta$ to $0.01\cdot\left\lVert A\right\rVert_{\infty}$ instead of $0.01$.
By making these two modifications, we manage to improve the average speedup factor from 4.77 to 12.9, and also improve the approximation from 2.55% to 0%. We summarize the results in Table 9.

**Table 9: Comparison of errors and speedups between Algorithm 5 and refined Algorithm 5.**
| Dataset | $d$ | $k$ | Algorithm 5 | Refined Algorithm 5 |  |  |
| --- | --- | --- | --- | --- | --- | --- |
|  |  |  | Error | Speedup | Error | Speedup |
| ArceneCov | 10000 | 3 | 0.0 | 0.2 | 0.0 | 0.6 |
| ArceneCov | 10000 | 5 | 0.0 | 0.6 | 0.0 | 1.8 |
| ArceneCov | 10000 | 10 | 0.0 | 9.6 | 0.0 | 27.7 |
| ArceneCov | 10000 | 15 | 0.7 | 8.7 | 0.0 | 21.4 |
| Overall |  |  | 0.2 | 4.8 | 0.0 | 12.9 |

### 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 4), in Table 10.
We also summarized the average approximation errors and average speedup factors for with respect to $k$ in Table 11.

In Table 10, we once again demonstrate the substantial speedup improvements our framework achieves when integrated with Chan’s algorithm.
The results show a median approximation error of only 0.38%, with the maximum approximation error not exceeding 2.73%, and a median speedup factor of 5.37.
We observe that the speedup factor grows (linearly) as the dimension of the dataset goes up, highlighting the efficacy of our framework, especially in large datasets.
Notably, in all instances within the DorotheaCov dataset, and the instance $k=1000$ in GLI85Cov dataset, our framework attains 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 2.

In Table 11, we do not observe an increase in speedup factors as
$k$ increases.
This phenomenon can be attributed to the runtime of Chan’s algorithm, which is $\mathcal{O}(d^{3})$ and independent of $k$.
Consequently, as predicted by Proposition 1, the speedup factors do not escalate with larger values of $k$.
Additionally, unlike the trends noted in Table 7, the average errors do not increase with $k$.
This deviation is likely due to Chan’s algorithm only providing an approximate solution with a certain multiplicative factor:
As shown in Theorem 3, our algorithm has the capability to improve this multiplicative factor at the expense of some additive error.
When both $k$ 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 are comparable to those obtained solely through Chan’s algorithm, maintaining an average error of less than 0.55% for each $k$.

Finally, we note that Chan’s algorithm is considerably more scalable than the Branch-and-Bound algorithm.
The statistics underscore our framework’s ability not only to enhance the speed of scalable algorithms like Chan’s, but also to maintain impressively low approximation errors in the process.

**Table 10: Comparison of runtime and objective values on real-world datasets, between Chan’s algorithm and Algorithm 5 integrated with Chan’s algorithm. $d$ is the dimension of Sparse PCA problem, and $k$ is the sparsity constant. We note that the objective value for GLI85Cov is scaled by $10^{9}$ and that for GLABRA180Cov is scaled by $10^{8}$. 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 | 186 | 91.3 | 0.44 | 1.23 |
| ArceneCov | 10000 | 500 | 239 | 140.4 | 187 | 137.5 | 2.07 | 1.28 |
| ArceneCov | 10000 | 1000 | 249 | 190.4 | 176 | 185.2 | 2.73 | 1.41 |
| ArceneCov | 10000 | 2000 | 256 | 233.3 | 199 | 232.9 | 0.17 | 1.29 |
| GLI85Cov | 22283 | 200 | 2836 | 35.8 | 847 | 34.9 | 2.51 | 3.35 |
| GLI85Cov | 22283 | 500 | 2728 | 38.8 | 923 | 38.8 | 0.00 | 2.96 |
| GLI85Cov | 22283 | 1000 | 2869 | 39.9 | 993 | 40.7 | -2.01 | 2.89 |
| GLI85Cov | 22283 | 2000 | 2680 | 42.3 | 1000 | 42.1 | 0.47 | 2.68 |
| GLABRA180Cov | 49151 | 200 | 35564 | 44.6 | 3796 | 44.1 | 1.12 | 9.37 |
| GLABRA180Cov | 49151 | 500 | 33225 | 60.4 | 4157 | 60.2 | 0.33 | 7.99 |
| GLABRA180Cov | 49151 | 1000 | 30864 | 70.5 | 4171 | 70.0 | 0.71 | 7.40 |
| GLABRA180Cov | 49151 | 2000 | 29675 | 78.1 | 3840 | 77.4 | 0.90 | 7.73 |
| DorotheaCov | 100000 | 200 | 246348 | 3.25 | 15864 | 3.31 | -1.85 | 15.53 |
| DorotheaCov | 100000 | 500 | 246539 | 4.94 | 17603 | 4.96 | -0.40 | 14.01 |
| DorotheaCov | 100000 | 1000 | 244196 | 6.23 | 16381 | 6.28 | -0.80 | 14.91 |
| DorotheaCov | 100000 | 2000 | 255281 | 7.18 | 17760 | 7.21 | -0.42 | 14.37 |

**Table 11: 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 parentheses.**
| k | Avg. Error (Std. Dev) | Avg. Speedup (Std. Dev) |
| --- | --- | --- |
| 200 | 0.55 (1.82) | 7.37 (6.44) |
| 500 | 0.50 (1.09) | 6.56 (5.73) |
| 1000 | 0.16 (2.04) | 6.65 (6.07) |
| 2000 | 0.28 (0.55) | 6.52 (5.92) |

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: Please address the concerns outlined in Weaknesses.
soundness: 3
strengths: 1. The framework significantly reduces runtime by decomposing the original problem into smaller, more manageable subproblems.

2. The paper provides approximation guarantees and time complexity analyses, ensuring that the method preserves solution quality.

3. Extensive experiments on diverse, large-scale datasets demonstrate the framework’s effectiveness in reducing computational time with minimal approximation errors.
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.
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.

2. 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.

[1] "Exact covariance thresholding into connected components for large-scale graphical lasso." Journal of Machine Learning Research, 13(27):781−794, 2012.

[2] "Graphical lasso and thresholding: Equivalence and closed-form solutions." Journal of Machine Learning Research, 20(10):1−44, 2019.

[3] "Learning large-scale MTP2 Gaussian graphical models via bridge-block decomposition." Advances in Neural Information Processing Systems 36:73211-73231, 2023.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
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? 
I was a little lost on the notation and didn't understand where the search over all permutations of rows and columns is addressed. 
Please clarify the computational complexity of this problem and how the algorithm gets around checking all row and column permutations. 


q2) 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. 

q3) Theorem 2 doesn't really use sparsity anywhere right? Its really a result about PCA for block-diagonal matrices, or do i misunderstand something?

Remark: 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: 
https://mathoverflow.net/questions/68041/showing-block-diagonal-structure-of-matrix-by-reordering
Another 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: 
Complexity Results for Bandwidth Minimization
M. R. Garey, R. L. Graham, D. S. Johnson, D. E. Knuth

q4) On experimental evaluation. 
Why 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  
the truncated power method (TPower) of (Yuan & Zhang, 2011). I believe Tpower is usually the fastest in practice. 
I 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. 

I'm willing to improve my score if the authors manage to convince me on these concerns.
soundness: 3
strengths: The paper has a cool observation and I think is well written. 
I have some concerns about the complexity of the approximation of a matrix by a block diagonal that I would really like clarified. 
The experimental section is reasonable some it misses additional experiments that I would like to see.
summary: This paper solves approximate Sparse PCA with matrix block-diagonalization. The proposed algorithm has 3 steps:
approximate a given matrix by a block diagonal, 
solve separate problems and then, combine them. 
 
The 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. 

The 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. 

The 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.

Edit after reading replies and discussion: I think the authors sort of addressed my concern, I give a small upgrade in my score
weaknesses: I have issues with the complexity of the block-diagonal matrix approximation and also with the experimental section, as I clarify below. 

Also some minor typos here:
Definition 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. 

typo: On the other hand, there also exist a number of algorithms that takes* polynomial runtime

typo: The axies* are indices of A; 

typo: th3, positive integer, and donote*

---

confidence: 2
contribution: 3
initial_rating: 8
presentation: 3
questions: See above.
soundness: 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.
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.
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?

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 2
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$)?
2. Could you explain why $|| E ||_\infty$ is the estimation in Model 1?
3. In algorithm 4, why $m=\lfloor (2C+2)d^{1+\alpha}  \rfloor$?
4. I'm curious about the inequality in Line 71. Can you prove it?
5. The proof of Theorem 1 is incomplete. Could you complete it?
6. In the proof of Theorem 3, how do you come up with $k || A-A^\epsilon ||_\infty \leq \epsilon $?
7. Could you prove why the g() function in Eq.(2) is convex?
8. Do you think the inequality in Line 41-45 is too relaxed?
9. Some descriptions are not rigorous. For example, "an absolute constant $c > 0$". Apparently, an absolute value is greater than 0 in this case.
10. 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.
11. 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?
12. Some sentences are too long for readers to understand. For example, Line 51-53. Could you rephrase them?
soundness: 3
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.  
2. The authors present theoretical analysis, providing guarantees for both approximation error and runtime complexity.  
3. This approach is evaluated across many real-world datasets, demonstrating consistent improvements in runtime and approximation error.
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.
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$).
2. There is not much discussion about the impact of hyperparameters.
3. In your notation, $p$ in $p$-norm should not represent the same value as $p$ in diag($A_1, \cdots, A_p$).
4. Some descriptions in your Theorem, Proposition, and Proof are not rigorous.
5. Some proofs are not completely correct.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 2
questions: I suggest that the authors take the above stated weaknesses into account to improve their manuscript. Apart from this suggestion:

Q1) How far from the bound implied by Theorem 1 are the measured suboptimality gaps in practice?

Q2) 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.

Q3) How hard is it to extend your results to the multi-component case?
soundness: 3
strengths: S1) The proposed scheme is simple and intuitively well-motivated.

S2) The theoretical results are also fairly simple and support the intuition.

S3) The analysis includes approximate algorithms, which have great practical relevance since the problem is computationally hard.

S4) 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.
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).
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.

W2) 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:
- Separate values for the different used $k$ should be reported (see point W2 below).
- 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. 
- 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).

W3) 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$. 

W4) 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).

W5) Though generally well-written, the paper lacks some clarity at times, and the notation is not always consistent/clear. In particular:
- 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.
- 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.
- 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.

Meta-review: This paper presents an efficient approach to Sparse Principal Component Analysis (SPCA) by leveraging a block-diagonal approximation of the covariance matrix, thereby reducing the original problem to a series of smaller subproblems. The paper is well-written, with the key ideas and contributions clearly presented. The authors provide  theoretical results, and the experimental section includes comparisons with other state-of-the-art methods on a variety of datasets, demonstrating advantages in both runtime efficiency and average approximation error. The reviewers acknowledged the contributions of this work, and the authors sufficiently addressed the minor concerns raised during the rebuttal phase. Overall, this paper makes a clear contribution to scalability of sparse PCA, which is a relevant topic of research.  Therefore, I recommend its acceptance to ICLR.

Decision: Accept (Poster)