[MOI] Add support for nonlinear problems without Hessian (#322)

ext/MadNLPMOI/MadNLPMOI.jl

@@ -99,6 +99,7 @@ MOI.eval_constraint(::_EmptyNLPEvaluator, g, x) = nothing
 MOI.jacobian_structure(::_EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
 MOI.hessian_lagrangian_structure(::_EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
 MOI.eval_constraint_jacobian(::_EmptyNLPEvaluator, J, x) = nothing
+MOI.eval_constraint_jacobian_transpose_product(::_EmptyNLPEvaluator, Jtv, x, v) = nothing
 MOI.eval_hessian_lagrangian(::_EmptyNLPEvaluator, H, x, σ, μ) = nothing
 
 function MOI.empty!(model::Optimizer)
@@ -123,6 +124,9 @@ function MOI.empty!(model::Optimizer)
     if haskey(model.options, :jacobian_constant)
         delete!(model.options, :jacobian_constant)
     end
+    if haskey(model.options, :hessian_approximation)
+        delete!(model.options, :hessian_approximation)
+    end
     return
 end
 
@@ -141,6 +145,16 @@ function MOI.copy_to(model::Optimizer, src::MOI.ModelLike)
     return MOI.Utilities.default_copy_to(model, src)
 end
 
+function _init_nlp_model(model)
+    if model.nlp_model === nothing
+        if !(model.nlp_data.evaluator isa _EmptyNLPEvaluator)
+            error("Cannot mix the new and legacy nonlinear APIs")
+        end
+        model.nlp_model = MOI.Nonlinear.Model()
+    end
+    return
+end
+
 MOI.get(::Optimizer, ::MOI.SolverName) = "MadNLP"
 
 function MOI.supports_add_constrained_variable(
@@ -487,6 +501,28 @@ function MOI.get(model::Optimizer, attr::MOI.ListOfSupportedNonlinearOperators)
     return MOI.get(model.nlp_model, attr)
 end
 
+### UserDefinedFunction
+
+MOI.supports(model::Optimizer, ::MOI.UserDefinedFunction) = true
+
+function MOI.set(model::Optimizer, attr::MOI.UserDefinedFunction, args)
+    _init_nlp_model(model)
+    MOI.Nonlinear.register_operator(
+        model.nlp_model,
+        attr.name,
+        attr.arity,
+        args...,
+    )
+    return
+end
+
+### ListOfSupportedNonlinearOperators
+
+function MOI.get(model::Optimizer, attr::MOI.ListOfSupportedNonlinearOperators)
+    _init_nlp_model(model)
+    return MOI.get(model.nlp_model, attr)
+end
+
 ### MOI.VariablePrimalStart
 
 function MOI.supports(
@@ -726,6 +762,13 @@ function MOI.eval_constraint_jacobian(model::Optimizer, values, x)
     return
 end
 
+function MOI.eval_constraint_jacobian_transpose_product(model::Optimizer, Jtv, x, v)
+    MOI.eval_constraint_jacobian_transpose_product(model.nlp_data.evaluator, Jtv, x, v)
+    # Evaluate QPBlockData after NLPEvaluator to ensure that Jtv is not reset.
+    MOI.eval_constraint_jacobian_transpose_product(model.qp_data, Jtv, x, v)
+    return
+end
+
 ### Eval_H_CB
 
 function MOI.hessian_lagrangian_structure(model::Optimizer)
@@ -751,23 +794,27 @@ end
 
 NLPModels.obj(nlp::MOIModel, x::AbstractVector{Float64}) = MOI.eval_objective(nlp.model,x)
 
-function NLPModels. grad!(nlp::MOIModel, x::AbstractVector{Float64}, g::AbstractVector{Float64})
+function NLPModels.grad!(nlp::MOIModel, x::AbstractVector{Float64}, g::AbstractVector{Float64})
     MOI.eval_objective_gradient(nlp.model, g, x)
 end
 
-function NLPModels. cons!(nlp::MOIModel, x::AbstractVector{Float64}, c::AbstractVector{Float64})
+function NLPModels.cons!(nlp::MOIModel, x::AbstractVector{Float64}, c::AbstractVector{Float64})
     MOI.eval_constraint(nlp.model, c, x)
 end
 
-function NLPModels. jac_coord!(nlp::MOIModel, x::AbstractVector{Float64}, jac::AbstractVector{Float64})
+function NLPModels.jac_coord!(nlp::MOIModel, x::AbstractVector{Float64}, jac::AbstractVector{Float64})
     MOI.eval_constraint_jacobian(nlp.model, jac, x)
 end
 
-function NLPModels. hess_coord!(nlp::MOIModel, x::AbstractVector{Float64}, l::AbstractVector{Float64}, hess::AbstractVector{Float64}; obj_weight::Float64=1.0)
+function NLPModels.jtprod!(nlp::MOIModel, x::AbstractVector{Float64}, v::Vector{Float64}, Jtv::AbstractVector{Float64})
+    MOI.eval_constraint_jacobian_transpose_product(nlp.model, Jtv, x, v)
+end
+
+function NLPModels.hess_coord!(nlp::MOIModel, x::AbstractVector{Float64}, l::AbstractVector{Float64}, hess::AbstractVector{Float64}; obj_weight::Float64=1.0)
     MOI.eval_hessian_lagrangian(nlp.model, hess, x, obj_weight, l)
 end
 
-function NLPModels. hess_structure!(nlp::MOIModel, I::AbstractVector{T}, J::AbstractVector{T}) where T
+function NLPModels.hess_structure!(nlp::MOIModel, I::AbstractVector{T}, J::AbstractVector{T}) where T
     @assert length(I) == length(J) == length(MOI.hessian_lagrangian_structure(nlp.model))
     cnt = 1
     for (row, col) in  MOI.hessian_lagrangian_structure(nlp.model)
@@ -804,11 +851,11 @@ function MOIModel(model::Optimizer)
         any(isequal(_kFunctionTypeScalarQuadratic), model.qp_data.function_type)
     has_nlp_constraints = !isempty(model.nlp_data.constraint_bounds)
     has_nlp_objective = model.nlp_data.has_objective
-    has_nlp_hessian = :Hess in MOI.features_available(model.nlp_data.evaluator)
+    has_hessian = :Hess in MOI.features_available(model.nlp_data.evaluator)
     is_nlp = has_nlp_constraints || has_nlp_objective
     # Initialize evaluator using model's structure.
     init_feat = [:Grad]
-    if has_nlp_hessian
+    if has_hessian
         push!(init_feat, :Hess)
     end
     if has_nlp_constraints
@@ -836,7 +883,11 @@ function MOIModel(model::Optimizer)
 
     # Sparsity
     jacobian_sparsity = MOI.jacobian_structure(model)
-    hessian_sparsity = MOI.hessian_lagrangian_structure(model)
+    hessian_sparsity = if has_hessian
+        MOI.hessian_lagrangian_structure(model)
+    else
+        Tuple{Int,Int}[]
+    end
     nnzh = length(hessian_sparsity)
     nnzj = length(jacobian_sparsity)
 
@@ -858,6 +909,9 @@ function MOIModel(model::Optimizer)
     if !has_nlp_constraints && !has_quadratic_constraints
         model.options[:jacobian_constant] = true
     end
+    if !has_hessian
+        model.options[:hessian_approximation] = MadNLP.CompactLBFGS
+    end
 
     return MOIModel(
         NLPModels.NLPModelMeta(
@@ -954,7 +1008,7 @@ function MOI.get(model::Optimizer, ::MOI.RawStatusString)
     elseif model.solver === nothing
         return "Optimize not called"
     end
-    return MadNLP.get_status_output(model.result.status, model.result.options) 
+    return MadNLP.get_status_output(model.result.status, model.result.options)
 end
 
 

ext/MadNLPMOI/utils.jl

@@ -201,24 +201,25 @@ function eval_sparse_gradient(
     f::MOI.ScalarQuadraticFunction{T},
     x::Vector{T},
     p::Dict{Int64,T},
+    adj::T,
 )::Int where {T}
     i = 0
     for term in f.affine_terms
         if !_is_parameter(term.variable)
             i += 1
-            ∇f[i] = term.coefficient
+            ∇f[i] += term.coefficient * adj
         end
     end
     for term in f.quadratic_terms
         if !_is_parameter(term.variable_1)
             v = _value(term.variable_2, x, p)
             i += 1
-            ∇f[i] = term.coefficient * v
+            ∇f[i] += term.coefficient * v * adj
         end
         if term.variable_1 != term.variable_2 && !_is_parameter(term.variable_2)
             v = _value(term.variable_1, x, p)
             i += 1
-            ∇f[i] = term.coefficient * v
+            ∇f[i] += term.coefficient * v * adj
         end
     end
     return i
@@ -229,12 +230,13 @@ function eval_sparse_gradient(
     f::MOI.ScalarAffineFunction{T},
     x::Vector{T},
     p::Dict{Int64,T},
+    adj::T,
 )::Int where {T}
     i = 0
     for term in f.terms
         if !_is_parameter(term.variable)
             i += 1
-            ∇f[i] = term.coefficient
+            ∇f[i] += term.coefficient * adj
         end
     end
     return i
@@ -508,13 +510,26 @@ function MOI.eval_constraint_jacobian(
     x::AbstractVector{T},
 ) where {T}
     i = 1
+    fill!(J, zero(T))
     for constraint in block.constraints
         ∇f = view(J, i:length(J))
-        i += eval_sparse_gradient(∇f, constraint, x, block.parameters)
+        i += eval_sparse_gradient(∇f, constraint, x, block.parameters, one(T))
     end
     return i
 end
 
+function MOI.eval_constraint_jacobian_transpose_product(
+    block::QPBlockData{T},
+    Jtv::AbstractVector{T},
+    x::AbstractVector{T},
+    v::AbstractVector{T},
+) where {T}
+    for (i, constraint) in enumerate(block.constraints)
+        eval_sparse_gradient(∇f, constraint, x, block.parameters, v[i])
+    end
+    return Jtv
+end
+
 function MOI.hessian_lagrangian_structure(block::QPBlockData)
     H = Tuple{Int,Int}[]
     append_sparse_hessian_structure!(block.objective, H)

lib/MadNLPTests/src/Instances/hs15nohessian.jl

@@ -0,0 +1,65 @@
+struct HS15NoHessianModel{T} <: NLPModels.AbstractNLPModel{T,Vector{T}}
+    meta::NLPModels.NLPModelMeta{T, Vector{T}}
+    counters::NLPModels.Counters
+end
+
+function HS15NoHessianModel(;T = Float64, x0=zeros(T,2), y0=zeros(T,2))
+    return HS15NoHessianModel(
+        NLPModels.NLPModelMeta(
+            2,     #nvar
+            ncon = 2,
+            nnzj = 4,
+            nnzh = 0,
+            x0 = x0,
+            y0 = y0,
+            lvar = T[-Inf, -Inf],
+            uvar = T[0.5, Inf],
+            lcon = T[1.0, 0.0],
+            ucon = T[Inf, Inf],
+            minimize = true
+        ),
+        NLPModels.Counters()
+    )
+end
+
+function NLPModels.obj(nlp::HS15NoHessianModel, x::AbstractVector)
+    return 100.0 * (x[2] - x[1]^2)^2 + (1.0 - x[1])^2
+end
+
+function NLPModels.grad!(nlp::HS15NoHessianModel, x::AbstractVector, g::AbstractVector)
+    z = x[2] - x[1]^2
+    g[1] = -400.0 * z * x[1] - 2.0 * (1.0 - x[1])
+    g[2] = 200.0 * z
+    return
+end
+
+function NLPModels.cons!(nlp::HS15NoHessianModel, x::AbstractVector, c::AbstractVector)
+    c[1] = x[1] * x[2]
+    c[2] = x[1] + x[2]^2
+end
+
+function NLPModels.jac_structure!(nlp::HS15NoHessianModel, I::AbstractVector{T}, J::AbstractVector{T}) where T
+    copyto!(I, [1, 1, 2, 2])
+    copyto!(J, [1, 2, 1, 2])
+end
+
+function NLPModels.jac_coord!(nlp::HS15NoHessianModel, x::AbstractVector, J::AbstractVector)
+    J[1] = x[2]    # (1, 1)
+    J[2] = x[1]    # (1, 2)
+    J[3] = 1.0     # (2, 1)
+    J[4] = 2*x[2]  # (2, 2)
+    return J
+end
+
+function NLPModels.jprod!(nlp::HS15NoHessianModel, x::AbstractVector, v::AbstractVector, jv::AbstractVector)
+    jv[1] = x[2] * v[1] + x[1] * v[2]
+    jv[2] = v[1] + 2 * x[2] * v[2]
+    return jv
+end
+
+function NLPModels.jtprod!(nlp::HS15NoHessianModel, x::AbstractVector, v::AbstractVector, jv::AbstractVector)
+    jv[1] = x[2] * v[1] + v[2]
+    jv[2] = x[1] * v[1] + 2 * x[2] * v[2]
+    return jv
+end
+

lib/MadNLPTests/src/MadNLPTests.jl

@@ -354,6 +354,7 @@ end
 
 include("Instances/dummy_qp.jl")
 include("Instances/hs15.jl")
+include("Instances/hs15nohessian.jl")
 include("Instances/nls.jl")
 include("wrapper.jl")
 

src/nlpmodels.jl

@@ -357,7 +357,9 @@ function create_callback(
     jac_scale = similar(jac_buffer, nnzj) ; fill!(jac_scale, one(T))
 
     NLPModels.jac_structure!(nlp, jac_I, jac_J)
-    NLPModels.hess_structure!(nlp, hess_I, hess_J)
+    if nnzh > 0
+        NLPModels.hess_structure!(nlp, hess_I, hess_J)
+    end
 
     fixed_handler, nnzj, nnzh = create_sparse_fixed_handler(
         fixed_variable_treatment,

test/MOI_interface_test.jl

@@ -50,8 +50,6 @@ function test_MOI_Test()
             # - Excluded because Hessian information is needed
             "test_nonlinear_hs071_hessian_vector_product",
             # - Excluded because Hessian information is needed
-            "test_nonlinear_hs071_no_hessian",
-            # - Excluded because Hessian information is needed
             "test_nonlinear_invalid",
 
             #  - Excluded because this test is optional
@@ -99,6 +97,20 @@ function test_invalid_number_in_hessian_lagrangian()
     return
 end
 
+# See issue #318
+function test_user_defined_function()
+    model = MadNLP.Optimizer()
+    MOI.set(model, MOI.Silent(), true)
+    # Define custom function.
+    f(a, b) = a^2 + b^2
+    x = MOI.add_variables(model, 2)
+    MOI.set(model, MOI.UserDefinedFunction(:f, 2), (f,))
+    obj_f = MOI.ScalarNonlinearFunction(:f, Any[x[1], x[2]])
+    MOI.set(model, MOI.ObjectiveFunction{typeof(obj_f)}(), obj_f)
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED
+end
+
 end
 
 TestMOIWrapper.runtests()

test/madnlp_quasi_newton.jl

@@ -34,7 +34,7 @@ end
 
 @testset "MadNLP: LBFGS" begin
     @testset "HS15" begin
-        nlp = MadNLPTests.HS15Model()
+        nlp = MadNLPTests.HS15NoHessianModel()
         solver_qn = MadNLP.MadNLPSolver(
             nlp;
             callback = MadNLP.SparseCallback,