A\b vs LinearSolve.jl at JuliaCon 2022

A\b is fine. It’s not going to go anywhere. Deprecating it would be silly. But, if you’re going to write a package that needs to allow customizability of the linear solver in order for the algorithm to be optimal in all possible cases, then A\b is a bad interface to build off of. At no time do I advocate for getting rid of A\b, in fact I say that it has quite a bit of smarts in it. But, it’s just not good enough for many applications, and it will never be a standard library since it would have to allow hooks for a ton of (very large) artifact dependencies and a lot of other cruft that doesn’t deserve to be in the standard library.

It actually goes pretty deep, and it takes a bit of implementing to really see why. One naïve thing you can do is say, okay, what about \(A,b;kwargs...) as the interface? But you end up with a lot of weird behavior really quickly. For example, if you do \(method(A),b;abstol,reltol), you might think that lu would work right off the bat, i.e. \(lu(A),b;abstol,reltol). But psych, if you do that today you’d get a MethodNotFound() error because the \ dispatch on LU does not match any keyword arguments. So you can go through and either

  1. Write out (and ignore) all keyword arguments that would make a common interface (pretty messy for Base).
  2. Have every \ on a factorization dispatch splat keyword arguments (this is an awful idea for the standard library because now it won’t throw errors when keyword arguments are misspelled, it’ll just ignore them :sweat_smile: You can then create separate error parsing and…).
  3. Target the standard interface at a different level to leave \ alone to avoid having to add more methods (and deal with ambiguities).

If we did go down this route (which is one I championed for about 5 years BTW), you still would have to have all downstream libraries adopt it, and it’s not necessarily easy. You’d have to create a lazy IterativeSolvers.cg(A) that doesn’t do anything until \(IterativeSolvers.cg(A),b;abstol,reltol). Cool, that’s just a bit of code churn. But then some of the harder to wrap libraries can have increasingly difficult issues with how to make that call nicely lazy because you’re also then trying to reuse that syntax for reuse.

What I mean is that luA = lu!(A); \(luA,b;abstol,reltol); \(luA,b2;abstol,reltol) reuses the LU-factorization of A. What is effectively going on now is that you’re mixing the notation of lazy call for method selection with the notation for factorization reuse. And if you actually want that to be efficient, then you need to make sure it’s actually a sufficient syntax for factorization reuse in all cases. In order for that to be the case with sparse matrices, you would need to know the sparsity pattern of the previous matrix in order to know whether to save the symbolic factorization. So instead the syntax would need to be something like \(lu(luA,A2),b;abstol,reltol) where luA is either the previous factorization (it cannot be the matrix itself because only the factorization object holds the symbolic factorization result, and it also needs to check the sparsity pattern).

Now I bet you $20 that if you had to go around telling people “we’re deprecating A\b for \(lu!(nothing,A),b;abstol,reltol) because it’s a cleaner syntax to generalize to iterative solvers and optimally using sparse matrices purely by dispatching on method and A!”, people would be pissed. No, that is absolutely worse for most cases!

That’s why it’s a separate interface, built off of (in many cases) tools first made in Base.LinearAlgebra, that is focused on a problem/solve syntax that works no matter the choice of A and method. The reason is because that’s exactly the kind of thing that a generic library needs to use: a generic library needs to handle any input someone throws at it with as minimal special casing as possible. LinearSolve.jl does that special casing.

But, I and hopefully many others will continue to use A\b in the REPL. It’s good enough for that, but not for everything else we need to build.

(There are other things I’m skipping too, like how for fully non-allocating lu!/ldiv! you need to allocate and store the workspace / pivot vectors somewhere, which is handled only when you get to the \(lu!(nothing,A),b;abstol,reltol) idea, so the sparse factorization stuff actually improves the dense case as well, but these are details)

5 Likes