Implementation roadmap

[jiahao]

A comprehensive listing of methods in the references.
Please remove any methods are impractical/obsolete and add methods worth implementing to the list.

Linear systems

• Stationary methods
  • Jacobi (LT)
  • Gauss-Seidel (LT)
  • Successive overrelaxation, SOR (LT)
  • Symmetric SOR (LT)
• Nonstationary methods
  • GMRES (LT, Saa)
  • MINRES (LT)
  • Conjugate Gradients, CG (LT)
  • LSQR
  • LSMR
  • BiCGStab(l) (Sle)
  • Chebyshev iteration (LT)
  • Quasi-Minimal Residual, QMR (LT)
  • DQGMRES

(Generalized) eigenproblems

• Without subspace
  • (Shift-and-inverted) power iteration (ET)
  • Rayleigh Quotient Iteration [1]
• Krylov subspace methods
  • Implicitly restarted Lanczos (ET)
  • Implicitly restarted Arnoldi (ET)
• Other subspace methods
  • Jacobi-Davidson method (ET)
  • LOPCG
  • Accelerated Rayleigh Quotient Iteration [1]
  • (Shift-and-inverted) subspace iteration (ET)

Singular value decomposition

• Golub-Kahan-Lanczos (ET)

References

• LT: Templates for the Solution of Linear Systems
• ET: Templates for the Solution of Algebraic Eigenvalue Problems
• Saa: GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems (pdf)
• Sle: BiCGstab(l) and other hybrid Bi-CG methods (pdf)
• TB: Numerical linear algebra
• HV: Iterative Krylov methods for large linear systems

[1] Low-priority methods: they are quite expensive per iteration.

[timholy]

An awesome list. Somewhere on my harddrive I have an lsqr.jl file I could contribute, if interested.

[stevengj]

I would say that all the stationary methods are obsolete (they are only useful now as a component of multigrid). CGN (conjugate-gradient on the normal equations) seems hardly used these days; the fact that it squares the condition number makes it useless for large systems. I don't see the point of the simple power method and friends; it seems subsumed by restarted Arnoldi (it is just the special case of restarting every iteration). I don't see the point of implementing non-restarted Arnoldi or Lanczos (these are just the special case of a large restart period).

[ViralBShah]

Great list. I guess the non-restarted methods are more starting points for the restarted methods.

[stevengj]

For the same reason, we only want restarted GMRES (of which there are several variants, if I recall correctly).

[stevengj]

Note that I added LOPCG to the list; this is one of the best iterative eigensolvers for Hermitian systems in my experience (and has the added benefit of supporting preconditioning).

[stevengj]

Also, the solvers for ordinary eigenproblems are merely special cases of the solvers for generalized eigenproblems, so we shouldn't have to implement both variants in most cases. With a little care in implementation, it should be possible to have negliglble performance penalty in the case where the right-hand side operator is the identity, since I(x) = x is essentially free and doesn't make a copy.

[StefanKarpinski]

This is an amazing issue ❤️

[JeffBezanson]

A bcgstab was posted here: JuliaLang/julia#4723

[ViralBShah]

We do not know about its license. For starters, the one from templates should be good.

[jiahao]

As @stevengj pointed out in the original discussion (JuliaLang/LinearAlgebra.jl#29), the state of the art for biconjugate gradients is the BiCGSTAB(l) (or ell?) variant. I purposely left the original BiCGSTAB algorithm off the list for that reason.

[jiahao]

I spent a few minutes transcribing symbol-for-symbol the BiCGSTAB(ell) algorithm in this gist. Unicode combining diacritics make it laughably easy to do this transcription. You can see clearly where when l=1 this reduces to BiCGSTAB.

Obviously this is subject to testing, etc.

[ghost]

A number of comments

1. Stationary methods are important and should be there. They can be used as smoothers in multigrid as well as
   first line preconditioners for other iterative methods.
2. Need to add least square solvers, in particular CGLS and LSQR
3. The solvers I added (bicgstab and gmres) are a direct translation from the book Templates so I do not see any license issue.
4. I plan to add cgls soon

[ViralBShah]

@timholy 's https://github.com/timholy/Grid.jl provides restriction and interpolation operators on grids.

[JasonPries]

With regard to the usefulness of CGN and CGS I would point to this paper:

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.210.62

[tkelman]

I'd be interested to see some of the symmetric indefinite methods MINRES, SYMMLQ, MINRES-QLP (http://www.stanford.edu/group/SOL/software.html)

[thraen]

I tried the gauss-seidel function of the package as smoother for multigrid. The implementation given doesn't account for sparse matrices and is thus not usable in such contexts (large, sparse).
Probably all the methods in stationary.jl have this issue.

[lruthotto]

@thraen: You might want to look into the function sor and other iterative solvers that are optimized for sparse matrices in https://github.com/lruthotto/KrylovMethods.jl. I use the sor a lot as a preconditioner and it is reasonably fast. It makes use of the sparse matrix format.

[thraen]

Thanks a lot for pointing that out!

On Wed, Jun 10, 2015 at 8:39 AM, Lars Ruthotto notifications@github.com
wrote:

@thraen https://github.com/thraen: You might want to look into the
function sor and other iterative solvers that are optimized for sparse
matrices in https://github.com/lruthotto/KrylovMethods.jl. I use the sor
a lot as a preconditioner and it is reasonably fast. It makes use of the
sparse matrix format.

—
Reply to this email directly or view it on GitHub
#1 (comment)
.

[Sisyphuss]

Is it worthwhile to add Aitken's delta-squared process into the list?

[klkeys]

I noticed that this package was listed for Julia GSOC 2016. Mike Innes suggested that I pipe up here if I wanted to help. Is anybody here able and willing to serve as a mentor?

[jiahao]

@klkeys Thanks for your interest. I'd be thrilled to have someone help out with this package!

One of the biggest challenges is to unify the various methods available in similar packages like https://github.com/lruthotto/KrylovMethods.jl and https://github.com/JuliaOptimizers/Krylov.jl. In some sense, this package has been more experimental in nature whereas the other two have been more focused on efficient implementations that work well for practical applications. Now that we (or at least, I) have some experience in writing these algorithms, I'd be very interested to figure out a common API across all the various solvers available.

cc: @lruthotto @dpo

[lruthotto]

I'd be glad to help work out a common API and a convergence of the packages in terms of efficiency. I'm not sure how much time I'll have to work on it myself during this summer, but I'd be happy to provide input as much as I can. Please keep me posted if you start a discussion somewhere.

[dpo]

I'm very interested in working out a good common API as well. I'm currently quite happy with the API in Krylov.jl but I'm of course open to discussing it. It's very unlikely that I'll have free time this summer to mentor though, but I'm happy to participate in discussions when I can.

[zimoun]

Will anyone work on Julia GSOC 2016 ?

It should be nice to have a common API to easily switch between the different implementations or to merge them. For example,

• the preconditioning scheme is different (e.g., M(x) in KrylovMethods.jl and M\x in IterativeSolvers.jl.
• an abstract layer to manage preconditioners and then implement classic ones (see there e.g., sec. 4.4 Preconditioners and Tab. 4 )
• write a "bridge" to be able to transparently use PETSc.jl, other GSOC 2016.

cc: @cortner related issue and JuliaLang/julia#71

[ChrisRackauckas]

I think MINRES should be added to this list.

[lruthotto]

@zimoun: We chose to provide the preconditioner as a function handle that computes the inverse, for example, M(x) = M\x to allow as much flexibility to the user as possible. This way you should be able to use routines from PETSc without an issue (you may have to write a wrapper). I'm not sure you need another abstract layer for preconditioners.

Something we could to in KrylovMethods is to allow M being a matrix and then building the function. We do something similar for A already. Do you think that's a good way to go?

@ChrisRackauckas : There is a basic version of minres in KrylovMethods (https://github.com/lruthotto/KrylovMethods.jl/blob/master/src/minres.jl) which might be a good starting point.

[cortner]

I'm not sure you need another abstract layer for preconditioners.

My own perspective: in optimisation you need more than just M(x). You need to provide M * x, M \ x, <x, My>, <x, M\y>.

But I would argue even in linear algebra you may want to create some abstraction. E.g., is the preconditioner right- or left-preconditioned? Maybe I am overcomplicating it. I do think for Optim.jl it will be useful to use dispatch instead of function handles. I won't get too stressed if this is incompatible with IterativeSolvers.jl, though it might be a shame.

[lruthotto]

I agree that we should distinguish between left and right preconditioners. How about having two variables, for example, Mr and Ml doing this?

Regarding the abstraction issue: Can you give me an example when you need to multiply with the preconditioner?

[cortner]

Regarding the abstraction issue: Can you give me an example when you need to multiply with the preconditioner?

Basically to compute distance, but I only need the inner products for that. So now that you pressure me to think about it, I am not sure whether it was just sloppiness on my part. I will think about it some more, but it is possible that only M \ x, <x, My>, <x, M \ y> is needed.

But don't you need M * x for right-preconditioning?

[stevengj]

For preconditioner (KSP) objects in PETSc.jl, we provide A_ldiv_B! (in-place M \ x), which corresponds to how preconditioners are defined in PETSc. Defining A_ldiv_B! for an object is sufficient to define \ as well, I believe, and it allows the solver to take advantage of in-place operations if it wants. So I would recommend this for IterativeSolvers too.

We also provide a function to do At_ldiv_B! (in-place M' \ x), but currently this is missing from Base. Actually, I just noticed that Base.LinAlg.At_ldiv_B! is defined.

[lruthotto]

@cortner: Thanks for your encouragement and suggestions. I just pushed some changes to the branch and switched more methods over to the new call to applyPrecond. I quite like it: It made the code not much longer and more flexible. Precondioners can now be provided as functions, arrays, and vectors.

What still remains is the option of left vs. right preconditioning in CG and other methods as well as preconditioning in minres (which @dpo has already implemented in his package)

[cortner]

Personally, I like applyPrecond, but I fear that it clashes with Julia naming conventions. Just applyprecond or apply_precond ok?

EDIT: Unless I hear otherwise I now think that the first bullet-point is crucial?! @stevengj: would you be willing to comment?

Although A_ldiv_B! hurts my eyes, this has some advantages:

• it is in Base, which means that if two modules implement it then they can import from Base and don't need to know of each other. Btw, this is an issue I never quite wrapped my head around - is this discussed somewhere in depth? Can anybody comment?
• it is widely understood what it means
• you don't need to implement boiler-plate code for the standard operators (even though this is minimal effort)

[ChrisRackauckas]

Would an Algebraic Multigrid be considered for this package?

[cortner]

I've written a simple wrapper for pyamg which I was hoping to integrate but if you have something better, I'd be very interest d.

P.S.: I now made it public at:
https://github.com/cortner/PyAMG.jl
(this is really just wrapping some documentation and tests around PyAMG + PyCall)

[ViralBShah]

PyAMG is a pretty solid solver.

[ChrisRackauckas]

I think we should get a native Julia version, but PyAMG will work really well in the meantime. Good job!

[cortner]

Thank you - I'll be happy to register the package then.

[lruthotto]

Talking about native julia implementation of mulitgrid methods. @erantreister has started putting out a package here: https://github.com/JuliaInv/Multigrid.jl

We're still working on some details, but it'll be in good shape soon. Always open to feedback/suggestions.

[cortner]

I am now almost certain that we should use A_ldiv_B! or simply \ for application of the Preconditioner. If we do not, then we would have to create a Preconditioning.jl module that just defines the interface without anything else. This is a real short-coming of the module system, which I believe traits are supposed to fix? But by using A_ldiv_B! or \ we can simply overload the operators that are in Base, so there is no problem anymore.

I would really appreciate to get this confirmed (or contradicted) by somebody who has a deeper understanding of Julia than I do.

[andreasnoack]

I really don't see the point in introducing a new a special linear solve function just for preconditioners. Applying the preconditioner is a linear solve and \,A_ldiv_B! are the functions for linear solves and they can be overloaded. Regarding diagonal preconditioning then we have Diagonal in Base.LinAlg for that. If some of the abstractions in Base.LinAlg are insufficient, then please raise issues.

[cortner]

@andreasnoack Thank you - I agree it would be good to stick with this.

[raminammour]

Hello all!

First thank you, I am a frequent user/abuser of this package! (I hope this is a good thread for this post)

1- I was wondering if you have any plans to support the LinearMaps.jl package as a general interface, instead of MatrixFcn which tries to accomplish the same, with less options. I have been using LinearMaps instead of MatrixFcn with IterativeSolvers, and for the most part, it works.

2- When the preconditioners are passed as a function (as in gmres for example), it is better to assume that the inverse is passed directly so Pinv*x replaces Pl\x in the code, or have the two options, maybe. For most of the things I do, the linear operators are large scale and beyond building explicitly.

3- I also cannot wait till this package is integrated with different optimization packages to replace all direct solvers. It will open many doors for large scale optimization.

Again, thanks for the developers,
Cheers!

[syntapy]

Hey, I'd like to work on this for GSOC 2017 if yalls would be fine with that. I think I'd have fun with it. I'm not sure yet whether I'd be able to do GSOC full time since I have work obligations, but if not, I'd still love to do this part time for free

[haampie]

Hi all,

I would like to work on this project as well for GSoC 2017. What are the current priorities? Implementing more iterative solvers? Integration with LinearMaps.jl?

Some ideas I have: implementing (eigen)solvers listed above such as Jacobi-Davidson, BiCGStab(l), IRAM (but there are already some old & open PRs on this); an maybe add some new flavours to existing methods like [thick restarted / augmented subspace / deflated] GMRES. All these latter methods are based on the idea that the Krylov subspace obtained at the end of a cycle contains very useful spectral data which can be exploited in the next cycle; potentially restoring superlinear convergence behaviour observed in full GMRES.

[barche]

From my point of view (as a potential user) I would be interested most in preconditioners (ILU, AMG).

[erantreister]

Hi Bart,
I have a package with some multigrid preconditioners in Multigrid.jl. In particular smoothed aggregation AMG. It's a part of the JuliaInv collection of packages.
Cheers,
Eran Treister.

[barche]

Aha, good to know, thanks!

[GregVernon]

How about the recently developed Alternating Anderson Richardson (AAR) method?

Alternating Anderson-Richardson method: An efficient alternative to preconditioned Krylov methods for large, sparse linear systems https://arxiv.org/pdf/1606.08740.pdf

[mohamed82008]

It would be nice to have an updated roadmap in light of the new developments in Julia's iterative solvers ecosystem.

[kiranshila]

Just submitted a PR for QMR