API for \

[ChrisRackauckas]

It would be nice to have IterativeSolvers.jl match the interface of other linear solvers so that way \ will solve using one of these methods. A discussion on this can be found here:

https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/5

If I'm not mistaken, all that would need to be done is the creation of types like K = CGSetup(A,opts...) (name pending of course) with a dispatch on \ which will call the correct iterative solver. This would then match using factorizations:

K = lufact(A)
x = K \ b

vs

K = CGSetup(A,args...)
x = K \ b

and will make it easy for package developers to incorporate IterativeSolvers.jl as a replacement for direct methods (without even requiring it as a dependency).

[lopezm94]

I like this, are you thinking of something like a new package with a purpose much like the one of RecipesBase? A unified API for all the linear algebra packages out there.

[ChrisRackauckas]

RecipesBase is not the best example because it has a very specific implementation for plotting. But yes, something like that, but more like MathProgBase (JuliaOpt/JuMP) and DiffEqBase (JuliaDiffEq/DifferentialEquations). I'll have a blog post explaining this all in some detail soon.

[andreasnoack]

I kind of like that \ would generally be "solve linear system" and that you could specify various methods to do that by wrapping the left-hand-side. However, in the case of direct solvers, this is also a matter of dividing computations into subcomponents where you can possibly reuse the expensive part of the computation (the factorization). This is less obvious for iterative solvers so it could be argued that

K = CGSetup(A,args...)
x = K \ b

is a complicated way of writing

x = CG(A,b,args...)

without compuational advantages.

I think that https://github.com/JuliaParallel/PETSc.jl has taken this approach, though, so it might be worth taking a look there.

[ChrisRackauckas]

without compuational advantages.

It's not about the computational advantages. I explained this in more detail in a blog post. If IterativeSolvers has a common interface that matches other interfaces, package authors would not be required to write a different set of code to use IterativeSolvers.jl's methods, and so not only users but also package authors would greatly benefit.

[jiahao]

I had a short discussion with @JaredCrean2 @stevengj last year about this specific API issue, but I don't think it was conclusive.

To summarize, the design problem as I view it, is really one about object design and what the left operand of \ could mean in a generic function system, and the fact that in iterative methods we want \ or solve to express the interaction between three [or even four] different types of data - A, b, algorithm specification, and possibly intermediate state of the algorithm.

a) Originally, I thought of as A\b as purely a way to spell "solve Ax = b". \ means "solve a matrix problem with its operands", in which case it is natural to interpret "A' \(A, b) is naturally extended with keyword arguments modifying the meaning of "solve", for example, i) \(A, b, method=:CG, atol=...) specifying an algorithm and roundoff tolerance, or even ii) \(A, b, method=CGSolver(atol=...), wrapping roundoff tolerance into the specification of the algorithm. CGSolver is now an object with its own fields and possibly even its own API. An advantage of this approach is being able to spell the propagation of intermediate state quite naturally, passing in a CGSolver(data=...), say, or even as a different.

The advantage of this spelling is that it keeps A, b, algorithm specification, and possibly intermediate state separate.

The disadvantage is that there is no meaningful infix expression and hence no natural way to "hijack" an existing solve written as A \ b to mean "do an iterative solve instead".

b) @JaredCrean2 presented to me last year an entirely different idea exposed in PETSc.jl. \ now means KSPSolve, i.e. "when A is a KSPSolve object, solve the system Ax = b using a Krylov subspace method". The first operand is now a Solver object and no longer a matrix or matrix-like function. The spelling then looks something like KSPSolver(A, CGSolver(options)) \ b, using a notation consistent with (a-ii).

The advantage of this spelling is that A \ b is still a meaningful infix expression: A can be substituted by a KSPSolver. Arguably, it is the only meaningful infix design - A \ Solver(b) is not sensible.

The disadvantage is that KSPSolver now takes on characteristics of a God object - all the details of how to calculate and its intermediate data all has to be stuffed into this object. In particular, KSPSolver \ b is no longer a pure operation (in the functional sense), since the internal state of the LHS will change (although IIRC it is not publicly accessible from the PETSc API).

My subjective opinion (which has changed since last year): the evidence points to (b) as the prevailing convention. To me, this looks like a pattern that is natural in class-based object-oriented programming, but not necessarily one that is natural in a generic function-based system. The convenient infix \ spelling doesn't work the moment you want any more complicated specification of what to do, requiring instead the prefix \(A, b, options=...) form. Since it is easy to define the necessary method to support (b) with the (a) design, say:

\(::Solver, ::b) = \ (Solver.A, b, solver=Solver(Solver.alltheotheroptions))

we could do it, and might as well do it.

[ChrisRackauckas]

Since it is easy to define the necessary method to support (b) with the (a) design, say:

That sounds like a good plan. I think the (b) design is necessary since otherwise it's really hard to build packages on top of. With the (b) design, a user just passes a function, and the package just needs to use that function. If it's done with a bunch of kwargs, the package then needs to keep those kwargs around, passing them through, maybe dividing them up if there were some other kwargs mixed in an IterativeSolvers.jl fails if giving an unknown kwarg: it just is a disaster.

So I have already setup part of DifferentialEquations.jl to work with the design (b). You can see from here that it gives a nice clean multi-package API:

http://docs.juliadiffeq.org/latest/features/linear_nonlinear.html#Linear-Solvers:-factorization-1

and it was dead simple to actually put it into affect (here's the only extra line of code in the actual method

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/rosenbrock_integrators.jl#L34

). I'll see if I get around to an IterativeSolvers.jl solver which implements at least style (b) to really get this rolling.