Which SciPy nonlinear solver when Jacobian is analytically known and sparse?

Question

I have a nonlinear function fun with n inputs and n outputs. I also have a function jac which calculates the Jacobian, which is sparse. Which SciPy method, if any, can I use to find the root of fun?

---

I am surprised that I can't find a method in scipy.optimize that takes a Jacobian in sparse matrix format. But I believe this can be simply solved by the Newton's method, and so I am confused if such a method does not exist. There are many quasi-Newton methods in this package that estimate the Jacobian, but they do not seem quite right. The only thing that looks similar is newton_krylov which documents something about sparse preconditioners. Is this what I should be looking at? If there is no such method, should I just do something like the following?

while not some_stopping_condition(x):
    x -= scipy.sparse.spsolve(jac(x), fun(x))

I would prefer to use an existing SciPy function in case I do something stupid. Notably, as it happens, I believe my particular problem is poorly conditioned, and I think I would benefit from being able to use scipy.sparse.lsqr or lsmr which takes a damping parameter.

Incidentally, the problem is solving large-displacement elasticity equations on a 2D membrane in 3D space with FEM. The problem is sparse because it is FEM. The problem is nonlinear because of the large displacements (and the corresponding Green-Lagrange strains). The poor conditioning comes from the fact that a small perturbation of any node in the direction perpendicular to the surface creates little change in strain energy (and the flatter the surface, the smaller the change).

Answer 1

I believe newton_krylov may be what you are looking for. The documentation is confusing but I think you can do what you want. There doesn't seem to be any way to give an analytic Jacobian for use in a Krylov solver in this software, but if you give it a function that takes in $x$ and returns a sparse matrix in the method argument then it looks like that should work. You may run into problems with the linear solve though because it looks like it's using a direct solve if you give it a sparse matrix, which can be a lot slower than Krylov methods depending on the problem.

Answer 2

You can consider your root finding problem as a minimisation problem by trying to minimize the least squares of your function :

$$\min_{x\in \mathbb{R}^n} \|F(x)\|_2 \, .$$

Then solve this problem with the optimize.least_squares solver of SciPy. With this solver you can provide the sparsity information with the 'jac_sparsity' input (and if you do so it will automattically use lsmr) or even give the analytical Jacobian as a callable to the 'jac' input.