Low-Rank Lyapunov ADI Methods

• Low-Rank Lyapunov ADI is a factorized approach for solving large-scale continuous-time Lyapunov equations using sparse shifted linear solves and low-rank approximations.
• It employs rational Krylov subspace methods and interpolation techniques to achieve efficient approximation and monitor convergence via exact low-rank residual factorizations.
• Modern extensions of LR-ADI include inexact solves, adaptive shift selection, warm starts, and mixed precision strategies that accelerate iterative performance on large-scale problems.

Searching arXiv for recent and foundational papers on LR-ADI, shift selection, inexact solves, projection interpretations, and recent extensions. Low-Rank Lyapunov ADI (LR-ADI) is the low-rank realization of the alternating direction implicit method for large-scale continuous-time Lyapunov equations, especially when the coefficient matrices are sparse and the constant term has low rank. Instead of forming the dense solution matrix explicitly, LR-ADI computes a factorized approximation, typically $X \approx ZZ^*$ in the positive semidefinite case or, more generally, an $LDL^T$-type representation, so that storage and arithmetic scale with the numerical rank rather than with $n^2$. Modern analysis places LR-ADI simultaneously in the traditions of stationary matrix iterations, rational Krylov projection, and interpolation-based model reduction, rather than treating it as only a fixed-point recurrence (Wolf et al., 2013, Flagg et al., 2012).

1. Problem class and algebraic setting

The standard Lyapunov equation treated by LR-ADI is

$AX+XA^T=-BB^T,$

with $A\in\mathbb R^{n\times n}$ stable and $B\in\mathbb R^{n\times r}$, $r\ll n$. A generalized form,

$A P E^T + E P A^T + B B^T = 0,$

is equally central in descriptor-system settings, under the assumptions that $E$ is nonsingular and the pencil $(A,E)$ is asymptotically stable (Wolf et al., 2013, Kürschner et al., 2018). In both cases the exact solution is typically dense, while the right-hand side is low rank and the solution often has rapidly decaying singular values, which motivates approximations of the form

$LDL^T$0

This factor viewpoint is the defining feature of LR-ADI. The method never needs to assemble the full dense solution, and its dominant operations are sparse shifted linear solves with a small number of right-hand sides. In positive semidefinite settings the factorization is usually Cholesky-like, but several extensions use symmetric-indefinite forms $LDL^T$1 or $LDL^T$2 when the initial value is nonzero, the right-hand side is indefinite, or the formulation is generalized (Schulze et al., 2024, Smith et al., 4 Dec 2025).

A source of notational confusion is sign convention. One common formulation uses shifts $LDL^T$3 and solves with $LDL^T$4; another uses $LDL^T$5 and solves with $LDL^T$6. These are equivalent parameterizations of the same basic mechanism, and the literature moves freely between them (Wolf et al., 2013, Kürschner et al., 2018).

2. Core LR-ADI recurrence and residual factors

In generalized form, a standard LR-ADI block recurrence is

$LDL^T$7

$LDL^T$8

after which the approximation is $LDL^T$9 with $n^2$0. In the alternative residual-factor form one solves

$n^2$1

and appends $n^2$2 to the low-rank factor (Wolf et al., 2013, Kürschner et al., 2018).

The residual factorization is one of LR-ADI’s most important computational properties. For the generalized Lyapunov residual

$n^2$3

the exact LR-ADI iterate satisfies

$n^2$4

so $n^2$5, where $n^2$6 is the number of columns of $n^2$7. Hence

$n^2$8

which reduces residual-norm evaluation to an $n^2$9 problem. In residual-factor notation this is the familiar identity $AX+XA^T=-BB^T,$0, again showing that the residual rank is bounded by the right-hand-side rank (Wolf et al., 2013, Kürschner et al., 2018).

This residual structure is often mistaken for a secondary implementation detail. It is not. It is the mechanism that makes LR-ADI viable at scale: residual monitoring, stopping criteria, and many adaptive strategies rely on the fact that both the approximation and the residual remain available in low-rank factored form.

3. Rational Krylov, projection, and $AX+XA^T=-BB^T,$1 pseudo-optimality

A central structural result is that LR-ADI iterates live in rational Krylov subspaces. For Sylvester and Lyapunov equations with low-rank right-hand side, the ADI factors satisfy

$AX+XA^T=-BB^T,$2

so the computed basis is not ad hoc: it is a rational block Krylov basis (Wolf et al., 2013).

This observation leads to a precise projection interpretation. For arbitrary shift sets closed under conjugation, there always exists an oblique test basis $AX+XA^T=-BB^T,$3 such that the LR-ADI approximation coincides with the solution of a projected Lyapunov equation,

$AX+XA^T=-BB^T,$4

with $AX+XA^T=-BB^T,$5 spanning the same rational block Krylov space as the ADI basis. Thus LR-ADI is always an implicitly oblique rational Krylov projection method, not only for special shifts. The associated reduced explicit system is unique, and the reduced model attached to the ADI iteration is $AX+XA^T=-BB^T,$6 pseudo-optimal within the class of systems having the same reduced pole structure (Wolf et al., 2013).

A sharper statement holds for special shifts. For rank-one Sylvester and Lyapunov equations, if the shifts are chosen so that the reduced poles are the mirror images of the interpolation points,

$AX+XA^T=-BB^T,$7

then the ADI approximation and the rational Krylov projection approximation coincide exactly. In the Lyapunov case, the corresponding residual satisfies the orthogonality condition

$AX+XA^T=-BB^T,$8

These shifts are called pseudo-$AX+XA^T=-BB^T,$9-optimal because they enforce the interpolation condition associated with $A\in\mathbb R^{n\times n}$0-optimal model reduction, but not the derivative-matching condition required for full $A\in\mathbb R^{n\times n}$1 optimality (Flagg et al., 2012).

One common misconception is that LR-ADI and projection methods are merely heuristically related. The projection interpretation is exact: always oblique in the general case, and orthogonal under the mirror-pole condition above (Wolf et al., 2013, Flagg et al., 2012).

4. Shift parameters and convergence behavior

Shift selection governs LR-ADI far more strongly than the outer recurrence itself. The same update formulas can converge rapidly or stagnate depending on the poles. Pseudo-$A\in\mathbb R^{n\times n}$2-optimal shifts are one theoretically distinguished class: they make ADI coincide with rational Krylov projection and, in the Lyapunov case, enforce residual orthogonality. Numerically, they also produce nearly best rank-$A\in\mathbb R^{n\times n}$3 approximations in spectral norm on several benchmark problems. In the EADY example, for instance, the best rank-20 benchmark was $A\in\mathbb R^{n\times n}$4, while the pseudo-$A\in\mathbb R^{n\times n}$5 approximation gave $A\in\mathbb R^{n\times n}$6; at rank 30 the values were $A\in\mathbb R^{n\times n}$7 and $A\in\mathbb R^{n\times n}$8, respectively (Flagg et al., 2012).

A different line of work selects each shift by minimizing the next residual norm. If $A\in\mathbb R^{n\times n}$9 is the current residual factor, the ideal next shift minimizes

$B\in\mathbb R^{n\times r}$0

but exact evaluation is too expensive because it would require repeated large shifted solves inside a nonlinear optimization loop. The practical solution is to optimize compressed objective functions built from small subspaces already generated during LR-ADI, such as recent columns of the low-rank factor or small extended Krylov spaces. On nonsymmetric examples these residual-minimizing shifts consistently outperform precomputed heuristics and usually outperform competing dynamic strategies, though at the price of more elaborate shift generation (Kürschner, 2018).

The paper on adaptive extended rational Krylov approximation adds a related reduced-model shift rule,

$B\in\mathbb R^{n\times r}$1

where $B\in\mathbb R^{n\times r}$2 is the projected matrix and $B\in\mathbb R^{n\times r}$3 is extracted from the last row of the reduced Lyapunov solution. Although this result is derived for a rational Krylov method rather than LR-ADI itself, it reinforces the broader principle that useful shifts can be extracted from the current reduced dynamics rather than from static spectral estimates alone (Kolesnikov et al., 2014).

Another misconception is that shift selection can be based solely on eigenvalue enclosure. For nonsymmetric problems, the literature explicitly criticizes purely spectral heuristics for ignoring $B\in\mathbb R^{n\times r}$4, eigenvector information, and the current residual state (Kürschner, 2018).

5. Inexact solves, warm starts, integrated Krylov-ADI, and mixed precision

For very large problems, the shifted systems inside LR-ADI are often solved iteratively rather than by sparse direct factorization. In that setting the outer iteration becomes inexact. The exact residual is no longer the computed surrogate $B\in\mathbb R^{n\times r}$5; instead, a structured residual gap appears, driven by the accumulated inner linear-system residuals. This leads to a precise inexactness model and to practical relaxation rules for the inner stopping tolerance. The preferred adaptive rule uses the computed residual and an accumulated residual-gap estimate to loosen the inner tolerances as the outer iteration progresses. Across the reported examples, relaxed inner solves reduced inner iteration counts by roughly $B\in\mathbb R^{n\times r}$6 to $B\in\mathbb R^{n\times r}$7, with the second practical strategy consistently better than the first (Kürschner et al., 2018).

The traditional requirement $B\in\mathbb R^{n\times r}$8 has also been removed. By factorizing the initial residual of an arbitrary low-rank initial guess $B\in\mathbb R^{n\times r}$9 as

$r\ll n$0

with fixed, generally indefinite $r\ll n$1, the standard low-rank increment and residual recurrences extend to nonzero starts. This matters in outer iterations such as Newton methods for algebraic Riccati equations and Rosenbrock schemes for differential Riccati equations, where warm starts from the previous outer iterate can drastically reduce the number of inner LR-ADI steps. Reported speed-ups were about $r\ll n$2 for the algebraic Riccati setting and up to $r\ll n$3 for the differential Riccati setting (Schulze et al., 2024).

A different acceleration route is to merge LR-ADI with an extended Krylov solver for the shifted systems. Because the right-hand sides $r\ll n$4 of all ADI inner solves stay in a single growing extended Krylov space $r\ll n$5, one can reuse one basis for all shifted systems and perform the LR-ADI iteration in projected coordinates. On a class of standard Lyapunov equations this integrated Krylov-ADI approach achieved up to about $r\ll n$6 speed-up over direct-solver-based LR-ADI, while retaining compatibility with modern shift strategies (Benner et al., 2022).

Mixed precision addresses memory rather than iteration count. A recent study split LR-ADI storage and arithmetic across three precision classes: the accumulated factor $r\ll n$7, the propagated factors $r\ll n$8, and the inner factors $r\ll n$9. The most successful compromise was to store only $A P E^T + E P A^T + B B^T = 0,$0 in single precision while keeping the shifted solves and residual propagation in double precision, denoted ADI(S,D,D). Lowering the precision of $A P E^T + E P A^T + B B^T = 0,$1 and the shifted solves caused much larger deterioration in explicit residuals, while keeping only $A P E^T + E P A^T + B B^T = 0,$2 in double precision provided essentially no benefit. A particularly important warning is that the implicit LR-ADI residual can remain excellent even when the explicit residual has stagnated badly in finite precision (Schulze et al., 1 Aug 2025).

6. Indefinite, tangential, and unified extensions

Classical block LR-ADI becomes expensive when the right-hand side has large rank or is indefinite. For equations of the form

$A P E^T + E P A^T + B B^T = 0,$3

with Hermitian indefinite $A P E^T + E P A^T + B B^T = 0,$4, a tangential reformulation replaces block updates by rank-1 updates

$A P E^T + E P A^T + B B^T = 0,$5

with tangential directions restricted to eigenvectors of $A P E^T + E P A^T + B B^T = 0,$6. This reduces per-step rank growth from the full right-hand-side rank to one column, while preserving the low-rank residual structure. The reported numerical results show especially large savings when the constant term has high rank: in one benchmark the tangential approximation had dimension about $A P E^T + E P A^T + B B^T = 0,$7 versus about $A P E^T + E P A^T + B B^T = 0,$8 for block ADI (Smith et al., 4 Dec 2025).

Recent work has also unified Lyapunov, Sylvester, and Riccati ADI methods under a single recursive Petrov-Galerkin interpolation framework. In this view, the expensive shifted linear solves are shared, while the difference between Lyapunov, Sylvester, and Riccati variants lies mainly in small-scale pole-placement steps. One concrete claim is that two shifted linear solves per iteration suffice to simultaneously solve six Lyapunov equations, one Sylvester equation, and ten Riccati equations; the same framework also exposes the reduced-order models implicitly constructed by the underlying ADI basis (Zulfiqar et al., 4 Dec 2025).

A broader conceptual extension comes from quadrature and Runge–Kutta formulations. By recasting Lyapunov and Sylvester equations as ODE-driven quadrature problems and imposing preservation of a low-rank residual invariant, one obtains structure-preserving Runge–Kutta schemes whose one-stage instance is exactly equivalent to standard low-rank ADI after the parameter transformation $A P E^T + E P A^T + B B^T = 0,$9. This does not change the algorithmic core of LR-ADI, but it clarifies why ADI is naturally compatible with low-rank residual factorizations (Bertram et al., 2019).

Taken together, these developments suggest that LR-ADI is best understood not as a single recurrence but as a family of low-rank, shift-driven, projection-compatible algorithms whose classical Lyapunov form remains the reference case. Its enduring strengths are sparse shifted solves, exact low-rank residual formulas, and a deep compatibility with rational approximation and reduced-order modeling; its persistent challenges are shift generation, large-scale robustness, and extensions to higher-rank, indefinite, or multiterm settings (Wolf et al., 2013, Zulfiqar et al., 4 Dec 2025).