pith. sign in

arxiv: 2510.02126 · v2 · pith:7RHEK5WTnew · submitted 2025-10-02 · 🧮 math.NA · cs.NA

Mixed-precision iterative refinement for low-rank Lyapunov equations

Peter Benner , Xiaobo Liu This is my paper

Pith reviewed 2026-05-18 10:17 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords low-rank Lyapunov equationsmixed precisioniterative refinementrounding error analysissign function methodresidual boundshalf precision
0
0 comments X

The pith

Mixed-precision iterative refinement lets half-precision inner solves handle low-rank Lyapunov equations with condition numbers up to 1/u without losing final accuracy.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper introduces a mixed-precision iterative refinement framework for low-rank Lyapunov equations of the form AX + XA^T + W = 0. Rounding error analysis yields conditions on attainable normwise residuals and rules for choosing algorithmic parameters across different precision combinations. These conditions hold for any inner solver that meets a prescribed residual target. The framework is illustrated with the sign-function Newton iteration, showing that half precision remains viable when the matrix condition number stays around the reciprocal of the half-precision unit roundoff.

Core claim

A mixed-precision iterative refinement scheme for low-rank Lyapunov equations supplies sufficient conditions, derived from rounding error analysis, under which the final residual accuracy is preserved even when inner linear solves run in reduced precision such as half precision; the conditions depend only on the chosen precisions and the requirement that each inner solve meets its own residual tolerance, not on the particular inner solver.

What carries the argument

Mixed-precision iterative refinement with derived residual bounds and explicit parameter-selection rules that remain valid whenever the inner solver attains the prescribed residual accuracy.

If this is right

  • Attainable normwise residuals can be bounded explicitly in terms of the working and inner precisions.
  • Algorithmic tolerances and iteration counts can be set in advance to guarantee a desired residual level.
  • The sign-function Newton iteration can safely use half precision when the Lyapunov operator has condition number of order 1/u_s.
  • The same residual guarantees apply to any inner solver that meets its local accuracy target.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same refinement strategy may transfer directly to low-rank Sylvester or Riccati equations.
  • Hardware that natively supports low-precision arithmetic could be used to solve larger-scale control problems without accuracy penalty.
  • An adaptive precision scheduler could be built by estimating the condition number once and then selecting inner precision accordingly.

Load-bearing premise

The inner solves must reach the prescribed residual accuracy for the derived overall residual bounds to remain valid.

What would settle it

An experiment in which an inner solver deliberately stops short of the target residual tolerance and the final computed solution shows a visibly larger residual or forward error than the predicted bound.

read the original 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. These conditions are independent of the choice of inner solver, provided that the prescribed residual accuracy is attained in the inner solves. Using the sign-function Newton iteration as the solver, we demonstrate that reduced precisions, such as half precision with unit roundoff $u_s$, can be used efficiently for Lyapunov equations with condition numbers of order $1/u_s$ without compromising the attainable solution quality. This provides an algorithmic framework towards exploiting native low-precision hardware to accelerate Lyapunov solvers without sacrificing accuracy.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 3 minor

Summary. The manuscript develops a mixed-precision iterative refinement framework for low-rank Lyapunov equations AX + XA^T + W = 0 (with W low-rank). Rounding error analysis derives sufficient conditions for attainable normwise residuals across precision settings and guidelines for algorithmic parameters. These conditions hold independently of the inner solver choice provided the inner solves attain their prescribed residual tolerances. The framework is demonstrated using the sign-function Newton iteration as inner solver, showing that reduced precisions (e.g., half precision with unit roundoff u_s) can be used for problems with condition numbers of order 1/u_s without loss of attainable solution quality.

Significance. If the analysis and the key assumption on inner-solve accuracy hold, the work offers a practical route to accelerate Lyapunov solvers on native low-precision hardware (e.g., GPUs or TPUs) for applications in control theory and model reduction, while preserving accuracy for moderately conditioned problems. The solver-independent formulation of the bounds is a notable strength that could extend to other inner methods.

major comments (2)
  1. [§3 (Rounding Error Analysis)] §3 (Rounding Error Analysis), around the derivation of residual bounds (e.g., the statements following Eq. (3.4) or (3.7)): The sufficient conditions for O(u) or O(u_s) residuals explicitly require that each inner solve attains its prescribed tolerance in the working precision. For the sign-function Newton iteration demonstrated in half precision, the Newton iterates themselves solve Sylvester equations whose conditioning mirrors that of the original Lyapunov operator; when cond ~ 1/u_s, accumulated rounding errors inside the low-precision inner loop may prevent the inner residual from being reached, rendering the outer mixed-precision bounds inapplicable precisely in the regime claimed to be advantageous. This assumption is load-bearing for the central claim and requires either an explicit bound on inner attainment or numerical verification that the tolerance is met.
  2. [§5 (Numerical Experiments)] §5 (Numerical Experiments), the half-precision runs with cond(A) ≈ 1/u_s: The reported residual histories and final accuracies should include explicit checks (or tabulated inner residuals) confirming that the prescribed inner tolerance was attained at each refinement step. Without this data, it remains unclear whether the observed outer accuracy validates the derived bounds or whether the inner solver benefited from implicit stabilization not captured by the analysis.
minor comments (3)
  1. [Introduction] The notation distinguishing full, single, and half precisions (u, u_s, etc.) and the precise definition of the low-rank factors L and S should be introduced earlier and used consistently from the abstract onward.
  2. [Figures] Figure captions for the convergence plots would benefit from stating the exact condition-number range and the inner tolerance used, to allow direct comparison with the theoretical thresholds.
  3. [Introduction] A short discussion of related mixed-precision work on Sylvester or Riccati equations would help situate the contribution; a few key references appear to be missing.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading and constructive feedback on our manuscript. We address each major comment below and describe the revisions planned for the next version.

read point-by-point responses

  1. Referee: §3 (Rounding Error Analysis), around the derivation of residual bounds (e.g., the statements following Eq. (3.4) or (3.7)): The sufficient conditions for O(u) or O(u_s) residuals explicitly require that each inner solve attains its prescribed tolerance in the working precision. For the sign-function Newton iteration demonstrated in half precision, the Newton iterates themselves solve Sylvester equations whose conditioning mirrors that of the original Lyapunov operator; when cond ~ 1/u_s, accumulated rounding errors inside the low-precision inner loop may prevent the inner residual from being reached, rendering the outer mixed-precision bounds inapplicable precisely in the regime claimed to be advantageous. This assumption is load-bearing for the central claim and requires either an explicit bound on inner attainment or numerical verification that the tolerance is met.

    Authors: We appreciate the referee highlighting the importance of the inner-solve attainment assumption. The analysis in Section 3 is explicitly conditional on the inner solver reaching the prescribed tolerance, as stated in the manuscript; the bounds are therefore solver-independent under this proviso. For the sign-function Newton iteration, a full a priori rounding-error bound on inner attainment in low precision would require a separate, detailed study of the iteration's behavior that lies outside the scope of the present work on the mixed-precision refinement framework. We will instead strengthen the manuscript by adding numerical verification of inner-residual attainment in the experiments (see response to the second comment). This directly addresses the referee's suggestion of numerical verification while preserving the general character of the analysis. revision: yes

  2. Referee: §5 (Numerical Experiments), the half-precision runs with cond(A) ≈ 1/u_s: The reported residual histories and final accuracies should include explicit checks (or tabulated inner residuals) confirming that the prescribed inner tolerance was attained at each refinement step. Without this data, it remains unclear whether the observed outer accuracy validates the derived bounds or whether the inner solver benefited from implicit stabilization not captured by the analysis.

    Authors: We agree that explicit confirmation of inner tolerance attainment is necessary to substantiate the numerical results. In the revised manuscript we will augment Section 5 with tabulated inner residuals (or additional plots) for all half-precision runs with cond(A) ≈ 1/u_s. These data will show the residuals achieved by the sign-function Newton iteration at each refinement step, allowing readers to verify that the prescribed tolerances were met and that the reported outer accuracies are consistent with the derived bounds. revision: yes

Circularity Check

0 steps flagged

No significant circularity; derivation relies on standard rounding-error analysis with explicit independent assumption

full rationale

The paper derives residual bounds and parameter choices via rounding error analysis for the mixed-precision iterative refinement applied to low-rank Lyapunov equations. The central claim that reduced precision (e.g., half precision with unit roundoff u_s) can be used for problems with cond ~ 1/u_s is conditioned on the explicit proviso that inner solves attain their prescribed residual tolerance. This assumption is stated as independent of the specific inner solver and is not derived from or equivalent to the target result; it is a standard precondition for the validity of the forward-error bounds. No self-definitional loops, fitted inputs renamed as predictions, or load-bearing self-citations that reduce the main result to its own inputs are present in the derivation chain. The analysis remains self-contained once the inner-solve accuracy condition is granted.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

Framework rests on standard numerical analysis assumptions about rounding errors and inner-solver accuracy; no explicit free parameters, invented entities, or non-standard axioms are stated in the abstract.

axioms (1)
  • domain assumption The inner solver attains the prescribed residual accuracy
    Explicitly required for the derived conditions and parameter choices to hold, independent of inner solver.

pith-pipeline@v0.9.0 · 5664 in / 1197 out tokens · 27465 ms · 2026-05-18T10:17:18.095978+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

38 extracted references · 38 canonical work pages

  1. [1]

    A. C. Antoulas, D. C. Sorensen, and Y. Zhou. On the decay rate of Hankel singular values and related issues.Sys. Control Lett., 46(5):323–342, 2002.doi:10.1016/S0167-6911(02) 00147-0

    work page doi:10.1016/s0167-6911(02 2002

  2. [2]

    Baur and P

    U. Baur and P. Benner. Gramian-based model reduction for data-sparse systems.SIAM J. Sci. Comput., 31(1):776–798, 2008.doi:10.1137/070711578

    work page doi:10.1137/070711578 2008

  3. [3]

    Benner, P

    P. Benner, P. Ezzatti, D. Kressner, E. S. Quintana-Ort´ ı, and A. Rem´ on. A mixed-precision algorithm for the solution of Lyapunov equations on hybrid CPU–GPU platforms.Parallel Comput., 37:439–450, 2011.doi:10.1016/j.parco.2010.12.002

    work page doi:10.1016/j.parco.2010.12.002 2011

  4. [4]

    Benner, J.-R

    P. Benner, J.-R. Li, and T. Penzl. Numerical solution of large-scale Lyapunov equations, Riccati equations, and linear-quadratic optimal control problems.Numer. Linear Algebra Appl., 15(9):755–777, 2008.doi:10.1002/nla.622

    work page doi:10.1002/nla.622 2008

  5. [5]

    Benner and E

    P. Benner and E. S. Quintana-Ort´ ı. Solving stable generalized Lyapunov equations with the matrix sign function.Numer. Algorithms, 20:75–100, 1999.doi:10.1023/A:1019191431273

    work page doi:10.1023/a:1019191431273 1999

  6. [6]

    Byers, C

    R. Byers, C. He, and V. Mehrmann. The matrix sign function method and the computation of invariant subspaces.SIAM J. Matrix. Anal. Appl., 18(3):615–632, 1997.doi:10.1137/ S0895479894277454

    work page 1997

  7. [7]

    Carson and N

    E. Carson and N. J. Higham. Accelerating the solution of linear systems by iterative refinement in three precisions.SIAM J. Sci. Comput., 40(2):A817–A847, 2018.doi: 10.1137/17M1140819

    work page doi:10.1137/17m1140819 2018

  8. [8]

    Chahlaoui and P

    Y. Chahlaoui and P. Van Dooren. A collection of benchmark examples for model reduction of linear time invariant dynamical systems. SLICOT Working Note 2002-2, Feb. 2002. URL: https://www.slicot.org/20-site/126-benchmark-examples-for-model-reduction

    work page 2002

  9. [9]

    J. W. Demmel.Applied Numerical Linear Algebra. Society for Industrial and Applied Math- ematics, Philadelphia, PA, USA, 1997.doi:10.1137/1.9781611971446

    work page doi:10.1137/1.9781611971446 1997

  10. [10]

    Dmytryshyn, M

    A. Dmytryshyn, M. Fasi, N. J. Higham, and X. Liu. Mixed-precision algorithms for solving the Sylvester matrix equation. ArXiv:2503.03456 [math.NA], Mar. 2025. URL:https://arxiv. org/abs/2503.03456

    work page arXiv 2025

  11. [11]

    Feitzinger, T

    F. Feitzinger, T. Hylla, and E. W. Sachs. Inexact Kleinman–Newton method for Riccati equations.SIAM J. Matrix. Anal. Appl., 31(2):272–288, 2009.doi:10.1137/070700978

    work page doi:10.1137/070700978 2009

  12. [12]

    G. H. Golub and C. F. Van Loan.Matrix Computations. Johns Hopkins University Press, Baltimore, MD, USA, 4th edition, 2013

    work page 2013

  13. [13]

    Gu and S

    M. Gu and S. C. Eisenstat. A divide-and-conquer algorithm for the symmetric tridiag- onal eigenproblem.SIAM J. Matrix. Anal. Appl., 16(1):172–191, 1995.doi:10.1137/ S0895479892241287. Preprint (Max Planck Institute for Dynamics of Complex Technical Systems, Magdeburg). 2025-10-03 P. Benner and X. Liu23

    work page 1995

  14. [14]

    Haidar, S

    A. Haidar, S. Tomov, and J. Dongarra. Towards half-precision computation for complex matrices: A case study for mixed precision solvers on GPUs. In2019 IEEE/ACM 10th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems (ScalA), pages 17–24, Denver, CO, USA, 2019. IEEE.doi:10.1109/ScalA49573.2019.00008

    work page doi:10.1109/scala49573.2019.00008 2019

  15. [15]

    In: 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp 409--419, doi:10.1109/CVPR.2018.00050

    A. Haidar, S. Tomov, J. Dongarra, and N. J. Higham. Harnessing GPU tensor cores for fast FP16 arithmetic to speed up mixed-precision iterative refinement solvers. InProceedings of the International Conference for High Performance Computing, Networking, Storage, and Analysis, SC18 (Dallas, TX), pages 47:1–47:11, Piscataway, NJ, USA, 2018. IEEE.doi: 10.1109...

    work page doi:10.1109/sc.2018.00050 2018

  16. [16]

    S. J. Hammarling. Numerical solution of the stable, non-negative definite Lyapunov equation. IMA J. Numer. Anal., 2(3):303–323, 1982.doi:10.1093/imanum/2.3.303

    work page doi:10.1093/imanum/2.3.303 1982

  17. [17]

    N. J. Higham. Computing the polar decomposition—with applications.SIAM J. Sci. Statist. Comput., 7(4):1160–1174, Oct. 1986.doi:10.1137/0907079

    work page doi:10.1137/0907079 1986

  18. [18]

    N. J. Higham.Accuracy and Stability of Numerical Algorithms. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, second edition, 2002.doi:10.1137/1. 9780898718027

    work page doi:10.1137/1 2002

  19. [19]

    N. J. Higham.Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2008.doi:10.1137/1.9780898717778

    work page doi:10.1137/1.9780898717778 2008

  20. [20]

    N. J. Higham and T. Mary. Mixed precision algorithms in numerical linear algebra.Acta Numerica, 31:347–414, May 2022.doi:10.1017/s0962492922000022

    work page doi:10.1017/s0962492922000022 2022

  21. [21]

    N. J. Higham and S. Pranesh. Simulating low precision floating-point arithmetic.SIAM J. Sci. Comput., 41(5):C585–C602, 2019.doi:10.1137/19M1251308

    work page doi:10.1137/19m1251308 2019

  22. [22]

    Huss–Lederman, E

    S. Huss–Lederman, E. S. Quintana-Ort´ ı, X. Sun, and Y. Y. Wu. Parallel spectral division using the matrix sign function for the generalized eigenproblem.Int. J. High Speed Comput., 11(1):1–14, 2000.doi:10.1142/S0129053300000084. [23]IEEE Standard for Floating-Point Arithmetic, IEEE Std754-2019 (Revision of IEEE754- 2008). The Institute of Electrical and ...

    work page doi:10.1142/s0129053300000084 2000

  23. [23]

    I. M. Jaimoukha and E. M. Kasenally. Krylov subspace methods for solving large Lyapunov equations.SIAM J. Numer. Anal., 31(1):227–251, 1994.doi:10.1137/0731012

    work page doi:10.1137/0731012 1994

  24. [24]

    Jbilou and A

    K. Jbilou and A. J. Riquet. Projection methods for large Lyapunov matrix equations.Linear Algebra Appl., 415(2):344–358, 2006. Special Issue on Order Reduction of Large-Scale Systems. doi:10.1016/j.laa.2004.11.004

    work page doi:10.1016/j.laa.2004.11.004 2006

  25. [25]

    Jing-Rebecca and J

    L. Jing-Rebecca and J. White. Low-rank solution of Lyapunov equations.SIAM Rev., 46(4):693–713, 2004.doi:10.1137/S0036144504443389

    work page doi:10.1137/s0036144504443389 2004

  26. [26]

    Komaroff

    N. Komaroff. Simultaneous eigenvalue lower bounds for the Lyapunov matrix equation.IEEE Trans. Automat. Control, 33(1):126–128, 1988.doi:10.1109/9.377

    work page doi:10.1109/9.377 1988

  27. [27]

    N. Lang, H. Mena, and J. Saak. On the benefits of theLDL T factorization for large-scale differential matrix equation solvers.Linear Algebra Appl., 480:44–71, 2015.doi:10.1016/j. laa.2015.04.006

    work page doi:10.1016/j 2015

  28. [28]

    V. B. Larin and F. A. Aliev. Construction of square root factor for solution of the Lyapunov matrix equation.Syst. Control Lett., 20(2):109–112, 1993.doi:10.1016/0167-6911(93) 90022-X

    work page doi:10.1016/0167-6911(93 1993

  29. [29]

    A. Laub, M. Heath, C. Paige, and R. Ward. Computation of system balancing transformations and other applications of simultaneous diagonalization algorithms.IEEE Trans. Automat. Control, 32(2):115–122, Feb. 1987.doi:10.1109/TAC.1987.1104549. Preprint (Max Planck Institute for Dynamics of Complex Technical Systems, Magdeburg). 2025-10-03 P. Benner and X. Liu24

    work page doi:10.1109/tac.1987.1104549 1987

  30. [30]

    Markidis, S

    S. Markidis, S. W. D. Chien, E. Laure, I. B. Peng, and J. S. Vetter. NVIDIA Tensor Core pro- grammability, performance & precision. In2018 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), pages 522–531, Los Alamitos, CA, USA, 2018. IEEE Computer Society.doi:10.1109/IPDPSW.2018.00091

    work page doi:10.1109/ipdpsw.2018.00091 2018

  31. [31]

    B. Moore. Principal component analysis in linear systems: Controllability, observability, and model reduction.IEEE Trans. Automat. Control, 26(1):17–32, 1981.doi:10.1109/TAC.1981. 1102568

    work page doi:10.1109/tac.1981 1981

  32. [32]

    T. Penzl. A cyclic low-rank smith method for large sparse Lyapunov equations.SIAM J. Sci. Comput., 21(4):1401–1418, 2000.doi:10.1137/S1064827598347666

    work page doi:10.1137/s1064827598347666 2000

  33. [33]

    Pernebo and L

    L. Pernebo and L. Silverman. Model reduction via balanced state space representations.IEEE Trans. Automat. Control, 27(2):382–387, 1982.doi:10.1109/TAC.1982.1102945

    work page doi:10.1109/tac.1982.1102945 1982

  34. [34]

    J. D. Roberts. Linear model reduction and solution of the algebraic Riccati equation by use of the sign function.Int. J. Control, 32(4):677–687, 1980.doi:10.1080/00207178008922881

    work page doi:10.1080/00207178008922881 1980

  35. [35]

    J. Schulze. A low-rank parareal solver for differential Riccati equations written in Julia. 2022. doi:10.5281/zenodo.7843198

    work page doi:10.5281/zenodo.7843198 2022

  36. [36]

    Sima and P

    V. Sima and P. Benner. Experimental evaluation of new SLICOT solvers for linear matrix equations based on the matrix sign function. In2008 IEEE Int Symposium on Computer- Aided Control System Design, Procceedings of the 2008 IEEE Multi-conference on Systems and Control, page 601–606. IEEE, 2008.doi:10.1109/CACSD.2008.4627361

    work page doi:10.1109/cacsd.2008.4627361 2008

  37. [37]

    Simoncini

    V. Simoncini. A new iterative method for solving large-scale Lyapunov matrix equations. SIAM J. Sci. Comput., 29(3):1268–1288, 2007.doi:10.1137/06066120X

    work page doi:10.1137/06066120x 2007

  38. [38]

    SIAM Rev.58, 377–441 (2016)

    V. Simoncini. Computational methods for linear matrix equations.SIAM Rev., 58(3):377–441, 2016.doi:10.1137/130912839. Preprint (Max Planck Institute for Dynamics of Complex Technical Systems, Magdeburg). 2025-10-03

    work page doi:10.1137/130912839 2016