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


2 Answers 2


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.

  • $\begingroup$ I would have loved to see a minimal working example of this answer. $\endgroup$ Feb 25, 2021 at 14:33

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.

  • 4
    $\begingroup$ This might create local minima that are not solution to the original problem. $\endgroup$
    – nicoguaro
    Sep 22, 2019 at 15:42
  • 3
    $\begingroup$ Nonlinear solvers are almost always a better choice than applying optimization algorithms to the residual $\endgroup$
    – whpowell96
    Feb 19, 2020 at 17:27

Your Answer

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

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