Solving Lyapunov equations

[jstac]

This is something we do a lot of. At present we use SciPy's solve_discrete_lyapunov. It might be worth comparing this against Tom's doublej routine.

[cc7768]

Hey @jstac.

So I spent a couple minutes looking at this (very contrived examples, nothing really interesting), and it looks like doublej(same as quantecon.quadsums.m_quadratic_sum) outperforms the scipy.linalg.solve_discrete_lyapunov for the cases where it applies (If we crank up the dimensions then the scipy version runs into memory errors because its algorithm uses a kronecker product while the doublej seems to work just fine). It only applies when both matrices have the modulus of their eigenvalues bounded by unit, will we be using it for other cases?

[jstac]

@cc7768: Thanks for this. Regarding the eigenvalues in different use cases, I'm not 100% sure. Would you mind to run that by Tom directly via email?

My preference would be to use doublej and perhaps put a note in the docstring indicating the existence of the SciPy option and a short rationale for using doublej (i.e., better performance).

Finally, I suggest we rename doublej to solve_discrete_lyapunov.

[cc7768]

@jstac It is currently named m_quadratic_sum in our files. I was thinking of making a utils.py file with some of these, I could just create a wrapper function called solve_discrete_lyapunov that simply calls m_quadratic_sum?

[albop]

If I understand well, the doubling algorithm implemented in doublej and m_quadratic_sum can be applied to solve Sylvester equations and also, as a special case, Lyapunov equations. Like @jstac said, it is indeed a strategic part and I would suggest a separate module for it (I mean a .py file) with at least both methods (Sylvester and Lyapunov). There are other linear methods for optimal controls that could naturally fit in it too (I will open an issue for that).
My bet would have been that the doubling algorithm works better for big systems, not necessarily for smaller ones. At any rate, one cannot discard, the case where some relevant eigenvalues are not strictly smaller than one. It contains the unit root case that is very useful for many models at least in the Sylvester case.
Extending cc7768's proposition, I would vote for two Sylvester/Lyapunov functions, each one taking amethod keyword argument whose value could be bartels-stewart for the Scipy implementation or doubling for Tom's implementation. In the latter case, the exception from m_quadratic could be caught, and transformed into a more explicit one asking to try a different method.
Looking at Dynare's manual, I see that they have kept 4 different implementations, which is a strong indication that one version will not cover all cases. They use the Bartels-Stewart as the default, so it's probably a safe bet too.

[cc7768]

@albop @jstac Thanks for the input. I think it is a good idea to put the Ricatti/Lyapunov/Sylvester solvers into the same file. I will write up the Lyapunov function and then submit a PR to close this issue and we can address Sylvester equations and other generalities in the #47.