# Libraries for solving Lyapunov's equation

The following matrix equation $$B\Sigma + \Sigma B^T + C = 0$$ in $\Sigma$ $-$ for given $B$ and $C$ matrices $-$ appears in my work as a characterization of a covariance matrix. I have learned that this equation is known, in particularly in continuous time control theory, as Lyapunov's equation, and that there are various well known algorithms for solving it that exploit the special nature of this linear equation.

From googling I have also learned that there exist Matlab and Fortran implementations. I have found SLICOT and RECSY. Due to licensing issues access to SLICOT source has been stopped, though.

Most of my work is implemented in R, and as I have been unable to find an R interface to a solver, I consider writing one myself. My question is then if SLICOT is the best available Fortran (or C) library with an implementation of a solver of Lyapunov's equation? I am also interested in implementations that can handle large sparse $B$ matrices.

• How big and how sparse? You might have to move away from R to deal with the bigger problems in reasonable time. Jun 26, 2012 at 16:45
• I probably shouldn't say this, but SLICOT is available here. Jun 26, 2012 at 18:11
• @BillBarth, dimensions in the order of 1000, $C$ diagonal and $B$ unstructured but potentially very sparse, 1% non-zero entries, say.
– NRH
Jun 27, 2012 at 6:24

SLICOT is the tool to use for dense problems.

For large but sparse system, there is the lyapack toolbox for MATLAB.

The algorithms in lyapack base on computing iteratively low-rank factors $Z_n$, so that $Z_n^HZ_n$ approaches $\Sigma$, where $\Sigma$ is the symmetric (positive or negative) definite solution of the Lyapunov equation. Computing only the factors in combination with sparsity of the coefficients makes this approach feasible for large scale equations.

There is vivid ongoing research at the Max-Planck Institute in Magdeburg, Germany, on sparse Lyapunov equations. However, the announcement of the upcoming release of the sucessor of lyapack - MESS - is quite a few years old. Nevertheless, it is worth checking the MESS's webpage and the publications of the contributing authors from time to time.

Disclaimer: My thesis supervisor is a major contributor both to SLICOT and lyapack and I am in regular contact with the developers of MESS.

You could connect to MATLAB using this.

Your matrices aren't too large : hand coding the algorithms shouldn't result in too much time loss, maybe it'll run for 1 hour. It may or may not be too long depending on various factors.

Though, coding it yourself may not be easy at all. I don't think I can, and I've been dealing with this for the past few months. But the SLICOT algorithm itself is here.

SLICOT's algorithm is not that complicated, it's a reduction to Schur form + some back-substitution. You can check the Bartels-Stewart paper http://dl.acm.org/citation.cfm?id=361582 which is reasonably readable and explains how it works. The paper is about the nonsymmetric case, but it shouldn't be hard to adapt it to the symmetric one --- you just need one Schur form instead of two.

You can probably also code it yourself in R if it already has a routine for the Schur form (I'd check myself, but it's always a mess to get meaningful results about R out of Google since due to their unfortunate naming choice).

This could settle the dense case. The large and sparse one is more technical.

There's a remarkably simple way to express solution in terms of eigenvalues/eigenvectors. If found this in Matlab implementation, even though none of the standard textbooks talk about this. Here's a proof that it works.

For $$AX+XA'=B$$ with Hermitian $$A$$ associated with a matrix of eigenvectors $$U$$ and column vector of eigenvalues $$s$$, you can write solution as follows:

$$X=U \left( \frac{U' BU}{s + s'} \right) U'$$

Here, $$s+s'$$ is a square matrix of pairwise sums of eigenvalues (this is what happens when you do column+row in numpy) and $$a/b$$ is Hadamard (pointwise) division.

Advantage of this formulation is that it's clear where this algorithm becomes numerically unstable -- the division by small values of $$s$$. To make it stable, simply skip division by near-0 entries. This essentially mirrors the logic of pseudo-inverse which ignores dimensions with near-0 singular values. This seems both more robust and faster than scipy's built-in scipy.linalg.solve_lyapunov, without impacting solution quality.

Here's one implementation of this skipping logic in numpy-like PyTorch syntax.

• Actually, if $A$ is not symmetric/Hermitian, you will need $U^{-1}$ rather than $U^T$, and this is another cause of major stability issues, since nothing guarantees that $U$ is well-conditioned. Jun 30 at 19:14
• Also, small nitpick: the algorithm is backward stable without "skipping logic" even in presence of small denominators: it's the problem that is ill-conditioned. Jun 30 at 21:45