A Quasi-Random Approach to Matrix Spectral Analysis
Lior Eldar†

Michael Ben-Or*

arXiv:1505.08126v3 [cs.DS] 6 Apr 2017

April 7, 2017

Abstract
Inspired by the quantum computing algorithms for Linear Algebra problems [HHL09, TS13]
we study how the simulation on a classical computer of this type of “Phase Estimation algorithms” performs when we apply it to solve the Eigen-Problem of Hermitian matrices. The
result is a completely new, efficient and stable, parallel algorithm to compute an approximate spectral decomposition of any Hermitian matrix. The algorithm can be implemented
by Boolean circuits in O(log2 n) parallel time with a total cost of O(nω+1 ) Boolean operations.
This Boolean complexity matches the best known rigorous O(log2 n) parallel time algorithms,
but unlike those algorithms our algorithm is (logarithmically) stable, so further improvements
may lead to practical implementations.
All previous efficient and rigorous approaches to solve the Eigen-Problem use randomization to avoid bad condition as we do too. Our algorithm makes further use of randomization in
a completely new way, taking random powers of a unitary matrix to randomize the phases of
its eigenvalues. Proving that a tiny Gaussian perturbation and a random polynomial power are
sufficient to ensure almost pairwise independence of the phases (mod 2π) is the main technical
contribution of this work. This randomization enables us, given a Hermitian matrix with well
separated eigenvalues, to sample a random eigenvalue and produce an approximate eigenvector in O(log2 n) parallel time and O(nω ) Boolean complexity. We conjecture that further
improvements of our method can provide a stable solution to the full approximate spectral
decomposition problem with complexity similar to the complexity (up to a logarithmic factor)
of sampling a single eigenvector.

1
1.1

Introduction
General

The eigen-problem of Hermitian matrices is the problem of computing the eigenvalues and
eigenvectors of a Hermitian matrix. This problem is ubiquitous in computer science and engineering, and because of its relatively high computational complexity imposes a high computational load on most modern information processing systems.
The Abel-Ruffini theorem implies there is no deterministic expression for the eigenvalues
of a general matrix A, and for this reason eigen-problem algorithms must be iterative. This
gives rise to a host of problems: these algorithms are hard to analyze rigorously, and often turn
out to be unstable, and hard to parallelize.
* School
†

of Computer Science and Engineering, The Hebrew University, Jerusalem, Israel
Center for Theoretical Physics, MIT, Email: leldar@mit.edu

1

As with many other problems in computer science, one typically considers an approximate
spectral decomposition of a matrix. Thus, given a matrix A, we are usually interested not in
its exact eigenvalues / eigenvectors, which may be very hard to compute, (and possibly very
long to describe once computed), but rather in an approximate decomposition:
Definition 1. Approximate Spectral Decomposition - ASD(A, δ)
Let A be some n × n Hermitian matrix. An approximate spectral decomposition of A, with accuracy
parameter δ = 1/poly(n) is a set of unit vectors {vi }ni=1 , kvi k = 1 such that each vi has kAvi −
λi vi k2 ≤ δ, for some number λi ∈ R, and
X

λi vi viT − A ≤ δ kAk ,

i∈[n]

where kXk is the operator norm of X.
By standard arguments the above can be generalized to an arbitrary n × n matrix A, by
considering the Hermitian matrix AH A, in which case ASD(AH A, δ) is an approximation of the
singular vectors (and singular values) of A. We note that the definition of ASD then corresponds
to a ”smooth analysis” of matrices: namely given input A, we do not find a spectral decomposition of A, but rather the decomposition of a matrix A0 , such that kA − A0 k ≤ δ. We also
point out, that the definition of ASD holds just as well in the case of nearly degenerate matrices: we do not require a one-to-one correspondence with the eigenvectors of A, which can be
extremely hard to achieve, but rather to find some set of approximate eigenvectors, such that
the corresponding weighted sum of rank-1 projections form an approximation of A.
When one considers an algorithm A for the ASD problem, one can examine its arithmetic
complexity or boolean complexity. The arithmetic complexity is the minimal size arithmetic circuit C (namely each node computes addition/multiplication/division to unbounded accuracy)
that implements A, whereas the boolean complexity counts the number of boolean AND/OR
gates of fan-in 2 required to implement A.
Given the definition above, and following Demmel et al. [DDH07] we consider an algorithm A to be log-stable (or stable for short), if there exists a circuit C that implements A on
n × n matrices, and a number t = O(log(n)), such that each arithmetic computation in C uses
at most t bits of precision, and output of the circuit deviates from the output of the arithmetic
circuit by at most 1/poly(n). We note that when an algorithm is stable then its boolean complexity is equal to its arithmetic complexity up to a factor O(log(n)). If, however, an algorithm
is unstable then its boolean complexity could be larger by a factor of up to n. In the study of
practical numerical linear algebra algorithms, one usually identifies algorithms that are stable
with ”practical”, and algorithms that are not stable to be impractical. This usually, because
the computing machines are restricted to representing numbers with a number of bits that is a
small fraction of the size of the input.
In terms of parallelism, we will refer to the complexity class NC(k) (see Definition 10) which
is the set of all computational problems that can be solved by uniform Boolean circuits of size
poly(n) in time O(logk (n)). Often, we will refer to the class RNC(k) , in which the parallel NC(k)
circuit is also allowed to accept uniform random bits. One would like an ASD algorithm to
have minimal arithmetic / boolean complexity, and minimal parallel time. Ideally, one would
also like this algorithm to be stable.

1.2

Main Contribution

Inspired by recent quantum computing algorithms [HHL09, TS13], we introduce a new perspective on the problem of computing the ASD that is based on low-discrepancy sequences.
Roughly speaking, low-discrepancy sequences are deterministic sequences which appear to be

2

random, because they ”visit” each small sub-cube the same number of times that a completely
random sequence would up to a small additive error.
Definition 2. Multi-dimensional Discrepancy
s
For integer s, put I s = [0, 1)s . Given a sequence x = (xn )N
n=1 , with xn ∈ I the discrepancy DM (x)
is defined as:
)
(
M
1 X
DM (x) = sup
χB (xn ) − vol(B) ,
M n=1
B∈B
where
Qs χB (xn ) is an indicator function which is 1 if xn ∈ B is the set of all s-products of intervals
i=1 [ui , vi ], with [ui , vi ](mod1) ⊆ [0, 1).
We recast the ASD problem as a question about the discrepancy of a certain sequence related
to the input matrix. Specifically, given a Hermitian matrix A with, say n unique eigenvalues
{λi }i∈[n] the central object of interest is the sequence comprised of n-dimensional vectors of
eigenvalue residuals:
S(A) = ({λ1 · 1}, . . . , {λn · 1}) , ({λ1 · 2}, . . . , {λn · 2}) , . . . ({λ1 · M }, . . . , {λn · M }) ,
where {x} is the fractional part of x ∈ R, and M = poly(n) is some large integer. S(A) is
hence a sequence of length M in [0, 1)n . We would like S(A) to have as small discrepancy as
possible. Hence, in sharp contrast to previous algorithms, instead of the computational effort
being concentrated on revealing ”structure” in the matrix, our algorithm is actually focused on
producing random-behaving dynamics.
The main application of our approach presented in this paper is a new stable and parallel
and stable algorithm for computing the ASD of any Hermitian matrix.
Theorem 3. For any Hermitian matrix 0  A  0.9I, and δ ≤ n−7 we have ASD(A, δ) ∈
RNC(2) , with total boolean complexity O(nω+1 ). The algorithm is log-stable.
The boolean complexity of our algorithm is O(nω+1 ). If however, one is interested in sampling
a uniformly random eigenvector, it can be achieved in complexity O(nω ). 1

1.3

Overview of the Algorithm

To compute the ASD of a given matrix A, we first consider a similar problem of sampling
uniformly an approximate eigenvector of A, where the eigenvalues of A are assumed to be
well-separated. Clearly, if one can sample from this distribution in RNC2 , then by the coupon
collector’s bound concatenating O(nlog(n)) many parallel copies of this routine, one can sample all eigenvectors quickly with high probability. To do this, we require a definition of a
Hermitian matrix that is δ-separated:
Definition 4. δ-separated
Let A be an n × n PSD matrix. A has a complete set of real eigenvalues λ1 > λ2 > . . . > λn ≥ 0. We
say that A is δ-separated if λj − λj+1 ≥ δ for all j < n, and λ1 ≤ 1 − δ.
Next, we introduce the notion of a separating integer w.r.t. a sequence of real numbers:
Definition 5. Separating Integer
Let λ̄ = (λ1 , . . . , λn ) ∈ [0, 1)n . For α > 4 define
Bout = [1 − 1/(4n), −1 + 1/(4n)] and

Bin (α) = [−1/(αn), 1/(αn)],

A positive integer m is said to separate the k-th element of λ̄ w.r.t. Bin , Bout if it satisfies:
1

ω signifies the infimum over all constants c such that one can multiply two matrices in time at most nc .

3

• {mλk } ∈ Bin (α)
• ∀j 6= k {mλj } ∈
/ Bout
and finally define the notion of a separating integer w.r.t. a δ-separated matrix.
Definition 6. A positive integer m is said to separate k in a δ-separated matrix A w.r.t. Bin , Bout , if
m separates the k-th element of L(A) w.r.t. Bin , Bout .
Following is a sketch of the main sampling routine. For complete details see Section 5.
The routine accepts a separating integer m of the i-th eigenvalue of a δ-separated matrix A, a
precision parameter δ and returns a δ approximation of the i-th eigenvector of A:
Algorithm 7. Filter(A, m, δ)
1. Compute parameters:
p = 24n2 dln(1/δ)e, ζ = δ 2 /(2pm).
2. Sample random unit vector:
Sample a standard complex Gaussian vector v, set w0 = v/kvk.
3. Approximate matrix exponent:
Compute a ζ Taylor approximation of eiA , denoted by Ũ .
4. Raise to power:
Compute Ũ m by repeated squaring.
5. Generate matrix polynomial:

p
m
Compute B = I+2Ũ
by repeated squaring.
6. Filter:
Compute w =

B·w0
kB·w0 k .

7. Decide:
Set z = A · w, i0 = arg maxi∈[n] |wi | and compute c = zi0 /wi0 . If
√
kA · w − c · wk ≤ 3δ n
return w, and otherwise reject.

In words - the algorithm samples a random vector and then multiplies it essentially by the
matrix B = ((I + eiAm )/2)p . After this ”filtering” step, it evaluates whether or not the resulting
vector is close to being an eigenvector of A, and keep it if it is. To understand the behavior of
the algorithm, it is insightful to consider the behavior in the eigenbasis of A.
X
w=
αi wi ,
i

where {wi }i∈[n] is an orthonormal basis for A corresponding to eigenvalues {λi }i∈[n] . If {mλi },
i.e. - the fractional part of mλi , is very close to 0, and {mλj } is ∼ 2 ln n/p far from 0 for all j 6= i,
then after multiplication by B and normalization, all eigenvectors wj for j 6= i are attenuated
by factor 1/n2 relative to wi , and hence the resulting vector is 1/n close to an eigenvector of λi .
Hence, a sufficient condition on the number m that would imply that w = F ilter(A, m, δ)
is an approximation of the i-th eigenvector is the following property: {mλi } is very close to 0,

4

and for all j 6= i {mλj } is bounded away from 0. This corresponds to the fact that m separates i
in A, as assumed.
So to sample uniformly an approximate eigenvector, we would like to call F ilter(A, m, δ)
for m ∼ U [M ] for M = poly(n) such that m separates i where i ∼ U [n]. The main observation
here, is that this condition is satisfied if the sequence of residuals of integer multiples of the
eigenvalues S(A) defined above has the aforementioned low discrepancy property.
Most of the work in this study is devoted to achieving this property. Computationally, we
achieve low-discrepancy of S(A) simply by additive Gaussian perturbation prior to calling the
sampling routine. We show that if we perturb a matrix using a Gaussian matrix E of variance
1/poly(n), then S(A + E) has discrepancy which is 1/poly(n). Showing this is non-trivial because arbitrary vectors of eigenvalues λ1 , . . . , λn do not generate low-discrepancy sequences
in general 2 , and on the other hand we are also severely limited in our ability to perturb the
eigenvalues without deviating too much from the original matrix. This is the subject of our
main technical theorem 32, which may be of independent interest:
Theorem. (sketch) Let A be an n × n Hermitian matrix, and E be a standard Gaussian matrix. For
any a > 0, b > 0 there exists M = M (a, b) = poly(n) such that w.p. at least 1 − n−b the sequence of
residuals of eigenvalue multiples of A + n−a · E of length M has discrepancy at most n−b .
Perturbing the input matrix has the additional benefit of making sure that A has a exactly
n unique eigenvalues with high probability. This follows from a breakthrough theorem by
Nguyen, Tao and Vu [NTV16] which has provided a resolution of this long-standing open
problem, which was considered unproven folklore until that point. This theorem allows us to
handle general Hermitian matrices without extra conditions on the conditioning number of A
/ its eigenvalue spacing.

1.4

Prior Art

There are numerous algorithms for computing the ASD of a matrix, relying most prominently
on the QR decomposition [TB97]. For specific types of matrices, like tridiagonal matrices much
faster algorithms are known [Rei05], but here we consider the most general Hermitian case.
We summarize the state of the art algorithms for this problems in terms of their complexity
(boolean / arithmetic, serial / parallel) and compare them to our own:

Csanky
Demmel et
al. [DDH07]
Bini et al.,
Reif [BP92,
Rei05]
New

Arithmetic
Complexity
Õ(nω+1 )
Õ(nω )

Boolean
Complexity
Õ(nω+2 )
Õ(nω )(∗)

Parallel Time

Log-Stable

log2 (n)
N/A

NO
YES

Õ(nω )

Õ(nω+1 )

O(log2 (n))

NO

Õ(nω+1 )

Õ(nω+1 )

log2 (n)

YES

Comments

∗ Conjectured by us, by
modifying the algorithm.
Working with Ω(n) bit Integers

Comparing our algorithm to the best known NC(2) algorithms, it is more efficient by a factor
of n compared with Csanky’s algorithm [Koz92]. Notably, our algorithm is completely disjoint
2
Consider for example the sequence of values 1, 1/2, . . . , 1/n. Then multiplying these numbers individually by a
random integer m and taking the residual would map to 0 all values 1/i for which i|m, and there is a logarithmic
number of these on average. This sequence of eigenvalues is well-separated, and at least potentially could arise from a
random matrix. We show, however, that this is not the typical case.

5

from Csanky’s techniques - which rely on computing explicitly high powers of the input matrix, and computes the characteristic polynomial of the matrix using the Newton identities on
the traces of those powers. This is an inherently unstable algorithm as it finds the eigenvalues by approximating the roots of the characteristic polynomial and small perturbation to the
coefficients of the polynomial may lead to large deviations of the roots.
The algorithms of Demmel et al., Bini et al. and Reif, rely on efficient implementation
of variants of the QR algorithm. Our asymptotic bounds are worse then Demmel et al. in
terms of total arithmetic/boolean complexity, though we conjecture that this is an artifact of
our proof strategy, and not an inherent problem (see the section on open problems), and in
fact, a variant of the algorithm could probably achieve a boolean complexity of O(nω ). We
note that the QR algorithm is not known to be parallelizable in a stable way, and hence the fast
parallel algorithms of Bini et al. and Reif are not stable and probably impractical. In fact the QR
decomposition has been shown, for standard implementations like the Given’s or Householder
method, to be P -complete [LMM99] assuming the real-RAM model. Thus, it is unlikely to be
stably-parallelizable unless P = NC. 3
Thus, to the best our knowledge, our algorithm is the first parallel algorithm for the ASD of
general Hermitian matrices that is both parallel and stable. In particular it achieves the smallest
bit-complexity of any RNC(2) algorithm to date. We conjecture that our approach may present
a practical and parallel alternative to computing the ASD. We dwell on this point a bit more in
Section 8.

1.4.1

Comparison to the power method / QR algorithm

An arguably natural benchmark by which to test the novelty of the proposed algorithm is the
iterative power-method for computing the eigenvalues of a Hermitian matrix. In this method,
one starts from some random vector b0 , and at each iteration k sets:
bk+1 =

Abk
.
kAbk k

This method can achieve polynomially good approximation of the top eigenvalue in time
which is logarithmic in n, for A with a constant spectral-gap.
Both the power method and our proposed scheme are similar in the sense that they attempt
to extract the eigenvectors of the input matrix directly. Also, if two eigenvalues are ε-close in
magnitude, for some ε > 0, then they require essentially the same exponent of A in the power
method, and of e2πiA in our scheme to distinguish between them. However, the similarity
stops here. We maintain, that the power method is both conceptually different, and for general Hermitian matrices performs much worse, in terms of running time, compared with our
proposed algorithm.
Conceptually, in the power method, we seek to leverage the difference in magnitude between
adjacent eigenvalues in order to extract the eigenvectors. On the other hand, in our proposed
scheme we recast the problem on the unit sphere S (1) , where we are interested in the spacing
of the residuals of integer multiples of the eigenvalues. Worded differently, our setting exploits
the additive group structure of the eigenvalues modulo 1, whereas the power method distinguishes between them multiplicatively. In the additive group setting, the advantage is that
we can consider the discrepancy of the sequence of residuals, and analyze how quickly these
residuals mimic a completely independent random distribution. Furthermore, in the additive
setting there is inherent symmetry between the eigenvalues, as no eigenvalue is more likely to
3

We point out that the algorithm of Reif [Rei05] achieves a QR factorization in parallel time O(log2 (n)) in the arithmetic model, thus showing that QR is indeed parallelizable, but it relies on computations modulo large integers and
therefore not stable and not practical.

6

be sampled than another. This allows for a natural parallelization of the algorithm to extract
simultaneously approximation of all eigenvectors.
The well-known QR algorithm for eigendecomposition [GVL96] is the de-facto standard for
computing the ASD, and is, in a sense, a parallel version of the power-method. That algorithm
applies an iterated sequence of of QR decompositions: At each step k we compute (where
A1 = A - the input matrix)
Ak = Qk Rk ,
and then set
Ak+1 = Rk Qk .
3

The algorithm runs in time Õ(n ), by applying several pre-processing steps [GVL96], and the
fast variant of Demmel et al. in time O(nω ). However, as stated above, the QR decomposition
which is at the core of these methods is not known to be stably parallel.

1.5

Open Questions

We outline several open questions that may be interesting to research following this work:
1. Is it possible to attain a serial run-time of O(nω ) for this algorithm? We conjecture that
this is possible based on numerical evidence for a variant of this algorithm, yet we do
not have a proof of this fact. While not directly improving on the best previously known
serial run-time ([DDH07]), it would reduce the overall work performed by our parallel
RNC(2) algorithm to match the work done by state-of-the-art serial run-time algorithm.
2. What other linear-algebra algorithms can be designed using our methods ? We would
like these algorithms to improve on previous algorithms in either the stability, boolean
complexity, parallel run-time, or all these parameters together.
3. Could one reduce the number of random bits required by the algorithm? Currently - we
show that using Õ(n2 ) random bits - i.e. applying additive Gaussian perturbation results
in a matrix whose eigenvalues seed a low-discrepancy sequence. However, can one do
away with only Õ(n) random bits - by applying a tri-diagonal perturbation to the matrix?
4. Is our algorithm practical ? What is the actual run-time of the algorithm on matrices
of ”reasonable” size, and does it compare with state-of-the art parallel algorithms? Our
numerical evidence suggests that in practice our algorithm may run much faster than the
analytical asymptotic bounds we provide here.

7

Table of contents

2

Preliminaries

9

3

Additive Perturbation

10

4

Low-Discrepancy Sequences

12

5

A Filtering Algorithm

16

6

Sampling Separating Integers

20

7

Parallel Algorithm for ASD

26

8

Numerical Experiments

30

8

2

Preliminaries

2.1

Notation

A random variable x distributed according to distribution D is denoted by x ∼ D. For a matrix
X, kXk signifies the operator norm of X. For a set S, U [S] is the uniform distribution on S.
For integer M > 0 the set [M ] is the set of integers {0, 1, . . . , M − 1}. For real number x, {x}
denotes the fractional part of x: {x} = x−bxc. For a Hermitian n×n matrix A, with eigenvalues
{λi }ni=1 , λ1 ≥ λ2 ≥ . . . ≥ λn L(A) = (λ1 , . . . , λn ) ∈ Rn denotes the vector of sorted eigenvalues
of A. For a measurable subset S ⊆ Rn vol(S) denotes the volume of S. Φ is the empty set. GUE
is the global unitary ensemble of random matrices. N, Z, C signify the natural, integer, and
complex numbers, respectively. For a matrix A, AH is the Hermitian conjugate-transpose of A.
For number n > 0 ln n denotes the natural logarithm, and logn denotes the binary logarithm.
µ(η, σ 2 ) is the standard Gaussian measure with mean η and variance σ 2 . U (n) is the set of n × n
unitary matrices.

2.2
2.2.1

Definitions
Hermitian matrices

We repeat again the definition of the ASD due to its importance:
Definition 8. Approximate Spectral Decomposition - ASD(A, δ)
Let A be some n × n δ-separated Hermitian matrix. An approximate spectral decomposition of A, with
accuracy parameter δ = 1/poly(n) is a set of unit vectors {vi }ni=1 , kvi k = 1 such that each vi has
kAvi − λi vi k2 ≤ δ, for some number λi ∈ C, and
X

λi vi viT − A ≤ δ,

i∈[n]

where kXk is the operator norm of X.
By standard arguments the above can be generalized to an arbitrary n × n matrix A, by
considering the Hermitian matrix AH A, in which case ASD(AH A, δ) is an approximation of
the singular vectors (and singular values) of A. We note that if {vi }i∈[n] is an ASD(A, δ) for some
matrix A, and B is a matrix such that kA − Bk ≤ δ, then by the triangle inequality {vi }i∈[n] is
√
an ASD(A, 2δ).

2.2.2

Complexity

Definition 9. Let ω denote the infimum over all t such that any two n × n matrices can be multiplied
using a number of products at most nt .
The current best upper-bound on ω is 2.372 due to Williams [Wil12].
Definition 10. Class NC
The class NC(k) is the set of problems computed by uniform boolean circuits, with a polynomial number
of gates, and depth at most O(logk n).
We will require the following known fact:
Fact 11. [Akl90] There exists an algorithm for sorting n numbers in time log(n), using n processors.
Definition 12. Class RNC
The class RNC(k) is the set of problems that can be computed by uniform boolean circuits, with a polynomial number of gates, accepting a polynomial number of random bits, and depth at most O(logk n).

9

For simplicity, we shall assume in this work that RNC circuits are allowed to accept t-bit
numbers, sampled from a truncated Gaussian distribution, and discretized to t-bits of precision.

2.2.3

Stable Computation

Following Demmel et al. [DDH07] we define the notion of log-stability as one where truncating
each binary arithmetic operation to O(log(n)) bits of precision doesn’t change the result by
much:
Definition 13. (t, δ)-stable randomized computation
Let C denote a randomized arithmetic circuit, and D be its output distribution supported on Rn . Let D
denote the discretization of C to t bits as follows: each infinite-precision arithmetic operation is followed
by rounding to t bits. Let D0 denote the output distribution of D. C is said to be (t, δ)-stable if
∀x ∃y, D(x) = D0 (y) and kx − yk ≤ δ.
Definition 14. Log-stable computation
Let C be a randomized arithmetic circuit that accepts n input numbers. C is said to be log-stable if for
any δ = 1/poly(n) it is (t, δ)-stable for some t = O(log(1/δ)).

3

Additive Perturbation

Matrix perturbation is a well-developed theory [SS90, GVL96] examining the behavior of eigenvalues and eigen-vectors under additive perturbation, usually much smaller compared to the
norm of the original matrix. While general eigenvalue problems are usually unstable against
perturbation, for Hermitian matrices the situation is much better: the Bauer-Fike theorem states
that the perturbed eigenvalues can only deviate from the original eigenvalues by an amount
corresponding to the relative strength of the perturbation.
In particular, when the perturbed matrix A is δ-separated one can compute an explicit estimate for the behavior of the perturbed eigenvalues. We use here a quantitative estimate by
[SS90]:
Fact 15. Rayleigh quotient for well-separated eigenvalues
Let A be a δ-separated n×n Hermitian matrix with eigenvalues λ1 > λ2 > . . . > λn , and corresponding
orthonormal basis {vi }∈[n] . Let E be an additive perturbation of A satisfying |Ei,j | ≤ ε for all i, j. Let
λ̃i denote the i-th eigenvalue of A + E. There exists a constant c > 0 satisfying:
∀i ∈ [n] λ̃i = λi + viH Evi + ζi ,

|ζi | ≤ cε2 /δ.

In fact, if the perturbation E is GUE a stronger characterization is readily available:
Corollary 16. Let A be a δ-separated n × n Hermitian matrix with eigenvalues {λi }i∈[n] , and corresponding orthonormal basis {vi }∈[n] . Let E be GUE. Then the eigenvalues {λ0i }i∈[n] of the perturbed
matrix A0 = A + ε · E are distributed as follows:
∀i ∈ [n] λ̃i = (1 − α) · µ(λi , ε2 ) + α · D + ζi ,

0 ≤ α ≤ 2−n , |ζi | ≤ 16cn · ε2 /δ

for some distribution D.
Proof. By Fact 15 the eigenvalues λ0i behave as
λ0i = λi + viH Evi + ζi ,

10

|ζi | ≤ c max |Ei,j |2 /δ.
i,j

The standard Gaussian satisfies:
√ 
Pµ |x| ≥ 4 n ≤ 2−2n .
√
Thus, by the union bound we have that |Ei,j | ≤ 4 n for all i, j w.p. at least 1 − 2−n . Hence,
w.p. at least 1 − 2−n we have:
∀i ∈ [n] |ζi | ≤ c max |Ei,j |2 /δ ≤ 16cn · ε2 /δ,
i,j

Suppose that this is the case, and let E 0 denote the GUE matrix, conditioned on having bounded
entries. We can write:
E 0 = (1 − α) · E + α · D,
where D is some distribution on n×n matrices and 0 ≤ α ≤ 2−n . The distribution E is invariant
under unitary conjugation, i.e.:
∀u ∈ U (n) uEuH = E
then

∀u ∈ U (n) uE 0 uH = (1 − α) · E + α · D0 ,

where D0 = uDuH is some distribution on n × n matrices. Hence, up to statistical distance at
most α we can assume w.l.o.g. that vi = ei . This implies that
λ0i = λi + (1 − α) · Ei,i + α · D0 + ζi ,
where |ζi | ≤ 16cn · ε2 /δ. Since Ei,i = µ(0, ε2 ) then
λ0i = (1 − α) · µ(λi , ε2 ) + α · D0 + ζi ,

|ζi | ≤ 16cn · ε2 /δ.

The stability of eigenvalues has also been generalized to eigenvectors and even general
invariant subspaces [SS90] in the following sense: if there is a cluster of eigenvalues that
is ”well-separated” from all other eigenvalues, then the orthogonal projection onto the subspace spanned by the corresponding eigenvectors is likewise stable under additive perturbation whose scale is negligible compared to the separation of the eigenvalue of this cluster from
the rest of the spectrum. In particular, if a matrix is δ-separated, then its eigenvectors are individually stable as follows:
Fact 17. [SS90] Rayleigh quotient for well-separated invariant subspaces
Let A be a δ-separated n × n Hermitian matrix with eigenbasis {vi }i∈[n] , and
A1 = A + E,

kEk ≤ ε.

There exists an orthonormal basis {vi0 }i∈[n] and a constant c > 0 satisfying
∀i ∈ [n] kvi − vi0 k ≤ cε2 /δ.
Our interest in additive perturbation, however, is not confined just to ”stability” arguments.
In fact, our main reason for using perturbation is to cause a scattering of the eigenvalues. The
first step of our algorithm in fact applies additive perturbation to provide a minimal spacing
between eigenvalues. Recently Nguyen et al. [NTV16] have provided the first proof that applying additive perturbation to any Hermitian matrix using a so-called Wigner ensemble, an
ensemble of random matrices that generalize GUE, in fact causes the eigenvalues of the perturbed matrix to achieve a minimal inverse polynomial separation. We state their result:

11

Lemma 18. [NTV16] Minimal eigenvalue spacing
Let Mn = Fn + ε · Xn , where Fn is a real symmetric matrix, kFn k2 ≤ 1, ε = n−γ for some constant
γ > 0, and Xn is GUE. Let λ1 ≥ λ2 ≥ . . . ≥ λn denote the eigenvalues of Mn , and put αi = λi − λi+1
for all i < n. Then for any fixed A > 0 there exists B = B(ε) > 0, such that

max P αi ≤ n−B = O(n−A ).
i∈[n]

In particular 4 for any A > 0 there exists B > 0 such that


P min αi ≥ n−B = 1 − O(n−A ).
i∈[n]

Using the lemma above we define:
Definition 19. For any δ = 1/poly(n), let B ∗ (δ) denote the largest number B > 0 such that for every
Fn the matrix Mn = Fn + δXn satisfies:


P min αi ≥ n−B ≥ 0.99
i∈[n]

4

Low-Discrepancy Sequences

4.1

Basic Introduction

Low discrepancy sequences (or ”quasi-random” sequences) are a powerful tool in random
sampling methods. Roughly speaking, these are deterministic sequences that visit any ”reasonable” subset B a number of times that is roughly proportional to the volume of B, up to
some small additive error, called the discrepancy.
Definition 20. Multi-dimensional discrepancy
s
For integer s, put I s = [0, 1)s . Given a sequence x = (xn )N
n=1 , with xn ∈ I the discrepancy DN (x) is
defined as:
)
(
N
1 X
χB (xn ) − vol(B) ,
DN (x) = sup
N n=1
B∈B
where χB (xn ) is an indicator function which is 1 if xn ∈ B and 0 otherwise, and B is a non-empty
family of Lebesgue-measurable subsets of I s .
In this work, we shall define B as the set of all s-products of intervals
s
Y

[ui , vi ],

[ui , vi ](mod1) ⊆ [0, 1).

i=1
∗
Often in literature, one considers instead the star-discrepancy
DN
(x) which is defined by optiQs
mizing over the set of all intervals of the form i=1 [0, ui ].
The definition of discrepancy naturally admits an interpretation in terms of probability: a
sequence x = {xn }n∈[N ] can be interpreted as a random variable x, that takes the value xn ,
as n is the random variable n ∼ U [N ]. Saying that DN (x) ≤ DN means that for any B ∈ B
the probability that xn is contained in B is equal to vol(B) + ε, where |ε| ≤ DN . In this work,
deviating somewhat from standard terminology, we shall often refer to the discrepancy of a
distribution x on length N , s-dimensional sequences. In this case DN (x) will denote

DN (x) = sup Px0 ∼x,n∼U [N ] (x0n ∈ B) − vol(B) ,
B∈B

4

applying the union bound over all eigenvalues

12

Low-discrepancy sequences have much in common with random sampling, or the MonteCarlo method, in the sense that they visit each cube a number of time that is roughly proportional to its volume, up to a small additive error. Yet, contrary to the Monte-Carlo method, such
sequences are not random, but only appear to be random in the sense above. Arguably such
sequences are most useful in the context of numerical integration: instead of computing an integral of a continuous function f over the s-dimensional cube [0, 1)s one can replace it with the
average value of the function over the points of a low-discrepancy sequence x = {xi }N
i=1 . The
well-known Koksma-Hlawka theorem, then connects the approximation error as a function of
the discrepancy DN (x) as follows:
Z
N
1 X
∗
f (xi ) −
(x),
f (z)dz ≤ V (f ) · DN
N i=1
[0,1)s
where V (f ) is the bounded-variation of f on [0, 1)s .
There are deterministic s-dimensional sequences x = {xi }N
i=1 with discrepancy as low as
DN (x) ≤ C ·

logs N
,
N

and matching lower-bounds (up to constant factors) on the smallest possible discrepancy are
known for s = 1 [Nie92]. Hence, usually one considers low-discrepancy sequences that are
very long (N ) compared to the dimension (s).
At this point, it may be insightful to consider an example of a low-discrepancy sequence:
the well-known van-der Corput sequence [Nie92]: consider the binary expansion of a positive
integer
b
X
b
∀n ∈ [N = 2 ] n =
αi 2i ,
i=0

then the van-der Corput sequence x = {xn }N
n=1 is defined as:
b

∀n ∈ [N = 2 ] xn =

b
X

αi 2−i−1 .

i=0

This sequence has
logN
.
N
We note that the discrepancy upper-bound decays asymptotically like O(1/N ) (assuming
small√dimension s) whereas the Monte-Carlo method converges more slowly, behaving as
O(1/ N ) - and hence these sequences are often preferred as a method of numerical integration
to Monte-Carlo. They are also advantageous compared with purely deterministic methods, like
defining a fine-resolution grid, because usually one can increase the length of the quasi-random
sequence, and reduce its discrepancy while making use of all previous points of the sequence.
∗
DN
(x) ≤ C

4.2

Some basic facts

We require a Lemma [2.5] due to Niederreiter [Nie92].
Lemma 21. [Nie92] Small point-wise distance implies similar discrepancy
Let x1 , . . . , xN , y1 , . . . , yN denote two s-dimensional sequences for which |xn,i − yn,i | ≤ ε, for all
n ∈ [N ], i ∈ [s]. Then the discrepancies of these sequences are related by:
|DN (x1 , . . . , xN ) − DN (y1 , . . . , yN )| ≤ s · ε.

13

(1)

We prove an additional fact:
Fact 22. Let x = {xn }n∈[N ] be a distribution on sequences with discrepancy at most DN (x), and let
y = {yn }n∈[N ] denote the following sequence:
yn = xn + zn ,
where z = {zn }n∈[N ] is some sequence chosen independently from x. Then
DN (y) = DN (x).
Proof. For each S ∈ B we have
Z
Py,n∼U [N ] (yn ∈ S) = Px,z,n∼U [N ] (xn +zn ∈ S) =

Px,n∼U [N ] (xn ∈ S−zn )·Pz,n∼U [N ] (zn )dzn
[0,1)s

By the discrepancy assumption on the sequence x we have that
∀S ∈ B

Px,n∼U [N ] (xn ∈ S) − vol(S) ≤ DN (x).

Therefore by the above
Z
Py,n∼U [N ] (yn ∈ S) =

(vol(S) + ε(zn )) · Pz,n (zn )dzn , |ε(zn )| ≤ DN (x).
[0,1)s

and hence by convexity:
Py,n∼U [N ] (yn ∈ S) = vol(S) + ε,

|ε| ≤ DN (x),

which implies the claim by the definition of discrepancy.

4.3

The Good Seed Problem

In our context we will be interested in sequences x = {xn }N
n=1 of the form
n gn o
xn =
,
N
where g ∈ [N ]s is some s-dimensional vector, called the seed of the sequence. Specifically, the
vector g/N will be the vector of eigenvalues of an n × n Hermitian matrix A whose spectrum
g = L(A) we would like to analyze. Since it is unreasonable to assume that the input matrix has
a spectrum that is a good seed, we would like to find a perturbation of the matrix A0 = A + E
such that g 0 = L(A0 ) has a corresponding sequence, defined as above, with low-discrepancy.
Niederreiter has shown [Nie92] that if g is sampled uniformly on [N ]s then it is a good seed
with high probability:
Lemma 23. Let s, N be an integers and g ∼ U ([N ]s ). Then


logs N
P DN (x) ≤
≥ 1 − 1/N.
N
For our application we require that N = poly(n), and s = 2, in which case the above
discrepancy is sufficiently low for our purposes. Yet, since it requires the normalized seed g/N
to be essentially uniform on [0, 1)n , it implies that the corresponding matrix perturbation E
added to A must be very strong - thereby loosing all connection to the input matrix.

14

4.4

Finding Reasonably-Good Seeds Locally

To bridge the gap between weak-perturbation and low-discrepancy we show a new lemma,
which may be of independent interest: it allows to trade-off the extent to which g is random,
and the discrepancy of the sequence generated by g. Specifically,
√ we will show that if g/N is
uniform on cubes of much√smaller side-length, i.e. at least 1/ N , then the resulting sequence
has discrepancy O(logs N/ N ). This is the subject of the following lemma:
Lemma 24. We are given integer N , with prime divisor M = Θ(N a ) for some constant a > 0, and an
integer s. Let g = (g1 , . . . , gs ) ∈ N s , such that each coordinate gi is independently chosen uniformly
on some interval Ii ⊆ [N ] of size M . Let x = x(g) = {xn }N
n=1 be the following s-dimensional sequence
of length N corresponding to residuals of g:
ng · no
.
xn =
N
Then


√
√ 
Pg DN (x) ≤ 2logs (M )/ M ≥ 1 − 1/ M .

Proof. For integers P, s put Cs∗ (P ) as the set of all vectors in Zs with entries in [−P/2, P/2) ∩ Z,
excluding the all-zero vector. Following Niederreiter [Nie92] we define for each h ∈ Zs
r(h) =

s
Y

(2)

r(hi ),

i=1

where r(hi ) = max(1, |hi |). For g = (g1 , . . . , gs ) ∈ Zs , we denote:
X
R(g, P ) =
r(h)−1 .

(3)

h·g=0(modP ),h∈Cs∗ (P )

By Theorem [5.10] of [Nie92] when each gi is randomly chosen on the entire interval [P ], for
prime P , then
 s

log (P )
Eg [R(g, P )] = O
.
(4)
P
Since R(g, P ) ≥ 0 for all vectors g, then


√
logs (P )
≤ 1/ P .
Pg R(g, P ) ≥ √
P
or


Pg

logs (P )
R(g, P ) ≤ √
P



(5)

√
≥ 1 − 1/ P .

(6)

Let us use the above equation to upper-bound the discrepancy of S(g). Recall that M is a prime
divisor of N , with M = Θ(N a ). We first observe that:
X
R(g, N ) =
r(h)−1 ≤
(7)
h·g=0(modN ),h∈Cs∗ (N )

X

X

r(h)−1 +

h·g=0(modM ),h∈Cs∗ (M )

h·g=0(modM ),h∈Cs∗ (N ),∃i,

r(h)−1
s.t.

(8)

|hi |≥M

We note that the second term is at most
s
M

X

r(h)−1 =

h∈Zs−1

15

s · logs−1 (N )
M

(9)

Regarding the first term, h · g = 0(modM ) if and only if h · (g(modM )) = 0(modM ), and since
each gi is uniform on some interval of size M , then gi (modM ) is uniform on [M ]. Since M is
prime, we can apply equation (6) to the first term with P = M . Using this equation and the
upper-bound on the second term we have:


√
2logs (M )
≥ 1 − 1/ M .
(10)
Pg R(g, N ) ≤ √
M
N −1

According to Theorem 5.6 of [Nie92], the discrepancy of the sequence {gn/N }n=0 is upperbounded by:
s
DN (S(g)) ≤
+ 2−s · R(g, N ).
(11)
N
Plugging equation (10) into equation (11) implies:


√
2logs (M )
(12)
≥ 1 − 1/ M .
Pg DN (S(g)) ≤ √
M

5

A Filtering Algorithm

In this section we provide the specification of the filtering algorithm, which is the main computational black box of our algorithm. This algorithm accepts an integer m that separates the i-th
eigenvalue of a Hermitian matrix A and computes an approximation for the i-th eigenvector,
with high probability:
Algorithm 25. Filter(A, m, δ)
1. Compute parameters:
p = 24n2 dln(1/δ)e, ζ = δ 2 /(2pm).
2. Sample random unit vector:
Sample a standard complex Gaussian vector v, set w0 = v/kvk.
3. Approximate matrix exponent:
Compute a ζ Taylor-series approximation of eiA , denoted by Ũ .
4. Raise to power:
Compute Ũ m by repeated squaring.
5. Generate matrix polynomial:

p
m
Compute B = I+2Ũ
by repeated squaring.
6. Filter:
Compute w =

B·w0
kB·w0 k .

7. Decide:
Set z = A · w, i0 = arg maxi∈[n] |wi | and compute c = zi0 /wi0 . If
√
kA · w − c · wk ≤ 3δ n
return w, and otherwise reject.

16

We now show that if the algorithm is provided with an integer m that separates the k-th eigenvalue of A in the sense defined in Definition 6, then the output is close to the k-th eigenvector
of A.
p
Theorem 26. Let δ ≤ n−10 and α = ln(1/δ). We are given an n × n Hermitian matrix A with
eigenvalues {λi }i∈[n] and a corresponding orthonormal eigenbasis {vi }i∈[n] . Additionally, we are
provided an integer m that separates k in A, w.r.t. Bin (α), Bout , in the sense of Definition 5. Let
w = F ilter(A, m, δ). Then
P (kw − vk k ≤ δ) ≥ 1 − 3n−3 .
The algorithm has boolean complexity at most 2cnω · log(2p2 m2 /δ 2 ), and runs in parallel time
O(log2 (n)).
Proof. Let {τ` }`∈[n] denote the set of eigenvalues of Ũ . Since Ũ is a polynomial in A (truncated
Taylor series) then {v` }`∈[n] is also an orthonormal basis for Ũ . Since in addition kŨ − eiA k ≤ ζ
then
∀` ∈ [n] τ` − eiλ` ≤ ζ.

(13)

Let w0 = B · w0 and denote
X

w0 =

β` v` ,

and

w0 =

`∈[n]

X

α` v` .

`∈[n]

Since A, Ũ share the same basis of eigenvectors, then by the definition of the matrix B the
coefficients α` , β` are related by:
|α` |2 = |β` |2 ·

1 + τ` m
2

2p

.

So by Equation 13
1 + eimλ`
|α` |2
≥
|β` |2
2

2p

− 2pmζ.

Since m separates k then {mλk } ∈ Bin , and for all ` 6= k we have {mλ` } ∈
/ Bout . Thus, for ` = k:


p
|αk |2
≥
1
+
cos(1/2n
ln(1/δ))
|βk |2

p

− 2pmζ

Using Claim 30

≥ 1−

1
2
4n ln(1/δ)

p
− 2pmζ ≥

1
.
2e6

On the other hand, for all ` 6= k we have:
|α` |2
1 + eimλ`
≤
2
|β` |
2
so since m separates k then:

2p

+ 2pmζ.

p

≤ |(1 + cos(1/2n))| + 2pmζ

17

(14)

Using Claim 30
≤ (1 − 1/12n2 )24n

2

ln(1/δ)

+ 2pmζ ≤ e−2 ln(1/δ) + 2pmζ ≤ 2δ 2 .

(15)

By Fact 29 for any ε = 1/poly(n) there exists a constant c > 0 such that
p
P(∀i, j |βj | ≤ c|βi | ln(1/ε)/ε) ≥ 1 − 3nε.
Choose ε = n−4 . Then by Equations 14 and 15:


|α` |2
2
2
6
8
≤ c (2δ ) · (2e ) · 4 ln n · n ≥ 1 − 3n−3 .
P ∀` 6= k
|αk |2
and so for δ ≤ n−10 there exists η ∈ C, |η| = 1 such that
w0
− η · vk
kw0 k

2

≤

1 X
|αj |2 ≤ 16c2 n9 ln nδ 2 e6 < δ.
|αk |2
j6=k

for sufficiently large n. Using Claim 28 we conclude that w.p. at least 1 − 3n−3 over choices
w0 , the criterion is met and the algorithm returns a vector w = w0 /kw0 k satisfying the equation
above.
Arithmetic run-time: The approximation of eiA by Ũ requires, using Fact 31 a time at most
cnω log(1/ζ) = cnω · log(2pm/δ 2 ).
Next, the repeated powering of Ũ to a power m requires time at most: cnω dlog(m)e and the
repeated powering of B to the power p requires time at most: cnω dlog(p)e Hence the total
complexity is : 4cnω · log(pm/δ 2 ).
Depth complexity: Each matrix product can be carried out in depth log(n). Each of steps 3 to 6
involves at most log(m) + log(p) sequential matrix multiplications. The sorting algorithm can
be computed in depth O(log(n)) by Fact 11. Hence the depth complexity of the entire circuit is
at most log(n) · (log(m) + log(p)) + O(log(n)) = O(log2 (n)).
We conclude the proof of the theorem by showing stability:
Claim 27. Under the assumption of Theorem 26 the algorithm is log-stable.
Proof: Consider the arithmetic operations involved in computing the filtering algorithm:
1. Generating an approximation Ũ of eiAm as a truncated Taylor series.
2. Raising Ũ to a power m ∈ [M ].
3. Computing ((I + Ũ )/2)p .
4. Normalizing Bw0 /kBw0 k.
Consider an arithmetic circuit C implementing the above, and the circuit D = D(C, t) the discretization of C to t bits of precision modeled as follows: after each arithmetic step,
the result is rounded to the nearest value of 2−t . Consider all steps except division. A is δseparated so in particular kAk ≤ 1. Thus, whenever we multiply two matrices at any of the
steps above
√ both have norm at most 1. Hence, at each rounding step the error is increased by
at most n2−t . Finally, considering the final division step, we observe that since m separates
k,
14 we have kBw0 k ≥ 1/(2e6 ). This implies that the total error is at most
√ then by Equation
−t
6
n(p + M ) · 2 · 2e . Since M, p are both polynomial in n then for any δ = 1/poly(n) the error
is at most δ for some t = O(log(1/δ)).

18

5.1

Supporting Claims

Claim 28. Suppose that kw − vk ≤ δ for some unit eigenvector v of A, and δ ≤ 1/4. Let z = A · w,
and i0 denote i0 = arg maxi∈[n] |wi |. Let c = zi0 /wi0 . Then
√
kA · w − c · wk ≤ 3δ n.
Proof. Write w = v + E, kEk ≤ δ. Since kAk ≤ 1 then z = Aw = λv · v + E 0 , where kE 0 k ≤ δ.
Therefore
zi0 = λv · vi0 + Ei00 = λv · (wi0 − Ei0 ) + Ei00 ,
So
√
Since |wi0 | ≥ 1/ n then:

zi0 = λv wi0 + E 00 , |E 00 | ≤ 2δ.
c=

√
zi0
= λv + ζ, |ζ| ≤ 2δ n.
wi0

Hence
kA · w − c · wk ≤ kλv · v + E 0 − (λv + ζ) · (v + E)k = kE 0 − (λv + ζ)E − ζvk
which by the triangle inequality is at most
√
√
√
δ + (1 + 2δ n)δ + 2δ n ≤ 3δ n,
for sufficiently large n.
Fact 29. Random unit vectors have well-balanced entries
Let {vi }i∈[n] be some orthonormal basis of Cn , 0 < ε = 1/poly(n), and v ∈ Cn a uniformly random
complex unit vector. For any i ∈ [n] let αi = |hv, vi i|. For any ε = 1/poly(n) there exists a number
c1 > 0 independent of n, such that


p
P ∀i, j |αi |/|αj | ≤ c1 ln(1/ε)/ε ≥ 1 − 3nε.
Proof. Sample a random unit vector by sampling a standard Gaussian vector and normalizing
to unity. Let x = (x1 , . . . , xn ) where xi ∼ µ(0, 1), and all are i.i.d. Since the Gaussian measure
is invariant under conjugation by a unitary matrix we can assume w.l.o.g. that vi = ei for each
i, in which case the αi ’s are i.i.d. αi ∼ µ(0, 1). By the Gaussian measure for every c2 > 0 there
exists c1 > 0 such that
√
P(|αi | ≥ c1 · ln n) ≤ n−c2 .
Therefore, by the union bound
√
P(∀i |αi | ≥ c1 · ln n) ≤ n−c2 +1 .
On the other hand, for any αi , 0 < ε < 1/n, we have
P(|αi | ≤ ε) ≤ 2ε.
So by the union bound
P(∀i ∈ [n] |αi | ≤ ε) ≤ 2nε.
Set c2 = lnn (1/ε). Since ε = n−k for some constant k, then c2 = k is also a constant. Then there
exists a constant c1 such that:


√
P ∀i, j |αi |/|αj | ≥ c1 ln n/ε ≤ 2nε + n−c2 +1 = 3nε,
Since normalization does not change these ratios, then the maximal ratio between the components of x, is the same as the maximal ratio between the components of x/ kxk. This implies
the proof.

19

Claim 30. ∀θ ∈ [−0.01, 0.01]

1−

θ2
2

≤

1+cos(θ)
2

≤1−

θ2
3 .

Proof. Follows by truncating the Taylor series of cos(x) to second order.
Fact 31. Efficient approximation of exponentiated matrix
Given a Hermitian n × n matrix A, kAk ≤ 1, and error parameter ε > 0, a Taylor approximation of
eiA , denoted by ŨA can be computed in time cnω log(1/ε) and satisfies eiA − ŨA ≤ ε.
Proof. Put s = log(1/ε) and consider the Taylor approximation of eiA up to the first s terms:
ŨA :=

s−1
X
(iA)m
m!
m=0

Then the approximation error can be bounded as:
eiA − ŨA ≤

∞
∞
∞
m
X
X
X
kAm k
kAk
1
≤
≤
≤ 2−s ≤ ε,
m!
m!
m!
m=s
m=s
m=s

(16)

where we have used the fact that kAk ≤ 1. The complexity of the approximation is comprised
of log(1/ε) matrix products for a total of cnω log(1/ε).

6

Sampling Separating Integers

In this section we show our main technical tool: which is that perturbing a δ-separated Hermitian matrix A by a Gaussian matrix of a carefully calibrated variance, results in a corresponding
sequence of residuals S(A) having low-discrepancy, at least for pair-wise variables, which in
turn implies that we can separate each eigenvalue of A almost uniformly:
Theorem 32. Let A be a δ-separated n × n PSD matrix, E GUE, ζ ≤ min{δ 13 , n−50 }, and
4 < α ≤ n. For any M ≥ ζ −1.6 we have:
∀k ∈ [n] PE,m∼U [M ] ( m separates k in A + ζ · E w.r.t. Bin (α), Bout ) ≥ 1/(5αn)

6.1

Additive Perturbation

Definition 33. (σ, ε)-normal vector
Let v = (v1 , . . . , vn ) be a vector of n random variables. v is said to be (σ, ε)-normal, if each vi ∼ xi + yi
with xi ∼ µ(λi , σ 2 ) for some λi ∈ R, x1 , . . . , xn are independent, and |yi | ≤ σε, for all i.
By our definitions above, Gaussian perturbation of a matrix with well-separated eigenvalues results in a (σ, ε)-normal vector as follows:
Fact 34. Perturbation of well-separated matrices
Let A be an n × n αmin -separated Hermitian matrix with eigenvalues λ1 ≥ λ2 . . . ≥ λn , where
αmin ≥ ε, and ε > 0 is some constant. Let E be GUE, and A0 = A + εL · E, where L ≥ 2. Then w.p. at
least 1 − n · 2−n the vector of eigenvalues of A0 (λ01 , . . . , λ0n ) is (εL , cnεL−1 )-normal, for some constant
c > 0.

20

Proof. We invoke Corollary 16 choosing ε as εL and δ as ε. We get:
∀i ∈ [n] λ0i = (1 − α) · µ(λi , ε2L ) + α · D + ζi ,

|ζi | ≤ 16cn · ε2L−1 , 0 ≤ α ≤ 2−n .

Choosing σ = εL implies the each λ0i w.p. at least 1 − 2−n is sampled according to:
λ0i = µ(λi , ε2L ) + ζi ,

|ζi | ≤ σ · 16cnεL−1 .

Taking the union bound over all i ∈ [n] then implies the proof.

6.2

Moderately Low-Discrepancy

Lemma 35. Low-discrepancy sequence from almost normal vectors
Let B > 0, and v = (v1 , . . . , vn ) be some (σ, ε)-normal vector, for σ = n−B , ε ≤ n−0.9B . There exists
M ≤ n1.6B such that for any S = {i1 , . . . , is } ⊆ [n] , |S| = s the distribution on s-dimensional
sequence of length M :
Vs ≡ {({m · vi1 } , . . . , {m · vis })}m∈[M ]
satisfies
DM (Vs ) ≤ 4logs (n) · n−0.1B .
Proof. Let P be the minimal prime at least n0.3B , and put M = P 5 . By Bertrand’s postulate, for
sufficiently large n we have that M = P 5 ≤ n1.51B ≤ n1.6B . For any z ∈ [0, 1) let z M be the
number closest to z in the grid m/M , m ∈ [M ].

Removal of non-independent component. Since v is (σ, ε)-normal then vi = Xi + Yi ,
where Xi ∼ (ηi , σ 2 ), |Yi | ≤ εσ, and the Xi ’s are independent. Let VSX denote the sequence
generated by taking only the X component of the seed vector v, i.e.:
VSX ≡ {({m · Xi1 } , . . . , {m · Xis })}m∈[M ]
Fact 36.

(17)

DM (VS ) ≤ DM (VSX ) + s · n−0.2B

Proof. Consider the r.v.’s Xi , Yi . By our assumption
∀i ∈ [n] |Yi | ≤ σε = n−1.9B .
Thus:
∀m ∈ [M ], i ∈ [n]

|{mvi } − {mXi }| ≤ m · n−1.9B ≤ M n−1.9B ≤ n−0.3B

(18)
(19)

By Lemma 21, we can conclude that the discrepancy of our target sequence VS follows tightly
the discrepancy of VSX :
DM (VS ) ≤ DM (VSX ) + s · n−0.3B
(20)

Reducing Gaussian measure to uniform measure Consider the vector derived by truncating each coordinate of the vector (Xi1 , . . . , Xis ) to the nearest point on the M -grid:
X M = (XiM
, . . . , XiM
) = (bM Xi1 e/M, . . . , bM Xis e/M ).
1
s
Consider the discrepancy of the distribution on s-dimensional sequences formed by taking
integer multiples of X M . We claim:

21

Fact 37.



Pv DM (VSX,M ) ≤ logs (n) · n−0.1B ≥ 1 − 3n−0.1B ,

Proof. In Fact 40 choose as parameter m = n0.2B+2 . We get that w.p. at least 1 − 2n2 /m =
1 − 2n−0.2B each Xi samples a convex mixture of variables {wj }j∈[m] where
wj ∼ U (Ij ), |Ij | = σ/m = n−1.2B−2

(21)

Hence, w.p. at least 1 − 2n−0.2B for all i ∈ [n], the variable M · {XiM } is a convex mixture of
uniform random variables on intervals Ij ⊆ [M ], where
|Ij | ≥

σM
≥ n1.5B · n−1.2B−2 ≥ M 0.2 .
m

(22)

We apply Lemma 24 to the sequence of residuals of integer multiples, with the seed X M :
VSX,M ≡





mX1M , . . . , mXsM m∈[M ] .

(23)

The lemma requires that each variable be distributed as:
M XiM ∼ U [I],
where I is some interval of [M ], for integer M > 1 satisfying:
|I| ≥ P, P prime , P ≥ N a , a > 0.
By our choice of parameters M has a prime divisor P equal to M 0.2 = P . Hence, by Equation
22 we can satisfy the assumption of the lemma by choosing the parameters N, M, a for Lemma
24 as follows:
N = M, M = P, a = 0.2.

(24)

Hence, by Lemma 24, and accounting for the Gaussian-to-uniform approximation error we get:


Pv DM (VSX,M ) ≤ 2logs (n) · n−0.1B ≥ 1 − n−0.1B − 2n−0.2B ≥ 1 − 3n−0.1B .
(25)

Treating the residual w.r.t. the M -grid Define: the truncation error
∀i ∈ [s] ri := Xi − XiM .
In Fact 36 we analyzed the error Yi whose magnitude is negligible even w.r.t. 1/M , and can
thus be disregarded for any element of the sequence VS . Unlike this, the residual error ri cannot
be disregarded because when multiplied by integers uniformly in [M ] it assumes magnitude
Ω(1). Thus, it requires a different treatment.
Corollary 38.


Pv DM (VSX ) ≤ 2logs (n) · n−0.1B ≥ 1 − 4n−0.1B

Proof. Express the i-th element of the sequence using ri :

∀i ∈ [s] {Xi · m} = (XiM + ri ) · m = {{mXiM } + {mri }}

22

(26)

Let E denote the event in which Xi is sampled according to wj ∼ U [Ij ] where wj is at distance
at least 1/M from either one of the edges of Ij . Conditioned on E, the random variables ri and
XiM are independent for all i ∈ [s] so conditioned on E we also have
∀i ∈ [s] {mXiM }, {mri } are independent.
VSX,M is a distribution on sequences formed by sampling the initial seed {XiM }, and VSX is a
distribution on sequences formed by adding to VSX,M a random sequence formed by sampling
the independently seed {ri }i∈[s] , and then generating the length-M sequence
{(mr1 , . . . , mrs )}m∈[M ] .
Then since conditioned on E {mXiM }, {mri } are independent for all i, then by Fact 22 we have:
DM (VSX |E) = DM (VSX,M )
and so by Fact 37


Pv DM (VSX |E) ≤ logs (n) · n−0.1B ≥ 1 − 3n−0.1B ,

(27)

By Equation 22 the probability of E is at least:
Pv (E) ≥ 1 − |Ij |/(2M ) ≥ 1 − M 0.2 /(2M ) ≥ 1 − n−B .
Thus:


Pv DM (VSX ) ≤ logs (n) · n−0.1B ≥ 1 − 3n−0.1B − P(E) ≥ 1 − 4n−0.1B .

(28)

Conclusion of proof:
By Corollary 38 we have

Pv DM (VSX ) ≤ 2logs (n) · n−0.1B ≥ 1 − 4n−0.1B
and by Fact 36 we have

DM (VS ) ≤ DM (VSX ) + s · n−0.2B

Thus by the union bound:

Pv DM (VS ) ≤ 2logs (n) · n−0.1B + s · n−0.2B ≥ 1 − 4n−0.1B
thus:


Pv DM (VS ) ≤ 3logs (n) · n−0.1B ≥ 1 − 4n−0.1B

Hence for all but a measure 4n−0.1B of sampled vectors v, the resulting sequence has discrepancy at most 3logs (n)n−0.1B . Since the discrepancy measures the worst-case additive error for
any set this implies that:
DM (VS ) ≤ 3logs (n)n−0.1B + 4n−0.1B ≤ 4logs (n)n−0.1B

23

6.3

Approximate Pairwise Independence

Lemma 39. Let λ̄ = (λ1 , . . . , λn ) ∈ [0, 1)n and M a positive integer that satisfy:

∀i 6= j DM {(mλi , mλj )}m∈[M ] ≤ ζ, ζ ≤ n−4
Let 4 < α ≤ n. For each k ∈ [n] w.p. at least 1/(5αn) over choices of m ∼ U [M ] the sampled sequence
m separates k w.r.t. Bin (α), Bout .
Proof. Let Ei denote the following event:
Ei := xi ∈ Bin ∧ ∀j 6= i xj ∈
/ Bout
We want to show that
∀i ∈ [n] P(Ei ) ≥

1
.
5αn

Let s denote the number of xj ’s not in Bout :
s = |{j j 6= i xj ∈
/ Bout }|
Then under this notation we have:
P(Ei ) = P(s = n − 1 ∧ xi ∈ Bin ).

(29)

Consider the conditional expectation:
E [s|xi ∈ Bin ]
By linearity of expectation:
E [s|xi ∈ Bin ] =

X

P [xj ∈
/ Bout |xi ∈ Bin ] .

j6=i

Considering each summand separately:
P(xj ∈
/ Bout |xi ∈ Bin ) =

P(xj ∈
/ Bout ∧ xi ∈ Bin )
P(xi ∈ Bin )

Using the pairwise discrepancy assumption, the above is at most:
(1 − |Bout |) · |Bin | + ζ
1
0.51
≤ 1 − |Bout | + 2ζαn =
+ 2αζn ≤
|Bin | − ζ
2n
n
and so by Equation 30
E [s|xi ∈ Bin ] = (n − 1) · P [xj ∈
/ Bout |xi ∈ Bin ] ≤ 0.51.
The variable s|xi ∈ Bin accepts only integral values, and by Markov’s inequality:
P(s ≥ 1|xi ∈ Bin ) ≤ 0.51
Therefore
P(s = 0|xi ∈ Bin ) ≥ 0.49.
Using again the 1-dimensional discrepancy we have
1
1
−ζ ≥
.
αn
2αn
Substituting the last two inequalities into Equation 29 yields:
P(xi ∈ Bin ) ≥

P(Ei ) ≥ 0.49 ·

24

1
1
≥
2αn
5αn

(30)

6.4

Proof of Theorem 32

Proof. By assumption A is δ-separated and ζ ≤ min{n−50 , δ 13 }. Consider the perturbed matrix
A0 = A + ζE.
Choose L = 13 and ε = ζ 1/13 . Then by Fact 34 w.p. at least 1 − n2−n all eigenvalues of A0 are
(σ, ε)-normal with parameters
σ = ζ ≤ n−50 , ε ≤ 16cnζ 12/13 ≤ ζ 0.9 = σ 0.9
where the last inequality follows because ζ ≤ n−50 . We assume that this is the case and account
for the negligible error at the end. Set now B = logn (1/σ). Then the eigenvalues of A0 are (σ, ε)normal for σ = n−B an ε ≤ n−0.9B . Since in addition 4 < α ≤ n then by Lemma 35 there exists
an integer M ≤ n1.6B satisfying:
∀S ⊆ [n] DM ({mλS }) ≤ 4logs (n)n−5 ≤ n−4 ,

(31)

for sufficiently large n. Hence, by Lemma 39 a random m ∼ U [M ] separates the k-th eigenvalue
of A + E w.r.t. Bin (α), Bout w.p at least 1/(5αn).

6.5

Technical approximations

Fact 40. Approximating a Gaussian by a convex sum of uniform distributions
Let g = (g1 , . . . , gn ) be standard Gaussian vector g ∈ Rn . Then g is equal to a convex combination
of two distributions Du , Dv as follows: (1 − p)DU + p√· DV , where DU is the n-fold distribution
of
Pm
independent variables z1 , . . . , zn , where each zi , |zi | ≤ n, and zi is the convex sum zi = j=1 tj wj
of m ≥ 2n2 i.i.d. uniformly distributed variables, with wj ∼ U [Ij ], with |Ij | = 1/m, and p ≤ 2n2 /m.
√
√
√
Proof. Partition the interval [− n/2, n/2] into m · n equal intervals Ij , each of size 1/m,
and let pj denote the point of minimal absolute value in the j-th interval. Set tj = Φ(pj ).
2
Consider the real Gaussian PDF Φ(x) = √12π e−x /2 . For any pair of neighboring intervals
pk , pk+1 , |pk+1 | > |pk |, we have
!
2
√
√
e−(pj+1 /m) /2
|Φ(pj ) − Φ(pj+1 )| ≤ Φ(pj ) · 1 −
≤ 1 − e−pj /m ≤ 1 − e− n/m ≤ n/m. (32)
2 /2
−(p
)
e j
Then by the above:
∀j ∈ [m]

max {|tj − tj−1 |, |tj − tj+1 |} ≤

√

n/m,

Let zi : R 7→ [0, 1] denote the following function, which is a sum of uniform distributions on
the intervals Ij :
√
m n
X
zi =
tj U [Ij ].
j=1

Consider the l1 -distance between gi and zi :
Z

√
m n

Z
|gi (x) − zi (x)| dx =

x∈R

gi (x) −
x∈R

25

X
j=1

Φ(pj )U [Ij ](x) dx

By Equation 32, and by bounding the Gaussian tail at

√

n we establish the upper bound:

√

m n
1 X
≤
max {|tj − tj+1 |, |tj − tj−1 |} + 2 · e−n ≤ 2n/m.
m j=1

On the other hand, by definition of the points tj , we have
gi (x) − zi (x) > 0, ∀x ∈ R.

(33)

Hence each gi may be written as a convex combination
gi = (1 − p) · ẑi + p · yi ,

ẑi (x) = R

zi (x)
z (x)dx
R i

p ≤ 2n/m.

Since the variables gi are independent then the n-fold distribution of g1 , . . . , gn , can be written
as a convex combination of the of the n-th fold distribution of such i.i.d. variables ẑ1 , . . . , ẑn ,
and another distribution, Dv occurring, by the union bound, w.p. at most n · p ≤ 2n2 /m.

7

Parallel Algorithm for ASD

The algorithm Filter(A, m, δ) described in Section 5 is given an integer m that separates the
i-th eigenvalue, and returns an approximation for the i-th eigenvector. In this section, we use
this algorithm in a black-box fashion and design a Las-Vegas algorithm for computing the full
ASD of a matrix. Essentially, it amounts to running sufficiently many copies of the filtering
algorithm in parallel so that all eigenvectors are (coupon) collected with high probability.

Algorithm 41. Input: n × n Hermitian matrix A, parameter δ ≤ 1/n.
1. Compute parameters:
√
B = min{δ, B ∗ (δ/(3 n))}, δ 0 = (min{δ, B})13 /4
p
α = ln(1/δ 0 ), M = (max{B −12 , n−50 })1.6 , T = 60nαlog(n)
Bout = [1 − 1/(2n), −1 + 1/(2n)] and Bin = [−1/(nα), 1/(nα)],
p
2. Perturb: A1 = A + (δ 0 )2 + δ 2 /(9n) · E, where E is GUE.
3. Run T parallel processes of the following procedure
(a) Sample m ∼ U [M ]
(b) Run Filter (A1 , m, δ 0 ).
4. For vector w = wk sampled at process i ∈ [T ], compute z = A · w, i0 = arg maxi∈[n] |wi |
and λ̃k = zi0 /wi0 .
5. Sort the values λ̃i : assume λ̃1 ≤ . . . ≤ λ̃T . Initialize: γ = λ̃1 , D = Φ. Iterate over all
i = 1, . . . , T . At each step i: if |γ − λ̃i | ≥ B/4 then add D → D ∪ {wi }, and set γ = λ̃i .

We now state our main theorem:

26

Theorem 42. For any Hermitian matrix 0  A  0.9I, and δ ≤ n−10 there exists an RNC(2)
algorithm computing ASD(A, δ), in boolean complexity Õ(nω+1 ). The algorithm is log-stable.
Proof.
Correctness:
Let E1 , E2 be independent GUE matrices, and set
√
A0 = A + (δ/(3 n)) · E1 .
By Lemma 18 and definition 19 the parameter B ≤ B ∗ above satisfies:
P (A0 is B separated) ≥ P (A0 is B ∗ separated) ≥ 0.99.
Assume that this is the case, and account for the error at the end. Consider now a small-scale
perturbation of A0 :
A1 = A0 + δ 0 · E 2 .

(34)

Since E1 , E2 are independent the matrix A1 is a perturbation of A as follows:
q
√
A1 = A + (δ 0 )2 + (δ/(3 n))2 · E,

(35)

for GUE matrix E. Hence the matrix A1 used by the algorithm starting at step 3 is a δ 0 perturbation of a B-separated matrix A0 .
We now invoke Theorem 32 w.r.t. A0 . A0 is separated with parameter B, and perturbed by
GU E with standard deviation ζ = δ 0 ≤ B 13 . We also have δ ≤ n−10 which means in particular
that the parameter
ζ used in the statement of Theorem 32 satisfies ζ ≤ n−50 . Finally we have
p
2 < α = O( ln(n)) = o(n). Therefore:
∀k ∈ [n] PE2 ,m∼U [M ] (m separates k in A1 = A0 + δ 0 E2 w.r.t. Bin (α), Bout ) ≥ 1/(5αn)
Let {vi }ni=1 denote an orthonormal basis for A1 , and let λi denote the eigenvalues of A1 .
Conditioned on sampling A1 , m such that m separates k in A1 , we invoke Theorem
p 26 for
w = Filter(A, m, δ 0 ). By our choice of parameters we have that δ 0 ≤ n−10 , and α = ln(1/δ 0 ).
Hence, the conditions of the theorem are met and so the output vector w satisfies:
Pw0 (kw − vk k ≤ δ 0 ) ≥ 1 − 3n−3 .
Thus, by the union bound we have
∀k ∈ [n] PE2 ,m,w0 (kw − vk k ≤ δ 0 ) ≥

1
1
− 3n−3 ≥
5nα
6nα

Therefore, by the coupon collector’s bound the probability that T = 60nαlog(n) parallel copies
sample an approximation for each k ∈ [n] is at least
1 − 6/60 − 0.01 = 0.89
Assume that this is the case. Let {wi }Ti=1 denote the output set of vectors from all T parallel
copies. We now want to show that each output vector is a valid approximate eigen-vector.
Each output vector satisfies by the stopping criterion:
∃c

√
λ̃i · wi − A1 · wi ≤ 3δ 0 n ≤ B/200,

27

where the last inequality is for sufficiently large n, and λ̃i is the algorithm’s approximation of
the i-th eigenvalue. By Fact 15 since A0 is B-separated then
P(A1 is B/2 separated) ≥ 1 − 2−n .
(36)
√
Assuming this is the case, δ 0 n is negligibly smaller than the separation of A1 - so by Fact 44
the value λ̃i above must satisfy:
∃k ∈ [n] λ̃i − λk ≤ B/10.
and there exists some eigenvector vk of λk such that
kwi − vk k ≤ B/30.
This implies by the triangle inequality that for any j 6= i, l 6= k for which λ̃j − λl ≤ B/10 we
have
λ̃i − λ̃j > B/2 − 2B/10 > B/4.
hence λ̃j and λ̃j are classified to different eigenvalue bins at step 5. We conclude that the set of
values {λ̃i }Ti=1 satisfies that
• For every k ∈ [n] there exists j ∈ T such that: λk − λ̃j ≤ B/10 and there exists a sampled
vector wj , and a unit eigenvector vk of λk such that kwj − vk k ≤ B/30.
• For every j ∈ T there exists a unique k ∈ [n] such that λk − λ̃j ≤ B/10
Therefore, at the end of step 5 we have |D| = n and
∀k ∈ [n] ∃w ∈ D kw − vk k ≤ B/10.
Hence, with probability at least 0.89 − 2−n ≥ 0.85 the algorithm returns a database D =
{w1 , . . . , wn }, such that there exists an orthonormal basis {v1 , . . . , vn } of A1 for which
∀k ∈ [n] kvk − wk k ≤ B/10 ≤ δ/10.
In that case D is an ASD(A1 , δ/10). On the other and, by Equation 35 we have that w.p. 1 − 2−n
A1 satisfies:
kA − A1 k2 ≤ (δ/3)2 .
p
Therefore w.p. at least 0.84 the database D is also ASD(A, (δ/10)2 + (δ/3)2 ) which is in particular an ASD(A, δ).
Run time:
By Theorem 26 the parallel time of each copy of Filter(A, m, δ 0 ) is at most O(log2 (n)), and
arithmetic complexity Õ(nω ). Hence, this is the depth complexity of all T parallel copies of
the circuit for Filter(A, m). The sorting step 5 can be implemented in depth O(log(n)) by Fact
11, using O(n) processors. Hence the complete algorithm runs in parallel time O(log2 (n)) in
arithmetic complexity Õ(nω+1 ).
Finally, we show that the entire algorithm can be implemented in log(n)-precision:
Proposition 43. The algorithm 41 is log-stable. In particular it runs in boolean complexity Õ(nω+1 ).
Proof. Each parallel process is log-stable by Claim 27. Hence, to show that the entire algorithm is log-stable, it is sufficient to show that for a Gaussian perturbation truncated to
t = O(log(n)) bits of precision, the statement of Theorem 32 still applies - namely a uniformly

28

random m ∼ U [M ] separates each eigenvalue, with high probability, for the parameter M set
in the algorithm. To do that, write
Et = E + D,
where E is the non-truncated GUE matrix, and D is some matrix sampled from a distribution
such that |Di,j | ≤ 2−t for all i, j. Following Equation 34 define
(t)

A1 ≡ A0 + E + D = A1 + D,
(t)

which is then used by the parallel process of Filter(A1 , m, δ 0 ). By Equation 36 w.p. at least
1 − 2−n the matrix A1 is separated with parameter at least B/2. Hence, for sufficiently large
t = O(log(n)), such that 2t > B and using again Fact 15 we have for each i ∈ [n]:
√
(t)
λi (A1 ) − λi (A1 ) ≤ 2 n · 2−t .
Hence by Lemma 21 if
({mλi1 (A1 )}, . . . , {mλi2 (A1 )})m∈[M ] ,
is a 2-dimensional sequence with discrepancy DM , then the discrepancy of the noisy sequence:
√
(t)
(t)
({mλi1 (A1 )}, . . . , {mλi2 (A1 )})m∈[M ] is at most DM + 2 · M · n · 2−t
√ . This implies that the
pairwise discrepancy error increases by an additive error at most 2M n2−t at Equation 31 in
the proof of theorem 32. Hence, there exists a choice of t = O(log(n)) so that the discrepancy
of the above sequence is at most n−4 , thereby satisfying Equation 31. Hence, the statement of
Theorem 32 holds for E = Et , for t = O(log(n)).
This implies that the algorithm can be computed in Õ(nω+1 ) arithmetic operation, where
each operation is a binary operation between two registers of size O(log(n)). Hence, the
boolean complexity of the algorithm is Õ(nω+1 ) · log(n) = Õ(nω+1 ).

7.1

Supporting Claims

Fact 44. Let A be a δ-separated matrix, w a unit vector satisfying
kAw − cwk ≤ B, B ≤ δ 2 /100,
for some |c| ≤ 1. Then there exists λ ∈ L(A) such that
|c − λ| ≤ δ/10,
and there exists some unit eigenvector v of λ such that
kv − wk ≤ δ/30.
Proof. Without loss of generality, we assume w = α1 v1 + α2 v2 where v1 , v2 are two distinct
elements of some orthonormal basis of A, corresponding to eigenvalues λ1 , λ2 . Then
kAw − cwk2 = kv1 α1 (c − λ1 ) + v2 α2 (c − λ2 )k2 ≤ B 2
By the fact that v1 , v2 are orthogonal:
|α1 |2 · |c − λ1 |2 + |α2 |2 · |c − λ2 |2 ≤ B 2 .
Since kwk = 1 then |α1 |2 + |α2 |2 = 1, so for at least one i ∈ {1, 2} we have
|c − λi |2 ≤ 2B 2 .

29

(37)

Suppose w.l.o.g. that i = 1. Since A is δ separated then |λ1 − λ2 | ≥ δ. Then together with the
triangle inequality we have
√
|λ2 − c| ≥ |λ2 − λ1 | − |c − λ1 | ≥ δ − 2B ≥ 0.9δ.
Hence by Equation 37 above:
|α2 |2 ≤ B 2 /(0.9δ)2 ≤ δ 2 /1000
This implies that for some unit eigenvector v of λ1 we have
kw − vk ≤ δ/30,
Using this inequality back in the assumption, and since kAk ≤ 1, |c| ≤ 1 then
kcw − Awk = kc(v + E) − A · (v + E)k ≤ B,

kEk ≤ δ/30

and since Av = λ1 v then
kcv − λ1 vk ≤ B + 2δ/30,
and since kvk = 1 we get:
|c − λ1 | ≤ B + 2δ/30 ≤ δ/10.

8

Numerical Experiments

In the work above we have outlined a new theoretical approach for solving the eigen-problem
of Hermitian matrices. Given the immense practical importance of this problem, we now provide in addition some numerical evidence about the performance of the algorithm. We run
the main algorithm for ASD of a random matrix A for n = 20, where the number of bits of
precision used is standard Matlab single precision of 32 bits, and require output precision of
δ = 10−4 . We then run 50 · n · logn ∼ 2500 copies of the filter procedure F ilter(A, m, δ). Then,
we iterate over a sequence of values for the integer M , i.e. the number such that such that by
Theorem 32 a uniformly random integer m ∼ U [M ] separates any eigenvalue of a perturbed
version of A w.h.p.
One can see that while the analytic behavior requires M to be at least n75 , here we see that
even M = n5 already suffices for the algorithm to return an approximate eigenvector for each
eigenvalue w.p. very close to 1. We conjecture, hence, that implementing the algorithm will
result in far better practical run time than the analytical bounds we provide here.

Acknowledgements
The authors thank Naomi Kirshner and Robin Kothari for very helpful comments, and Yosi
Atia for pointing to us an error in an earlier version of this paper. We also thank anonymous
reviewers for their helpful comments and suggestions. This research project was supported in
part by the Israeli Centers of Research Excellence (I-CORE) program (Center No. 4/11), by the
Israeli Science Foundation (ISF) research grant 1446/09, by an EU FP7 ERC grant (no.280157),
and by the EU FP7-ICT project QALGO (FET-Proactive Scheme). LE is thankful to the Templeton Foundation for their support of this work.

30

Figure 1: Algorithm Behavior

References
[Akl90]

Selim G. Akl. Parallel Sorting Algorithms. Academic Press, Inc., Orlando, FL, USA,
1990.

[BP92]

D. Bini and V. Pan. Practical improvement of the divide-and-conquer eigenvalue
algorithms. Computing, 48(1):109–123, 1992.

[DDH07] James Demmel, Ioana Dumitriu, and Olga Holtz. Fast linear algebra is stable. Numerische Mathematik, 108(1):59–91, 2007.
[GVL96] Gene H. Golub and Charles F. Van Loan. Matrix Computations (3rd Ed.). Johns Hopkins University Press, Baltimore, MD, USA, 1996.
[HHL09] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd. Quantum algorithm for
linear systems of equations. Phys. Rev. Lett., 103:150502, Oct 2009.
[Koz92]

Dexter C. Kozen. The Design and Analysis of Algorithms. Springer-Verlag New York,
Inc., New York, NY, USA, 1992.

[LMM99] Mauro Leoncini, Giovanni Manzini, and Luciano Margara. Parallel complexity of
numerically accurate linear system solvers. SIAM Journal on Computing, 28(6):2030–
2058, 1999.
[Nie92]

H. Niederreiter. Random Number Generation and Quasi-Monte Carlo Methods. Society
for Industrial and Applied Mathematics, 1992.

[NTV16] Hoi Nguyen, Terence Tao, and Van Vu. Random matrices: tail bounds for gaps
between eigenvalues. Probability Theory and Related Fields, pages 1–40, 2016.
[Rei05]

John H. Reif. Efficient parallel factorization and solution of structured and unstructured linear systems. Journal of Computer and System Sciences, 71(1):86 – 143, 2005.

31

[SS90]

Gilbert W. Stewart and Jiguang Sun. Matrix perturbation theory. Computer science
and scientific computing. Academic Press, Boston, 1990.

[TB97]

Lloyd N. Trefethen and David Bau. Numerical linear algebra. Society for Industrial
and Applied Mathematics, Philadelphia, 1997.

[TS13]

Amnon Ta-Shma. Inverting well conditioned matrices in quantum logspace. In
Proceedings of the Forty-fifth Annual ACM Symposium on Theory of Computing, STOC
’13, pages 881–890, New York, NY, USA, 2013. ACM.

[Wil12]

Virginia Vassilevska Williams. Multiplying matrices faster than coppersmithwinograd. In Proceedings of the Forty-fourth Annual ACM Symposium on Theory of
Computing, STOC ’12, pages 887–898, New York, NY, USA, 2012. ACM.

32

