objcons with AD

[tpapp]

Consider the problem

$$\max_x g(f(x))$$

subject to

$$h(f(x)) = 0$$

and

$$l \le x \le u$$

where $x \in \mathbb{R}^n$ is a vector, $f: \mathbb{R}^n \to \mathbb{R}^m$ is expensive to calculate compared to $g$ and $h$, $g: \mathbb{R}^m \to \mathbb{R}$ and $h: \mathbb{R}^m \to \mathbb{R}^k$.

The natural solution would be using objcons, but I need help figuring out how to integrate AD with that.

My (admittedly very superficial) understanding of the API is that it allows value & gradient/jacobian for the objective and the constraint separately, but not together. But perhaps I missed something.

For the purposes of an MWE, a walk-through using

x0 = zeros(3)
f(x) = vcat(abs2.(x), x)        # we just pretend this is expensive
g(y) = sum(@view y[1:3])        # pretent problem is not separable
h(y) = y[4:6] .- 1              # still pretending this is not trivial
# optimum is ones(3)

would help me a lot. Remember, the exercise is to evaluate f once for each x, and use AD (ForwardDiff is fine).

[abelsiqueira]

Since you have a custom function, there is no easy way to do it, so you have to define an AbstractNLPModel and implement some of its values.
Something like

using NLPModels

mutable struct YourProblem{T, S} <: AbstractNLPModel{T, S}
  meta::NLPModelMeta{T, S}
  counters::Counters
end

function YourProblem(::Type{T}) where {T}
  meta = NLPModelMeta{T, Vector{T}}(
    2, # number of variables,
    ncon = 2, # number of constraints,
    x0 = ...
    lvar = ...
    uvar = ...
    lcon = ...
    ucon = ...
    name = ...
  )
end
YourProblem() = YourProblem(Float64)

function NLPModels.objcons(nlp::YourProblem, x::AbstractVector)
  @lencheck 2 x
  increment!(nlp, :neval_obj)
  increment!(nlp, :neval_cons)

  fx = ...
  obj_value = g(fx)
  cons_value = h(fx) 
 
  return obj_value, cons_value
end

You will need to define other functions as well, and you need to handle AD yourself.
Which functions you'll need to implement depend on the solver.

Since the model is completely new, you can do other things that you might deem necessary, for instance, storing fx inside the struct itself. It might be useful to pass it around functions.

There are a few examples to get inspiration, such as ManualNLPModels.jl, or hs10.jl from NLPModelsTest.jl.
There is a Guidelines for creating NLPModels, but some things might be outdated (hopefully not).

Let me know if this is helpful.

[tpapp]

Thanks for your answer. The only part unclear is the AD. You say that

You will need to define other functions as well, and you need to handle AD yourself.

but I could not figure out which part of the API handles the objective, its gradient, the constraint, and the Jacobian at once. Is there even such a thing?

If that is not available, I could do the following: for each x, whenever either objcons, objgrad, jac_op, etc are called, evaluate $f$, $\nabla f$, $c$, $J$, using AD, and save it (only keeping the last $n$, how large should $n$ be with linearsearch etc?) Then I could just look these up and serve it in the relevant functions.

In any case, I am still in the planning phase, comments are appreciated.

[abelsiqueira]

Hi again,

Thanks for your answer. The only part unclear is the AD. You say that

You will need to define other functions as well, and you need to handle AD yourself.

but I could not figure out which part of the API handles the objective, its gradient, the constraint, and the Jacobian at once. Is there even such a thing?

No.

If that is not available, I could do the following: for each x, whenever either objcons, objgrad, jac_op, etc are called, evaluate f, ∇f, c, J, using AD, and save it (only keeping the last n, how large should n be with linearsearch etc?) Then I could just look these up and serve it in the relevant functions.

In any case, I am still in the planning phase, comments are appreciated.

This is the way. The functions that your model needs to implement are the ones called here: https://github.com/JuliaSmoothOptimizers/Percival.jl/blob/main/src/AugLagModel.jl

For instance, update_cx! calls cons!, so you need to implement cons!. etc.
I haven't really checked, but I think the list is obj, grad!, cons!, objcons!, jprod!, jtprod! and hprod!.

The internal steps are done by tron from JSOSolvers.jl. I don't have a good recollection of when it needs to re-evaluate any of the functions, maybe @dpo or @tmigot know better? My initial guess would be the maximum is given by max_cgiter, which defaults to 50.

[tmigot]

Hi @tpapp ! I hope @abelsiqueira answer "solved" your issue. It seems the right choice to store the computations to be reused.
I did something similar here https://github.com/JuliaSmoothOptimizers/FletcherPenaltySolver.jl/blob/57980538c5581264e323a16f17038ad0540a11e0/src/model-Fletcherpenaltynlp.jl#L234 in another solver where I compute:

shahx = hash(x)

and if shahx is different than the stored hash, then we update.

[tpapp]

Forgot to mention: I am probably going to be using Percival.jl, I find it very robust.

[tpapp]

I started an implementation at https://github.com/tpapp/ObjConsNLPModels.jl.

But I am running into a dimension error with Percival:

julia> using ObjConsNLPModels

julia> using Percival

julia> model = objcons_nlpmodel(x -> [sum(abs2, x), sum(x) - 1]; x0 = [2.0, 2.0])
ObjConsNLPModels.ObjConsNLPModel{Float64, Vector{Float64}, Float64, var"#24#25"}
  Problem name: Generic
                  All variables: ████████████████████ 2                     All constraints: ████████████████████ 1                                free: ████████████████████ 2                                free: ████████████████████ 1                               lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                               lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                               upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                               upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                             low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                             low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                               fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                               fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                              infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                              infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 nnzh: (  0.00% sparsity)   3                              linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
                                                                                  nonlinear: ████████████████████ 1                                                                             nnzj: (  0.00% sparsity)   2     

  Counters:
                            obj: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                                grad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                                cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
                       cons_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                            cons_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                                jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
                          jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                                 jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                             jac_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
                        jac_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                               jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                           jprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
                      jprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                              jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                          jtprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
                     jtprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                                hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                               hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
                          jhess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                              jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     

julia> output = percival(model)
ERROR: DimensionError: Input lvar should have length 3 not 2

What am I doing wrong?

[tmigot]

Hi @tpapp !

In your example, since there is no lcon and ucon in the NLPModelMeta, it means that you are modeling the constraint:

-Inf <= c(x) <= Inf

which is not handled everywhere. That's a nice one, though. I think we need to return an error when this happens.

Can you try with lcon and ucon to see if this is the bug?

[tpapp]

Yes, it is a minor bug (non-informative error message), I opened an issue. That said, not providing lcon and ucon is of course user error.

[tpapp]

Thanks for all the help so far. I updated the repository https://github.com/tpapp/ObjConsNLPModels.jl with new methods and fixed the ucon and lcon. Currently I am stuck with

julia> using ObjConsNLPModels

julia> using Percival

julia> model = objcons_nlpmodel(x -> [sum(abs2, x), sum(x) - 1]; x0 = [2.0, 2.0])

julia> percival(model)

giving me a DimensionMismatch in one of the updating functions. Maybe I am still violating some implicit assumption about dimensions, however it seems not to be caught by the @lenchecks. Help would be appreciated.

[tmigot]

I think the error comes from your objcons! function.
https://github.com/tpapp/ObjConsNLPModels.jl/blob/ea697b41e69699f7d989d10390bc6b6895a635f0/src/ObjConsNLPModels.jl#L169

function NLPModels.objcons!(model::ObjConsNLPModel, x::AbstractVector, buffer::AbstractVector)
    copy!(buffer, objcons_at_point(model, x).objcons)
end

If I understand correctly objcons_at_point(model, x).objcons return an array whose first component is the objective value. So, buffer is of size ncon + 1 while it should be of size ncon.

Additionally, the function objcons! must return a tuple as shown here https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/283ed28037d3ddb32366b65b3e001d7f546903bb/src/nlp/api.jl#L131 .

Let me know if it helps. Could you also provide more details on the error for next time? Thanks!

[tpapp]

This helped a lot. I fixed the code, and submitted a trivial docstring PR to clarify that a tuple needs to be returned (JuliaSmoothOptimizers/NLPModels.jl#436).