see title.

An algorithm is 'good' if it is able to distinguish between zero Eigenvalues and nonzero Eigenvalues.

  • $\begingroup$ I'm would also be very interested to hear an answer to this question, if "best" is interpreted as "best in terms of theoretical guarantees" and not "best in practice." $\endgroup$ – alex May 11 '10 at 20:20
  • 3
    $\begingroup$ To make things my inquiry more specific, suppose the matrix is $n \times n$ and entries are rational numbers, taking $K$ bits to specify the numerator and denominator. How many operations do you need to give me: (i) the first $m$ bits of the smallest eigenvalue (ii) whether the smallest eigenvalue equals zero or not. $\endgroup$ – alex May 11 '10 at 20:22
  • $\begingroup$ for (ii), you just have to compute the determinant and check if it is zero. The traditional bound is $2/3n^3$ operations, but much depends on what is an "operation" if you are dealing with rational numbers. You may want to compute everything exactly, and in this case you need to operate with very large integers, or rely on an approximation. In the former case, the time is probably going to be exponential, since the worst-case bounds on the coefficient growth in Gaussian elimination are exponential. In the latter, it's just $2/3n^3 \log \epsilon$, with $\epsilon$ the desired precision. $\endgroup$ – Federico Poloni Oct 1 '10 at 12:24
  • 1
    $\begingroup$ You may also rely on an algorithm of this kind: reduce to integer coefficients; then compute the determinant $\mod p$ for several values of $p$, and check if it is zero. If you are fine with a probabilistic result, check just a few primes. If you want to be 100% sure of your computation, I am afraid you will need exponentially many primes (or exponentially large primes, which is the same). $\endgroup$ – Federico Poloni Oct 1 '10 at 12:28

One method is to reduce the computation to that of computing matrix multiplication of $n \times n$ matrices. In particular, the determinant of a symbolic matrix can be computed in $O(n^{\omega})$ arithmetic operations, where $\omega < 2.376$ is the matrix multiplication exponent, and from a symbolic determinant of course one can recover all eigenvalues. However, since the operations here will be over polynomials of degree $n$ with coefficients in $m$ bits, this method would take about $O(n^{1+\omega} m)$ time to get $m$ bits of the eigenvalues.

More complex methods can get you the eigenvalues in $O(n^3 + n^2 \log^2 n \log b)$ time, where the eigenvalues are approximated to within $2^{-b}$. For some structured matrices you can get about $O(n^{\omega})$. See

Victor Y. Pan, Zhao Q. Chen: The Complexity of the Matrix Eigenproblem. STOC 1999: 507-516

(Actually it appears this paper never appeared in a journal form, so study it very carefully if you are serious about this problem.)

I don't see a simple way to exploit the fact that (a) it is symmetric and (b) you just want to find the smallest nonzero eigenvalue. It seems doubtful to me that you could do this much faster than $O(n^{\omega})$ (without finding all nonzero eigenvalues faster than this), but this is just based on intuition, not fact.

  • 2
    $\begingroup$ In my opinion, this is not an answer to the question. Because complexity has nothing to do here. What is required is a method with stability properties under roundoff errors. Especially because the smallest eigenvalue is likely to be close to other ones. $\endgroup$ – Denis Serre Nov 14 '10 at 8:29
  • 1
    $\begingroup$ An examination of the above reference yields that the stated $O(n^3)$ complexity is only arithmetic complexity, and the authors leave as an open problem, what is the bit-complexity of their scheme, suggesting the use of the well-known Chinese remainder theorem. $\endgroup$ – Lior Eldar Jul 1 '13 at 8:39
  • $\begingroup$ Also, the complexity of computing the eigenvalues themselves is only (arithmetic) $O(n^{\omega})$. The $O(n^3)$ complexity is the result of computing the {\it eigenvectors} - all $n$ of them at once, instead of $n$ iterations of Gaussian elimination. $\endgroup$ – Lior Eldar Jul 1 '13 at 8:44

Tao and Vu (1) have shown that the distribution of the smallest singular value of a random matrix is "universal", i.e. independent of the particular random variable populating the matrix. They are interested in analyzing distributions, not particular matrices, but it appears that it might be possible to use their machinery for addressing this problem.

First, several caveats.
1. Their analysis assumes there are no zero eigenvalues. With random matrices, this isn't much of a restriction, but it might be for your application. If there were at most a small number of zero eigenvalues, one could pick them off somewhat efficiently and then reduce to the invertible case. 2. The results hold with high probability, not certainty. 3. The results recover estimates of eigenvalues, not the exact values, although they can be made arbitrarily accurate. 4. Tao and Vu's analysis is not aimed at this particular application, so there very well may be some issue translating them to this domain.

Now that the pussy-footing is out of the way, the basic idea is as follows: We wish to find the smallest eigenvalue $\lambda_n$ of a matrix $A$. Suppose that $A$ is invertible. Then the largest eigenvalue of $A^{-1}$ is $1/\lambda_n$. Since finding the largest eigenvalue is a much easier problem (by, e.g. the power method), we might be in a better position.

But computing $A^{-1}$ is (generically) as difficult as computing all the eigenvalues! Tao and Vu dodge this problem by taking random subsets of the columns of $A$ and considering the orthogonal complement; this is inspired by "property testing" arguments in complexity theory. Then with high probability one can estimate the largest eigenvalue of $A^{-1}$ from these restrictions, and we are done.

(1) Terence Tao and Van Vu. "Random matrices: the distribution of the smallest singular values". March, 2009. http://arxiv.org/abs/0903.0614v1


The QR algorithm gives quite rapidly a good approximation of the eigenvalues of a real symmetric (or complexHermitian). But overall, it gives first the smallest eigenvalues. The reason is that the ratio $\lambda_{n-1}/\lambda_n$ is the convergence rate.

  • $\begingroup$ "it gives first the smallest eigenvalues" - not necessarily true if the Wilkinson shift is used, which is often the case in modern QR implementations. $\endgroup$ – J. M. is not a mathematician Oct 1 '10 at 12:05
  • $\begingroup$ @J.M. The Wilkinson shift uses the behaviour I mentionned above. You get a first approximation $\mu$ of the smallest e.v. $\lambda_n$, then shift by this $\mu$. This accelerate the convergence to $\lambda_n$ because $(\lambda_{n-1}-\mu)/(\lambda_n-\mu)$ is very large. Once you like the approximate e.v., you store it and delete the last row and column of the matrix before continuing. Incredibly fast ! $\endgroup$ – Denis Serre Oct 1 '10 at 13:32
  • $\begingroup$ @Denis: What I meant was that the first eigenvalue you get using tridiagonal QR/QL+Wilkinson shift isn't necessarily the tiniest (in magnitude) eigenvalue; an example I have offhand is the positive definite tridiagonal Toeplitz matrix with 2 on the diagonal and 1 on the off-diagonals. QR gives the eigenvalue 3 first; QL gives the eigenvalue 1 first; neither return the tiniest eigenvalue (0.267949...) first. $\endgroup$ – J. M. is not a mathematician Oct 1 '10 at 14:36
  • $\begingroup$ ...though I suspect whatever is returned first by either version is machine-dependent. $\endgroup$ – J. M. is not a mathematician Oct 1 '10 at 14:37

First, a disclaimer: I know absolutely nothing about numerical algorithms for finding the eigenvalue of a matrix, symmetric or not. So my feelings will not be hurt if my answer gets downvoted into oblivion.

It seems to me that an obvious but perhaps overly naive approach is the following:

Let $A$ be the symmetric matrix in question.

a) Use some standard minimization algorithm (maybe the conjugate gradient method?) to minimize

$$\frac{x\cdot Ax}{x\cdot x}$$

over nonzero $x$ (an obvious thing to do is to restrict to $|x| = 1$ but you might save some arithmetic if you don't bother with this normalization).

b) See what eigenvalue you get. If it's the eigenvalue you want, then you're done. If not, save the eigenvector you found and proceed to c)

c) Repeat a), except restrict to the subspace orthogonal to all of the eigenvectors you've found so far. See what eigenvalue you get. If it's the one you want, you're done. Otherwise, save the eigenvector and repeat this step again.

Eventually, you'll have all of the eigenvalues and eigenvectors. Depending on what "smallest" means, you may or may not be able to stop before you have found all of the eigenvectors. Actually, if "smallest" means "eigenvalue with the smallest nonzero absolute value", then just do the steps above with $A^2$ instead of $A$.

For small matrices this seems like a practical approach to me. But, as I said, I don't know anything about this stuff.


This is the method they use with LaPACK, which is usually the fastest for general problems (fastest noncommercial anyway)


and here's a discussion regarding this computation


I honestly am not sure if you can hunt down the smallest eigenvalue without finding all of them.

However, I would not under any circumstances do the symbolic determinant. So far as I know if you're doing a numerical calculation, an introduction of symbols will give you a major slowdown. I don't have the time at the moment to look up the exact computation, but I think that outside of quantum computing symbolic factorization is NP-hard/NP-complete. In the multivariate case (using Groebner Bases) it's doubly exponential.


Probably still your best bet, after of course reducing your original symmetric matrix to tridiagonal form, would be either bisection (with the help of Gerschgorin bounds) or an appropriate modification of the dqd/MRRR algorithm of Parlett, Fernando, and Dhillon.

If your matrix has additional structure apart from symmetry (e.g. it is a Toeplitz or arrowhead matrix), there of course may be even more slick apporaches. I suggest looking at the references given in LAPACK and other numerical linear algebra books.


If smallest means closest to zero, than Rayleigh quotient iterations gives cubic convergence, but still requires matrix inverse, or actually to solve the system of linear equations.


see also http://en.wikipedia.org/wiki/Inverse_iteration


A quick search led me to this paper, which deals specifically with sparse symmetric matrices, although some of its references might be useful.

Jang, Ho-Jong, and Lee, Sung-Ho, "NUMERICAL STABILITY OF UPDATE METHOD FOR SYMMETRIC EIGENVALUE PROBLEM," J. Appl. Math. & Computing Vol. 22(2006), No. 1 - 2, pp. 467 - 474.

A PDF copy is available here: http://www.mathnet.or.kr/mathnet/kms_tex/986075.pdf

I should also mention that "best" is a difficult superlative to qualify without knowing the structure of your matrices. Probably the best algorithm for a sparse symmetric matrix is not the best algorithm for a symmetric Toeplitz matrix.


Here is a simple idea, with no complexity analysis. Compute a basis $v_1$, ..., $v_{n-r}$ for the kernel of $A$; this can be done with exact arithmetic in $n^3$ operations by Gaussian elimination. Compute a basis $w_1$, ..., $w_{r}$ for the othogonal complement: Again, doable in exact arithmetic with $n^3$ operations. $\mathrm{Span}(w_1, ..., w_r)$ is the orthogonal complement to the $0$-eigenspace of $A$ and hence, since $A$ is symmetric, it is the span of the nonzero eigenspaces of $A$. So $A$ maps $\mathrm{Span}(w_1, \ldots, w_r)$ to itself. Let $X$ be the $r \times r$ matrix of this map. This is again a matrix of rational numbers, computable with exact arithmetic in reasonable time, whose eigenvalues are the same as the nonzero eigenvalues of $A$. Note, however, that $X$ is not symmetric.

Invert $X$ and find its largest eigenvalue by one of the standard methods.

What I am gambling here is that the advantages of working in exact arithmetic are greater than the disadvantages of passing to a nonsymmetric matrix.

  • $\begingroup$ Is $X$ symmetric if $w_1,...,w_r$ are chosen to form an orthogonal basis ? $\endgroup$ – KBuck Jun 26 '12 at 15:35
  • $\begingroup$ Yeah, but to find an orthonormal basis you have to take square roots, which either requires passing to floating point or working in number fields. $\endgroup$ – David E Speyer Jun 26 '12 at 15:50
  • $\begingroup$ One can imagine a compromise where one first computes $X$ exactly, then computes the transition matrix to an orthogonal basis $w'_i$ numerically, and thus changes to a symmetric floating point matrix $X'$ which one finds the eigenvalues of. $\endgroup$ – David E Speyer Jun 26 '12 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.