# Diagonalization of Dense Ill Conditioned Matrices

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

-
What routine are you using to diagonalize your real symmetric matrices? And in what sense is the eigenvalue decomposition inaccurate? – Jack Poulson Jul 8 '12 at 23:27
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. – J. M. Jul 12 '12 at 5:34
@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. – Arnold Neumaier Jul 12 '12 at 10:48
@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... – J. M. Jul 12 '12 at 10:51
@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. – Arnold Neumaier Jul 12 '12 at 10:54

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.

-
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. – Jack Poulson Jul 9 '12 at 14:19
@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. – Arnold Neumaier Jul 9 '12 at 15:15
Reference? Jacobi's method is known to be highly accurate (albeit slow). – Jack Poulson Jul 9 '12 at 15:21
@JackPoulson: Try and see. Demmel & Kahan, Accurate Singular Values of Bidiagonal Matrices, 202.38.126.65/oldmirrors/ftp.netlib.org/lapack/lawnspdf/… – Arnold Neumaier Jul 9 '12 at 15:31
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. – Jack Poulson Jul 9 '12 at 16:05