5 ms·

Nonlinearsolve.jl: Fast and Robust Solvers for Nonlinear Equations in Julia

• wenc 3y ago
  Would be interested to see this compared to IPOPT. https://coin-or.github.io/Ipopt/ 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).
  • adgjlsfhk1 3y ago
    it should be pretty easy to add to the benchmarks. From a first glance, it looks like they have a lot of the necessary machinery to do well.
  • jpfr 3y ago
    This article is about “solving” differential equations and not convex optimization.
    • wenc 3y ago
      > This article is about “solving” differential equations and not convex optimization. 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)
      • FreakLegion 3y ago
        > all nonlinear equation systems are nonconvex 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.
        • Bimos 3y ago
          The feasible region {x | f(x) = 0} is nonconvex no matter whether f is convex.
          • graycat 3y ago
            Consider claim: > 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.
          • ubj 3y ago
            I'm assuming you're referring to nonlinear f(x) because this statement is trivially false for linear f(x). 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.
        • wenc 3y ago
          I was talking in terms of convex optimization. The criteria for convex optimization is convex objective, convex inequality constraints and convex feasible region, and linear (not just convex) equality constraints. I’m aware that ipopt is a local solver.
          • FreakLegion 3y ago
            > all nonlinear equation systems are nonconvex > 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.
  • akutlay 3y ago
    I dont know much about this but the fact that they did not mention IPOPT or Highs (which I would think them as two of the most popular software used for these problems) made me question this too.
    • affinepplan 3y ago
      those are indeed very popular tools but for a slightly different kind of problem than this package is attempting to solve
      • ChrisRackauckas 3y ago
        Yes, especially Highs is completely unrelated because it's for (mixed integer) linear programming problems and quadratic programming problems, so it's not able to solve the types of problems described in this manuscript. 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.
        • akutlay 3y ago
          Thank you very much for your answer!
    • avikpal1410 3y ago
      IPOPT as mentioned in other posts solves Nonlinear Optimization Problems. Working around it to make it solve nonlinear equations or nonlinear least squares problems typically doesn't go well. We will add benchmarks for the NLLS part (which we don't talk about in the paper) and that will contain IPOPT comparisons. For some comparison (this is not a benchmark me or anyone on this paper have written), see https://juliapackagecomparisons.github.io/comparisons/math/nonlinear_solvers/index.html#performance_evaluation https://juliapackagecomparisons.github.io/comparisons/math/n.... Specialized solvers for NLLS pretty much always beat general optimization solvers.
  • ChrisRackauckas 3y ago
    Definitely an interesting response I didn't see coming. This is one of the reasons to share preprints! 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 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/benchmarks/NonlinearProblem/nonlinear_solver_23_tests.jmd 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.
    • sobellian 3y ago
      You can also put it as a minimization problem to minimize ||g(u)||, but I assume that does not magically make anything more efficient.
      • 3y ago
• Bimos 3y ago
  On GPU, it can only solve small instances which one single GPU thread can handle, and can only benefit from GPU by solving multiple instances, right?
  • eigenspace 3y ago
    That's the GPU parallelization they offer, but if you had a problem where say the objective function itself was really big and cumbersome and benefited from GPU parallelization, then you could put that objective function on the GPU just fine. 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.
  • ChrisRackauckas 3y ago
    No, there are two types of GPU usage: the ensemble form and the array form. In the array form, you simply make your state variable a GPU and can run all internal operations on the GPU. This is used for GPU acceleration of `f` functions where you have a large state and a relatively "regular" function in `f`, such as neural networks or PDE discretizations. Thus you can see this in action for example on things like the DeepEquilibriumNetworks.jl https://docs.sciml.ai/DeepEquilibriumNetworks/stable/tutorials/basic_mnist_deq/ https://docs.sciml.ai/DeepEquilibriumNetworks/stable/tutoria... use cases, where it's big neural networks with nonlinear solvers in them. 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/S0045782523007156 https://www.sciencedirect.com/science/article/abs/pii/S00457..., or the arxiv version https://arxiv.org/abs/2304.06835 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/ 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.
• deleted 3y ago
  [deleted]