Add Lyapunov / Sylvester solver

[dlfivefifty]

There is no Lyapunov / Sylvester solver built-in to Julia, for solving

AXB + CXD = F

This could be fixed by wrapping the LAPack xTRSYL routine.

[lindahua]

+1

[stevengj]

Seems worthwhile and low-cost. What would the interface look like?

[jiahao]

I'm in favor of just putting in the LAPACK wrapper for now and deferring a higher-level interface to a symbolic rewrite rule in the same vein as #48. Or we could just do the simple thing and define a new solvesylvester function.

[dlfivefifty]

You could call the function lyap to match matlab

Sent from my iPad

On 15 Feb 2014, at 3:47 am, Jiahao Chen [email protected] wrote:

I'm in favor of just putting in the LAPACK wrapper for now and deferring a higher-level interface to a symbolic rewrite rule in the same vein as #48. Or we could just do the simple thing and define a new solvesylvester function.

—
Reply to this email directly or view it on GitHub.

[lsorber]

@dlfivefifty lyap also solves Sylvester equations. Imho, maybe lyap tries to do too much and its functionality should be split into smaller pieces in Julia.

[atisharma]

Would a wrapper to SLICOT be a suitable response?

[mschauer]

As I understood, SLICOT is not really free, http://slicot.org/obtaining-slicot

[atisharma]

SLICOT 4.5 is GPL v3.0. SLICOT 5.0 does indeed look non-free from their website. But strangely it is available in the Debian main repo, and Debian are pretty strict about licencing.

Writing good Lyapunov / Riccati solvers is not that easy and if the SLICOT licence is acceptable it would be preferable to a reimplementation. Perhaps using SLICOT 4.5?

https://github.com/jcrist/Slicot.jl is attempting to create a Julia SLICOT wrapper.

[mschauer]

Ah, tricky. Could be that they had 5.0 for a limited time on the website under GPL and then removed the link, so that version 5.0 would be free as well, http://www.dynare.org/pipermail/dev/2011-July/001230.html .

[johnmyleswhite]

I think we'd generally try to avoid pulling something GPL3 into Base. Seems better suited for a package.

[dlfivefifty]

I've implemented my own lyap in Julia for the time being. If anyone is interested you can find it as part of ApproxFun:

https://github.com/dlfivefifty/ApproxFun/blob/master/src/PDE/lyap.jl

[mschauer]

@dlfivefifty Cool, thank you for implementing it. I would like to use your implementation for the time. Could you put your lyap.jl public domain or under a licence?

[dlfivefifty]

As far as I'm concerned it's public domain, but if you require an official license I'm happy to give it to you

Sent from my iPhone

On Jun 26, 2014, at 10:54 PM, "M. Schauer" [email protected] wrote:

@dlfivefifty Cool, than you for implementing it. I would like to use your implementation for the time. Could you put your lyap.jl public domain or under a licence?

—
Reply to this email directly or view it on GitHub.

[mschauer]

Ok, cool. I also have a look at the xTRSYL

[mschauer]

I gave it a try, if somebody familiar with Julia-LAPACK interaction could have a look, that would be nice, https://gist.github.com/mschauer/e967a2f1e7dc745c23e7 including naive implementation for comparison.
With some guidance I can make a pull request, or somebody else takes over.

[ViralBShah]

You should definitely make a PR, add tests and docs.

[dlfivefifty]

Looking forward to changing my code to use this!  I didn’t see a wrapper for solving

    AXB^t + CXD^t = F

Is a routine for solving this generalized Sylvester equation included in LAPack? Is it easy to access via your code?

Cheers,

Sheehan

On 27 Jun 2014, at 8:03 pm, M. Schauer [email protected] wrote:

I gave it a try, if somebody familiar with Julia-LAPACK interaction could have a look, that would be nice, https://gist.github.com/mschauer/e967a2f1e7dc745c23e7 including naive implementation for comparison.
With some guidance I can make a pull request, or somebody else takes over.

—
Reply to this email directly or view it on GitHub.

[mschauer]

I think there is only Generalized coupled Sylvester in LAPACK, but one can add a wrapper for that as well

DTGSYL solves the generalized Sylvester equation
              A * R - L * B = scale * C                 
              D * R - L * E = scale * F
  where R and L are unknown m-by-n matrices

[dlfivefifty]

Interesting, Mathematica has a LyapunovSolve that does AXB+CXD=F. I'd be surprised if Wolfram bothered to implement it themselves, so I wonder what they use

Sent from my iPad

On 28 Jun 2014, at 10:59 pm, "M. Schauer" [email protected] wrote:

I think there is only Generalized coupled Sylvester in LAPACK, but one can add a wrapper for that as well

DTGSYL solves the generalized Sylvester equation
A * R - L * B = scale * C
D * R - L * E = scale * F
where R and L are unknown m-by-n matrices
—
Reply to this email directly or view it on GitHub.

[andreasnoack]

MATLAB uses SLICOT but we cannot do that because of the license.

[tkelman]

There is a slicot wrapper under-development at https://github.com/jcrist/slicot.jl, it may not have all of the functions exposed yet though.

[JeffBezanson]

implemented by JuliaLang/julia#7435

[atisharma]

It seems the SLICOT library is dual licence. It will be forked if the licence of future versions is closed. The current version is in the Debian repo, and it seems that there is a commitment to maintain the fork. So perhaps it would be possible to implement as a package?
See https://lists.debian.org/debian-science/2012/06/msg00154.html for the copyright status of the library, and https://lists.debian.org/debian-science/2012/06/msg00154.html for the debian mailing list thread.

[mschauer]

Should we reopen? We also still need the discrete time Lyapunov equation

A X A' - X + Q = 0 

[tkelman]

@mschauer I could absolutely use that (also algebraic Riccati if someone wants to be awesome), but if the functionality isn't easily doable via lapack wrappers that would already be in openblas then probably best as a package.

[mschauer]

I fear you are right. Maybe we can ping each other if something starts moving.

[mschauer]

Perhaps not the canonical way, but https://github.com/mschauer/IDRsSolver.jl solves Stein/discrete Lyapunov equation and other linear equations.

[zouhairm]

+1 for discrete time Lyapunov equation and algebraic Riccati...
FWIW, the QuantEcon package https://github.com/QuantEcon/QuantEcon.jl implements solve_discrete_lyapunov.