# Nested optimization - matrix multiplication and internal functions

**URL:** <https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [August 12, 2024, 3:15pm UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098 "2024-08-12T15:15:07Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 12, 2024, 3:15pm UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098/1 "2024-08-12T15:15:07Z")

</div>

I have a nested optimization problem in JuMP that fails to run because of Dual going through. In my original problem the lower level problem is a convex QP so it is possible to [compute correct gradients](https://jump.dev/DiffOpt.jl/dev/examples/matrix-inversion-manual/) but first the code needs to run. For simplicity I have set a scalar outer variable, but the problem fails in an even stranger way:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/e/7e7e959ee72e73adfaf01e4cab746fc8da28b577.png)

The problem for the MWE:  
\min\limits\_x 2x^2+0.5V(x) where x\geq 0 and V(x)=\max\limits\_y \:y^T \begin{bmatrix} x&0\\0&x \end{bmatrix} y s.t. (y\_1-1)^2+(y\_2-1)^2\leq 2, y\_i-x\leq 100, y\_i-x\geq -100, i=1,2

Case 1: outer variable x is a vector. Then the code runs, but using internal functions fails and matrix multiplication fails:

> **Code case 1**
>
> ```julia
> using JuMP, Ipopt
> 
> function some_function(a,x...)
> return a*x
> end
> 
> function solve_lower_level(x...)
> sm_fcn = x -> 1.0*x
> 
> model = Model(Ipopt.Optimizer)
> set_silent(model)
> @variable(model, y[1:2])
> # @expression(model,testexp,[sm_fcn(x...)[1] 0.0; 0.0 x[2]]') ####Doesn't run
> # @expression(model,testexp,[some_function(1.0,x...) 0.0; 0.0 x[2]]') ####Doesn't run
> @expression(model,testexp,[1.0*x[1] 0.0; 0.0 x[2]]') ####Runs
> 
> @objective(
> model,
> Max,
> y'*[1.0*x[1] 0.0; 0.0 x[2]]'*y
> )
> 
> Cm = [1 0;0 -1]
> Dm = [100;100]
> 
> @constraint(model, (y[1] - 1)^2 + (y[2] - 1)^2 <= 2)
> @constraint(model, Cm[1,1]*(y[1]-x[1])<=Dm[1]) ######If use this, the code runs
> @constraint(model, Cm[2,2]*(y[2]-x[2])<=Dm[2])
> 
> # @constraint(model, Cm*(y-x)<=Dm) #####Supposed to be -100<=y[i]-x[i]<=100, i=1,2, but this doesn't work
> 
> optimize!(model)
> @assert is_solved_and_feasible(model)
> 
> # @show value.(y)
> return objective_value(model), value.(y)
> end
> 
> function V(x...)
> f, _ = solve_lower_level(x...)
> return f
> end
> 
> function ∇V(g::AbstractVector, x...)
> _, y = solve_lower_level(x...)
> g[1] = -y[1]^2
> g[2] = -y[1]^2
> return
> end
> 
> function ∇²V(H::AbstractMatrix, x...)
> _, y = solve_lower_level(x...)
> H[1, 1] = 0.0
> H[2, 2] = 0.0
> return
> end
> 
> model = Model(Ipopt.Optimizer)
> @variable(model, x[1:2] >= 0)
> @operator(model, op_V, 2, V, ∇V, ∇²V)
> @objective(model, Min, x[1]^2 + x[2]^2 + 0.5*op_V(x[1], x[2]))
> optimize!(model)
> @assert is_solved_and_feasible(model)
> solution_summary(model)
> 
> ```

Case one: outer variable x is a scalar. Nothing runs, so I must be making a JuMP error:

> **Code case 2**
>
> ```julia
> using JuMP, Ipopt
> 
> function some_function(a,x)
> return a*x
> end
> 
> function solve_lower_level(x)
> sm_fcn = x -> 1.0*x
> 
> model = Model(Ipopt.Optimizer)
> set_silent(model)
> @variable(model, y[1:2])
> # @expression(model,testexp,[sm_fcn(x) 0.0; 0.0 x]') ####Doesn't run
> # @expression(model,testexp,[some_function(1.0,x) 0.0; 0.0 x]') ####Doesn't run
> # @expression(model,testexp,[1.0*x 0.0; 0.0 x]') ####Doesn't run
> # @expression(model,testexp,[x 0.0; 0.0 x]') ####Doesn't run
> 
> @objective(
> model,
> Max,
> y[1]*y[1]+ y[2]*y[2]
> )
> 
> Cm = [1 0;0 -1]
> Dm = [100;100]
> 
> @constraint(model, (y[1] - 1)^2 + (y[2] - 1)^2 <= 2)
> # @constraint(model, Cm[1,1]*(y[1].-x)<=Dm[1]) ######this doesn't work
> # @constraint(model, Cm[2,2]*(y[2].-x)<=Dm[2])
> 
> # @constraint(model, Cm*(y-x)<=Dm) #####Supposed to be -100<=y[i]-x[i]<=100, i=1,2, but this doesn't work
> 
> optimize!(model)
> @assert is_solved_and_feasible(model)
> 
> # @show value.(y)
> return objective_value(model), value.(y)
> end
> 
> function V(x...)
> f, _ = solve_lower_level(x)
> return f
> end
> 
> function ∇V(g::AbstractVector, x)
> _, y = solve_lower_level(x)
> g[1] = -y[1]^2
> return
> end
> 
> function ∇²V(H::AbstractMatrix, x...)
> _, y = solve_lower_level(x...)
> H[1, 1] = 0.0
> return
> end
> 
> model = Model(Ipopt.Optimizer)
> @variable(model, x >= 0)
> @operator(model, op_V, 1, V, ∇V, ∇²V)
> @objective(model, Min, x^2 + x^2 + 0.5*op_V(x))
> optimize!(model)
> @assert is_solved_and_feasible(model)
> solution_summary(model)
> 
> ```

Thanks for any inputs!

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [August 12, 2024, 3:30pm UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098/2 "2024-08-12T15:30:13Z")

</div>

As usual, please give the complete error message, in a Markdown code block (with ````julia-repl`), instead of a screenshot. See the PSA for more such advice:

> [@Please read: make it easier to help you](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757):
>
> Welcome to the Julia Discourse! We are enthusiastic about helping Julia programmers, both beginner and experienced. This public service announcement (PSA) outlines best practices when asking for help. Following these points makes it easier for us to help you and more likely you’ll get a prompt, useful answer. Keywords are highlighted to make it easier to refer to specific points. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of …

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 12, 2024, 7:01pm UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098/3 "2024-08-12T19:01:25Z")

</div>

The syntax for univariate functions is different to multivariate.

See [Nonlinear Modeling · JuMP](https://jump.dev/JuMP.jl/stable/manual/nonlinear/#Univariate-functions)

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 13, 2024, 7:02am UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098/4 "2024-08-13T07:02:45Z")

</div>

Thanks! I still cannot figure it out even for a single variable. The code I have runs if I don’t use the upper level x in the lower level problem. If I use x in the lower level problem, for instance for multiplication, then I get an error.

```julia
using JuMP, Ipopt

function solve_lower_level(x)

    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model, y[1:2])

    # @objective(
    # model,
    # Min,
    # y[1]*y[1]+ y[2]*y[2]
    # ) ####This works

    @expression(model,testExp, y[1]*x*y[1]+ y[2]*x*y[2]) ###This doesn't

    @objective(
        model,
        Min,
        testExp
    )
    @constraint(model, (y[1] - 1)^2 + (y[2] - 1)^2 <= 2)

    optimize!(model)
    @assert is_solved_and_feasible(model)

    return objective_value(model), value.(y)
end

function V(x)
    f, _ = solve_lower_level(x)
    return f
end

function ∇V(x)
    _, y = solve_lower_level(x)
    return -y[1]^2
end

model = Model(Ipopt.Optimizer)
@variable(model, x >= 0)
@operator(model, op_V, 1, V, ∇V)
@objective(model, Min, x^2 + x^2 + 0.5*op_V(x))
optimize!(model)
@assert is_solved_and_feasible(model)
solution_summary(model)

```

> **Error message**
>
> ```julia
> ERROR: Unable to register the function :op_V.
> 
> Common reasons for this include:
> 
> * The function takes `f(x::Vector)` as input, instead of the splatted
> `f(x...)`.
> * The function assumes `Float64` will be passed as input, it must work for any
> generic `Real` type.
> * The function allocates temporary storage using `zeros(3)` or similar. This
> defaults to `Float64`, so use `zeros(T, 3)` instead.
> 
> As examples, instead of:
> julia
> my_function(x::Vector) = sum(x.^2)
> 
> use:
> julia
> my_function(x::T...) where {T<:Real} = sum(x[i]^2 for i in 1:length(x))
> 
> Instead of:
> julia
> function my_function(x::Float64...)
> y = zeros(length(x))
> for i in 1:length(x)
> y[i] = x[i]^2
> end
> return sum(y)
> end
> 
> use:
> julia
> function my_function(x::T...) where {T<:Real}
> y = zeros(T, length(x))
> for i in 1:length(x)
> y[i] = x[i]^2
> end
> return sum(y)
> end
> 
> Review the stacktrace below for more information, but it can often be hard to
> understand why and where your function is failing. You can also debug this
> outside of JuMP as follows:
> julia
> import ForwardDiff
> 
> # If the input dimension is 1
> x = 1.0
> my_function(a) = a^2
> ForwardDiff.derivative(my_function, x)
> 
> # If the input dimension is more than 1
> x = [1.0, 2.0]
> my_function(a, b) = a^2 + b^2
> ForwardDiff.gradient(x -> my_function(x...), x)
> 
> Stacktrace:
> [1] error(s::String)
> @ Base .\error.jl:35
> [2] _validate_register_assumptions(f::typeof(∇V), name::Symbol, dimension::Int64)
> @ MathOptInterface.Nonlinear C:\.julia\packages\MathOptInterface\aJZbq\src\Nonlinear\operators.jl:327
> [3] MathOptInterface.Nonlinear._UnivariateOperator(op::Symbol, f::Function, f′::Function)
> @ MathOptInterface.Nonlinear C:\.julia\packages\MathOptInterface\aJZbq\src\Nonlinear\operators.jl:347
> [4] register_operator(::MathOptInterface.Nonlinear.OperatorRegistry, ::Symbol, ::Int64, ::Function, ::Vararg{Function})
> @ MathOptInterface.Nonlinear C:\.julia\packages\MathOptInterface\aJZbq\src\Nonlinear\operators.jl:399
> [5] register_operator
> @ C:\.julia\packages\MathOptInterface\aJZbq\src\Nonlinear\model.jl:237 [inlined]
> [6] set(model::Ipopt.Optimizer, attr::MathOptInterface.UserDefinedFunction, args::Tuple{typeof(V), typeof(∇V)})
> @ Ipopt C:\.julia\packages\Ipopt\bqp63\src\MOI_wrapper.jl:490
> [7] set(b::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, attr::MathOptInterface.UserDefinedFunction, value::Tuple{typeof(V), typeof(∇V)})    
> @ MathOptInterface.Bridges C:\.julia\packages\MathOptInterface\aJZbq\src\Bridges\bridge_optimizer.jl:942
> [8] _pass_attribute(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, index_map::MathOptInterface.Utilities.IndexMap, attr::MathOptInterface.UserDefinedFunction)
> @ MathOptInterface.Utilities C:\.julia\packages\MathOptInterface\aJZbq\src\Utilities\copy.jl:51
> [9] pass_attributes(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, index_map::MathOptInterface.Utilities.IndexMap)
> @ MathOptInterface.Utilities C:\.julia\packages\MathOptInterface\aJZbq\src\Utilities\copy.jl:38
> [10] default_copy_to(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}})
> @ MathOptInterface.Utilities C:\.julia\packages\MathOptInterface\aJZbq\src\Utilities\copy.jl:391
> [11] copy_to
> @ C:\.julia\packages\MathOptInterface\aJZbq\src\Bridges\bridge_optimizer.jl:442 [inlined]
> [12] optimize!
> @ C:\.julia\packages\MathOptInterface\aJZbq\src\MathOptInterface.jl:121 [inlined]
> [13] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
> @ MathOptInterface.Utilities C:\.julia\packages\MathOptInterface\aJZbq\src\Utilities\cachingoptimizer.jl:321
> [14] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
> @ JuMP C:\.julia\packages\JuMP\7rBNn\src\optimizer_interface.jl:595
> [15] optimize!(model::Model)
> @ JuMP C:\.julia\packages\JuMP\7rBNn\src\optimizer_interface.jl:546
> [16] top-level scope
> @ c:\Tests_diff3.jl:44
> 
> caused by: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(∇V), Float64}, Float64, 1})
> 
> Closest candidates are:
> (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
> @ Base rounding.jl:207
> (::Type{T})(::T) where T<:Number
> @ Core boot.jl:792
> (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number}
> @ Base char.jl:50
> ...
> 
> Stacktrace:
> [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(∇V), Float64}, Float64, 1})
> @ Base .\number.jl:7
> [2] _complex_convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(∇V), Float64}, Float64, 1})
> @ JuMP C:\.julia\packages\JuMP\7rBNn\src\operators.jl:22
> [3] *(lhs::ForwardDiff.Dual{ForwardDiff.Tag{typeof(∇V), Float64}, Float64, 1}, rhs::VariableRef)
> @ JuMP C:\.julia\packages\JuMP\7rBNn\src\operators.jl:48
> [4] *(lhs::VariableRef, rhs::ForwardDiff.Dual{ForwardDiff.Tag{typeof(∇V), Float64}, Float64, 1})
> @ JuMP C:\.julia\packages\JuMP\7rBNn\src\operators.jl:105
> [5] operate(::typeof(*), ::VariableRef, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(∇V), Float64}, Float64, 1})
> @ MutableArithmetics C:\.julia\packages\MutableArithmetics\SXYDN\src\interface.jl:210
> [6] macro expansion
> @ C:\.julia\packages\MutableArithmetics\SXYDN\src\rewrite.jl:340 [inlined]
> [7] macro expansion
> @ C:\.julia\packages\JuMP\7rBNn\src\macros.jl:257 [inlined]
> [8] macro expansion
> @ C:\.julia\packages\JuMP\7rBNn\src\macros\@expression.jl:86 [inlined]
> [9] macro expansion
> @ C:\.julia\packages\JuMP\7rBNn\src\macros.jl:393 [inlined]
> [10] solve_lower_level(x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(∇V), Float64}, Float64, 1})
> @ Main c:\Tests_diff3.jl:14
> [11] ∇V(x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(∇V), Float64}, Float64, 1})
> @ Main c:\Tests_diff3.jl:36
> [12] derivative
> @ C:\.julia\packages\ForwardDiff\PcZ48\src\derivative.jl:14 [inlined]
> [13] _validate_register_assumptions(f::typeof(∇V), name::Symbol, dimension::Int64)
> @ MathOptInterface.Nonlinear C:\.julia\packages\MathOptInterface\aJZbq\src\Nonlinear\operators.jl:321
> [14] MathOptInterface.Nonlinear._UnivariateOperator(op::Symbol, f::Function, f′::Function)
> @ MathOptInterface.Nonlinear C:\.julia\packages\MathOptInterface\aJZbq\src\Nonlinear\operators.jl:347
> [15] register_operator(::MathOptInterface.Nonlinear.OperatorRegistry, ::Symbol, ::Int64, ::Function, ::Vararg{Function})
> @ MathOptInterface.Nonlinear C:\.julia\packages\MathOptInterface\aJZbq\src\Nonlinear\operators.jl:399
> [16] register_operator
> @ C:\.julia\packages\MathOptInterface\aJZbq\src\Nonlinear\model.jl:237 [inlined]
> [17] set(model::Ipopt.Optimizer, attr::MathOptInterface.UserDefinedFunction, args::Tuple{typeof(V), typeof(∇V)})
> @ Ipopt C:\.julia\packages\Ipopt\bqp63\src\MOI_wrapper.jl:490
> [18] set(b::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, attr::MathOptInterface.UserDefinedFunction, value::Tuple{typeof(V), typeof(∇V)})    
> @ MathOptInterface.Bridges C:\.julia\packages\MathOptInterface\aJZbq\src\Bridges\bridge_optimizer.jl:942
> [19] _pass_attribute(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, index_map::MathOptInterface.Utilities.IndexMap, attr::MathOptInterface.UserDefinedFunction)
> @ MathOptInterface.Utilities C:\.julia\packages\MathOptInterface\aJZbq\src\Utilities\copy.jl:51
> [20] pass_attributes(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, index_map::MathOptInterface.Utilities.IndexMap)
> @ MathOptInterface.Utilities C:\.julia\packages\MathOptInterface\aJZbq\src\Utilities\copy.jl:38
> [21] default_copy_to(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}})
> @ MathOptInterface.Utilities C:\.julia\packages\MathOptInterface\aJZbq\src\Utilities\copy.jl:391
> [22] copy_to
> @ C:\.julia\packages\MathOptInterface\aJZbq\src\Bridges\bridge_optimizer.jl:442 [inlined]
> [23] optimize!
> @ C:\.julia\packages\MathOptInterface\aJZbq\src\MathOptInterface.jl:121 [inlined]
> [24] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
> @ MathOptInterface.Utilities C:\.julia\packages\MathOptInterface\aJZbq\src\Utilities\cachingoptimizer.jl:321
> [25] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
> @ JuMP C:\.julia\packages\JuMP\7rBNn\src\optimizer_interface.jl:595
> [26] optimize!(model::Model)
> @ JuMP C:\.julia\packages\JuMP\7rBNn\src\optimizer_interface.jl:546
> [27] top-level scope
> @ c:\Tests_diff3.jl:44
> 
> ```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 13, 2024, 7:19am UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098/5 "2024-08-13T07:19:08Z")

</div>

Hmm. I think this is because you haven’t provided the second derivative, so we try to automatically complete it using FowardDiff.jl. I should clarify the documentation.

I think you need something like:

```julia
using JuMP, Ipopt

function solve_lower_level(x)
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model, y[1:2])
    @objective(model, Min, y[1] * x * y[1]+ y[2] * x * y[2])
    @constraint(model, (y[1] - 1)^2 + (y[2] - 1)^2 <= 2)
    optimize!(model)
    @assert is_solved_and_feasible(model)
    return objective_value(model), value.(y)
end

function V(x)
    f, _ = solve_lower_level(x)
    return f
end

function ∇V(x)
    _, y = solve_lower_level(x)
    return -(y[1]^2 + y[2]^2)
end

function ∇²V(x)
    _, y = solve_lower_level(x)
    return -2 * y[1] - 2 * y[2]
end

model = Model(Ipopt.Optimizer)
@variable(model, x >= 0)
@operator(model, op_V, 1, V, ∇V, ∇²V)
@objective(model, Min, x^2 + x^2 + 0.5 * op_V(x))
optimize!(model)
@assert is_solved_and_feasible(model)
solution_summary(model)

```

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 13, 2024, 7:29am UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098/6 "2024-08-13T07:29:04Z")

</div>

Yes, that runs, thanks. The thing is that I don’t have second derivatives to my original problem so I thought I could get away without it ☹

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 13, 2024, 8:04am UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098/7 "2024-08-13T08:04:47Z")

</div>

For future reference, I found the documentation bit where it says that ForwardDiff cannot go through a JuMP problem on its own: [Nonlinear Modeling · JuMP](https://jump.dev/JuMP.jl/stable/manual/nonlinear/#Operator-calls-something-unsupported-by-ForwardDiff)

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 13, 2024, 9:04am UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098/8 "2024-08-13T09:04:04Z")

</div>

Perhaps you could try:

```Julia
function V(x, tmp)
    f, _ = solve_lower_level(x)
    return f
end

function ∇V(g, x, tmp)
    _, y = solve_lower_level(x)
    g[1] = -(y[1]^2 + y[2]^2)
    g[2] = 0.0
    return
end

model = Model(Ipopt.Optimizer)
@variable(model, x >= 0)
@operator(model, op_V, 2, V, ∇V)
@objective(model, Min, x^2 + x^2 + 0.5 * op_V(x, 0))

```

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 13, 2024, 9:23am UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098/9 "2024-08-13T09:23:09Z")

</div>

That works, thank you! Why it works, that’s a different question 😅 From all I understand, we shouldn’t even be able to pass two arguments to op\_V without `x...`.

Thanks again!

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 13, 2024, 9:25am UTC](https://discourse.julialang.org/t/nested-optimization-matrix-multiplication-and-internal-functions/118098/10 "2024-08-13T09:25:41Z")

</div>

> Why it works

We have a different implementation for uni-variate and multi-variate functions. This is confusing. We should probably be able to opt-in to the multivariate implementation even for `N = 1`.

Edit: I’ve opened an issue [Univariate user-defined functions with the multivariate signature · Issue #2534 · jump-dev/MathOptInterface.jl · GitHub](https://github.com/jump-dev/MathOptInterface.jl/issues/2534)

> we shouldn’t even be able to pass two arguments

This is just a regular Julia thing. These two are the same:

```julia
foo(x, y) = x[1] + x[2]
foo(x...) = x[1] + x[2]

```

Although the `x...` doesn’t give as nice an error if you call it with the wrong number of arguments.