Invariant subspaces and PCA in nearly matrix multiplication time
Authors: Aleksandros Sobczyk
Marko Mladenovic
Mathieu Luisier
Venue: NeurIPS 2024 poster (2024)
Keywords: Invariant subspace
Generalized eigenvalue problem
PCA
Spectral projector
Spectral gap
Matrix multiplication
Bit complexity
Primary area: learning_theory
TL;DR: Invariant subspaces and PCA embeddings can be provably approximated in nearly matrix multiplication time in finite precision

Abstract: Approximating invariant subspaces of generalized eigenvalue problems (GEPs) is a fundamental computational problem at the core of machine learning and scientific computing. It is, for example, the root of Principal Component Analysis (PCA) for dimensionality reduction, data visualization, and noise filtering, and of Density Functional Theory (DFT), arguably the most popular method to calculate the electronic structure of materials. 
Given Hermitian $H,S\in\mathbb{C}^{n\times n}$, where $S$ is positive-definite,  let $\Pi_k$ be the true spectral projector on the invariant subspace that is associated with the $k$ smallest (or largest) eigenvalues of the GEP $HC=SC\Lambda$, for some $k\in[n]$. 
We show that we can compute a matrix $\widetilde\Pi_k$ such that $\lVert\Pi_k-\widetilde\Pi_k\rVert_2\leq \epsilon$, in $O\left( n^{\omega+\eta}\mathrm{polylog}(n,\epsilon^{-1},\kappa(S),\mathrm{gap}_k^{-1}) \right)$ bit operations in the floating point model, for some $\epsilon\in(0,1)$, with probability $1-1/n$. Here, $\eta>0$ is arbitrarily small, $\omega\lesssim 2.372$ is the matrix multiplication exponent, $\kappa(S)=\lVert S\rVert_2\lVert S^{-1}\rVert_2$, and $\mathrm{gap}_k$ is the gap between eigenvalues $k$ and $k+1$. 
To achieve such provable "forward-error" guarantees, our methods rely on a new $O(n^{\omega+\eta})$ stability analysis for the Cholesky factorization, and a smoothed analysis for computing spectral gaps, which can be of independent interest.
Ultimately, we obtain new matrix multiplication-type bit complexity upper bounds for PCA problems, including classical PCA and (randomized) low-rank approximation.

Invariant subspaces and PCA in nearly matrix multiplication time
Aleksandros Sobezyk IBM Research and ETH Zurich obc0zurich ibm. com
Marko Mladenović
Mathieu Luisier
ETH Zurich mmladeeovvcc@eth..coo
ETH Zurich
mluisier@iis.cm www..eeth....
Abstract
Approximating invariant subspaces generalized eigenvalue problems (GEPs) is fundamental computational problem the core of machine learning and sci- entific computing. is, for example, the root Principal Component Analysis (PCA) for dimensionality reduction, data visualization, and noise filtering, and of Density Functional Theory (DFT), arguably the most popular method to calcu- late the electronic structure of materials. Given Hermitian H, S C**n, where S is positive-definite let Ilg be the true spectral projector on the invariant sub- space that is associated with the & smallest (or largest) eigenvalues of the GEP HC SCA, for some k [n]. We show that we can compute a matrix II, such that | &l S c, in O (n=** polylog(n, e c_ k()), g)pp bit oper- ations in the floating point model, for some e E (0, 1), with probability 1 - 1/n. Here, 17 > 0 is arbitrarily small, S 2.372 is the matrix multiplication expo- nent, x(S) = ||S || and gap; is the gap between eigenvalues k and k ++1 1. achieve such provable "forward-error" guarantees, our methods rely on a new O(n*+)) stability analysis for the Cholesky factorization, and smoothed analy- sis for computing spectral gaps, which can be independent interest. Ultimately, we obtain matrix multiplication-tyi bit complexity upper bounds for PCA problems, including classical PCA and (randomized) low-rank approximation.
Introduction
Generalized eigenvalue problems (GEPs) arise naturally a plethora of applications in machine learning, scientific computing, and engineering. Given a pair of matrices H and S, often referred to as matrix pencil, the problem of interest has the following form
HC = SCA,
(I)
where C and A are the unknown eigenvector and eigenvalue matrices, respectively. Of particular importance are the so-called "Hermitian definite" simply "definite" GEPs/pencils, in which case H is Hermitian and is Hermitian and positive-definite In many important applications, the quan- tity of interest is an (arbitrarily large) subset of eigenvectors, defining invariant subspace, rather than the entire C and A solutions of the GEP.
In data science and machine learning, invariant subspaces play central role in many problems, including Spectral Clustering [111, 124], Language Models [73], Image Processing [118, 11], Rec- ommendation Systems [46], Principal Components Analysis (PCA) [45, 79, 122], Support Vector Machines [99], and many others [50, 35, 96, 12]. We particularly focus on PCA applications, which take the form of GEP as in Eq. (1) where H is the sample covariance and S the identity. In more advanced settings, H and can be defined over kernel [122, 20]: See Section and Appendix G for more details. Another closely related application comes from Density Functional
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Theory [88] (DFT), which is not a machine learning problem per se, but it is probably the most com- monly used method (it was awarded the Nobel prize Chemistry in 1998) compute the electronic and structural properties of materials. In this case, H is the Hamiltonian and S the overlap matrix the smallest
(cf. Appendix F). The spectral projector on the invariant subspace corresponding generalized eigenvalues (occupied energies) directly provides the density matrix and the electron density. Obtaining them often presents a challenge from the computational point of view.
1.1 Problem definition
The main focus of this work is the computation of spectral projectors on invariant subspaces that are associated with a subset of the spectrum of Hermitian definite GEPs. As the Abel-Ruffini theorem excludes exact computation, even in exact arithmetic, we seek for approximate computations, as described in the following Problem 1.1.
Problem 1.1 (Spectral projector). Given Hermitian definite GEP HC SCA of size n, an integer S k S n 1, and accuracy e E (0, 1), compute a matrix I1g E C such that
f
Ix <.
(2)
where II, is the true spectral projector on the invariant subspace associated with the smallest or largest eigenvalues.
Before proposing algorithms to solve Problem 1.1, we first make some clarifications and define useful concepts.
Type of approximation:
The approximation of the form of Equation (2) is commonly called a "forward approximation" or "forward error" in numerical analysis. It quantifies the distance between the true solution of the problem and the one returned by an algorithm. It is a stronger and harder to achieve notion of approximation than the related "backward error." For details see Appendix A.1.
Model of computation: While many finite precision models of computation exist in the litera- ture, all algorithms in this work are analyzed in the floating point model of computation, which is also the prominent model implemented in existing computers. Each real number a is rounded to floating point number fl(()) = (1 + 0)a, where 0 E C satisfies [0] u Rz0- The machine preci- sion bounds also the errors introduced by arithmetic operations {+, -, /)/ and the expression log(1 /)) gives the number of bits required to achieve the desired precision. More details can be found in Appendix A.2.
The complexity of numerical algorithms is often measured in terms the arith- metic operations executed, commonly referred to as arithmetic complexity. A more realistic notion the bit complexity, which bounds the number of boolean operations. In the floating point model, straightforward to translate the arithmetic to the bit complexity if we have an upper bound on the number of bits. For instance, arithmetic operations on b bits can be typically carried out O66)) bit operations, or even faster by using more advanced algorithms [123, 55, 68].
Bit complexity:
Matrix multiplication time: two seminal works [39, 40], it was demonstrated that matrix mul- tiplication and other fundamental problems in Numerical Linear Algebra can be solved in the float- ing point model with nearly (n**t)) bit-complexity (up to polylogarithmic factors), where is an arbitrarily small positive number and is the matrix multiplication exponent, to-date bounded by w S 2.372 [47, 141, 1, 6). Hereafter, we will use the notation (n) = =n**.
1.2 Existing algorithms
Here give brief overview of existing algorithms. We refer Appendix A.7 for more details. GEPs in general can be solved using classic eigensolvers and related techniques in Ö(n²) floating operations, e.g., by reducing the matrix (or pencil) to tridiagonal form with similarity transforma- tions and applying the shifted QR algorithm on the tridiagonal matrix, or by using a divide-and- conquer method (see [37, 109, 75, 8, 9, 64, 42, 10, 113] and references therein). Significant pro- gresses beyond the (n²²) bit complexity barrier have been made [40, 39, 15, 95, 25, 41, 107, 121].

---

Regarding the computation of eigenvalues, two notable examples are the (T)M((n)) algorithm of [95] for the largest eigenvalye, and the O(n²) algorithm of [107] for the spectral norm.
The first have addressed the problem of computing invariant subspaces in nearly O(TMM (n)) in floating point is [39] (see also [14, 13J). The authors described an iterative algorithm for the Schur decomposition, and showed that each individual step is numerically stable, and it takes ((TMa(()) operations. An end-to-end bound on the number of iterations to achieve a backward approximate solution was left open. More recently, the seminal work of [15] extended the analysis to obtain an end-to-end (T)M complexity to approximately diagonalize a matrix, and [41, 121] provided rigorous analysis for the generalized eigenproblem case. In Corollary 1.7 and Proposition 1.1 of [15], it was also outlined how to translate the backward diagonalization error to a forward error for the eigenvectors. The reported bound, however, has two main limitations: it relies on simplicity of the spectrum, which is a strict assumption, and it requires as input an over-estimate on the eigenvec- tor condition number of the problem. which is unknown, and [15] does not describe how to compute it (see Appendix B.1 for more details). In this work we describe how to overcome these limitations and provide novel, end-to-end, provably accurate analysis (in the sense of Eq. (2)) for arbitrary invariant subspaces of definite GEPs with ((TMm(n)) boolean complexity.
1.3 Contributions and methods
Our main contribution, summarized in the following Theorem 1.1 and Algorithm 1, an end-to-end analysis to solve Problem 1.1 in nearly ((TM(())) time.
Theorem 1.1. Let (H, be a Hermitian definite pencil of size n with |H] IS ||- s 1, A s 22 < < An its eigenvalues, gap; = X X and k(S) = - |||||||| Algorithm
I PROJECTOR(H, S,k,
takes as imputs H, S. an integer & [n - 1}, an error parameter E (0, and returns matrix ITg such that
Pr [. - H. - sc
where Il is the true spectral projector on the invariant subspace that is associated with the k smallest (or largest) eigenvalues. The algorithm executes
gaps O (Thna(n) (log( [ي ( + log((x((S)) (og(n(S)) log (log(( ksS) egaps
floating point operations with O (log(n) (1og't (log'( kppp) ) log' (nk(()) + log
(ERAny gaDy log( gAps
)))
bits of precision. Internally, the algorithm needs to generate a total of at most Ö(n) standard normal floating point numbers using additional O(log(og(())) bits.
To achieve the results of Theorem 1.1, we provide a novel (n)Mm())) complexity analysis of several problems in numerical linear algebra that can be of independent interest.
In brief, our methodology is as follows. We first observe that we can determine reasonable "guesses" for the spectral gap (gap;) and for the midpoint (7²) between the X₃ and λ4+1 eigen- values then we can efficiently compute the spectral projector by approximating the sign function
sgn(āx - S S-H
using the analysis of [15] for the Newton iteration. The matrix ([ + sgn(āu: SH) indeed transforms in exact arithmetic all eigenvalues that are smaller than 7 to and the ones larger than j&g to zero. As will be proved in Proposition 2.1, in Section 2, this approach is sufficient to provide an accurate spectral projector II in floating point. As consequence, the problem reduces to approximating the aforementioned midpoint and gap. As a baseline, in Appendix B. we prove that this can be done in nearly similarly, by iteratively calling generalized diagonalization [41]. However, this approach presents two drawbacks: does not take advantage of the inherent symmetry of the problem, and, at the same time, it performs full diagonalization when only interested in the gap between two
with iterative inversion [39] and diagonalization [15] or,
3

---

specific eigenvalues, which is seemingly redundant. We formally prove this claim by designing a novel approach that achieves better complexity, typically by a factor of (I)g(n)) (cf. Section 3.4). Importantly, no explicit diagonalization is necessary.
To minimize the complexity of our algorithm, it is crucial leverage symmetry. To that end the Cholesky factorization of S in the spirit of the Cholesky-QR algorithm [37]. We highlight that, while other factorizations have been solved in ((TMM(n)) in floating point, an end-to-end analysis for Cholesky remains open. In exact arithmetic, for example, the LU of a Hermitian definite matrix directly provides its Cholesky and [39] showed that the LU factorization of non-symmetric matri- ces can be obtained in O (n))) However, when considering arithmetic errors, the relationship between LU and Cholesky does not hold in floating point, demonstrated by the counter-example of Appendix C.5. Other fast Cholesky algorithms have been proposed for special classes of ma- trices, e.g., for matrices with well-defined separators [58, 94, 61] and graph Laplacians [91, 92]. However, they do not generalize to arbitrary dense matrices. Our analysis is the first to improve the classic (nn²) floating point Cholesky algorithms [84, 71] for the general case, with provable error bounds. In the following Theorem 1.2 we summarize our new analysis of the Cholesky factorization Algorithm (see also Appendix C). We note that the algorithm itself is not new, only its analysis.
Theorem 1.2. Given a Hermitian positive-definite matrix M. there exists algorithm L &- (MOLESKY(M)), listed in Algorithum 2, which requires ([TMM(n)) arithmetic operations. This algorithm is logarithmically stable, in a sense that, there exist global constants C2 Cg2 Cg2 such that for all c E (0, 1), if executed in a floating point machine with precision
u S UCCOLLSKY == "1 (M) log
which translates into O (log(n) log(k(M)) - log(t)) required bits of precision, then it does not
break down due to arithmetic errors, and the solution returned satisfies |LLL" M|| S <||M]].
This stand-alone result fulfills the definition of "logarithmic-stabbliyy"" a notion of numerical stability that is commonly used in the related literature [39, 15]. Given this new Cholesky analysis, the following transformation of the GEP to regular Hermitian eigenvalue problem:
HC = SCA = L LL'LLLL = (L-IC)A,
can be carried out accurately in (T)) in floating point, with provable forward-error bounds for all eigenvalues of the transformed problem. Here, L is the Cholesky factor of s instead of S. Specifically, in Proposition C.3 in Appendix C.4, we prove that the corresponding Algorithm 4, H +- REDUCE(H, e)). returns Hermitian matrix H such that, for any given accuracy e (0, 1),
( ( H)) S e, Vi E [n.]
The symmetry induced by the Cholesky transformation is crucial to design an efficient algorithm for the spectral gap. As described in Section 3, and analyzed in Appendices D and E,, any spectral gap eigenvalue of Hermitian definite pencil can be approximated by iterative algorithm that uses only "counting-quere querries" i.e., queries that ask how many eigenvalues are smaller than a given threshold. This way we completely avoid diagonalization, thus leading lower complexity.
To perform the counting queries efficiently, the transformed matrix H must be regularized with small random perturbations, in the spirit of smoothed analysis [128), which has recently drawn attention in the context of matrix algorithms [31, 25, 15, 41, 104] (see Appendix D for the analysis). These aforementioned works typically require a guarantee on the minimum eigenvalue gap of the perturbed matrix, e.g., [104] uses a Minami-type bound [105], while in [15, 41] the entire pseudospectrum of the perturbed matrix must be shattered with respect to a grid. The latter is even more challenging to achieve than a minimum gap and it requires Ō(n²) random bits. Our algorithm is significantly less demanding in terms of randomness: All we need is the Wegner estimate [138, 4] for the density-of- states of random Hermitian operators, and only Ö(n) random bits in total.
Finally, in Section 4, we apply our main results to prove the first matrix multiplication-yyye upper bounds for the bit complexity of PCA algorithms. Specifically, for the standard PCA formulation, we show that we can first compute the spectral projector and then use deflation to obtain a basis for the desired diwedinesli embedding in nearly matrix multiplication time. We then apply similar arguments to the seminal Bhock-Krylov PCA algorithm of [106].

---

PROJECTOR
Input: Hermitian definite pencil H E H"" E + gap index accuracy 6.
Requires:
|||] 1, ||- k, E (n- 1].
Algorithm: II PROJECTOR(H.SSKK c).
1: M gapp <- GAP(H, S, g: ).
2: K < COND(S, 2 )-
b This is skipped when = I.
3: Íî PURIFY(H, S, Pk, gap,, R).
4: return II.
Output: Approximate projector II the invariant subspace associated with the & smallest eigenvalues.
Ensures:
| II] I with probability at least 1/n.
Algorithm 1: PROJECTOR.
CHOLESKY
Input: Matrix M = (B
C E H
Requires: is positive-definite u S "M)))=((M)) log(n) for some constants cı, c2, ca and e E (0, 1).
Algorithm: L < CM)LEEKY(M)..
if
2:
return L MM,1.
3: else
4:
L1 + (AOLESKY(A))
b LIi Li
+ E
5:
BAi + MM (B, INV (A))
BAi = B (A -1 + NNV + MM
BA - ERA
6:
L21 <- MM (.),,L))
L21
(BA EBA:) L1 + EM
BA - L1 E E1
7:
5 C HEMMMM
S C ((BAi) B EMM) E
Es
where Es BET B EM B + +E
E
8: L2a < CHOLESKY
LoL = S E
Lu
9: return L
Lai
L2a
Output:
Lower triangular Cholesky factor L E C
Ensures:
LL'"
M S e|[M], L is lower triangular, LL* Hermitian and positive-definite
Algorithm 2: CHOLESKY(M))
1.4 Notation
Matrices denoted by bold capital letters and vectors by bold small letters. For real complex constants we typically use Greek letters, or the Latin letters C. The vector e; denotes the i-th column of the standard basis. A is the conjugate transpose of A and A denotes the pseudoinverse. The 2-norm is the default for matrices and vectors. <(A) = ||||||| is the two-norm condition number of A. For the етror analysis of the various algorithms, we use E.P to denote the error matrices that are introduced by the floating point errors of the i-th operation OP. The letters c and 8 typically denote (scalar) eтrоr quantities and failure probabilities, respectively. [n] is the set {1, 2..,,... n1), ). We denote by H! C IR" the set of Hermitian matrices of size x n, HH the set of Hermitian positive semi-definite matrices and H3 the set of Hermitian positive definite matrices. For matrix and scalar z we write z A a shorthand for 2I A. A(A) and A(A, B) denote the spectrum of matrix A and matrix pencil (A, B). respectively. The eigenvalues and

---

singular values are always sorted in ascending order by default: X S A2 .... Ax. Ac(A) is the e-pseudospectruum of A (see Definition A.2).
2 Computing spectral projectors with the sign function
Given Hermitian definite pencil (H, S). our ultimate goal is compute forward error approx- imation of the spectral projector that is associated with the k smallest eigenvalues, as described in Problem 1.1. Algorithm 3 solves this problem provably and efficiently, but it requires that we al- ready have a suitable approximation of the eigenvalue gap that separates the desired subspace from the rest of the eigenspace. The algorithm called PURIFY since it is inspired by "purification" techniques in DFT, referring to the removal of the unoccupied orbitals. The computation of the gap and the midpoint is in fact the bottleneck of our main algorithm, however, we still show that they can be computed efficiently in Section 3, and, importantly, without diagonalizing any matrices. The properties of Algorithm are stated in Proposition 2.1.
Proposition 2.1. Lei E IHT with ||]]| S 1, S H with || 1 k E [n-1] and E (0,1). 1).
Let l4 A and gap; A Akt1, where À1 .... An are the generalized eigenvalues of the Hermitian definite pencil (H. S) and assume that we want to compute II, which is the true spectral projector associated with the k smallest eigenvalues. we have access to for some constant C' > 1, then Algorithm 3 computes II + PURIFY S, jık, gapk, R. such
Th E Pb + gaps
gapę E (1 5 gapk:
K E [N(S), Ck()),
that | - 1+ll I in (n) (log( gaps
log(log( ((ss..))) floating point operations
using O (log(n) log"( gaps
1 log(
<(8) egupp bits of precision.
Proof. The full proof of Proposition 2.1 can be found in Appendix B. We briefly summarize it here. The main idea is S S-H If we already know that Pa is reasonable approximation of j4k- and that it is located inside the correct eigenvalue gap, then, in exact arithmetic, our problem is equivalent to computing sgn(px - H). The result can be used to filter the desired spectral projector, often referred as "purification" in the context of DFT. The main challenge is to ensure that all propagated numerical errors, success probabilities, and input parameters for all algorithms are well-defined and bounded. To obtain the final forward errors must rely matrix similarity arguments, the properties of the pseudospectrum, the eigenvalue bounds of Weyl and Kahan from Fact 1, and Lemma B.1, which gives explicit bounds on the sign function under small floating point perturbations.
the sign function algorithm from [15], SGN, approximate sgn(
The rest of the paper is devoted to the analysis of our new algorithm for the spectral gap and the mid- point based on eigenvalue counting queries, described in Theorem 3.1. For comparison purposes, in Appendix B.1 we analyze a diagonalization based algorithm for the same task (which is a new result itself), specifically, using the of-te-ofthhe-aee EIG algorithm of [15]. We compare the two algorithms in Section 3.4, demonstrating that our counting- based algorithm is indeed faster.
3 Fast spectral gaps with counting queries
Our core algorithm efficiently approximates spectral gaps based on "eigenvalue counting queries" only, thus avoiding explicit (and expensive) diagonalization. To give some intuition the main idea, consider the following simplified version of the problem.
Problem 3.1 (Gap finder). Let >1 S... & Xn in [-, be n (unknown) real values (e.g., they can be the eigenvalues of the original matrix pencil) μi xx and gap; = >x++1 - > for some k [n 1). Given and some error parameter e 10 as input, we want to approximate pa and gap; up to additive gaap- i.e. we look for F P4 t gap,, and gap; E (1 + c) gaps- Only queries of the following form can be performed: We fix parameter E (0. 1/2), which distorts all 7,, >;
7). We can then choose any value h E (-1- ask how many values X; are smaller than h. For each Y, we can query arbitrarily many different
X; to some (unknown) X; E [>i-7,
+n) and
values for b, and each h-query costs q(1/>) O(polylog(1//) (1/7)).

---

The query cost is arbitrary to avoid trivial solutions by setting = 0. The following proposition is proved in Appendix E:
Proposition 3.1. Problem 3.1 can be solved iteratively by execuling a total of O(log( Egapp) iter-
ations and (1) queries per iteration, where each query costs at most q gap,
3.1 Smoothed analysis of eigenvalue counting
To use the counting query model of Problem 3.1 and Proposition 3.1 compute the spectral gap of a matrix pencil, we need a "black-box" method to count eigenvalues that are smaller than threshold. We first describe straightforward, deterministic algorithm COUNT(X, h, E) for this task (see Lemma E.I), which takes as input Hermitian matrix X, a scalar h. and parameter 5, X) E. returns the precise number of eigenvalues of X
with the requirement that min that are smaller than h. The runtime of the algorithm depends on log(1/e), and must therefore be minimized. For this we resort to smoothed analysis: We apply a random perturbation to ensure that e is at least polynomial in 1/n, up some other factors detailed in Appendix D.
To build random "regularizer." in Definition D.I we introduce a random oracle that samples num- bers from standard normal distribution and returns their floating point representation using pre- specified number of bits. Based on this simple oracle, we can design a floating point algorithm X f REGULARIZE(A 7, 7) which has the following properties described Proposition 3.2:
Proposition 3.2. Ler with ||||| S 1 be a Hermitian matrix, 7, 8 E (0, 1/4) two given parameters, and X + REGULARIZE(A, 7, Let g be an arbitrary (but fixed) grid of points in [-2, 2] with cardinality gl = T. For every element húg & consider the matrices M, h4 X and M
h4 X Eg, where E, denote the diagonal floating point error matrices induced by the shift. All the following hold simultaneously with probability 1 - 28 if we use O(log( 76 )) bits of precision:
X < 4/3,
1 M(AA)
S
9 16
mmin (M.)
78
4n 14T
Proof. The main result that we use the proof can be traced back to the Wegner estimate for the ofnsity-of-states of Hermitian operators under random diagonal disorder [138]. See Appendix D and in particular D.2 for more details.
3.2 Computing the gap and the midpoint
We can describe the algorithm GAP in Theorem 3.1, that computes the k-th gap and the mid- point a Hermitian definite pencil. The same methodology can be extended to approximate any singular value, described in Proposition E.2 in Appendix E.3.
Theorem 3.1 (GAP). Let H E H", E H and |H]] ||||||| S 1, which define a Hermitian definite pencil (H,S). Given &s E [n - 1), accuracy E (0, 1), and failure probability 8 E (0, 1/2), there exists an algorithm
Pk. gap + GAP(H,S,k,er
which returns juk pk + egap, and gap (1 + E) gapk, where pa
As A and gap,
d Akt1. The algorithm requires
(TTMM(1) log( de gegm log( (ggge)
(aegap. (xr log (H(()) bits, where X; are the eigen- values of (H, If ((S) is unknown, additional ( (n) log( mm(S)) log(n(S))) floating point operations and (log(n) log* ()n(s)) bits are sufficient to compute it with Corollary E.I.
arithmetic operations using (log(n) (log"
Proof. The full proof builds upon the results that are detailed in Appendices D and E. A summary is the following. We first fix our initial error parameter 70 1/8 and call H REDUCE (H, S, 30
7

---

(Algorithm 4) which internally uses CHOLESKY to reduce the GEP to a regular Hermitian one. From Proposition C.3, Hermitian matrix H is returned such that [^, (H) X; (H. S 4
Next, we use the same counting query model as in Proposition 3.1. We first regularize H using X + REGULARIZE(H, Iin 2 60), where 6. = 6/2 is the initial failure probability, Conditioning on success of Proposition 3.2 (with probability 1 50), for all E [n]. it holds that (X, (H) 9 Summing the two eigenvalue error bounds, we conclude that all eigenvalues of (H. S), which initially lie in [. are distorted by at most "ro in X. We now have all necessary tools to go back to the counting query model of Proposition 3.1: In the first step we construct a grid g = 1.9/8). Clearly, igl (()). Since we conditioned on the success of
- -7/8, 7/8, .6,
.....
E) 2 éo Proposition 3.2, the regularization ensures that for every h2 E g it holds that m (hjn X +
with Eo = 8 a (n)
7006
This allows us to efficiently execute COUNT(X, h2, Eu) for every hj-
At the end of the first iteration, we have computed two intervals Ig and I++1, where I contains A and 1 contains Ak+1, and each interval has size at most 1, i.e., half the size of [-1, 1]. We continue by halving at each step both 7 and 6, constructing the corresponding grids as per the proof of Proposition 3.1, and counting eigenvalues over the grid. In each iteration after the first one, we keep track of two intervals I and fx+1, and two corresponding grids go and gk-1 with size ge = gk+1| = (1). We therefore only need to execute constant number of COUNT queries, and in each iteration the size of the intervals I and 1x+1 is halved. The algorithm terminates
after a total of m O (log(= ) iterations and finally provides the advertised complexity, bit requirements, failure probability, and approximation guarantees.
3.3 Sketch proof Theorem 1.1
The proof of our main Theorem 1.1 directly follows from Theorem 3.1 together with Proposition 2.1 as well as the algorithm SIGMAK (described in Appendix E.3) which used to compute the condition number of S. The full proof can be found in Appendix E.4.
3.4 Comparison with diagonalization
We can now compare Theorem 3.1 with diagonalization- -hased approach that detailed in Propo- sition B.2. We fix 8 O(1/n) so that both algorithms succeed with the same probability.
((poly(L/n)) the total arithmetic complexity of the algorithm of Theo- rem 3.1 is OTTM (n) lo))(n)) using O(log bits. For the same parameters, Proposition B.2 requires need a total of (TTMM (n) log'(()) arithmetic operations, and Ollog (n)) bits. Thus, in total, Algorithm 3.1 is faster by factor of (llg(n))).
For e, gap K (S)
In the extreme case where c, gapa, R(S) (11, Theorem 3.1 counts ( (n) log(n)) arithmetic operations and O(lgg" (n)) bits, while Proposition B.2 requires (TTMa(nn) log2(n) (n)) operations, and Ollog((m) (n)) bits. Thus, Proposition B.2 is again slower by factor of Ollog(n)). Interestingly, in this case Theorem 3.1 is faster than even single call to EIG, which requires (TTi((n) log (n)) arithmetic operations. We conclude that, at least based on the currently existing algorithms, diago- nalization is redundant for the computation of spectral gaps and invariant subspaces.
3.5 Application in DFT
In Appendix F we demonstrate how our main results can be directly applied to approximate density matrices and electron densities of atomic systems DFT. Even though is not a machine learning problem per se, we decided to dedicate section in the Appendix due to its importance: DFT calculations persistently occupy supercomputing clusters and the corresponding software libraries and literature receive tens of thousands of citations annually at the time of this writing [89, 59, 60, 127, 74]. Our work the first analysis to provide forward-error guarantees in finite precision for these problems in nearly matrix multiplication time.
8

---

4 PCA
Since its introduction in the early twentieth century [116, 72], Principal Component Analysis is one of the most important tools in statistics, data science, and machine learning. It can be used, for example, to visualize data, to reduce dimensionality, or to remove noise from data; cf. [79, 45] for reviews on the vast bibliography. In its simplest formulation, given (centered) data matrix X Rmxx, the goal is to learn a k-dimensional embedding Cs. where & < n, that maximizes the sample variance, which can be written as optimization problem
C+ arg max tr(C' HC),
(3)
C C=
where H = X E 1Tx" is the sample covariance. It can be shown that the solution C₁ corre- sponds to the principal k singular vectors of H. i.c. the ones that correspond to the largest singular values. Evidently, since the sample covariance is always symmetric and positive semi-definite, this can be written as a Hermitian eigenvalue problem
HC CA,
(which is indeed definite GEP in Equation (1) with S = 1). By solving for Ck, can project the data in h dimensions by computing XC2, preserving as much of the variance in & dimensions main results. However, the solution of Equa-
as possible. To compute C₂ we can directly tion (3) is an actual orthonormal basis for the invariant subspace rather than the spectral projector that Theorem 1.1 returns. This can be addressed with deflation: Once we have the spectral projec- tor I, assuming that the approximation is sufficiently tight, we can apply a subsequent deflation step based on rank-revealing QR factorization to obtain k-dimensional basis. This can be done deterministically in O(n²) time [66] or in randomized (n*) [39, 15).
In Appendix G.I we prove the following Theorem 4.1 for Algorithm 7, which builds upon our main Theorem 1.1, the algorithm of Proposition E.2 approximate || XxIl, and the DEFLATE algorithm of [19], to solve the standard PCA problem of Eq. (3). Following the existing literature, the result is stated for real matrices, but it can be trivially adapted the complex case as well.
Theorem 4.1 (PCA). Let X E Rmxn be centered data matrix, H the n x n symmetric sample covariance matrix, i.e., H = x H 1, k E [ni] target rank, and E (0,1) an accuracy
parameter. Given H. we compute matrix C₁ with lo columns such that |X - XC C2 s
(1 XXCCCK |l, where Cj E Ban contains the top-k (right) singular vectors of X o (Tmm((n) (log(
n log( 1 log(log( co (app)) n own log( oe log( n gap,
arithmetic operations using O (log()) (log"(.gm,)) n sion, with probability at least 1 - O(1/n).
(egip) ) + +log (om1)) + log( Log(roin)) bits of preci-
4.1 Block-Krylov PCA
In some applications, the target dimension k might be small, i.c., k < n. This condition has driven a whole area of research in so-called low-rank approximation algorithms for PCA [53, 120, 33, 67, 100, 106, 30, 29, 5]. Such approaches are also suitable for kernel PCA, since they rely on matrix-vector products and therefore the kernel matrix does not need to be explicitly formed. The techniques from the previous section can be directly applied to obtain new bit complexity upper bounds for existing state-of-the-art algorithms, which are typically analyzed in exact arithmetic. They internally rely on the computation of the principal singular vectors of submatrices, which can improved with our methods. Specifically, in Appendix G.2 we summarize a floating point Iteration algorithm (see Algorithm 8), essentially, providing matrix
analysis of the Block- multiplication type upper bound on the bit complexity with only a polylogarithmic dependence on the singular value gap. In nutshell, we directly obtain the following result:
Theorem 4.2 (Bit complexity analysis of Block- Block-Krrlov PCA). Let X be data matrix X E RRmx log(n) VePCA
| < 1, Ae E [n] a target rank, éPCA E (0, an accuracy parameter, and =
Let
TMMx(k) denote the complexity stably multiply X X with dense matrix with k columns

---

from the right (see Def. G.I). Using the Steps 1-6 that are detailed in Appendix G.3 as floating point implementation of Algorithm 8, compute matrix Zy E mRx that satisfies
Z Z S O(epca);
with high probability, where Zy is an approximate basis for the top-k principal components of X, returned by Algorithm 8 in exact arithmetic. The total cost is at most
o o (TTTaax())) log(
x(K) (MP) + mm(gk)) log( gapslMy + (gk)) polylog( (MP)
ole
floating point operations, using O (polylog( mgx(() EPCA
) bits of precision. K, M are in Alg. 8.
Proof. The full proof can be found in Thm. G.2, Appendix G.3. The main idea is to apply the count- ing query methodology compute the condition number of the Block- Krylovv matrix K, as well as the k-th spectral gap and the midpoint of the reduced matrix M in Line 5 of Alg. 8. Thereafter, we can compute spectral projector and an approximate basis for the top-k singular vectors of M using PURIFY and DEFLATE, similar to the analysis of classical PCA in the previous section.
5 Open problems
We mention some open problems and interesting future directions.
(i) Bit requirement of SGN: The major bottleneck for the bit requirements of our main algo- rithms comes from the SGN algorithm of [15]. An inverse-free Newton-Schultz iteration [83], or the implicit repeated squaring of [41, 121] can potentially give significant improvements.
(ii) Sparse algorithms: In applications like DFT it commonly occurs that the matrices have spe- cial structure, i.e., they are banded and/or sparse. remains open whether Problem 1.1 can be provably solved faster than reported results in finite precision for these special cases (recall that the tridiagonal QR algorithm requires O(n²) operations to return the eigenvectors). An end-to-end stability analysis of existing fast eigensolvers would be the place start [65, 137].
(iii) Distributed PCA:
The techniques for Block-Krylov PCA be potentially applied dis- tributed or streaming PCA algorithms, which are also based on randomized low-rank approx- imations. E.g., the distributed PCA algorithm of [31], straightforward replace the SVD computation on the server by counting query iteration. The full analysis of such an approach is left future work.
Conclusion
In this work we provided an end-to-end analysis to approximate spectral projectors on k- dimensional invariant subspaces of Hermitian definite matrix pencils (H,S) that require at most O (Z) polylog(n, - <(S), gap) ") bit operations in the floating point model of computa- tion. This is the first end-to-end analysis that improves the Ö(n3) complexity of classic eigensolvers for both the regular and the generalized case. To achieve this result introduced new method to approximate spectral gaps by querying the number of eigenvalues that are smaller than a threshold, and therefore completely avoid explicit diagonalization of any matrix pencil. This approach quired proving that the Cholesky factorization can be stably computed in O(TMM (n)) floating point operations, novel result per se. Our results have direct implications PCA problems, providing matrix multiplication type upper bounds for the bit complexity of classical and Block-Krylov PCA.
Acknowledgements
We would like to thank Daniel Kressner, Nian Shao, Ryan Schneider, Nicolas Deutschmann, Nikhil Srivastava, Ilse Ipsen, Cameron Musco, David Woodruff, Uria Mor, and Haim Avron for helpful discussions.
10

---

References
[1] Dimitris Achlioptas. Database-friendly random projections. In Proc. Twentieth ACM SIGMOD-SIGACT-SICANN Symposium on Principles of Database Systems, pages 274-281, 2001. 61
Satoru Adachi, Satoru Iwata, Yuji Nakatsukasa, and Akiko Takeda. Solving the trust-region subproblem by generalized eigenvalue problem. SIAM Journal on Optimization, 27(1):269-
291, 2017.
[3] Nir Ailon and Bernard Chazelle. Approximate nearest neighbors and the fast Johnson- Lindenstrauss transform. In Proc. Thirty-Fighth Annual ACM Symposium on Theory of Com-
puting, pages 557-563, 2006. 62
[4] Michael Aizemman, Ron Peled, Jeffrey Schenker, Mira Shamis, and Sasha Sodin. Matrix reg- ularizing effects of Gaussian perturbations. Communications in Contemporary Mathematics,
19(03):1750028, 2017. 4, 47
[5] Zeyuan Allen-Zhu and Yuanzhi Li. LazySVD: Even faster SVD decomposition yet without agonizing pain. Advances Neural Information Processing Systems, 29, 2016. 9
[6] Josh Alman, Ran Duan, Virginia Vassilevska Williams, Yinzhan Xu, Zixuan Xu, and Renfei Zhou. More asymmetry yields faster matrix multiplication. arXiv preprint 2404 2024. 2
[7] Diego Armentano, Carlos Beltrán, Peter Bürgisser, Felipe Cucker, and Michael Shub, A stable, polynomial-time algorithm for the eigenpair problem. Journal of the European Math- ematical Sociery, 20(6):1375-143) 2018. 26
[8] Zhaojun Bai and James Demmel. On a direct algorithm for computing invariant subspaces with specified eigenvalues. University of Tennessee, Computer Science Department, 1991. 2,
26
[9] Zhaojun Bai and James Demmel. Using the matrix sign function to compute invariant sub- spaces. SIAM Journal Mairix Analysis and Applications, 19(1):005-225, 1998. 2, 26,
[10] Zhaojun Bai, James Demmel, and Ming Gu. An inverse free parallel spectral divide and conquer algorithm for nonsymmetric eigenproblems. Numerische Mathematik, 76(3):279- 308, 1997. 2, 26
[11] Gökhan H Bakur, Jason Weston, and Bernhard Schölkopf. Learning to find pre-images. Ad-
vances in Neural Information Processing Systems, 16:449-456, 2004. 1
[12] Pierre Baldi and Kurt Hornik. Neural networks and principal component analysis: Learning from examples without local minima. Neural networks, 2(1):53-58, 1989.
[13] Grey Ballard, James Demmel, and Ioana Dumitriu. Minimizing Communication for Eigen- problems and the Singular Value Decomposition. Technical Report UCB/EECS-2011-14, February 2011.
[14] Grey Ballard, James Demmel, Olga Holtz, and Oded Schwartz. Minimizing Communication in Numerical Linear Algebra. SIAM Journal on Matrix Analysis and Applications, 32(3):866~
901, 2.11.
[15] Jess Banks, Jorge Garza- Varrza, Archit Kulkarni, and Nikhil Srivastava. Pseudospectral Shat- tering, the Sign Function, and Diagonalization in Nearly Matrix Multiplication Time. Foun- dations of Computational Mathematics, pages 1-89, 2022. 2, 3, 4, 6, 9, 10, 19, 21, 22, 23,
26, 27, 31, 32, 33
[16] Jess Banks, Jorge Garza- Garza-Vargas, and Nikhil Srivastava. Global Convergence of Hessenberg Shifted QR II: Numerical Stability. arXiv preprint (iXX::2205.06//1 2022. 26
[17] Jess Banks, Jorge Garza- Garza-Vargas, and Nikhil Srivastava. Global Convergence of Hessen- berg Shifted QR Approximate Ritz Values via Shifted Inverse Iteration. arXiv preprint 0rXiv:2205.06804 2022. 26
[18] Jess Banks, Jorge Vargas, and Nikhil Srivastava. Global Convergence of Hessenberg Shifted QR I: Exact Arithmetic. Foundations of Computational Mathematics, pages 1-34, 2024. 26
11

---

[19] Jess Banks, Archit Kulkarni, Satyaki Mukherjee, and Nikhil Srivastava. Gaussian Regulariza- tion of the Pseudospectrum and Davies' Conjecture. Communications on Pure and Applied Mathematics, 74(10):2114-2133 2021. 9, 62
[20] Elnaz Barshan, Ali Ghodsi, Zohreh Azimifar, and Mansoor Zolghadri Jahromi. Supervised principal component analysis: Visualization, classification and regression on subspaces and submanifolds. Pattern Recognition, 44(7):1357-1371, 2011.
[21] Yair Bartal, Ben Recht, and Leonard J Schulman. Dimensionality reduction: beyond the Johnson-Lindenstrauuss bound. In Proc. 2011 Annual ACM-SIAM Symposium on Discrete Algorithmss, pages 868-887 SIAM, 2011. 62
[22] Harold Basch, C. J. Hornback, J. W. Moskowit.. Gaussian-Orbital Basis Sets for the First- Row Transition-Mett Atoms. The Journal of Chemical Physics, 51(4):1311-1318, 1969. 58
[23] Robert A Baston and Yuji Nakatsukasa. Stochastic diagonal estimation: probabilistic bounds and an improved algorithm. arXiv preprint wrxiv:221..16644 2022. 62
[24] Friedrich L Bauer and Charles T Fike. Norms and exclusion theorems. Numerische Mathe-
matik, 2:137-141, 1960. 20
Proc. 9th Innovations in Theoretical Computer Science Conference, pages 6:1-6:22. Schloss [25] Michael Ben-Or and Lior Eldar. A Quasi-Random Approach to Matrix Spectral Analysis. In
Daggtuhl-Laibniz-Zeenttu fuer Informatik, 2018. 2, 4, 26
[26] Peter Benner, Volker Mehrmann, and Hongguo Xu. A new method for computing the sta- ble invariant subspace of real Hamiltonian matrix. Journal of computational and applied mathematics, 86(1):17-43, 1997. 26
[27] Rajendra Bhatia. Pinching, trimming, truncating, and averaging of matrices. The American Mathematical Monthly, 107(7):602 608, 2000. 25
[28] Rajendra Bhatia. Perturbation bounds for matrix eigenvalues. SIAM, 2007. 20
[29] Srinadh Bhojanapalli,, Prateek Jain, and Sujay Sanghavi. Tighter low-rank approximation via sampling the leveraged element. In Proc. Twenty-Sixth Annual ACM-SIAM Symposium on Discrete Algorithms,, pages 902-920. SIAM, 2014.
[30] Christos Boutsidis, Petros Drineas, and Malik Magdon-Ismail. Near-optimal column-based matrix reconstruction. SIAM Journal on Computing, 43(2):687-717, 2014. 9
[31] Christos Boutsidis, David P Woodruff, and Peilin Zhong. Optimal principal component anal- ysis in distributed and streaming models. In Proc. Forty-Eighth Anmual ACM Symposium on Theory of Computing, pages 236-249, 2016. 4, 10
[32] Tim Clark and Rainer Koch. Linear Combination of Atomic Orbitals, pages 5-22. Springer Berlin Heidelberg, 1999. 58
[33] Kenneth I. Clarkson and David P Woodruff. Numerical linear algebra in the streaming model. In Proc. Forty-First Annual ACM Symposium on Theory of Computing, pages 205-214, 2009.
9, 64
[34] Kenneth L Clarkson and David P Woodruff. Low-rank approximation and regression in input sparsity time. Journal of the ACM, 63(6):1-45, 2017. 64
[35] Michael B Cohen, Sam Elder, Cameron Musco, Christopher Musco, and Madalina Persu. Dimensionality reduction for k-means clustering and low rank approximation. In Proc. Forty-
Seventh Annual ACM Symposium on Theory of Computing, pages 163-172, 2015.
[36] Charles R. Crawford. Reduction of Band-symmetric generalized eigenvalue problem. Com- munications of the ACM, 16(1):41-44 44, 1973. 26
[37] Philip I Davies, Nicholas J Higham, and Françoise Tisseur. Analysis of the cholesky method with iterative refinement for solving the symmetric definite generalized eigenproblem. SIAM Journal on Matrix Analysis and Applications, (72)):7224493 2001. 2, 4, 26
[38] Theodorus J Dekker and Joseph F Traub. The shifted QR algorithm for Hermitian matrices. Linear Algebru and its Applications, 4(3):137-154, 1971. 26
[39] James Demmel, Ioana Dumitriu, and Olga Holtz, Fast linear algebra is stable. Numerische Mathematik, 108(1):59-9)) 2007. 2, 3, 4, 9, 21, 22, 23, 34, 47, 64
12

---

[40] James Demmel, Ioana Dumitriu, Olga Holtz, and Robert Kleinberg. Fast matrix multiplica- tion is stable. Numerische Mathematik, 106(2):199-22) 2007. 2, 22,
[41] James Demmel, Ioana Dumitriu, and Ryan Schneider. Generalized Pseudospectral Shattering and Inverse-Free Matrix Pencil Diagonalization. arXiv preprint WrXiv:2306.037000 2023. 2,
3,4, 4, 10, 23
[42] James W Demmel. Applied numerical linear algebra. SIAM, 1997. 2,26
[43] James Weldon Demmel. Three methods for refining estimates of invariant subspaces. Com-
puting, :3355,, 1987. 26
[44] James Weldon Demmel and Bo Kägström. Computing stable cigendecompositions of matrix pencils. Linear Algebra and its Applications, 88: 18:1391186, 1987. 26
[45] Konstantinos I Diamantaras and Sun Yuan Kung. Principal component neural networks: theory and applications. John Wiley & Sons, Inc., 1996. 1, 9, 62
[46] Petros Drineas, Iordanis Kerenidis, and Prabhakar Raghavan. Competitive recommendation systems. In Proc. Thiry-Fourth Anmual ACM Symposium on Theory of Computing, pages 82-90, 2002. 1
[47] Ran Duan, Hongxun Wu, and Renfei Zhou. Faster matrix multiplication via asymmetric hashing. In 2023 IEEE 64th Annual Symposium Foundations of Computer Science, pages 2129-2138. IEEE, 2023. 2, 21
[48] Stanley C Eisenstat and Ilse CF Ipsen. Relative perturbation results for eigenvalues and eigenvectors of diagonalisable matrices. BIT Numerical Mathematics, 38(3):502-509, 1998.
20
[49] Ethan N Epperly, Joel A Tropp, and Robert J Webber. XTrace: Making the Most of Every Sample in Stochastic Trace Estimation. SIAM Journal on Matrix Analysis and Applications, 45(1):1-23, 2024. 62
[50] Dan Feldman, Melanie Schmidt, and Christian Sohler. Turning big data into tiny data: Constant-size coresets for k-means, PCA, and projective clustering. SIAM Journal on Com- puting, 49(3)x601-657, 2020.
[51] John GF Francis. The QR transformation a unitary analogue to the LR transformation Part
1. The Computer Journal, 4(3):265-271, 1961. 26
[52] John GF Francis. The QR transformation Part 2. The Computer Journal, 4(4):332-345,
1962. 26
[53] Alan Frieze, Ravi Kannan, and Santosh Vempala. Fast monte-carlo algorithms for finding low-rank approximations. In Proc. 39th Annual Symposium on Foundations of Computer Science, page 370, 1998. 9, 64
[54] Alan Frieze, Ravi Kannan, and Santosh Vempala. Fast Monte-Carlo algorithms for finding low-rank approximations. Journal of the ACM, (1)) 2004. 64
[55] Martin Fürer. Faster integer multiplication. Theory of Computing, pages 57-66, 2007. 2
Proc. Thirty-Ninth Annual ACM Symposium
[56] Giulia Galli. Linear scaling methods for electronic structure calculations and quantum molec- ular dynamics simulations. Current Opinion in Solid State and Materials Science, 1(6):864- 874, 1996. 59
[57] Giulia Galli and Michele Parrinello. Large scale electronic structure calculations. Physical Review Letters, 69(24)3554,, 1992. 59
[58] Alan George. Nested dissection of regular finite element mesh. SIAM Journal on Numerical Analysis, (0(2))345-366 1973.
[59] P Giannozzi, O Andreussi, T Brumme, O Bunau, M Buongiorno Nardelli, M Calandra, R Car, C Cavazzoni, Ceresoli, M Cococcioni, N Colonna, Camimeo, A Dal Corso, S de Giron- coli, P Delugas, R A DiStasio, A Ferretti, A Floris, G Fratesi, G Fugallo, R Gebauer, U Gerst- mann, Giustino, Gorni, J Jia, M Kawamura, H-Y Ko, A Kokalj, E Küçükbenli, M Lazzeri, M Marsili, N Marzari, F Mauri, N L Nguyen, H-V Nguyen, A Otero de-la Roza, L. Paulatto, S Poncé, D Rocca, Sabatini, B Santra, M Schlipf, A Seitsonen, A Smogunov, Tim- rov, T Thonhauser, P Umari, N Vast, X Wu, and S Baroni. Advanced capabilities for
13

---

materials modelling with Quantum ESPRESSO. Journal of Physics: Condensed Matter, 29446::46590, oct 2017. 58
[60] Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli, Guido L Chiarotti, Matteo Cococcioni, Ismaila Dabo, An- drea Dal Corso, Stefano de Gironcoli, Stefano Fabris, Guido Fratesi, Ralph Gebauer, Uwe Gerstmann, Christos Gougoussis, Anton Kokalj, Michele Lazzeri, Layla Martin-Samos, Nicola Marzari, Francesco Mauri, Riccardo Mazzarello, Stefano Paolini, Alfredo Pasquarello, Lorenzo Paulatto, Carlo Sbraccia, Sandro Scandolo, Gabriele Sclauzero, Ari P Seitsonen, Alexander Smogunov, Paolo Umari, and Renata M Wentzcovitch. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. Journal of Physies: Condensed Matter, 21(39):395502, Sep 2009. 8
[61] John R Gilbert and Robert Endre Tarjan. The analysis of nested dissection algorithm. Numerische Mathematik, 50(4):377-404, 1986. 4
[62] Stefan Goedecker. Linear scaling electronic structure methods. Reviews of Modern Physics,
71 114): 1999. 59
[63] Stefan Goedecker and GE Scuserza. Linear scaling electronic structure methods in chemistry and physics. Computing in Science & Engineering, 5(4):14-21, 2003. 59
[64] Gene H Golub and Charles F Van Loan. Matrix Computations. Johns Hopkins University
Press, 2013. 2. 19, 26, 64
[65] Ming Gu and Stanley C Fisenstat. A divide-and-comquer algorithm for the symmetric tridi- agonal eigenproblem. SIAM Journal on Matrix Analysis and Applications, 16(1):172-191, 1995. 10, 59
[66] Ming Gu and Stanley C Fisenstat. Efficient algorithms for computing strong rank-revealing QR factorization. SIAM Journal on Scientific Computing. 17(4):848-888 1996. 9
[67] Nathan Halko, Ber-Gunnar Martinsson, and Joel A Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM Review,
53(2):217-288, 2011. 9, 64
[68] David Harvey and Joris Van Der Hoeven, Integer multiplication in time o(nlog) .n). Annals of Mathematics, 193(2):563-611) 2021. 2
[69] Qiu He, Bin Yu, Zhaohuai Li, and Yan Zhao. Density Functional Theory for Battery Materi- als. Energy & Environmental Materials, 2(4):264- 279.,, 2019. 58
[70] Nicholas J Higham. Accuracy and stability of numerical algorithms. SIAM, 2002. 19, 20,
61, 64
[71] Nicholas J Higham. Cholesky factorization. Wiley Interdisciplinary Reviews: Computational Statisties, 1(2):251-22) 2009.
[72] Harold Hotelling. Analysis of complex of statistical variables into principal components.
Journal of educational psychology, 24(6):417, 1933. 9, 62
[73] Edward J Hu, yelong shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Lu Wang, and Weizhu Chen. LoRA: Low-rank adaptation of large language models. In International Conference on Learning Representattonss 2022.
[74] Jürg Hutter, Marcella lannuzzi, Florian Schiffmann, and Joost Vondelo cp2k: atomistic simulations of condensed matter systems. Wiley Interdisciplinary Reviews: Computational
Molecular Science, 4(1):15-25, 2014. 8, 58
[75] Ilse CF Ipsen. Computing an eigenvector with inverse iteration. SIAM Review, 39(2):254
291, 1997. 2, 26
[76] Ilse CF Ipsen. Absolute and relative perturbation bounds for invariant subspaces of matrices.
Linear Algebra and its Applications, (09(1-3):45-6 2000. 26
[77] Zhongxiao Jia and GW Stewart. An analysis of the rayleigh-rit method for approximating
eigenspaces. Mathematics of computation, 70(234):637-647, 2001. 26
[78] William B Johnson and Joram Lindenstrauss. Extensions of Lipschitz mappings into a Hilbert
space. Contemp. Math., 26(1):189-206, 1984. 61
[79] Ian T Jolliffe, Principal component analysis for special types of data. Springer, 2002. 1, 9,
62
14

---

[80] Praneeth Kacham and David P Woodruff. Faster Algorithms for Schatten-p Low Rank Ap- proximation. Approximation, Randomization, and Combinatorial Optimization. Algorithmss and Techniques, 2024. 64
[81] William Kahan. Spectra of nearly Hermitian matrices. Proc. American Mathematical Society,
48(1):11-17, 1975. 20
[82]
Daniel M Kane and Jelani Nelson. Sparser Johnson-LLnndenstraus transforms. Journal of the ACM, (11) 2014, 62
[83] Charles S Kenney and Alan J Laub. The matrix sign function. IEEE transactions on automatic
control, 4003):133--1348, 1995. 10
[84] Andrzej Kielbasinskki. A note rounding- error analysis of Cholesky factorization. Linear Algebra and its Applications, 88:487 1987. 4
[85] Cedric Klinkert, Aron Szabó, Christian Stieger, Davide Campi, Nicola Marzari, and Mathieu Luisier. 2-d Materials for Ultrascaled Field-Effect Transistors: One Hundred Candidates under the Ab Initio Microscope. ACS Nano, 14(7):8005-615,, 2020. 58
[86] Walter Kohn. Density functional/Wanmier function theory for systems of very many atoms. Chemical Physics Letters,
1993. 59
[87] Walter Kohn. Density functional and density matrix method scaling linearly with the number of atoms. Physical Review Letters, 76(17):3168, 1996. 59
[88] Walter Kohn and Lu Jeu Sham. Self-consistent equations including exchange and correlation effects. Physical Review, 140(4A):A1133, 1965. 2,
[89] Georg Kresse and Jürgen Furthmüller. Efficient iterative schemes for ab initio total-eneryy calculations using a plane-wave basis set. Physical Review B, 1666))111199111111 1996. 8, 58
[90] Vera N Kublanovskaya. On some algorithms for the solution of the complete eigenvalue problem. USSR Computational Mathematics and Mathematical Physics, 1(3):637-657, 1962. 26
[91] Rasmus Kyng, Yin Tat Lee, Richard Peng, Sushant Sachdeva, and Daniel Spielman. Spar- sified Cholesky and Multigrid Solvers for Connection Laplacians. In Proc. 48th Anmual ACM
Symposium on Theory of Computing, pages 842-850, 2016. 4
[92] Rasmus Kyng and Sushant Sachdeva. Approximate Gaussian Elimination for Laplacians Fast, Sparse, and Simple. In Proc. IEEE 57th Annual Symposium on Foundations of Computer Science, pages 573-582, 2016. 4
[93] Kasper Green Larsen and Jelani Nelson. Optimality of the Johnson-Lindeessraaus lemma. In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science, pages 633-638. IEEE, 2017. 62
[94] Richard J Lipton, Donald J Rose, and Robert Endre Tarjan. Generalized nested dissection. SIAM Journal on Numerical Analysis, 16(2):346-335 1979. 4
[95] Anand Louis and Santosh S Vempala. Accelcrated Newton iteration for roots of black box polynomials. In Proc. IEEE 57th Annual Symposium on Foundations of Computer Science, pages 732-740. IEEE, 2016. 2,3
[96] Amaz Malhi and Robert X Gao. PCA-based feature selection scheme for machine defect classification. IEEE transactions on instrumentation and measurement, 5316:1517-1555, 2004.
[97] Alexander N Malyshev. Computing invariant subspaces of regular linear pencil of matrices.
Siberian Mathematical Journal, 30(4):559 1989. 26
[98] Alexander N Malyshev. Parallel algorithm for solving some spectral problems of linear alge-
bra. Linear Algebra and its Applications, 188:489-520, 1993.
[99] Olvi I. Mangasarian and Edward W Wild. Multisurface proximal support vector machine classification via generalized eigenvalues. IEEE transactions on pattern analysis and machine intelligence, 28(1):69-74, 2005.
[100] Per-Gunnar Martinsson, Vladimir Rokhlin, and Mark Tygert. A randomized algorithm for the decomposition of matrices. Applied and Computational Harmonic Analysis, 30(1):47-
68, 2011. 9, 64
15

---

[101] Nicola Marzari, Arash A. Mostofi, Jonathan R. Yates, Ivo Souza, and David Vanderbilt. Max- imally localized Wannier functions: Theory and applications. Reviews of Modern Physies, 84(4):1419, 2012. 58
[102) Reinhard J Maurer, Christoph Freysoldt, Anthony M Reilly, Jan Gerit Brandenburg, Oliver Hofmann, Torbjörn Björkman, Sébastien Lebègue, and Alexandre Tkatchenko. Advances in Density-Functien Calculations for Materials Modeling. Annual Review of Materials Re- search, 49(1):1-30, 2019. 58
[103] István Mayer. Simple theorems, proofs, and derivations in quanium chemistry. Springer Science & Business Media, 2003. 59
[104] Raphael Meyer, Cameron Musco, and Christopher Musco. On the unreasonable effectiveness of single vector Krylov methods for low-rank approximation. In Proc. 2024 Annual ACM- SIAM Symposium on Discrete Algorithms, pages 81 1-85.. SIAM, 2024. 4
[105] Nariyuki Minami. Local fluctuation of the spectrum of multidimensional Anderson tight
binding model. Communications in mathematical physics, 177:709-777 -250099---,117 1996. 4
[106] Cameron Musco and Christopher Musco. Randomized block Krylov methods for stronger and faster approximate singular value decomposition. Advances in Neural Information Processing Systems, 28, 2015. 4, 9, 64, 65
[107] Cameron Musco, Christopher Musco, and Aaron Sidford. Stability of the lanczos method for matrix function approximation. In Proc. Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1605-1624. SIAM, 2018. 2, 3, 24, 64
[108] Yuji Nakatsukasa. On the condition numbers of a multiple eigenvalue a generalized eigen- value problem. Numerische Mathemarik, 121:531-520 -444,, 2012. 26
[109] Yuji Nakatsukasa and Nicholas Higham. Stable and efficient spectral divide and conquer algorithms for the symmetric eigenvalue decomposition and the SVD. SIAM Journal on Scientific Computing, 3335) A149, 2013. 2
[110] Jelani Nelson and Huy L Nguyên. OSNAP: Faster numerical linear algebra algorithms via sparser subspace embeddings. In Proc. IEEE 54th Annual Symposium on Foundations of Computer Science, pages 117-126. IEEE, 2013. 62
[111] Andrew Ng. Michael Jordan, and Yair Weiss. On spectral clustering: Analysis and an algo- rithm. Advances in Neural Information Processing Systems, 14, 2001. 1
[112] Christopher C Paige. Error analysis of the lanczos algorithm for tridiagonalizing symmetric matrix. IMA Journal of Applied Mathematics, 18(3):34--34, 1976. 26
[113] Christopher Conway Paige. The computation of eigenvalues and eigenvectors of very large
sparse malrices. PhD thesis, University of London, 1971. 2, 26
[114] Victor Y Pan and Zhao Q Chen. The complexity of the matrix eigenproblem. In Proc. 31st Annual ACM Symposium on Theory of Computing, pages 507-516, 1999. 26
[115] Beresford N Parlett. The Symmetric Eigenvalue Problem. SIAM, 1998. 26
[116] Karl Pearson. Liii on lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin philosophical magazine and journal of science, 2(11):559-572, 1901.
9, 62
[117] Andrzej Pokrzywa. Spectra of operators with fixed imaginary parts. Pruc. American Mathe- matical Society, 81(3):359-364 1981. 20
[118] Awwal Mohammed Rufai, Gholamreza Anbarjafari, and Hasan Demirel. Lossy image com- pression using singular value decomposition and wavelet difference reduction. Digital signal processing, 24:117-123, 2014.
[119] Yousef Saad. Numerical methods for large eigenvalue problems: revised edition. SIAM,
2011 1. 26
[120] Tamás Sarlós. Improved approximation algorithms for large matrices via random projections. In Proc. 47th Annual Symposium on Foundations of Computer Science, pages 143-152. IEEE, 2006. 9, 64
[121] Ryan Schneider. When fast, implicit squaring of A-1B stable? arXiv:2310.001990 2023. 2, 3, 10
arXiv preprint
16

---

[122] Bernhard Schölkopf, Alexander Smola, and Klaus-Robert Müller. Nonlinear component anal- ysis as a kemnl eigenvalue problem. Neural computation, (15):1299-1319 1998. 1
[123] Amold Schönhage and Volker Strassen. Fast multiplication of large numbers. Computing, 7:281-292, 1971. 2
[124] Jianbo Shi and Jitendra Malik. Normalized cuts and image segmentation. IEEE Transactions on pattern analysis and machine intelligence, (288 905, 2000.
[125] Max Simchowitz, Ahmed El Alaoui, and Benjamin Recht. Tight query complexity lower bounds for PCA via finite sample deformed wigner law. In Proc. 50th Annual ACM Sympo- sium on Theory of Computing, pages 1249-1259, 2018. 64
[126] Aleksandros Sobczyk and Mathieu Luisier. Approximate Euclidean lengths and dis- tances beyond
Advances in Neural Information Processing Systems,
35:19357-19369, 2022. 62
[127] José M Soler, Emilio Artacho, Julian D Gale, Alberto García, Javier Junquera, Pablo Ordejón, and Daniel Sánchez-Portal. The SIESTA method for ab initio order-N materials simulation.
Journal of Physies: Condensed Matter, 14(1 14(11):2745, 2002. 8, 58
[128] Daniel Spielman and Shang-Hua Teng. Smoothed analysis of algorithms: Why the simplex algorithm usually takes polynomial time. Journal of the ACM, 51(3):885-463, 2004. 4
[129] Nikhil Srivastava. The complexity of diagonalization. In Proc. 2023 International Symposium
on Symbolic and Algebraic Computation, pages 1-6, 2023. 26
[130] Gilbert W Stewart. On the sensitivity of the eigenvalue problem ax=Abx. SIAM Journal on Numerical Analysis, 9(4):669-686, 1972. 26
[131] Gilbert W Stewart. Error and perturbation bounds for subspaces associated with certain eigen- value problems. SIAM Review, 15(4):727-764, 1973. 26
[132] Gilbert W Stewart. On the perturbation of pseudo-inverses, projections and linear least squares problems. SIAM Review, 19(4):6346662, 1977. 26
[133] Gilbert Stewart. Pertubation bounds for the definite generalized eigenvalue problem. Lin-
ear Algebra and its Applications, 33:69-85,1 1979. 26
[134] Volker Strassen. Gaussian elimination is not optimal. Numerische Mathematik, 13(4):354-
356, 1969. 21
[135] Lloyd N Trefethen and Mark Embree. Spectra and Pseudospectra: The Behavior of Nonnor- mal Matrices and Operators. Princeton University Press, 2005. 20, 21, 26
[136] Joost Vondele, Urban Borstnik, and Jurg Hutter. Linear scaling self-consistent field calculations with millions of atoms in the condensed phase. Journal of Chemical Theory and Computation, 8(10):356--335 3573, 2012. 59
[137] James Vogel, Jianlin Xia, Stephen Caulcy, and Venkataramanan Balakrishnan. Superfast divide-and-conquer method and perturbation analysis for structured eigenvalue solutions.
SIAM Journal Scientific Computing, 38(3):A1358-A) (1588-13322, 2016. 10, 59
[138] Franz Wegner. Bounds the density of states in disordered systems. Zeitschrift für Physik
B Condensed Matter, 44(1):9-15, 1981. 4, 4.4
[139] Hermann Weyl. Das asymptotische Verteilungsgesetz der Eigenwerte linearer partieller Dif- ferentialgleichungen (mit einer Anwendung auf die Theorie der Hohiraumstrahlung) Math-
ematische Annalen, 71(4):441-47) 1912. 20
[140] James Hardy Wilkinson. Global convergene of tridiagonal QR algorithm with origin shifts. Linear Algebra and its Applications, 1(3):409 4200 1968. 26
[141] Virginia Vassilevska Williams, Yinzhan Xu, Zixuan Xu, and Renfei Zhou. New bounds for matrix multiplication: from alpha to omega. In Proc. 2024 Annual ACM-SIAM Symposium Discrete Algorithms, pages 3792-3835. SIAM, 2024. 2.
[142] Zewen Xiao and Yanfa Yan. Progress in Theoretical Study of Metal Halide Perovskite Solar Cell Materials. Advanced Energy Materials, 7(22): 2017. 58
[143] Weitao Yang. Direct calculation of electron density in dunsity-l theory. Physical Review Letters, 66(1 1 1): 1438, 1991. 59
17

---

[144] Yunkai Zhou, Yousef Saad, Murilo L Tiago, and James R Chelikowsky. Parallel self- consistent-field calculations via fhebyshev-filteres subspace acceleration. Physical Review E, 74(6):066704, 2006. 59
[145] Yunkai Zhou, Yousef Saad, Murilo L. Tiago, and James R Chelikowsky. Self-consistent-fiele caleulations using Chebyshev-filteee subspace iteration. Journal of Computational Physics, 2 2006. 59
18

---

A Preliminaries
In this section the model of computation is defined and basic linear algebra principles are summa- rized. Further details can be found in standard textbooks such as [70, 64].
A.1 Forward and backward approximation
In numerical analysis, a "forward-error" often measures the distance between the true solution to problem and the solution returned by a numerical algorithm. related notion is the so-called "back- ward error" or "backward approximation." For more details we can refer to the standard textbook of Higham [70]. In this case, the solution returned the exact solution of "nearby problem," " and the backward error quantifies the distance of the original problem to this nearby problem. In PCA, for example, backward- aawraaadpppooiittioonnn could be defined by C arg T tr(C' НC).
The backward error would be given || H| for some norm. An algorithm is backward stable if the backward error is always well-defined and bounded. A forward error type approximation often harder to achieve. A common rule-of-thumb states that
forward error backward emror condition number of the problem.
However, this does not generally hold for eigenvalue problems. For example, Proposition 1.1 in [15] can be used to transformed a backward approximate diagonalization error to a forward error for the eigenvectors, but the bound depends on the minimum eigenvalue gap between any eigenvalue pair. As the authors point in Remark 1.2, special treatment is needed in terms of the invariant subspaces in the presence of multiple eigenvalues, which are thoroughly analyzed in this work.
A.2 Floating point model
We assume standard floating point model of computation and borrow its axioms from [70, Chapter 2]. There is fixed number of bits to represent floating point numbers, specifically, one bit reserved for the sign s, bits used for the exponent e, and & bits are used for the significand m.
A real number a E I is rounded floating point number
fl(a) s 2 m.
The sign s is + if the corresponding bit is one, and the bit is zero. The exponent is stored as a binary number in the so-called biased form, and its range is e E [-M,M, M). where M = 2n-1, The significand mn is an integer that satisfies 2 S m S 2 1, where the lower bound is enforced to ensure that the system is normalized, i.e. the first bit of is always 1. We can therefore write fl((x) more intuitive representation
((aa) -2
122
where the first bit 1 of m is always equal to one for normalized numbers. The range of normalized numbers is therefore [2-M 2 (2-2-1)) Numbers that smaller than 2 and they will be ignored for simplicity, since we can either add more bits in the exponent, or account for them in the failure probability when the numbers random (the latter strategy is used for example in Lemma D.I). Similarly, numbers that are larger than 2 M (2 - 2- are assumed to be numerically equal infinity, denoted by INF.
called subnormal
From [70, Theorem 2.2], for real numbers a the normalized range holds that
((()) (1 + 0)a,
where 0 E IR satisfies |0| 2- u. where is the machine precision. Clearly, f Ollgg I/u)), which case we can always obtain a bound for the number of required bits a numerical algorithm if we have upper bound for the precision u. We will write the same for complex numbers which are represented as a pair of normalized floating point numbers.
The floating point implementation of each arithmetic operation @ E [+, x./ also satisfies
fl(a + B) = (1 + (a ② 3),
|0| < u.
(4)
19

---

Divisions and multiplications with 1 and 2 do not introduce errors (for the latter we simply in- crease/decrease the exponent). We assume that we also have an implementation of v such that fl(va) = (1 + 0lya where [ s u. From [70, Lemma 3.1] we can bound products of errors as
II + (1) 1++10 Th.
where P: = +1 and |7 1 I-nu
The above can be extended also for complex arithmetic (see [70, Lemma 3.5J), where the bound becomes [0] < Ou)) but will ignore the constant prefactor for simplicity.
Operations on matrices can be analyzed in a similar manner. Let @ denote the element-wise multi- plication between two matrices and the element-wise division. The floating point representation of a matrix A satisfies
(I(A) = A AA
S u.
It can be shown that |||| S uu/n|AA
For any operation {+, -, 0, 2 and matrices A and holds that
fl A B) A O B
3 B), A. u, A & (A o e S
Suvn & BlJ. (5)
A.3 Spectral decomposition, pseudospectrum, and eigenvalue bounds
We first recall the definition of the spectral decomposition of diagonalizable matrix.
Definition A.1 (Spectrum and spectral decomposition). A matrix A is diagonalizable if there exist invertüble matrix V and diagonal matrix A such that A VAV-1 This is called the spectral decomposition of A. The set A(A) {A,,li 1, n.) is the spectrum of
The spectral theorem states that Hermitian matrices (or, generally, normal matrices) can be always diagonalized by unitary transformations.
Theorem A.1 (Spectral theorem). A E C is Hermitian, then there exists orthogonal matrix Q C such that A- QAQ", where A diagonal matrix with real diagonal elements.
In the sections that follow we need to bound the (forward) errors on the eigenvalues of matrices under perturbations that introduced due to the finite precision arithmetic. Such bounds can be derived by the classic Bauer-Fike theorem the following (cf. [24, 48] for more details).
There many bounds in the literature describing the effect of perturbations the eigenvalues. We summarize some classic results in the following proposition and refer to Bhatia's 's monograph for detailed overview [28].
Fact A.1. Let A and B be two eigenvalues (A) I (2) < < (A) of and (B < ... & An (B) of B:
matrices. The following bounds are known between the
bound
reference [139]
Hermitian
Hermitian
[A.(A)
(B S B
Hermitian
Non-Hermitian
[A((A) A. - (llo(())] Bl|
[81, 117]
Pseudospectral analysis is useful when the aforementioned bounds are not applicable.
Definition A.2 (Pseudospectrum), For some c >> the c-pseudospectrum of a matrix M Cox
is defined as:
(M) = (AEC: A(M E) for some El| <
=:EC: [(AI M) 1/e}
where A(M) is the spectrum of M.
Recall some useful properties from the seminal book of Trefethen and Embree [135].
20

---

Proposition A.1 (Collective results from [135]), Let D(2,r) denote the open disk of radius in the complex plane centered at 2 C. M and E be two n. X n matrices, and E be a positive real number. The following hold:
(i) A, (-]| (M) S (M C ([EE) (M) ([135, Thm 52.4]).
(ii) Any bounded connected component of A, (M) has a nonempty intersection with A(M,, i... contains least one eigenvalue (I135, Thm 4.3]),
(ii) U, c)), S Ac(M) S U, D(a, ckty (M) ([135, Thmss 4.3 and 52.2)).
A.4 Eigenvector condition number of definite pencils
Some approximation bounds throughout the paper depend on the eigenvector condition number of the generalized eigenproblem (1). For arbitrary GEPs a bound for this quantity might not always exist, but a straightforward bound exists for the Hermitian definite case. obtain such a bound we define the eigenvector condition number of diagonalizable matrix as follows:
Kv (A) ==
inf
|VI|||
VDV-I=A D: diagonal
Proposition A.2. Let H E IHI and S H and consider the definite GEP HC SCA, Then
kv S (K(S).
Proof. Since Hermitian positive-definite, can be written as LL" for some matrix L. Then we can transform the GEP to Hermitian eigenproblem, specifically
L HL 'L'C = L'CA.
Since LL" S we have that |[]] LL'I| || = |||| and similarly ||- = |||- 1| which means that <L) = /S(). Since L 'LL L-HLL- is Hermitian, it can be diagonalized by unitary matrix, i.e. and L diagonalizes 'HL" For Ö we
there exists C such that [ Cll = [[(=
have that
ट = L L < L L L' = ||
Similarly,
11111110010101100000001111111111 |L*||
Undoing the transformation, we can see that Ò also satisfies:
HÒ = SČA,
ie. Ò diagonalizes ' since C Č's s-'s-ppe 'HC = A. We conclude that (y I x(O) <
[|||
A.5 Imported floating point algorithms for fast linear algebra
Before dive into the details, recall the concept of "fast matrix multiplication", pioneered by Strassen [134], who showed that two square matrices can be multiplied O(n"), where w = log, 2.807 3, in real arithmetic. Since then, the matrix multiplication exponent w has been significantly reduced, the record to date being w < 2.371552 [47, 141]. In two seminal works [40, 39], Demmel, Dumitriu, Holtz, and Kleinberg, proved that any fast matrix multiplication algorithm can be executed numerically stably in floating point machine with almost the same arithmetic complexity as in real arithmetic. They also showed that other problems in Numerical Linear Algebra can be reduced to such matrix multiplications, including inversion, LU and QR factorizations, and solving linear systems of equations. Our algorithms build on the existing results for stable fast matrix multiplication [40], inversion [39], as well as backward-approximatt diagonalization and the matrix sign function [15]. We import the corresponding results from the aforementioned works in the following theorems.
21

---

Theorem A.2 (MM, stable fast matrix multiplication [40]). For every 1) > 0, there exists fast ma- trix multiplication algorithm MM which takes as input two matrices A, B C* and a machine precision u >> O and returns C <- MM(A, A, B) such that
|CC ABII S (MM(N))
u||
on floating point machine with precision u, where (MM(()) n,",, for constant C, independent of n. Such an algorithm is called (MM(() (n)--table, and it requires TMM(n) O(n arithmetic operations, where is the exponent of matrix multiplication in real arithmetic.
Theorem A.3 (INV, logarithmically-tahe
fast inversion [39]). For every > 0, there exists
fast inversion algorithm INV which takes as input an invertible matrix E C and a machine precision u > O and returns C <- INV(A) such that
IC-A (TN u- R((A(A)) Civy log(n) || ||
on floating point machine with precision u, where (Nnv(n) = Onn for some constants
ANV S 8 and Cn independent of n. Such algorithm is called (MNN((n), CNV -stable and requires Tiwv(n) = (TTNma(() arithmetic operations, where (Mm(()) is the same as in Theorem A.2.
Sign function and deflation. The next result that import is floating point algorithm for the matrix sign function, which, for diagonalizable matrix A VAV such that A(A) does not interesect with the imaginary axis, is defined as sgn(A) (A) = V sgn(A) (A)V- where sgn(A) A) is diagonal matrix and each diagonal entry contains the sign of the real part of the corresponding eigenvalue. To state the main result, we first recall the definition of the Circles of Apollonius.
Definition A.3 (Circles of Apollonius, imported from Section 4.1 of [15]), Let E (0, The Circles of Apollonius C. are defined as C., = C+ U C..
C = C: |m(()] S sa), =
a).
where m(z) =
1-s is the Möbius transformation taking the right half-plane to the unit disk. The disk 1+2
C is centered at 1+a² and has radius 1-a2. and Cā is its reflection with respect to the imaginary axis.
2a
Theorem A.4 (SGN, imported Theorem 4.9 from [15]). There is deterministic algorithm SGN(A, a, n, e) which takes as imput a matrix A E C"x", an accuracy parameter E (0. 1011, and parameters (0, 1), 0 < 1 1/100, such that it is guaranteed that A, (A) Ca- The algorithm returns a matrix such that
S sgn(()) Se
as long as the machine precision satisfies
a N (enov log(n))
USGN
(INv(N)
corresponding to most
(log(n)log) log (1-2) (log(?) + +log(t))
bits of precision. Here
- [log( [log(1-a) 3log(log( + log(log( 1) + 7.79]
denotes the number of iterations that the algorithm executes. The arithmetic complexity is
O(NTMM()))
In PCA, we will need to be able to compute a basis for the column space of an approximate spectral projector. For this, we recall the following algorithm.
22

---

Theorem A.5 (DEFLATE, imported Theorem 5,3 of [15]). There exists randomized algorithm C₁ + DE k 3, which takes as input a matrix II. rank parameter k [n], parameter B such that |[I II|I I for some projector matrix II of rank-k, and a desired accuracy e, and returns a (complex) matrix C₁ such that, there exists a matrix C₁ with k columns that forms an orthonormal basis for im ( and
C C.
The algorithm succeeds with probability at least 1 (20m) The number of arithmetic operations is most (T)))) and it requires (Ilog(n) {oog(})) bits of precision. Internally, the algorithm generates (n²) random numbers to form n x complex Ginibre matrix.
Computing an orthonormal basis and the spectral norm. Lastly, we recall the following two results. The first one is for the QR factorization from [39], which adapt suitably for our analysis, specifically, to compute an orthonormal basis for the column space rectangular matrix.
Theorem A.6 (Basis computation, follows from Section 4 in [39]). Let A E Rmxx m n. There exists an algorithm Q R = QR(A), which returns a matrix Q E Rm and an upper triangular matrix R _x in O(mn²=1) floating point operations using Olog( mx(A) )) bits, where CQR (0, is given accuracy. The matrix Q satisfies
|1Q Q S CQR,
for some orthogonal matrix & E Rax", where Q has orthonormal columns and A =QR is the true economy- QR of A.
Proof. We first scale A' + A/M, where M is the smallest power of two that is larger than nllA max- This ensures that ((1/n) S |||l|| < 1. Since we scale by a power 2, then there are no floating point errors and the orthonormal basis from the QR factorization remains the same (only the upper triangular factor scaled).
We now use the corresponding algorithm of [39] on A'. returrss three matrices: an upper triangular matrix R R=*x matrix W Rmxn with |[ E O(n), and matrix Y Raxm with || Onn). The matrix L I WY + Eg exactly satisfies D T I, and the enror matrix
satisfies |E || (ppolynn)u), Moreover, L
R
= À, where Ā- A' + EA and |EAl <
(0-nXx
((poly(n)u) A| S O(poly(n)u). Note that |R|| S
|| S ||'I EA S 1 + (ppoyy
The total cost of the algorithm is 0(mm"-1) floating point operations (ignoring the negligible term r > in the exponent).
Using W and Y, we can construct an approximate basis Q as follows. Let Y, E Ra contain the MM(W, Y, in O(mm(-1)
first n columns of Y. We can compute the matrix Q I
(m-nxx
floating point operations by performing the multiplication in blocks of size n x n. We can write Q Q Eg where Q contains the first n columns of &p and || E (p)oly)n)u).
Note that R satisfies QR = A' + Ex. Since Q has orthonormal columns, the singular values of R satisfy o (R) = (QQ) = ai + EA) E (A') + ||EAll where the last comes from the stability of singular values which is consequence of Weyl's inequality in Fact A.1. Now, if ||All < 1 min (A'), which is achieved by setting u S (Pol((n)(AA) for some constant c, then R is full rank, and holds that | || 1/amin (R)-
3o (A') & <O(n){(A).
Since R remains full rank, we have that QR = À Q = A'R + EAR If we write the true economy-QR QR of A' as A' = QR, then we can see that |RR = |QRR = |
We highlight that our Gaussian sampling oracle in Definition D.I is slightly different than the one described in [15]. but it should be straightforward to adapt the analysis using either of the two definitions.
23

---

Q RAR J. This implies that R.R. 11 E | 11,1 + |EAR II].
With similar
arguments obtain
|RR 11 [1 E R 11, EARRIII.
(())n(A). Thus, if we fur- We previously argued that |R]] R S O(1), and R | 11 I
ther enforce |Eal S 61 poly(a)<(A} for some constant C1. then max BAR 1l. EAR EAR
(,, for some E (0, which means that all the singular values of R.R are inside the inter- val [1 - 2,,, + 2e1]. In other words, RR is approximately orthogonal. Indeed, we can write RR = Φ Egg where $ T = I and all the singular values of Eg. are inside [-26], 261] In particular, = UV where U, V come from the SVD R... -
UEV.
We can now go back to Q. From the above, Q can be written as
Q AR
EAR
Q + QE EAR
which ultimately gives
Q + Q = QΦ + Q + EAR + Eq.
This means that Q is just rotation of the true Q plus some additive error terms that we can control with u. In particular,
I Q
$ QI + EAR + E +|+ (poly(nn))
e + O(poly(n))) < ((e)),
where in the last we applied the required bound for u v 61 poly(a)a(A) Combining everything, we conclude that if the machine precision satisfies
S CQ (oly())x(A)
then we compute a matrix Q in 0(mm" floating point operations such that there exists orthogonal matrix T satisfying
I - QΦ S EQR:
using ((0o((m)(A)/eqq)) bits of precision.
The second result is for the computation of the spectral norm of the matrix from [107] using the Lanczos algorithm.
Theorem A.7 (Imported variant of Theorem 18 from the full version [107]. Let A Rm and 5 E (0, 1/2) a failure probability parameter, We can compute vector such that, with probability at least 8.
|
S
AAyl TyT
|||||
TO
in O (mn log(n) log(1/6)) floating point operations using O(log(n)) bits of precision. Internally the algorithm generates On log(1/6)) random bits.
We can obtain the following corollary.
Corollary A.1 (Spectral norm), Let B E Cm We can compute value E E O[[ with probability at least 1 - 20 in O (mn log(n) log(1/8)) floating point operations using O(log(n) bits of precision.
Proof. Let R. be the real part of B and be the imaginary part. We can approximate R E 9 R||
1o R Z 0.9 R|
[|] and Z
10.
9
Z
[Z]]] using the algorithm of Theorem A.7. We then set Σ [Z]]] Using the triangle inequality we have that
Z
R|
|
B
< 0.9(
R| + [|]]) < S < R + ZZ| < B
or, in other words, E E [B
B]]]
24

---

A.6 Symmetrization
There are certain cases where is crucial to ensure that the result floating point multiplication inversion remains Hermitian. The following proposition states that always "symmetrize" a floating point matrix with small additional errors.
Proposition A.3 (HERM). Let C be a Hermitian matrix and C' C+E for some (non- ) matrix E, such that ||C C'Il = ||B|||- Consider the matrix HERM(C"), where HERM( A) is algorithm that replaces the strictly lower triangular part of A with the strictly upper triangular part of A. Then (ER)(C") Hermitian and it holds that:
||C HERM(C (C)) S CHERM log(n) ||||
for some constant CHERM-
For any matrix A, the upper triangular part of A, denoted as Au (A(A), satisfies the following inequality:
Proof.
| (A)ll S clog(n) ||A||,
for some constant c. The same holds for the lower triangular part A (A). For the diagonal part
diag(A ()) ) it holds that diagg(A A]l
|| (see [27] for proofs). Let C" HERM(C'). By the
definition of HERM(-). the matrix C" be written as C" (C(C") AC')" diag(C')
C (Ar (E) + Au (E)* diag E)). Therefore
Au (E) + AU (E)* diag(E)]l S 2clog(n)
E + El| S CHERM Iog(n)
[B]]
for some constant CHERM
We can directly this to derive bounds for matrix multiplication where the result is Hermitian matrix and for Hermitian matrix inversion.
Corollary A.2 (Symmetrized matrix multiplication). Let A and B be two matrices and C = AB be a Hermitian matrix. If MM is a ptm-stable matrix multiplication algorithm then the matrix
C' = HERM(MM(A, B))
is Hermitian and it satisfies
||C' Cll CHERM log(n) u| | B|| (MM
Proof. Straightforward combination of Theorem A.2 and Proposition A.3.
Corollary A.3 (Symmetrized inversion). Let A be an invertible Hermitian matrix. If INV (UNN),(iv)-stabble inversion algorithm then the matrix
C HERM(INV(A))
is Hermitian and it satisfies
|C A -- S CHERM log(n) MIN((n) u x(A) JNv 1i A A-1
In addition, if u satisfies
S
CHERM log(n). (INv(n)) k( A caN log(n)+I
for some E (0,1/2), then both of the following hold:
|C - -C-A--10 < <|IA -111.
x(A) <c(C) & 2x(A).
2
25

---

Proof. The first part is straightforward combination of Theorem A.3 and Proposition A.3. For the second part, we directly bound
( (1-+)l (1
I|A 1 - ||||| < (1 + R(A) 11 s (1+ + "1l
R[A]
A 1 | ||||| - | '11.
Note that |C--A A S S<AA A 'Il = T
1 E0 min (A Since C and A are both
Hermitian then from Weyl's inequality
e)o (A A min (C) S (1 e)) (A 1-
7 m (A (A-1) S min (C) S 20 (A 11)
where we conclude that
}<((A) S (C) & 4k(A).
Remark A.1. Throughout this paper, to simplify our bounds, we assume that
1 15 pMM((T) (1), CNV log(n), CHERM log(n).
Moreover, we define a global upper bound j such that
Huna(n), CHERM log(n) (n), (uNv (n) S p(n) S O(n%)).
where, from Theorem A.2, cn a constant that does not depend
A.7 Existing algorithms for invariant subspaces and the (generalized) eigenproblem
We give a brief overview of some results related to the computation of spectral projectors, eigen- vectors, and invariant subspaces for the Hermitian definite generalized eigenproblem. Many more details can be found in standard textbooks and references therein [42, 119, 64, 135, 115]. as well the recent overview of [129] for the complexity of the eigenproblem. For further reading can highlight some influential works from the enormous bibliography on eigenvalue algorithms and perturbation bounds [44, 43, 8, 130, 131, 132, 133, 75, 108, 9, 26, 76, 77, 37, 2, 113, 112).
The prominent classic method to solve the Hermitian (generalized or regular) eigenproblem is to use numerically stable Householder transformations to reduce Hermitian matrix to tridiagonal form, and then apply a tridiagonal shifted QR algorithm to diagonalize the tridiagonal matrix, which gives both the eigenvalues and the eigenvectors. The paramount QR algorithm was originally proposed in [51, 52, 90]. its convergence for the symmetric case was analyzed in [140, 38] in exact arithmetic, and the non-symmetric was recently analyzed in floating point in [18, 16, 17]. Other classical algorithms for computing invariant subspaces include [97, 98, 10, 9].
With the aforementioned procedure, the floating point complexity is about Ö(r²) floating point operations, potentially up to polylog(1/e, 1/ gap;) factors. It is known that for Hermitian definite pencils (H, S) one can extend this procedure such that, simultantcouslyy H will become (Hermitian) tridiagonal and S the identity [36]. The Hermitian tridiagonal matrix can be further reduced with similarity transformations to a symmetric tridiagonal matrix, in which case one can apply again tridiagonal symmetric QR. Using standard perturbation bounds the ones used in this work one can potentially obtain provable forward errors invariant subspaces.
While there exist works that overcome the Ö(n²) worst-case complexity barrier for some compu- tations related to the cigenproblem, none of them provide end-to-end forward errors for individual invariant subspaces. 14] showed that the eigenvalues of matrix can be computed in O(z)) arith- metic operations, however, it becomes O(n*+1) boolean operations in rational arithmetic (up to
some other omitted factors), and therefore slower than standard Ö(n²) floating point eigensolvers. [25] proposed a quantun-inspired method for approximate diagonalization in backward-error sense, with (n bit complexity. [15] proved that hapkward- diagonalization can be solved in Ö(n²) bit complexity in floating point by using smoothed analysis, improving both the Ō(n²) classic algorithms for the Hermitian case and the (ppoy(n//e) algorithm of [7] for
26

---

the non-Hermitian case. This was extended for the generalized eigenproblem case in [41]. As already mentioned in the introduction, obtain forward errors for eigenvectors from the backward- approximate solution, the corresponding bounds that are reported in [15, 41] require the existence of minimum eigenvalue gap in the original matrix, while the analysis for invariant subspaces was left as open problem. Our work completes this analysis, and takes step even further, by show- ing that (at least based on existing results), explicit diagonalization is redundant for approximating individual invariant subspaces.
B Proof of Proposition 2.1
PURIFY
Input: Hermitian definite pencil H E HI" S E + approximate condition number ñ of S, approximate
values Pe. gapg. accuracy (0, 1/12).
Requires: ||]|| < 1, |SS' S 1, R ((x(S)), jaa E μa + gapu. and gap, E (1 + ( gap,. where A A and gap; = A 1 for some k E [n 1). Algorithm: P + PURIFY(H, S, Paa gapg: c).
1: SINV < ISV(S))
b SINV = S-1
1--1- E
2: H MM(SINv. H).
H #SI=Ii SINv
EM
3: M jiu H E
4: C 5GN
M, 16-199 5: II - (1- + E
SaP
C sgn(M) + ESC
6: return II.
Output: Approximate spectral projector IT.
Ensures:
[T II| S.
Algorithm PURIFY.
For the analysis of Algorithm 3 first recall that, in exact arithmetic, small perturbations do not have a severe influence in the sign function of a matrix.
Lemma B.1. Let H E H S E H where is a scalar, and M = - S--S and M (Amis(Mii for some E (0,1), where Mmii((M is the smallest 128kx (M)
- S S E. If [ S
eigenvalue of M in absolute value, then [M) sgn(M)il S.
Proof. The proof uses standard techniques from holomorphic functional calculus and the properties of the pseudospectrum, but it is stated for completeness. A similar proof, for example, can be found in [9] for general matrices A (i.e. in that proof A can have complex eigenvalues).
Recall that for a matrix that has no eigenvalues on the imaginary axis, holds that
sgn(A) (A) P+(A) P-(A),
where P- (-(A) is the spectral projector on the subspace spanned by the eigenvectors corresponding eigenvalues with negative real part, and P+(A) is the corresponding spectral projector to the positive halfplane. Then
sggnnMM M) sgn(M) (M)l IP+ (M) + (M) P+ + E) P_ + E)]l
< IP+ (M) - P+ (M + E) + |P_(M) P P-(M + E)]l
S- H has real eigenvalues because it is similar to a Hermitian matrix, i.e. A(S H = A(L'HL, where L the lower triangular Cholesky factor of S S-1. Thus M has also real eigenvalues. From Proposition A.1 (i) and (iii)
A Tij M) S A2 BI (M) S
UD
A (M), EM)
27

---

If |B| S E1
Aratn (M)
M for some E (0, 1) then Amiin (M) 2 Amin (M)) and Xmmax (M)
Mm [M) Aii (M)//2.
Let R be rectangle whose bottom-left comner is located at -i] Amin (M) and its top-right corner at 2 + Amin (M)) + / (M)). Any point the boundary of R has distance of at least Amina (M)
from ((M) and at least at least [Amin (M))/2 from A(M). From holomorphic functional calculus
we have that
L (3-M)-
P+(A):
JoR
where in this case 0R denotes positively oriented rectifiable curve over the boundary of the rect- angle. Using the resolvent identity (z (M) (2 (2- = (z (z-MM) M) (M- M') M for P+ we have
(2 M) d2 1 $ (2 M + E)
(P+(MM P+ M E)
dz
2Ti 8R
2Ti JaR
(2 - M) -1 - (2 M + E) da
2Ti 0R
$ (2 M) -1 (2 M (z M E))(2 M+ E) dz
2Ti 0R
(2 M) -1 EE(2 M + E) dz
2Ti aa
1
[1 M)- M E) 10-1 dz
2T
JaR
1
length(At)
|
sup =EaR
[1 M)
sup I1(2 - M E1 =EδR
2m
1
Amin (M)
2
(6)
Amin (M) + 4)
2T
4Ny (M)
Amiin (M)
Dmmin [M)
4
é1
TKV (M) Amin (M)
Since (M)) then we set 61 S
Amin(M)] 16
(which satisfies 61 < 1), where is the
desired final accuracy, this gives
(M) P+ + E) S 5/22.
The same can be obtained for P_. Putting everything together, EE|
Pmmin (M) for some 128xV (M)
E (0, 1) ensures that ||sgn(M EE) -ggl(M)] s c. From Proposition A.2 (M S /S)), which gives the final bound.
We can now prove Proposition 2.1, which we restate for readability.
Proposition B.1 (Restatement Proposition 2.1). Let H E H with |H| < 1, H with [n 1] and
1). Let = duthe+1 and gapı = A ARt1. where A1 S... Xn are the generalized eigenvalues of the Hermitian definite pencil (H, S) and assume that we want to compute II which is the true spectral projector associated with the i smallest eigenvalues. we have access to
|S I
Pa E Pa gapp gap, E (1 gaPEs
R (n(S),Cn(SS)
for some constant C then Algorithm 3 computes I1 + PURIFY S, Jie, gap Rú E) such
that [ n&ll S 6, in (1) (log( gape + log(log( ( 11) floating point operations
using (log(n) log"( Baps log( (()) bits of precision.
28

---

Proof. Using the notation of Algorithm 3, we have that H = S S initially set
INV
H
+
Ez
Let us
11
u(n) ((ln)nnwv (1g
u(n)x(S)aw logtn]"
(6)
for some eo E 1/8) (to be determined later). From Theorems A.3 and A.2, and from the assump- tion that |[]] ||- S 1 we have that
|NV| || - (
1og(n) || ||- S 60.
E2 MM
S (
H [ SiNv|l & wpmm (1z) H (IS- '1 [E.'/]) 2 co(1 + 60) S 2Eg
We next have M iu H E
F S H E H E EAM + EnJ juk -5 -S'H 'H+B.
where Eg is diagonal error matrix with
E3 < u
H|| < u
INV H
E
EM
< u( + co + 2e0) S 4u < 4ea,
and B = E H EM + E We can bound the norm of B as:
[B S TEp.
Now we need to apply Lemma B.1 to argue that |ssgn(u - SS sgn(i - S S <
é1, for some E (0, 1). To satisfy the requirements of the lemma, we need to ensure that B 5
n A (6p- S- S-'H) This can be achieved by setting
61 128 x(S)
T T| min (TN S S-H
60
7. 128 R(S)
/R(S)
However, we do not know Amin (āx S S'H
H)| and <(S). To circumvent this, we use the
assumption that üig is well-placed close to the center of gap which implies that
Amin (it-
H))
2010 gapk
gapk:
and also use the assumption that k{S) $ R. This means that we can set
60 7. 128V K
gaps
which implies the desired bound for B:
T122
min (Mu S S-HH)
| S 7 7 128
gapt 61
128 ((())
We can now apply Lemma B.1 to argue that
[sg S H) sgn( - S B) S1.
By assumption we know that 74u = ph gap, which means that sgn (itx - -S-' = =sgn(4)
S- S-'H), and therefore
||sgu(uu- -S -S - S B) S11-
To not interrupt the flow, we will skip for now the analysis of SGN and we will leave it for the end of the proof. From Theorem A.4 we know that the result returned by SGN satisfies ||C sgn(M) <
<SGN, where we have full control over <sGN since is passed as an argument to SGN, and therefore
it can be set to ESGN = ej. Denoting C = sgn(us - S -S-'H) and Ec = C - C. we have
| sgn(st S S-'H) = |Eel S + ESGN = 2
29

---

We now proceed the line of Algorithm 3, which computes II
1 2 (1+
+ E where
Eg is diagonal error matrix with norm bounded by
B || - ul| u||C|| - u(1 + 61 + eSGN) S 3u < 361.
Then we can write IT = (1 + + Ec) + E +) = + Ec 2 + E (+)
1I ETI, where
II [1 + C) and Eµ 2 Ec E Combining with the above:
|En I 461,
in which case we can set 61 4. where is the desired accuracy, guarantee forward error of most c for the spectral projector II.
Gathering all the requirements for the machine precision (except for SGN, which detailed below) from Equations (6) and (7), and the assumption on K and gap; it suffices to set
47 128VV
ggapk
p(n (n))oowwwWe log(n)
which translates
(log(1/u)) = (log) (kjuln)zsw/bate)) bante)
)) =0 (1o((n)
log (ggpp,
(8)
rgaps
bits of precision.
We finally proceed with the analysis SGN. use Theorem A.4 for the convergence of SGN, we need to find the appropriate parameters OSGN, 1/SGN call it. These parameters must be such that 100 OSGN 1. From the properties of the pseudospectrum, specifically,
Aysa (M) C Coscy- and 99
from Proposition A. (i), we know that for any > 0
(M) An (īR S
H B) C A (i S S S Aq (in S --'
where Buax is the upper bound for |B] that we obtained above:
Tgapk
B S
Binax
128V R
For 17SGN = Bmax and from Proposition A.I (iii) we obtain
Byac (M) C A2Bmar (Tt S S-'H H
SUD (it A, 2B 2B By (īq S --)
C L D ju - H- 221
mgap2
9-128
where we used the fact that Rv (S
H S (S)
We can get a rough upper bound
ngapt 2e1
S H 6l
gr
((S)
gaps S
gap T28
128VR
9-128.8' v<(S)
which gives
A Brase M C
U (īo
A H), El
gap
suffices to find the appropriate asGN such that the Apollonian circles Casax will contain the disks of the pseudospectrum. Since the smallest singular value of jág - S is lower bounded by gap,, the lefumost positive point where any pseudospectral disk intersects with the real axis is
[ gap; -1 gap2
> ( 1zs
gapk 128
gaPa Similarly for the negative halfplane, the
rightmost negative point where disk intersects with the real axis is gapa Denoting by Ç gap.
30

---

it suffices to set OSGN 1 8-gape
1-5 s-gape Clearly, OSSNN < 1 for all x gaps > O. For the lower
bound, to ensure OSGN > 99/100 we need that gap; <
8
To achieve this we can simply scale
199
the imput matrices by constant, since, by assumption, ||S H| || conclude, we have the required parameters call SGN.
1, and therefore gapı _ 2. To
We can now go back to the complexity analysis. We are running
Č f SGN M aSeN, 7SGN, ESGN).
where
8 gapp
Tgapa
aSSN
TSGN -max E1
ESGN E1
8 + gapp
9 128
4
We now invoke Theorem A.4. The number of iterations is bounded by
N = 0 (log( (asso))) + log(log( VISNNESON
O [log( (z++,) + log (x
O [og b + log(log( x(S)
gap, ))
The number of required bits is
O (log(n) log (T log( (manesso))
0 (log(n) log
1111 log
gape
cgaps
o (log(n) (2) log'(
log
x(S)
gape
(gap,,
which dominates the previous bound of Equation
It remains to bound the arithmetic complexity. There is constant number of inversions, matrix multiplications, and scalar-mate matrix operations. which total take most (1)))(nz) floating point operations. Thus, the dominant factor is the call to SGN, which amounts to
O (TMM(()) (log( gAp, log(log(, ((S) 1))
arithetic operations using the aforementioned number of bits.
B.1 Spectral gaps with diagonalization
The straightforward approach to compute the desired spectral gaps is to iteratively compute S-'H and diagonalize it, until the eigenvalues are well approximated. If we compute A <- MM(INV(S)) H), then we can write A = S S + E, for some eπror matrix E. The next step is to approximate (S- and (s For this one could use the recent state-of-the-art approxattaprooiimate diagonalization algorithm of [15].
Theorem B.1 (EIG, imported Theorem 1.6 from [15]). There exists randomized algorithm EIG(A, EEIG) which takes any matrix A E C with ||A]] < 1 and a desired accuracy param- eter EEIG inputs and returns a diagonal D and an invertible matrix V such that
A VDV <EIG and Kl x(V) < 32n 2.5 /EEIG,
in
o (M) ( arithmetic operations on a floating point machine with
Ollog* (n/gg((n/eee (n/ceI))
bits of precision, with probability at least 1 14/n.
31

---

Applying EIG on A we can obtain a bocckwrd-approxinnne diagonallization. But we are not fin- ished yet, since we are interested in each individual eigenvalue. To translate the backward error to a forward error for the eigenvalues, and, ultimately, the spectral gap, one can try to use Corollary 1.7 and Proposition 1.1 of [15]. However, this approach has two main limitations. First, relies on simplicity of the spectrum, i.e., it assumes that the minimum gap between any pair of eigenvalues is larger than zero. This assumption is quite restrictive, since the desired gap might be well-defined even at the presence of other multiple eigenvalues. For example, in DFT applications it is not un- common to have eigenenergies with algebraic multiplicity larger than one, and at the same time have large band-gap gan that separates the occupied from the unoccupied orbitals. The second limitation is that the aforementioned Corollary 1.7 requires as an imput parameter an actual over-estimate for the minimum eigenvalue gap. Even if such gap exists, it is not described how to estimate it.
For diagonalizable matrices, we can leverage the following Corollary B.1, which an immediate consequence of Kahan's ss inequality (Fact A.1), and it overcomes the aforementioned limitations.
Corollary B.1. If X is diagonalizable and it has real eigenvalues then for any Z the following bound holds for the eigenvalues of X:
(xx(X) (Z) S O(log(n))xv(X)1x> -Z]l-
WAW-1 where W diagonalizes X and is chosen such that K s(W) xv(X) This is always possible since if W is any matrix that diagonalizes X then A is similar to X and since X has real eigenvalues then has to be real, and therefore symmetric. Then
Proof. Write X
write
A W Zzw| & || w |||X
Z
K [[] ZII-
Since A is symmetric, and since W ZW is similar to Z then from Kahan's inequality (Fact A.I)
(A(X) - (z) = (A) A (W 'ZW)
s Ollog(n))lIA w ZW]
s (()))( - Z|
((log()))xv((IIII ZI|-
Now we need to use this to get a bound for the computed generalized eigenvalues after INV and EIG. Proposition B.2 (EIG-based gap). Given a definite pencil (H,S) with |H],, ||- < 1, we can compute gap; (1 e) gapı and ju pa te gapı, by iteratively calling INV, MM, and EIG, using
o (T) log( 1
log2?
pm((S) egap,
floating point operations using
O (log' 4 (x((8) log(n)
bits of precision with probability at least 1 - 0( log l/e gaps) ).
Proof. Let D, V EIG(A, CEIG) be the solution returned by EIG when applied A -S -S-H 'H+E E.
Note that
||S-'H VDV |S S H E VDV- |S- + E VDV- E S <EIG + |E|
From Corollary 1. since S HH is diagonalizable with real cigenvalues then we conclude that
|A:(S-HH) D.2l |A (VDV=1) < ((0og(n)xxx(5~11 S- H VDV-11
O(lg((n) (oog(nn)nv(S-11)) (EEIG + |[B|])-
for some chosen e', and We can tune the machine precision such that ||||| = EEIG = clog(n) va(S)
some global constant then finally
( ( - \:VVVV SE.
32

---

We can now consider an iterative scheme, where we call EIG on S S''H, H, halving e' at each step. We need keep halving e' until it reaches e' ( gap,). in which case we have total of Ollog( ( calls to EIG. In the worst case every call costs
o (Thaa( o log ((e)))
=0 (TTham(()) log (N(S))
gap, log(n)) bits (note that the bits required by MM and INV to achieve accuracy e' are dominated by those required for EIG, and we therefore ignore them). Since we do not know gapk, we can set the termination criterion to be c' = ((ggap)). where gapı the approximate gap that we obtain from EIG.
arithmetic operations using log ( (
Each iteration succeeds with high probability 1 - 1/n, in which case a union bound gives 1 Ollog(1/(e gapx))/n, which can potentially be improved complexity, but we do not expand further.
O(1/n) without impacting the
A subtle detail in the analysis above is that need an estimate for k(S) in order to be able to INV and to set the machine precision. Since x(S) is generally unknown, we need to compute it. We will later show in Appendix E how to compute the condition number quickly.
When S I, we can avoid the inversion and the computation of ((S). The arithmetic complexity and the bit requirement are the same by setting k(S) 1.
O
Remark B.1. We highlight that a similar result can be obtained by using the more recent pencil diagonalization algorithm of [41]. The authors the latter mention that their algorithm should be "more numerically stable" than EIG of [15], since it uses an inverse-free iteration internally, and they do provide strong theoretical and experimental evidence for this statement. However, there is no formal, end-to-end proof for the bit complexity at the time of this writing, we choose to com- pare against EIG of [15]. Note that, any improvements on the bit requirements for the matrix sign function, directly provide the same improvements for our main Theorems 1.1 and 3.1.
C
Analysis of Cholesky
This section devoted the analysis of Algorithm 2 and the proof of Theorem 1.2. Let M
A
B"
be Hermitian matrix. It can be factorized in the form
B
C
BA-1 () s) ( A A )
I
) WYW",
M =
(9)
I
where S = C -B- -B-BBB B" is the Schur complement and
w BA-1 -A-+1) 1). Y=
A
$ 8)
A
Proposition C.I. Let M E H and consider the partitioning M
B
C
The following
hold
(i) A, C and the Schur complement
=C-BA are all positive definite;
(ii) For all X E (A, C,S) we have that
X <
|M| and X -I1|| S ||M
(ii)) B < /
Proof. It is easy
that A is positive definite: since M 0. then it must also hold that A > since the quadratic form x" Ax can be written as y My for some vector y. For the recall the factorized form of Equation (9). Consider the quadratic form --A- (-A---(A-BBB x
Schur complement
Then y My = x'(C - BA -B''''),, which means that every quadratic form for the Schur complement can be written as a quadratic form for the matrix M, and therefore they are both positive definite.
(C - BA - B' Let y yety=
33

---

For the spectral norm bounds, since A and C are both positive definite then their norms are equal the largest absolute eigenvalues. From the variational characterizaiton of eigenvalues and the discussion above it is easy to see that ||||| < |M For the Schur complement, let z be the eigenvector such that ||C BA B x' ( BA 'B")z. Then
z (C- - -BA " B² =2= Cz - z BA 5 M
where we used the fact that BA "B" is positive semi-definite (since A -1 is positive definite) and therefore for all
x BA BBA-BB'X2 0,
where equality with zero holds only when z E ter(B").
We finally prove the bound for |B|| Let be the top left singular vector of B and be the top right singular vector. Specifically, B UEV* the economy SVD of B, then u is the first column of U and v² is the first of V Then u Bv m VBv'vv (B)u"u |B|]. Consider the
vector z
Then
A B
3 - A'v
[ 2 z Mz = (v*
u²)
C
2u Bv B|I
where used the fact that A,C positive-definite. We then obtain bounds for and ||S by observing that minn (X) min Xz and lower bound by Gmin (M) using similar argumnents.
Since A and S are Hermitian and positive-definite, we can recursively compute their Cholesky factors. This gives rise to the recursive Algorithm 2. The complexity is as follows. Proposition C.2. Algorithm 2 requires (TTmM(n)) arithmetic operations.
Proof. Let (nn) be the number of operations executed for matrix of size require time (2) Step 5 requires time TMM (n/2) + TiNv (n/2). 2), and step 6 requires TMM (n/2). The Schur complement in Step 7 requires TMM (n/2) + THERM(n/2) + (n/2)² TmM(n//) + 2(n/2)2. This becomes
Steps 4 and 8
[(n) (T()) 3TMM (Tnv(T))+2 2T
+ ((TMM O(TMM (n)))
C.1 Error analysis
The goal is to show that CHOLESKY(M)) will return backward-approximate Cholesky factor L such that the error
ILL - All S f(u, II.M)
where f is some function, and use this bound to argue about the number of bits that required to achieve it. Recall that
Bd) )).
L1 L
Lu
L L Lij + L2eL
M
B
L
L2z)
LL
L
Li
is easy to see that
Li L La Li
L L L + L
A B
B'
ILL Ml|
LL
All L B
L22L2
ClI.
(10)
which means that it suffices to bound the individual terms in the sum. The first tenm is the error of the first recursive call. Following the notation of [39], let us denote err(n) the norm-wise error of the algorithm for size We can then write
L L All <
E < err(n) (nr(n/2).
(11)
34

---

To simplify the various inequalities in the proofs we use the global bound μ from Remark A.1, where all the terms PMM (n) and (INv(n) are bounded by p(n). Finally, we also define
E (M, n) μ(() <(MM,
EI n) == p(n) (M(M) (w log(n))
(12)
We can now get first expression for the bounds of the terms in Equation (10). Lemma C.1. In Algorithm 2, the blocks L21, Lo2 of the returned matrix L satisfy:
Li All - er enr(m//)),
|m B|| < (x(M) + EBAil)) emr(1/2) + |EAAill [M]
| |
||
و + L - C|I < 2x err(n/2) + 2x( | + |L I1?
+ Eg
Proof. The enror for LnLii already described. For the second inequality, first expand Ln as follows:
Li
(BA -1 + EBAi E Li E Li B BA BA -A- + EBA.A + EBA E + LM
[BA -1 (AAi)(A EXAt)
Then
L L
B
BA -EC EBA. E E E Li
B
|| + |BA|||| + EBA
E
+ | EM
|Li
kr M err (n/2) + |BBiill M| + EBA: err (n/2)
With similar procedure we derive the third bound. We first write L21 where E contains several eтror matrices from previous computations. Then
BA -' Ln
ELg
Lg Lij [BA- L1 EL.1)
BA
L1 ELL))
= (BA -1 LLn + EL,1)
BA Ln + EL,))
BA -1 L A 'B* + EL BA 'Ln + Li B ELG E
-BA '(A + (H)
A B3* EL, BA LL + LiA -'B' 'B'E,, + EL EL..
= BA - 'B" + BA CCH A B" + ET, BA TA-'L, LL + A A + ELg Ei
Similarly
L =5 + E
=S + Es ECH
= C - BA BA--B'
Es + E
Then
|L L + L C||
BA CH A -BB" + Er,, BA -A-L
LLi + L; A - B'E Ei,, + E E + Es + ECH
BA | || + E,, BA || [n
+ E
Es
CH
S
x(M)² 2 e + 2x(M) ET. Ln| + Er,, 1| + Es + er (n.2))
2x(M)² err
(rr((/2)) (n/2) + 2x(M) EL,
| + EL,, | |E||
With this, the task reduces to finding appropriate bounds for the norms of the matrices EM EBA., Es and EL We first derive bounds for the norms of the error matrices E(.)-
35

---

Lemma C.2. In Algorithm 2, the following bounds hold for the eroor matrices
(i)
_ u. &
(2 M) |M
(i)
S u & [ M) (1 + u &INV| (2: M))
(iii)
EBAiL
S u εiNv| - M) x: <MM) + & ( M) (1 + u EiNV (5. ())
(iv)
MM
S E
1. M) (1 + u ENNv (.) Lull
(v)
E
S u- & ( M) (1 u &NNV M))² |M]]
(vi)
SUB
u vn/2 |M]] (1+ (M(M) + & (5. M)) (1 + u- EINN (3. ())))
(vii)
E
< EiNv / M) ((M)) +(2 (2+ EIN (3 M)) |Lnil|:
(vili)
[s|| S u
EiNv (5. M) (2)(22(M) (1 (1+u u- EINV (5. M)) (3x(M)++ + E3 (1, M))
Proof. Each term is bounded as follows.
(i) For |EM we have:
| S MNN((n/2) (s(A) Jan log 2 I|A |||
(..frrom Thm. A.3)
S m(n/2) u x(A) SINV log ||M-1|
u EINV (1 M) |M~
iii) Similarly, for |E ||:
E EEMM
< MMM (MM (n/2)
B||||| -1 E (...from Thm. A.2)
- PMM (n/2)
B||B|||| || -I + PMM (n/2) ul B||
EENY [|.
< MMM (n/2)
u- |M]] M + (1/2 |M] EINV (g, IM
I u E1 (g. M) +u² E1 (3 M) ENNV (H.
& ( M) (1 ENN (g M))
(ii) Using these two inequalities
bound the norm of EBA::
EBAi
BBEN
EM
B
E INV
+ E MM
M ENN [ M) |M |M +u E1 (g. M) (1 EiNv (H,
$ 6( (M) EINV
( M) & ( () + u EINV M)).
(iv) Next is the error term the matrix multiplication between the result of MM(B, INV(A)) and Lu Expanding Theorem A.2 gives
|| I pMm (n/2) |n BAill
suffices to bound |/Bij] = +
ENV
INV + MM
< BB -I + BE +
|BMM]] [|.
Using Proposition C.I. the triangle inequality, and the previous bounds, the two
36

---

unknown norms in the sum are bounded as follows
[BA
1
[ ||| A -1 S |M]] x(M),
BE NV
- [ I
I |M] u- & ( M) |M
=u uu-EIN (3. M) (MM)
This gives
|Bil| & ((M) + u ENN (2. M) {(M) + u &1 ( M (1 (1+u u. ENN (2.M)
= x(MM (1 + EINN 2 M) (1 + &1 ( M))
- x(MM (1 (1+u u SiNv M)
2
where in the last inequality we simplified (1 + 3. ( M) S (1 + u EiNN (5. M)). Then
MM
S (2M(2/22)- u- ||11||| x(M) (1 11+00 u ENN (9) (1 + u-uuu E (H.
S. Ei (2.M) (1 u ENN M))² ||,|]]]
concluding the bound for E EM
(v) The norm of MM is very similar to that of E since they both involve multiplication with BAi:
[BM < CHERM log(n/2)- PMM (MM((/22 u ||BAi]|| (..from Corollary A.2)
CHERM log(n/2) (MMM(n/2) u (MM (1 ENN (5. ()) ||M]
& - M) (1 + u εNN (g. |M]]
(vi) Next is ESUB
From Equation (5):
SUB
< u n/2
C (BAi) + EMM)
< n//
|C]|+|[[BA)B"* (|BAII(BBBBB]) JEM MM
$ nnn/2
[M ||M]+|[[BA]]]] (BAA) + ||EM
S n//
(IM +MMM (1 + u EIN (9. ()) |M +u & (9. M) (1 u EINV (8. |[M]]
u//n2. |MM] (1 <M) (1 u E ()) u & (5. M) (1 (1+u EIN (3. (2)
uv/n2 ||M] (1 (M((M +u & (G, (1 (1+ u &NNv (.MM) )))
(vii) The next term EL which we can bound as follows:
EL
EBALLI + EMMI
< EBAil [Ln|] + E
< [K(M). u EiNv (g. M) + u Er (8. M) (1 + u EINV (2.M)) [Ln +...
+ u & (5, M) (1 + u ENN (ă. M))² |L
[MM) u- EINV 2. M) + El (g. M) (1 + u EINV (5. ()) (2 + u EINV ( M)) [LL
S u EINv (5, M) k(M) + (2 + u ENN (3. (5.M)) ||L|]]]-
(viii) The final and most involved error term is the norm of Eg. Recall that from Line 7 of Algorithm the matrix can be written as
Es = BEI B + MM B" + EMM
E SUB
37

---

Once more, the norm of each term is bounded separately.
[BE] BB|
< [ E INV [||
- |M]] u εINV M) |M
<MM) ||M]] EINV (B.M).
Similarly,
E MMM
S u & (g. M) (1 + u SIN (5. ()) [|M]]
The final bound for Eg is given by the sum of the four bounds that were derived.
EEsl = |BE]] IVE + EZ B + E E
EMM
< [BET INV B*1 + E MM B" || + | MM
+ E
SUB
S u EiNv (2. M) <(M)- |M]]
E (3, M) (1 + ENN (5 M)) ||
& (2. M) (1 + u- ENN (%. M) 2 |M
n/2 || (1 (x(M) + u El (3. M)) (1+ + u- EINV (5 M))
-u-||M|- EiN ( M) .<(M)...
3 (2 M) (1 + u. EINV (g. (2.M)
1 (5, M) (1 + u EIN (2)
nn2 [ + (M(M) + uu-a E1 (). (1 u. EIN (f.)))²
=u- uM EINv G M) [((M) + -1+ u. EIN (2.M)) + EiNv (2)
[1 + (x(M) + u- & M)) (1+ u- EINV (5. (B))))))))))))))})))
= =u- u.|[M] ENN. (E.M) [M(M) + (1+ u. EIN M))
(1 (1+1 u ENN ())
+ [1 + (x(M) + u E (3) M)) (1 + u EINV (.. M))²)
< u- |M| EIN (g.M)
[<(M) + 1 + (1 + u EIN (5. M)) 2 (1 1+ (s(M) + E1 (g) (g,M)) M)))
< u M ENNv (g.M) M) 2x(M) + (1 u EINv ())² (3x(M)++-B8
+ + M))
C.2 Maintaining pesitive-definitenesss, norms, and condition numbers throughout the recursion
Given the bounds of Lemma C.2, we can now calculate the appropriate machine precision u and the corresponding number of bits such that the algorithm will not break down due to loss of positive- definiteness of the submatrices. Hereafter, we will denote by M the matrix that is passed as an argument in any of the recursive calls of Algorithm 2, and M, will denote the original matrix that needs to be factorized.
38

---

Lemma C.3. Let Mo Hi be matrix that is factorized using Algorithm 2 and k. be its condition number. Then there exist constants Ct C2. Cz. C4 1, such that n. > C1 and
C2 1 K Ing
then every matrix that is constructed in Line during the recursion, that is, every Schur complement and every upper-left block A, will be Hermitian and positive-defiite. Moreover, for each such matrix X it holds that ||X|| s 2 Moll and x(X) S 2M(MM). This value of translates to
log(1/u) - ((log(n) log(s))
required bits of precision.
Proof. At each recursive step, two matrices need to remain positive-definit: A and Š. M is positive-definite then is A. It remains ensure the same for the Schur complement. To prove this we first fix u be bounded by the value
U 2ane EiN (5. Mo) 72x((M,22
(13)
Assume that in the first level of recursion, i.e. M Mo. Per Lemma C.2, write S S Es, where
EEs <
|Moll EINV Mo) (2n((Mo) + (1+ u. ENN (o)) (3n(M))) - u & (3 Mo))
For the assumed value of u, both E1 (1.MM) S and ENN 22. Mo) I 1, and also know that 1 I <(Ma), in which case the bound simplifies to
EE| I
[M || EINV (.. Mo) (2k(ML) + (1 (3x))) + +1)
M || εIN (g. Mo) (2x(MM) + 4 (3x(M)) <(Mo)))
Moll ENN (2. Mo) 18x(M²)
|Mo EINN (g. 188(Mo)
n 20n EiNV M.) 72n(MM)²
20we Amin (Mo).
(14)
Using this, we derive the following useful inequalities
Amax (S) < ||| |Es! || I 1 2cwx |Mull.
Amin (S) Z Amin (S) - | Es 2 Amin(S) 2ow Amin (Mo) 2 (1 1/n22m0
2en0 Ammi (Ma),
naw +
x(S) < now -1 x(Mo),
(
= ((nn4))x(S) AwN (n/(4/4)
2ane
aw log(n/4)
S u(n/4)
now
-1
#( (Mo) Jonv log(n/4)
< n anw p(n/4) x(M)))
ONV log(n/4)
< n nw EINN (2.
In the above, we used the fact that for CINV 2 1 we have that -- < 2, for all 2 3, and
2oIN +1
therefore 20nv log(n/4) (n/4)anw < sww, Since the smallest eigenvalue of S is larger than positive value, we can conclude that it is Hermitian and positive-definite, and its condition number and spectral norm are appropriately bounded.
39

---

In similar manner, we can now get bound for the Schur complement of the Schur complement. In this recursive step call Algorithm 2 with M S. Let
A Bi B
$
C
Let S be the true Schur complement of S, i.e. S Ci - and let S be the approximate Schur complement of S that is constructed in line 7 of Algorithm 2 when executed on S. Let Es, be the corresponding error matrix, i.c. S1 = S + Es, Using Lemma C.2 and simplifying the terms
(1 + 1/n s 2 and PONY +1 < 2 we have that
2
S Mo)
SallMall
(2n(Mo)
<2x(Mo)
[Es, S
|||
εNN 1 S
K(S)
1+u-
n S ENNV 4.
3 k(S) +u & 11
Sww Enevl (3t .M..
Sn-ww Envv MM))
S u Moll
n nw &INV (2,
(AK(MM)) + (1 -.uu MINVEINV (o) (Gk(M)) +u 2. (g. Mo))
< u Moll nw EINV (o) (AR(MM)) + (1 1)² (6x(MM) + 2))
|Mo n CNV EINV (a. Mo) (2x(Mo)
$ "1201 EIN ( Mo) 72x ||Mol| n2 cwN EiNv (H. EiNV Mo)
72x((Mo)
Amin (Mo).
noW
The corresponding bounds S can be derived:
Amax ($))
S + [Es| ||
[S] + EEll + EE, II
|Moll
| ||
+
|Moll
(1 + 2/T Moll.
Amin (S1) 2 (1 2/n (NV) Amii (Mo),
n Fwv +2
x(Si)
n nww (M)),
(g.5,) u(n (8)) cav log(n188
EiNV
now +22
ONV log(n/8)
S u(n/8).
x((M) cINv log(n/8)
now -2
S no m(n/8)) -x(Mo) log(n/8)
nIN M,)
We can therefore identify that S is also Hermitian desitive-definite and its extremal eigenvalues similarly bounded with those of S.
we keep applying the same analysis for all log(n) iterations, as long as n is greater than some constant such that Swv + log (n log n) < 2, then for each i = 1, ..., log n it holds that
1
Es. <
Amia (Mo) S FIN
|M
now
40

---

which implies that
Amaax (8t) (1 i/n mM ||Mbl
Amin (S:) 2 (1 mm Ammin (Mo).
nww
x(S,)) - nww - x(Mo),
εiNV 27. S.) - n n EINV (# Mo).
The same bounds hold for the matrices A;, since they always originate from the top-left corner of a matrix Sj. where j < i, i.e., A; either has no errors or it inherits the errors from a matrix Sy- We can therefore conclude that, for the value of u in Inequality (13), every matrix A; and S; that constructed during the recursion of Algorithm 2 will be Hermitian and positive-defintt and its condition number and spectral norm will be at most 2x(Mo) and 2 M| respectively. The required number of bits is
(1/4+++)+ log (n 2aw EIN Mo) -7xx(M)) O (Iog(n) log(x(Mo))).
C.3 Final lackward-approxxiationn bounds and proof of Theorem 1.2
Having safeguarded the possibility of breakdown due to loss of positive-deffmiteness, we now revisit the bounds of Lemma C. and finalize the proof of Theorem 1.2.
Theorem C.1 (Restatement of Theorem 1.2). Given a Hermitian positive-definte matrix M. there exists an algorithm L + CHOLESKY (M), listed in Algoritim 2, which requires MMTM (n)) arith- metic operations. This algorithm is logarithmically stable, sense that, there exist global con- stants C2 Cz such that for all E (0, 1), if executed in a floating point machine with precision
UCHOLESKY == C1T2 ((M) logn'
which translates into (log(n) log(n(M)) log(t)) required bits of precision, then does not break down due to arithmetic errors, and the solution returned satisfies [LL MM] S c||M]]
Proof. The arithmetic complexity of the algorithm was already bounded in Proposition C.2.
For the error proof assume the bound for S uu++ from Inequality (13). Combining the in- equalities for single recursion step from Lemma C.2 with the results of the previous section and the bound for u. we have the following inequalities (recall once more that M is the input of the algorithm for the specific recursive call, while Mo is the matrix in the original call):
MM E
&
22 M) (1 + u EINV ({. |Lnl||
2
8 El
n Mo) Lnlll-
2
EBAi ||
< EIN
(3, M)
x< <(M) + u E (3. M) (1 +u n. EIN M))
o2n EINv
(3) Mo) <(MM) + 4& (g. Mo))
6n ANV EINv (g. M,) .<MM)),
EL
I εINY (3 M)
[s(M) + (2 12+u u- E ( M)) |لل|| u. ni SINV EINY (g. Mo) (2x(M)) + 9) -)- |LHll
W
SNV EINO (2 Mo) -11(MM) |Lnl.
Es|
M ENNV (# M) (2R(M) + (1+0) u.uu EINV (2, M))² (3k(M) + u. E (g. ()MM))
u 2 Moll n AN EINV (H Mo) (AK(MO) + (1+ 1)² (6n(Mo) + +1)
EINN (3, Mo) 64x(MM) -|Moll-
41

---

Using those inequalities, the first term that we need to bound from Lemma C.I is the following
وا Bll - (m(M) + EBAil) ( /2) + [EBAill [M] +
| ||
- (2n(Mo) + u 6n AN EINV (2. Mo) <MM))) (rr(n/2)
u 6n SNy EIN (3. Mo) <MMM)) Moll
u. 8 &1 (2. M,)
< ~| ||
|| Mo or((n22
< (on + u 6n ANY EINV
(2 Mo) (M))) err (1/2)
u 1m 1ww &NN
M,) <MM))
+ u 16 & (3: Mo) ||M
+ u 8 E1 (3. Mo) err(n//2)
S (Mn( + u 6n DNV εiN3 (3. Mo) is (M)) + u 8 8. (3. Mo)) err(n/22
+ u (12n SwN WNN ( Mo) xMM))) + u 16 16-1 El (5. (3) ||Mo]l.
S 22 . MM)) err(n/ (n/2) + u 28
Sow EINV
(H) Mo) -<(M)) |Moll
where in the last inequality we used the assumed bound for u, namely Inequality (13). Similarly
Li + L2gLing - CII s 2n (M) err(n/2) + 2x( EL,.
اللس + | 1 + |Es|
I 8n² (Mo) err (n/2)
4m TNN EINV ( Mo) (Lx(MM) L1
+ (u n naw EINV (2. Mo) 11m(MM) |
+ 1n nwW EIN (,) (64x(Ma) |Moll
< 8x² (Mo) (nr()
44 )2 u n NN ENN (M) (2| + (nr(n/2)
(u n aw EINN (3. Mo) 11x(M)))² (2| (nr(n/2)) ANV EINV Mo) 644((M)) |Moll
8x 8m²(Mo))
44ks (Mo)) n n RNV EIN (g. Mo)
121 (u n2 nn (5. Mo) x(M)))) err(n/2)
+ (8M)(M)) u n RNN EINv (z. (,,M,, MD)
242 (u Eiv (g. x(M)))
ANY EINs (5, Mo) (M))
|Mall
+ (44 + 12lu n aw
EINv ( Mo)) u nm EiN (g.Ma) x(Mo)> err(n/2)
|Mall <(Mo) any
EINV EINv (g. Mo) (88((M))+ +4
242x(Mo)) n CWN ENN (g. Mo)
< 173 x(,))² err(n/2) + u 394k(MM)) nwv EINV (g. |Mooll
where in the last inequality we used again (13).
42

---

We can now derive recursive formula for the main error bounds, Inequality (10).
err(n) (1)
ILL - M||
| All Li -B| L Leg CII
err(n/2) +...
22 <<(M)) ( 28n nw EiNv Mo) (<(M)) ||Moll +...
+ 173 x(M))² 2 err((1 (M4
2 nww EINV (3. Mo) |MMol
S 196) err(n + u 422k (M) nww EiNV ( Mo) ||Moll-
(15)
The base case is trivial since M is just positive real number and we can compute L = fl(MM) = (1 + M M with |6| I which means that LL = L = (1 + 6) > 0 and
M LL [M 8)1
M s 6u||Mol
we denote a = 1 196k 2 and / = 422k 422n(MM)) 2 EiNv (5 Mo). then the solution of the recursion can be written
loga-1
err(n) - a log " eer((1) + M
I los 6u M
log log Mal
= un xMM)) log Mell + log(n)n x(M) 2log
4422k(M2) n. EINV SiNV (2,Mo)
|Mool
u cm (Mo) log ||Mo||-
We can now absorb the remaining terms inside u, by tuning the constants Inequality (13), to argue that achieve any desired (multiplicatiie) backward-accuracy (0, 1) and if we set
E Cn log
for some constants C4, dg, dg. This translates to the advertised
log(1/u) O (Igg(n<<< ((M,) log " /0) = O (log(n) (og(x(M)) + +log(
bits of precision.
C.4 Reducing the GEP to a regular Hermitian eigenproblem
We can now use CHOLESKY reduce the definite GEP to regular Hermitian eigenproblem For the rest of the paper, we assume that || |S- || 1. This is a limitation since approximate the norms of H and S-1 in floating point using the algorithm SIGMAK, which is described later, and then scale accordingly. Formally, let 7j Z rewrite the generalized eigenproblem
||]]] and a 2 ||S Then we can
HC = SCA
(H (OS)C(A
H'C S'CA',
ic. it is the same generalized eigenproblem only with scaled eigenvalues. We thus safely make the unit-norms assumption. Proposition C.3 we prove the properties of the reduc- tion, which, in brief, states that can compute in ([Tmm((T) floating point operations us- ing O[log(n) log(n(S)) + log(1/e)) bits the matrix + L HL,, where is returned by CHOLESKY(INV(SS) The eigenvalues of H provably approximate the pencil eigenvalues:
( (H. <
(16)
where (0, 1) a given accuracy. Algorithm details the reduction.
It has to be real otherwise complex 1-by-1 matrix not Hermitian.
43

---

REDUCE
Input: Matrix H E H", matrix E HQ--. accuracy parameter E (0, 1).
Requires:
Hermitian and [H] < Hermitian and positive definite and S 1--1 S 1.
Algorithm: H REDUCE(H,,))
1: SINV HERM(INV(S)).
b Sinv 5 + E
2: L CHOLESKY(SIw).
LL Snnv E E
3: H <- HERM(MM(MM(L') H)).
(L'H+
H EMM
HL Es-
4: return H.
Output: Hermitian matrix H.
Ensures: (А:(H) A in (T)mM(n)) floating point operations using (log(n) log(x(S))
log(1/c)) bits.
Algorithm 4: REDUCE.
Proposition C.3. Given H H S H9 and (0, 1), Algorithm 4 executes O(TMa (n))
floating point operations and returns a matrix H. + REDUCE(H, c), such that, if the machine precision satisfies
UREDUCE
PinVa PLnVa (4x(S)) log(n))
for some constants P1, P2, Ps. and E (0, 1). which translates to
O (log(1/ (/UREDCCE) O (log(n) log(x(S)) + log(1/e)
required bits of precision, then for all i [n]
(H) (H, S) S e s |||S)
where S), are the eigenvalues the Hermitian definite pencil (H,S) S).
Proof. Let be the (exact) lower triangular Cholesky factor of S-1 i.e. LL = S-1 (which is
unique up column phases). Expanding the equation in line of Algorithm we have
A(H) = A L HI Es
A LL H LE
A (IS EINV + E CH H LE,L
A (s-'H E H ECHH LE
[L'HL + LE HL L HL L L -.ELEE
At the same time, (SS-') AL''H, Combining this with Kahan's bound from Fact A. 1 and denoting
B - -E- HL L -EE BECH L LL L,,
we have that
(H,)) (HH (A.L'LL)) ALL + B)] S C| log(n),
(17)
for some constant C. It remains bound for B||
44

---

(i E INV originates from the inversion of SINV- Let SINV = HERM(INV((S) = s'' + E
From Corollary A.3, we know that as long as the machine precision satisfies
UREDUCE S
CHERM log (wn( (x) log(n)++11
for some 8 (0, then all of the following hold: ||
I ø||S
8||s_1|
2/15-11 Σ
|Sivvll & S| S21. I1. ((()) & k: ((Siv) & 4k(().
(ii) As second step, the Cholesky factor L = CHOLESKY((INv) is computed, such that LL SINV ECH, From Theorem 1.2 we know that as long as
1
UREDUCE S 6
C)n ( )as log(n)
for some constants C12 C22 C3. then in (i), the following is sufficient:
| 8|[Sinv|l S 261|s I From the bound of <(SINv)
UREDUCE S 6 C17e2 (4x(S)) log(n)
HERM (MM(MM(I, H), HL E EMML +
(iii) Finally, we form the matrix
EMM Eg. Using Theorem A.2 and Corollary A.2, as long as
BREDUCE S ((n)
CHERM log(n (M)MM)))
[E, = H - L HL||
L HL + E L EMM L HI
IEEMI + EMM]
(n [[|]]]]
I CCHER (og(n) (x) + |-I||
EEM
MEM E
(MM(1))
+ CHERM log(n) UA (n) (p((n)
L H + EM
u ((())
+ u/((1)
Eg
MM
< 2u p(n)
u)((n) u [MM
((x)
(2 (1M(1))
2|1 + ||
3u p(n)
||]||
H
< 3u p(n)
(||S
N + [ I) |
< 3u p(n)
(2 H
< [1
<
(iv) For the final desired error need to bound
[L]]
[
|]. Taking the square:
1) " |ا
LL
LL
- = IV + E ECH
(S
E -11
= IS CC
< 1 < For the spectral norm we have that [|SIN++ E < || + || D
3||s~1|| ||. For the spectral norm of the inverse need to bound the smallest singular value of
45

---

S IN
E2 from below:
min INV E2 2 minn ( (S - | 2 Fiiu (S- S
CH
1
28|
S
2/
1 1-46x(S)
2
Combining the last two bounds, that we have
3 3||-~1
6x(S)
x(S
E2
1-1öw(S)
4öx((S
|||S|||
We now have bounds for all the required quantities to bound |B
B
L
E INV
HL
+
L
"EL"HL
HL
+
LE L
L
S (K(S)
E
INV
E
H
LE
S (*(S)
(IE)
E2
Eg
1I)
Gr((S)
S (K(S)
8||S 8||8 H|
+
- || +
128
||||||
1- 45x(S)
(*((S)
S 8 5 A(S) | H
3 + 12
45k
Setting 8 = 64c lug(n)x(S)) where C is the constant from Inequality (17) and E (0, 1), this finally
gives
6x(S)
BB| S
A5 (s(()) S S-
[|
3 12
64C (o()))))))
4e 640 Tog (M) (C())
6R(SS)
S ||]
3 + 12
64C'log(n) k(S)
c 16
1 64C log()))()) / k(S)
96
| |||
3 + 12
R(S)
15
< 64C log(n) / x(S)
||| + /k())
S < (ogg(n) ||
which we can use in Inequality (17) to obtain
(H, ( - < E
We can now gather all the requirements for UREDUCE. From the above, UREDUCE needs to satisfy
min UREDUCE I (4C)oo(())xS))
CHERM (og(n)))nvv()) lo log(n)+11
1
1
cnn (**=((x(())) co log(n)
CHERM log(n)
<e PinPa (4x(S)) Ps log(n)
for some suitably chosen constants Pl, P2, Ps- This translates to
O (O/UREEC)))) = O (log(n) log(x(S)) + log(l/e)) required bits of precision. Having established the eigenvalue bounds, note that (H, (L I). Since H and L'HL are both Hermitian, the eigenvalue bound holds also for the singular values. For scaled matrices ||S~ ||| < eigenvalues.
we obtain maximum additive error for the
D
46

---

C.5
Counterexample for floating point LU
In this section provide a counter-example which illustrates why the backward-stable LU factoriza- tion algorithm of [39] is not sufficient to obtain a Cholesky factor for Hermitian postive-definite matrix. We apply one-by-one the steps the LUR algorithm of [39] on the matrix:
3 1 3
A
(a) In the first step the matrices Lt and Uy are computed by calling LUR(A1.211) to compute the LU factorization of the first column of A. It easy to see that this returns L [1 (1(1+
where [01] u. and Ur 3 (the Uz factor is 1-by-1 matrix in this case).
(b) The second step updates the upper right corner of A, by multiplying from the left with the inverse of the top-left block of Lz. Assuming that multiplication and inversion with the number does not incur errrr,, then for this specific example A1,2 is not modified in this step, i.e. remains equal to 1, without errors.
(c) The third step updates the Schur complement, which becomes 22. + fl - [1 + 01)) =
3 (1 + 01) + 02), where |02| < u is the error term incurred by subtraction.
(d) The next step computes the Lg, UR f- LUR(A2.2) of the (updated) Schur complement. this
specific case trivially returns LR = and U = A2.2 = (3 [1 + 01)) + 02).
(e) The final steps combine together the left and right LU factors to finally return L
(1 0 1)
(f and U = 0 (3 {1 + 01)) 0b))
(1 3
1 (1+ 01 3(1 + 02) + #22(1 + +01))) which is not symmetrie.
Then, in exact arithmetic, LU
D Regularization with diagonal disorder
We now arrive the interesting part of the analysis. In Section 2 argued that if we know the gap and the midpoint then the sign function can yield the desired spectral projector. To compute the gap and midpoint efficiently, reduction, and use the Wegner estimate [138] for Hermitian diagonal perturbations to regularize the problem in the spirit of smoothed analysis. We use the following variant of the Wegner estimate.
will take advantage of the symmetry induced by the Cholesky-based
Proposition D.1. Let G be random diagonal matrix with independent diagonal elements sampled from (0, for 22/21((1n/8) and some E (0, 1/2). Then for any interval C
Pr||A + G) 1J Intn//4/0mn]
where H is fixed Hermitian matrix and I R is fixed interval.
Proof. This directly comes from the well-known result of Wegner [138], which states that that for any interval C
20 n
where [0
V2naz the supremum the probability density.
Similar results can be obtained for other classes of random matrices G. Such example is the Gaussian Unitary Ensemble (GUE) [4], which case G dense, but it is invariant under rotations which might be useful for other applications.
47

---

D.1 Sampling Gaussians and perturbing in floating point
The next step describe how to use such random diagonal perturbations in floating point. Briefly, we first assume that we have fixed number of bits for the floating point exponent and for the mantissa. If we sample standard normal Gaussian (in infinite precision). can happen that it is too large to fit in the given number of bits. This is accounted for in the failure probability, which is very small due to the decay of the normal distribution. We describe our sampling method in the following definition, which can be thought of as sampling only the most significant bits of the floating point representation of a Gaussian.
Definition D.1. A floating point standard normal sampler N(p,) t) takes as imput the number of exponent bits p. and the number of mantissa bits t, and returns a floating point number g t- fl(g)
where g is sampled from (O. 1). Following the definitions of the floating point model in Appendix
A.2, if [g) E [2-M 2M (2 )- where M 20-1, then [ gl S 2- =
uigl- If [g] >
2d (2 2 2-²) the sampler returns +INF, i.e. the floating point representation of infinity, and if
lgl < 2 d returns zero.
Lemma D.1 (Diagonal Gaussian sampler). Letry E (0, and a > 0, and assume that we want to obtain a floating representation of the matrix G, where is a diagonal matrix with independent diagonal entries from (0,o²)) Let N(p, t) be a floating point standard normal sampler as in Definition D.I, and u 2-t. Let ō be such that lō ol I MeT, where Mel S 1-cu for some (integer) constant e 1, and
=ddiagãnn .2..... Jn).
where 9: N(p, 1). 1f Cllo(loo({)) for some global constant C. then the diagonal matrix V Vffl(y- - G) G + E, has all the following properties with probability at least 1 5:
(i) V.. & (0. INF).
(ii E is diagonal and | S Mo+20 V 2n llnn///)),
(ii) [M < (1 No+2) 2ln(4n/
Proof. Let us first analyze the conditions that are necessary such that all the sampled numbers lie within the floating point bounds from Appendix A.2, specifically, Igil E [2-M 2- 2 (2 2-))],
i [n]. where M = 2p11 For simplicity we consider 19i1 E [2 M 2,2"
For the lower bound, we have that
1 2-M
Pr[la. S 2
< 2 2
V2T
where we naively upper bounded the standard normal probability density in the interval [0,2-M] by this implies that there are no subnormal gi for all i [n] simultaneously. This value of M translates to po= = (log(log(3) ())) bits.
V2= the constant function I For M logg
222
2
0
For the upper bound, a standard normal tail bound gives
Pr|lss| 2 2 2*] < 2exp(-222M /2).
(18)
For the aforementioned value of M, the probability that each one of the gi is larger than 2 is exponentially small 8/n, therefore we can conclude that the sampler does not return any INF values with exponentially high probability.
This is already enough to argue that the sampler returns floating point numbers that are not subnor- mal and not INF. However, we will need a tighter bound for the magnitude of the to bound the norm of the diagonal random matrix. Analyzing Eq (18) for M _ 1log2(2 In( 4n )). and taking a 2
union bound over all i E [n]] we have that 1gal < / In( with probability at least 5/2.
2n( An << 2 holds for all simultaneously
Conditioning on the event that all the gi are in the correct range, from Definition D.1, 7h- gil S u|gi holds for all i, Then
ü = fl(ō g.) (1 + 0jōg. = ag. (1 + N 1):
48

---

where >> is constant and Me++l -
[-(++]j This means that every diagonal element of fl(ōG)
(c+1)u
(1 + 1/0+1 )ffar from random variable that is sampled from (O,, Multiplying each diagonal element of fl(āG) with to form V (similar to the Step 4 of Algorithm 5) simply increases the (relative) error + Me+22.. Then can write
V=
where V is diagonal with independent diagonal entries from (N and E is diagonal matrix with diagonal clements E bounded in magnitude by
(c + 2)u
E S
(c + 2)u Val
It also holds that || S ov21n(4n/0) since the diagonal elements Vi, since they are just
the gs's scaled by a, and we already bounded [9s1 /2ln(4n/0), which implies that |E| s
(c+2)u 7 1- V 22n(4n//)) Finally
I
>
| 7 (1 (c+2)n (ct2)
2 2n(4n/6)).
We conditioned only on two random events, where each holds with probability at least 8/2. giving the final success probability least 1 8.
REGULARIZE
Input: Hermitian matrix A E H", perturbation scale factor E (0, 1/4), failure probability parameter E (0, 1/4).
Requires:
||A|| 1.
Algorithm:
X + REGULARIZE(A, 6).).
1: J
2v 222
2: p [Clog(log( ( [.
3: G diag(01... ................................
..n). kn).ā. +- N(p, log(1/u)
4: V f1(7 5 G) =oG + E.
> Where G diag(gn....... <- -N(0, 1).
5: X A + V B(+)
Output: Hermitian perturbed matrix X.
Ensures: See Proposition D.2.
Algorithm REGULARIZE.
We can now use the Wegner estimate to get a minimum singular value bound for diagonal shifts over grid, defined as follows.
Definition D.2 (Grid). A 1-d grid in the real line is a set of s points defined as
grid(l, r,,r, h) {l+jhlj Zzo, jh < <r.
The cardinality of a grid g grid(l, h) is denoted as ggl
The following Proposition 3.2 summarizes the properties of Algorithm 5, REGULARIZE, that we will perform efficient eigenvalue counting queries to compute the spectral gap.
Proposition D.2. Ler A with ||| S 1 be a Hermitian matrix, 8 E (0. 1/4) two given parameters, and X REGULARIZE( 7, 8). Let g be arbitrary (but fixed) grid of points in [- 2] with size ig T. For every element b; g consider the matrices =;= h, X and M -i; X + E,, where E, denote the diagonal flouting point error matrices induced by the shift. If the machine precision u satisfies
78
32(c+))nT (m(4n/6)
49

---

which translates to O 4) bits, then all the following hold simultaneously with probability 1 20:
4/3, |1(X) M((A) - 16 min MM) X
76
4n Ar
Proof. Let X A + YoG, where oG is a diagonal matrix with independent elements from (( where a = 2v22m(4n/5)) Let B = X/7 = A/7 + oG. Denote by EG the event
|| & $. From a tail bound as Lemma D.I. Pr[Ea] 2 1- Conditioning on Eg. AA/r oGI| S 1/7 1/2 8. This directly implies that A(B) C
obtain [ (-1/7 - 1/2,1/7+ 1/2). Let h4 he/r for all i 1,...,,
From the Wegner estimate (Proposition D.I), for every fixed interval I,
Pr|A(H - B) 2 1] S V 4x In(4n/8)m)
Setting 4; to be an interval centered at hi with size [4,] = 4T (m(4n/5) we have that the probabil-
ity that contains any eigenvalues at most 8/T. Since there are no eigenvalues 7t, this means
that the smallest singular value of hy B is at least half the size of It, i.e. mmin (hg B) 2 |////2
6
Hence, mmin (h, X) "Y (h² B)
76
holds well.
2T da (o(4n/5)
nT 4n In(1n/5)
From Lemma D.1. if the diagonal Gaussian sampler is called with p (log(loo ()) and log(1 1 /u) then
1
E || S 77
| S 7(1 + Nc+2) 2'
where |M+2|
S
The elements of the diagonal error matrix E(+) in Algorithm are (c(22)
(c+2)u
bounded by
S SuA + Va.l S u( + MVI s u (1 11 + net2))
Then
X x = |E (+) + E|
$ uu|A V 7
S u 7 (1 + no+2)) + Mo+2.
If S 2[c+2] for some E (0, 1) then Mc+2| < and since < 1/2 we also have that ||
| S E.
Recall that ||| A + G < 1 + 7/2 < 5/4 and thus
|E < u [/s + X < u h1| + X || + c) S u + (1 + 7/2) + +e) < 5u < 2(c+2)
Se,
assuming that
1. We now collect all the error bounds and apply Weyl's inequality to argue that
1y X) 2y x + E.) < X X E S 2c.
The same holds for the singular values ay since the matrices are Hermitian. Thus if
78
<
S Omin (h; - -X)/4
8Tn 4T In(n//)
then
Omin (14 - x + E.) - amin (h: - X) S amin (lis - -X)/2.
which implies
76
Omin (h4 X E,) 2 Tmin (h; X)/2
4T'm 4T Iu in(4n/0)
50

---

The aforementioned bound on e also gives two (loose) bounds that are useful for the analysis of algorithms later:
X
< X + c 5/44 5/4 1/32 4/3,
and also from Weyl's inequality
(:(A) (:(X) S XII | -
2 +
9
The requirement for the machine precision is
76
u S 2(c+ 2)
$
32(c+ 2)Tn / π ln(4n/8)
which translates to
o [og 78
Tn
bits of precision (recall that also required p (log((og((n/))) for the exponent of the Gaussian sampler, but this does not affect the total number bits).
For the complexity, we have n calls to ((p, t) and On)) scalar additions and multiplications there-
after, which accumulate to (n)) total floating point operations using O (log( Tn bits.
For the success probability, there are two types random events that took into account: ÆG. which fails with probability at most 6, and the Wegner estimate for 7 points h2 in the grid g. where each one fails with probability at most 6/T. A union bound over all random events gives success probability of at least 5 T = 26.
E
Fast spectral gaps and eigenvalues with counting queries
The core of our methods is algorithm efficiently approximate spectral gaps using only "eigen- value counting queries," avoiding an explicit (and expensive) diagonalization.
E.1 Counting eigenvalues
The main subroutine of our algorithm is used to count the cigenvalues that smaller than given threshold efficiently.
Lemma E.1. Let X be Hermitian, || < 2, and h E [- is fixed point that can be exactly reprexented in floating point, i.e. fl(h) h. There exists algorithm COUNT (X. h,e) which takes as input h, X, and e E (0, 4/199) with the guarantee that a (h- - -) Z e, and returns an integer n()) which is the precise number of eigenvalues of X that are smaller than h. The algorithum requires
O (TMM(n)) (log(?) + log(log(")))
arithmetic operations using
O (log(m) (og (2) log("))
bits of precision.
Proof. Let M h X and M h X E be its floating point counterpart. As the shift only distorts the diagonal elements, ||]|| max, E uumax, 1h X < 4u. By assumption
omin (M) 2 &, hence if u < e/8 then Weyl's inequality implies omiu (M+ + E) 2 amin (M) e/2 e/2.
Assume that we want to approximate n(-) by the expression
ブ() COUNT X h, round (((csSGN(M, VISGN, (SGN)))
51

---

for some aSGN, 7sGN, ESGN (to be determined). For simplicity we will use (GN(M) to denote SGN(M, OSGN, TIGNN, ESGN). Assuming for now that SGN converged the requested error ESGN, we know that (M) SGN(M) S ESGN. To proceed recall that for any matrix
[l(t((A)) tr(A)l log(n) log (n) |||||
which can be achieved with binary tree-type addition.
Now consider the error
(m(-) (-))
(r(sgn(M)) tr (r(sgn(M)) + tr(sgn (M)) (r)) + tr(SGN(M)) (((r((GN(M))
tr(sgn(M)) tr (sgm(( [M) + (r(ssn(MM) tr(SGN(M)) + tr(SGN(M)) f((r(SGN(MM)
(r(sgn(M)) (r(sgu( (M)) + TIESGN +
log(n)j
(1 + ESGN)
[ssgn(M) (M)
NESGG +
log(n)u I-log(nju
+ ESGN).
(19)
It remains to bound (Mg (ggn(M) For this we can directly use Lemma B.1 where we set [Amin x) then || |sggnh -
S = 1, H = X, and j h. The lemma states that if || |E| S
Eshift
128
X) sg X E I Shiff for some Chiff E (0, 1) that we can choose. By assumption
|mmia (h- X 2 E, which means that
E²'T
4u Cshift 128n
sufficient to guarantee that
sgn(MM sgn sgn(MM I Eshift.
For SGN, we set 175NN = 4' 0SGN = 1+5
1-ē and ESGN = 8n
1. Note that aSeN satisfies the requirement
99/100 < aSGN < as long as e < 4/199). Then based on Theorem A.4 SGN requires
N 0 (lo((- log(log( esow TIGN 1) O (log( log(log( iterations and
o (log(n) (T (log(
( essa
+ log( (pddaw)) - O (log(n) (agenn)oo@( (o) log(())
bits of precision.
If we set Eshift 1/8, then the bound for u becomes
2
32 128n'
in which case
log(n)u (1 ESGN)
1
<
1 log(n)u
Then (19) becomes
(m(-) カ(ー) <
2 = 3/8 < 1/2,
which means that we round to the closest integer then (ound((())) = n()) Note that log(1/u) O (log(2)) which is dominated by the bit requirements of SGN.
52

---

E.2 Approximate midpoint and gap
Having an efficient way to count eigenvalues smaller than a threshold, we show how to approximate spectral gaps. We can now prove Proposition 3.1, which we restate below for readability.
Proposition E.1. Problem 3.1 can be solved in time O (log( gap2 ) (esappe)
Proof. We start with To = 1/4, and make a grid g = {-1, 2/44 1/44, 1/4, 2/4, 2/4, 1, 1,444 the grid can be exactly represented in floating point). Let count(hy) denote the number of yu- distorted values X that (hy) Set jæ = arg min, (hyl count(hy) 2 k} and j++1 = arg miny [hyl count(hy) 2 k k+ Now X E [1j.-1, [iju] and X [ (. The points in the grid are equispaced each interval has length y. Since by definition X [Ax Yo, A 70) then X [hj 70 h2. + [hje-2, +1). Similarly for As+1- We have now restricted both 0 and Ax+1 inside some intervals of length 3-yo- Let us denote those intervals I [h 12.--22- h and Ix+1 [h - h+++++1]- In the next step we halve To to 71 = 1/8. We now create two grids, one inside Is and one inside TA+1: gk: = grid(hj, h/ 71) and similar for g Each of the new grids has exactly points and each point can be exactly represented in floating point. As in the previous iteration, we query all the 8 points hy of BA for count(hy) and pick new jr: arg min(hy count(hy) 2 k). and we do the same for Jk+1. There are a total of 16 queries. As before, we have now restricted X and 1 inside two new intervals 43 and 4 where [s] = | 371, which half the size of the corresponding intervals of the previous iterations. we set d equal to the midpoint of 4g, then A - Au < /2/ The same for AR+1- We keep repeating the same procedure by halving % in each step
7) which consists of the points (a useful note for later is that all the numbers in
smaller than h,, for some hg g.j 1, ................................ 10. Query all h, E g for
remains to find proper termination condition. Let us denote 7i = & and gap; = d - 2
where d and A4+1 are as above. Then, after m steps.
e = 1 + 3mm /2,
gap gap; +m
Thus, T'e x 6 egaps /3 is a sufficient terminating criterion. We don't know gapę, but we can use gap; instead, i..., rion 7'm I egap,/6. /6. But since gap, 2 gap; - -3 then the terminating criterion will be reached in the worst case when
- is also sufficient, which gives quantifiable termination crite-
gap
-37m
6 ggap
Y' S
7m
meaning that the algorithm will halt and return an e-accurate Tua gap, in at most O (oo((()))
O (logl, (og((gAp))) iterations.
Each qy-query at iteration costs ((polylog(//x)) by assumption, and there are exactly 16 hy- queries in each iteration, which gives total query cost of ((polylog(1/7)) per iteration, which is maximized for the smallest Ym O(cgap,))
We now describe the algorithm GAP that computes the k-th gap and the midpoint of a Hermitian definite pencil.
Theorem E. (GAP, Restatement of Theorem 3.1). Let H E H S E Hy and H ||
which define a Hermitian definite pencil (H. Given [n - 1], accuracy (0, 1). and failure probability 8 E (0, 1/2), there exists an algorithm
Pkk gaps + S,
which returns Ph Pk + gapı and gaps (1 + gapp, where No Act and gaps = 2
A A1+11 1. The algorithm requires
(T) log( de gapp log( s
53

---

arithmetic operations using O (log(n) (log*) (se gape + log(x(5))) bits, where X are the eigen- values of (H, S). If (() is unknown, additional (TTMm (n) log( = log( log(k(()))) floating point operations and (llg((n) log" (x(()) bits are suficient to compute it with Corollary E.I.
Proof. We first set 70 1/8 and call H REDUCE (H, }). From Proposition C.3, REDUCE requires O (T) floating point operations using ((lg((n) log(s(S)) + 1og((170)) bits, and it returns a matrix H that satisfies | (H) (H, S) S 7.
We then call X < REGULARIZE(H, ,, 6n), where 60 6/2 is the initial failure probability. From
Proposition 3.2, for all it holds that (X, X (H) - 9 /6. Therefore, all the eigenvalues, which initially lie in
[-1, 1], are distorted by at most %ao
Next, recall the counting query model of Proposition 3.1: in the first step construct grid g grid(-1, 1 + 70: 70), where igl = 10. From Proposition 3.2, the regularization ensures that for every hg E it holds that o min (hg X + E) 2 for Eo =
Iu/4n/ (this is included in
the success probability of - 60 of Proposition 3.2). Then execute COUNT (X, &o) 6n) for every hy. We continue by halving at each step both and 6, constructing the corresponding grids as per Proposition 3.1, and counting eigenvalues over the grid. In each iteration after the first one, keep track of two grids g with size gl ((1) each, and therefore we can ignore [g] in the complexity.
After total of m O (log( ) iterations, we have that im ((pp)))) and 8m
(oge gapş). Invoking Lemma E.1, and plugging in the value of em each call COUNT dur- ing the entire algorithm costs at most
o (Thma(() (l08(=#)) + log(lo((=")) (0)) - (Tima (n) log( log(< mm))m))
arithmetic operations using O (log(n) ) log" ( log( m)) = O (log(n) log (yn²
bits. The last iteration where "Ym. = ( gap;) and 8m = Oōge gapx) is the most expensive, which gives total of
(TMM(1) (n) log( dg n
gaps )mm))mmmm))))) m = O (IMM (n) log( d Sap
log( gups ))
arithmetic operations and
O (log(n) log"
Vdge gaps
Nogegap
bits. The failure probability 8n in the first iteration and halved each subsequent iteration, in which case a union bound converges to 26g =
To be able to call REDUCE in every step, we need to approximate x(S in order to set the machine precision appropriately. For this we can use Corollary E. which is detailed in the next section, and it returns E (x(S)) with probability 5- It requires O(T log( RN(SS) log(x(S))) floating point operations using (log(n) log" 4( (x(S)) bits.
E.3 Singular values, singular gaps, and condition number
The gap finder can be extended compute eigenvalues and singular values and singular gaps or arbitrary matrices instead of eigenvalue gaps of Hermitian matrices. The following algorithm in Proposition E.2 computes the k-th singular value of matrix, for arbitrary [n]. using counting queries.
Proposition E.2 (SIGMAK). Given matrix A E Cm 2 m with ||||| < 1, an integer k E [n]] such that rank(A) (A) 2 k, an accuracy (0, 1), and failure probability 8 (0, 1). there exists an algorithm
0 + SIGMAK(A, c,),
54

---

which returns value ox (1 = (A) with probability at least 1 - 8. The algorithm executes
o (Tum (n))
m (m+ Iogg(( log( 1.
floating point operations and requires O (log(n) log"( (x)))) bits.
Notes: If m > n. then we can apply the algorithm on A A. |f 1, we can use Theorem A.7 to approximate ||A]] and scale accordingly. The rank(A) 2 k assumption can potentially be omitted by a more sophisticated counting query strategy.
Proof. We will solve the problem with counting queries, similar to Theorem 3.1, halving and ó in every iteration. We start with To = 1/8 and 6g = 8/2. First, we divide by two, to ensure that 1/2. Then we construct A = HERM HMRM(MM(A)) A, A*)) = AA" /4 + EM To perform this rectangular multiplication practice, we partition A in m/n blocks of size n n each, perform the individual multiplications and sum the results. We can ensure that |M || I 7//4 if we use
< ||
Ollog(20 (2) bits. Then from Weyl's inequality (Ax(A) Ax /4) & 70/4. We then call
X C REGULARIZE(A A, 70, 6o). Conditioning on success of Proposition 3.2, for all it holds that (X((X) X.(A)) z 9 176/16, and therefore |A (X) A:(AA" /4) S 70-
Let X = (X) and A, /4AA"/4). Since all the X are non-negative, then all the X [-70 1+ ro]. We thus construct grid g = grid(0, 1 + 70: 70) = {0. 1/8, 3/.. 9/88. We now need to perform counting queries similar to Theorem 3.1. Since we conditioned on success of Proposition 3.2, REGULARIZE ensures that for every h2 E g it holds that min (h2 X Eng)) 2 so for Eu /
Then we execute C₁ + COUNT(X, h3, for every hy E Eg (in total
(1)) calls to COUNT). We set jk be the smallest such that C3 2 k. Now X1 E [boer - Since Xe E [A 70. d + 7ol. then d E [h -o. hy. + yYl = [h hje+1l := Ix. where |4| 3% (we can always ignore the negative portion Ig. if there is any).
1/8, and we create another grid inside Ig. g1 = m) The new grid has at most points. We query all the points h of g1 for C, + COUNT(X, h1,61), and we pick a new j1 = arg min [C, 2 k}. Now 0 lies in- side new interval I with L4l = 3y1, which is half the size of the interval of the previous iteration.
In the next step we halve To to 71 grid hi
we set A equal to the midpoint of I (which is always positive), then A Au S 3//2. We keep repeating the same procedure by halving 7/ and &e in each step i.
After steps have that
d A + 3mm
A sufficient termination criterion is 7m edg. Instead of d (which is unknown) use de instead. Then ? < e A--3ym is also a sufficient, which gives a termination criterion 7m S t that we can actually calculate. In turn, this termination criterion will be reached the worst case when
A 3/m/2
7'
7m
9
meaning that the algorithm will halt and return an e-accurate A in at most O (log( (log(;2)) (-m)
O (log( iterations.
Note that in each iteration we make total of @(1) calls to COUNT. From Lemma E. if we plug in the corresponding E,, each call to COUNT during the entire algorithm costs at most
o (1) log( (,66
arithmetic operations using O (1og(n) ) log* (g:0)) bits. The cost is maximized in the last iteration where 7m ((cAs) and 6,. O6pe))s) which gives
o (Tum log( -m)m O (Thm (m)) log( oyeAo log( 1og(=t)))
55

---

arithmetic operations and
o (log(n) log (B,EXe)
bits. As in Theorem 3.1, since the failure probability is halved in each iteration, a union bound converges to 26n- =ô. 6.
Now we have that | 1 Akl S cdg. Recall that (A) = AA(AAA) 2 since we computed
A and divided by two. If set 0z 2/ Am then
k
2 (1+ c) of - (1 (1-6) -)0 S ok (1+c)ox
(1-c)
4 4
where the last holds since vI E < 1 1+e and 1 1 vI- for E (0, 1). Note that we lazily assumed that we computed 0い = 2 / exactly, which does not hold, the square root introduces a machine precision error. The calibration is left as an exercise. If we replace D with of and 60 6/2 the complexity bounds we get the desired
O (T)m(()) log( Teoa log( edd))
arithmetic operations and O (log(n) log"( (x bits.
We now have a tool for the condition number.
Corollary E.1 (COND), Let A Hi, Given 8 (0. 1/2). we can compute such that r((A) S R S 32x(A), with an algorithm K + COND(A, 5) in
(TmM(N) log( (a(((((( log(ss())))))
using (log(n) log"( (mngA))) bits of precision with probability at least 1 5.
Proof. We first compute Σ E [0.9|||], A| using Corollary A.1 with parameter 8/3 in O(n" log(n) log(1/8)) floating point operations using O(lo((n) bits of precision. The algo- rithm succeeds with probability at least 1
3 We then scale A' + A/M where M is |A'll
the smallest power of two that is larger than 4E.
This implies
I and also 3.6
Tmin (A') (min(A) (A) A|
Thus, it suffices approximate amin (A')
IoN(A 3.0A(A) and scale it to obtain the desired approximation for x(A).
3.6||AT|
We now call J'in + SIGMAK(A', 1, 8/3,3 which succeeds probability 1 5/3 and returns o Fiii (1 + (iimni')))) It requires
O (Ihma((n) (10g( (oss (AT) log( (mmn(()))))))))))))))))) =0 (IiMm(1)))g) mx(A) log(m<(A))
arithmetic operations and
O (log(n))ogg log (5
( (5)))((A)))) O (log(n) log log((nx(A))
bits. Then we can
E [i((A), 32k)))22(A) The success probability is =1-8. min
E.4 Proof of Theorem 1.1
We now have all the prerequisites to prove that Algorithm 1 provides the guarantees of our main Theorem 1.1, which we restate for readability.
Theorem E.2 (Restatement of Theorem 1.1), Let S) be a Hermitian definite pencil of size n, with H|] || ||= < 1. and \ < X2 S & An its eigenvalues. There exists an algorithm
I + PROJECTOR(H, S, c).
56

---

which takes as input H. S, an integer 1 S k S 1, an enorr parameter E (0, 1) and returns matrix I1 such that
Pr n H1 Se 2 121-1111 1/n,
where II is the true spectral projector on the invariant subspace that is associated with the k smallest eigenvalues. The algorithm executes
o (Tum((n) (log( gaps. log( gaps + log(nx log(k(S)) log (log( m )))
floating point operations, using
O (loB(n) (log" gapp + log (nn((S) - ( log(
K(S) egap,
egap
bits of precision, where x(S) = ||||||S* and gapk AR+1
An. Internally, the algorithm
needs to generate total of at most (n) standard normal floating point numbers using additional (log(log(())) bits. In the regular case, when L, the ((log(vx(SS) (og())())) term in the arithmetic complexity and the O(log" (nx(S))) term in the number of bits are removed.
Proof. We first compute PA. gapg.k which satisfy the requirements of Proposition 2.1 using COND and GAP. Afterwards we use them to call PURIFY to obtain the spectral projector II, such that III II S e. We first call COND(S, ). which, from Corollary E. requires
O (TMM(N) (n) log(nx(S)) log(((S)))
arithmetic operations and
O (log(n) log (ogg(nx(SS)
bits. Then, given the result of COND, we call GAP(H, S. k, 8 2n) which requires
o TMM (n) log( a_
log( gapp )))
gap,
arithmetic operations and
O (log(n) log gap,)) ( log(x(()))
bits. They both succeed the same time with 1/n probability. PURIFY gap, K. requires
o (T))(n) (log( gaps + log(log( x(S) gaps
floating point operations and
(log(r) log
log( x(S) gap,
gap,
bits. The total arithmetic complexity is therefore
O (Thma(n) (Log( n gaps
log( s.) + log(nn(())) log(k(S)) log (log( m(9)
egaps
and the bit requirement is
o (log(n) (log* (1og((=a)) EAP,
+ log* (nk(S)) log" log( -app
Notice that the term appears only inside log log in the arithmetic complexity and in only one the log factors of the number of bits.
D
57

---

F DFT background
Density Functional Theory (DFT) [88], which awarded the Nobel prize in Chemistry in 1998, is considered as one of the most common methods to perform electronic structure calculations thanks its ab initio character. That means that no input besides the initial atomic structure is required to determine the relaxed geometry, energy levels, and the electron structure of a given material, in their bulk or nanostructured form. The main idea behind DFT is to describe the system and its properties by the electron density only. DFT calculations are widely used in industry and in academia to predict, for example, the properties of novel materials [102] or to optimize the performance of batteries [69], solar cells [142], or nanoclectronic devices [85]. Many scientific libraries implementing DFT algorithms have been developed, which persistently occupy supercomputing clusters and receive up to tens of thousands of citations annually at the time of this writing [89, 59, 127, 74]. Despite the success of DFT calculations, important theoretical aspects of their underlying algorithms remain unclear.
Assume a system that is composed of a set of na E N atoms positioned inside a fixed three- dimensional domain. The Kohn-Sham equation in real-space describes the electronic wave function (r) of the system
H(r)v(r) Ev(r).
(20)
where r (z, y, 2) is any position in the domain. Here H(r) is the so-called tonian operator of the system, which can be written in the following form:
Hamil-
1 v2 Fiin(r) + H(r) -
/ n (rrr) r'l dr' + (r)e(r),
|r
where the first term is the kinetic term, the second term accounts for electron-ion interactions, the third term represents the electron-electron interactions, which are solved at the Hartree level through Poisson's ss equation, and the fourth term is the so-called exchange-correlation term. The electron density n(r) the main observable of interest in DFT. It is related to the solutions E, and ()( of Eq. (20) through
n(r) = Σ
(21)
where EF is the Fermi level of the system, E, is the energy level (eigenvalue) corresponding to the wave function (r). and f(E; EF) is the occupation term. At zero temperature (T' = OK), [EE; Ep) is equal to one for occupied states (E, < Ep)) and zero otherwise." The Ep, E, and J/ are the unknown quantities for which there exists no analytical solution, except for some spe- cial cases. Approximate numerical solutions can be sought after expanding the wave functions into a suitable basis. In this case, the v are approximated by a set of n = (m,,) basis functions (r), E [n). typically localized around atomic positions.' Common basis sets that lead to effi- cient algorithms to solve the eigenvalue problem of Eq. (20) include Linear Combination of Atomic Orbitals (L.CAO) [32], Gaussian-Type Orbitals (GTO) [22], or Maximally Localized Wannier Func- tions (MLWF) [101]. In this context, the wave functions take the form of linear combinations of the localized orbitals. In particular, for all E [n]]
(b) Σ (r))
(22)
where Gj are the (unknown) complex coefficients. This expansion can also be written in matrix notation as
Un (r) (r)
C11
C12
Cin
X1 (1)
C21
C22
Can
x2(r)
(())
= Cx(r),
(23)
(r))
Cm1
C22
Con
'The definition of f'E; Er) can vary to allow for partial occupations, iee. at nonzero temperature 0 can be replaced by the Fermi-Diral distribution I+exp[[E, -EF)/RT)) where ki is Boltzmann's 's constant. This case is not considered this work. A common alternative is to approximate the electronic wave function using Plane Waves. In this work we focus on localized basis sets.
58

---

where y(r) E C" and x(r E C" are (complex) vectors and C n x n matrix. By expanding the wave functions into LCAO, GTO, or MLWF basis, the Hamiltonian operator becomes matrix
H (xx|H)))) (t)())())
Similarly, an overlap matrix S must be introduced
S. = (xllx)) - xi (t),
which is equal to the identity matrix I if the basis elements are orthogonal (e.g. S = I for MI.WF). The entries of H and S represent the interaction between two basis elements that can be located on the same or on different atoms. Furthermore, due to the complex conjugate property of the inner product, both H and S are Hermitian matrices. Moreover, S is positive-definite ([103, Chapter 3)). Putting everything together, the expansion coefficients Gú are obtained by solving the following generalized eigenvalue problem already formulated in Eq. (1):
HC = SCA.
The so-called density matrix P can be derived from the C and A solutions
P Cf(A; CF)
(24)
The density matrix possesses interesting properties. The matrix PS a projector matrix and there- fore is idempotent, i.e. PS (PS)². Moreover, it can be used derive the electron density at any position
n(r) = Σ
-x'))xx)
(25)
LJEn
where x(r) C is vector such that (r) xx(r), i.c. x(r) contains the values of all the basis functions at position r. Given these basic definitions, we next give an overview of existing algorithms.
F1 Linear scaling methods
Because of the strong localization of the (r)'ss, the (i.j) interactions rapidly decay with the dis- tance between the atoms and can be discarded beyond pre-defined cut-off radius. As a conse- quence, H and S are often banded matrices with sometimes a sparse band. It should also be noted that at equilibrium, the number of negatively charged electrons present in the system (which corre- sponds to the integral of n(r) over space) must exactly compensate the number of positively charged protons Zp. This allows for the computation of the Fermi level: from Eq. (25) and the definition of S. it follows that trace = Zp. which means that one can perform a binary search on the Fermi level Er until trace (PS) Zyp is satisfied. These observation has led the development of so-called linear scaling methods, which aim to exploit the sparsity of the matrices to obtain faster solutions [143, 57, 86, 87, 56, 62, 63]. Typically, the density matrix and the Fermi level it- cratively computed. until all required constraints are satisfied, i.c., the density matrix should be idempotent and trace(PS) function using Newton iteration [136]. To optimize performance, heuristics such as truncating ma- trix elements with negligible magnitude are applied, the cost of decreasing the solution accuracy. Empirical evidence shows that such methods tend to scale nearly linearly to the system size, how- ever, they do not exhibit provable theoretical guarantees. Similar methods have been studied based on finite-diffet miite-differences. which aim to approximate the Fermi distribution with polynomial expansions [145, 144]. Finally, the work of [65, 137] offer the possibility to achieve lower complexities than Onn)) for both the eigenvalues and eigenvectors of structured matrices, by leveraging the fast mul- tipole method. However, an end-to-end stability analysis of those algorithms remains open. Indeed, are not aware of any end-to-end analysis with provable approximation guarantees that has lower than O worst-case complexity for any of the aforementioned algorithms.
Zy. This can be done, for example, by approximating the matrix sign
F.2
Density matrix
We can directly apply Theorem 1.1 to compute density matrices in DFT.
The following Theorem summarizes the analysis of Algorithm 6 for the density matrix.
59

---

DENSITY
Input: Hermitian definite pencil H E H"" E 4 number of occupied orbitals k.
Requires: |||| 1, ||S- 1, E [n- 1].
Algorithm: P DENSITY(H. S,
1: II < PROJECTOR(H. S, )
2: SINV f INV(S).
D SiNV = 8-1 + !NV
3: P + HERM (MM(î, MM(S )))
b P = fi ii(ii +E!M) + EMM
4: return P.
Output: Approximate density matrix P.
Ensures: |P PP-Pl Pll S with probability least 1 O(1/n).
Algorithm 6: DENSITY.
Theorem F.1 (DENSITY). Given Hermitian definite pencil (H, S) with |H]|, S 1. an integer k [n- 1] denoting the number of occupied states in the system, and (0, Algorithm 6 returns a matrix +- DENSITY(H, S, k c) such that
P PP-
where P is the true density matrix of the system and succeeds with probability at least O (1/n). The floating point arithmetic complexity is
" (Tmm(()) (log( gaps
log( gaps
log(ns(S)))og((((()))
and it requires
o (log(n) (log"( gapp + log (nn(S)) - -log* -oog"(
(egap, log( k(S)
gaa,,
bits of precision.
Proof. In the first step of the algorithm compute IT such that
= éTT:
where is the spectral projector on the invariant subspace of the definite pencil associated with the A smallest eigenvalues.
Then we invert S, and the final step of the algorithm computes P. Unrolling all the computations, P can be written as -IS-III" + Ep where Ep is Hermitian error matrix (yet be detailed), is not hard to verify that IIS 'I" is equal to P. i.c the true density matrix.
It thus remains to bound |Epl We can directly force ||E|N|| s <n|[S-11| s cn by setting S Then from line of Algorithm
TI (n)()))(S)) (ww log(n)-
P II ПE
MM
MM + E
= (II En) (s
INV
(II En) + (II + En MM + E
IIS IIS-'T + HE INV
+ S- + ENN Ej En (s- 1 + ENN II +...
En (s-1 + E E + (II En EM + MM
bound each term as follows.
i IIE INV II"
<
|E|V < en.
60

---

(
II s E.M' Ei
- (1+
en 2 2em
(iii)
En (s* " + E En
- (1 + m क - 20
(iv) EEM || S (u((n) IT
S up(n) (IS | (1 [nll)) ( <<(++m)
(n) S 4cm.
(v) E MM I ( || |SNv
E
< (n(1 + + [(1 +n + + + 4c)) S cn
2. (4 + 4 I 1641-
Putting everything together we have that
P PP-Pl S en + 2E + 2e + Δеπ 16Eл 25e
25
=32 S
32
Note that the aforementioned value of u is already dominated by the requirement of PROJECTOR. The same holds for the arithmetic complexity. Therefore, the arithmetic complexity, the success probability, and the number of required bits, similar to those of Theorem 1.1 up constant factor.
F.3 Electron density queries
Having good approximation of the density matrix, at any point of interest. A direct application of Theorem F.1 gives the following.
use it to query for the electron density
the atomic domain, let x be vector x Xn(r)), where Xi are the basis functions that were used to construct the Hamiltonian and overlap matrices H.S. Let P DENSITY c, k) where k is the number of occupied states in the atomic system, and let P be the true density matrix. Let n(r) x Px be the true electron density at and ñ(r)= = Px. Then long DENSITY succeeds.
Corollary F.1. Given position (xi(r) (r)
[n(r) m(r) s 8c|[x|12
Proof. First we recall that E (0, 1/12) and note that |P||| S PI + S 1 + e. The quadratic form x" Px is computed in two steps: y fl(Px) and 7()) fl(x"y). From [70], Eq. (3.12)
and Lemma 6.6, we know that the matrix-vector product satisfies y Px ||[xx - yll S Vrnral P |||x||, and Tn. := nu
where ||e|
i- For the chosen u we have vTnn I 2e. From
the same equations, the subsequent dot product satisfies 7i(r) = fl(x"y) = x"y + where |e| S
mh ||x||||||| - 7n ||x|||
(IIPx|| + ||el] S mm || (1 + von..) S 4e|[x|12 Then (n(r) 7(r) =
xx (x"x y + c)
xx Px x Px x'e cl Px x P x |c| 2c |x|
2e| + 4e |c| 8c|
This be generalized for many different points r. If N is the number of those points and N = (n) (which is typically the case in applications), then one can stack the vectors x(ri), x(r2), ... x(rx) the rows of a matrix X and use fast matrix multiplication compute all (rr)) simultancously, instead of querying each n(r;)) one-by-one.
If N is asymptotically larger than the system size, i.e., N (ppo)y(n)). the query complexity can be reduced by applying the Johnson-Lindenstrun Lindenstrauss (JL) lemma [78], at the cost of a polynomial dependence on the accuracy. In brief, since the matrix PS is an orthogonal projector, and therefore idempotent, can be deduced that
x Px = x x'PSP x 'L L Px = |x'PL 1
where L is the lower triangular Cholesky factor of S Then |L PL can be approximated by ||x PT] where T is n x sparse random sign matrix scaled by 1/vF with = Olog((N) (N) columns. The representation of T requires 0(nxr) = O(n (oe)) random bits. Alternative random matrix distributions satisfying the JL property, with varying sparsity, runtime,
61

---

and dimension requirements exist in the literature [3, 82. 110]. The value of r is worst-case optimal [93], however, tighter bounds have been studied in the literature [21. 126, 23, 49], thus potentially leading to more practical algorithms. The full analysis is omitted.
G Deflation and bit complexity of PCA
Since its introduction in the early twentieth century [116, 72], Principal Component Analysis is one of the most important tools statistics, data science, and machine learning. It can be used, for example, to visualize data, to reduce dimensionality, or to remove noise in measurements; cf. [79, 45] for reviews on the vast bibliography. In its simplest formulation, given a (centered) data matrix X E IR*x,, the goal is to find k-ddimnnsional embedding Cg, where l < n, that maximizes the sample variance, which can be written an optimization problem
C₁ = arg
max C=Txxk
tr(C HC),
(26)
where H X E Rax= is the sample covariance. It can be shown that the solution C₁ corre- sponds the principal & singular vectors of H. i.e. the ones that correspond the largest singular values. Evidently, since the sample covariance is always symmetric and positive semi-defnnit, this can be written as Hermitian eigenvalue problem
HC =C=CA.
(which is indeed definite GEP in Equation (1) with S I). This way can project the data in dimensions by computing XCR, preserving as much of the variance in k dimensions as possible. Classically, PCA can be solved by approximating the SVD of H via diagonalization.
G.1 Vanilla PCA
With this short introduction, we can now describe how to enhance the analysis of PCA algorithms by using our results for spectral projectors obtain forward-error guarantees for all the SVD- related computations. Using our Theorem 1.1, the SIGMAK algorithm of Proposition E.2, and the DEFLATE algorithm of [19], as described in Algorithm 7, we can state the following result for the classic ("vanilla") PCA, which provides forward error guarantees in matrix multiplication time regardless how large k is. This can serve as the backbone for more advanced algorithms that are analyzed in the next section, which typically first use random projections to reduce the matrix size to approximately ( x k) and then perform the SVD computation on the smaller matrix.
Remark G.1. For simplicity, we have assumed that the covariance matrix contains no errors. This might not be true if H is computed numerically from the data matrix X. It is not hard to extend the algorithum to take this into account as well, and nothing changes in the analysis (except some negligible constant factors). To handle errors in the input matrix one can adapt the analysis from the next Section G.2.
Theorem G.1 (PCA). Let H be an n. X n. symmetric sample covariance matrix of a centered data 1, k E [n] is a target rank, and c E (001) an accuracy parameter. Then we can compute a matrix C with is columns such that X (1 + +X XC&C |. where C₁ E BRnk contains the top-k (right) singular vectors of X in
matrix X E mRx, i.e. H = X X H <
o (Thma(()) (108( log( outt 1) log( gaps
log( (1 + log(log( ca2+1 gape
gaps
arithmetic operations using O (log(n) (log'( (egip)) log log'(( (mmm)) + log log(.a).)) bits of preci- sion, with probability at least ((1nn).
Proof. We first compute ō SIGMAK ks 11,, 1)) such that, from Proposition E.2, 0k+1 (1 2) (H) with probability 1/n. requires
O (ToM(()) (n) log( o a_ log (+4
62

---

PCA
Input: Sample Covariance Matrix H Rax of a centered data matrix X E R=x*, target rank k accuracy (0, 1).
Requires: ||H] I LI,E in 1].
Algorithm: C2= t PCA(H, k,c)
1: 0 (IGMAK(H, h; 1: 3
2: é' -.///
3: II; PROJECTOR(H,I 20 )²)
4: C2 <- DEFLATE(II, k, 120 e.).
5: return Ca.
Output: Approximate principal component matrix Cs.
Ensures:
X 1-XCC,, s (1 +1+66) + Cf with probability at least O(1/n)
Algorithm 7: PCA.
floating point operations and O (log(n) log (a bits. Then we set + c0x+1/4. In the next step, compute II PROJECTOR(H, I, k, 20 )²), which, from Theorem 1.1, returns
2²
spectral projector II that satisfies I I& I 20 )a with probability 1/n. It requires
O (Tama(n) (logí gapy
log( gaps
+ log(log gaps
n =0 (Tim((n) (log(
log( = + log(log( car-1 gups
gaps
arithmetic operations and O (log(m) log"
O (log(n) log" (cors" ssap.) bits. Then we
gapp
use DEFLATE to compute the matrix C₁ DEFLATE k, (20 )2,c'). From Theorem A.5, DEFLATE succeeds with probability
(20n) (20
1-1114
and it requires ((TMM (n)) floating point operations using Olog(n//e) O(log(n/(ox+1e)) bits, and it internally generates (n²) random complex normal variables. On success returns C² E Caxk such that | - C/ S c', where Cg matrix whose columns form an orthonormal basis for 11. We can then keep only the real part of C²2, and set C2 <- Re (e(C).). The spectral norm can only decrease by removing the imaginary part, therefore |C C& < é' holds well. Now write C₁- Cs E where |E|
S e', which means that
x C-X
= X X(Cr + EC + +E)
$ X XC
+ X 1
+ XC,
||
X
< x XC
+
2
+ +
- C.C. +3
4
(1+ ()) X XC
where in the last we used the fact that
T H)) < 3 ( = 30 (X)/2 30 (X)22 - XC,C
The result follows by rescaling by a constant and by summing together the individual algorithm complexities and bits.
63

---

G.2 Low-rank PCA
In many applications of PCA the parameter k is chosen to be small, i.. k < d. This has led to an extensive research area on the so-called low-rank approximation algorithms for PCA. One of the carliest works which introduced this type of randomized low-rank approximations is [53, 54]. Some other landmark works in the field include the analysis of randomized PCA and low-rank approximation [100, 67, 120, 33, 34], and the pioneering Block-Krylon PCA of [106], which is essentially optimal in the matrix vector query model [125]. The approximation accuracy of low-rank approximation-basee PCA methods is often measured with respect to the spectral or the Frobenius norm enror, i.c. the matrix C₁ that is returned should satisfy Cj C2=
Lxx and:
X C
-- XC
(27)
(2.F)
(2.F)
Most of the aforementioned low-rank approximation type algorithms assume exact arithmetic. A subtle detail is that they often rely on the computation of an (exact) SVD of some smaller submatrix, which is not realistic due to the Abel-Ruffini theorem. This can be justified for practical reasons since SVD can be approximated in polynomial time to arbitrary accuracy using classical solvers [64]. Quoting [67]:
"Techniques for computing the SVD are iterative by necessity. but they converge so fast that we can treat them as finite for practical purposes."
However, if one wants to rigorously prove forward error approximations and end-to-end complexity upper bounds, there necessarily a dependence on the singular value gap that separates the principal invariant subspace from the rest of the spectrum (consider a small perturbation of the 2-by-2 identity matrix and k: 1 as a straightforward example).
Assuming an exact SVD algorithm, the seminal analysis of the Block-Krylov PCA of [106] can in fact provide per-singular vector guarantees, which are much stronger than the classical norm-wise bounds. Then arithmetic complexity of Krylov PCA depends on poly( (.t))) It is favorable for coarse accuracy (i.e. when 6 = =(1)) and small rank k. not suitable, however, for larger k. F.g., it can be the case in applications where we need to keep k = n/20 of the original dimensions. The same holds when higher accuracy is required, c.g. when = 1/ poly(n). Then the complexity of these methods is already higher than standard eigensolvers.
In finite precision, the landscape is even less clear. At the time of this submission, the only work related to low-rank approximation PCA that we are aware of with end-to-end bit complexity bounds [107] for the approximation of matrix functions applied on vectors. They prove that the Lanezos method can be stably applied to compute a vector u such that the quantity ||Aul approximates || This is aapraadd-approximate solution for the top-1 singular vector. Concurrently with this work, and independently. [80] analyzed the bit complexity of Block-Krylon PCA and achieved similar bounds as ours, albeit with different techniques. However, [80] did not describe how to compute the condition number which is required in order to adjust the machine precision. Our analysis covers the computation of all the involved parameters.
For general invariant subspaces (k > 1), consider one of the simplest randomized low-rank ap- proximation algorithms, often referred to as "subspace iteration" [67] or "simultaneous iteration" [106]. The algorithm first samples a matrix G E Rxx with i.i.d standard normal elements, where I == (kk), and computes Y + (AAT) A where typically s log(min m,, n))) It then returns the Q factor from the economy QR factorization of Y. It can be shown (see e.g. Corollary 10.10 of [67]) that
[[ 1 )AI S - Apll-
The returned matrix has slightly more than & columns, often referred to as "'oversampling." but, importantly, there no explicit SVD involved in the computation. Even in this case, floating point arithmetic already spoils the approximation guarantees: we cannot compute Q without rounding errors. Stability analysis of QR factorization is typically carried out the backward-errov sense [70, 39], which makes the analysis of low-rank approximation PCA algorithms even more complicated.
64

---

G.3 Bit complexity analysis of Block-Krylov PCA
We now analyze the seminal Block- Block-Krylov iteration algorithm of [106], which is listed in Algorithm 8 for convenience (in exact arithmetic). For the floating point analysis the following methodology is used. The main result is stated Theorem 4.2. We shall denote by (X. q) the cost of multiply- ing X with dense matrix with y columns with a numerically stable multiplication algorithm, like the one of Theorem A.2. Formally, for the rest of this section, we assume the following subroutine.
Definition G.1 (MMX). Let X E Rmxx be the input matrix of Algorithm 8. We assume subroutine (M), which takes as input a matrix B with k columns. It returns a matrix E RRxr such that
|CC - XBI| S poly(mm, X B
in total of (i)x(()) floating point operations, using most (log(1/u)) bits of precision. We also assume MMX (B) which approximates the product X XBB with the same cost and approximation bounds.
The standard inner-product based algorithm as well as block variant of the MM algorithm of Theorem A.2 satisfy this definition.
Block-Krylov Iteration (exact arithmetic)
Input: Data matrix X E Rmxn target rank k accuracy e E (0, 1).
Requires: Exact arithmetic.
Algorithm: Ze Block-Krylov Iteration (X, e)
1: e log(m)).
2: G - -V, 1)nx*
3: K < XG, (Xx JG,...... .XX.
(XX yx
4: Q R < Fconomy- (KK).
b Q E mogk
5: M Q xx Q E Rokxok
6: U2 top-k singular vectors of M.
7: return Z₁ QU.
Output:
Approximate principal component matrix Z.
Ensures: | Z.Z2 X|| $ (1+ e) XxIl with high probability (in exact arithmetic).
Algorithm 8: Block-Krylov Iteration (Alg. 2 of [106))
Step 1: Constructing the Block-Kryloo matrix. We first observe that we can scale the matrix G in Algorithm 8 to have norm at most one, since this does not affect the Krylov basis. We assume that the (scaled) Gaussian matrix G in Algorithm greatly simplifies the analysis and it can be addressed by using the floating point Gaussian sampler of Definition D.1 obtain small error with high probability.
given exactly. This not realistic, but
Given G and X, we can construct the Krylov matrix K using O(1) calls to MMX(G). The first call X, MMX(G) returns matrix X, XG E, where |E < u (mpk)) since we assumed constant norms. we keep performing multiplications re- cursively to build the Krylov matrix ||K|], after O(g) multiplications it holds that || K| < ((poly(mkqq). The total cost is ((qTTMxx()) floating point operations. If X sparse and we the standard inner product-based algorithm. TMMxx (k) O O(knzz(X), If X is dense and use Theorem A.2 then Oggmk" 2) floating point operations are sufficient.
< ((polymmk)) |||
Step 2: Condition number of the Block-Kryle matrix. So far we have approximated K by K K Ex.
In the following steps we will need an approximation for the condition number We can get such approximation using variant of the COND algorithm (Corollary E. which internally uses
65

---

the SIGMAK algorithm of Proposition E.2 based on counting queries to approximate the condition number. In that case, we start with some arbitrary error |E ||. and we keep keep halving until we have a sufficiently good approximation. At each iteration I, the algorithm executes
O (FTTmmx(k))
floating point operations to construct the Krylov matrix using Oppolylog(m, a bits. At most O(log(((K)) iterations are required. The most expensive iteration is the last one, when et x 1/x(K), which gives a total cost of at most
O (q (q]lo((((K))TTmxx(())
floating point operations using at most (ppolylog(m, x(K))) bits. It returns value R (k))),, and succeeds with high probability.
Step 3: Computing a basis. The computation of the basis in Step 4 of Algorithm 8 is arguably the hardest part the analysis of the Block- Block-Krylov iteration. With similar arguments as in the proof of
Theorem A.6, if we write the economy-QU factorizations K -x R3 and K = QR. there exists an orthogonal matrix 4 such that Q Qg 9 E, and || S e, for some e'. This requires
that |K K||
S E' poly(mngk)x(K) which implies a requirement
u S (K))(
for some constant c. Instead of kK) we can use from the previous Step We can now use Theorem A.6 (which computes a basis via QR) on K to approximate Qg by the matrix
Q Og 2 + EQR = (Q+1 + E + EQR = Q 42 + E 42 EQR.
Let Eq E, 2 + EQR. We have that |Eqll |[ 4b EQR I é' EQR- From Theorem
O(log( ( )) bits. Thus, Q is approximate orthogonal basis for the range of the true Block-Krylov matrix K.
A.6 we ensure that || S e' if
To summarize, as long u é' (mg)(s)) which can be achieved by replacing <(K) with ğ above, then can compute Q Q + Eq. where Q is the true orthonormal basis from the QR factorization of K, 4 is a gl x gk orthogonal matrix, and ||Eqll S <'. The total cost is (mm(g))=- floating point operations for Theorem A.6, using O llog( mgg((K) bits.
Step 4: Computing the reduced matrix. In Step 5, Algorithm 8 forms the matrix M to compute its top-k singular vectors. To analyze the computation of M, the first observation is that Q in lines 4- 7 can be replaced by any basis for the column space of K. In particular, replace it by Q. We then perform the multiplication in two steps: M MMX (Q), which returns M1
X EM E u O(poly(gm)) Then we compute M + MM(M, M, M,M, 1 + E where || S u (:ppoll(gkk)) M E u- (ppoly(gkn))
where |EM < u Oppoly(gkmm)
|||
where the last implied if u is sufficiently smaller than 1/ poly(gm). Putting everything together we can write
M 2 + Q XX Q 22 + EM,
where EM contains the errors in Q and also the ones from the multiplication in floating point, and |M || < E[Il E O(e'). e' is as in Step 3. So far we have assumed Ologg( mqq((K) )) bits.
Step 5: Spectral gap and midpoint. Now we have written M = Q Xx Q+++x+++e where + 4 is orthogonal and |EM]l= < O(e') for some é' (0, 1). We can the counting queries in the spirit of Theorem 3.1 to compute the spectral gap and the midpoint of M.
In particular, we start with some desired bound é E (0, 1) for |M by setting the number of bits to Olog( gm)((K) )). At each iteration we repeat Steps 1, 3, and 4 (we do not need to compute
66

---

again the condition number of K in Step 2). After at most m = (log( gap,, (M)) iterations, the error en satisfies <m (gga;;(M) and we obtain two quantities gap, (1 + 8 gap; (M) and
Ne E (MM + & gap; (M) with high probability. In each iteration need to construct M using Steps 3, and 4 with the specified number of bits.
Step 1 costs O(qTmmx(()) operations. Step 3 requires O(m(gk)u-) operations. Step 4 m(gk))=11 operations. The result the ok x qk matrix M on which we call COUNT on each iteration. From Lemma E.I, assuming that we have regular- ized M appropriately using REGULARIZE, COUNT costs at most O((gk) polylog(gk/<m.)) (((a))" polylog(ak/ gap, (M))). The maximum number of bits in the last iteration, which is (gs(()) equal O (pllylog (MP)(M)
executes O(TMMx
Step 6: Principal singular vectors. Given suitable approximation for the gap and the midpoint,
we next use Lemma B. to prove forward error bounds between the true k-spectral projector of M
T T Q xx QΦ and the one of M = M + EM. From Lemma B. if we set μ = jík. we have that
Damio (āg -M]'x eSGN (gan (M)²).
sggn
M jis) sgn (M
|1 S ESGN if |Ev|l E5GN
128
The same holds for the spectral projectors, i.e. | (M) II, (M)l S CSGN- The bound for |E
is achieved by setting u
escweap2
which translates to (log( escoggapp bits.
poly(mgin
It remains to approximate I1, (M), denoted II for simplicity, and then use deflation, similar to tica the vanilla PCA Algorithm 7. We first set CPURIFY 22
) and then use Algorithm 3,
11 < PURIFY (URIFY(M,,L) 1, Pk. gapk EPURIFY).
followed by U < DEFLATE(TT,, k (PURIFY, From Proposition 2.1, PURIFY costs (((gk)" log(ak/ gapx)) floating point operations using
(polylog
)) = Oppolylog(
gl
g&PE
EPCA ge/P
bits. From Theorem A.5, DEFLATE costs (((gk)²) operations using Q(polylog(gk/epca)) bits, and succeeds with high probability. The returned matrix Uμ satisfies ||Ù Unll S EPCA, where U is matrix whose columns form orthonormal basis for the span of the top-k singular vectors
of M. This also implies that Ü UU || < 3CPCA. Then
(M) Ü U < (M) II, [M) + (M) Ü cseN 3CPCA-
We finally return Ze MM( Üz). It holds that where [2 - QU,l S O(poly(mg)). We can then write
Qū, Q EZ 2 T
(MI) + EnJ + Ez
(M) QENQ Q T + Ez
=Q + (M) + EQ) TT)TT QEnQ Ez
= -QIII,MM Q + E' Q + Ez,
where QMI QQ = Z is equivalent to the true, exact arithmetic Z₁ of Algorithm
8. The enror matrices satisfy EEz < Oppo)y(m4),, [|QEn < O(ESGN + EPCA), and [E'| < O[[ M([] E (ESGN gap, (M))).
Putting everything together. We can now state the main result by summarizing the Steps 1-6.
Theorem G.2 (Restatement of Theorem 4.2), Let X be a data matrix X E mm ||||| S 1, & E [n]] a target rank, éPCA E (0, 1) an accuracy parameter, and g =
log(n) Let TMMx (k) denote the
/cpcA
67

---

complexity to stably multiply X or X with dense matrix with &o columns from the right (see Def. G.1). Using the Steps 1-6 that are detailed in Appendix G.3 as a floating point implementation of Algorithm 8, we can compute a matrix Z mm that satisfies
Z ZiZZ 5
with high probability, where Zy is an approximate basis for the top-k principal components of X. returned by Algorithm 8 in exact arithmetic. The total cost is at most
( (k) log(
x(K) gap, (Mi m(gk)) log( (M)) + (gk)* polylog(
ql gap, (M)
floating point operations, using O polyyoog( mgre(K) ePca gap, bits of precision. K, M are the same as in
Alg. 8. The only parameters that are initially required are and EPCA-
Proof. We first compute v. E (K)(K)) in Steps and 2. This costs
O lo (M)))xxx(())))
floating point operations using Oppolylog(qmxz) (K))) bits.
Then, in Step we iteratively use Steps 1, 3, and 4 to compute the midpoint and the gap, in a total of (log(1/ gap,)) iterations. The total cost is
gk o ((TTMmxx()) (k) + m(gk)"=1) log( gap, (M)) (gk)* polylog( gap, (M)
floating point operations using at most
magre (K) (polylog( gaps (M3
bits of precision. The only parameter that we require to know beforehand, except for the matrix sizes, K from the previous step.
In the last Step 6, assuming that CSGN = CPCA- we need to run Steps 1, 3, and 4, using O polylog( mgx(K)) bits. The number of arithmetic operations does exceed the one from EPCA gap," the previous step. Thereafter the costs of DEFLATE, PURIFY, and MM are negligible compared to the costs of Steps 1, 3, and 4. The final matrix Z₁ satisfies
Z.Z. Z S (((PCA)).
where Zμ is the true, exact arithmetic projector in Algorithm 8.
Summarizing everything, compute Z₁ advertised above in total of at most
(TTMMx()) log( <(K) (M) + m(gk) - log( (Mp) + (gk)" polylog( gk gap, (M))
floating point operations, using (polylog( mgg((K) cocA gaps
bits of precision.
68

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
All claims accurately follow the paper contributions.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
Limitations have been clearly discussed in the paper.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The au- thors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the ap- proach. For example, a facial recognition algorithm may perform poorly when image resolution is low images taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to ad- dress problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Ys]]
69

---

Justification: All theoretical results include an informal proof in the main paper and reference the full, rigorous proof in the appendix. All definitions and assumptions are clearly stated.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should clearly stated or referenced in the statement of any theo- rems.
The proofs can either appear in the main paper the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be comple- mented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclu- sions of the paper (regardless of whether the code and data are provided or not)?
Answer: [NA]
Justification:
There are no experiments, only theoretical results.
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is a dataset and/or model, the authors should describe the steps taken make their results reproducible verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, other means that are appropriate the research performed.
While NeurIPS does not require releasing code, the conference does require all sub- missions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results way produce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which thors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing verifying the results.
5. Open access to data and code
70

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [NA]
Justification:
There are no experiments and no code associated with the paper.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines lhttps://nipp) cc/
mublic//guides/CoeeSuuuui sos.eeeessummisisooooooo.o... for more details.
While we encourage the release code and data, we understand that this might not be possible, so "No" is acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: wwtpcc/pppppp..iiiieeee.
Codeeubm ssi onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only a subset of experiments are reproducible, they should state which ones omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [NA]
Justification: The paper does not include experiments.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Question:
Does the paper report error bars suitably and correctly defined or other appropri- ate information about the statistical significance of the experiments?
Answer: [NA]
Justification:
There are no experiments in the paper.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
71

---

The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report I-sigma error bars, but one should state The authors should prefer- ably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [NA]
Justification:
The paper does not include experiments.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics httppp///oorrrs..ccpppp.cahaaaas hctppp/bccctthhcccctueeee.e..eee ines?
Answer: [Yes]
Justification: The research conducted the paper conform, in every respect, with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: There is no apparant societal impact of the work performed as it only includes mathematical results.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact spe- cific groups), privacy considerations, and security considerations.
72

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitiga- tion strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility MI).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
There is no data or models involved in the paper.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by re- quiring that users adhere usage guidelines or restrictions access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Lieenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification:
The paper does not use existing assets.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name the license (c.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
assets are released, the license, copyright information, and terms of use the pack- age should be provided. For popular datasets, paperswi tecodhhodee com/ dattasts has curated licenses for some datasets. Their licensing guide can help determine the li- cense of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
73

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documenta- tion provided alongside the assets?
Answer: [NA]
Justification:
The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the per include the full text instructions given to participants and sereenshots, applicable, well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contri- bution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, cura- tion, or other labor should be paid at least the minimum wage the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equiva- lent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
74


Reviews:
[{"review_id": "EPL4jpdMA7", "reviewer": "Reviewer_ZDvn", "summary": "This paper studies the computational cost of invariant subspaces for generalized eigenvalue problems (GEPs) which is a fundamental computational problem with applications in machine learning and scientific computing. The authors propose a novel method that approximates the spectral projectors of GEPs, and give the first $\\tilde{O}(n^{\\omega})$ bit complexity result for forward-error approximation in the floating point model, improve upon the previous $\\tilde{O}(n^3)$ result. Based on this result, the authors also give new matrix multiplication-type upper bounds for PCA problems in terms of bit complexity.", "questions": "1. In the statement of Theorem 1.1 and Proposition 2.1, is $||H|| \\leq 1$ a proper assumption?\n2. Also in Theorem 1.1, what if $gap_k$ is very small for the given $k$ (say, $gap_k = 2^{-n}$)?\n2. 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?\n3. 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?", "limitations": "The authors have adequately discussed the limitations of their work.", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. This paper is well-written with very solid theoretical results and proofs. The problem, the model, the results are stated clearly. The existing works and previous results are discussed thoroughly.\n2. The idea of reducing the problem to approximating the average and the gap among two adjacent eigenvalues is interesting. The result of approximating these two quantities (Theorem 3.1) can have potential other applications.\n3. This paper gives a new stability analysis of the Cholesky factorization under floating point model and improves upon the previous $\\tilde{O}(n^3)$ result, which is of independent interest (especially in the TCS community).", "weaknesses": "The main concern is that this paper may be too dense with technical details and rigorous proofs as a NeurIPS submission (instead of a TCS conference). It might be better if the authors consider presenting this paper in a way that is easier to follow (e.g., delay more proof details to appendix, while adding some examples of what the results look like for specific schemes of spectral decay, and more detailed discussion about the potential applications / examples of GEPs), so that this problem can be brought to more people in this community.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their positive feedback. Below we elaborate on the specific questions that were raised:\n1. **Question:**\n    In the statement of Theorem 1.1 and Proposition 2.1, is $||H||\\leq 1$ a proper assumption?\n    \n    **Answer:**\n    $||\\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. \n\n2. **Question:**\n    Also in Theorem 1.1, what if \n    $gap_k$ is very small for the given $k$ (say, $k=2^{-n}$)?\n    \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.\n\n3. **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?\n\n    **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.\n\n4. **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?\n\n    **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. \n\n\nThank you again for the feedback and for the interest! We are happy to discuss these points more in depth, if needed.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "yWd8YYmlcZ", "reviewer": "Reviewer_5csW", "summary": "The paper presents a novel approach to approximating invariant subspaces of generalized eigenvalue problems (GEPs), which are fundamental in many applications such as Principal Component Analysis (PCA) and Density Functional Theory (DFT). The authors introduce an algorithm that computes a spectral projector $ \\Pi_k $ with a forward-error approximation guarantee in near matrix multiplication time $ O(n^{\\omega + \\eta}) $, where $ \\omega $ is the matrix multiplication exponent. The approach advances a new analysis for Cholesky factorization and a smoothed analysis for computing spectral gaps, which are key innovations applied to obtain the desired bound with high probability.\n\nThe paper's technical claims are well-supported by rigorous mathematical proofs and thorough analysis. The use of Cholesky factorization and smoothed analysis for spectral gaps is innovative and effectively addresses the computational challenges. The proposed algorithm’s performance is theoretically grounded. However, some sections could benefit from additional explanations, empirical computations and examples to enhance understanding, especially for readers less familiar with the theoretical work developed in the paper, and potentially reaching a larger audience.", "questions": "1) 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?\n2) 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?", "limitations": "The presented paper discusses in great depth the theoretical underpinning of the proposed algorithm, and we understand this already has a significant impact. Nevertheless, the lack of practical implementation, empirical validation, and analysis of some aspects of the algorithms hinders the potentially higher impact of the proposed method. Even a limited empirical section, including investigations into the role of the hyperparameters and covering varied structures of matrices, would highly improve the paper's potential impact and readership.", "rating": 8, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "- **Originality:** The paper introduces a novel approach to solving GEPs and PCA with nearly matrix multiplication time complexity, which is a significant advancement over classical methods. The novel analysis for Cholesky factorization and a smoothed analysis for computing spectral gaps could have a wider application too.\n- **Quality:** The mathematical rigor and comprehensive analysis in terms of bit complexity ensure that the proposed methods are both theoretically sound and practically applicable.\n- **Clarity:** The paper is generally well-written, with detailed explanations of the methods and thorough proofs.\n- **Significance:** The results have broad implications for many applications in machine learning and scientific computing, providing a more efficient computational framework.", "weaknesses": "- **Complexity of Implementation:** The proposed methods, while theoretically sound, may be complex to implement in practice. Detailed implementation guidelines for the algorithms 1, 2, 3, and 4. A standard implementation also would be of value, showing how some of the factors hidden in the $O(.)$ complexity could play a significant role in practice. For example, as it is, it remains open how, for example, the significance of the factor defining the bound on the approximation error $\\epsilon$ (which appears inside a $\\text{log} \\text{log}$ expression) when working with varied sizes of matrices $n$ and approximation error $\\epsilon$. Theorem 1 also should include some comments about this issue.\n- **Hyperparameter Sensitivity:** The algorithm's performance depends on several parameters (e.g., spectral gap, condition number). A more detailed discussion on selecting these parameters and their impact on performance would be helpful, along with empirical implementations and discussion.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "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:\n\n1. **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?\n\n    **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 \\textbf{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. \n2. **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?\n\n    **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.\n\n- **Question:** Complexity of Implementation.\n\n  **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).\n\nTo 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.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "mZhLypqRGU", "reviewer": "Reviewer_FrmC", "summary": "The paper tackles the fundamental problem of GEP subspace approximation (with forward error approximation). The problem is very relevant to different areas of Machine learning. This paper improves the time complexity of this approximation from cubic in n to matrix multiplication time, which is a significant improvement. \n\nI have read the first nine pages and Appendix B, and the idea seems reasonable, and the paper is well explained. However, given the length of the paper (which exceeds 50 pages), I have not been able to verify the correctness of the solution. Therefore, I will not comment on the strengths and weaknesses of the paper and give a low confidence score. \n\nIn my opinion, a more theoretical venue such as STOC/FOCS would be more appropriate for this paper (both because of the emphasis on theory and because the reviewers have significantly lower paper load). I wish the authors the best.", "questions": "1. 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. \n\n2. Is a general dependence on the term gap_k necessary for analysis like this?", "limitations": "The authors have adequately addressed the limitations of this work.", "rating": 5, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "-", "weaknesses": "-", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "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:\n1. **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.\n\n    **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.\n2. **Question:** Is a general dependence on the term gap\\_k necessary for analysis like this?\n\n    **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:\n    $\n    A = \n    [\n        1+\\epsilon, 0;\n        0, 1   \n    ].\n    $\n    The gap is equal to $\\epsilon$, and the spectral projector on the top-1 invariant subspace is $[\n        1, 0; 0, 0].$\n    Now consider the following $\\epsilon$-perturbation of $A$: $\\widetilde A = [1, 0; 0, 1+\\epsilon]$.\n    The spectral projector of the top-1 invariant subspace is completely different than before. In fact, it is orthogonal to the first one.\n\n\nWe hope that these responses clarify the reviewer's concerns. Let us know if the answers are not sufficient so we can elaborate further.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "X0nExWLs2v", "reviewer": "Reviewer_PrHp", "summary": "This paper considers the following optimization problem: given Hermitian $H$ and Hermitian, positive definite $S$, find matrices $C$ and $\\Lambda$ such that $HC=SC\\Lambda$ where $C$ is the eigenvectors and $\\Lambda$ is the eigenvalues. The important application is when the eigenvectors of the interest form an invariant subspace, which covers applications such as DFT and PCA. So the goal is to compute a rank-$k$ projection onto the top or bottom-$k$ eigenvectors. The main contribution of this paper is a stable algorithm that computes an approximate projector $\\tilde \\Pi_k$ such that $\\|\\tilde \\Pi_k - \\Pi_k\\|\\leq \\epsilon$ where $\\Pi_k$ is the projector onto top/bottom-$k$ principal components. The algorithm runs in the current matrix multiplication time up to log factors in $n$, $\\kappa(S)$, $1/\\epsilon$  and the gap between $\\lambda_k, \\lambda_{k+1}$. It uses polylog bits of precision in the above parameters. It also provides error bounds for Cholesky beyond $O(n^3)$ classical algorithm. As corollaries, the algorithm leads to improvement for DFT, PCA and Block-PCA.", "questions": "For your open question on rectangular FMM: is the current runtime bottleneck due to block square FMM? I didn't dig too deep into the algorithm details, but I would imagine the $n^\\omega$ bottleneck is due to certain linear algebraic operations on $n\\times n$ matrices, instead of a bunch of small rectangular FMMs?", "limitations": "Authors have discussed limitations.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper designs a novel algorithm for invariant subspace projection and PCA in the current matrix multiplication time. The algorithm is a good combination of existing techniques (such as approximating the sign function via inverse Newton) and new insights, such as reducing gap and midpoint of eigenvalues computation to eigenvalue threshold counting, and solve the counting via smoothed analysis. \n\nThe paper is well-written and results are presented clearly in the main body, some intuitions and proof sketches are provided for better understanding of the algorithms and analysis. Overall, I think this is a good paper with solid results.", "weaknesses": "The only complaint I have with this paper (might be a bit unfair) is due to the sheer amount of contents and its very numerical linear algebra nature, this paper might be better suited for conferences such as ISSAC or journals such as SIAM journal on matrix analysis. Otherwise, I think this paper is definitely strong enough to be published in NeurIPS.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "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.\n\nIf more information is needed, we would be pleased to further expand on the rectangular FMM or other topics.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "V7pkgHb2on", "reviewer": "Reviewer_KBna", "summary": "The paper analyses complexity of generalized symmetric eigenvalue problem computation.", "questions": "My suggestion is to better emphasize the computational aspect and make the paper more suitable for a conference. For example, \n- Move the main Algorithms to the main body as it gives a more visual impact, which is of higher importance for a conference paper\n- Use tables to compare accuracy and complexity of different algorithms, such as baselines like Lanczos, Randomized SVD, etc.\n- It would be helpful if the authors could show numerical results w.r.t. accuracy, complexity, or relevant quantities to support their theoretical analysis.", "limitations": "Addressed.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- Solid analysis of the computations involved in generalized eigenvalue problems.\n- Excellent survey of related works.\n- Relevant applications in ML (PCA).", "weaknesses": "- The assumption of H being symmetric needs to be in the Abstract as well, as it is quite significant.\n- Missing citation to kernel PCA [1]\n- While the analysis is strong, the computational aspect is not as emphasized. For example, Algorithm 1 is not shown in main body and only in Appendix, and the paper has no experimental evaluation.\n\n[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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "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:\n1. **Question:** The assumption of $H$ being symmetric needs to be in the Abstract as well, as it is quite significant.\n\n    **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.\n\n2. **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. \n\n    **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. \n\n3. **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.\n\n    **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.\n\n4. **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.\n\n    **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.\n\nWe 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.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The reviewers generally agreed that this paper provides strong theoretical results for the central problem of solving generalized eigenvalue problems in finite precision. The paper is well-written, appears sound, and while it is more of a numerical analysis paper than an ML paper, should be of interest to some in the broader NeurIPS community.

Beyond the various helpful suggestions made by the reviewers, I'd like to strongly suggest to the reviewers that in the final version they should tone down/clarify some of the claims about what the work contributes to the special case of the Hermitian eigenvalue problem.

Cor 1.7 of https://arxiv.org/pdf/1912.08805 gives a forward stable solution to this problem in matrix multiplication time. A subtle difference is that it returns a factorization V D V^-1 rather than VD V^T although given that V will be nearly orthogonal due to the accuracy guarantee, I suspect that this is easily resolved. See Remark 6.1 of the same paper for some discussion. 

Thus I feel that the claim in the intro of the submitted paper: "Still, to the best of our knowledge, the state-of-the-art algorithms providing an end-to-end, provably accurate analysis (in the sense of Eq. (2)) for arbitrary invariant subspaces of definite GEPs or even regular Hermitian eigenproblems, in any finite precision model, remain the classic aforementioned O􏰙(n3) algorithms" is either incorrect or at least very misleading/overstated. It should be more carefully discussed how the author's solution improves on or buttons up prior work in the Hermitian case.

I apologize for not brining this issue up during the response period. I don't feel that it detracts from the main contribution of the paper on generalized eigenvalue problems enough to change my recommendation. If it is helpful and possible to discuss with the authors post decision I would be happy to.
Decision: Accept (poster)