4 ms·

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

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.
- ChrisRackauckas 3y ago
  That makes your Newton method require calculating the Hessian of the objective function, which is a second derivative. A normal Newton method for nonlinear systems uses the Jacobian, which is the first derivative (the connection is that the Hessian is the Jacobian of the gradient). But indeed, formulating it so that you have discontinuities (absolute values) and require second derivatives instead of the first for the same algorithm is not the good kind of magic :-/. You could then try using only gradient-based methods, like Adam, but you can see that is actually using less problem structure. BFGS for example is related to Broyden's method and is thus using simple Hessian approximations, so for nonlinear system solving this is effectively less stable than just using Newton. 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.
  - kxyvr 3y ago
    This is not entirely true. A typical trick for accomplishing this is to simply use the Hessian approximation `g'(x)*g'(x)` (star denotes adjoint), thus cutting off the second-derivative evaluation of g, and then use a modified CG method like Newton-Krylov for line-search methods or Steihaug-Toint for trust-region methods. It's very robust, matrix-free, and doesn't require an excessive amount of code to implement. Further, the method can be preconditioned with a positive definite preconditioner. 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.
wenc 3y ago
Yes tribal knowledge says NLP solvers aren’t tailored to NLE solvers as NLE solvers are. There have been papers in this past showing this but in my experience, it’s not obvious in practice. Nonlinear systems can often surprise us. The initial point (and randomness) has more influence on the solution time than almost anything else. 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.
- adgjlsfhk1 3y ago
  isn't MA57 just a standard good sparse linear solver? For sparse problems, NonlinearSolve.jl (via LinearSolve.jl) will either use KLU or Cholmod both of which are also quite fast. Also, even if the linear solver is the bottleneck, that doesn't mean different algorithms won't make a differenc.e Different nonlinear solve algorithms can lead to needing fewer iterations (e.g. fewer linear solves).
  - wenc 3y ago
    I don’t know. I’m curious is all. 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.
    - kxyvr 3y ago
      MA57 is a sparse symmetric solver. IPOPT uses this because it solves a modified KKT system that includes both the Hessian and constraint information simultaneously. This differs from composite step methods, which split the optimization step into a step for feasibility (quasi-normal step) and a step for optimality (tangential step). An algorithm like NITRO uses a composite step method and as a result uses a different kind of linear system solver.
- guyomes 3y ago
  > The initial point (and randomness) has more influence on the solution time than almost anything else. 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/mpsolve/ https://numpi.dm.unipi.it/scientific-computing-libraries/mps... [2]: https://woelen.homescience.net/science/math/exps/polynomials/mpsolve3_theory_paper.pdf https://woelen.homescience.net/science/math/exps/polynomials...

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.

ChrisRackauckas 3y ago
That makes your Newton method require calculating the Hessian of the objective function, which is a second derivative. A normal Newton method for nonlinear systems uses the Jacobian, which is the first derivative (the connection is that the Hessian is the Jacobian of the gradient). But indeed, formulating it so that you have discontinuities (absolute values) and require second derivatives instead of the first for the same algorithm is not the good kind of magic :-/. You could then try using only gradient-based methods, like Adam, but you can see that is actually using less problem structure. BFGS for example is related to Broyden's method and is thus using simple Hessian approximations, so for nonlinear system solving this is effectively less stable than just using Newton. 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.
- kxyvr 3y ago
  This is not entirely true. A typical trick for accomplishing this is to simply use the Hessian approximation `g'(x)*g'(x)` (star denotes adjoint), thus cutting off the second-derivative evaluation of g, and then use a modified CG method like Newton-Krylov for line-search methods or Steihaug-Toint for trust-region methods. It's very robust, matrix-free, and doesn't require an excessive amount of code to implement. Further, the method can be preconditioned with a positive definite preconditioner. 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.

wenc 3y ago

Yes tribal knowledge says NLP solvers aren’t tailored to NLE solvers as NLE solvers are. There have been papers in this past showing this but in my experience, it’s not obvious in practice. Nonlinear systems can often surprise us. The initial point (and randomness) has more influence on the solution time than almost anything else. 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.

adgjlsfhk1 3y ago
isn't MA57 just a standard good sparse linear solver? For sparse problems, NonlinearSolve.jl (via LinearSolve.jl) will either use KLU or Cholmod both of which are also quite fast. Also, even if the linear solver is the bottleneck, that doesn't mean different algorithms won't make a differenc.e Different nonlinear solve algorithms can lead to needing fewer iterations (e.g. fewer linear solves).
- wenc 3y ago
  I don’t know. I’m curious is all. 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.
  - kxyvr 3y ago
    MA57 is a sparse symmetric solver. IPOPT uses this because it solves a modified KKT system that includes both the Hessian and constraint information simultaneously. This differs from composite step methods, which split the optimization step into a step for feasibility (quasi-normal step) and a step for optimality (tangential step). An algorithm like NITRO uses a composite step method and as a result uses a different kind of linear system solver.
guyomes 3y ago
> The initial point (and randomness) has more influence on the solution time than almost anything else. 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/mpsolve/ https://numpi.dm.unipi.it/scientific-computing-libraries/mps... [2]: https://woelen.homescience.net/science/math/exps/polynomials/mpsolve3_theory_paper.pdf https://woelen.homescience.net/science/math/exps/polynomials...