dlyapchol

Square-root solver for discrete-time Lyapunov equations

Syntax

R = dlyapchol(A,B)

R = dlyapchol(A,B,E)

R = dlyapchol(___,Name=Value)

[R,res] = dlyapchol(___)

Description

R = dlyapchol(A,B) returns a Cholesky factorization X = R^TR of the solution X to the Lyapunov matrix equation:

$$AX{A}^T-X+B{B}^T=0$$

All eigenvalues of A matrix must lie in the open unit disk for R to exist.

If A is a sparse matrix, the function returns a low rank-approximation X ≈ R^TR. (since R2026a)

R = dlyapchol(A,B,E) computes a Cholesky factorization X = R^TR of X solving the generalized Lyapunov equation

$$AX{A}^T-EX{E}^T+B{B}^T=0$$

All generalized eigenvalues of (A,E) must lie in the open unit disk for R to exist.

If both A and E are sparse matrices, the function returns a low rank-approximation X ≈ R^TR. (since R2026a)

R = dlyapchol(___,Name=Value) specifies additional options. All options except Scaling apply only when A or A and E are sparse matrices. (since R2026a)

[R,res] = dlyapchol(___) also returns the relative residual in the Frobenius norm. (since R2026a)

Examples

Compute Cholesky Factorization of Lyapunov Solution

This example shows how to directly compute the Cholesky factor of the solution of discrete Lyapunov equations.

Consider a randomly generated state-space model and obtain its data.

rng(0)
sysd = drss(8);
[A,B,~,~] = ssdata(sysd);

Obtain the Cholesky factorization $X=R^{T}R$ of the solution of Lyapunov equation $AX{A}^T-X+B{B}^T=0$.

[r,resF] = dlyapchol(A,B)

r = 8×8

    0.4172    1.0934   -0.3481   -0.3331    0.6403   -0.5475   -0.5600   -0.9851
         0    0.5844    0.6294   -0.6418   -1.2987    0.6227    0.2539    1.0850
         0         0    0.8293   -0.2011   -1.5942    1.3404    0.1181    0.6564
         0         0         0    0.2972   -0.1477    0.0388   -0.2957    0.0511
         0         0         0         0    0.2342   -0.2879   -0.1138   -0.1482
         0         0         0         0         0    0.0337   -0.0086   -0.0273
         0         0         0         0         0         0    0.0276    0.0082
         0         0         0         0         0         0         0    0.0004

resF = 
2.6404e-14

Compute Low-Rank Cholesky Factor of Discrete Lyapunov Equation with Sparse Matrices

Since R2026a

This example shows how to directly compute a low-rank Cholesky factor of the solution of discrete-time Lyapunov equations for sparse matrices.

For this example, consider flowmeterSparse.mat which contains a continuous-time sparss model sys. Load the model sys and convert it to discrete time.

load flowmeterSparse.mat
Ts = 0.1;
sysd = c2d(sys,Ts,"tustin")

Sparse discrete-time state-space model with 5 outputs, 1 inputs, and 9669 states.
Model Properties

Use "spy" and "showStateInfo" to inspect model structure. 
Type "help sparssOptions" for available solver options for this model.

Use sparssdata to extract the sparse matrices A, B, and E

[A,B,~,~,E] = sparssdata(sysd);

Obtain the low-rank approximation of the Cholesky factorization $X\approx R^{T}R$ of the solution of Lyapunov equation $AX{A}^T-EX{E}^T+B{B}^T=0$.

[R,res] = dlyapchol(A,B,E);

Initializing...
Running ADI with built-in shifts.....
Running ADI with adaptive shifts.............
Solved Lyapunov equations to desired accuracy.

size(R)

ans = 1×2

          54        9669

res

res = 
1.9702e-12

dlyapchol uses the low-rank alternating directions implicit (LR-ADI) algorithm to compute a low-rank approximation of the Cholesky factorization.

You can also specify additional options, such as increase RankTol to reduce the rank of the factorization.

[R,res] = dlyapchol(A,B,E,RankTol=1e-5);

Initializing...
Running ADI with built-in shifts.....
Running ADI with adaptive shifts.............
Solved Lyapunov equations to desired accuracy.

size(R)

ans = 1×2

          35        9669

res

res = 
1.9702e-12

The algorithm now returns a lower-rank approximation while maintaining a low residual.

Input Arguments

A, B, E — Input matrices

matrices

Input matrices, specified as matrices of the following sizes.

• A and E — m-by-m square matrices

• B — m-by-n matrix

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Example: [R,res] = dlyapchol(A,B,E,RankTol=1e-5);

Focus — Dynamic range

[0 Inf] (default) | vector

Since R2026a

Dynamic range, specified as a vector of form [fmix,fmax]. Use this option to speed up the LR-ADI algorithm when you know that all poles are in this range. When fmin = 0, the smallest natural frequency is estimated by inverse power iterations. When fmax = Inf, the largest natural frequency is estimated by power iterations.

Because dlyapchol does not take the sample time Ts, you must specify the frequencies as normalized values w*Ts, with pi corresponding to the Nyquist frequency (that is, dlyapchol assumes Ts = 1).

CustomShift — Shifts to accelerate convergence

zeros(0,1) (default) | scalar | vector

Since R2026a

Custom shifts to accelerate the convergence of the algorithm, specified as a scalar or an r-by-1 vector.

Specify shifts based on prior knowledge of pole locations. The algorithm applies these shifts in addition to the default shifts. See [4] and [5] for details.

MaxRank — Maximum rank of Cholesky factors

5000 (default) | positive integer

Since R2026a

Maximum rank of Cholesky factors, specified as a positive integer. The LR-ADI algorithm terminates when the row size of R exceeds this target. The final row size of R can be larger.

LyapTol — Relative tolerance for Lyapunov residuals

1e-8 (default) | positive scalar

Since R2026a

Maximum value for the relative residual, specified as a positive scalar. The LR-ADI iterations stop when the relative residual falls below LyapTol. Increasing LyapTol helps speed up computation at the expense of accuracy.

RankTol — Relative tolerance for rank decisions

1e-8 (default) | positive scalar

Since R2026a

Relative tolerance for rank decisions, specified as a positive scalar. Increasing this value reduces the rank of R but can also increase the residual by dropping parts of R.

Scaling — Automatic scaling of A and E

"on" (default) | "off"

Since R2026a

Use this option to enable or disable scaling of A and E matrices. This option does a form of balancing for full matrices and uses equilibrate for sparse matrices. Scaling can improve accuracy by compressing the numerical range but can sometimes make things worse when better scaling for (A,E) results in worse scaling for B.

Display — Show or hide progress report

"on" (default) | "off"

Since R2026a

Show or hide progress report for LR-ADI iterations, specified as either "off" or "on".

Output Arguments

R — Cholesky factor of Lyapunov solution

matrix

Cholesky factor of the Lyapunov solution, returned as matrix of size:

• m-by-m square matrix when A, E, or both are full matrices.

• r-by-m matrix when both A and E are sparse matrices. Typically, r ≪ m. (since R2026a)

res — Relative residual

scalar

Since R2026a

Relative residual, returned as a scalar. The relative residual is defined as follows:

$$\frac{{\Vert AX{A}^T-EX{E}^T+B{B}^T\Vert }_F}{{\Vert B{B}^T\Vert }_F}.$$

Algorithms

For full matrices, dlyapchol uses SLICOT routines SB03OD and SG03BD.

For sparse matrices, dlyapchol uses the low-rank alternating directions implicit (LR-ADI) algorithm to compute a low-rank factorization X ≈ R^TR where R is wide and skinny. For singular E, dlyapchol solves the projected Lyapunov equations described in [6]. (since R2026a)

References

[1] Bartels, R.H. and G.W. Stewart, "Solution of the Matrix Equation AX + XB = C," Comm. of the ACM, Vol. 15, No. 9, 1972.

[2] Hammarling, S.J., “Numerical solution of the stable, non-negative definite Lyapunov equation,” IMA J. Num. Anal., Vol. 2, pp. 303-325, 1982.

[3] Penzl, T., ”Numerical solution of generalized Lyapunov equations,” Advances in Comp. Math., Vol. 8, pp. 33-48, 1998.

[4] Benner, Peter, Jing-Rebecca Li, and Thilo Penzl. “Numerical Solution of Large-Scale Lyapunov Equations, Riccati Equations, and Linear-Quadratic Optimal Control Problems.” Numerical Linear Algebra with Applications 15, no. 9 (November 2008): 755–77. https://doi.org/10.1002/nla.622.

[5] Benner, Peter, Martin Köhler, and Jens Saak. “Matrix Equations, Sparse Solvers: M-M.E.S.S.-2.0.1—Philosophy, Features, and Application for (Parametric) Model Order Reduction.” In Model Reduction of Complex Dynamical Systems, edited by Peter Benner, Tobias Breiten, Heike Faßbender, Michael Hinze, Tatjana Stykel, and Ralf Zimmermann, 171:369–92. Cham: Springer International Publishing, 2021. https://doi.org/10.1007/978-3-030-72983-7_18.

[6] Benner, Peter, and Tatjana Stykel. “Model Order Reduction for Differential-Algebraic Equations: A Survey.” In Surveys in Differential-Algebraic Equations IV, edited by Achim Ilchmann and Timo Reis. Springer International Publishing, 2017. https://doi.org/10.1007/978-3-319-46618-7_3.

Version History

Introduced before R2006a

R2026a: Support for LR-ADI algorithm for solving equations with sparse matrices

dlyapchol now uses the low-rank alternating directions implicit (LR-ADI) algorithm to compute a low-rank Cholesky factorization when both A and E inputs are sparse matrices. Previously, the function converted A and E to full matrices before computing the factorization. To get the old behavior, you can call dlyapchol with full(A) instead of A as input.