Matrix-Free/Function-Only Linear Solves

Suppose I want to solve a system of equations A*x = b, and I have a function A_func such that A_func(x) = A*x. How could I go about solving the system of equations without constructing the matrix A?

I could construct the matrix A if I really wanted to, but I have written enough in terms of A_func that I would like to avoid it if possible.

I saw some similar topics, but they seemed to be focused on the nuances of specific package implementations.

One option is to use an iterative method such as those from Krylov.jl. Have a look at the docs to choose a method based on the properties of A (if any).

https://docs.sciml.ai/SciMLOperators/stable/

SciMLOperators.jl isn’t perfect yet, but it’s pretty close. It has the most integration with the Julia package system by far too, and it’s only growing. So I’d recommend it.

https://docs.sciml.ai/LinearSolve/stable/

LinearSolve.jl takes AbstractSciMLOperators and will use that with the various Krylov methods. This also means any SciML solver (ODEs, nonlinear solvers, etc.) will also take AbstractSciMLOperators for Jacobian/Hessian definitions.

3 Likes

@leespen1 The relevant documentation page for Matrix-free operators is here if want you to use it in Krylov.jl.

We did two examples with LinearOperators.jl but you can use SciMLOperators.jl, LinearMaps.jl or your own operator if you prefer.

For information, LinearSolve.jl mainly uses Krylov.jl for Krylov methods under the hood.

1 Like

Thanks!

Thanks! LinearMaps.jl works very well for me!