WIP: Solver for Lyapunov and Sylvester equations

[olof3]

Straight-forward Implementation of Bartels--Stewart's algorithm.

It seems to perform quite well and the code is reasonably concise compared to e.g., the implementation in MatrixEquations. The intention is not to pull this into ControlSystems.jl (at least not in the long run), but to see if anyone has some comments @baggepinnen , @mfalt

Here are some benchmarks, the number of (view) allocations are very high but that doesn't seem to hurt the performance (I tried an alternative version with a lot less allocations, but the execution time did not seem to improve, if anything it got a bit worse). The number of allocations is reduced by more than half in julia 1.5 (1.6?).

The existing implementations of sylvc and lyapc in LinearAlgebra (only continuous time) and MatrixEquations use LAPACK's trsyl, which is a bit hard to compete with.

1. sylvc, Complex{Float64}
┌─────┬──────────────────────────────────┬───────────────────────────────────┬──────────────────────────────────┐
│   N │                         LyapTest │                   MatrixEquations │                    LinearAlgebra │
├─────┼──────────────────────────────────┼───────────────────────────────────┼──────────────────────────────────┤
│   5 │ 35.522 μs (19 allocs, 10.08 KiB) │  68.492 μs (39 allocs, 18.66 KiB) │ 34.180 μs (23 allocs, 10.80 KiB) │
│  50 │ 6.060 ms (27 allocs, 451.48 KiB) │ 12.201 ms (51 allocs, 746.84 KiB) │ 6.756 ms (32 allocs, 490.86 KiB) │
│ 500 │   2.507 s (27 allocs, 31.05 MiB) │    4.240 s (51 allocs, 46.83 MiB) │   2.470 s (32 allocs, 34.86 MiB) │
└─────┴──────────────────────────────────┴───────────────────────────────────┴──────────────────────────────────┘
1. sylvc, Float64
┌─────┬───────────────────────────────────────┬──────────────────────────────────┬───────────────────────────────────┐
│   N │                              LyapTest │                  MatrixEquations │                     LinearAlgebra │
├─────┼───────────────────────────────────────┼──────────────────────────────────┼───────────────────────────────────┤
│   5 │     19.701 μs (123 allocs, 12.98 KiB) │ 27.519 μs (48 allocs, 12.17 KiB) │   15.506 μs (27 allocs, 6.94 KiB) │
│  50 │    2.220 ms (5404 allocs, 568.05 KiB) │ 3.753 ms (60 allocs, 381.42 KiB) │  2.084 ms (36 allocs, 249.64 KiB) │
│ 500 │ 870.797 ms (468566 allocs, 44.17 MiB) │   1.372 s (60 allocs, 23.47 MiB) │ 891.372 ms (36 allocs, 17.46 MiB) │
└─────┴───────────────────────────────────────┴──────────────────────────────────┴───────────────────────────────────┘
2. lyapc, Complex{Float64}
┌─────┬──────────────────────────────────┬──────────────────────────────────┬──────────────────────────────────┐
│   N │                         LyapTest │                  MatrixEquations │                    LinearAlgebra │
├─────┼──────────────────────────────────┼──────────────────────────────────┼──────────────────────────────────┤
│   5 │  18.227 μs (12 allocs, 6.02 KiB) │  17.968 μs (20 allocs, 8.11 KiB) │  20.011 μs (16 allocs, 6.73 KiB) │
│  50 │ 3.222 ms (18 allocs, 304.03 KiB) │ 6.852 ms (29 allocs, 424.97 KiB) │ 3.908 ms (23 allocs, 343.41 KiB) │
│ 500 │   1.281 s (18 allocs, 23.15 MiB) │   1.187 s (29 allocs, 34.63 MiB) │   1.411 s (23 allocs, 26.97 MiB) │
└─────┴──────────────────────────────────┴──────────────────────────────────┴──────────────────────────────────┘
2. lyapc, Float64
┌─────┬───────────────────────────────────────┬────────────────────────────────────────┬───────────────────────────────────┐
│   N │                              LyapTest │                        MatrixEquations │                     LinearAlgebra │
├─────┼───────────────────────────────────────┼────────────────────────────────────────┼───────────────────────────────────┤
│   5 │       12.416 μs (78 allocs, 8.09 KiB) │      33.488 μs (218 allocs, 23.08 KiB) │   12.523 μs (18 allocs, 4.30 KiB) │
│  50 │    1.132 ms (2614 allocs, 318.77 KiB) │       3.480 ms (7903 allocs, 1.23 MiB) │  1.108 ms (25 allocs, 174.05 KiB) │
│ 500 │ 520.334 ms (235194 allocs, 25.97 MiB) │ 601.178 ms (706382 allocs, 263.91 MiB) │ 607.320 ms (25 allocs, 13.50 MiB) │
└─────┴───────────────────────────────────────┴────────────────────────────────────────┴───────────────────────────────────┘
3. sylvd, Complex{Float64}
┌─────┬──────────────────────────────────┬────────────────────────────────────┐
│   N │                         LyapTest │                    MatrixEquations │
├─────┼──────────────────────────────────┼────────────────────────────────────┤
│   5 │ 35.822 μs (20 allocs, 10.56 KiB) │  45.797 μs (204 allocs, 30.56 KiB) │
│  50 │ 6.170 ms (29 allocs, 490.63 KiB) │  9.911 ms (22082 allocs, 8.08 MiB) │
│ 500 │   2.186 s (29 allocs, 34.86 MiB) │ 4.695 s (2245532 allocs, 5.88 GiB) │
└─────┴──────────────────────────────────┴────────────────────────────────────┘
3. sylvd, Float64
┌─────┬───────────────────────────────────────┬────────────────────────────────────┐
│   N │                              LyapTest │                    MatrixEquations │
├─────┼───────────────────────────────────────┼────────────────────────────────────┤
│   5 │     23.834 μs (148 allocs, 14.77 KiB) │  37.485 μs (376 allocs, 39.27 KiB) │
│  50 │    2.352 ms (6243 allocs, 640.00 KiB) │  3.619 ms (19225 allocs, 3.71 MiB) │
│ 500 │ 901.859 ms (536167 allocs, 50.20 MiB) │ 1.652 s (1673379 allocs, 1.65 GiB) │
└─────┴───────────────────────────────────────┴────────────────────────────────────┘
4. lyapd, Complex{Float64}
┌─────┬──────────────────────────────────┬──────────────────────────────────┐
│   N │                         LyapTest │                  MatrixEquations │
├─────┼──────────────────────────────────┼──────────────────────────────────┤
│   5 │  18.391 μs (13 allocs, 6.50 KiB) │  21.959 μs (20 allocs, 8.11 KiB) │
│  50 │ 2.839 ms (20 allocs, 343.17 KiB) │ 6.858 ms (29 allocs, 424.97 KiB) │
│ 500 │   1.185 s (20 allocs, 26.97 MiB) │   1.276 s (29 allocs, 34.63 MiB) │
└─────┴──────────────────────────────────┴──────────────────────────────────┘
4. lyapd, Float64
┌─────┬───────────────────────────────────────┬────────────────────────────────────────┐
│   N │                              LyapTest │                        MatrixEquations │
├─────┼───────────────────────────────────────┼────────────────────────────────────────┤
│   5 │       14.759 μs (97 allocs, 9.50 KiB) │      23.062 μs (253 allocs, 25.86 KiB) │
│  50 │    1.265 ms (3048 allocs, 365.41 KiB) │       4.024 ms (9477 allocs, 1.68 MiB) │
│ 500 │ 463.967 ms (269384 allocs, 29.96 MiB) │ 724.206 ms (841792 allocs, 604.26 MiB) │
└─────┴───────────────────────────────────────┴────────────────────────────────────────┘

[coveralls]

Coverage decreased (-8.4%) to 50.765% when pulling 67af15a on olof3:sylvester_work into b059312 on JuliaControl:master.

[coveralls]

Coverage decreased (-8.8%) to 50.35% when pulling eb4e945 on olof3:sylvester_work into b059312 on JuliaControl:master.

[baggepinnen]

On my working branch, I have replaced usage of lyap with the corresponding function from MatrixEquations.jl, they seem to perform better and error less often. Do you think it would be a good idea to take up a dependency on ME.jl and make use of it more extensively?

[imciner2]

Do you think it would be a good idea to take up a dependency on ME.jl and make use of it more extensively?

Yes, I think we should definitely switch over to using the MatrixEquations package instead of having custom Solvers. One thing that also immediately gets us is the ability to have the cross-term in the LQR design functions.

[olof3]

The current implementation of dlyap use the "naive" version and is really worthless, so anything would be better than that. MatrixEquations uses vanilla Bartels-Stewart, which is quite straight-forward, but the implementation is thousands of lines ported Fortran code (like most of the other code in that package).

Some time ago I made an attempt a shot to at a more appropriate Julia implementation, which made it a lot more readable IMO. https://github.com/olof3/SylvesterEquations.jl. For most testcases it performs on par with ME, sometimes some 30-50% faster, (the numbers in the first post of this thread are a bit outdated).

For the Riccati equations it is less trivial, but the version in MatrixEquations seems quite vanilla. I have made some attempts to get something together for that as well.

On my todo list is to get something together and investigate if a joint effort can be made MatrixEquations (which I guess would be ideal). The alternative is a new package under JuliaControl, perhaps ControlMatrixEquations, there are advantages with that as well. With the current state of MatrixEquations I don't think we should include it as a dependency too hastily. Give me a week and I should be able to get something together.

What about the naming? I really like the approach in MatrixEquations with lyapc, lyapd, sylvc? I think we should start using that. Also for lqg, etc. #307 (comment)_
possibly keeping care/dare since they really are quite standard, although I think the consistency of arec/ared (as in MatrixEquations) is very nice.

[baggepinnen]

What about the naming? I really like the approach in MatrixEquations with lyapc, lyapd, sylvc? I think we should start using that. Also for lqg, etc. #307 (comment)_

For methods accepting LTISystem we no longer need to add the c/d pre/postfix since this info is contained in the type of the system, hence, we can just call the corresponding methods lqr/kalman. For methods accepting matrices, I'm in favor of having the distinction c/d in the end. We can include deprecations of care/dare to help the user find the right function if we want.

[olof3]

For methods accepting LTISystem we no longer need to add the c/d pre/postfix since this info is contained in the type of the system, hence, we can just call the corresponding methods lqr/kalman.

Right, for lqg we should only accept statespace arguments, so there will be only lqg and possibly (if we added) lqg_sd.
But I I think there should be both lqrc/lqrd (in order not having to come up with some C/D matrices), in addition to lqr/kalman (for state-space arguments).

For methods accepting matrices, I'm in favor of having the distinction c/d in the end. We can include deprecations of care/dare to help the user find the right function if we want.

Excellent, I think we should keep the the reference from care/dare indefinitely in that case. I'm not sure if that makes it a deprecation. care/dare really is very established terminology, so although I personally prefer arec/ared I'm really a bit hesitant