MatrixEquations.jl

MatrixEquations.jl copied to clipboard

We have a quick implementation of ADI:

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

This could live in this package.

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

Whoops! Its public now

Could you provide an example for using this function?

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

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.

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.

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.)

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

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!

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