| Summary: | <p>This report may be considered as a non-trivial extension of an unpublished report by William Kahan (<em>Accurate Eigenvalues of a symmetric tri-diagonal matrix, Technical Report CS 41, Computer Science Department, Stanford University, 1966</em>). His interplay between matrix theory and computer arithmetic led to the development of algorithms for computing accurate eigenvalues and singular values. His report is generally considered as the precursor for the development of IEEE standard 754 for binary arithmetic. This standard has been universally adopted by virtually all PC, workstation and midrange hardware manufactures and tens of billions of such machines have been produced. Now we use the features in this standard to improve the original algorithm.</p> <p>In this paper, we describe an algorithm in floating-point arithmetic to compute the exact inertia of a real symmetric (shifted) tridiagonal matrix. The inertia, denoted by the integer triplet (<em>π, ν, ζ</em>), is defined as the number of positive, negative and zero eigenvalues of a real symmetric (or complex Hermitian) matrix and the adjective exact refers to the eigenvalues computed in exact arithmetic. This requires the floating-point computation of the diagonal matrix <em>D</em> of the <em>LDL</em><sup>t</sup> factorization of the shifted tridiagonal matrix <em>T</em> − <em>τI</em> with +∞ and −∞ rounding modes defined in IEEE 754 standard. We are not aware of any other algorithm which gives the exact answer to a numerical problem when implemented in floating-point arithmetic in standard working precisions. The guaranteed intervals for eigenvalues are obtained by bisection or multisection with this exact inertia information. Similarly, using the Golub–Kahan form, guaranteed intervals for singular values of bidiagonal matrices can be computed. The diameter of the eigenvalue (singular value) intervals depends on the number of shifts with inconsistent inertia in two rounding modes. Our algorithm not only guarantees the accuracy of the solutions but is also consistent across different IEEE 754 standard compliant architectures. The unprecedented accuracy provided by our algorithms could be also used to debug and validate standard floating-point algorithms for computation of eigenvalues (singular values). Accurate eigenvalues (singular values) are also required by certain algorithms to compute accurate eigenvectors (singular vectors).</p> <p>We demonstrate the accuracy of our algorithms by using standard matrix examples. For the Wilkinson <em>W</em><sup>+</sup><sub style="position:relative;left:-.7em;">21</sub>, the eigenvalues (in IEEE double precision) are very accurate with an (open) interval diameter of 6 ulps (units of the last place held of the mantissa) for one of the eigenvalues and lesser (down to 2 ulps) for others. These results are consistent across many architectures including Intel, AMD, SGI and DEC Alpha. However, by enabling IEEE double extended precision arithmetic in Intel/AMD 32-bit architectures at no extra computational cost, the (open) interval diameters were reduced to one ulp, which is the best possible solution for this problem. We have also computed the eigenvalues of a tridiagonal matrix which manifests in Gauss–Laguerre quadrature and the results are extremely good in double extended precision but less so in double precision. To demonstrate the accuracy of computed singular values, we have also computed the eigenvalues of the <em>Kac<sub>30</sub></em> matrix, which is the Golub–Kahan form of a bidiagonal matrix. The tridiagonal matrix has known integer eigenvalues. The bidiagonal Cholesky factor of the Gauss–Laguerre tridiagonal is also included in the singular value study.</p>
|
|---|