Define MathOptNLSModel when Hessian is not available

[DanielVandH]

The following problem

using JuMP                  
using NLPModelsJuMP         
using JSOSolvers             
model = Model()
g(x...) = 10*length(x) + sum(x[i].^2 .- 10*cos.(2π*x[i]) for i in 1:length(x))
register(model, :g, 5, g, autodiff = true)
x₀ = [1.0, 0.1, 0.2, -0.5, 1.0]
@variable(model, x[i=1:5], start = x₀[i])
@NLexpression(model, res[i in 1:5], g(x...))
nls = MathOptNLSModel(model, res, name = "NL")

returns an error LoadError: Unsupported feature Hess. Providing a Hessian through register does not work as JuMP does not support Hessians for multivariate functions. I've seen in #68 that this issue was solved for MathOptNLPModel.

Has this been considered for MathOptNLSModel?

[amontoison]

Hi @DanielVandH !
Why do you need to register g?
The following model works:

using JuMP                  
using NLPModelsJuMP
using JSOSolvers
model = Model()
n = 5
x₀ = [1.0, 0.1, 0.2, -0.5, 1.0]
@variable(model, x[i=1:5], start = x₀[i])
@NLexpression(model, res[i=1:5], 10*n + sum(x[j]^2 - 10*cos(2π*x[j]) for j=1:n))
nls = MathOptNLSModel(model, res, name = "NL")

Is it what you want ?

[DanielVandH]

Hi @amontoison,

Sorry, I should clarify that this is just a much simpler example than what I'm really trying to do where registering functions really is necessary. My real problem involves five other functions and a much more complicated expression that changes based on what functions a user provides, but this simpler case is representative of all the parts needed.

[amontoison]

I work on a fix @DanielVandH (#95)
If you have time, can you try my branch hessian-nls

pkg> add https://github.com/amontoison/NLPModelsJuMP.jl.git#hessian-nls

with your more complex nls model ?

[DanielVandH]

Thanks @amontoison . I've tried this out, and it does work for the construction of the model. When trying to solve it, I get the following error. Is this a related problem?

julia> output = trunk(nls)
ERROR: LoadError: BoundsError: attempt to access 0-element Vector{Float64} at index [1:3]
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] forward_eval(storage::Vector{Float64}, partials_storage::Vector{Float64}, nd::Vector{JuMP._Derivatives.NodeData}, adj::SparseArrays.SparseMatrixCSC{Bool, Int64}, const_values::Vector{Float64}, parameter_values::Vector{Float64}, x_values::Vector{Float64}, subexpression_values::Vector{Float64}, user_input_buffer::Vector{Float64}, user_output_buffer::Vector{Float64}, user_operators::JuMP._Derivatives.UserOperatorRegistry)
    @ JuMP._Derivatives C:\Users\licer\.julia\packages\JuMP\lnUbA\src\_Derivatives\forward.jl:165
  [5] _forward_eval_all(d::NLPEvaluator, x::Vector{Float64})
    @ JuMP C:\Users\licer\.julia\packages\JuMP\lnUbA\src\nlp.jl:709
  [6] macro expansion
    @ C:\Users\licer\.julia\packages\JuMP\lnUbA\src\nlp.jl:801 [inlined]
  [7] macro expansion
    @ .\timing.jl:299 [inlined]
  [8] eval_constraint(d::NLPEvaluator, g::SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false}, x::Vector{Float64})
    @ JuMP C:\Users\licer\.julia\packages\JuMP\lnUbA\src\nlp.jl:799
  [9] residual!(nls::MathOptNLSModel, x::Vector{Float64}, Fx::Vector{Float64})
    @ NLPModelsJuMP C:\Users\licer\.julia\packages\NLPModelsJuMP\eKJIH\src\moi_nls_model.jl:96
 [10] residual(nls::MathOptNLSModel, x::Vector{Float64})
    @ NLPModels C:\Users\licer\.julia\packages\NLPModels\s4zXc\src\nls\api.jl:16
 [11] trunk(::Val{:GaussNewton}, nlp::MathOptNLSModel; x::Vector{Float64}, subsolver::Symbol, atol::Float64, rtol::Float64, max_eval::Int64, max_time::Float64, bk_max::Int64, monotone::Bool, nm_itmax::Int64, trsolver_args::Dict{Symbol, Any})
    @ JSOSolvers C:\Users\licer\.julia\packages\JSOSolvers\NLU9e\src\trunkls.jl:57
 [12] trunk
    @ C:\Users\licer\.julia\packages\JSOSolvers\NLU9e\src\trunkls.jl:35 [inlined]
 [13] trunk(nlp::MathOptNLSModel; variant::Symbol, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ JSOSolvers C:\Users\licer\.julia\packages\JSOSolvers\NLU9e\src\trunkls.jl:3
 [14] trunk(nlp::MathOptNLSModel)
    @ JSOSolvers C:\Users\licer\.julia\packages\JSOSolvers\NLU9e\src\trunkls.jl:3
 [15] top-level scope
    @ c:\Users\licer\TestIssue.jl:48
in expression starting at c:\Users\licer\TestIssue.jl:48

Here's the code I used, it's essentially what I'm trying to do (with some generated data for the purpose of this issue).

## Define functions
using Pkg, JuMP
#]add https://github.com/amontoison/NLPModelsJuMP.jl.git#hessian-nls
function T(α...)
    t = α[3]
    1 / (1 + exp(-α[1] - α[2] * t))
end
function D(β...)
    u = β[end]
    β[1] + β[2] * u
end
function D′(β...)
    u = β[end]
    β[2]
end
function R(γ...)
    u = γ[end]
    γ[1] * u + γ[2] * u^2
end
## Setup the model 
model = Model()
α₀ = [-3.0, 0.2]
β₀ = [300.0, 300.0]
γ₀ = [0.05, 0.1]
tt = length(α₀)
d = length(β₀)
r = length(γ₀)
register(model, :T, tt + 1, T, autodiff = true)
register(model, :D, d + 1, D, autodiff = true)
register(model, :D′, d + 1, D′, autodiff = true)
register(model, :R, r + 1, R, autodiff = true)
@variable(model, α[i = 1:tt], start = α₀[i])
@variable(model, β[i = 1:d], start = β₀[i])
@variable(model, γ[i = 1:r], start = γ₀[i])
## Generating some data for issue 
using Random
Random.seed!(123)
fₜ = abs.(randn(100))
t = 1:100
f = abs.(randn(100))
fₓ = abs.(randn(100))
fₓₓ = abs.(randn(100))
inIdx = 1:100
## Define and optimise
@NLexpression(model, res[i in inIdx], f[i] - T(α..., t[i]) * (D′(β..., f[i]) + D(β..., f[i]) * fₓₓ[i] * fₜ[i] + R(γ..., f[i])))
nls = MathOptNLSModel(model, res, name = "NonlinearEQL", hessian = false)
using JSOSolvers
output = trunk(nls)

[amontoison]

@DanielVandH
I finally found the culprit !
An auxiliary vector user_output_buffer is required for user defined functions and must be manually copied (https://github.com/JuliaSmoothOptimizers/NLPModelsJuMP.jl/pull/95/files#diff-47c27891e951c8cd946b850dc2df31082624afdf57446c21cb6992f5f4b74aa2R404).
The model of your previous message works fine on my computer.

Can you try again please with my updated fork?

[DanielVandH]

@amontoison Awesome. Works great! I modified my test function a bit and also tested to see if it works on problems with bounds or nonlinear constraints. It doesn't work on problems with nonlinear constraints, giving the following issue with CaNNOLeS:

ERROR: LoadError: Hessian computations were not requested on the call to initialize!.
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:33
 [2] hessian_lagrangian_structure(d::NLPEvaluator)
   @ JuMP C:\Users\licer\.julia\packages\JuMP\lnUbA\src\nlp.jl:1359
 [3] hess_structure_residual!(nls::MathOptNLSModel, rows::Vector{Int64}, cols::Vector{Int64})
   @ NLPModelsJuMP C:\Users\licer\.julia\packages\NLPModelsJuMP\zt7JO\src\moi_nls_model.jl:210
 [4] hess_structure_residual
   @ C:\Users\licer\.julia\packages\NLPModels\s4zXc\src\nls\api.jl:286 [inlined]
 [5] cannoles(nls::MathOptNLSModel; x::Vector{Float64}, λ::Vector{Float64}, method::Symbol, merit::Symbol, linsolve::Symbol, max_f::Int64, max_time::Float64, max_inner::Int64, ϵtol::Float64, check_small_residual::Bool, always_accept_extrapolation::Bool, ϵkchoice::Symbol, δdec::Float64)
   @ CaNNOLeS C:\Users\licer\.julia\packages\CaNNOLeS\CGoHr\src\CaNNOLeS.jl:97
 [6] cannoles(nls::MathOptNLSModel)
   @ CaNNOLeS C:\Users\licer\.julia\packages\CaNNOLeS\CGoHr\src\CaNNOLeS.jl:50
 [7] top-level scope
   @ c:\Users\licer\TestIssue.jl:73
in expression starting at c:\Users\licer\TestIssue.jl:73

This could be a separate problem, though. Maybe I'm reading the CaNNOLeS code wrong but it seems purely second-order so might not be a way around it. Here is my modified script for testing.

## Define functions
using Pkg, JuMP
#]add https://github.com/amontoison/NLPModelsJuMP.jl.git#hessian-nls
function T(α...)
    t = α[3]
    1 / (1 + exp(-α[1] - α[2] * t))
end
function D(β...)
    u = β[end]
    β[1] + β[2] * u
end
function D′(β...)
    u = β[end]
    β[2]
end
function R(γ...)
    u = γ[end]
    γ[1] * u + γ[2] * u^2
end
## Setup the model 
model = Model()
α₀ = [-3.0, 0.2]
β₀ = [5.0, 5.0]
γ₀ = [0.05, 0.1]
tt = length(α₀)
d = length(β₀)
r = length(γ₀)
register(model, :T, tt + 1, T, autodiff = true)
register(model, :D, d + 1, D, autodiff = true)
register(model, :D′, d + 1, D′, autodiff = true)
register(model, :R, r + 1, R, autodiff = true)
@variable(model, α[i = 1:tt], start = α₀[i])
@variable(model, β[i = 1:d], start = β₀[i])
@variable(model, γ[i = 1:r], start = γ₀[i])
## Generating some data for issue 
using Random
Random.seed!(123)
fₜ = abs.(randn(100))
t = 1:100
f = abs.(randn(100))
fₓ = abs.(randn(100))
fₓₓ = abs.(randn(100))
inIdx = 1:100
## Define and optimise
@NLexpression(model, res[i in inIdx], f[i] - T(α..., t[i]) * (D′(β..., fₓ[i]) + D(β..., fₓ[i]) * fₓₓ[i] * fₜ[i] + R(γ..., fₜ[i])))
nls = MathOptNLSModel(model, res, name = "NonlinearEQL", hessian = false)
using JSOSolvers

# Is the fix solver dependent?
output = trunk(nls)
output2 = tron(nls)

# How does it compare to an independent solver?
using LsqFit 
xdata = hcat(t, fₜ, fₓ, fₓₓ)
model_inplace(f, x, θ) = (@. f = (T(θ[1:2]..., xdata[:, 1]) * (D′(θ[3:4]..., xdata[:, 3]) + D(θ[3:4]..., xdata[:, 3]) * xdata[:, 2] * xdata[:, 4] + R(θ[5:6]..., xdata[:, 4]))))
output3 = curve_fit(model_inplace, xdata, f, [α₀;β₀;γ₀]; inplace = true)

soln1 = output.solution 
soln2 = output2.solution 
soln3 = output3.param # Looks good 

# Does it change based on bounds? 
set_lower_bound(γ[1], 0.0)
nls = MathOptNLSModel(model, res, name = "Bounded", hessian = false)
soln4 = tron(nls).solution # Good 

# Try a nonlinear constraint 
using CaNNOLeS
@NLconstraint(model, α[1]^2 + α[2]^2 == 80.0)
delete_lower_bound(model[:γ][1])
nls = MathOptNLSModel(model, res, name = "NonlinearConstraint ", hessian = false)
soln5 = cannoles(nls).solution

[amontoison]

@DanielVandH, CaNNOLeS requires the hessian :(

@tmigot, @abelsiqueira
Do we have a solver for the previous problem that do not use the hessian ?
Do we have a ADNLSModel to avoid the JuMP / MathOptNLSModel model ?

[tmigot]

Thanks @DanielVandH for the issue!

One way I see is you can penalize the constraint in the JuMP model and use trunk, however, there are some drawbacks.

As @amontoison mentioned, you can use ADNLSModel, see docstring in ADNLPModels.jl for your problem and then use CaNNOLeS. Here is a tutorial for unconstrained problems juliasmoothoptimizers.github.io/CaNNOLeS.jl/dev/tutorial/.
The drawback of this approach right now is that we don't have the sparse hessian and jacobian just yet (but that's a work in progress).