Nonlinear root solving libraries which accept a Jacobian in band-storage

Question

I'm in search for a library for solving large systems of non-linear equations, similar to MINPACK, but unlike MINPACK, can accept a Jacobian in band-storage.

My Jacobian is sometimes not invertible, so a method like the Levenberg-Marquardt algorithm in MINPACK, which can deal with singular Jacobians is important. Also, it would be nice if the package was in Fortran, but C/C++ is ok too.

I'm solving a system of ~10,000 nonlinear equations, with a non-zero diagonal about ~200 wide, so using band-storage is necessary. Thanks for any suggestions.

Answer 1

A few options are

• KINSOL (Fortran and C interfaces are available; page 18 in the KINSOL documentation says you can use the banded solvers from LAPACK)
• scipy.optimize.root (particularly the large scale nonlinear solvers)
• fsolve in MATLAB allows you to specify the Jacobian as a sparse matrix.
• MATLAB codes from the book "Iterative Methods for Linear and Nonlinear Equations"
• Newton-Krylov method in Fortran or MATLAB

Alternatively you could try and use a nonlinear minimizer like the Ceres solver from Google. This library supports sparse Jacobian matrices, but their Levenberg-Marquardt implementation is primarily aimed at nonlinear least squares problems. An answer over at Stack Overflow explains how you might still be able to use it, but it involves the risk of converging to a minimum instead of the root.

Apart from KINSOL, none of these codes appear to be particularly adapted to banded matrices. In general if you store your Jacobian as a sparse matrix (e.g. CSR format) a back-of-the-envelope calculation says you should be able to reduce your memory consumption from ~100 MB to ~2 MB (assuming you are using 64-bit reals). Compared to a general sparse array, the LAPACK banded storage does not require the integer pointer arrays, so you would have some additional savings (~1 MB). In both cases I think you could expect to see a significant performance increase.