gempy
gempy copied to clipboard
gempy-project
Improve sparse linear system solver
Is your feature request related to a problem? Please describe.
@Leguark said that, based on a conversation with @denise91, improving the performance of sparse linear system solves might be as simple as choosing a different solver in scipy. My guess is that it is by default using a direct method such as LU decomposition and that, depending on the properties of the linear systems, we may be able to apply something more efficient.
Describe the solution you'd like
Consider the options in scipy and document our choice of linear solver.
Describe alternatives you've considered
Another option would be uisng the linear solver suite in PETSc via petsc4py. I think this is only worth considering when we need to scale on large distributed memory systems (MPI). scipy should be the right place to start.
Additional context I am missing plenty of background from @Leguark who I think has already done some work on this.
Related issues #323
#252 The current sparse implementation is broken. It is due to how I am passing the fault matrices now (after the octtree prototype)
@Leguark, can you get us up to speed on what you've already done? Also has @denise91 already looked into this?
I am looking for a test case. I found a single occurence of using a sparse solver: https://github.com/cgre-aachen/gempy/blob/dce8f68fea46fddf4014f8765311a0d82785043c/gempy/core/theano/theano_graph_pro.py
That file contains this class:
class SolveSparse(T.Op):
#itypes = [T.dvector]
#otypes = [T.dvector]
def make_node(self, x, y):
x, y = as_sparse_variable(x), as_sparse_variable(y)
assert x.format in ["csr", "csc"]
assert y.format in ["csr", "csc"]
out_dtype = theano.scalar.upcast(x.type.dtype, y.type.dtype)
return theano.gof.Apply(self,
[x, y], [T.tensor(out_dtype, broadcastable=(False,))])
# [theano.sparse.type.SparseType(dtype=out_dtype,
# format=x.type.format)()])
def perform(self, node, inputs, outputs):
from scipy.sparse.linalg import spsolve
(C, b) = inputs
b.ndim =1
weights = spsolve(C, b)
outputs[0][0] = np.array(weights)
Also perhaps related, here is another scipy linear solve call (but not sparse): https://github.com/cgre-aachen/gempy/blob/0204d522230a9970a4f3da1202dade7c75f920bd/notebooks/experimental/GemFe.ipynb
I am still waiting on background info from at least @Leguark , maybe also @denise91 before I dig too much into this.
Dear Alex, The default solver is a direct solver and from past experiences, they are much slower than iterative solver, in case you have large matrices. You can find here a list of alternative solvers: https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html If you have a symmetric matrix it is usually the best to start with the CG solver (scipy.sparse.linalg.cg()).
Yes, iterative solvers will always be faster (at the cost of only approximately solving the system). The problems are with robustness, which very much depends on the properties of the matrix. CG will work if the matrices are SPD. It's more likely that we will be able to use a more robust Krylov method such as GMRES out of the box.
For the issue at hand, I need information about the systems GemPy needs to solve. If there are some representative test cases, then I can probably figure out how to analyze the matrices. The first thing I will try is probably just throwing GMRES at them before getting any more fancy. Also we will eventually have to talk about preconditioning, which is always problem-dependent.
Are @javoha or @Leguark on the hook for providing a test case? @Leguark had mentioned some reasons to wait, but I don't quite understand (and then I'm not sure if this issue is still slated for GemPy 2.2) and it'd be good to write it down here.
Please re-assign me once there is a test case. As we discussed offline, I don't have that expertise.
Also if my memory tells me right, then we are also waiting until some other big changes are done before working this issue. It would help to have those documented here. On my side, I'm just waiting on the test case.
We are going to use symbolic solver instead sparse at least for the mid term. #554
