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

I am developing 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 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 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:

``````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, 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.

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.

1 Like

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?

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

``````using Pkg
Pkg.activate(; temp=true)

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)
``````
1 Like

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?

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

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

``````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:

``````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

``````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

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

yielding

``````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})
β3     @Base/broadcast.jl:868; materialize!(dest::SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false}, bc::Base....
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)
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 @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!
β 912   @ForwardDiff/src/dual.jl:240; log
56β 56    @MathOptInterface/src/Nonlinear/operators.jl:755; eval_multivariate_hessian(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, H::...
78β 78    @MathOptInterface/src/Nonlinear/operators.jl:571; eval_univariate_hessian(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::Fl...
β 7     @Base/array.jl:1058; push!
β 9     @Base/array.jl:1058; push!
β 2     @Base/dict.jl:497; getindex
11β 11    @NaNMath/src/NaNMath.jl:21; pow
320β 320   @NaNMath/src/NaNMath.jl:9; log
76β 76    @MathOptInterface/src/Nonlinear/operators.jl:517; eval_univariate_function(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::F...
99β 99    @MathOptInterface/src/Nonlinear/operators.jl:544; eval_univariate_gradient(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, x::F...
β 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:...
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)
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
``````
1 Like

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

``````using Optim
``````

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

``````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 *= 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))
``````

Can you specify exactly what doesnβt work? Does it just compile forever?

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.

1 Like

So the error I get comes from

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

it is:

``````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
[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

`````` 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.

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.

1 Like

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

1 Like