pith. sign in

arxiv: 2604.06556 · v3 · submitted 2026-04-08 · 🧮 math.NA · cs.NA· cs.SY· eess.SY

LDL^top Factorization-based Generalized Low-rank ADI Algorithm for Solving Large-scale Algebraic Riccati Equations

Umair Zulfiqar This is my paper

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

classification 🧮 math.NA cs.NAcs.SYeess.SY
keywords algebraic Riccati equationslow-rank approximationADI methodLDL^T factorizationlarge-scale systemsnumerical algorithmscontrol theory
0
0 comments X

The pith

The generalized RADI algorithm with LDL^T factorization efficiently computes low-rank solutions to general algebraic Riccati equations.

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

This paper develops a new solver for large-scale general-form algebraic Riccati equations by extending the low-rank ADI method with an LDL^T factorization step. The approach allows handling equations that arise in state estimation and controller design, which previous RADI methods could not address directly. It bases the computation on a low-rank Cholesky factor ADI method to avoid certain matrix operations and includes a strategy for generating ADI shifts automatically. Tests on systems with millions of variables show it produces accurate results quickly.

Core claim

By integrating LDL^T factorization into the low-rank alternating direction implicit iteration and using the low-rank Cholesky factor ADI method as the underlying solver, the algorithm computes low-rank factors for the solutions of general Riccati equations without relying on the Sherman-Morrison-Woodbury formula.

What carries the argument

LDL^T factorization combined with low-rank Cholesky factor ADI iteration, which computes the low-rank solution factors for general-form Riccati equations.

If this is right

  • It handles general Riccati equations from applications like state estimation and controller design.
  • It avoids the Sherman-Morrison-Woodbury formula through use of the low-rank Cholesky ADI base.
  • Automatic generation of ADI shifts improves efficiency.
  • Numerical performance is demonstrated on Riccati equations of orders from 10^6 to 10^7.

Where Pith is reading between the lines

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

  • The method could extend to other large-scale matrix equations in control theory.
  • Integration with modern computing architectures might further speed up real-time applications.
  • Comparison with other low-rank solvers could reveal trade-offs in accuracy for specific problem classes.

Load-bearing premise

The general-form Riccati equations must admit low-rank solutions that the combined LDL^T and ADI procedure can capture accurately without numerical instability.

What would settle it

Running the algorithm on a general Riccati equation with a known low-rank solution and observing if the computed factors yield a residual norm much larger than machine epsilon would disprove its reliability.

Figures

Figures reproduced from arXiv: 2604.06556 by Umair Zulfiqar.

Figure 1
Figure 1. Figure 1: Normalized residual for X 3.2. RLC Ladder This example considers the RLC network model from [24]. It consists of 5 × 106 segments of the RLC ladder described in [24], resulting in a large-scale state-space model with the following dimensions: E ∈ R 107×107 , A ∈ R 107×107 , B ∈ R 107×1 , C ∈ R 1×107 , and D ∈ R 1×1 . For implicit restart in the proposed shift generation strategy, the maximum number of colu… view at source ↗
Figure 2
Figure 2. Figure 2: Normalized residual for Qricc [PITH_FULL_IMAGE:figures/full_fig_p012_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Normalized residual for Q′ ricc This can be solved using UG-RADI by setting B1 = B, B2 = [ ], R1 = −(I − D⊤D), R2 = [ ], C1 = C, Z = I, and C2 = D⊤C. UG-RADI converged in 15 iterations in 108.3384 seconds. The decay in normalized residual is plotted in [PITH_FULL_IMAGE:figures/full_fig_p012_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Normalized residual for Qpr [PITH_FULL_IMAGE:figures/full_fig_p013_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Normalized residual for Qbr Let QLQG solve the following CARE A ⊤QLQGE + E ⊤QLQGA − (B ⊤QLQGE + D⊤C) ⊤(R + D⊤D) −1 (B ⊤QLQGE + D⊤C) + C ⊤QC = 0. 13 [PITH_FULL_IMAGE:figures/full_fig_p013_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Normalized residual for QLQG Let Q∞ solve the following CARE A ⊤Q∞E + E ⊤Q∞A − E ⊤Q∞  BR−1B ⊤ − 1 γ 2 BB⊤  Q∞E + C ⊤QC = 0. Set γ = 1.5. This can be solved using G-RADI by setting B1 = B, B2 = q 1 γ2 B, R1 = R, R2 = I, C1 = C, Z = Q, and C2 = [ ]. UG-RADI converged in 16 iterations in 82.2031 seconds. The decay in normalized residual is plotted in [PITH_FULL_IMAGE:figures/full_fig_p014_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Normalized residual for Q∞ Appendix A: Proof of Theorem 2.2 Proof. 1. Post-multiplying (11) by X˜(i) yields −(S (i) w ) ⊤ − (X(i) ) −1S (i) w X˜(i) + (L (i) w ) ⊤ZˆCˆ(i) r + Bˆ(i) r Rˆ−1 (Bˆ(i) r ) ⊤X˜(i) = 0 Aˆ r = −(X(i) ) −1S (i) w X˜(i) + Bˆ(i) r Rˆ−1 (Bˆ(i) r ) ⊤X˜(i) Aˆ (i) cl = −(X(i) ) −1S (i) w X˜(i) . Hence, Aˆ (i) cl has eigenvalues α1, . . . , αk, each with multiplicity p + m1. Now consider (Aˆ… view at source ↗
read the original abstract

The low-rank alternating direction implicit (ADI) method is an efficient and effective solver for large-scale standard continuous-time algebraic Riccati equations that admit low-rank solutions. However, the existing low-rank ADI algorithm for Riccati equations (RADI) cannot be directly applied to general-form Riccati equations. This paper introduces a generalized RADI algorithm based on an $LDL^\top$ factorization, which efficiently handles the general Riccati equations arising in important applications like state estimation and controller design. An efficient implementation is presented that avoids the Sherman-Morrison-Woodbury formula and instead uses a low-rank Cholesky factor ADI method as the base algorithm to compute low-rank factors of general-form Riccati equations. Sample MATLAB-based implementations of the proposed algorithm are also provided. An approach for automatically and efficiently generating ADI shifts is discussed. Numerical examples solving several Riccati equations of orders ranging from $10^6$ to $10^7$ accurately and efficiently are presented, demonstrating the effectiveness of the proposed algorithm.

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 / 2 minor

Summary. The paper introduces a generalized low-rank ADI (RADI) algorithm for large-scale general-form continuous-time algebraic Riccati equations. It employs an LDL^T factorization step to enable application of an existing low-rank Cholesky-factor ADI solver as the base method, avoiding the Sherman-Morrison-Woodbury formula. The work includes an automatic ADI shift generation strategy, provides sample MATLAB implementations, and reports numerical results on Riccati equations of orders 10^6 to 10^7.

Significance. If the stability and efficiency claims hold, the method would extend low-rank ADI techniques to general Riccati equations arising in state estimation and controller design, offering a practical tool for large-scale problems in control theory. The provision of reproducible MATLAB code and an automatic shift strategy are strengths that support usability and verification.

major comments (2)
  1. [Algorithm section] Algorithm derivation (likely §3): the central construction assumes that the LDL^T factorization step preserves the low-rank structure and does not amplify rounding errors when combined with the low-rank Cholesky ADI iteration, but no explicit a-priori bound on the perturbation to the ADI residual or shift computation is provided, especially for cases where the diagonal factor D has negative eigenvalues.
  2. [Numerical examples] Numerical examples section: the reported results for matrix orders 10^6–10^7 cite accurate and efficient solves but omit condition numbers of the factors, residual growth under controlled rounding, or direct comparisons against SMW-based alternatives, leaving the stability claim for indefinite or ill-conditioned general Riccati equations unverified.
minor comments (2)
  1. [Abstract] Abstract: the phrasing 'accurately and efficiently' would benefit from a brief quantitative summary of achieved residuals or timings to ground the claim.
  2. [Shift generation] The description of the automatic ADI shift generation could include a short pseudocode or reference to the specific heuristic used for clarity.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We are grateful to the referee for the thorough review and valuable comments on our manuscript. The feedback highlights important aspects regarding the theoretical stability of the proposed algorithm and the verification in numerical experiments. We address each major comment below and indicate the revisions made to the manuscript.

read point-by-point responses

  1. Referee: [Algorithm section] Algorithm derivation (likely §3): the central construction assumes that the LDL^T factorization step preserves the low-rank structure and does not amplify rounding errors when combined with the low-rank Cholesky ADI iteration, but no explicit a-priori bound on the perturbation to the ADI residual or shift computation is provided, especially for cases where the diagonal factor D has negative eigenvalues.

    Authors: We thank the referee for pointing out this aspect of the theoretical analysis. The LDL^T factorization is performed exactly in exact arithmetic and the low-rank structure is maintained by the subsequent low-rank Cholesky-factor ADI method as described in Section 3. However, we agree that no explicit a-priori perturbation bound is derived for the combined process, particularly when the diagonal matrix D has negative entries. Such a bound would require a detailed rounding error analysis that goes beyond the scope of the current work. In the revised manuscript, we have added a paragraph in the algorithm section discussing the numerical stability considerations and advising users to monitor the residual norm for cases with potentially indefinite factors. We believe this addresses the concern practically while noting the limitation in the theoretical guarantees. revision: partial

  2. Referee: [Numerical examples] Numerical examples section: the reported results for matrix orders 10^6–10^7 cite accurate and efficient solves but omit condition numbers of the factors, residual growth under controlled rounding, or direct comparisons against SMW-based alternatives, leaving the stability claim for indefinite or ill-conditioned general Riccati equations unverified.

    Authors: The referee correctly identifies that the numerical section focuses on demonstrating efficiency and accuracy for large-scale problems but does not include condition numbers of the factors or explicit residual growth analysis under controlled rounding. Direct comparisons with SMW-based methods are challenging for the largest problems due to computational cost, but we have performed such comparisons on medium-scale instances in the revised manuscript. We have also added tables reporting the condition numbers of the computed low-rank factors and included additional figures showing the evolution of the residual norm. For the very large problems, a full controlled rounding experiment is not feasible within reasonable computational resources, but the results in double precision show consistent convergence without apparent instability. We have updated the numerical examples section accordingly. revision: partial

Circularity Check

0 steps flagged

Direct algorithmic construction from LDL^T factorization and prior ADI solver; no reduction to self-inputs

full rationale

The paper describes an explicit algorithmic extension: it applies an LDL^T factorization to rewrite the general Riccati equation so that the existing low-rank Cholesky-factor ADI iteration (RADI) can be invoked as a black-box subroutine. No equation in the derivation equates the output low-rank factors to a fitted parameter, a self-referential definition, or a quantity computed only from the paper's own fitted values. The cited RADI base algorithm is treated as an independent, previously published method rather than a result derived inside this manuscript. Numerical examples and shift-selection heuristics are presented as implementation details, not as load-bearing premises that collapse back into the algorithm definition. Consequently the claimed efficiency and applicability rest on standard matrix-algebraic identities and external solver behavior, not on circular re-labeling of inputs.

Axiom & Free-Parameter Ledger

1 free parameters · 2 axioms · 0 invented entities

The algorithm relies on standard assumptions from numerical linear algebra (existence of low-rank solutions, stability of ADI iteration under chosen shifts) and introduces no new postulated entities. The automatic shift generation is described as an approach but its exact parameters are not detailed in the abstract.

free parameters (1)
  • ADI shifts
    The paper discusses an approach for automatically generating ADI shifts; these are algorithmic parameters chosen during execution rather than fitted constants.
axioms (2)
  • domain assumption The Riccati equation admits a low-rank solution that can be represented by low-rank factors.
    Stated in the abstract as the setting in which the low-rank ADI method is efficient.
  • domain assumption The low-rank Cholesky factor ADI method can serve as a stable base solver when combined with LDL^T factorization.
    Implicit in the claim that the generalized algorithm reuses the base method without the Sherman-Morrison-Woodbury formula.

pith-pipeline@v0.9.0 · 5496 in / 1518 out tokens · 40935 ms · 2026-05-10T18:35:42.595819+00:00 · methodology

discussion (0)

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

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score.

  1. A Low-rank ADI Algorithm for Solving Large-scale Non-symmetric Algebraic Riccati Equations

    math.NA 2026-04 unverdicted novelty 7.0

    Develops the first low-rank ADI algorithm for non-symmetric algebraic Riccati equations, with autonomous shift generation, and demonstrates it on a benchmark problem of order 10^6.

Reference graph

Works this paper leans on

24 extracted references · 24 canonical work pages · cited by 1 Pith paper

  1. [1]

    Y. Lin, V. Simoncini, A new subspace iteration method for the algebraic Riccati equation, Numerical Linear Algebra with Applications 22 (1) (2015) 26–47

    work page 2015

  2. [2]

    Benner, Z

    P. Benner, Z. Bujanović, P. Kürschner, J. Saak, RADI: A low-rank ADI-type algorithm for large scale algebraic Riccati equations, Numerische Mathematik 138 (2) (2018) 301–330

    work page 2018

  3. [3]

    Bertram, H

    C. Bertram, H. Faßbender, On a family of low-rank algorithms for large-scale algebraic Riccati equa- tions, Linear Algebra and Its Applications 687 (2024) 38–67. 30

    work page 2024

  4. [4]

    Benner, J

    P. Benner, J. Heiland, S. W. Werner, A low-rank solution method for Riccati equations with indefinite quadratic terms, Numerical Algorithms 92 (2) (2023) 1083–1103

    work page 2023

  5. [5]

    J. Saak, S. W. Werner, UsingLDL⊤ factorizations in Newton’s method for solving general large-scale algebraic Riccati equations, Electronic Transactions on Numerical Analysis 62 (2024) 95–118

    work page 2024

  6. [6]

    Zulfiqar, Z.-Y

    U. Zulfiqar, Z.-Y. Huang, Q.-Y. Song, Z.-Y. Gao, A unified low-rank ADI framework with shared linear solves for simultaneously solving multiple Lyapunov, Sylvester, and Riccati equations, arXiv preprint arXiv:2512.04676 (2025)

    work page arXiv 2025

  7. [7]

    K.Gallivan, A.Vandendorpe, P.VanDooren, Sylvesterequationsandprojection-basedmodelreduction, Journal of Computational and Applied Mathematics 162 (1) (2004) 213–229

    work page 2004

  8. [8]

    Wolf,H 2 pseudo-optimal model order reduction, Ph.D

    T. Wolf,H 2 pseudo-optimal model order reduction, Ph.D. thesis, Technische Universität München (2014)

    work page 2014

  9. [9]

    H. K. Panzer, Model order reduction by Krylov subspace methods with global error bounds and auto- matic choice of parameters, Ph.D. thesis, Technische Universität München (2014)

    work page 2014

  10. [10]

    Astolfi, Model reduction by moment matching for linear and nonlinear systems, IEEE Transactions on Automatic Control 55 (10) (2010) 2321–2336

    A. Astolfi, Model reduction by moment matching for linear and nonlinear systems, IEEE Transactions on Automatic Control 55 (10) (2010) 2321–2336

    work page 2010

  11. [11]

    G. H. Golub, C. F. Van Loan, Matrix computations, JHU press, 2013

    work page 2013

  12. [12]

    Benner, P

    P. Benner, P. Kürschner, J. Saak, A reformulated low-rank ADI iteration with explicit residual factors, PAMM 13 (1) (2013) 585–586

    work page 2013

  13. [13]

    Gugercin, A

    S. Gugercin, A. C. Antoulas, C. Beattie,H2 model reduction for large-scale linear dynamical systems, SIAM Journal on Matrix Analysis and Applications 30 (2) (2008) 609–638

    work page 2008

  14. [14]

    Saak, Efficient numerical solution of large scale algebraic matrix equations in PDE control and model order reduction, Ph.D

    J. Saak, Efficient numerical solution of large scale algebraic matrix equations in PDE control and model order reduction, Ph.D. thesis, TU Chemnitz (2009)

    work page 2009

  15. [15]

    Benner, P

    P. Benner, P. Kürschner, J. Saak, Self-generating and efficient shift parameters in ADI methods for large Lyapunov and Sylvester equations, Electronic Transactions on Numerical Analysis 43 (2014) 142–162

    work page 2014

  16. [16]

    Rommes, Methods for eigenvalue problems with applications in model order reduction, Ph.D

    J. Rommes, Methods for eigenvalue problems with applications in model order reduction, Ph.D. thesis, Utrecht University (2007)

    work page 2007

  17. [17]

    Rommes, Modal approximation and computation of dominant poles, in: Model Order Reduction: Theory, Research Aspects and Applications, Springer, 2008, pp

    J. Rommes, Modal approximation and computation of dominant poles, in: Model Order Reduction: Theory, Research Aspects and Applications, Springer, 2008, pp. 177–193

    work page 2008

  18. [18]

    Martins, L

    N. Martins, L. T. Lima, H. J. Pinto, Computing dominant poles of power system transfer functions, IEEE Transactions on Power Systems 11 (1) (1996) 162–170

    work page 1996

  19. [19]

    E.Mengi, Large-scaleestimationofdominantpolesofatransferfunctionbyaninterpolatoryframework, SIAM Journal on Scientific Computing 44 (4) (2022) A2412–A2438

    work page 2022

  20. [20]

    Saad, Numerical methods for large eigenvalue problems: Revised edition, SIAM, 2011

    Y. Saad, Numerical methods for large eigenvalue problems: Revised edition, SIAM, 2011

    work page 2011

  21. [21]

    U. Zulfiqar, MATLAB codes forLDL ⊤ factorization-based generalized low-rank ADI algorithm for solving large-scale algebraic riccati equations, https://doi.org/10.5281/zenodo.19568180 (2026). URLhttps://doi.org/10.5281/zenodo.19568180

    work page doi:10.5281/zenodo.19568180 2026

  22. [22]

    Benner, J

    P. Benner, J. Saak, A semi-discretized heat transfer model for optimal cooling of steel profiles, in: Di- mension Reduction of Large-Scale Systems: Proceedings of a Workshop held in Oberwolfach, Germany, October 19–25, 2003, Springer, 2005, pp. 353–356. 31

    work page 2003

  23. [23]

    J. Saak, M. Köhler, P. Benner, M-M.E.S.S-2.0 - The Matrix Equation Sparse Solver Library, Zenodo (2021)

    work page 2021

  24. [24]

    Gugercin, A

    S. Gugercin, A. C. Antoulas, A survey of model reduction by balanced truncation and some new results, International Journal of Control 77 (8) (2004) 748–766. 32

    work page 2004