Nonlinearsolve.jl: Fast and Robust Solvers for Nonlinear Equations in Julia
arxiv.org
arxiv.org
https://coin-or.github.io/Ipopt/
IPOPT isn't tailored to NLEs specifically, but it does solve NLPs well. It also uses some amazing Fortran linear algebra routines from Harwell (MA57).
This article is about solving nonlinear equations (not differential equations, not sure where you got that from). All NLP optimizers can solve nonlinear equations — it’s a special case where the objective is constant.
Ipopt is not a convex solver so am not sure what convex optimization you are referring to. It is a general nonlinear solver, which covers nonconvex problems as well (I worked on nonconvex nonlinear programs for a decade and it was my primary solver)
Also all nonlinear equation systems are nonconvex. (A convex program requires equality constraints to be linear)
Maybe you have something more particular in mind when you say "systems", but not all nonlinear functions are non-convex. Least squares, for example, is nonlinear and convex.
Also note that IPOPT, while wonderful, is a local solver. It may not be limited to convex problems, but those are the only ones it's guaranteed to solve to optimality.
> The feasible region {x | f(x) = 0} is nonconvex no matter whether f is convex.
For some positive integer n and for the set of real numbers R, consider closed (in the usual topology for R^n), convex set A a subset of R^n. Define function f: R^n --> R so that f(x) = 0 for all x in A, f(x) > 0 for all x not in A, and for all x in R^n f infinitely differentiable at x. It is a theorem that such an f exists.
Then { x | f(x) = 0 } = A and is both closed and convex in contradiction to the claim.
But consider the function f(x) = max(0, g(x) - c) where the following holds:
* g(x) is nonlinear, positive definite in x, and convex.
* c > 0
Then f(x) is nonlinear and convex (it's the pointwise maximum of convex functions), and the set {x : f(x) = 0} is a convex set.
I’m aware that ipopt is a local solver.
> I was talking in terms of convex optimization
I'm not sure what subtlety you're pointing to here. As far as non-convex problems, though, my point is that IPOPT isn't special in this regard. Any convex solver can be a non-convex solver if you don't care about global optimality.
Aside: Structurally I’m not sure how this would be true.
Convex solvers have very specific forms. For instance a QP solver requires a very particular form and does not admit any arbitrary non convex form except for one: the non-PSD Hessian which is the concave problem.
My point is that all NLEs power inside an optimization problem gives rise to a non convex optimization problem with no guarantee of a global solution. So convex optimization is not applicable here.
IPOPT is a bit of a stretch but I do know that there are some people in the mathematical programming space that use IPOPT's constraint satisfaction piece for solving a nonlinear system in a pinch, but it's not the main use case of IPOPT which is nonlinear optimization with nonlinear (in)equality constraints. This case has no optimization and is just nonlinear equality constraints, which you can then specialize quite a bit on. But while it is considered common knowledge in numerical circles to not use a nonlinear optimizer to solve nonlinear systems, I don't have a good canonical benchmark to point to on that, so we might as well make this it.
For some comparison (this is not a benchmark me or anyone on this paper have written), see https://juliapackagecomparisons.github.io/comparisons/math/n.... Specialized solvers for NLLS pretty much always beat general optimization solvers.
IPOPT is a nonlinear optimization tool, not a nonlinear solver. Generally from what I've seen it's not a good idea to solve nonlinear systems with a nonlinear optimizer, but you can phrase it as a "constraint satisfaction problem", i.e. a nonlinear optimization with a trivial `maximize 0` loss function but with equality constraints to be satisfied `f(u) = 0` and allow it to do that. IIRC the sparse symmetric linear solver is only used in the optimization process as it's that part where it can guarantee a sparse symmetric linear system, the constraints themselves `f` would have to be symmetric to reuse that in the Newton method which isn't true for any of the benchmarks. As such, I don't think you'd get most of IPOPT's advantages even showing up in a constraint satisfaction problem at all. Don't get me wrong, IPOPT is an amazing software for nonlinear optimization, but when not doing optimization it would lose a lot of what makes it amazing.
However, perfectly agreed that I cannot lay this to rest right now because I don't have a benchmark to point to. The best thing to do here would be to add it to the benchmarks and fully describe why in the paper. We'll definitely follow up with this in a revision.
In the meantime, if you have more requests for the benchmarks, please feel free to open issues at https://github.com/SciML/SciMLBenchmarks.jl so we can track them. All of our benchmarks run on this open platform and anyone can add things via a PR. In particular, the 23 benchmarks diagram in the paper for example is simply just the result of this script https://github.com/SciML/SciMLBenchmarks.jl/blob/master/benc... which runs automatically on any PR (for new folks it requires we click yes for security reasons though). So please feel free to send any benchmark requests. We do plan to do a lot more comprehensive nonlinear optimization benchmarks in the near future. For this case with IPOPT though, we're happy to add that in ourselves hopefully in the next few weeks.
That's not to say that there's no reason to use this trick though. I have seen that some optimization tools can improve robustness in some cases, and this is a way to make use of heuristic methods (genetic algorithms, differential evolution, etc.) for globalizing the process, but it's probably not the first trick you'd want to do in most cases.
While raising derivative order is one way to understand why this should be less efficient, but another is to understand that it loses some specific structure. In the nonlinear system problem f(u)=0, any u that solves this is a solution. All solutions are the same, if you find one you're good. In the optimization problem, you want to find "the" u that minimizes f(u). Now if you are minimizing f(u)=||g(u)|| and you know there is a u that causes g(u) to be zero, then all u that causes g(u) to be zero is a solution, all match the global optima. But solvers don't necessarily know that. Solvers need to use second order information to know that they have not hit a saddle point or other behavior. They don't know that "near zero = done", and they can't specialize on this information. This is another way in which you're just giving sophisticated more solvers more work to do where you already know the answer by problem design.
And, to be clear, there's no guarantee that a zero residual be found, only that `g'(x)*g(x)=0`, but that tends to be the nature of nonlinear equations.
It is not robust, in fact it's well-known as a numerically unstable method. g'(x)*g'(x)` is numerically unstable because it squares the condition number of the matrix. If you have a condition number of 1e-10, which is true for many examples in scientific problems (like the DFN battery model in the paper), the condition number of J'J is 1e-20, which effectively means a linear solve is unable to retrieve any digits of accuracy with 64-bit floating point numbers. So while I agree that you can use symmetric linear solvers as a small performance improvement in some cases, you have to be very careful when you do this as it is a very numerically unstable method. For some details, see https://math.stackexchange.com/a/2874902
Also, if you look at Figure 5 in the paper, you will see that there is a way to choose operator conditioning https://docs.sciml.ai/LinearSolve/stable/basics/OperatorAssu..., and if you tell it to assume the operator is well-conditioned it will actually use this trick. But we do not default to that because of the ill-conditioning.
https://help.juliahub.com/batteries/stable/api/#JuliaSimBatt...
https://docs.sciml.ai/SciMLBenchmarksOutput/dev/NonlinearPro...
But as Avik states, this same kind of thing is done as part of the trust region methods, and so the trust region stuff with some ways to avoid J'J is effectively doing the same thing, though of course that cannot target a symmetric linear solver because J is not generally symmetric and it needs to do the L2 solve so it needs a QR/GMRES.
But indeed, let's get IPOPT in this interface and benchmarks and show it in the plots.
Note that we are preparing an extension to easily use the optimizers https://github.com/SciML/NonlinearSolve.jl/pull/395 to add it to the benchmarks.
Note that we don't however solve the normal form J'J system as that is generally bad (unless user opts in to doing so) and instead use least squares formulation which is more numeraically stable
I forget the exact internal formulation of the problem in IPOPT but for an NLE it’s likely the efficiency of the linear algebra (MA57) algorithm will dominate since it is as you say, constraint satisfaction.
MA57 seems generic but I’ve tried different linear solvers with IPOPT even ones that were supposed to be better, but MA57 still performed the best.
Nonlinear numerical problems involve a bit of art. Theory doesn’t always pan out in practice.
On that remark, a theoretical tool to choose the initial point is the Newton Polygon. For univariate polynomials, this leads to some of the fastest solver [1]. The idea to choose the initial point in this case is notably described in Section 6 of this paper [2].
A very rough intuition of the idea is that a polynomial p(z) doesn't vanish if one of its monomials is significantly larger than all the other. Reciprocally, the polynomial has more chances to vanish for values z such that two monomials have the same order of magnitude. Hence this is a good value to choose your initial point. Finally, deciding if two monomials have the same order of magnitude can be done by taking their log, and then it becomes geometry with polygons, hence the name Newton Polygon.
[1]: https://numpi.dm.unipi.it/scientific-computing-libraries/mps...
[2]: https://woelen.homescience.net/science/math/exps/polynomials...
I think this use case isn't really discussed in the paper because the solver doesn't need to specialize on that case, it should 'just work' without any intervention needed.
The other use case is if you're solving tons of f(u,p)=0 equations for many different p, but all the same `f`. This shows up for example when doing parameter sweeps. You might think that you can just use the former for this: if you originally had `u` a vector, then you can make `u` a matrix and solve `f(u[:,i],p[:,i])` all simultaneously. This is in fact what machine learning libraries tend to do with their `vmap` approach. However, that hides a choice: how do you change `f` to build the composite `f`? You can either do the matrix operations on each operation of `f`, or you can build a specific GPU kernel for `f`. Of course you can do the former, that is exactly the same as the use case above, reusing the same machinery. However, we also offer the latter. And the reason for that is because calling GPU kernels has a lot of overhead and inhibits interprocedural optimizations, and therefore the former is a lot slower than the latter when you can exploit this known repeated structure. We were able to demonstrate this in a recent publication which showed using the latter approach for ODEs was about 20x-100x faster than PyTorch, Jax, and our own implementation of the former, while our approach to the latter is just about as fast as MPGOS which is a CUDA kernel, indicating that it's some limitation of the design (see https://www.sciencedirect.com/science/article/abs/pii/S00457..., or the arxiv version https://arxiv.org/abs/2304.06835). The kernel building approach is thus simply more efficient for things like parameter sweeps, and it's a pretty unique way of using the GPU as most machine learning libraries are not built around kernel generation but instead only give you the linear algebra kernel calling approaches.
tl;dr, DiffEqGPU.jl's documentation is very clear on the two types of GPU usage https://docs.sciml.ai/DiffEqGPU/stable/getting_started/, but is specifically for ODEs. We need to write some similar tutorials for GPU usage on the nonlinear solvers to this effect.