The dual of VectorNonlinearOracle

[odow]

See https://discourse.julialang.org/t/lagrange-multipliers-of-a-vectornonlinearoracle-after-solve/134236

julia> using JuMP, Ipopt

julia> set = MOI.VectorNonlinearOracle(;
           dimension = 2,
           l = [-Inf],
           u = [1.0],
           eval_f = (ret, x) -> (ret[1] = x[1]^2 + x[2]^2),
           jacobian_structure = [(1, 1), (1, 2)],
           eval_jacobian = (ret, x) -> ret .= 2.0 .* x,
           hessian_lagrangian_structure = [(1, 1), (2, 2)],
           eval_hessian_lagrangian = (ret, x, u) -> ret .= 2.0 .* u[1],
       )
VectorNonlinearOracle{Float64}(;
    dimension = 2,
    l = [-Inf],
    u = [1.0],
    ...,
)

julia> model = Model(Ipopt.Optimizer)
A JuMP Model
├ solver: Ipopt
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

julia> set_silent(model)

julia> @variable(model, x[1:2])
2-element Vector{VariableRef}:
 x[1]
 x[2]

julia> @objective(model, Max, sum(x))
x[1] + x[2]

julia> @constraint(model, c, x in set)
c : [x[1], x[2]] ∈ VectorNonlinearOracle{Float64}(;
    dimension = 2,
    l = [-Inf],
    u = [1.0],
    ...,
)

julia> optimize!(model)

julia> value(x)
2-element Vector{Float64}:
 0.7071067834685273
 0.7071067834685273

julia> dual(c)
2-element Vector{Float64}:
 -0.9999999999982961
 -0.9999999999982961

What should dual(c) here be?

I based it on:

julia> using JuMP, Ipopt

julia> model = Model(Ipopt.Optimizer)
A JuMP Model
├ solver: Ipopt
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

julia> set_silent(model)

julia> @variable(model, x[1:2])
2-element Vector{VariableRef}:
 x[1]
 x[2]

julia> @objective(model, Max, sum(x))
x[1] + x[2]

julia> @constraint(model, c, sum(x.^2) <= 1)
c : x[1]² + x[2]² ≤ 1

julia> optimize!(model)

julia> value(x)
2-element Vector{Float64}:
 0.7071067834685273
 0.7071067834685273

julia> dual(c)
-0.7071067789033629

julia> dual(c)' * (2 * value(x))
2-element Vector{Float64}:
 -0.9999999999982961
 -0.9999999999982961

but the user doesn't have an easy way to access the multipliers on the individual rows.

cc @Robbybp

[Robbybp]

dual should give the multipliers, not the multipliers times the Jacobian. They should have the same dimension as the output of the function. That you don't have to handle this constraint specially when constructing the (gradient of the) Lagrangian:
$$\nabla \mathcal{L} = \nabla f + J^T \lambda$$

[franckgaga]

So, from my understanding, if dual returns $v = \mu^\top J(x)$, I can compute back the Lagrange multiplier with:

$$\mu^\top = v J^{-1}(x)$$

This is a relatively cheap computation, so I'm okay with this solution. But is it always possible ? Presumably not if $J(x)$ has more rows than columns i.e.: an inconsistent system. This is an issue since we really need the $\mu$ vector to compute the Hessian of the Lagrangian at optimum. Or can it be computed from $v$ ?

[Robbybp]

I don't think that is well-defined. We should just return the multipliers.

[odow]

I think @blegat needs to chime in here.

I get that the Lagrange multipliers are a metric that folks might want, but it isn't obvious to me that they are the MOI.ConstraintDual of this set. I need to think about this a bit more.

[odow]

Our issue is that we have (from MOI's view) the set $x \in S \subseteq \mathbb{R}^n$. I think the dual vector needs to have the same dimension, so it doesn't make sense to return $ y \in \mathbb{R}^m$.

Here's the conic equivalent of our example.

julia> using JuMP, SCS

julia> begin
           model = Model(SCS.Optimizer)
           set_silent(model)
           @variable(model, x[1:2])
           @objective(model, Max, sum(x))
           @constraint(model, c, [1; 0.5; x] in RotatedSecondOrderCone())
           optimize!(model)
           value(x), dual(c)
       end
([0.7071113337601804, 0.7071113337601804], [0.7071208604303918, 1.414186043597682, -1.0000002259786969, -1.0000002259786969])

The "dual" on the x part of the RSOC is [-1, -1], which matches the "dual" we return at present.

The question is really, how do we get access the the sqrt(2) components of the multipliers.

One answer is to make the set [t; x] in S where t - f(x) == 0. Now the duals all migrate their way out of the set and onto the variable bounds of t:

julia> using JuMP, Ipopt

julia> set = MOI.VectorNonlinearOracle(;
           dimension = 3,
           l = [0.0],
           u = [0.0],
           eval_f = (ret, x) -> (ret[1] = x[1] - x[2]^2 - x[3]^2),
           jacobian_structure = [(1, 1), (1, 2), (1, 3)],
           eval_jacobian = (ret, x) -> (ret[1] = 1.0; ret[2:3] .= -2.0 .* x[2:3]),
           hessian_lagrangian_structure = [(2, 2), (3, 3)],
           eval_hessian_lagrangian = (ret, x, u) -> ret .= -2.0 .* u[1],
       )
VectorNonlinearOracle{Float64}(;
    dimension = 3,
    l = [0.0],
    u = [0.0],
    ...,
)

julia> model = Model(Ipopt.Optimizer)
A JuMP Model
├ solver: Ipopt
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

julia> set_silent(model)

julia> @variable(model, x[1:2])
2-element Vector{VariableRef}:
 x[1]
 x[2]

julia> @variable(model, t <= 1.0)
t

julia> @objective(model, Max, sum(x))
x[1] + x[2]

julia> @constraint(model, c, [t; x] in set)
c : [t, x[1], x[2]] ∈ VectorNonlinearOracle{Float64}(;
    dimension = 3,
    l = [0.0],
    u = [0.0],
    ...,
)

julia> optimize!(model)

julia> value(t)
1.0000000064527088

julia> dual(UpperBoundRef(t))
-0.7071067789033879

julia> dual(c)
3-element Vector{Float64}:
  0.7071067789033629
 -0.9999999999982961
 -0.9999999999982961

[Robbybp]

I guess the way I think about it, in MOI-speak, is more like we have VectorNonlinearOperator(x)$\in [l, u]$, so the duals would be in $\mathbb{R}^m$. But we have "x-in-some weird set" rather than "some weird function-in-interval".

[odow]

But we have "x-in-some weird set" rather than "some weird function-in-interval".

Yes, precisely

[franckgaga]

Our issue is that we have (from MOI's view) the set x ∈ S ⊆ R n . I think the dual vector needs to have the same dimension, so it doesn't make sense to return $ y \in \mathbb{R}^m$.

Why does it needs to have the same dimension ? Is it because:

1. you use dual(oracle_constraint) internally and the code expects the same dimension ? or
2. to be consistent with the MOI convention used throughout all other kind of constraint ?

If the answer is 2: I'm the first one to fully agree that conventions are important in a complex codebase like MOI. But what would be the consequence of breaking this convention only for VectorNonlinearOracle ?

Otherwise, could we could define another method instead of dual to fetch the Lagrange multipliers ? My knowledge on LP is very poor so I may suggest something dumb, but, maybe shadow_prices or something similar ?

One answer is to make the set [t; x] in S where t - f(x) == 0. Now the duals all migrate their way out of the set and onto the variable bounds of t:

It's good to know that's possible but that's really not an intuitive way of fetching the multipliers. My use case is a getinfo function that can be called the compute additional information about the optimum, to help troubleshooting/debugging. This function should not be called in "normal" operation. I would be very reluctant the "augment" the optimization problem just to be able to fetch the multipliers, in cases that the user calls my debugging function. Especially knowing that it can impact the performances.

[odow]

It's 2. We're not breaking conventions by having duals of different dimensions. It's arguably also 1, but there are probably other issues in the code.

Otherwise, could we could define another method instead of dual to fetch the Lagrange multipliers ?

Yeah, we could have some other attribute. That's not a problem. There's also something else to consider: we don't even define that the set will have a dual or a Lagrange multiplier for the rows; the Hessian is optional.

The underlying issue is that we've defined a new set (because that's easy to do in MOI), rather than a new function (because that's hard to do in MOI). If we had done MOI.VectorNonlinearOracleFunction -in- MOI.Interval, there would be no issue. But although it's conceptually easy to think we could have done that, implementing it in practice in JuMP, MOI, and Ipopt is quite tricky.

[Robbybp]

For my part, I often use some code like this (where nlp is a thin wrapper around the NLEvaluator):

function eval_lagrangian_gradient(nlp, x, λ)
    # x comes from JuMP.value.(all_variables)
    # λ comes from JuMP.dual.(all_constraints)
    grad_obj = eval_objective_gradient(nlp, x)
    # jac::SparseMatrixCSC
    jac = eval_constraint_jacobian(nlp, x)
    grad_lagrangian = grad_obj * nlp.lagrangian_objective_factor + jac' * λ
    return grad_lagrangian
end

to evaluate $\nabla \mathcal{L}$. I guess this doesn't work if I have VectorNonlinearOracle constraints. I will need something like:

grad_lagrangian = grad_obj * nlp.lagrangian_objective_factor + jac[normal_indices, :]' * λ[normal_indices] + λ[vno_indices]

[odow]

@Robbybp: isn't this an exact argument for why our MOI.ConstraintDual is correct? You can just use λ[vno_indices] directly.

[Robbybp]

You can just use λ[vno_indices] directly.

Personally, I would rather it be consistent with all the other constraints that nonlinear solvers support, where I need jac' * λ.

[odow]

In our case, the function is now x, so the Jacobian is the identity matrix.

I need to get @blegat's opinion.