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

Abstract: This paper considers large-scale nonsymmetric continuous-time algebraic Riccati equations (NAREs) that admit low-rank solutions. Low-rank alternating direction implicit (ADI) methods have proven to be an efficient approach for solving several matrix equations, including Lyapunov equations, Sylvester equations, and symmetric Riccati equations. Although a low-rank algorithm for the Sylvester equation has been used as an inner loop in computing low-rank solutions of NAREs, no low-rank ADI algorithm currently exists for NAREs themselves. This paper fills this gap by developing a low-rank ADI algorithm for large-scale NAREs that admit a low-rank solution. Since Lyapunov equations, Sylvester equations, and symmetric Riccati equations are special cases of the NARE, the existing low-rank ADI methods in the literature are special cases of the more general low-rank ADI method proposed here. An automatic and computationally efficient method for shift generation is also discussed, and a subspace-accelerated projection approach is presented to generate shifts for subsequent iterations without user intervention. Once initialized with arbitrary shifts, the proposed algorithm solves large-scale NAREs autonomously, generating its own shifts. Numerical results are presented using benchmark example of order $10^6$, demonstrating the computational efficiency and accuracy of the proposed algorithm.

• The paper introduces a low-rank ADI framework that unifies and advances solvers for large-scale non-symmetric algebraic Riccati equations.
• It employs projection-based techniques and automatic shift generation to achieve rapid convergence via adaptive, block Krylov subspace formulations.
• Numerical benchmarks on systems exceeding one million order demonstrate the method's scalability, efficiency, and applicability in control and model reduction.

Low-rank ADI Algorithm for Large-scale Non-symmetric Algebraic Riccati Equations

Problem Formulation and Context

The paper "A Low-rank ADI Algorithm for Solving Large-scale Non-symmetric Algebraic Riccati Equations" (2604.23208) develops specialized iterative techniques for computing low-rank solutions of large-scale NAREs of the form

$$A X \hat{E} + E X \hat{A} - E X \hat{B} C X \hat{E} + B \hat{C} = 0,$$

where $A$, $E$, $\hat{A}$, $\hat{E}$ are large sparse matrices, and the solution $X$ is expected to have low numerical rank. The NARE generalizes the class of matrix equations including Lyapunov, Sylvester, and symmetric algebraic Riccati equations, thus offering broad relevance across control, signal processing, and model reduction domains.

Existing low-rank ADI algorithms efficiently handle Lyapunov, Sylvester, and symmetric Riccati equations but do not directly apply to general NAREs. Previous strategies employed low-rank Sylvester solvers as inner iterations within alternate approaches for NAREs, but a unified ADI algorithm for NAREs had not been available. This paper closes that gap by developing the first low-rank ADI method for large-scale NAREs.

Theoretical Framework

The proposed algorithm synthesizes projection-based interpolatory ADI formulations, leveraging block Krylov subspaces and rational interpolation insights to construct low-rank approximations. At each iteration, the solution is factored as

$$X \approx \tilde{X}^{(i)} = V^{(i)} \bar{X}^{(i)} (\hat{W}^{(i)})^\top,$$

where $V^{(i)}$ and $\hat{W}^{(i)}$ are basis matrices generated via shifted linear solves, and $\bar{X}^{(i)}$ is recursively updated. The alternating direction implicit mechanism admits selection of complex-valued shifts $A$0, enabling accelerated convergence for stiff or poorly conditioned problems.

The algorithm constructs a sequence of subspaces such that

$A$1

with analogous definitions for $A$2. The residual is characterized by explicit closed forms, and block structures in $A$3, $A$4 allow simultaneous handling of real and complex shifts, including conjugate pairs.

A key result is the recursive expression for the residual:

$A$5

where residual factors are incrementally updated, providing a rigorous means to monitor convergence in the large-scale setting.

Algorithmic Developments

Two implementations are presented:

• The first uses the Sherman–Morrison–Woodbury (SMW) formula for efficient inversion of shifted linear systems with low-rank updates, thus reducing the computational cost associated with large $A$6.
• The second, termed UN-RADI, avoids SMW and exploits shared solves for Lyapunov equations, extracting low-rank factors for NAREs from the factors of Lyapunov problems. This enables simultaneous solving of NAREs and related equations, facilitating efficient computation in multi-query or multi-equation scenarios.

Both are equipped with adaptive shift strategies that eliminate manual tuning: the shifts are generated automatically using projected dominant pole estimation in Krylov subspaces, with implicit restart and subspace acceleration to maintain a manageable basis size. This parallels techniques from eigenvalue solvers and dominant pole algorithms.

The projection conditions and the parameterizations ensure that stabilizing solutions are obtained when required (i.e., closed-loop matrices are Hurwitz), and the method is flexible to produce general solutions when stabilization is not necessary.

Automatic Shift Generation

A crucial contribution is the automated shift generation. The ADI shifts are chosen according to controllability and observability metrics in projected subspaces:

• Projected eigen-decompositions identify poles with maximal $A$7, targeting modes most relevant for convergence.
• Shifts are restricted to $A$8 for improved interpolatory behavior, as demonstrated for related Sylvester equations.

The process draws from subspace-accelerated dominant pole estimation (SADPA) and ensures that deflation and implicit restart avoid unnecessary expansion of subspaces, keeping computational effort and storage in check.

Numerical Performance

Applying the algorithm to large-scale rail cooling models benchmark problems (with matrix order exceeding $A$9), the residual decays rapidly, achieving relative tolerance of $E$0 within 57 iterations and completing in 274.7s (for N-RADI) on commodity hardware. UN-RADI (SMW-free) converged in 376.7s. Both methods demonstrate strong scalability and accuracy; the automatic shift generation eliminates manual tuning and adapts to spectral features of the evolving residual.

Implications and Outlook

The low-rank ADI framework extends the reach of efficient solvers to large non-symmetric Riccati problems, encompassing Lyapunov, Sylvester, and symmetric Riccati cases as special instances. The theoretical results for residuals, projection conditions, and recursive gain updates provide robust footholds for further analysis in model reduction, control synthesis, and scalable numerical linear algebra.

Practically, the algorithm enables solution of previously intractable problems (e.g., optimal control for complex PDE systems) and offers a drop-in replacement for regular Riccati solvers in software frameworks. The SMW-free variant promises improved scaling for scenarios with large output dimensions.

Automated shift selection reduces reliance on expert tuning and supports batch solution of related equations. The method is amenable to further extensions, including block ADI, parallelization, and integration into unified frameworks for multi-equation model reduction.

Conclusion

The paper introduces a mathematically rigorous, computationally efficient low-rank ADI algorithm for large-scale NAREs. The method generalizes existing ADI schemes, provides automated shift generation, and demonstrates strong numerical performance on benchmark problems. Theoretical properties are validated, and MATLAB implementations are made available, supporting reproducibility and further experimentation. The developments pave the way for scalable solution of high-dimensional Riccati-type equations across modern applications in control, design, and data-driven modeling.