Refactor matrix equations (Ricatti, Lyapunov, Sylvester)

[albop]

Currently these equations live in different modules and some of them are not easily discoverable.

For all of them there are also several possible implementations:

• in scipy (calling linpack)
• in quantecon (i.e. doublej)
• in Slycot, which is based on the last opensource version of SLICOT

Slycot, used to be complicated to build, but this has changed, and it is now pip-installable. It is supposed to be efficient and AFAIK is the only option for Generalized Sylvester equations.
The python-control library has a mateqn.py file (https://github.com/python-control/python-control/blob/master/src/mateqn.py) that is also quite complete and that we could partly reuse or depend on. It seems to do extensive error handling and input checks.

My point is that we should deal with all of these matrix equations for optimal control together in a systematic way. For instance we could

• move them to a common module (like matrix_equations.py)
• provide wrapper functions (solve_lyapunov, solve_generalized_sylvester, solve_riccati, etc...) that can switch between different implementations
• investigate which implementation is the best.

[jstac]

@albop Thanks for these suggestions. It makes sense for matrix equations to be treated systematically and located in the same module. @cc7768 is currently working on this.

[sglyon]

2764a8c added a matrix_eqn.py file with lyapunov and riccati solvers.

Is this satisfactory enough to close the issue, @albop? We didn't do anything with slycot or any investigation for which algorithm works best.

[cc7768]

Still missing Sylvester equations. Needs to stay open.
On Aug 28, 2014 8:00 PM, "Spencer Lyon" notifications@github.com wrote:

2764a8c
2764a8c
added a matrix_eqn.py file with lyapunov and riccati solvers.

Is this satisfactory enough to close the issue, @albop
https://github.com/albop? We didn't do anything with slycot or any
investigation for which algorithm works best.

—
Reply to this email directly or view it on GitHub
#47 (comment).

[mmcky]

@albop @cc7768 @spencerlyon2 Are we still looking to implement a solver for Sylvester equations? Is this just requiring a wrapper on: http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.linalg.solve_sylvester.html

[albop]

There are several types of Sylvester equations. the one in scipy is
AX+XB=Q, but there is also the so-called generalized sylvester form AX + B.X.C=D which is useful for higher order perturbations, and is implemented
in slycot.

Actually, in these higher order perturbation problems, the system has a
more precise structure AX + B.X.[C,...C]=D where X is tensor
(multidimensional matrix). There is a better way to do that than slycot,
documented in Kamenik 2005.
In https://github.com/EconForge/dolo/blob/master/dolo/numeric/matrix_equations.py
wrote a function which solve the third problem, by reformulating it as a
generalized sylvester and using slycot on it, so it' s not optimal but very
useful (well, for my workflow, at least).

If you want we can move this function to quantecon, and think later about
improving it. The cost is an additional dependency on slycot, which is
another non-trivial piece of fortran software.

The one I use, for higher order discrete time approximation is sometimes
called generalized sylvester equation in the form AX + B.X.C=0 while
Actually the useful

On Tue, Aug 4, 2015 at 4:12 PM, mmcky notifications@github.com wrote:

@albop https://github.com/albop @cc7768 https://github.com/cc7768
@spencerlyon2 https://github.com/spencerlyon2 Are we still looking to
implement a solver for Sylvester equations? Is this just requiring a
wrapper on:
http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.linalg.solve_sylvester.html

—
Reply to this email directly or view it on GitHub
#47 (comment)
.

[cc7768]

I think this would still be nice to have. I don't know that it is an immediate priority -- Hooking into slycot might be problematic though. What do others think about having dependencies that some might be unable (due to OS or other reasons) to compile? Hooking into the scipy code also seems like a reasonable solution.

@albop or @jstac : Do you know if there is a simple algorithm (simple in the sense of the doubling algorithm for solving the Ricatti equation) that we could write ourselves?

[jstac]

I would like to avoid adding a dependency such as slycot unless it's really critical.

Here's a discussion of algorithms for generalized Sylvester equations with a link to a paper:

http://math.stackexchange.com/questions/568405/efficient-algorithm-for-generalized-sylvesters-equation