Using preconditioners efficiently

Question

I am trying to numerically solve a linear system of equations of the form

A x = a

where A is really ill-conditioned and a a vector. Both A and a are known and I need to determine x.

Such systems of linear equations arise when considering discretisations of PDEs using spectral collocation methods, for instance. In order to solve this problem, I am trying to use arbitrary precision and iterative methods that use preconditioners.

One such preconditioner, is to find a finite difference approximation to A and a, which I will denote by B and b respectively. I would like to know how to take advantage of this using Mathematica.

So far, my code looks like the following:

g = LinearSolve[B];
f = x \[Function] g[x];
solF = f[b];
LinearSolve[A, a, 
 Method -> {"Krylov", "Method" -> "BiCGSTAB", 
   "Preconditioner" -> f, "StartingVector" -> solF}]

Is there a more efficient way of doing this?

Answer 1

 SparseArray`KrylovLinearSolve[A, a, 
   "Method" -> "BiCGSTAB", "Preconditioner" -> f, "StartingVector" -> solF
   ]

might get you rid of some calling overhead. But improving the quality of your preconditioner will have a considerably larger effect.

If LinearSolve[B] works as a preconditioner, then possibly its incomplete LU preconditioner might work as well. Since it is cheaper to apply then LinearSolve[B], you should give it a try. If I recall it correctly, you have to set up this preconditioner as follows:

data = SparseArray`SparseMatrixILU[B];
f = x \[Function] SparseArray`SparseMatrixApplyILU[data, x];

Note that this will work only if B is a sparse matrix.

If you manage to implement the action of A on vectors as a function fA so that fA[x] == A.x, then you can use SparseArray`KrylovLinearSolve as follows

 SparseArray`KrylovLinearSolve[fA, a, 
   "Method" -> "BiCGSTAB", "Preconditioner" -> f, "StartingVector" -> solF
   ]

This is useful for matrices whose fA action can be implemented without assembling A. If A is a dense matrix of size $N \times N$, then its assembly cost $O(N^2)$ time and memory and A.x costs $O(N^2)$ time. But if A stems from a convolution with a function, i.e. $A \, x = x * u$, then the fast Fourier transform can be employed to implement fA that computes fA[x] in $O(N \, \log(N))$ time.