lyapchol

Square-root solver for continuous-time Lyapunov equation

Syntax

R = lyapchol(A,B)

R = lyapchol(A,B,E)

R = lyapchol(___,Name=Value)

[R,res] = lyapchol(___)

Description

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

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

All eigenvalues of matrix A must lie in the open left half-plane for R to exist.

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

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

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

All generalized eigenvalues of (A,E) must lie in the open left half-plane 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 = lyapchol(___,Name=Value) specifies additional options. All options except Scaling apply only when A or A and E are sparse matrices. (since R2026a)

[R,res] = lyapchol(___) 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 Lyapunov equations.

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

rng(0)
sys = rss(8);
[A,B,~,~] = ssdata(sys);

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

[R,res] = lyapchol(A,B)

R = 8×8

    0.4089    0.9042   -0.3494   -0.4282    1.1370   -0.8339   -0.6044   -1.4592
         0    0.2074    0.2725   -0.2124   -0.5418    0.2698    0.0737    0.4182
         0         0    0.5366   -0.3523   -0.6570    0.5135    0.1761    0.2253
         0         0         0    0.1297   -0.1899    0.1940   -0.0609    0.0885
         0         0         0         0    0.0354   -0.0115   -0.0248   -0.0433
         0         0         0         0         0    0.0223    0.0070   -0.0115
         0         0         0         0         0         0    0.0058    0.0017
         0         0         0         0         0         0         0    0.0001

res = 
7.1792e-14

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

Since R2026a

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

For this example, consider flowmeterSparse.mat which contains a continuous-time sparss model sys. Load the model sys to the workspace and use sparssdata to extract the sparse matrices A, B, and E.

load flowmeterSparse.mat   
[A,B,~,~,E] = sparssdata(sys);

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

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

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

size(R)

ans = 1×2

          36        9669

res

res = 
9.0288e-09

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

Specify Frequency Focus and Custom Shifts in LR-ADI for Lyapunov Equation

Since R2026a

This example shows how specifying Focus and CustomShift options affects LR-ADI convergence when computing the low-rank Cholesky factor for sparse Lyapunov equations.

Load the sparse model and extract the sparse matrices.

load flowmeterSparse.mat   
[A,B,~,~,E] = sparssdata(sys);

Compute the Cholesky factor. By default, if you do not specify any name-value options, the algorithm uses their default values to compute the low-rank factor.

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

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

size(R)

ans = 1×2

          36        9669

res

res = 
9.0288e-09

You can use additional options, such as Focus, to specify a dynamic range and speed up LR-ADI when all poles are in this range. For example, the dynamic range of this system is approximately [0 4.64e5], as the largest eigenvalue is about -4.6e5.

eigs(A,E,1)

ans = 
-4.6449e+05

If you set the Focus range to cover this entire interval, you get roughly the same solution as before.

[R2,res2]=lyapchol(A,B,E,Focus=[0 4.65e5]);

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

size(R2)

ans = 1×2

          36        9669

res2

res2 = 
9.3884e-09

The LR-ADI algorithm selects shifts within the specified frequency range. If the focus range is too narrow, the algorithm may fail to converge. For example, specify a frequency range up to 1000 rad/s.

[R3,res3]=lyapchol(A,B,E,Focus=[0 1e3]);

Initializing...
Running ADI with built-in shifts.....
Running ADI with adaptive shifts........................................
........................................................................
.......
Could not solve Lyapunov equations to desired accuracy. 
This is typically due to A or (A,E) having undamped or poorly damped eigenvalues.

The algorithm fails to converge because this range excludes some key poles. To restore the convergence, you can use the CustomShift option to introduce a shift near -4.6e5.

[R4,res4]=lyapchol(A,B,E,Focus=[0 1e3],CustomShift=-5e5);

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

size(R4)

ans = 1×2

         116        9669

res4

res4 = 
9.6902e-09

The algorithm now converges to desired accuracy.

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]=lyapchol(A,B,E,Focus=[0 1e3],CustomShift=-4.64e5)

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.

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{E}^T+EX{A}^T+B{B}^T\Vert }_F}{{\Vert B{B}^T\Vert }_F}.$$

Algorithms

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

For sparse matrices, lyapchol 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, lyapchol 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

lyapchol now uses the low-rank alternating directions implicit (LR-ADI) algorithm to compute a low-rank Cholesky factorization when both A and E 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 lyapchol with full(A) instead of A as input.