# Efficient Dual Problems in JuMP? Or: Avoiding Symbolic Overheads with Complex Objectives?

**URL:** <https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, performance\
**Created:** [May 7, 2024, 3:55pm UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953 "2024-05-07T15:55:32Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![SebKrantz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebkrantz/32/50256_2.png) [@SebKrantz](https://discourse.julialang.org/u/SebKrantz)\
**Post date:** [May 7, 2024, 3:55pm UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/1 "2024-05-07T15:55:32Z")

</div>

I am developing [OptimalTransportNetworks,jl](https://github.com/SebKrantz/OptimalTransportNetworks.jl), a Julia library to compute optimal transport networks in spatial equilibirium originally implemented in MATLAB.

Having succesfully implemented all model cases in JuMP, I am now examining their performance where I observe that the [Primal Solution](https://github.com/SebKrantz/OptimalTransportNetworks.jl/blob/main/src/models/model_fixed.jl) to the simplest model case without population mobility and cross-good congestion, optimizing across 2 x J x N + ndeg x N variables (J = number of nodes and ndeg = number of edges of the underlying graph and N = number of goods traded), is significantly slower than the identical [Dual Solution](https://github.com/SebKrantz/OptimalTransportNetworks.jl/blob/main/src/models/model_fixed_duality.jl) optimizing over just J x N variables (the prices of each good in each location which are the shadow values of the constraints in the primal problem). For example executing [Example 1](https://github.com/SebKrantz/OptimalTransportNetworks.jl/blob/main/examples/example01.jl):

```julia
using OptimalTransportNetworks
# Model Parameters
param = init_parameters(labor_mobility = false, K = 10, gamma = 1, beta = 1, verbose = true,
                        N = 1, tol = 1e-5, cross_good_congestion = false, nu = 1, 
                        duality = false) # Set duality = true/false

# Init network
param, graph = create_graph(param, 11, 11, type = "map") # create a map network of 11x11 nodes located in [0,10]x[0,10]
# note: by default, productivity and population are equalized everywhere

# Customize graph
param[:Zjn] = fill(0.1, param[:J], param[:N]) # set most places to low productivity
Ni = find_node(graph, 6, 6) # Find index of the central node at (6,6)
param[:Zjn][Ni, :] .= 1 # central node more productive

# RESOLUTION
@time results = optimal_network(param, graph) # solve the network

# Plot the optimal network with the optimal transportation plan
plot_graph(graph, results[:Ijk], node_sizes = results[:Cj])
# The size/shade of the nodes reflects the total consumption at each node

```

Gives a total execution time of 8.263035 seconds on my M1 Mac (using the MA57 linear solver) in the primal case and 195.564940 seconds in the dual case (on a small network using a single good).

The opposite is the case in MATLAB, here the [Dual Solution](https://github.com/SebKrantz/OptimalTransportNetworks.jl/blob/main/OptimalTransportNetworkToolbox_v.1.0.4b/Code/objective_duality.m), written as a single vectorized function, is much faster - as expected - than the primal solution.

My diagnosis of this is that the way the problem is written down in JuMP, using ` @expression`’s, seems to cause much unnecessary overhead. Another possibility is the number of Nonzeros in the Hessian (which is 7381 in Julia and 542 in MATLAB, using ADiGator for symbolic autodifferentiation, but since the primal problem has 51062 nonzeros in Julia and is much faster I think this is not the issue).

I am thus seeking ways to implement this dual problem in a more efficient way (having a fast dual solution is crucial for large networks). Ideally I would just translate the MATLAB objective, which takes the price array and returns a scalar (the value of the lagrangian), to Julia, and somehow pass it to Ipopt (including the necessary autodifferentiation). Using JuMP I presume this function would internally create even more symbolic bloat. Having succesfully translated the bulk of the libaray to JuMP, a JuMP solution would be preferable, also given the effective sparse autodifferentiation. So the question I guess is: is there a way to avoid symbolic computations in JuMP to optimize a complex objective without constraints, or any other way to avoid symbolic overheads? Many thanks for your answers.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 7, 2024, 4:18pm UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/2 "2024-05-07T16:18:01Z")

</div>

Disclaimer: I don’t know how to improve things on the JuMP side (luckily Oscar will).  
However, if providing sparse hessians to JuMP from outside can be useful, I’d be thrilled to help you try out the new combo SparseConnectivityTracer.jl + DifferentiationInterface.jl for sparse autodiff with automatic pattern detection.

---

<div class="post-metadata">

**Author:** ![SebKrantz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebkrantz/32/50256_2.png) [@SebKrantz](https://discourse.julialang.org/u/SebKrantz)\
**Post date:** [May 7, 2024, 7:58pm UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/3 "2024-05-07T19:58:08Z")

</div>

Thanks @gdalle, I thinks this would in any case be of interest for performance optimization. Can you provide a minimal example or outline how this can be specified?

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 7, 2024, 9:34pm UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/4 "2024-05-07T21:34:18Z")

</div>

Here’s a minimal example. It’s still experimental so I don’t guarantee optimal performance, but it should scale much better than dense AD in any case. If you run into any trouble, feel free to ping and I or @hill will lend a hand

```julia
using Pkg
Pkg.activate(; temp=true)
Pkg.add(["DifferentiationInterface", "SparseConnectivityTracer", "ForwardDiff", "Zygote"])

using DifferentiationInterface
using SparseConnectivityTracer
using ForwardDiff: ForwardDiff
using Zygote: Zygote

# define your basic second order backend by using a reverse mode inner backend to compute gradients, and a forward mode outer backend to compute pushforwards
inner_backend = AutoZygote()
outer_backend = AutoForwardDiff()
second_order_backend = SecondOrder(outer_backend, inner_backend)

# bring together the three ingredients of sparse AD: sparsity detection + coloring + AD
sparsity_detector = TracerSparsityDetector()
coloring_algorithm = DifferentiationInterface.GreedyColoringAlgorithm()
sparse_second_order_backend = AutoSparse(
    second_order_backend; sparsity_detector, coloring_algorithm
)

# define your function and input
f(x) = sum(abs2, diff(x))
x = collect(1.0:5.0);

# prepare once, reuse the `extras` for different values of x
extras = prepare_hessian(f, sparse_second_order_backend, x)

# compute as many hessians as you want, and enjoy the sparsity
hessian(f, sparse_second_order_backend, x, extras)

```

---

<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:** [May 7, 2024, 9:59pm UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/5 "2024-05-07T21:59:22Z")

</div>

Can you make a minimal reproducible example of just the problematic JuMP model?

It looks like your dual problem has one complicated objective and no constraints?

> <https://github.com/SebKrantz/OptimalTransportNetworks.jl/blob/e6451e80584f94cd3fb7b5fbf9590df6aca35d08/src/models/model_fixed_duality.jl#L24-L81>

Do you have the log from Ipopt? Have you profiled the 195 seconds? Where is the time being spent? In JuMP? Or in Ipopt?

---

<div class="post-metadata">

**Author:** ![SebKrantz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebkrantz/32/50256_2.png) [@SebKrantz](https://discourse.julialang.org/u/SebKrantz)\
**Post date:** [May 8, 2024, 7:50am UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/6 "2024-05-08T07:50:03Z")

</div>

Hi @odow, thanks! So the 195 seconds are for the whole run, which involves solving the model (which computes trade flows, prices and welfare) many times until convergence on the transport network. To perform a single run I can

```julia
using OptimalTransportNetworks
# Model Parameters
param = init_parameters(labor_mobility = false, K = 10, gamma = 1, beta = 1, verbose = true,
                        N = 1, tol = 1e-5, cross_good_congestion = false, nu = 1) 

# Init network
param, graph = create_graph(param, 11, 11, type = "map") # create a map network of 11x11 nodes located in [0,10]x[0,10]
# note: by default, productivity and population are equalized everywhere

# Customize graph
param[:Zjn] = fill(0.1, param[:J], param[:N]) # set most places to low productivity
Ni = find_node(graph, 6, 6) # Find index of the central node at (6,6)
param[:Zjn][Ni, :] .= 1 # central node more productive

# Get model and speficy Ipopt output
param[:duality] = true # Change to see dual/primal models
model, _ = optimal_network(param, graph, verbose = true, return_model = true)

# Just the initial solve
using JuMP
optimize!(model)
# Doing this repeatedly strongly suggests that most time is spent in Ipopt

```

For the dual model I get:

```julia
This is Ipopt version 3.14.4, running with linear solver ma57.

Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 7381

Total number of variables............................: 121
                     variables with only lower bounds: 121
                variables with lower and upper bounds: 0
                     variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 0
        inequality constraints with only lower bounds: 0
   inequality constraints with lower and upper bounds: 0
        inequality constraints with only upper bounds: 0

iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
   0 5.0295207e-01 0.00e+00 9.99e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
   1 5.2422604e-01 0.00e+00 9.31e-03 -1.0 9.86e-02 - 9.91e-01 1.00e+00f 1
   2 4.4131102e-01 0.00e+00 6.97e-03 -2.5 1.12e+00 - 1.00e+00 1.00e+00f 1
   3 -6.9396595e-02 0.00e+00 1.60e-02 -3.8 3.17e+00 - 1.00e+00 1.00e+00f 1
   4 -7.9029511e-02 0.00e+00 4.97e-02 -3.8 2.59e-01 - 1.00e+00 1.00e+00f 1
   5 -8.2938275e-02 0.00e+00 9.04e-03 -3.8 4.41e-02 - 1.00e+00 1.00e+00f 1
   6 -8.2796379e-02 0.00e+00 1.63e-03 -3.8 1.95e-02 - 1.00e+00 1.00e+00f 1
   7 -8.2763843e-02 0.00e+00 5.31e-05 -3.8 3.83e-03 - 1.00e+00 1.00e+00f 1
   8 -8.6539460e-02 0.00e+00 2.51e-04 -5.7 8.97e-02 - 1.00e+00 1.00e+00f 1
   9 -8.6557002e-02 0.00e+00 1.90e-06 -5.7 4.28e-03 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
  10 -8.6558622e-02 0.00e+00 2.38e-07 -8.6 1.73e-03 - 1.00e+00 1.00e+00f 1
  11 -8.6558622e-02 0.00e+00 1.30e-12 -8.6 3.93e-06 - 1.00e+00 1.00e+00f 1

Number of Iterations....: 11

                                   (scaled) (unscaled)
Objective...............: -8.6558622239778449e-02 -8.6558622239778449e-02
Dual infeasibility......: 1.3013525949420977e-12 1.3013525949420977e-12
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5065925976051084e-09 2.5065925976051084e-09
Overall NLP error.......: 2.5065925976051084e-09 2.5065925976051084e-09

Number of objective function evaluations = 12
Number of objective gradient evaluations = 12
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 11
Total seconds in IPOPT = 9.386

EXIT: Optimal Solution Found.

```

And for the primal version

```julia
This is Ipopt version 3.14.4, running with linear solver ma57.

Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 51183
Number of nonzeros in Lagrangian Hessian.............: 51062

Total number of variables............................: 662
                     variables with only lower bounds: 242
                variables with lower and upper bounds: 0
                     variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 242
        inequality constraints with only lower bounds: 0
   inequality constraints with lower and upper bounds: 0
        inequality constraints with only upper bounds: 242

iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
   0 -4.1322314e-02 7.49e-03 8.07e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
   1 -5.2189229e-02 6.69e-03 1.47e+00 -1.7 2.35e-02 - 1.00e+00 4.22e-01h 1
   2 -1.2066922e-01 2.24e-03 1.59e+00 -1.7 1.19e-02 - 1.00e+00 9.27e-01h 1
   3 -1.2688148e-01 2.14e-04 7.87e+00 -1.7 8.53e-03 - 1.00e+00 4.83e-01h 1
   4 -1.3898200e-01 0.00e+00 5.75e-01 -1.7 1.14e-03 - 1.00e+00 1.00e+00h 1
   5 -1.3907991e-01 0.00e+00 1.80e-02 -1.7 2.08e-04 - 1.00e+00 1.00e+00h 1
   6 -1.3907550e-01 0.00e+00 2.71e-03 -3.8 8.38e-07 - 1.00e+00 1.00e+00h 1
   7 -1.3823317e-01 0.00e+00 3.25e-06 -3.8 1.83e-05 - 1.00e+00 1.00e+00h 1
   8 -1.3741170e-01 0.00e+00 2.94e-05 -5.7 1.28e-05 - 1.00e+00 1.00e+00h 1
   9 -1.1489670e-01 0.00e+00 9.17e-04 -5.7 4.97e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
  10 -1.1267338e-01 0.00e+00 6.76e-06 -5.7 6.42e-05 - 1.00e+00 1.00e+00h 1
  11 -1.0276234e-01 0.00e+00 5.31e-04 -8.6 3.97e-04 - 8.31e-01 1.00e+00h 1
  12 -9.2027305e-02 0.00e+00 1.43e-04 -8.6 8.50e-04 - 8.44e-01 1.00e+00h 1
  13 -8.7321303e-02 0.00e+00 2.34e-05 -8.6 6.04e-04 - 1.00e+00 1.00e+00h 1
  14 -8.6660709e-02 0.00e+00 4.18e-07 -8.6 1.05e-04 - 1.00e+00 1.00e+00h 1
  15 -8.6661544e-02 0.00e+00 4.99e-12 -8.6 2.20e-07 - 1.00e+00 1.00e+00h 1

Number of Iterations....: 15

                                   (scaled) (unscaled)
Objective...............: 5.0752466587507044e-04 -8.6661543930583784e-02
Dual infeasibility......: 4.9852007424451824e-12 8.5123979619650069e-10
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5059099014031048e-09 4.2789254514771992e-07
Overall NLP error.......: 2.5059099014031048e-09 4.2789254514771992e-07

Number of objective function evaluations = 16
Number of objective gradient evaluations = 16
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 16
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 16
Number of Lagrangian Hessian evaluations = 15
Total seconds in IPOPT = 0.279

EXIT: Optimal Solution Found.

```

I have also profiled the dual problem using

```julia
using JuMP, Profile
@profile optimize!(model)
Profile.print(combine = true, maxdepth = 1) #, C = true)

```

yielding

```julia
Overhead ╎ [+additional indent] Count File:Line; Function
=========================================================
     ╎1 @Base/array.jl:797; collect(itr::Base.Generator{Base.LogicalIndex{Int64, BitMatrix}, typeof(identity)})
     ╎ 1 @Base/array.jl:823; collect_to_with_first!
    2╎2 @Base/array.jl:2534; filter(f::Core.Compiler.var"#389#392"{Vector{Int64}}, a::Vector{Int64})
    1╎1 @Base/array.jl:1055; push!(a::Vector{Tuple{NonlinearExpr, Int64, NonlinearExpr}}, item::Tuple{NonlinearExpr, Int64, N...
    1╎1 @Base/bitarray.jl:405; falses(dims::Tuple{Int64})
    1╎2 @Base/broadcast.jl:860; materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(...
     ╎ 1 @Base/broadcast.jl:885; copy
     ╎3 @Base/broadcast.jl:868; materialize!(dest::SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false}, bc::Base....
     ╎ 3 @Base/broadcast.jl:871; materialize!
    1╎1 @Base/compiler/abstractinterpretation.jl:1751; abstract_call(interp::Core.Compiler.NativeInterpreter, arginfo::Core.Compiler.ArgInfo, sv::Core....
     ╎1 @Base/compiler/abstractinterpretation.jl:234; abstract_call_gf_by_type(interp::Core.Compiler.NativeInterpreter, f::Any, arginfo::Core.Compiler...
     ╎ 1 @Base/idset.jl:18; isempty
     ╎1 @Base/compiler/abstractinterpretation.jl:2416; typeinf_local(interp::Core.Compiler.NativeInterpreter, frame::Core.Compiler.InferenceState)
     ╎ 1 @Base/idset.jl:18; isempty
    1╎1 @Base/compiler/inferencestate.jl:175; compute_trycatch(code::Vector{Any}, ip::Core.Compiler.BitSet)
    2╎2 @Base/compiler/optimize.jl:342; argextype(x::Any, src::Core.Compiler.IncrementalCompact, sptypes::Vector{Any}, slottypes::Vector...
    1╎1 @Base/compiler/optimize.jl:809; renumber_ir_elements!(body::Vector{Any}, ssachangemap::Vector{Int64}, labelchangemap::Vector{Int...
    1╎1 @Base/compiler/ssair/domtree.jl:252; SNCA!(domtree::Core.Compiler.DomTree, blocks::Vector{Core.Compiler.BasicBlock}, max_pre::Int64)
    1╎1 @Base/compiler/ssair/inlining.jl:948; call_sig(ir::Core.Compiler.IRCode, stmt::Expr)
     ╎1 @Base/compiler/ssair/inlining.jl:321; ir_inline_item!(compact::Core.Compiler.IncrementalCompact, idx::Int64, argexprs::Vector{Any}, li...
    1╎ 1 @Base/compiler/utilities.jl:363; coverage_enabled
    1╎1 @Base/compiler/ssair/inlining.jl:1107; process_simple!(ir::Core.Compiler.IRCode, idx::Int64, state::Core.Compiler.InliningState{Core.Co...
    1╎1 @Base/compiler/ssair/ir.jl:36; basic_blocks_starts(stmts::Vector{Any})
    1╎1 @Base/compiler/ssair/ir.jl:1149; finish_current_bb!(compact::Core.Compiler.IncrementalCompact, active_bb::Int64, old_result_idx::...
   20╎20 @Base/compiler/ssair/ir.jl:460; iterate(it::Core.Compiler.UseRefIterator, #unused#::Nothing)
    2╎2 @Base/compiler/ssair/ir.jl:1325; maybe_erase_unused!(extra_worklist::Vector{Int64}, compact::Core.Compiler.IncrementalCompact, id...
    3╎3 @Base/compiler/ssair/ir.jl:989; process_node!(compact::Core.Compiler.IncrementalCompact, result_idx::Int64, inst::Core.Compiler....
    1╎1 @Base/compiler/ssair/ir.jl:853; process_phinode_values(old_values::Vector{Any}, late_fixup::Vector{Int64}, processed_idx::Int64,...
    3╎3 @Base/compiler/ssair/ir.jl:903; renumber_ssa2!(stmt::Any, ssanums::Vector{Any}, used_ssas::Vector{Int64}, late_fixup::Vector{Int...
   10╎10 @Base/compiler/ssair/ir.jl:410; setindex!(x::Core.Compiler.UseRef, v::Any)
    1╎1 @Base/compiler/ssair/ir.jl:797; setindex!(compact::Core.Compiler.IncrementalCompact, v::Any, idx::Core.SSAValue)
    1╎1 @Base/compiler/ssair/ir.jl:498; ssamap(f::Core.Compiler.var"#359#360"{Vector{Core.SSAValue}, Bool}, stmt::Any)
    4╎4 @Base/compiler/ssair/passes.jl:3; is_known_call(x::Any, func::Any, ir::Core.Compiler.IncrementalCompact)
    1╎1 @Base/compiler/ssair/passes.jl:344; lift_leaves(compact::Core.Compiler.IncrementalCompact, result_t::Any, field::Int64, leaves::Vect...
    2╎2 @Base/compiler/tfuncs.jl:597; typeassert_tfunc(v::Any, t::Any)
    1╎1 @Base/compiler/typeinfer.jl:408; cycle_fix_limited(typ::Any, sv::Core.Compiler.InferenceState)
    1╎1 @Base/compiler/typeinfer.jl:593; visit_slot_load!(sl::Core.SlotNumber, vtypes::Vector{Core.Compiler.VarState}, sv::Core.Compiler....
    2╎2 @Base/compiler/typelattice.jl:359; widenconditional(typ::Any)
    2╎2 @Base/compiler/typelattice.jl:148; ⊑(a::Any, b::Any)
    1╎1 @Base/compiler/utilities.jl:193; specialize_method(method::Method, atype::Any, sparams::Core.SimpleVector; preexisting::Bool, com...
    2╎2 @Base/compiler/utilities.jl:43; anymap(f::Core.Compiler.var"#261#262", a::Vector{Any})
    3╎3 @Base/compiler/utilities.jl:270; find_ssavalue_uses(e::Expr, uses::Vector{Core.Compiler.BitSet}, line::Int64)
    2╎2 @Base/compiler/utilities.jl:235; singleton_type(ft::Any)
     ╎1 @Base/dict.jl:104; Dict{Symbol, Int64}(kv::Base.Generator{Base.Iterators.Enumerate{Vector{Symbol}}, MathOptInterfac...
    1╎ 1 @Base/dict.jl:381; setindex!(h::Dict{Symbol, Int64}, v0::Int64, key::Symbol)
    1╎1 @Base/dict.jl:175; rehash!(h::Dict{Tuple{Int64, Int64}, Nothing}, newsz::Int64)
    1╎1 @Base/expr.jl:37; copy(e::Expr)
    5╎5 @Base/expr.jl:40; copy_exprs(x::Any)
    1╎1 @Base/int.jl:87; +(x::Int64, y::Int64)
    7╎7 @Base/range.jl:321; steprange_last(start::Int64, step::Int64, stop::Int64)
    1╎1 @Base/reflection.jl:1128; may_invoke_generator(method::Method, atype::Any, sparams::Core.SimpleVector)
     ╎1 @Base/sort.jl:571; sort!(v::Vector{Int64}, lo::Int64, hi::Int64, a::Base.Sort.QuickSortAlg, o::Base.Order.ForwardOr...
     ╎ 1 @Base/sort.jl:504; sort!(v::Vector{Int64}, lo::Int64, hi::Int64, #unused#::Base.Sort.InsertionSortAlg, o::Base.Orde...
     ╎48773 @Base/task.jl:484; (::VSCodeServer.var"#64#65")()
     ╎ 48773 @VSCodeServer/src/eval.jl:34; macro expansion
     ╎1 @Base/tuple.jl:434; hash(t::Tuple{DataType, DataType}, h::UInt64)
     ╎ 1 @Base/tuple.jl:434; hash
    5╎5 @Ipopt/src/MOI_wrapper.jl:181; _replace_parameters(model::Ipopt.Optimizer, f::MathOptInterface.VariableIndex)
     ╎2 @Ipopt/src/MOI_wrapper.jl:183; _replace_parameters(model::Ipopt.Optimizer, f::MathOptInterface.VariableIndex)
     ╎ 2 @Base/dict.jl:497; getindex
    2╎2 @Ipopt/src/MOI_wrapper.jl:188; _replace_parameters(model::Ipopt.Optimizer, f::MathOptInterface.ScalarAffineFunction{Float64})
   15╎15 @Ipopt/src/MOI_wrapper.jl:205; _replace_parameters(model::Ipopt.Optimizer, f::MathOptInterface.ScalarNonlinearFunction)
     ╎3 @Ipopt/src/MOI_wrapper.jl:998; get(model::Ipopt.Optimizer, #unused#::MathOptInterface.TerminationStatus)
     ╎ 3 @Base/dict.jl:497; getindex
    1╎1 @Ipopt/src/MOI_wrapper.jl:1110; get(model::Ipopt.Optimizer, attr::MathOptInterface.VariablePrimal, vi::MathOptInterface.Variable...
     ╎1 @JuMP/src/Containers/macro.jl:539; (::OptimalTransportNetworks.var"#90#98"{Matrix{NonlinearExpr}, Vector{VariableRef}, Matrix{Varia...
     ╎ 1 @JuMP/src/macros/@constraint.jl:131; macro expansion
    2╎2 @JuMP/src/nlp_expr.jl:505; check_belongs_to_model(expr::NonlinearExpr, model::Model)
    2╎2 @JuMP/src/nlp_expr.jl:510; check_belongs_to_model(expr::NonlinearExpr, model::Model)
    1╎1 @JuMP/src/optimizer_interface.jl:723; _moi_get_result(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridg...
     ╎1 @JuMP/src/variables.jl:1703; value(v::VariableRef; result::Int64)
    1╎ 1 @JuMP/src/optimizer_interface.jl:752; get(model::Model, attr::MathOptInterface.VariablePrimal, v::VariableRef)
    1╎1 @JuMP/src/variables.jl:338; check_belongs_to_model(v::VariableRef, model::Model)
    1╎1 @MathOptInterface/src/Bridges/bridge_optimizer.jl:1257; get(b::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, attr::MathOptInterface.Var...
     ╎2 @MathOptInterface/src/Bridges/bridge_optimizer.jl:1265; get(b::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer}, attr::MathOptInterface.Var...
    2╎ 2 @Ipopt/src/MOI_wrapper.jl:1110; get(model::Ipopt.Optimizer, attr::MathOptInterface.VariablePrimal, vi::MathOptInterface.Variabl...
     ╎1 @MathOptInterface/src/Nonlinear/ReverseAD/Coloring/Coloring.jl:453; hessian_color_preprocess(edgelist::Set{Tuple{Int64, Int64}}, num_total_var::Int64, seen_idx::Mat...
     ╎ 1 @Base/array.jl:1058; push!
     ╎1 @MathOptInterface/src/Nonlinear/ReverseAD/forward_over_reverse.jl:35; _eval_hessian(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, f::MathOptInterface.Nonlinea...
    1╎ 1 @MathOptInterface/src/Nonlinear/ReverseAD/forward_over_reverse.jl:41; _eval_hessian_inner(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, ex::MathOptInterface.N...
     ╎912 @MathOptInterface/src/Nonlinear/ReverseAD/forward_over_reverse.jl:319; _forward_eval_ϵ(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, ex::MathOptInterface.Nonli...
     ╎ 912 @ForwardDiff/src/dual.jl:240; log
     ╎8 @MathOptInterface/src/Nonlinear/ReverseAD/forward_over_reverse.jl:342; _forward_eval_ϵ(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, ex::MathOptInterface.Nonli...
    8╎ 8 @MathOptInterface/src/Nonlinear/ReverseAD/utils.jl:69; MathOptInterface.Nonlinear.ReverseAD._UnsafeVectorView(x::Vector{Float64}, N::Int64)
     ╎25 @MathOptInterface/src/Nonlinear/ReverseAD/forward_over_reverse.jl:346; _forward_eval_ϵ(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, ex::MathOptInterface.Nonli...
   25╎ 25 @MathOptInterface/src/Nonlinear/ReverseAD/utils.jl:144; MathOptInterface.Nonlinear.ReverseAD._UnsafeLowerTriangularMatrixView(x::Vector{Float64}, N::Int64)
     ╎56 @MathOptInterface/src/Nonlinear/ReverseAD/forward_over_reverse.jl:350; _forward_eval_ϵ(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, ex::MathOptInterface.Nonli...
   56╎ 56 @MathOptInterface/src/Nonlinear/operators.jl:755; eval_multivariate_hessian(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, H::...
     ╎78 @MathOptInterface/src/Nonlinear/ReverseAD/forward_over_reverse.jl:377; _forward_eval_ϵ(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, ex::MathOptInterface.Nonli...
   78╎ 78 @MathOptInterface/src/Nonlinear/operators.jl:571; eval_univariate_hessian(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::Fl...
    2╎2 @MathOptInterface/src/Nonlinear/ReverseAD/forward_over_reverse.jl:130; _hessian_slice_inner(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, ex::MathOptInterface....
     ╎3 @MathOptInterface/src/Nonlinear/ReverseAD/forward_over_reverse.jl:161; _hessian_slice_inner(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, ex::MathOptInterface....
    3╎ 3 @MathOptInterface/src/Nonlinear/ReverseAD/forward_over_reverse.jl:391; _reverse_eval_ϵ(output_ϵ::MathOptInterface.Nonlinear.ReverseAD._UnsafeVectorView{ForwardDiff.Par...
     ╎7 @MathOptInterface/src/Nonlinear/ReverseAD/graph_tools.jl:312; _compute_hessian_sparsity(nodes::Vector{MathOptInterface.Nonlinear.Node}, adj::SparseArrays.Spar...
     ╎ 7 @Base/array.jl:1058; push!
     ╎9 @MathOptInterface/src/Nonlinear/ReverseAD/graph_tools.jl:313; _compute_hessian_sparsity(nodes::Vector{MathOptInterface.Nonlinear.Node}, adj::SparseArrays.Spar...
     ╎ 9 @Base/array.jl:1058; push!
     ╎2 @MathOptInterface/src/Nonlinear/ReverseAD/graph_tools.jl:24; _replace_moi_variables(nodes::Vector{MathOptInterface.Nonlinear.Node}, moi_index_to_consecutive_...
     ╎ 2 @Base/dict.jl:497; getindex
     ╎11 @MathOptInterface/src/Nonlinear/ReverseAD/reverse_mode.jl:178; _forward_eval(f::MathOptInterface.Nonlinear.ReverseAD._FunctionStorage, d::MathOptInterface.Nonl...
   11╎ 11 @NaNMath/src/NaNMath.jl:21; pow
     ╎320 @MathOptInterface/src/Nonlinear/ReverseAD/reverse_mode.jl:181; _forward_eval(f::MathOptInterface.Nonlinear.ReverseAD._FunctionStorage, d::MathOptInterface.Nonl...
  320╎ 320 @NaNMath/src/NaNMath.jl:9; log
     ╎76 @MathOptInterface/src/Nonlinear/ReverseAD/reverse_mode.jl:230; _forward_eval(f::MathOptInterface.Nonlinear.ReverseAD._FunctionStorage, d::MathOptInterface.Nonl...
   76╎ 76 @MathOptInterface/src/Nonlinear/operators.jl:517; eval_univariate_function(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::F...
     ╎99 @MathOptInterface/src/Nonlinear/ReverseAD/reverse_mode.jl:235; _forward_eval(f::MathOptInterface.Nonlinear.ReverseAD._FunctionStorage, d::MathOptInterface.Nonl...
   99╎ 99 @MathOptInterface/src/Nonlinear/operators.jl:544; eval_univariate_gradient(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::F...
     ╎1 @MathOptInterface/src/Nonlinear/ReverseAD/reverse_mode.jl:32; _reverse_mode(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, x::Vector{Float64})
     ╎ 1 @Base/abstractarray.jl:2626; ==(A::Vector{Float64}, B::Vector{Float64})
   67╎67 @MathOptInterface/src/Nonlinear/operators.jl:607; eval_multivariate_function(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x:...
   79╎79 @MathOptInterface/src/Nonlinear/operators.jl:658; eval_multivariate_gradient(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, g:...
  851╎851 @MathOptInterface/src/Nonlinear/operators.jl:755; eval_multivariate_hessian(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, H::...
  185╎185 @MathOptInterface/src/Nonlinear/operators.jl:517; eval_univariate_function(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::F...
     ╎62 @MathOptInterface/src/Nonlinear/operators.jl:522; eval_univariate_function(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::F...
   62╎ 62 @Base/dict.jl:496; getindex(h::Dict{Symbol, Int64}, key::Symbol)
  159╎159 @MathOptInterface/src/Nonlinear/operators.jl:544; eval_univariate_gradient(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::F...
     ╎95 @MathOptInterface/src/Nonlinear/operators.jl:549; eval_univariate_gradient(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::F...
   95╎ 95 @Base/dict.jl:496; getindex(h::Dict{Symbol, Int64}, key::Symbol)
  789╎789 @MathOptInterface/src/Nonlinear/operators.jl:571; eval_univariate_hessian(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::Fl...
     ╎41 @MathOptInterface/src/Nonlinear/operators.jl:576; eval_univariate_hessian(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::Fl...
   41╎ 41 @Base/dict.jl:496; getindex(h::Dict{Symbol, Int64}, key::Symbol)
    9╎9 @MathOptInterface/src/Nonlinear/parse.jl:59; _get_node_type(data::MathOptInterface.Nonlinear.Model, x::MathOptInterface.ScalarNonlinearFunction)
   16╎16 @MathOptInterface/src/Nonlinear/parse.jl:79; _parse_without_recursion_inner(stack::Vector{Tuple{Int64, Any}}, data::MathOptInterface.Nonlinea...
     ╎8 @MathOptInterface/src/Nonlinear/parse.jl:81; _parse_without_recursion_inner(stack::Vector{Tuple{Int64, Any}}, data::MathOptInterface.Nonlinea...
     ╎ 8 @Base/array.jl:1058; push!
     ╎13 @MathOptInterface/src/Nonlinear/parse.jl:86; _parse_without_recursion_inner(stack::Vector{Tuple{Int64, Any}}, data::MathOptInterface.Nonlinea...
     ╎ 13 @Base/array.jl:1058; push!
     ╎3 @MathOptInterface/src/Nonlinear/parse.jl:261; parse_expression(#unused#::MathOptInterface.Nonlinear.Model, expr::MathOptInterface.Nonlinear.Ex...
     ╎ 3 @Base/array.jl:1058; push!
     ╎4 @MathOptInterface/src/Nonlinear/parse.jl:364; parse_expression(#unused#::MathOptInterface.Nonlinear.Model, expr::MathOptInterface.Nonlinear.Ex...
     ╎ 4 @Base/array.jl:1058; push!
     ╎8 @MathOptInterface/src/Nonlinear/parse.jl:365; parse_expression(#unused#::MathOptInterface.Nonlinear.Model, expr::MathOptInterface.Nonlinear.Ex...
     ╎ 8 @Base/array.jl:1058; push!
    1╎1 @MathOptInterface/src/Nonlinear/parse.jl:366; parse_expression(#unused#::MathOptInterface.Nonlinear.Model, expr::MathOptInterface.Nonlinear.Ex...
     ╎1 @MathOptInterface/src/Nonlinear/parse.jl:375; parse_expression(#unused#::MathOptInterface.Nonlinear.Model, expr::MathOptInterface.Nonlinear.Ex...
     ╎ 1 @Base/array.jl:1058; push!
     ╎1 @MathOptInterface/src/Utilities/cachingoptimizer.jl:876; _throw_if_get_attribute_not_allowed(model::MathOptInterface.Utilities.CachingOptimizer{MathOptIn...
     ╎ 1 @Base/essentials.jl:788; isempty
    1╎1 @MathOptInterface/src/Utilities/cachingoptimizer.jl:584; _replace_constraint_function_or_set(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterf...
    1╎1 @MathOptInterface/src/Utilities/cachingoptimizer.jl:898; get(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimi...
     ╎1 @MathOptInterface/src/Utilities/cachingoptimizer.jl:905; get(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimi...
     ╎ 1 @MathOptInterface/src/Bridges/bridge_optimizer.jl:902; get
    1╎1 @MathOptInterface/src/Utilities/cachingoptimizer.jl:920; get(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimi...
     ╎4 @MathOptInterface/src/Utilities/cachingoptimizer.jl:932; get(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimi...
     ╎ 4 @MathOptInterface/src/Utilities/copy/index_map.jl:64; getindex
    1╎1 @MathOptInterface/src/Utilities/cachingoptimizer.jl:780; set(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{...
     ╎2 @MathOptInterface/src/Utilities/functions.jl:952; is_canonical(f::MathOptInterface.ScalarNonlinearFunction)
     ╎ 2 @Base/array.jl:1064; push!
    2╎2 @MathOptInterface/src/Utilities/functions.jl:246; map_indices(index_map::Base.Fix1{typeof(getindex), MathOptInterface.Utilities.IndexMap}, vi::Mat...
     ╎5 @MathOptInterface/src/Utilities/functions.jl:247; map_indices(index_map::Base.Fix1{typeof(getindex), MathOptInterface.Utilities.IndexMap}, vi::Mat...
     ╎ 5 @Base/operators.jl:1096; Fix1
    1╎1 @MathOptInterface/src/Utilities/functions.jl:254; map_indices(index_map::Base.Fix1{typeof(getindex), MathOptInterface.Utilities.IndexMap}, x::Vect...
     ╎8 @MathOptInterface/src/Utilities/functions.jl:365; map_indices(index_map::Base.Fix1{typeof(getindex), MathOptInterface.Utilities.IndexMap}, f::Math...
     ╎ 8 @Base/array.jl:1058; push!
     ╎1 @MathOptInterface/src/Utilities/universalfallback.jl:463; get(uf::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}},...
     ╎ 1 @MathOptInterface/src/Utilities/model.jl:378; get
    1╎1 @MathOptInterface/src/functions.jl:1221; convert(F::Type{MathOptInterface.ScalarNonlinearFunction}, f::MathOptInterface.ScalarAffineFunct...
    1╎1 @OrderedCollections/src/ordered_dict.jl:235; ht_keyindex(h::OrderedCollections.OrderedDict{MathOptInterface.Nonlinear.ConstraintIndex, MathOp...
     ╎1 @OrderedCollections/src/ordered_dict.jl:240; ht_keyindex(h::OrderedCollections.OrderedDict{MathOptInterface.Nonlinear.ConstraintIndex, MathOp...
     ╎ 1 @OrderedCollections/src/dict_support.jl:6; hashindex
Total snapshots: 66771. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task

```

---

<div class="post-metadata">

**Author:** ![SebKrantz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebkrantz/32/50256_2.png) [@SebKrantz](https://discourse.julialang.org/u/SebKrantz)\
**Post date:** [May 8, 2024, 5:02pm UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/7 "2024-05-08T17:02:18Z")

</div>

Unfortunately this doesn’t work with my rather complex function. I’ve tried a host of libraries now (including `Nonconvex.jl` and `ADNLPModels.jl`), and the only thing that does work is

```julia
using Optim
Optim.optimize(objective, lb, ub, x0, Fminbox(ConjugateGradient()))

```

Takes about 5 seconds though which is only slightly faster than the current JuMP solution. It would be really great if for this dual case I could have a way to solve the problem with `Ipopt` (which takes much less than 1 sec for the 121 prices here). Is there another (safer) way to do sparse autodifferentiation? In that case I may be able to set the problem up using the Ipopt C API.

FYI the objective is

```julia
using OptimalTransportNetworks
import OptimalTransportNetworks as otn

# Model Parameters
param = init_parameters(labor_mobility = false, K = 10, gamma = 1, beta = 1, verbose = true,
                        N = 1, tol = 1e-5, cross_good_congestion = false, nu = 1) 

# Init network
param, graph = create_graph(param, 11, 11, type = "map") # create a map network of 11x11 nodes located in [0,10]x[0,10]

# Customize graph
param[:Zjn] = fill(0.1, param[:J], param[:N]) # set most places to low productivity
Ni = find_node(graph, 6, 6) # Find index of the central node at (6,6)
param[:Zjn][Ni, :] .= 1 # central node more productive

# Initial infrastructure
edges = otn.represent_edges(otn.dict_to_namedtuple(graph))
I0 = Float64.(graph[:adjacency])
I0 *= param[:K] / sum(graph[:delta_i] .* I0) 

# Auxiliary data
auxdata = otn.create_auxdata(otn.dict_to_namedtuple(param), otn.dict_to_namedtuple(graph), edges, I0)

# Now the objective
function objective_duality(x, auxdata)

    # Extract parameters
    param = auxdata.param
    graph = auxdata.graph
    kappa_ex = auxdata.kappa_ex
    A = auxdata.edges.A
    Apos = auxdata.edges.Apos
    Aneg = auxdata.edges.Aneg
    es = auxdata.edges.edge_start
    ee = auxdata.edges.edge_end
    omegaj = param.omegaj
    Lj = param.Lj
    alpha = param.alpha
    sigma = param.sigma
    rho = param.rho
    hj = param.hj
    hj1malpha = (hj / (1-alpha)) .^ (1-alpha)
    Zjn = param.Zjn
    beta = param.beta
    a = param.a
    onesN = ones(param.N)

    # Extract price vector
    Pjn = reshape(x, (graph.J, param.N))

    # Calculate consumption
    cj = zeros(graph.J)
    for j=1:graph.J
       cj[j] = alpha * (sum(Pjn[j, n]^(1-sigma) for n=1:param.N)^(1/(1-sigma)) / omegaj[j])^(-1/(1+alpha*(rho-1))) * hj1malpha[j]^(-((rho-1)/(1+alpha*(rho-1))))
    end
    uj = zeros(graph.J)
    for j=1:graph.J
       uj[j] = ((cj[j]/alpha)^alpha * hj1malpha[j])^(1-rho)/(1-rho)
    end
    zeta = zeros(graph.J)
    for j=1:graph.J
       zeta[j] = omegaj[j] * uj[j] * (1-rho) * alpha / cj[j]
    end
    cjn = zeros(graph.J, param.N)
    for j=1:graph.J, n=1:param.N
        cjn[j, n] = (Pjn[j, n] / zeta[j])^(-sigma) * cj[j]
    end

    # Calculate Q, Qin_direct (which is the flow in the direction of the edge)
    # and Qin_indirect (flow in edge opposite direction)
    Qin_direct = zeros(graph.ndeg, param.N)
    Qin_indirect = zeros(graph.ndeg, param.N)
    for i=1:graph.ndeg, n=1:param.N
        Qin_direct[i,n] = max(1/(1+beta) * kappa_ex[i] * (Pjn[ee[i],n]/Pjn[es[i],n] - 1), 0)^(1/beta)
    end
    for i=1:graph.ndeg, n=1:param.N
        Qin_indirect[i,n] = max(1/(1+beta) * kappa_ex[i] * (Pjn[es[i],n]/Pjn[ee[i],n] - 1), 0)^(1/beta)
    end

    # Calculate labor allocation
    PZ = (Pjn .* Zjn) .^ (1/(1-a))
    Ljn = zeros(graph.J, param.N)
    for j=1:graph.J, n=1:param.N
        Ljn[j, n] = ifelse(Zjn[j,n] == 0, 0, PZ[j, n] / sum(PZ[j, n] for n=1:param.N) * Lj[j])
    end

    # Calculate production Yjn
    Yjn = Zjn .* Ljn .^ a

    # Create flow constraint
    B_direct = ((Qin_direct .^ (1+beta)) * onesN) ./ kappa_ex
    B_indirect = ((Qin_indirect .^ (1+beta)) * onesN) ./ kappa_ex
    consraw = cjn .* repeat(Lj, 1, param.N) + A * Qin_direct - A * Qin_indirect - Yjn + Apos * B_direct + Aneg * B_indirect
    cons = sum(Pjn .* consraw, dims = 2)

    # Lagrangian
    f = sum(omegaj .* Lj .* uj - cons)
    return f
end

# Simplified version
function objective(x)
    return objective_duality(x, auxdata)
end

# Initial values
v1 = range(1, 2, length=graph[:J])
v2 = param[:N] == 1 ? 2.0 : range(1, 2, length=param[:N])
x0 = collect(v1 * v2') 
lb = repeat([1e-6], length(x0)) 
ub = repeat([Inf], length(x0))

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 8, 2024, 5:21pm UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/8 "2024-05-08T17:21:07Z")

</div>

> [@SebKrantz](#):
>
> Unfortunately this doesn’t work with my rather complex function.

Can you specify exactly what doesn’t work? Does it just compile forever?

---

<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:** [May 8, 2024, 8:24pm UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/9 "2024-05-08T20:24:40Z")

</div>

9 seconds seems far too long for such a small problem. Let me take a look. This is likely a performance “bug” somewhere in the AD code.

---

<div class="post-metadata">

**Author:** ![SebKrantz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebkrantz/32/50256_2.png) [@SebKrantz](https://discourse.julialang.org/u/SebKrantz)\
**Post date:** [May 8, 2024, 8:58pm UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/10 "2024-05-08T20:58:38Z")

</div>

So the error I get comes from

```julia
extras = prepare_hessian(objective, sparse_second_order_backend, x0)

```

it is:

```julia
ERROR: TypeError: in typeassert, expected Float64, got a value of type HessianTracer{BitSet}
Stacktrace:
 [1] setindex!(A::Vector{Float64}, x::HessianTracer{BitSet}, i1::Int64)
   @ Base ./array.jl:966
 [2] objective_duality(x::Vector{HessianTracer{BitSet}}, auxdata::NamedTuple{(:param, :graph, :edges, :kappa, :kappa_ex), Tuple{NamedTuple{(:alpha, :kappa_min, :gamma, :annealing, :duality, :rho, :F, :uprimeinv, :cong, :usecond, :min_iter, :nu, :K, :warm_start, :sigma, :N, :verbose, :tol, :max_iter, :Fprime, :omegaj, :u, :Zjn, :mobility, :a, :m, :beta, :hj, :J, :uprime, :Lj, :Hj, :optimizer_attr), Tuple{Float64, Float64, Int64, Bool, Bool, Int64, OptimalTransportNetworks.var"#24#30", OptimalTransportNetworks.var"#23#29"{Float64, Int64}, Bool, OptimalTransportNetworks.var"#22#28"{Float64, Int64}, Int64, Int64, Int64, Bool, Int64, Int64, Bool, Float64, Int64, OptimalTransportNetworks.var"#25#31", Vector{Float64}, OptimalTransportNetworks.var"#20#26"{Float64, Int64}, Matrix{Float64}, Bool, Float64, Vector{Float64}, Int64, Vector{Float64}, Int64, OptimalTransportNetworks.var"#21#27"{Float64, Int64}, Vector{Float64}, Vector{Float64}, Dict{Symbol, String}}}, NamedTuple{(:J, :y, :nodes, :delta_i, :adjacency, :ndeg, :delta_tau, :x), Tuple{Int64, Vector{Int64}, Vector{Vector{Int64}}, Matrix{Float64}, BitMatrix, Int64, Matrix{Float64}, Vector{Int64}}}, NamedTuple{(:A, :Apos, :Aneg, :edge_start, :edge_end), Tuple{Matrix{Int64}, Matrix{Int64}, Matrix{Int64}, Vector{Int64}, Vector{Int64}}}, Matrix{Float64}, Vector{Float64}}})
   @ Main ~/Documents/Julia/OptimalTransportNetworks.jl/misc/experimental/dual_solution_manual.jl:66
 [3] objective(x::Vector{HessianTracer{BitSet}})
   @ Main ~/Documents/Julia/OptimalTransportNetworks.jl/misc/experimental/dual_solution_manual.jl:114
 [4] trace_function(#unused#::Type{HessianTracer{BitSet}}, f::typeof(objective), x::Vector{Float64})
   @ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/KDSIC/src/pattern.jl:22
 [5] hessian_pattern(f::Function, x::Vector{Float64}, ::Type{BitSet})
   @ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/KDSIC/src/pattern.jl:192
 [6] hessian_sparsity
   @ ~/.julia/packages/SparseConnectivityTracer/KDSIC/src/adtypes.jl:44 [inlined]
 [7] prepare_hessian(f::typeof(objective), backend::AutoSparse{SecondOrder{AutoForwardDiff{nothing, Nothing}, AutoZygote}, TracerSparsityDetector{BitSet}, DifferentiationInterface.GreedyColoringAlgorithm}, x::Vector{Float64})
   @ DifferentiationInterface ~/.julia/packages/DifferentiationInterface/zNsrZ/src/sparse/hessian.jl:13
 [8] top-level scope
   @ ~/Documents/Julia/OptimalTransportNetworks.jl/misc/experimental/dual_solution_manual.jl:205

```

and highlighting the first loop in the objective

```julia
 cj = zeros(graph.J)
  for j=1:graph.J
     cj[j] = alpha * (sum(Pjn[j, n]^(1-sigma) for n=1:param.N)^(1/(1-sigma)) / omegaj[j])^(-1/(1+alpha*(rho-1))) * hj1malpha[j]^(-((rho-1)/(1+alpha*(rho-1))))
  end

```

I presume that indexing assignments are simply not allowed in an objective. Potentially I could rewrite this without indexing assignments by vectorizing everything, but this would complicate certain expressions and potentially yield efficiency losses. Overall my impression is that this solution, though neat, seems a bit complicated (also dependency wise for the package), and the sparse hessian is of course just one part of it, I need a suitable optimizer (ideally Ipopt) to go with it. I appreciate the effort though! If you can get it to run and it is indeed significantly faster than the simple `Optim.optimize(objective, lb, ub, x0, Fminbox(ConjugateGradient()))` solution, I’ll gladly adopt it.

---

<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:** [May 9, 2024, 1:43am UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/11 "2024-05-09T01:43:30Z")

</div>

Let’s chat over here: [Debug performance issue in Nonlinear submodule · Issue #2496 · jump-dev/MathOptInterface.jl · GitHub](https://github.com/jump-dev/MathOptInterface.jl/issues/2496)

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 9, 2024, 5:32am UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/12 "2024-05-09T05:32:09Z")

</div>

> [@SebKrantz](#):
>
> `ERROR: TypeError: in typeassert, expected Float64, got a value of type HessianTracer{BitSet}`

That’s not because of the indexing, that’s because your code does not allow non-floating point number types. The way SparseConnectivityTracer.jl works (like ForwardDiff.jl) is operator overloading: replace the usual number formats with something more fancy to trace dependencies. This alternative number type is called `HessianTracer`.

When you create `zeros(graph.J)`, it falls back on `zeros(Float64, graph.J)`. Then, when you give a new value to `cj[j]`, it has to be cast to a `Float64` to fit into the array. Except it’s a `HessianTracer`.  
The solution is to try and make every quantity that you allocate in the algorithm have the same `eltype` as the initial vector.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 9, 2024, 5:34am UTC](https://discourse.julialang.org/t/efficient-dual-problems-in-jump-or-avoiding-symbolic-overheads-with-complex-objectives/113953/13 "2024-05-09T05:34:53Z")

</div>

> [@SebKrantz](#):
>
> If you can get it to run and it is indeed significantly faster than the simple `Optim.optimize(objective, lb, ub, x0, Fminbox(ConjugateGradient()))` solution, I’ll gladly adopt it.

I’ll give it a shot! Scalar indexing is not a bug but it will make things slow though, since you need a reverse mode backend in addition to forward mode for the hessian autodiff. I’ll see if I can get Enzyme to work instead of Zygote