Scipy riccati solvers are too naive

[ilayn]

It basically implements the textbook definition. Implement the Newton + line search for both care and dare (with proper function names if possible). Slyvester and Lyapunov solvers use LAPACK versions hence they are OK on paper.

Reminder : For control purposes, speed is irrelevant! So being implemented in Fortran doesn't mean much. The residuals and conditioning of the variables decide whether certain synthesis algorithms work or not. This needs to be state-of-art without relying on anything off-the-shelf.

[ilayn]

This has led to a turtle-all-the-way-down situation. Hence, I'm trying to replicate the Sylvester solver which is needed for Lyapunov solvers which is needed for Riccati solvers. Current implementation is related to SLICOT 2014 report by Köhler and Saak and also Kågström and Westin and then Penzl. Need to figure out a way to make incremental advances. Full blown implementation is too much work and might resort back to LAPACK and implement an error bound warning.

[ilayn]

Lyapunov solvers are in place but slow (~10x) but more precise. For problems up to 300x300 the difference is bearable. There are obvious improvements waiting but main bottleneck is np.linalg.solve doing too many unnecessary checks at each call. Might switch to the low-level gufunc instead.

[ilayn]

Commits 96814038b2f68ba379e133e24f2a27cfddeab9be and a39813b4d80f70a9341a00ff3c4d91a520aacf1d for the regular Lyapunov equations brought down the speed difference to x3. The precision for not-ill-conditioned problems are at least equivalent. For discrete Lyapunov equations the precision gets better with problem size (with respect to max residual) if compared to scipy and matlab.

[ilayn]

Lyapunov solvers are working pretty OK. For Riccati solvers I'm stuck with the balancing algorithm of LAPACK's GEBAL not performing as good as matlab's balance.m. I surely appreciate any pointers if anybody wants to chime in.

[ilayn]

With, commit ba4e9731d1214e0e394d7293af45128e3fda0e16 scipy solvers are dropped and custom ones are in place. There is still this difference between GEBAL and balance.m which is beyond me.

[ilayn]

Changed my mind an started to commit riccati solvers to official scipy repo directly.

[ilayn]

SciPy 0.19 will have the new Riccati solvers via https://github.com/scipy/scipy/pull/6616. Lyapunov solvers are on the way either to that or to the next release of SciPy. It actually makes more sense to commit to numpy and scipy as much as possible since it is maintained and reviewed more meticulously.