What does LinearSolve Preconditioner do?

Question

I have looked at a few places in order to find an answer to my question - specifically, I have looked at this old Mathematica documentation page, at the Preconditioner option for NDSolve, as explained here, also at this older question and at this unanswered post here, but I still don't understand what the Preconditioner option does. Does it multiply the matrix and the vector in the linear solve by something (or by its inverse)? How and from where does it take arguments? How can I specify a user-defined function and what should be the form of that?

Answer 1

TL;DR Basically, a preconditioner is meant to speed up the convergence of iterative linear solvers. A preconditioner in Mathematica is the inverse of what is called preconditioner in numerical Mathematics.

---

The computational backend for LinearSolve[A,b,Method->"Krylov" ] is SparseArray`KrylovLinearSolve.

More or less,

SparseArray`KrylovLinearSolve[A,b, "Method" -> "BiCGSTAB", "Preconditioner" -> f]

should be equivalent to

LinearSolve[A,b, Method -> {
  "Krylov", 
  "Method" -> "BiCGSTAB", "Preconditioner" -> "ILU0"
  }]

However, calling SparseArray`KrylovLinearSolve directly is usually a bit faster because LinearSolve seems to have some overhead.

Other supported Krylov methods are "ConjugateGradient" (only for symmetric positive-definite matrices) and "GMRES". See the documentation of LinearSolve, section Options, subsection Methods, subsubsection "Krylov".

You can use an arbitrary function f as preconditioner. Here an example for a built-in preconditioner (incomplete LU-factorization without fill-in):

precdata = SparseArray`SparseMatrixILU[A, "Method" -> "ILU0"]
f = x \[Function] SparseArray`SparseMatrixApplyILU[precdata, x]

Other accepted values for the option "Method" in SparseArray`SparseMatrixILU are "ILUT" and "ILUTP".

Another example is this:

f = With[{invdiag = 1/Normal[Diagonal[A]]},
   x |-> invdiag x
   ];

It is the well-known Jacobi preconditioner -- or rather its inverse.

In constrast to the mathematical nomenclature, a preconditioner itself gets applied, not its inverse, in Mathematica. That's very plausible from an algorithmic point of view.

See this and the standard literature about numerical linear equations for more information about what a preconditioner is. Basically, a preconditioner $f$ is meant to speed up the convergence of iterative linear solvers. Sloppily formulated, $f$ is a good preconditioner for the matrix $A$ if 1. $f(x)$ is easy to compute for each vector $x$ and 2. $A \,f(x)$ is close to $x$. So, a good preconditioner $f$ should best be almost an inverse of $A$ - but orders of magnitude faster than computing $A^{-1} x$ with direct methods.