10
$\begingroup$

I am trying to diagonalize some dense, ill-conditioned matrices. In machine precision, results are inaccurate (returning negative eigenvalues, eigenvectors do not have the expected symmetries). I switched over to Mathematica's Eigensystem[] function to take advantage of arbitrary precision, but computations are extremely slow. I am open to any number of solutions. Are there packages/algorithms that are well suited to ill-conditioned problems? I am not an expert on preconditioning, so I am not sure how much this could help. Otherwise, all I can think of are parallelized arbitrary precision eigenvalue solvers, but I am not familiar with anything beyond Mathematica, MATLAB and C++.

To give some background on the problem, the matrices are large, but not huge (4096x4096 to 32768x32768 at the most). They are real, symmetric, and the eigenvalues are bounded between 0 and 1 (exclusive), with many eigenvalues being very close to 0 and none close to 1. The matrix is essentially a convolution operator. I do not need to diagonalize all of my matrices, but the larger I can go, the better. I have access to computing clusters with many processors and distributed computing capabilities.

Thank you

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ What routine are you using to diagonalize your real symmetric matrices? And in what sense is the eigenvalue decomposition inaccurate? $\endgroup$
    – Jack Poulson
    Jul 8, 2012 at 23:27
  • $\begingroup$ Here's an idea related to Arnold's answer: perform a Cholesky decomposition of your SPD matrix, and then find the singular values of the Cholesky triangle you have just obtained, possibly using a dqd-type algorithm to preserve accuracy. $\endgroup$
    – J. M.
    Jul 12, 2012 at 5:34
  • 1
    $\begingroup$ @J.M.: Forming the Cholesky decompositon of a numerically singular positive definite matrix is numerically unstable with the usual method, as one likely encounters negative pivots. (E.g., Matlab's chol(A) typically fails.) One would have to set them to zero and annihilate the corresponding rows of the factors. Doing this gives a way to reliably get the numerical null space. $\endgroup$
    – Arnold Neumaier
    Jul 12, 2012 at 10:48
  • $\begingroup$ @Arnold, if memory serves, there are adaptations of Cholesky that use symmetric pivoting for those cases where the matrix is positive semi-definite (or nearly so). Maybe those could be used... $\endgroup$
    – J. M.
    Jul 12, 2012 at 10:51
  • $\begingroup$ @J.M.: One doesn't need pivoting to resolve the semidefinite case; the recipe I gave is enough. I just wanted to point out that one cannot use the standard canned programs but has to modify them oneself. $\endgroup$
    – Arnold Neumaier
    Jul 12, 2012 at 10:54

2 Answers 2

Reset to default
7
$\begingroup$

Calculate the SVD in place of the spectral decomposition. The results are the same in exact arithmetic, as your matrix is symmetric positve definite, but in finite precision arithmetic, you'll get the small eigenvalues with much more accuracy.

Edit: See Demmel & Kahan, Accurate Singular Values of Bidiagonal Matrices, SIAM J. Sci. Stat. Comput. 11 (1990), 873-912.
ftp://netlib2.cs.utk.edu/lapack/lawnspdf/lawn03.pdf

Edit2; Note that no method will be able to resolve eigenvalues smaller than about the norm times the machine accuracy used, as changing a single entry by one ulp may already change a small eigenvalue by this much. Thus getting zero eigenvalues in place of very tiny ones is appropriate, and no method (except working with higher precision) will disentangle the corresponding eigenvectors, but just return a basis for the common numerical null space.

Improve this answer
$\endgroup$
7
  • $\begingroup$ I'm not sure I believe this, as most SVD implementations start with a unitary reduction to real bidiagonal form and then essentially compute the Hermitian EVD of a related matrix, such as $[0, B^T; B, 0]$, which can easily be permuted to real symmetric tridiagonal form. The relative accuracy is highly dependent on which method is used to solve the condensed EVP/SVD, and I don't see where the SVD has an advantage...I'm sure this is discussed in one or more of Demmel's papers. $\endgroup$
    – Jack Poulson
    Jul 9, 2012 at 14:19
  • 2
    $\begingroup$ @JackPoulson: The point is that the bidiagonal form determines small singular values much better. The associated symmetric tridiagonal form has zeros on the diagonal, which are preserved by the bidiagonal reduction to diagonal, but not by QR applied to the tridiagonal. $\endgroup$
    – Arnold Neumaier
    Jul 9, 2012 at 15:15
  • 1
    $\begingroup$ Reference? Jacobi's method is known to be highly accurate (albeit slow). $\endgroup$
    – Jack Poulson
    Jul 9, 2012 at 15:21
  • $\begingroup$ @JackPoulson: Try and see. Demmel & Kahan, Accurate Singular Values of Bidiagonal Matrices, 202.38.126.65/oldmirrors/ftp.netlib.org/lapack/lawnspdf/… $\endgroup$
    – Arnold Neumaier
    Jul 9, 2012 at 15:31
  • $\begingroup$ I see what you're getting at: the QR algorithm for the bidiagonal SVD can be made relatively accurate because of the zero diagonal of the permuted $[0,B^T; B,0]$ matrix, whereas this "zero-shift" technique does not work for arbitrary tridiagonal matrices, so QR-algorithm based tridiagonal EVP's will be less accurate for small eigenvalues. The catch is that this reasoning assumes an eigensolver based on using the QR algorithm on the condensed form; Jacobi is a notable exception, but perhaps it is slower than a QR-algorithm based SVD. MRRR can sometimes also achieve high relative accuracy. $\endgroup$
    – Jack Poulson
    Jul 9, 2012 at 16:05
1
$\begingroup$

Thank you for this suggestion. I tried Mathematica's SVD command, but I get no noticeable improvement (still missing appropriate symmetries, 'eigenvalues' are incorrectly zero where they were incorrectly coming out negative before). Maybe I would need to implement one of the algorithms you describe above instead of a built-in function? I would probably want to avoid going to the trouble of using a specific method like this unless I was sure ahead of time that it would offer a significant improvement.

@JackPoulson, I skimmed the paper on Jacobi's method you referenced, and it looks promising. Can you or anyone recommend a good way to implement Jacobi's method for finding eigensystems? I am guessing that if I coded it up myself (in MATLAB), it would be extremely slow.

Improve this answer
$\endgroup$
6
  • $\begingroup$ I haven't tested it, but there is a MATLAB implementation here: groups.google.com/forum/?fromgroups#!msg/sci.math.num-analysis/… $\endgroup$
    – Jack Poulson
    Jul 9, 2012 at 22:20
  • $\begingroup$ Note that no method will be able to resolve eigenvalues smaller than about the norm times the machine accuracy used, as changing a single entry by one ulp may already change a small eigenvalue by this much. Thus getting zero eigenvalues in place of very tiny ones is appropriate, and no method (except working with higher precision) will disentangle the corresponding eigenvectors, but just return a basis for the common numerical null space. What do you need the eigenvalues for? $\endgroup$
    – Arnold Neumaier
    Jul 11, 2012 at 7:17
  • $\begingroup$ @ArnoldNeumaier: I ran some tests in MATLAB with eigenvalues in the range of [0,1], with one eigenvalue manually set to values such as 6.3e-16, and Octave's SVD routine (based on dgesvd, which uses reduction to bidiagonal and then QR) does pick up these values much more accurately than Octave's eig. The linked Jacobi code seems to be too slow to use, even on modest-sized matrices. $\endgroup$
    – Jack Poulson
    Jul 11, 2012 at 13:36
  • $\begingroup$ @JackPoulson: Yes. But Leigh seems to complain about multiple very tiny eigenvalues, and their eigenvectors will rarely be those designed but will mix freely, no matter which method is used. And positive very tiny positive values (smaller than 1e-16) will of course be found zero. $\endgroup$
    – Arnold Neumaier
    Jul 11, 2012 at 14:07
  • $\begingroup$ @ArnoldNeumaier is right that I am finding multiple very small eigenvalues, which I am guessing exacerbates the problem. I did not realize (though in retrospect it is obvious) that eigenvalues less than 1e-16 will be zero in floating point. I guess although the number can be stored, rounding error occurs when adding it to a larger number. The eigenvectors tell me if a certain problem is solvable. The eigenvector allow decomposition of the problem into solvable and non-solvable parts. If I am fundamentally limited by precision, then can you recommend any packages for faster solution? $\endgroup$
    – Leigh
    Jul 11, 2012 at 17:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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