ampl / gsl

GNU Scientific Library with CMake build support and AMPL bindings
https://gsl.ampl.com
GNU General Public License v3.0
555 stars 213 forks source link

Wrong hessian order in amplgsl_hypot3 #76

Closed HunterTracer closed 1 year ago

HunterTracer commented 1 year ago

I use the following code to build a model via pyomo and check hessian via IPOPT. import pyomo.environ as pyo import pyomo.core.base.external as external import pyomo.common.gsl as gsl

dll = gsl.find_GSL() model = pyo.ConcreteModel() model.test_indices = pyo.RangeSet(0, 2) x0 = [-0.2, 0.3, 0.1] model.x = pyo.Var(model.test_indices, domain=pyo.Reals, initialize=lambda m, i: x0[i]) model.gsl_hypot3 = external.ExternalFunction(library=dll, function="gsl_hypot3") model.cons = pyo.Constraint(rule=lambda m: (m.gsl_hypot3(m.x[0], m.x[1], m.x[2]) <= 100)) model.obj = pyo.Objective(rule=lambda m: m.x[0] + m.x[1] + m.x[2], sense=pyo.minimize) optimizer = pyo.SolverFactory("ipopt") optimizer.options["derivative_test"] = "second-order" results = optimizer.solve(model, tee=True)

Then, IPOPT will report warnings like:

Starting derivative checker for second derivatives with obj_factor or lambda[i] set to 1.5.

\* 1-th constr_hess[ 1, 3] = 1.2351616171153007e-01 v ~ -5.4190095373535634e-02 [ 3.279e+00] \* 1-th constr_hess[ 2, 2] = -5.4190094682647579e-02 v ~ 1.2351616144287625e-01 [ 1.439e+00] \* 1-th constr_hess[ 3, 1] = 1.2351616171153007e-01 v ~ -5.4190087606698210e-02 [ 3.279e+00] After checking the source code of amplgsl_hypot3 function in amplgsl.cc, I found that the original hessian order is incorrect, which does not accord with the instruction in funcadd.h: In this case, the function should set al->hes[i + j*(j+1)/2] to the partial derivative of the function with respect to al->ra[i] and al->ra[j] for all 0 <= i <= j < al->nr. In amplgsl_hypot3, the terms in hes[2] and hes[3] should be swaped. I also submit new code at https://github.com/ampl/gsl/pull/75. May the code be merged? Thanks! @4er4er4er @fdabrandao @mapgccv

HunterTracer commented 1 year ago

The fix about the hessian order, which is provided at https://github.com/ampl/gsl/pull/77, has been merged so far. Therefore, the issue has been resolved.