Legacy and standard constraint/objective macros represent `^2` differently in nonlinear expressions

[DimitriAlston]

@objective and @constraint seem to be reformulating x^2 to x*x in nonlinear expressions, which leads to poorer McCormick relaxations. This isn't an issue when using the legacy macros of @NLobjective and @NLconstraint. I suspect this happens somewhere in

EAGO.jl/src/eago_optimizer/moi_wrapper.jl

Lines 74 to 100 in 15bb6cb

	function MOI.add_constraint(m::Optimizer, f::F, s::S) where {F<:Union{SAF,SQF,MOI.ScalarNonlinearFunction},S<:INEQ_SETS}
	check_inbounds!(m, f)
	ci = CI{F, S}(m._input_problem._constraint_count += 1)
	if f isa MOI.ScalarNonlinearFunction
	if isnothing(m._input_problem._nlp_data)
	model = MOI.Nonlinear.Model()
	backend = MOI.Nonlinear.SparseReverseMode()
	vars = MOI.get(m, MOI.ListOfVariableIndices())
	evaluator = MOI.Nonlinear.Evaluator(model, backend, vars)
	m._input_problem._nlp_data = MOI.NLPBlockData(evaluator)
	end
	MOI.Nonlinear.add_constraint(m._input_problem._nlp_data.evaluator.model, f, s)
	constraint_bounds = m._input_problem._nlp_data.constraint_bounds
	has_objective = m._input_problem._nlp_data.has_objective
	if s isa MOI.LessThan
	push!(constraint_bounds, MOI.NLPBoundsPair(-Inf, s.upper))
	elseif s isa MOI.GreaterThan
	push!(constraint_bounds, MOI.NLPBoundsPair(s.lower, Inf))
	else
	push!(constraint_bounds, MOI.NLPBoundsPair(s.value, s.value))
	end
	m._input_problem._nlp_data = MOI.NLPBlockData(constraint_bounds, m._input_problem._nlp_data.evaluator, has_objective)
	else
	_constraints(m, F, S)[ci] = (f, s)
	end
	return ci
	end

and

EAGO.jl/src/eago_optimizer/moi_wrapper.jl

Lines 277 to 294 in 15bb6cb

	function MOI.set(m::Optimizer, ::MOI.ObjectiveFunction{T}, f::T) where T <: Union{VI,SAF,SQF,MOI.ScalarNonlinearFunction}
	check_inbounds!(m, f)
	if f isa MOI.ScalarNonlinearFunction
	if isnothing(m._input_problem._nlp_data)
	model = MOI.Nonlinear.Model()
	MOI.Nonlinear.set_objective(model, f)
	backend = MOI.Nonlinear.SparseReverseMode()
	vars = MOI.get(m, MOI.ListOfVariableIndices())
	evaluator = MOI.Nonlinear.Evaluator(model, backend, vars)
	m._input_problem._nlp_data = MOI.NLPBlockData(evaluator)
	else
	MOI.Nonlinear.set_objective(m._input_problem._nlp_data.evaluator.model, f)
	m._input_problem._nlp_data = MOI.NLPBlockData(m._input_problem._nlp_data.evaluator)
	end
	else
	m._input_problem._objective = f
	end
	end

Example:

model = JuMP.Model(EAGO.Optimizer)
@variable(model, x[1:2])
@objective(model, Min, x[1]^2 + x[2]^3)
JuMP.optimize!(model)

returns

MathOptInterface.Nonlinear.Expression(
    MathOptInterface.Nonlinear.Node[
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_CALL_MULTIVARIATE, 2, -1),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_CALL_MULTIVARIATE, 1, 1),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_CALL_MULTIVARIATE, 3, 2),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_MOI_VARIABLE, 1, 3),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_MOI_VARIABLE, 1, 3),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_CALL_MULTIVARIATE, 4, 2),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_MOI_VARIABLE, 2, 6),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_VALUE, 1, 6),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_MOI_VARIABLE, 3, 1)
    ],
    [3.0]
)

model = JuMP.Model(EAGO.Optimizer)
@variable(model, x[1:2])
@NLobjective(model, Min, x[1]^2 + x[2]^3)
JuMP.optimize!(model)

returns

MathOptInterface.Nonlinear.Expression(
    MathOptInterface.Nonlinear.Node[
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_CALL_MULTIVARIATE, 2, -1),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_CALL_MULTIVARIATE, 1, 1),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_CALL_MULTIVARIATE, 4, 2),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_MOI_VARIABLE, 1, 3),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_VALUE, 1, 3),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_CALL_MULTIVARIATE, 4, 2),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_MOI_VARIABLE, 2, 6),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_VALUE, 2, 6),
        MathOptInterface.Nonlinear.Node(MathOptInterface.Nonlinear.NODE_MOI_VARIABLE, 3, 1)
    ],
    [2.0, 3.0]
)

(i.e., x^2 + x^3, as expected).

[DimitriAlston]

@odow can you tell if this is an EAGO or a JuMP issue?

[odow]

This is coming from MOI:

https://github.com/jump-dev/MathOptInterface.jl/blob/406f7cd9b0e681f1934121055b7619ab3d666202/src/Nonlinear/parse.jl#L342-L361

The reason is that x[1]^2 is treated as a QuadExpr, not ^(x[1], 2).

julia> using JuMP

julia> model = Model();

julia> @variable(model, x[1:2]);

julia> @objective(model, Min, x[1]^2 + x[2]^3);

julia> f = objective_function(model)
(x[1]²) + (x[2] ^ 3.0)

julia> f.head
:+

julia> f.args
2-element Vector{Any}:
 x[1]²
 x[2] ^ 3.0

julia> f.args[1]
x[1]²

julia> typeof(f.args[1])
QuadExpr (alias for GenericQuadExpr{Float64, GenericVariableRef{Float64}})

Isn't there scope in EAGO to recognize this expression type? If someone writes x * x you don't currently exploit that as x^2?

[odow]

It's not really a fix, but a work-around is

julia> using JuMP

julia> model = Model();

julia> @variable(model, x[1:2]);

julia> @objective(model, Min, x[1]^2 + x[2]^3)
(x[1]²) + (x[2] ^ 3)

julia> @force_nonlinear @objective(model, Min, x[1]^2 + x[2]^3)
(x[1] ^ 2) + (x[2] ^ 3)

[DimitriAlston]

Isn't there scope in EAGO to recognize this expression type? If someone writes x * x you don't currently exploit that as x^2?

I believe EAGO used to provide a functionality like this before v0.7 for user-defined functions, but it hasn't been updated after a lot of other internal changes. After the MOI changes for handling nonlinear functions, I'd written a recursive function to convert ScalarNonlinearFunctions to a DAG in EAGO, which would have checked for things like this, but after some discussion at the last JuMP-dev I switched to using a workaround in MOI to convert back to the old DAG format. So currently we rely on MOI's expression representation, though this might change if we get more time to develop out this functionality in EAGO.

[odow]

I've opened a PR here: jump-dev/MathOptInterface.jl#2799

It's not guaranteed that I'll merge it, but I'll see if there are any other downstream impacts.

[odow]

This should be fixed in the next release of MOI