Computing accurate eigenvalues using a mixed-precision Jacobi algorithm

TL;DR

This work develops a mixed-precision framework for computing eigenvalues of symmetric matrices by preconditioning A with a low-precision-based matrix Q, applying the preconditioner at high precision, and solving with the Jacobi method at working precision. The authors prove forward-error bounds that are tighter than classical Jacobi bounds and, under small off-diagonal preconditioning, can be independent of the original condition number. They present two concrete preconditioner constructions exploiting low precision, derive detailed error and condition-number bounds, and confirm the theory with comprehensive numerical experiments showing substantial accuracy gains, albeit with higher cost dominated by high-precision matrix multiplications. The results suggest practical pathways to high-accuracy spectral decompositions, including dynamic precision selection and potential extensions to singular values and one-sided problems as hardware support for high-precision arithmetic improves.

Abstract

We provide a rounding error analysis of a mixed-precision preconditioned Jacobi algorithm, which uses low precision to compute the preconditioner, applies it at high precision (amounting to two matrix-matrix multiplications) and solves the eigenproblem using the Jacobi algorithm at working precision. Our analysis yields meaningfully smaller relative forward error bounds for the computed eigenvalues compared with those of the Jacobi algorithm. We further prove that, after preconditioning, if the off-diagonal entries of the preconditioned matrix are sufficiently small relative to its smallest diagonal entry, the relative forward error bound is independent of the condition number of the original matrix. We present two constructions for the preconditioner that exploit low precision, along with their error analyses. Our numerical experiments confirm our theoretical results and compare the relative forward error of the proposed algorithm with the Jacobi algorithm, a preconditioned Jacobi algorithm, and MATLAB's $\texttt{eig}$ function. Timings using Julia suggest that the dominant cost of obtaining this level of accuracy comes from the high precision matrix-matrix multiplies; if support in software or hardware for this were improved, then this would become a negligible cost.

Figures (7)

• Figure 1: Size of ${\operatorname{off}}({\widetilde{A}})/\|{A}\|_F$ against the matrix size $n$ of ${\widetilde{A}}=\widetilde{Q}^{T} A\widetilde{Q}$ with $\widetilde{Q}$ computed by four different algorithms specified in the legend. The subcaptions state the MODE used in . In all cases $\kappa_{2}(A) = 10^{6}$.
• Figure 2: Maximum relative forward error, $\max_{k}{{\varepsilon}_{fwd}^{(k)}}$, against $\kappa_{2}(A)$ for four different algorithms and $A\in\mathbb{R}^{100\times 100}$ generated by with different values of MODE.
• Figure 3: Maximum relative forward error, $\max_{k}{{\varepsilon}_{fwd}^{(k)}}$, for four different algorithms against the size $n$ of matrices $A$ generated by with different values of MODE and $\kappa_{2}(A) = 10^{10}$.
• Figure 4: Relative forward error for the eigenvalues ordered from largest to smallest of the test matrices in section \ref{['sec.evals-special-mtx']}. Subfigure (d) only displays ${{\varepsilon}_{fwd}^{(k)}}$ for every tenth eigenvalue, and the data missing for MP3Jacobi and MATLAB is because the error was smaller than $10^{-16}$.
• Figure 5: Ratio of the scaled condition number of $A$ and ${\widetilde{A}}$, and $\kappa_{2}(A)$. The subcaptions state the singular value distributions for the test matrices $A$ which have fixed size $n=100$.

Theorems & Definitions (22)

• Theorem 2.2
• Lemma 2.3
• Proof 1
• Remark 2.4
• Lemma 2.5
• Proof 2
• Theorem 2.6: vesl93
• Theorem 2.7
• Theorem 2.8: chu95
• Corollary 2.9