# Generic Iterative Solvers: Krylov's methods example

**URL:** <https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086>\
**Category:** General Usage\
**Tags:** question, package\
**Created:** [June 1, 2022, 9:36pm UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086 "2022-06-01T21:36:21Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [June 1, 2022, 9:36pm UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086/1 "2022-06-01T21:36:21Z")

</div>

Hi,

This is related to a discussion within the lengthy tread that spun off Yuri’s blog post. I don’t want to make that thread even longer with off-topic material.

In that thread, when discussing Julia’s ecosystem @juliohm asked why we have 3 packages for iterative solvers, and he cited IterativeSolvers.jl, Krylov.jl, and KrylovKit.jl.

As I’m very interested in iterative solvers and this kind of “package ecosystem” issues, I dug a little bit into this and I think I found the answer in KrylovKit docs:

> KrylovKit accepts general functions to represent the linear map or operator that defines the problem, without having to wrap them in a [`LinearMap`](https://github.com/Jutho/LinearMaps.jl) or [`LinearOperator`](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl) type. Of course, subtypes of `AbstractMatrix` are also supported. If the linear map (always the first argument) is a subtype of `AbstractMatrix`, matrix-vector multiplication is used, otherwise, it is applied as a function call.

The magic is, of course, achieved via multiple-dispatch.

See the extremely simple and elegant [apply.jl](https://github.com/Jutho/KrylovKit.jl/blob/master/src/apply.jl) code in that package, which contains lines like:

```julia
apply(A::AbstractMatrix, x::AbstractVector) = A * x
apply(f, x) = f(x)

```

which are then used:

```julia
function GMRES(operator, ...)
    y₀ = apply(operator, x₀)
    (...)
end

```

I believe this is both a pretty elegant example of multiple dispatch (actually, I’m going to make this my example of multiple-dispatch in [my blog’s tutorial](https://www.matecdev.com/posts/julia-tutorial-science-engineering.html) (I still don’t have one) and that this is also the most “Julian” way to write an Iterative solver, so I guess I’m going to start recommending this package as well in my [numerical computing tutorial](https://www.matecdev.com/posts/julia-numerical-computing.html) (after testing it, of course).

Now, a few thoughts about the ecosystem.

As Julia developed, many great package names were available, like “IterativeSolvers”, which will win Google rankings in a trivial way (by direct keyword hit). This package is also maintained by an organization called JuliaLinearAlgebra, which lead one to think that, somehow, IterativeSolvers.jl is the “official” Julia package for iterative solvers.

However, I believe InterativeSolvers.jl is simply the oldest package, and, for some reason, the creators of both Krylov.jl and KrylovKit.jl decided to create a separate package instead of creating PRs to the existing one, maybe because that was easier.

As for Krylov.jl, they list in their docs:

> All solvers in Krylov.jl have an in-place version, are compatible with **GPU** , and work in any floating-point data type.

So I guess it was one of those features that made them deviate from the other packages.

I guess in situations like this, the different package authors can learn from each other, and yes, eventually it would be desirable that someone merges all the contributions into a single package, to avoid duplicated efforts, and to give all users the best possible experience.

However, I’m not really sure how this could be accomplished. Maybe it’s just a matter of letting the individual contributors decide to which package they are willing to contribute to, and let the more motivated/friendlier package maintainer summarize the best efforts.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [June 1, 2022, 10:12pm UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086/2 "2022-06-01T22:12:47Z")

</div>

The advantage of KrylovKit is that it allows vectors not `<: AbstractArray` plus being compatible with GPUs. This was my main motivation for [using it](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials2b/#d-Swift-Hohenberg-equation-(non-local)-on-the-GPU,-periodic-BC-(Advanced)) in BifurcationKit.jl. The con is that it is out of place.

Now, we have [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [June 2, 2022, 1:54am UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086/3 "2022-06-02T01:54:43Z")

</div>

Cool, it looks like you’ve wrapped every package for solving linear systems out there!

What I truly like about KrylovKit.jl is that you accept “anything” that has an associated “apply” function to it, and you then provide such apply functions for matrices or functions. A user could even define it’s own type, code an “apply” function, and forward the type directly to your solver. I guess that is truly generic.

For example, I can do:

```julia
using KrylovKit

iter_A(x) = [-2x[1] + x[2]; x[1:end-2] - 2x[2:end-1] + x[3:end] ; -2x[end] + x[end-1] ]
xt = rand(1000)
b = iter_A(xt)

sol = linsolve(iter_A,b)

```

Now, this runs, but it doesn’t converge (I’m not sure what’s your default algorithm here, but that was not my main point anyway)

At first sight, this kind of interface doesn’t seem to be possible with LinearSolve:

```julia
using LinearSolve

iter_A(x) = [-2x[1] + x[2]; x[1:end-2] - 2x[2:end-1] + x[3:end] ; -2x[end] + x[end-1] ]
xt = rand(1000)
b = iter_A(xt)

prob = LinearProblem(iter_A, b) 

ERROR: MethodError: no method matching size(::typeof(iter_A), ::Int64)

```

So I guess that in LinearSolve you still need to set up something like LinearMaps somehow?

Anyway, I’m sticking to whatever gives me unrestarted GMRES with a custom number of iterations, and it could be KrylovKit if I get to configure it correctly, or IterativeSolvers as I’m already using it.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [June 2, 2022, 6:39am UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086/4 "2022-06-02T06:39:21Z")

</div>

I did no write that package

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 2, 2022, 10:26am UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086/5 "2022-06-02T10:26:14Z")

</div>

> [@martin.d.maas](#):
>
> At first sight, this kind of interface doesn’t seem to be possible with LinearSolve:
> 
> ```julia
> using LinearSolve
> 
> iter_A(x) = [-2x[1] + x[2]; x[1:end-2] - 2x[2:end-1] + x[3:end] ; -2x[end] + x[end-1] ]
> xt = rand(1000)
> b = iter_A(xt)
> 
> prob = LinearProblem(iter_A, b) 
> 
> ERROR: MethodError: no method matching size(::typeof(iter_A), ::Int64)
> 
> ```
> 
> So I guess that in LinearSolve you still need to set up something like LinearMaps somehow?

You need to use a SciMLOperator. The interface package is being separated out soon into [GitHub - SciML/SciMLOperators.jl: SciMLOperators.jl: Matrix-Free Operators for the SciML Scientific Machine Learning Common Interface in Julia](https://github.com/SciML/SciMLOperators.jl), and the documentation won’t be up until the end of the week but will be mostly the same as [https://diffeq.sciml.ai/stable/features/diffeq\_operator/](https://diffeq.sciml.ai/stable/features/diffeq_operator/) . The problem is that LinearMaps ended up being useless in most solver contexts because it is not an expressive enough API, so we just couldn’t use it in differential equation contexts and necessitated creating a whole interface, which LinearSolve.jl expects.

Anyways, we planned to have it finalized by the end of the week and are on track so we’ll have an update rather soon.

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [June 3, 2022, 9:27am UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086/6 "2022-06-03T09:27:57Z")

</div>

Thank you, Chris. I’ll keep an eye on LinearSolve then and wait for these updates. It looks great as a reference package implementing every possible method out there.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 3, 2022, 1:09pm UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086/7 "2022-06-03T13:09:07Z")

</div>

> [@martin.d.maas](#):
>
> What I truly like about KrylovKit.jl is that you accept “anything” that has an associated “apply” function to it, and you then provide such apply functions for matrices or functions. A user could even define it’s own type, code an “apply” function, and forward the type directly to your solver. I guess that is truly generic.

Doesn’t IterativeSolvers.jl equivalently accept any object that defines `LinearAlgebra.mul!` (which has the advantage that it gives you an in-place interface)? Similarly for Krylov.jl.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 3, 2022, 1:41pm UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086/8 "2022-06-03T13:41:17Z")

</div>

> [@stevengj](#):
>
> Doesn’t IterativeSolvers.jl equivalently accept any object that defines `LinearAlgebra.mul!` (which has the advantage that it gives you an in-place interface)? Similarly for Krylov.jl.

Yes. The bigger issue is that the preconditioner interfaces are not aligned well though.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [June 3, 2022, 2:41pm UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086/9 "2022-06-03T14:41:15Z")

</div>

The above answers show you that the real bonus of KK is the state space which is not an AbstractArray.  
Another pro of this library is the eigensolver.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 3, 2022, 2:45pm UTC](https://discourse.julialang.org/t/generic-iterative-solvers-krylovs-methods-example/82086/10 "2022-06-03T14:45:49Z")

</div>

> [@rveltz](#):
>
> Another pro of this library is the eigensolver which, to me, has no analog apart from [NonlinearEigenproblems](https://nep-pack.github.io/NonlinearEigenproblems.jl/). The other iterative solvers focus on hermitian problems.

The other main packages for non-Hermitian iterative eigensolvers are [ArnoldiMethod.jl](https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl) and [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl)