A Mixed Precision Eigensolver Based on the Jacobi Algorithm

The classic method for computing the spectral decomposition of a real symmetric matrix, the Jacobi algorithm, can be accelerated by using mixed precision arithmetic. The Jacobi algorithm is aiming to reduce the off-diagonal entries iteratively using Givens rotations. We investigate how to use the low precision to speed up this algorithm based on the approximate spectral decomposition in low precision. We first study two different index choosing techniques, classical and cyclic-by-row, for the Jacobi algorithm. Numerical testing suggests that cyclic-by-row is more efficient. Then we discuss two different methods of orthogonalizing an almost orthogonal matrix: the QR factorization and the polar decomposition. For polar decomposition, we speed up the Newton iteration by using the one-step Schulz iteration. Based on numerical testing, using the polar decomposition approach (Newton--Schulz iteration) is not only faster but also more accurate than using the QR factorization. A mixed precision algorithm for computing the spectral decomposition of a real symmetric matrix at double precision is provided. In doing so we compute the approximate eigenvector matrix $Q_\ell$ of $A$ in single precision using $\texttt{eig}$ and $\texttt{single}$ in MATLAB. We then use the Newton--Schulz iteration to orthogonalize the eigenvector matrix $Q_\ell$ into an orthogonal matrix $Q_d$ in double precision. Finally, we apply the cyclic-by-row Jacobi algorithm on $Q_d^TAQ_d$ and obtain the spectral decomposition of $A$. At this stage, we will see, from the testings, the cyclic-by-row Jacobi algorithm only need less than 10 iterations to converge by utilizing the quadratic convergence. The new mixed precision algorithm requires roughly 30\% of the time used by the Jacobi algorithm on its own.