Efficient Sparse PCA via Block-Diagonalization

Dear Reviewer,

Thank you for taking the time to review our paper and provide thoughtful comments. We have outlined our responses and clarifications below. We hope that our efforts have addressed your concerns, and we kindly request that you consider increasing your scores.

Responses to Weaknesses:

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

Response: We thank you for pointing this out. In fact, our framework requires the input matrix $A$ to be symmetric. We have made the clarification on Page 3 in our revised manuscript. Other than that, our framework accommodates general symmetric input matrices, provided the subroutine used is compatible with such matrices. For PSD input matrices, our framework can be modified to achieve better approximation error. Specifically, if we define $d^\star = intdim(A, \epsilon)$, the additive error term in Theorem 3 decreases from $a\left(k, d^\star\right) + (1 + \frac{1}{m\left(k, d^\star\right)}) \cdot k\epsilon$ to $a\left(k, d^\star\right) + \frac{2k\epsilon}{m\left(k, d^\star\right)},$ using a similar proof to the one presented for Theorem 3. The modification is on Algorithm 2: instead of outputting the thresholded matrix blocks, we output the corresponding matrix blocks from the original input matrix $A$, which ensures the blocks are PSD. We are happy to add more discussions on this during the revision.

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

1. Separate values for the different $k$ used should be reported (see point below).
2. It would be interesting to report the threshold value $\epsilon$ 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.
3. 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).

Responses: We appreciate your suggestions. Regarding the values for separate $k$, we have reported and discussed the average speedup factors and average approximation errors in Appendices C.2 and C.3 of our original manuscript. The summary of results across different $k$ values was not included in the main text due to space constraints. We next report the relative threshold value $\epsilon / |A|$, the percentage of zero entries after the thresholding, and Jaccard index for selected datasets. We first present the results on the LymphomaCov1 dataset, comparing our Algorithm 5 with Branch-and-Bound algorithm and the vanilla Branch-and-Bound algorithm:

From the table above, it is evident that the change of Jaccard index aligns with the change of the approximation error. Additionally, the relative threshold value $\epsilon / |A|$ and the percentage of zeros remain relatively stable across $k$, as the largest block size $d_0$ is set to 40, and the best solutions are typically found in blocks near this size. We next present results on the GLABRA180Cov dataset, comparing our Algorithm 5 with Chan’s algortihm and the vanilla Chan’s algorithm:

For Chan’s algorithm, the Jaccard index does not vary significantly with $k$, nor does its change align with the change of the approximation error. A possible reason is that Chan’s algorithm is an approximation algorithm, and the Jaccard index may not fully capture the similarity between the solutions obtained and the true optimal solutions. As expected, since $d_0 = 2k$, the relative threshold value $\epsilon / |A|$ and the percentage of zeros decrease as $k$ increases. We are happy to conduct more experiments and add discussions on this during the revision.

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

Response: We have made clarifications in Theorem4, as well as included the additional details in the proof of Theorem 4 in our revised manuscript. The high level idea of obtaining the running time is that, the running time for each iteration is upper bounded by $\mathcal{O}\left( \left\lceil \frac{d}{d_0} \right\rceil \cdot g\left(k,d_0\right) + d^2\right)$. This is clear due to the fact that $g$ is convex and the trivial fact that $g(k, 0) = 0$, and the fact that the call of an algorithm $\mathcal{A}$ would be executed only if the intrinsic dimension is less than or equal to $d_0$. Since the total number of iterations is upper bounded by $\mathcal{O}(\log(|A| / \delta))$, we can obtain the desired running time. For the approximation error, we directly make use of Theorem 3.

10. I noticed that you frequently use expressions like $\max\max\max$ or $\inf\inf\inf$ in your proof, which appears unprofessional. Could you consider rephrasing these using constraints for clarity?

Response: Our intention is to remain consistent with the expressions used in Comminges et al. (2021). This is a common mathematical notation in the literature, as each $\max$ or $\inf$ corresponds to taking the maximum or infimum with respect to different variables, indices, or distributions.

11. Some sentences are too long for readers to understand. For example, Line 51-53. Could you rephrase them?

Response: We are happy to rephrase the sentence in our revised manuscript. The intention is to introduce that there are three types of algorithms: Some algorithms are fast but yield sub-optimal solutions, others provide high-quality solutions at the cost of greater computational complexity, and a few achieve both efficiency and accuracy but only under specific statistical assumptions. We have modified the sentence in our revised manuscript, as highlighted on Page 1.

Dear Reviewer,

We are grateful for your detailed feedback and helpful suggestions. Please find our responses and clarifications below.

Responses to 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$).

Response: We thank the reviewer for highlighting this concern. However, in our paper, it is not necessary to know $p$ or the specific values of $d_i$ prior to running the algorithms, as these quantities become clear during execution. More specifically, $p$ simply represents the number of blocks obtained after running Algorithm 2, and $d_i$ denotes the dimensions of these blocks. Both $p$ and $d_i$ are determined automatically and become evident upon completion of Algorithm 2.

2. There is not much discussion about the impact of hyperparameters.

Response: Except for Algorithm 5, our framework does not include any hyperparameters (we assume that the parameters are known beforehand in Algorithm 4). We have discussed the hyperparameters $d_0$ and $\delta$ in Algorithm 5 in Appendix C, on Page 22 (Page 24 in our revised manuscript). In summary, increasing $d_0$ yields a better solution but also increases the runtime of our framework. Setting $\delta$ to a moderate value that scales with $|A|_{\infty}$ has minimal impact on solution quality and can speed up the framework compared to using a very small universal constant for $\delta$.

3. In your notation, $p$ in $p$-norm should not represent the same value as in $diag(A_1, A_2, \dots, A_p)$. Some descriptions in your Theorem, Proposition, and Proof are not rigorous. Some proofs are not completely correct.

Response: In fact, $p$ in $p$-norm is used only for vectors, while in $diag(A_1, A_2, \dots, A_p)$ the integer $p$ denotes the total number of blocks in matrices. To avoid further confusion, we have changed $p$-norm to $q$-norm on Page 4 in our revised manuscript.

We believe that we have now provided rigorous statements and proofs in our revised manuscript, addressing your comments below. Additionally, we are happy to clarify any points that may still be unclear to you.

Responses to questions:

1. Could you explain why $|E|_{\infty}$ is the estimation in Model 1?

Response: The intuition for estimating $|E|_{\infty}$ is that we want to obtain an estimate of $A$ in Model 1, and then by our denoising algorithm Algorithm 1, we can get rid of the impact of $E$ and find a high quality approximate solution to SPCA with input $(A, k)$.

2. In algorithm 4, why $m = \lceil (2C+2) d^{1+\alpha} \rceil$?

Response: We choose this particular value as we make use of the analysis proposed in Comminges et al. (2021), where the idea is to make sure there are sufficient random variables in each block to make accurate estimation. The details are left in Proof of Proposition 3 on Page 17 - 18.

3. I'm curious about the inequality in Line 71. Can you prove it?

Response: We cannot find any inequality on Line 71, as it is Figure 1 in our paper introducing our proposed framework. However, we are happy to prove it if you could clarify which inequality you are referring to.

4. The proof of Theorem 1 is incomplete. Could you complete it?

Response: We respectfully disagree with this assessment. We believe that we have proved the desired result stated in Theorem 1. However, we are happy to provide further clarification if you could indicate which part of the proof is not clear.

5. In the proof of Theorem 3, how do you come up with $k|A - A^\epsilon| \le \epsilon$?

Response: We thank the reviewer for spotting this typo. What we meant is $\epsilon \cdot k$ on the right-hand-side, and we have corrected it on Page 16 in our revised manuscript. Note that we indeed use the correct inequality in the last inequality in this proof.

6. Could you prove why the g() function in Eq.(2) is convex?

Response: We assume in Proposition 1 that the function $g(k,d)$ is convex with respect to $d$.

7. Do you think the inequality in Line 41-45 is too relaxed?

Response: We cannot find any inequality on lines 41–45, as this part is focused on the literature review. However, we are happy to provide further clarification if the reviewer could specify which inequality they are referring to.

8. Some descriptions are not rigorous. For example, "an absolute constant $c>0$". Apparently, an absolute value is greater than 0 in this case.

Response: What we intended to emphasize is that such a constant is a universal constant, independent of other factors, such as the dimension of the problem. It does not refer to absolute values.

Thank you for your response. Most of my questions are answered.

I’m curious about the inequality in Line 771. Can you provide a proof for it? (Question 3)

For the proof of Theorem 1, I think it would be clearer if the author could state $||A-A^\epsilon||_\infty \le 2\epsilon$.

Dear Reviewer,

We are grateful for your detailed feedback. Please find our responses and clarifications below.

Responses to Weaknesses:

1. 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?

Response: Thank you for raising this insightful question. We believe that the issue you mentioned does not affect any of our theorems. In essence, Theorems 1 and 2 do not rely on finding an optimal solution; rather, they are universally true properties of the optimal solutions to SPCA, independent of the specific method used to solve SPCA. For Theorems 3 and 4, we require only that an approximate algorithm achieves a certain approximation ratio. As long as the algorithm satisfies the stated approximation bound for the input, the theorems hold. For instance, Chan’s algorithm can provide such an approximate solution for any input matrix $A$, including the scenario you described with nearly identical leading eigenvalues.

Dear reviewer,

Thank you very much for your comments and suggestions. Please find below our responses and clarifications. We believe we have addressed all your concerns and would appreciate it if you could consider raising your scores.

Responses to Weaknesses:

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

Response: We thank you for raising this question. While finding the best block-diagonal approximation to a matrix $A$ is a hard problem, we don’t need to find the best one. Instead, we only need to find a matrix approximation that is $\epsilon$ away (after permutation) under L_inf distance, for some moderate $\epsilon > 0$. It can be identified in $O(d^2)$ using Algorithm 1 and Algorithm 2. Although such approximation might be far away from the best block-diagonal approximation, applying Algorithm 3 to this approximation matrix yields desired computational speedups and approximation guarantees.

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

Response: We thank you for highlighting this insightful high-level intuition behind the result. To the best of our knowledge, this specific result has not been previously developed in the literature. We conduct a rigorous proof by contradiction on pages 14–15 (and on Page 16 in our revised manuscript) in the Appendix to prove it formally. We also believe the implications of this result are non-trivial, as it plays a key role in the development of our framework, as outlined in Remark 1 on Page 6.

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

Response: In Theorem 2, we develop a more general result that applies to both PCA and SPCA, depending on the values of $k$ and $d$. In the special case where $k = d$, we show that there exists an optimal solution to PCA that lies within one of the block sub-matrices. For the case where $k < d$, we demonstrate that finding a $k$-sparse optimal solution to SPCA is equivalent to solving SPCA within each block.

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

Response: We thank you for the feedback. Our framework is mainly designed to speed up SPCA algorithms that come with approximation guarantees. TPower, while efficient in practice, lacks such guarantees and may not converge to global optima. This was the primary reason for not including it in our initial comparisons.

To address your concerns, we conducted additional experiments with FullPath and TPower. Specifically, we run three sets of experiments on the DorotheaCov dataset: (1) vanilla TPower, (2) Algorithm 5 with FullPath, and (3) Algorithm 5 with Chan’s algorithm. Table below shows the objective values identified by each approach for different values of $k$. As you can see, TPower only identifies sub-optimal solutions with objective values ~ 0.2, yet both Alg 5 w/ FullPath and Alg 5 w/ Chan’s obtain solutions with much higher objective values. We didn’t include vanilla FullPath and Chan’s algorithms in this table since they have not outputted a solution after 24 hours. In contrast, Alg 5 w/ FullPath and Alg 5 w/ Chan’s output solutions in 4 - 5 hours.

Dear Reviewer,

We greatly appreciate your valuable feedback and constructive suggestions. Below, we have provided our responses and clarifications. We believe we have addressed your concerns and would be grateful if you could consider increasing your scores.

Responses to 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.

Response: Thank you for your suggestion. We have added additional discussion on this point and position our contribution to the broader context of optimization techniques, in Appendex A on Page 14 in our revised manuscript.

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.

Response: We appreciate it for raising this important point. Our framework (Algorithm 3) actually ensures that the final solution remains $k$-sparse, as long as it calls a valid SPCA algorithm $\mathcal{A}$. Indeed, the final output is guaranteed to be $k$-sparse thanks to the construction on lines 5--8 of Algorithm 3, provided each $y_i$ is $k$-sparse. Examples of such SPCA algorithms $\mathcal{A}$ include Chan's algorithm or Branch-and-Bound algorithm, both of which are discussed in detail in our paper.