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

Specific Domains Optimization (Mathematical)
,
1

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.

2

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
3

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?

4

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

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?

6

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

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

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))
8

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

9

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
10

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
   @ ~/.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

 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.

11

Let’s chat over here: Debug performance issue in Nonlinear submodule Β· Issue #2496 Β· jump-dev/MathOptInterface.jl Β· GitHub

12

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
13

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