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?

  • $\begingroup$ 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? $\endgroup$ Jan 9, 2020 at 16:26
  • $\begingroup$ A is really dense, but B is rather sparse. $\endgroup$
    – user12588
    Jan 9, 2020 at 16:28
  • $\begingroup$ Is there a way to implement the action of A on vectors without assembling it? $\endgroup$ Jan 9, 2020 at 16:31
  • $\begingroup$ Unfortunately, no, there is no way that I know of implement the action of A on vectors... $\endgroup$
    – user12588
    Jan 9, 2020 at 16:56

1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.