symbolify models cause BoundsError

[odow]

julia> import Nonconvex

julia> import NonconvexUtils

julia> Nonconvex.@load Ipopt
[ Info: Attempting to load the package NonconvexIpopt.
[ Info: Loading succesful.

julia> model = Nonconvex.Model(x -> x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]);

julia> addvar!(model, fill(1.0, 4), fill(5.0, 4));

julia> add_ineq_constraint!(model, x -> 25.0 - x[1] * x[2] * x[3] * x[4]);

julia> add_eq_constraint!(model, x -> 40.0 - x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2);

julia> sym_model = NonconvexUtils.symbolify(
           model;
           hessian = true,
           sparse = true,
           simplify = true,
       );

julia> result = Nonconvex.optimize(
           sym_model,
           IpoptAlg(),
           [1.0, 5.0, 5.0, 1.0];
           options = IpoptOptions(; first_order = false, sparse = true),
       )
This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.

Number of nonzeros in equality constraint Jacobian...:        4
Number of nonzeros in inequality constraint Jacobian.:        4
Number of nonzeros in Lagrangian Hessian.............:        9

Total number of variables............................:        4
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        4
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        1

ERROR: BoundsError: attempt to access 9-element Vector{Float64} at index [10:10]
Stacktrace:
  [1] throw_boundserror(A::Vector{Float64}, I::Tuple{UnitRange{Int64}})
    @ Base ./abstractarray.jl:691
  [2] checkbounds
    @ ./abstractarray.jl:656 [inlined]
  [3] view
    @ ./subarray.jl:177 [inlined]
  [4] maybeview
    @ ./views.jl:146 [inlined]
  [5] dotview
    @ ./broadcast.jl:1200 [inlined]
  [6] add_values!(values::Vector{Float64}, HL::LowerTriangular{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}}; factor::Int64, offset::Int64)
    @ NonconvexCore ~/.julia/packages/NonconvexCore/nK7uX/src/models/moi.jl:299
  [7] add_values!(values::Vector{Float64}, HL::LowerTriangular{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}})
    @ NonconvexCore ~/.julia/packages/NonconvexCore/nK7uX/src/models/moi.jl:295
  [8] (::NonconvexIpopt.var"#8#16"{Bool, NonconvexIpopt.var"#lag#10"{NonconvexCore.CountingFunction{NonconvexCore.Objective{NonconvexCore.var"#144#151"{Model{Vector{Any}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Base.RefValue{Float64}, Set{Symbol}}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.IneqConstraint{NonconvexCore.var"#148#155"{NonconvexCore.IneqConstraint{NonconvexUtils.var"#99#100"{NonconvexUtils.SymbolicFunction{NonconvexUtils.var"#5#6"{NonconvexCore.IneqConstraint{FunctionWrapper{var"#299#300", Int64}, Float64, Int64, Set{Symbol}}, NonconvexCore.Unflatten{Tuple{Vector{Any}}, NonconvexCore.var"#unflatten_to_Tuple#70"{Tuple{Int64}, Tuple{Int64}, Tuple{NonconvexCore.var"#Vector_from_vec#62"{Vector{Any}, Vector{NonconvexCore.var"#unflatten_to_Real#58"}, Vector{Vector{Float64}}}}}}}, NonconvexUtils.var"#66#74"{NonconvexUtils.var"#64#72"{DataType}}, NonconvexUtils.var"#68#76"{NonconvexUtils.var"#64#72"{DataType}}, Vector{Float64}}, typeof(identity)}, Float64, Int64, Set{Symbol}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Float64, Int64, Set{Symbol}}}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.EqConstraint{NonconvexCore.var"#146#153"{NonconvexCore.EqConstraint{NonconvexUtils.var"#99#100"{NonconvexUtils.SymbolicFunction{NonconvexUtils.var"#5#6"{NonconvexCore.EqConstraint{FunctionWrapper{var"#301#302", Int64}, Float64, Int64, Set{Symbol}}, NonconvexCore.Unflatten{Tuple{Vector{Any}}, NonconvexCore.var"#unflatten_to_Tuple#70"{Tuple{Int64}, Tuple{Int64}, Tuple{NonconvexCore.var"#Vector_from_vec#62"{Vector{Any}, Vector{NonconvexCore.var"#unflatten_to_Real#58"}, Vector{Vector{Float64}}}}}}}, NonconvexUtils.var"#66#74"{NonconvexUtils.var"#64#72"{DataType}}, NonconvexUtils.var"#68#76"{NonconvexUtils.var"#64#72"{DataType}}, Vector{Float64}}, typeof(identity)}, Float64, Int64, Set{Symbol}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Float64, Int64, Set{Symbol}}}}}})(x::Vector{Float64}, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
    @ NonconvexIpopt ~/.julia/packages/NonconvexIpopt/zl6RS/src/ipopt.jl:211
  [9] _Eval_H_CB(n::Int32, x_ptr::Ptr{Float64}, #unused#::Int32, obj_factor::Float64, m::Int32, lambda_ptr::Ptr{Float64}, #unused#::Int32, nele_hess::Int32, iRow::Ptr{Int32}, jCol::Ptr{Int32}, values_ptr::Ptr{Float64}, user_data::Ptr{Nothing})
    @ Ipopt ~/.julia/packages/Ipopt/M2QE8/src/C_wrapper.jl:122
 [10] IpoptSolve(prob::Ipopt.IpoptProblem)
    @ Ipopt ~/.julia/packages/Ipopt/M2QE8/src/C_wrapper.jl:433
 [11] optimize!(workspace::NonconvexIpopt.IpoptWorkspace{NonconvexCore.VecModel{NonconvexCore.Objective{NonconvexCore.var"#144#151"{Model{Vector{Any}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Base.RefValue{Float64}, Set{Symbol}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.EqConstraint{NonconvexCore.var"#146#153"{NonconvexCore.EqConstraint{NonconvexUtils.var"#99#100"{NonconvexUtils.SymbolicFunction{NonconvexUtils.var"#5#6"{NonconvexCore.EqConstraint{FunctionWrapper{var"#301#302", Int64}, Float64, Int64, Set{Symbol}}, NonconvexCore.Unflatten{Tuple{Vector{Any}}, NonconvexCore.var"#unflatten_to_Tuple#70"{Tuple{Int64}, Tuple{Int64}, Tuple{NonconvexCore.var"#Vector_from_vec#62"{Vector{Any}, Vector{NonconvexCore.var"#unflatten_to_Real#58"}, Vector{Vector{Float64}}}}}}}, NonconvexUtils.var"#66#74"{NonconvexUtils.var"#64#72"{DataType}}, NonconvexUtils.var"#68#76"{NonconvexUtils.var"#64#72"{DataType}}, Vector{Float64}}, typeof(identity)}, Float64, Int64, Set{Symbol}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Float64, Int64, Set{Symbol}}}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.IneqConstraint{NonconvexCore.var"#148#155"{NonconvexCore.IneqConstraint{NonconvexUtils.var"#99#100"{NonconvexUtils.SymbolicFunction{NonconvexUtils.var"#5#6"{NonconvexCore.IneqConstraint{FunctionWrapper{var"#299#300", Int64}, Float64, Int64, Set{Symbol}}, NonconvexCore.Unflatten{Tuple{Vector{Any}}, NonconvexCore.var"#unflatten_to_Tuple#70"{Tuple{Int64}, Tuple{Int64}, Tuple{NonconvexCore.var"#Vector_from_vec#62"{Vector{Any}, Vector{NonconvexCore.var"#unflatten_to_Real#58"}, Vector{Vector{Float64}}}}}}}, NonconvexUtils.var"#66#74"{NonconvexUtils.var"#64#72"{DataType}}, NonconvexUtils.var"#68#76"{NonconvexUtils.var"#64#72"{DataType}}, Vector{Float64}}, typeof(identity)}, Float64, Int64, Set{Symbol}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Float64, Int64, Set{Symbol}}}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.SDConstraint}}, Vector{Float64}, Vector{Float64}, Vector{Float64}}, Ipopt.IpoptProblem, Vector{Float64}, IpoptOptions{NamedTuple{(:hessian_approximation, :jac_c_constant, :jac_d_constant, :sparse), Tuple{String, String, String, Bool}}}, Base.RefValue{Int64}})
    @ NonconvexIpopt ~/.julia/packages/NonconvexIpopt/zl6RS/src/ipopt.jl:56
 [12] #optimize#142
    @ ~/.julia/packages/NonconvexCore/nK7uX/src/models/vec_model.jl:63 [inlined]
 [13] optimize(::Model{Vector{Any}}, ::IpoptAlg, ::Vector{Float64}; kwargs::Base.Pairs{Symbol, IpoptOptions{NamedTuple{(:hessian_approximation, :jac_c_constant, :jac_d_constant, :sparse), Tuple{String, String, String, Bool}}}, Tuple{Symbol}, NamedTuple{(:options,), Tuple{IpoptOptions{NamedTuple{(:hessian_approximation, :jac_c_constant, :jac_d_constant, :sparse), Tuple{String, String, String, Bool}}}}}})
    @ NonconvexCore ~/.julia/packages/NonconvexCore/nK7uX/src/common.jl:233
 [14] top-level scope
    @ REPL[103]:1

[mohamed82008]

This particular example is fixed in the latest release JuliaRegistries/General#62556. However, a bigger problem remains.

The problem is that sometimes the initial Hessian has structural non-zeros that have the value of 0. The way I get Hessians is using ForwardDiff over Zygote. ForwardDiff causes a loss in sparsity in the Hessian in the process. So I workaround this issue by using sparse on the Hessian result. This is problematic because a non-structural zeros with value 0 get assumed to be structural zeros. So if the initial Hessian gets mistaken for being more sparse than it really is, the Hessians in the later iterations would cause the above error to be thrown.

The reason why we can't use the sparsity pattern returned by Symbolics here is that it is hidden in the Lagrangian function which calls the constraint/objective functions which call the symbolified functions. So the SymbolicFunction type is actually not directly accessible from the Lagrangian anonymous function defined in the Ipopt wrapper.

There are 2 ways to workaround this:

1. Make ForwardDiff correctly return a sparse Jacobian (requires more ForwardDiff hacking). Nothing I tried seems to work. I may need to do more Dual number hacking to get this working. I successfully hacked around Zygote to make it return a sparse Jacobian.
2. Use Zygote over Zygote but that currently runs into Zygote higher order differentiation bugs.
3. Use sparsify/symbolify directly on the Lagrangian function in NonconvexIpopt and use the sparse Hessian function defined in the returned types. This should have the correct sparsity pattern that Symbolics figures out and it is possible now but it requires having NonconvexIpopt depend on NonconvexUtils which is mostly unnecessary.

I might have to go with the third option. For now, you can just make sure the initial Hessian is non-zero for any structural non-zeros at the initial solution. I worked around the individual example in the issue reported by ensuring the diagonal of the Hessian is always treated as structural non-zeros. The 0 value in that example happened to be on the diagonal so that solved the problem in this case.

[odow]

For now, you can just make sure the initial Hessian is non-zero for any structural non-zeros at the initial solution

How would one do this?

[mohamed82008]

By using a random initial solution that's unlikely to lead to term cancellation.

[ccoffrin]

By using a random initial solution that's unlikely to lead to term cancellation.

I am not sure I follow what this means, but I will remark that the starting point for your NLP solver is quite important to convergence and solution quality, users will need the ability to set specific starting points that are connected to domain knowledge about the problem that is being solved.

[ccoffrin]

Based on my testing for lanl-ansi/rosetta-opf#23 I do believe that the selected starting point impacts the correctness of the symbolic AD system.

I would recommend that this behavior be avoided if possible as it will be quite unexpected to users of Nonconvex, who have experience using other NLP modeling frameworks.

[mohamed82008]

NonconvexIpopt 0.4 should have fixed this issue.

[mohamed82008]

I will go ahead and close this issue since I don't think the issue exists anymore with NonconvexIpopt 0.4. If you think otherwise, please report with a new MWE.

[mohamed82008]

For the record, this is the new, simpler working code with NonconvexIpopt 0.4:

import Nonconvex
import NonconvexUtils
Nonconvex.@load Ipopt

model = Nonconvex.Model(x -> x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]);
addvar!(model, fill(1.0, 4), fill(5.0, 4));
add_ineq_constraint!(model, x -> 25.0 - x[1] * x[2] * x[3] * x[4]);
add_eq_constraint!(model, x -> 40.0 - x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2);
result = Nonconvex.optimize(
   model,
   IpoptAlg(),
   [1.0, 5.0, 5.0, 1.0];
   options = IpoptOptions(; first_order = false, symbolic = true, sparse = true),
);

[ccoffrin]

Based on my testing for lanl-ansi/rosetta-opf#23 I do believe that the selected starting point impacts the correctness of the symbolic AD system.

I would recommend that this behavior be avoided if possible as it will be quite unexpected to users of Nonconvex, who have experience using other NLP modeling frameworks.

Ok, I'll open a new issue for this point, with an example on how to reproduce.

[ccoffrin]

You might want to re-open this issue. This appears to happen when testing the nonconvex model in rosetta-opf on the pglib_opf_case14_ieee data file from this repo, https://github.com/power-grid-lib/pglib-opf

[odow]

using a random initial solution that's unlikely to lead to term cancellation

Presumably Nonconvex could do this internally? Why do you need to use the user's primal starting point? Even then, it's still a probabilistic approach. It wouldn't work for something like f(x) = ifelse(x <= 5, 0, (x-5)^2).

[mohamed82008]

Symbolic differentiation will not work for any function that's not Symbolics-compatible. A function with if statements (I believe ifesle is an exception) may violate the assumptions in Symbolics. I believe the issue with the reported example has now been resolved so I don't think it falls under the Symbolics-incompatible class of functions.

[ccoffrin]

The pglib_opf_case14_ieee appears to be working now, thanks!