I just end watching LinearSolve.jl: because A\b is not good enough by Chris Rackauckas himself and I’m a bit confused. When I get basic points why such packages was developed, somethings bugs me about A\b being “deprecated”.
I understand why LinearSolve.jl was created, factor 2 in solving differential equation is something that I cannot ignore, but I just like A\b because of its syntax.
Also, it seems to me that even when LinearSolve.jl is package that unified dozens of packages under the hood, simple numerical computation in Julia become to complicated for my taste, because you have to many options to choose. Maybe I’m just to stupid to learn Julia and its ecosystem well at this point in time.
A \ b has never been optimal. It has always been the case that there are better alternatives. In particular, you can often use knowledge about the specific matrix (structures, size, sparsity, …) to pick a better solver.
However, in some cases A \ b is good enough – not all programs spend the majority of time in the linear solver. A \ b is not being deprecated, and you can still use it.
One thing worth noting is that we will probably at some point make A\b smarter. There’s no reason A\b needs to be slower than the default LinearSolve solver. The big advantage of LinearSolve is that it provides an interface for customization, but that’s a fairly advanced use case. For 90% of users, we can make A\b as good.
From what I understand form Rackauckas talk, the biggest stumbling block is that some algorithms need more information that knowing A and b in A\b, like relative error of the method (result?). People with much more knowledge of Julia internals are need to answer how serious problem it is.
But, when in the cases that all information is in A and b making some multiple dispatch magic to make algorithm smarter doesn’t seems impossible. Are they some catches here?
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
Write out (and ignore) all keyword arguments that would make a common interface (pretty messy for Base).
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 You can then create separate error parsing and…).
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)
You can then create separate error parsing and…).