search

andreasvarga / MatrixEquations.jl

Solution of Lyapunov, Sylvester and Riccati matrix equations using Julia
MIT License
79 stars 8 forks source link

Move Alternating-direction Implicit (ADI) Method solver here? #25

Open dlfivefifty opened 11 months ago

dlfivefifty commented 11 months ago

We have a quick implementation of ADI:

https://github.com/dlfivefifty/AlternatingDirectionImplicit.jl

This could live in this package.

baggepinnen commented 11 months ago

The repo link you posted appears to be private, I can't access it

dlfivefifty commented 11 months ago

Whoops! Its public now

andreasvarga commented 11 months ago

Could you provide an example for using this function?

baggepinnen commented 11 months ago

The ADI package comes with a number of extra dependencies, something that is worth considering before integrating the functionality here

dlfivefifty commented 11 months ago

I don't have a MWE and at the moment its specialised to my use case but maybe the Wikipedia article is a better description of what the algorithm does:

https://en.wikipedia.org/wiki/Alternating-direction_implicit_method

the numbers a,b,c,d are used to embed the spectrum of the matrices which allow one to deduce apriori the number of iterations needed to converge within the specified tolerance.

andreasvarga commented 11 months ago

It is a good idea to add some solvers for large Sylvester and Lyapunov equations. The provided implementation has in my opinion some issues, as for example, (1) the lack of a meaningful default setting for the parameters a,b,c,d; (2) it is not clear which equation is actually solved, becase the role of C matrix is not explained. At this moment, it appears that only C = I or C a square nxn matrix are acceptable, which restricts A and B to be of the same dimension, which is clearly not desirable. Finally, the solver should be compatible with the already existing solvers, thus should solve AX+XB = C (and not AX-XB = C) (this is however a minor issue at this stage). In this context, an interesting possibility (especially for Lyapunov solvers) is to allow C to have a factored form (e.g., Cholesky) , which is then maintained also by the solution.

I will invest some time to study existing ADI algorithms, especially those which already have good implementations in MATLAB.

I would like very much to know which algorithm has been used to generate the shift parameters in the provided code. The Wikipedia article provides no help in this direction.

dlfivefifty commented 11 months ago

The implementation was just a quick-and-dirty thing, I'm not proposing it be moved as-is.

The parameters can be chosen by just calling eigmin and eigmax which could be to added as the default but In some applications they can be deduced analytically. In my case these are banded matrices so the complexity is O(n) I think, but even O(n^2) is fine as thats optimal. (If A and B are dense its O(n^3) which is probably the same complexity as your existing solver.)

dlfivefifty commented 11 months ago

(Sorry, I forgot to mention we are solving a more general equation AXC - C'XB = F)

andreasvarga commented 11 months ago

The ADI package comes with a number of extra dependencies, something that is worth considering before integrating the functionality here

Do you mean the shift generation stuff or the "factorization" related suff? So, a separate dedicated package for large dimension problems would be also a possibility!

baggepinnen commented 11 months ago

I'm talking about these two dependencies being listed in the project file https://github.com/dlfivefifty/AlternatingDirectionImplicit.jl/blob/main/Project.toml#L7-L8 These are not (direct) dependencies of MatrixEquations, which is a very lightweight package in terms of dependencies