11
$\begingroup$

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.

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

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.

Improve this answer
$\endgroup$
1
3
$\begingroup$

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.

Improve this answer
$\endgroup$
3
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.