Add four functions in the NLP API (#2162)

Author: amontoison

docs/src/reference/nonlinear.md

@@ -32,10 +32,14 @@ eval_objective
 eval_constraint
 eval_objective_gradient
 jacobian_structure
-hessian_lagrangian_structure
 eval_constraint_jacobian
 eval_constraint_jacobian_product
 eval_constraint_jacobian_transpose_product
+hessian_lagrangian_structure
+hessian_objective_structure
+hessian_constraint_structure
+eval_hessian_objective
+eval_hessian_constraint
 eval_hessian_lagrangian
 eval_hessian_lagrangian_product
 objective_expr

src/Nonlinear/ReverseAD/mathoptinterface_api.jl

@@ -156,7 +156,16 @@ function MOI.initialize(d::NLPEvaluator, requested_features::Vector{Symbol})
         end
         d.max_chunk = max_chunk
         if d.want_hess
-            _compute_hessian_lagrangian_structure(d)
+            d.hessian_sparsity = Tuple{Int64,Int64}[]
+            if d.objective !== nothing
+                append!(
+                    d.hessian_sparsity,
+                    zip(d.objective.hess_I, d.objective.hess_J),
+                )
+            end
+            for c in d.constraints
+                append!(d.hessian_sparsity, zip(c.hess_I, c.hess_J))
+            end
         end
     end
     return
@@ -250,22 +259,41 @@ function MOI.eval_constraint_jacobian_transpose_product(
     return
 end
 
+function MOI.hessian_objective_structure(d::NLPEvaluator)
+    @assert d.want_hess
+    obj = d.objective
+    return Tuple{Int64,Int64}[(i, j) for (i, j) in zip(obj.hess_I, obj.hess_J)]
+end
+
+function MOI.hessian_constraint_structure(d::NLPEvaluator, c::Integer)
+    @assert d.want_hess
+    con = d.constraints[c]
+    return Tuple{Int64,Int64}[(i, j) for (i, j) in zip(con.hess_I, con.hess_J)]
+end
+
 function MOI.hessian_lagrangian_structure(d::NLPEvaluator)
     @assert d.want_hess
     return d.hessian_sparsity
 end
 
-function _compute_hessian_lagrangian_structure(d::NLPEvaluator)
-    d.hessian_sparsity = Tuple{Int64,Int64}[]
+function MOI.eval_hessian_objective(d::NLPEvaluator, H, x)
+    @assert d.want_hess
+    _reverse_mode(d, x)
+    fill!(d.input_ϵ, 0.0)
     if d.objective !== nothing
-        append!(d.hessian_sparsity, zip(d.objective.hess_I, d.objective.hess_J))
-    end
-    for c in d.constraints
-        append!(d.hessian_sparsity, zip(c.hess_I, c.hess_J))
+        _eval_hessian(d, d.objective, H, 1.0, 0)
     end
     return
 end
 
+function MOI.eval_hessian_constraint(d::NLPEvaluator, H, x, i)
+    @assert d.want_hess
+    _reverse_mode(d, x)
+    fill!(d.input_ϵ, 0.0)
+    _eval_hessian(d, d.constraints[i], H, 1.0, 0)
+    return
+end
+
 function MOI.eval_hessian_lagrangian(d::NLPEvaluator, H, x, σ, μ)
     @assert d.want_hess
     _reverse_mode(d, x)

src/Nonlinear/ReverseAD/precompile.jl

@@ -6,6 +6,23 @@
 
 function _precompile_()
     ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
+    Base.precompile(
+        Tuple{
+            typeof(MOI.eval_hessian_objective),
+            NLPEvaluator,
+            SubArray{Float64,1,Vector{Float64},Tuple{UnitRange{Int64}},true},
+            Vector{Float64},
+        },
+    )
+    Base.precompile(
+        Tuple{
+            typeof(MOI.eval_hessian_constraint),
+            NLPEvaluator,
+            SubArray{Float64,1,Vector{Float64},Tuple{UnitRange{Int64}},true},
+            Vector{Float64},
+            Int64,
+        },
+    )
     Base.precompile(
         Tuple{
             typeof(MOI.eval_hessian_lagrangian),
@@ -42,6 +59,13 @@ function _precompile_()
             Vector{Float64},
         },
     )
+    Base.precompile(Tuple{typeof(MOI.hessian_objective_structure),NLPEvaluator})
+    Base.precompile(
+        Tuple{typeof(MOI.hessian_constraint_structure),NLPEvaluator,Int64},
+    )
+    Base.precompile(
+        Tuple{typeof(MOI.hessian_lagrangian_structure),NLPEvaluator},
+    )
     Base.precompile(Tuple{typeof(MOI.jacobian_structure),NLPEvaluator})
     Base.precompile(
         Tuple{typeof(MOI.eval_objective),NLPEvaluator,Vector{Float64}},

src/Nonlinear/evaluator.jl

@@ -69,6 +69,8 @@ function MOI.initialize(evaluator::Evaluator, features::Vector{Symbol})
     evaluator.eval_objective_gradient_timer = 0.0
     evaluator.eval_constraint_timer = 0.0
     evaluator.eval_constraint_jacobian_timer = 0.0
+    evaluator.eval_hessian_objective_timer = 0.0
+    evaluator.eval_hessian_constraint_timer = 0.0
     evaluator.eval_hessian_lagrangian_timer = 0.0
     append!(evaluator.ordered_constraints, keys(evaluator.model.constraints))
     filter!(f -> f != :ExprGraph, features)
@@ -143,10 +145,32 @@ function MOI.eval_constraint_jacobian(evaluator::Evaluator, J, x)
     return
 end
 
+function MOI.hessian_objective_structure(evaluator::Evaluator)
+    return MOI.hessian_objective_structure(evaluator.backend)
+end
+
+function MOI.hessian_constraint_structure(evaluator::Evaluator, i)
+    return MOI.hessian_constraint_structure(evaluator.backend, i)
+end
+
 function MOI.hessian_lagrangian_structure(evaluator::Evaluator)
     return MOI.hessian_lagrangian_structure(evaluator.backend)
 end
 
+function MOI.eval_hessian_objective(evaluator::Evaluator, H, x)
+    start = time()
+    MOI.eval_hessian_objective(evaluator.backend, H, x)
+    evaluator.eval_hessian_objective_timer += time() - start
+    return
+end
+
+function MOI.eval_hessian_constraint(evaluator::Evaluator, H, x, i)
+    start = time()
+    MOI.eval_hessian_constraint(evaluator.backend, H, x, i)
+    evaluator.eval_hessian_constraint_timer += time() - start
+    return
+end
+
 function MOI.eval_hessian_lagrangian(evaluator::Evaluator, H, x, σ, μ)
     start = time()
     MOI.eval_hessian_lagrangian(evaluator.backend, H, x, σ, μ)

src/Nonlinear/types.jl

@@ -209,6 +209,8 @@ mutable struct Evaluator{B} <: MOI.AbstractNLPEvaluator
     eval_constraint_timer::Float64
     eval_objective_gradient_timer::Float64
     eval_constraint_jacobian_timer::Float64
+    eval_hessian_objective_timer::Float64
+    eval_hessian_constraint_timer::Float64
     eval_hessian_lagrangian_timer::Float64
 
     function Evaluator(
@@ -225,6 +227,8 @@ mutable struct Evaluator{B} <: MOI.AbstractNLPEvaluator
             0.0,
             0.0,
             0.0,
+            0.0,
+            0.0,
         )
     end
 end

src/nlp.jl

@@ -203,6 +203,47 @@ The sparsity structure is assumed to be independent of the point ``x``.
 """
 function jacobian_structure end
 
+"""
+    hessian_objective_structure(
+        d::AbstractNLPEvaluator,
+    )::Vector{Tuple{Int64,Int64}}
+
+Returns a vector of tuples, `(row, column)`, where each indicates the position
+of a structurally nonzero element in the Hessian matrix:
+``\\nabla^2 f(x)``.
+
+The indices are not required to be sorted and can contain duplicates, in which
+case the solver should combine the corresponding elements by adding them
+together.
+
+Any mix of lower and upper-triangular indices is valid. Elements `(i,j)` and
+`(j,i)`, if both present, should be treated as duplicates.
+
+The sparsity structure is assumed to be independent of the point ``x``.
+"""
+function hessian_objective_structure end
+
+"""
+    hessian_constraint_structure(
+        d::AbstractNLPEvaluator,
+        i::Int64,
+    )::Vector{Tuple{Int64,Int64}}
+
+Returns a vector of tuples, `(row, column)`, where each indicates the position
+of a structurally nonzero element in the Hessian matrix:
+``nabla^2 g_i(x)``.
+
+The indices are not required to be sorted and can contain duplicates, in which
+case the solver should combine the corresponding elements by adding them
+together.
+
+Any mix of lower and upper-triangular indices is valid. Elements `(i,j)` and
+`(j,i)`, if both present, should be treated as duplicates.
+
+The sparsity structure is assumed to be independent of the point ``x``.
+"""
+function hessian_constraint_structure end
+
 """
     hessian_lagrangian_structure(
         d::AbstractNLPEvaluator,
@@ -312,6 +353,47 @@ For example, it may be the `view` of a vector.
 """
 function eval_hessian_lagrangian_product end
 
+"""
+    eval_hessian_objective(
+        d::AbstractNLPEvaluator,
+        H::AbstractVector{T},
+        x::AbstractVector{T},
+    )::Nothing where {T}
+
+This function computes the sparse Hessian matrix:
+``\\nabla^2 f(x)``,
+storing the result in the vector `H` in the same order as the indices
+returned by [`hessian_objective_structure`](@ref).
+
+## Implementation notes
+
+When implementing this method, you must not assume that `H` is
+`Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
+For example, it may be the `view` of a vector.
+"""
+function eval_hessian_objective end
+
+"""
+    eval_hessian_constraint(
+        d::AbstractNLPEvaluator,
+        H::AbstractVector{T},
+        x::AbstractVector{T},
+        i::Int64,
+    )::Nothing where {T}
+
+This function computes the sparse Hessian matrix:
+``\\nabla^2 g_i(x)``,
+storing the result in the vector `H` in the same order as the indices
+returned by [`hessian_constraint_structure`](@ref).
+
+## Implementation notes
+
+When implementing this method, you must not assume that `H` is
+`Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
+For example, it may be the `view` of a vector.
+"""
+function eval_hessian_constraint end
+
 """
     eval_hessian_lagrangian(
         d::AbstractNLPEvaluator,

test/Nonlinear/ReverseAD.jl

@@ -36,6 +36,10 @@ function test_objective_quadratic_univariate()
     g = [NaN]
     MOI.eval_objective_gradient(evaluator, g, [1.2])
     @test g == [2.4]
+    @test MOI.hessian_objective_structure(evaluator) == [(1, 1)]
+    H = [NaN]
+    MOI.eval_hessian_objective(evaluator, H, [1.2])
+    @test H == [2.0]
     @test MOI.hessian_lagrangian_structure(evaluator) == [(1, 1)]
     H = [NaN]
     MOI.eval_hessian_lagrangian(evaluator, H, [1.2], 1.5, Float64[])
@@ -63,6 +67,14 @@ function test_objective_and_constraints_quadratic_univariate()
     g = [NaN]
     MOI.eval_objective_gradient(evaluator, g, [1.2])
     @test g == [2.4]
+    @test MOI.hessian_objective_structure(evaluator) == [(1, 1)]
+    H = [NaN]
+    MOI.eval_hessian_objective(evaluator, H, [1.2])
+    @test H == [2.0]
+    @test MOI.hessian_constraint_structure(evaluator, 1) == [(1, 1)]
+    H = [NaN]
+    MOI.eval_hessian_constraint(evaluator, H, [1.3], 1)
+    @test H == [2.0]
     @test MOI.hessian_lagrangian_structure(evaluator) == [(1, 1), (1, 1)]
     H = [NaN, NaN]
     MOI.eval_hessian_lagrangian(evaluator, H, [1.2], 1.5, Float64[1.3])
@@ -92,6 +104,10 @@ function test_objective_quadratic_multivariate()
     g = [NaN, NaN]
     MOI.eval_objective_gradient(evaluator, g, [1.2, 2.3])
     @test g == [2 * 1.2 + 2.3, 1.2 + 2 * 2.3]
+    @test MOI.hessian_objective_structure(evaluator) == [(1, 1), (2, 2), (2, 1)]
+    H = [NaN, NaN, NaN]
+    MOI.eval_hessian_objective(evaluator, H, [1.2, 2.3])
+    @test H == [2.0, 2.0, 1.0]
     @test MOI.hessian_lagrangian_structure(evaluator) ==
           [(1, 1), (2, 2), (2, 1)]
     H = [NaN, NaN, NaN]
@@ -126,6 +142,10 @@ function test_objective_quadratic_multivariate_subexpressions()
     g = [NaN, NaN]
     MOI.eval_objective_gradient(evaluator, g, [1.2, 2.3])
     @test g == [2 * 1.2 + 2.3, 1.2 + 2 * 2.3]
+    @test MOI.hessian_objective_structure(evaluator) == [(1, 1), (2, 2), (2, 1)]
+    H = [NaN, NaN, NaN]
+    MOI.eval_hessian_objective(evaluator, H, [1.2, 2.3])
+    @test H == [2.0, 2.0, 1.0]
     @test MOI.hessian_lagrangian_structure(evaluator) ==
           [(1, 1), (2, 2), (2, 1)]
     H = [NaN, NaN, NaN]
@@ -225,6 +245,10 @@ function test_constraint_quadratic_univariate()
     J = [NaN]
     MOI.eval_constraint_jacobian(evaluator, J, x_val)
     @test J == 2 .* x_val
+    @test MOI.hessian_constraint_structure(evaluator, 1) == [(1, 1)]
+    H = [NaN]
+    MOI.eval_hessian_constraint(evaluator, H, x_val, 1)
+    @test H == [2.0]
     @test MOI.hessian_lagrangian_structure(evaluator) == [(1, 1)]
     H = [NaN]
     MOI.eval_hessian_lagrangian(evaluator, H, x_val, 0.0, [1.5])