Inaccurate eigenvectors with dstedc

[caliarim]

Dear all,

I think dstedc produces inaccurate eigenvectors for a specific tridiagonal symmetric matrix. In order to show it, I assembled the attached test. It computes eigenpairs of the Jacobi matrix related to Gauss-Hermite quadrature formulas (Golub-Welsch algorithm). Then, it applies the quadrature formula to int_{-infty}^{infty} x^34 * exp(-x^2) dx. The quadrature formula should be exact for 18 or more points. It is numerically exact if the eigenpairs are computed with dsteqr. But it diverges (as the number of points increases) with dstedc. I tested lapack 3.7.1. tsteig.f.zip

[langou]

Hi @caliarim

Jim Demmel followed up on this. Here is what Jim says:

It seems that Matlab (using dstevd) gets O(macheps) values for norm(VV'-eye) and norm(JV-V*D) where [V,D]=eig(J), which is all we promise. The user prefers the result of V from dstev, in particular the entries of V much smaller than macheps, which is something that we can't make any general guarantees about of course, it depends on the special structure of your matrix. Allowing an option to use dstev instead of dstevd sounds like the best option for those expert users who know they have such special structure. We can also refer the user to some literature about what special structures we have studied, but that certainly won't cover all cases.

One more comment on the special structure of this problem. Since the matrix is tridiagonal with 0 diagonal, it is really equivalent to the bidiagonal SVD, which we have extensively studied. This is an "easy" problem in the sense that the condition number (for m=100) is 121, and the minimum relative gap between eigenvalues/singular values is .031, so all the eigen/singular values are determined to high relative accuracy, and so are the eigen/singular vectors, in a norm sense. But the the user wants componentwise high relative accuracy, even for the tiniest components, which we at least did not analyze/promise.

Jim

Cheers, Julien.