Why using JuMP or AMPL gives different number of nonzeros in Lagrangian Hessian?

slmontes · January 27, 2025, 12:25pm

Hello,

I was comparing NLP optimization between JuMP and AMPL and noticed that even when the problems are identical and the presolve option is turned off in AMPL, the number of nonzeros in Lagrangian Hessian varies.

This also occurs when I use amplNLWriter instead of using JuMP alone.

Does anyone know why this might be happening?

Here is a simple example using the clnlbeam problem from JuMP.jl/stable/tutorials/nonlinear/simple_examples/

In JuMP:

using JuMP
using Ipopt
using Test

function example_clnlbeam()
 N = 1000
 h = 1 / N
 alpha = 350
 model = Model(Ipopt.Optimizer)
 @variables(model, begin
 -1 <= t[1:(N+1)] <= 1
 -0.05 <= x[1:(N+1)] <= 0.05
 u[1:(N+1)]
    end)
 @objective(
 model,
 Min,
 sum(
            0.5 * h * (u[i+1]^2 + u[i]^2) +
            0.5 * alpha * h * (cos(t[i+1]) + cos(t[i])) for i in 1:N
 ),
 )
 @constraint(
 model,
 [i = 1:N],
 x[i+1] - x[i] - 0.5 * h * (sin(t[i+1]) + sin(t[i])) == 0,
 )
 @constraint(
 model,
 [i = 1:N],
 t[i+1] - t[i] - 0.5 * h * u[i+1] - 0.5 * h * u[i] == 0,
 )
 optimize!(model)
 println("""
 termination_status = $(termination_status(model))
 primal_status      = $(primal_status(model))
 objective_value    = $(objective_value(model))
 """)
 Test.@test is_solved_and_feasible(model)
    return
end

example_clnlbeam()

Output:

This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

Number of nonzeros in equality constraint Jacobian...:     8000
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     4002

Total number of variables............................:     3003
                     variables with only lower bounds:        0
                variables with lower and upper bounds:     2002
                     variables with only upper bounds:        0
Total number of equality constraints.................:     2000
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  3.5000000e+02 0.00e+00 0.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  3.5000000e+02 0.00e+00 0.00e+00  -1.7 0.00e+00    -  1.00e+00 1.00e+00   0
   2  3.5000000e+02 0.00e+00 0.00e+00  -3.8 0.00e+00  -2.0 1.00e+00 1.00e+00T  0
   3  3.5000000e+02 0.00e+00 0.00e+00  -5.7 0.00e+00   0.2 1.00e+00 1.00e+00T  0
   4  3.5000000e+02 0.00e+00 0.00e+00  -8.6 0.00e+00  -0.2 1.00e+00 1.00e+00T  0

Number of Iterations....: 4

                                   (scaled)                 (unscaled)
Objective...............:   3.5000000000000318e+02    3.5000000000000318e+02
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5059035596802450e-09    2.5059035596802450e-09
Overall NLP error.......:   2.5059035596802450e-09    2.5059035596802450e-09


Number of objective function evaluations             = 5
Number of objective gradient evaluations             = 5
Number of equality constraint evaluations            = 5
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 5
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 4
Total seconds in IPOPT                               = 0.037

EXIT: Optimal Solution Found.
termination_status = LOCALLY_SOLVED
primal_status      = FEASIBLE_POINT
objective_value    = 350.0000000000032

Using AmplNLWriter:

using JuMP
using AmplNLWriter
using Ipopt
using Test

model = Model(() -> AmplNLWriter.Optimizer("ipopt"))

function example_clnlbeam()
 N = 1000
 h = 1 / N
 alpha = 350
 @variables(model, begin
 -1 <= t[1:(N+1)] <= 1
 -0.05 <= x[1:(N+1)] <= 0.05
 u[1:(N+1)]
    end)
 @objective(
 model,
 Min,
 sum(
            0.5 * h * (u[i+1]^2 + u[i]^2) +
            0.5 * alpha * h * (cos(t[i+1]) + cos(t[i])) for i in 1:N
 ),
 )
 @constraint(
 model,
 [i = 1:N],
 x[i+1] - x[i] - 0.5 * h * (sin(t[i+1]) + sin(t[i])) == 0,
 )
 @constraint(
 model,
 [i = 1:N],
 t[i+1] - t[i] - 0.5 * h * u[i+1] - 0.5 * h * u[i] == 0,
 )
 optimize!(model)
 println("""
 termination_status = $(termination_status(model))
 primal_status      = $(primal_status(model))
 objective_value    = $(objective_value(model))
 """)
 Test.@test is_solved_and_feasible(model)
    return
end

example_clnlbeam()

Output:

This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

Number of nonzeros in equality constraint Jacobian...:     8000
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     2002

Total number of variables............................:     3003
                     variables with only lower bounds:        0
                variables with lower and upper bounds:     2002
                     variables with only upper bounds:        0
Total number of equality constraints.................:     2000
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  3.5000000e+02 0.00e+00 0.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  3.5000000e+02 0.00e+00 0.00e+00  -1.7 0.00e+00    -  1.00e+00 1.00e+00   0
   2  3.5000000e+02 0.00e+00 0.00e+00  -3.8 0.00e+00  -2.0 1.00e+00 1.00e+00T  0
   3  3.5000000e+02 0.00e+00 0.00e+00  -5.7 0.00e+00   0.2 1.00e+00 1.00e+00T  0
   4  3.5000000e+02 0.00e+00 0.00e+00  -8.6 0.00e+00  -0.2 1.00e+00 1.00e+00T  0

Number of Iterations....: 4

                                   (scaled)                 (unscaled)
Objective...............:   3.5000000000000318e+02    3.5000000000000318e+02
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5059035596802450e-09    2.5059035596802450e-09
Overall NLP error.......:   2.5059035596802450e-09    2.5059035596802450e-09


Number of objective function evaluations             = 5
Number of objective gradient evaluations             = 5
Number of equality constraint evaluations            = 5
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 5
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 4
Total seconds in IPOPT                               = 0.049

EXIT: Optimal Solution Found.
termination_status = LOCALLY_SOLVED
primal_status      = FEASIBLE_POINT
objective_value    = 350.00000000000074

gdalle · January 27, 2025, 2:36pm

Hi @slmontes,
Since the number of zeros for JuMP is twice that of AMPL, could it be that AMPL is counting only the upper or lower triangle of the Hessian?
It may also be due to JuMP’s sparsity detection algorithm, which is sometimes a bit conservative, see e.g. Nonzero elements in the Hessian that should be zero · Issue #1811 · jump-dev/JuMP.jl · GitHub

odow · January 27, 2025, 9:18pm

Hi @slmontes, welcome to the forum

You can safely ignore the differences in the non-zero count between Ipopt.jl and AmplNLWriter.jl. As you can see from the logs, Ipopt takes identical sequences of iterations.

The differences are due to two reasons:

Multiple non-zeros can be specified for the same element in the Hessian. These are summed together and have no effect on the output.
JuMP’s sparsity detection is conservative. It may have “non-zeros” in the Hessian which are in fact always 0.0.

In some rare cases there might be a different sequence of iterations due to small numerical round off differences in how the Hessian matrix is evaluated.

odow · January 27, 2025, 9:20pm

It may also be due to JuMP’s sparsity detection algorithm, which is sometimes a bit conservative, see e.g. Nonzero elements in the Hessian that should be zero · Issue #1811 · jump-dev/JuMP.jl · GitHub

I think you meant instead [Nonlinear] sparsity pattern of Hessian with :(x * y) · Issue #2527 · jump-dev/MathOptInterface.jl · GitHub

slmontes · January 28, 2025, 1:39pm

Thank you! This is helpful to know. I also thought it may have been due to using symmetry assumptions for simplification. However, it does not always yield double the number of non-zeros, so this makes sense.

Topic		Replies	Views
Too many nonzeroes for NLP in JuMP and Examodels Optimization (Mathematical) jump	8	98	December 18, 2024
Error in JuMP with nonlinear objective Optimization (Mathematical)	6	563	February 28, 2020
Debugging in Ipopt without JuMP Optimization (Mathematical) ipopt	5	939	June 29, 2022
How to use non-linear solvers with JuMP using AmplNLWriter Optimization (Mathematical) jump , optimization	2	238	February 8, 2024
Difference in sparsity pattern caused by @NLconstraint and @constraint in JuMP Optimization (Mathematical)	4	1299	July 30, 2019

Why using JuMP or AMPL gives different number of nonzeros in Lagrangian Hessian?

Related topics