Invariant subspaces and PCA in nearly matrix multiplication time

We thank the reviewer for their useful feedback and for their suggestions to improve the readability of our manuscript. Below we provide replies to the specific questions that were raised:

1. Question: The assumption of $H$ being symmetric needs to be in the Abstract as well, as it is quite significant.

   Answer: We will fix the missing mention in the abstract that $H$ must be Hermitian and $S$ Hermitian and positive-definite. Indeed, these are crucial assumptions that must be emphasized.

2. Question: Move the main Algorithms to the main body as it gives a more visual impact, which is of higher importance for a conference paper.

   Answer: We will move the most relevant algorithms from the appendices into the main body of the manuscript, taking advantage of the additional page of content that is allowed in the final version.

3. Question: It would be helpful if the authors could show numerical results w.r.t. accuracy, complexity, or relevant quantities to support their theoretical analysis.

   Answer: Due to the limited time, we well not be able to incorporate numerical results in our manuscript, but will provide indicative examples in the final version.

4. Question: Missing citation to kernel PCA [1] Schölkopf, Bernhard, Alexander Smola, and Klaus-Robert Müller. "Nonlinear component analysis as a kernel eigenvalue problem." Neural computation 10.5 (1998): 1299-1319.

   Answer: We are aware of the paper by Schölkopf et al. and cited it twice in our manuscript (Ref. [116]). If other relevant references are missing, we would gladly add them.

We hope that our responses clarify the concerns raised by the reviewer. If (other) major concerns remain, we will address them to the best of our capabilities. Again, we would like to thank the reviewer for their useful recommendations.

We would like to thank the reviewer for their constructive feedback and for taking the time to review our manuscript. The reviewer is absolutely right about the rectangular FMM. Actually, our comment about rectangular FMM only applies to the analysis of the Block-Krylov method in Section 4.2 and in Appendix G, where we used blocked-square FMM to perform Square$\times$Rectangular matrix multiplication. Using Rectangular FMM instead of blocked-square FMM might provide minor (but not substantial) improvements for specific regimes of the input parameters.

If more information is needed, we would be pleased to further expand on the rectangular FMM or other topics.

We would like to thank the reviewer for their comments and carefully examining our manuscript as well as Appendix B. Our replies to their comments can be found below:

1. Question: Could the authors please elaborate on the difference between their result and the forward error approximation in Theorem 5.2 of [1] (which is the [37] citation of the authors)? The authors of [1] write that their Theorem 5.2 is restrictive as it needs $S$ to be invertible and $H$ to have a simple eigenspectrum (no repeated eigenvalues). It seems that Theorem 1.1 by the authors also requires $S$ to be invertible but can allow a non-simple spectrum.

   Answer: This is indeed a very interesting point that deserves further investigations. We have been in close contact with the authors of [37] the past months. After discussing and comparing our results in detail, we concluded that our works complement each other. The major difference comes from the fact that our Theorem 1.1 is stated specifically for Hermitian-definite pencils: $H$ must be Hermitian, while $S$ must be Hermitian and positive-definite, which is stricter than just invertible. However, we do not need simplicity in the spectrum because of Weyl's inequality (or the slightly more general Kahan's bound): For Hermitian matrices, a backward-error directly translates into a forward-error in the eigenvalues (see also our Lemma B.1 in Appendix B). This property does not hold for general matrices, as the ones that are assumed in Theorem 5.2 of [37]. In that case, one would resort to either the Bauer-Fike theorem, or to the properties of the pseudospectrum, in order to obtain forward-error bounds from the backward-error. If the reviewer is interested, we can also refer to the thesis of the last author of [37] for additional remarks.

2. Question: Is a general dependence on the term gap_k necessary for analysis like this?

   Answer: We have thoroughly discussed the $gap_k$ dependence in the bounds with other authors and experts in the field. It appears to be unavoidable. Intuitively, if we want to distinguish between two invariant subspaces, the error must be smaller than the gap. Consider the following arbitrary $2\times 2$ example: $A = [ 1+\epsilon, 0; 0, 1 ].$ The gap is equal to $\epsilon$, and the spectral projector on the top-1 invariant subspace is $[ 1, 0; 0, 0].$ Now consider the following $\epsilon$-perturbation of $A$: $\widetilde A = [1, 0; 0, 1+\epsilon]$. The spectral projector of the top-1 invariant subspace is completely different than before. In fact, it is orthogonal to the first one.

We hope that these responses clarify the reviewer's concerns. Let us know if the answers are not sufficient so we can elaborate further.

We would like to thank the reviewer for their positive feedback and for carefully reviewing our manuscript! Below we provide answers to the specific questions that were asked:

1. Question: Parameter Selection: How should the hyperparameters (e.g., spectral gap, condition number) be chosen in practice? Are there heuristic methods or guidelines for selecting these parameters to ensure optimal performance?

   Answer: We would like to stress out that all required parameters (gaps, condition numbers, norms, etc) are determined by the algorithm, i.e., the user does not need to provide them. The only inputs are the matrices to be handled, their size $n$, and the desired accuracy $\epsilon\in(0,1)$. A reduced set of input parameters is an additional contribution compared to previous algorithms, where methods to determine the required parameters are typically not provided. In our case, the complexity of determining the parameters is included in the main theorems.

2. Question: Extensions to Sparse Matrices: How does the proposed method handle sparse matrices common in large-scale applications? Are there any modifications or optimizations specific to sparse matrix structures?

   Answer: Extending these methods to sparse matrices is a very intriguing, but rather difficult open problem. To-date, for large-scale sparse problems, a Krylov-based method is probably the algorithm of choice. It takes advantage of matrix-vector products and exhibits a somewhat sparsity-time complexity. We thoroughly described the bit-complexity analysis of the celebrated Block-Krylov PCA method (see Section 4 and Ref. [100]). The problem with such methods is that $k$ (the target rank) must be very small to work well. As mentioned in Section 5 (iii), large-scale banded problems can probably be solved with so-called ``superfast eigensolvers'' (e.g., refs [61, 131]), but we are not aware of any end-to-end complexity/stability analysis of these algorithms. For general sparse matrices (non-banded) and large $k$ we believe that $O(n^\omega)$ is probably the best that can be achieved.

• Question: Complexity of Implementation.

  Answer: We expect that the implementation of the proposed methods does not involve substantial difficulties (although it remains to be demonstrated). The main steps are the diagonal perturbations, the Newton iteration, and the bisection-like algorithm to compute the spectral gap. The level of difficulty should definitely not exceed that of, for example, Ref. [11,37] (for which open source implementations exist).

To conclude, we would like to thank again the reviewer for the detailed feedback and for their interest, they are both highly appreciated! Let us know if further clarifications of our replies are required.

We would like to thank the reviewer for their positive feedback. Below we elaborate on the specific questions that were raised:

1. Question: In the statement of Theorem 1.1 and Proposition 2.1, is $||H||\leq 1$ a proper assumption?

   Answer: $||\cdot||\leq 1$ is typically a standard assumption in eigenvalue algorithms. For Theorem 1.1., in particular, $||H||$ can be approximated up to a constant factor in floating point using Lanczos, e.g., by applying Theorem 18 in the ArXiv version of Ref. [101]. We can then scale $H$ by the approximated norm. Similarly with $S^{-1}$, we can use the $\mathsf{SIGMAK}$ algorithm in Appendix E to compute the minimum singular value of $S$, and scale accordingly.

2. Question: Also in Theorem 1.1, what if $gap_k$ is very small for the given $k$ (say, $k=2^{-n}$)?

   Answer: This is an excellent question that we discussed with other authors and experts in the field. A dependence on the gap should be unavoidable. Otherwise, it is not possible to distinguish between the two subspaces and obtain forward-error bounds (backward-error bounds are easier). If $gap=2^{-n}$, then the problem is rather ill-posed to begin with, but the algorithm will still terminate in polynomial time.

3. Question: How does the result of Theorem 3.1 compare with previous works (if any)? How hard is it to approximate $\mu_k$ and $gap_k$ for a given $k$, is there any previous work that also achieve the $\tilde O(n^\omega)$ bit complexity result?

   Answer: For general $k$ we are not aware of any other $\widetilde O(n^\omega)$ algorithms, except for Banks et al. (Ref. [11]). We briefly described in Proposition B.2 how to use [11] to compute the gap in $\widetilde O(n^\omega)$. However, [11] leads to a slightly slower algorithm than our main $\mathsf{GAP}$ algorithm, which is diagonalization-free. We are also not aware of any lower bounds to compute a single spectral gap. We believe that it should be $\Omega(n^\omega)$-hard, but we have no formal proof for that.

4. Question: Can the authors provide a more detailed discussion of how the proposed algorithm can (or cannot) be applied to real world problems such as solving large scale PCA tasks?

   Answer: When it comes to practical applications, the main benefit of the proposed ``matrix multiplication reduction'' is not the asymptotic complexity of the resulting algorithm, but rather the fact that we can exploit matrix multiplication in the computations. It is fully parallelizable and lends itself to cache-optimization, which can provide significant speedups (up to orders of magnitude) in practice. We also think that random perturbations (diagonal, Ginibre, ...) might actually become very useful in the near future, due to their eigenvalue repulsion properties, but this field is currently under-explored.

Thank you again for the feedback and for the interest! We are happy to discuss these points more in depth, if needed.