Mixed-precision algorithms for solving the Sylvester matrix equation

We consider the solution of the general Sylvester equation $AX+XB=C$ in mixed precision. First, we investigate the use of GMRES-based iterative refinement (GMRES-IR) to solve the equation using implicitly its Kronecker product form: we propose an efficient scheme to use the Schur factors of the coefficient matrices as preconditioners, but we demonstrate that this approach is not suitable in the case of the Sylvester equation. By revisiting a stationary iteration for linear systems, we therefore derive a new iterative refinement scheme for the quasi-triangular Sylvester equation, and our rounding error analysis provides sufficient conditions for convergence and a bound on the attainable relative residual. We leverage this iterative scheme to solve the general Sylvester equation in mixed precision. The new algorithms compute the Schur decomposition of the matrix coefficients in low precision, use the low-precision Schur factors to obtain an approximate solution to the quasi-triangular equation, and iteratively refine it to obtain a working-precision solution to the quasi-triangular equation. However, being only orthonormal to low precision, the unitary Schur factors of $A$ and $B$ cannot be used to recover the solution to the original equation. We propose two effective approaches to address this issue: one is based on re-orthonormalization in the working precision, and the other on explicit inversion of the almost-unitary factors. We test these mixed-precision algorithms on various Sylvester and Lyapunov equations from the literature. Our numerical experiments show that, for both classes of equations, the new algorithms are at least as accurate as existing ones. Our cost analysis, on the other hand, suggests that they would typically be faster than mono-precision alternatives if implemented on hardware that natively supports low precision.