Repository navigation
The dual of VectorNonlinearOracle #2887
Description
Activity
dualshould 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$$ Reacted by Francis GagnonSo, from my understanding, if
dualreturns$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$ ?I don't think that is well-defined. We should just return the multipliers.
Reacted by Francis GagnonOur 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
xpart 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 Swheret - f(x) == 0. Now the duals all migrate their way out of the set and onto the variable bounds oft: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
Reacted by Benoît LegatI 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".But we have "x-in-some weird set" rather than "some weird function-in-interval".
Yes, precisely
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:
- you use
dual(oracle_constraint)internally and the code expects the same dimension ? or - 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
dualto fetch the Lagrange multipliers ? My knowledge on LP is very poor so I may suggest something dumb, but, maybeshadow_pricesor 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
getinfofunction 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 to "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.- you use
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.For my part, I often use some code like this (where
nlpis 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]
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' * λ.But what would be the consequence of breaking this convention only for
VectorNonlinearOracle?The ability to work with abstract sets and functions in MOI requires us to be quite strict on this so I prefer not making an exception.
My use case is a
getinfofunction that can be called the compute additional information about the optimum, to help troubleshooting/debugging. This function should not be called in "normal" operation.Then this can be solved by a solver-specific constraint attributed defined by Ipopt, something like jump-dev/Ipopt.jl#521. We could standardize it at some point into MOI if we converge on this.
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".I plan to introduce array operators in https://github.com/blegat/ArrayDiff.jl/ so that might be a possible direction then. But it would be similar to the user-defined operators in
ScalarNonlinearFunction, not a new MOI function typeVectorNonlinearOperator.Reacted by Robert ParkerReacted by Francis GagnonThen this can be solved by a solver-specific constraint attributed defined by Ipopt, something like jump-dev/Ipopt.jl#521. We could standardize it at some point into MOI if we converge on this.
That would be a good solution for me! I guess it would need to be manually implemented in all NLP solver package e.g. KNITRO.jl, UnoSolver.jl, etc. ? As long as a clear exception is thrown if the attribute is missing, I'm happy with this solution.
See https://discourse.julialang.org/t/lagrange-multipliers-of-a-vectornonlinearoracle-after-solve/134236
What should
dual(c)here be?I based it on:
but the user doesn't have an easy way to access the multipliers on the individual rows.
cc @Robbybp