# Using preconditioners efficiently

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?

• Information is a bit rare here. In what form is A given? Is A a dense matrix, a sparse array or only a function that implements matrix-vector multiplication? – Henrik Schumacher Jan 9 at 16:26
• A is really dense, but B is rather sparse. – user12588 Jan 9 at 16:28
• Is there a way to implement the action of A on vectors without assembling it? – Henrik Schumacher Jan 9 at 16:31
• Unfortunately, no, there is no way that I know of implement the action of A on vectors... – user12588 Jan 9 at 16:56

 SparseArrayKrylovLinearSolve[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 = SparseArraySparseMatrixILU[B];
f = x \[Function] SparseArraySparseMatrixApplyILU[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 SparseArrayKrylovLinearSolve as follows

 SparseArrayKrylovLinearSolve[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.