feature request to pass in user-defined matvec, or add simple demo of such to tutorial

[ahbarnett]

Dear Devs,

Thanks for implementing GMRES in Krylov.jl. This is a suggestion: any user is looking for the simplest way to use their own matrix-vector multiply, ie function f(x) that implements Ax. (This is probably the 2nd most important use case, sparse matrices being the 1st.) Your documentation
https://docs.juliahub.com/Krylov/0fcC3/0.7.11/factorization-free/
highly recommends using LinearOperators.jl to create a LinearOperator object. However, to do this for is extremely messy, involving making a 5-argument mul! operator first:
https://jso.dev/tutorials/introduction-to-linear-operators/
It needs another package and is simply too confusing/general; I gave up!

I did get this route to work (still needs another pkg):

A = rand(4,4)
f(x) = A*x              # my dummy matvec, toy example of user function
using LinearMaps
opA = LinearMap(f, 4)
using Krylov
b = rand(4)
gmres(opA, b)

However, there should be a simpler route, just passing in the function f to gmres. No operator dimension is needed since b determines it (and in the rectangular case an N arg could be added). This is how MATLAB does it:

>> A = rand(4); f = @(x) A*x;
>> b = rand(4,1);
>> gmres(f, b)

Could this interface feature be added, to reach the MathWorks level of good design and usability?
Failing that, could a simple example be added to your tutorial of how to achieve the above without installing another pkg?
Otherwise, it is hard to attract new users...

Thanks so much! Alex

[amontoison]

Hi @ahbarnett!

Your link of the documentation is very old (version 0.7.11), please check the last version of documentation with the link available in the README.
You just need create a 3-args mul! with LinearOperators.jl, not a 5-args mul!.
Can you check again with the last version of the documentation?

Are you not allocating a new vector each time that you are performing a matrix-vector product with our operator? It's something that we want to avoid because you can't solve large linear problems if you allocate a vector of 1 000 000 coefficients at each iteration.

[amontoison]

@ahbarnett can you also explain what kind of linear operator you want to implement?
With your example, you can directly call gmres(A, b) and it works.

[ahbarnett]

Thanks for quick reply. I hope this clarifies my question/request:

1. The docs are hard to find, since you have to spot a tiny unlabelled 3-blue-bar icon in the top-right of this page to access the menu. A suggestion is to have the menu clearly visible on the left, as in, eg, https://docs.readthedocs.io/
   Anyway, I checked that your current docs are just the same as I indicated, and point to the same tutorial page
   https://jso.dev/tutorials/introduction-to-linear-operators/
   which has a complicated explanation involving a 5-arg mul! op (which only has 4 args; go figure). Most users will be lost trying to figure out how to call GMRES on their matvec function.

2. Calling a user-defined matvec function is a standard use case for Krylov methods. I have to disagree with you that, unless the operator is extremely simple, allocation time in most scientific applications will be negligible (even in Julia :). Eg alloc size 1e6 vector is <<1ms, but FFT on this size is 20ms. And the fast multipole method (FMM, an algorithm commonly combined with GMRES in my world) will be >1sec.

3. In my comment I explained that the f(x) I gave is just a dummy example to be replaced by a real-world user-supplied matvec function. In the real-world situation there is no "A". Eg, imagine you wrapped some Fortran code to apply A. (I would not be asking the question if it the dense or sparse A were available.) Put another way: this is a question of interface design. How can a user most simply call your solver using their matvec function? (without having to read complicated docs from other pkgs?)

Addressing this would open up the user base in scientific computing. Thanks, Alex

[amontoison]

Thanks @ahbarnett for all your comments, give me a few days and I will refactor the documentation about linear operators. 😃 A few remark about your answer:

1. I am not to understand what is the tiny unlabelled 3-blue-bar icon in the top-right of this page that you are talking about.
   The only place where where we can add a link for the documentation of a Julia package is in the README.md.
   
   The documentation is displayed as follow when I click on the link in the README:
   
   Can you elaborate your suggestion?

2. I am not talking about allocation time but storage. Each Krylov method was highly optimized (see here) to avoid unnecessary allocations. With a sparse linear systems of 10 000 000 variables, one n-vector is 762.939 MiB and on GPU you will saturate your VRAM after a few iterations. Allocations on GPU are also much slower than CPU.
   Note that I understand your point, I just explain why I said this comment initially. 👍

3. Ok I clearly understand this part and I will on it asap. A quick question about your application, is it possible to do y ← A * x in practice or you can't provide the vector that will store the result of A * x?

[ahbarnett]

Dear Alexis,
Great - just a minimal example of calling a "matvec" function (like I did) would be great.
Re 1) the 3-bars menu button is on the page you get to when you click on the "Docs stable" badge. Ie, here: https://jso.dev/Krylov.jl/stable/
It's hard to find that menu. Maybe there is a config for your HTML generator that adds a menu/contents on the left, like the readthedocs page I gave?
Re 2) It's true storage is a concern, but that's because eg GMRES needs to orthogonalize against all prior vectors (up to "restart").
Re 3) I guess there's two possibilities:
i) the user wants to use matvec!(x) which replaces x with Ax.
ii) the user wants to use y = matvec(x) which allocates a new vector.
In MATLAB-land the 2nd is assumed, so that's a decent starting point.
Having both types of interface would be ideal.

Any documented interface adjustments that make using such matvec functions easy with GMRES, etc, would be amazing, and really improve the package accessibility.
(I am currently using the LinearMaps route but it feels a clunky).
Thanks, Best, Alex

[amontoison]

@ahbarnett I added a matrix-free example with a FFT a few days ago.
I don't know if it's still relevant now but I share the link:
https://jso.dev/Krylov.jl/dev/matrix_free/#Example-3:-Solving-the-Poisson-equation-with-FFT-and-IFFT