Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: Solver for Lyapunov and Sylvester equations #303

Open
wants to merge 16 commits into
base: master
Choose a base branch
from

Conversation

olof3
Copy link
Contributor

@olof3 olof3 commented May 13, 2020

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 │
├─────┼──────────────────────────────────┼───────────────────────────────────┼──────────────────────────────────┤
│   535.522 μs (19 allocs, 10.08 KiB) │  68.492 μs (39 allocs, 18.66 KiB) │ 34.180 μs (23 allocs, 10.80 KiB) │
│  506.060 ms (27 allocs, 451.48 KiB) │ 12.201 ms (51 allocs, 746.84 KiB) │ 6.756 ms (32 allocs, 490.86 KiB) │
│ 5002.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 │
├─────┼───────────────────────────────────────┼──────────────────────────────────┼───────────────────────────────────┤
│   519.701 μs (123 allocs, 12.98 KiB) │ 27.519 μs (48 allocs, 12.17 KiB) │   15.506 μs (27 allocs, 6.94 KiB) │
│  502.220 ms (5404 allocs, 568.05 KiB) │ 3.753 ms (60 allocs, 381.42 KiB) │  2.084 ms (36 allocs, 249.64 KiB) │
│ 500870.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 │
├─────┼──────────────────────────────────┼──────────────────────────────────┼──────────────────────────────────┤
│   518.227 μs (12 allocs, 6.02 KiB) │  17.968 μs (20 allocs, 8.11 KiB) │  20.011 μs (16 allocs, 6.73 KiB) │
│  503.222 ms (18 allocs, 304.03 KiB) │ 6.852 ms (29 allocs, 424.97 KiB) │ 3.908 ms (23 allocs, 343.41 KiB) │
│ 5001.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 │
├─────┼───────────────────────────────────────┼────────────────────────────────────────┼───────────────────────────────────┤
│   512.416 μs (78 allocs, 8.09 KiB) │      33.488 μs (218 allocs, 23.08 KiB) │   12.523 μs (18 allocs, 4.30 KiB) │
│  501.132 ms (2614 allocs, 318.77 KiB) │       3.480 ms (7903 allocs, 1.23 MiB) │  1.108 ms (25 allocs, 174.05 KiB) │
│ 500520.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 │
├─────┼──────────────────────────────────┼────────────────────────────────────┤
│   535.822 μs (20 allocs, 10.56 KiB) │  45.797 μs (204 allocs, 30.56 KiB) │
│  506.170 ms (29 allocs, 490.63 KiB) │  9.911 ms (22082 allocs, 8.08 MiB) │
│ 5002.186 s (29 allocs, 34.86 MiB) │ 4.695 s (2245532 allocs, 5.88 GiB) │
└─────┴──────────────────────────────────┴────────────────────────────────────┘
3. sylvd, Float64
┌─────┬───────────────────────────────────────┬────────────────────────────────────┐
│   N │                              LyapTest │                    MatrixEquations │
├─────┼───────────────────────────────────────┼────────────────────────────────────┤
│   523.834 μs (148 allocs, 14.77 KiB) │  37.485 μs (376 allocs, 39.27 KiB) │
│  502.352 ms (6243 allocs, 640.00 KiB) │  3.619 ms (19225 allocs, 3.71 MiB) │
│ 500901.859 ms (536167 allocs, 50.20 MiB) │ 1.652 s (1673379 allocs, 1.65 GiB) │
└─────┴───────────────────────────────────────┴────────────────────────────────────┘
4. lyapd, Complex{Float64}
┌─────┬──────────────────────────────────┬──────────────────────────────────┐
│   N │                         LyapTest │                  MatrixEquations │
├─────┼──────────────────────────────────┼──────────────────────────────────┤
│   518.391 μs (13 allocs, 6.50 KiB) │  21.959 μs (20 allocs, 8.11 KiB) │
│  502.839 ms (20 allocs, 343.17 KiB) │ 6.858 ms (29 allocs, 424.97 KiB) │
│ 5001.185 s (20 allocs, 26.97 MiB) │   1.276 s (29 allocs, 34.63 MiB) │
└─────┴──────────────────────────────────┴──────────────────────────────────┘
4. lyapd, Float64
┌─────┬───────────────────────────────────────┬────────────────────────────────────────┐
│   N │                              LyapTest │                        MatrixEquations │
├─────┼───────────────────────────────────────┼────────────────────────────────────────┤
│   514.759 μs (97 allocs, 9.50 KiB) │      23.062 μs (253 allocs, 25.86 KiB) │
│  501.265 ms (3048 allocs, 365.41 KiB) │       4.024 ms (9477 allocs, 1.68 MiB) │
│ 500463.967 ms (269384 allocs, 29.96 MiB) │ 724.206 ms (841792 allocs, 604.26 MiB) │
└─────┴───────────────────────────────────────┴────────────────────────────────────────┘

@coveralls
Copy link

Coverage Status

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

@coveralls
Copy link

coveralls commented May 13, 2020

Coverage Status

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

src/LyapTest.jl Outdated Show resolved Hide resolved
src/LyapTest.jl Outdated Show resolved Hide resolved
src/LyapTest.jl Outdated Show resolved Hide resolved
@baggepinnen
Copy link
Member

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
Copy link
Contributor

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
Copy link
Contributor Author

olof3 commented Jan 17, 2021

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
Copy link
Member

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
Copy link
Contributor Author

olof3 commented Jan 17, 2021

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants