What is the difference between nonlinear "functions with vector outputs" and nonlinear "operators with vector outputs"?

Specific Domains Optimization (Mathematical)
,
1

The JuMP documentation and this answer suggests that nonlinear user-defined functions with vector outputs need to be defined using many @operators. The documentation actually says “User-defined operators …”.

This post suggest that nonlinear functions can be used directly. The date on that post and syntax suggest an old version of JuMP.

The example below also suggest that nonlinear functions can be used directly.

This works

import Pkg; Pkg.activate(temp=true)
Pkg.add(name="JuMP", version="1.15.1");
Pkg.add("Ipopt");
Pkg.instantiate();

import JuMP
import Ipopt

model = JuMP.Model(Ipopt.Optimizer)

function dy(x)
    return [cos(x[1])*x[2]; sin(x[1]) + sin(x[2])]
end

x0 = [.1;.1]
xc = [1;0.5]

JuMP.@variable(model, x[i=1:2], start = x0[i])
JuMP.@constraint(model, dy(x) == xc)
JuMP.optimize!(model)

print("x: ", JuMP.value.(x))
print("dy(x): ", dy(JuMP.value.(x)))

...
EXIT: Optimal Solution Found.
x: [-0.39325899095511446, 1.0826434588266454]dy(x): [0.9999999999993564, 0.5000000000001374]

So does this

import Pkg; Pkg.activate(temp=true)
Pkg.add(name="JuMP", version="1.15.1");
Pkg.add("Ipopt");
Pkg.instantiate();

import JuMP
import Ipopt

function dy(x)
    return [cos(x[1])*x[2]; sin(x[1]) + sin(x[2])]
end

x0 = [.1;.1]
xc = [1;0.5]

model = JuMP.Model(Ipopt.Optimizer)
JuMP.@variable(model, x[i=1:2], start = x0[i])

JuMP.@operator(model, op_f1, 2, (x...) -> dy(collect(x))[1])
JuMP.@operator(model, op_f2, 2, (x...) -> dy(collect(x))[2])

JuMP.@constraint(model, op_f1(x...) == xc[1])
JuMP.@constraint(model, op_f2(x...) == xc[2])

JuMP.optimize!(model)

print("x: ", JuMP.value.(x))
print("dy(x): ", dy(JuMP.value.(x)))

...
EXIT: Optimal Solution Found.
x: [-0.39325899095511446, 1.0826434588266454]dy(x): [0.9999999999993564, 0.5000000000001374]

Is there a preferred syntax?

Is there a reason to prefer the @operator macro over the user-defined function?

If we have a nice conclusion I will open a PR to add to the documentation :smiley:

PS:
This for nonlinear and this for linear are possibly related.

2

It is possible that my confusion lies entirely with when to use an operator and is unrelated to “vector outputs”. The below example illustrates the need (?) for an @operator in the case of taking a derivative.

This does not work

import Pkg; Pkg.activate(temp=true)
Pkg.add(name="JuMP", version="1.15.1");
Pkg.add("Ipopt");
Pkg.add("ForwardDiff");
Pkg.instantiate();

import JuMP
import Ipopt
import ForwardDiff as FD

f(x) = FD.derivative(sin,x)

model = JuMP.Model(Ipopt.Optimizer)
JuMP.@variable(model, x, start = 0.4)
JuMP.@constraint(model, f(x) == 0.5)
JuMP.optimize!(model)

print("x: ", JuMP.value.(x))

MethodError: no method matching derivative(::typeof(sin), ::JuMP.VariableRef)

Closest candidates are:
  derivative(::F, ::R) where {F, R<:Real}
   @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/derivative.jl:12
  derivative(::Any, ::AbstractArray)
   @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/derivative.jl:72
  derivative(::Any, ::AbstractArray, ::Real)
   @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/derivative.jl:25
  ...


Stacktrace:
 [1] f(x::JuMP.VariableRef)
   @ Main ./In[2]:11
 [2] macro expansion
   @ ~/.julia/packages/MutableArithmetics/K9YPJ/src/rewrite.jl:321 [inlined]
 [3] macro expansion
   @ ~/.julia/packages/JuMP/ToPd2/src/macros.jl:717 [inlined]
 [4] top-level scope
   @ In[2]:15

This works

import Pkg; Pkg.activate(temp=true)
Pkg.add(name="JuMP", version="1.15.1");
Pkg.add("Ipopt");
Pkg.add("ForwardDiff");
Pkg.instantiate();

import JuMP
import Ipopt
import ForwardDiff as FD

f(x) = FD.derivative(sin,x)

model = JuMP.Model(Ipopt.Optimizer)
JuMP.@operator(model, op_f, 1, f)

JuMP.@variable(model, x, start = 0.4)
JuMP.@constraint(model, op_f(x) == 0.5)
JuMP.optimize!(model)

print("x: ", JuMP.value.(x))

...
EXIT: Optimal Solution Found.
x: 1.0471975562573468

The macro and struct suggests we need to use an operator " When called with non-JuMP types, the struct returns the evaluation of func(args...) .".

So the answer to question might be.

  1. vector outputs are fine to use as is for JuMP types.
  2. use @operators when multiple dispatch has not been implemented for JuMP types.

Or I could be wrong again :smiley:

1 Like
3

You’re on the right track.

For some simple algebraic functions, like your dy, you can call it directly, including vector valued outputs. Under the hood, JuMP traces the function to create a vector of nonlinear expressions. (Try calling the function in the REPL.)

For functions that we cannot trace, like your gradient function, you must define it as an operator. We either need to be able to automatically differentiate it, or you need to provide the gradient. User defined operators must return a scalar.

If you have any suggestions for how we could improve the documentation please let me know. Its always hard for me to know what bits are confusing :slight_smile:

1 Like