Mixed-precision iterative refinement for low-rank Lyapunov equations
Peter Benner, Xiaobo Liu
TL;DR
This work develops a mixed-precision iterative refinement framework for solving large-scale, low-rank Lyapunov equations of the form $AX + XA^T + W = 0$ with $W=LL^T$ or $W=LSL^T$, where $A$ is Hurwitz and $W$ is low rank. It analyzes two factored IR variants, Cholesky-type and LDL^T-type, deriving rounding-error bounds that guide precision choices ($u$, $u_s$, $u_r$, $u_c$) and rank truncation strategies, and demonstrates the use of the sign function Newton iteration as the solver within the IR loop. The paper shows that reduced precisions such as bf16 or fp16 can accelerate the solution for moderately conditioned problems while preserving accuracy, provided the precision and truncation parameters satisfy the derived bounds. Numerical experiments on synthetic and benchmark problems validate the analysis and reveal substantial potential speedups in cache-bound scenarios, especially for well-conditioned instances, while also highlighting conditions under which low-precision solvers may fail to converge. These results open avenues for hardware-accelerated, mixed-precision Lyapunov solvers and motivate exploring other solver backends within the IR framework.
Abstract
We develop a mixed-precision iterative refinement framework for solving low-rank Lyapunov matrix equations $AX + XA^T + W =0$, where $W=LL^T$ or $W=LSL^T$. Via rounding error analysis of the algorithms we derive sufficient conditions for the attainable normwise residuals in different precision settings and show how the algorithmic parameters should be chosen. Using the sign function Newton iteration as the solver, we show that reduced precisions, such as the half precision, can be used as the solver precision (with unit roundoff $u_s$) to accelerate the solution of Lyapunov equations of condition number up to $1/u_s$ without compromising its quality.