Automatic differentiation and Trapz

Specific Domains Optimization (Mathematical)
1

It seems that Nonconvex is not able to deal with the trapz function, but maybe, I’ve done something wrong. The code (that is a nonsense toy code) :

using Nonconvex
Nonconvex.@load Ipopt
using Trapz

function obj(hV,N,M)
    h=reshape(hV,(N,M))
    vx=range(0,1,length=N)
    vy=range(0,1,length=M)
    output = trapz((vx,vy),h)
    return tr(h)
end

function cstr(hV,N,M,b)
    h=reshape(hV,(N,M))
    return det(h)-b
end

N=2
M=2

model = Model(x->obj(x,N,M))
addvar!(model, vec(zeros(N,M)), 4.0*vec(ones(N,M)))
add_ineq_constraint!(model, x -> cstr(x, N, M, 2.0))

alg = IpoptAlg()
options = IpoptOptions(first_order = true, tol = 1e-7)
r = optimize(model, alg, vec(rand(N,M)), options = options)
2

Any idea ?

3

Cc @mohamed82008

What’s the error?

4 
ERROR: LoadError: Mutating arrays is not supported -- called setindex!(::Array{Float64, 0}, _...)
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] (::Zygote.var"#437#438"{Array{Float64, 0}})(#unused#::Nothing)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/lib/array.jl:71
  [3] (::Zygote.var"#2337#back#439"{Zygote.var"#437#438"{Array{Float64, 0}}})(Δ::Nothing)
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
  [4] macro expansion
    @ ~/.julia/packages/Trapz/UoNri/src/kernels.jl:42 [inlined]
  [5] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [6] Pullback
    @ ~/.julia/packages/Trapz/UoNri/src/kernels.jl:41 [inlined]
  [7] macro expansion
    @ ./none:0 [inlined]
  [8] Pullback
    @ ./none:0 [inlined]
  [9] Pullback
    @ ~/.julia/packages/Trapz/UoNri/src/Trapz.jl:7 [inlined]
 [10] (::typeof(∂(trapz)))(Δ::Nothing)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [11] Pullback
    @ ~/Nextcloud/Dauphine/OptiForm/Julia/Nonconvex/testTrapz.jl:9 [inlined]
 [12] (::typeof(∂(obj)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [13] Pullback
    @ ~/Nextcloud/Dauphine/OptiForm/Julia/Nonconvex/testTrapz.jl:21 [inlined]
 [14] (::typeof(∂(#1)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [15] #208
    @ ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:207 [inlined]
 [16] (::Zygote.var"#1750#back#210"{Zygote.var"#208#209"{Tuple{Tuple{Nothing}}, typeof(∂(#1))}})(Δ::Float64)
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [17] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/functions/functions.jl:170 [inlined]
 [18] (::typeof(∂(#_#8)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [19] (::Zygote.var"#208#209"{Tuple{Tuple{Nothing, Nothing}, Tuple{Nothing}}, typeof(∂(#_#8))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:207
 [20] #1750#back
    @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67 [inlined]
 [21] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/functions/functions.jl:170 [inlined]
 [22] (::typeof(∂(λ)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [23] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/models/vec_model.jl:79 [inlined]
 [24] (::typeof(∂(λ)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [25] #208
    @ ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:207 [inlined]
 [26] (::Zygote.var"#1750#back#210"{Zygote.var"#208#209"{Tuple{Tuple{Nothing}}, typeof(∂(λ))}})(Δ::Float64)
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [27] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/functions/functions.jl:170 [inlined]
 [28] (::typeof(∂(#_#8)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [29] (::Zygote.var"#208#209"{Tuple{Tuple{Nothing, Nothing}, Tuple{Nothing}}, typeof(∂(#_#8))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:207
 [30] #1750#back
    @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67 [inlined]
 [31] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/functions/functions.jl:170 [inlined]
 [32] (::typeof(∂(λ)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [33] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/functions/counting_function.jl:9 [inlined]
 [34] (::typeof(∂(λ)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [35] (::Zygote.var"#52#53"{typeof(∂(λ))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface.jl:41
 [36] gradient(f::Function, args::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface.jl:76
 [37] (::NonconvexIpopt.var"#eval_grad_f#13"{NonconvexCore.CountingFunction{NonconvexCore.Objective{NonconvexCore.var"#133#140"{Model{Vector{Any}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Base.RefValue{Float64}, Set{Symbol}}}})(x::Vector{Float64}, grad_f::Vector{Float64})
    @ NonconvexIpopt ~/.julia/packages/NonconvexIpopt/jkxtN/src/ipopt.jl:135
 [38] _Eval_Grad_F_CB(n::Int32, x_ptr::Ptr{Float64}, #unused#::Int32, grad_f::Ptr{Float64}, user_data::Ptr{Nothing})
    @ Ipopt ~/.julia/packages/Ipopt/M2QE8/src/C_wrapper.jl:49
 [39] IpoptSolve(prob::Ipopt.IpoptProblem)
    @ Ipopt ~/.julia/packages/Ipopt/M2QE8/src/C_wrapper.jl:433
 [40] optimize!(workspace::NonconvexIpopt.IpoptWorkspace{NonconvexCore.VecModel{NonconvexCore.Objective{NonconvexCore.var"#133#140"{Model{Vector{Any}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Base.RefValue{Float64}, Set{Symbol}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.EqConstraint}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.IneqConstraint{NonconvexCore.var"#137#144"{NonconvexCore.IneqConstraint{FunctionWrapper{var"#3#4", Int64}, Float64, Int64, Set{Symbol}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Float64, Int64, Set{Symbol}}}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.SDConstraint}}, Vector{Float64}, Vector{Float64}, Vector{Float64}}, Ipopt.IpoptProblem, Vector{Float64}, IpoptOptions{NamedTuple{(:tol, :hessian_approximation, :jac_c_constant, :jac_d_constant), Tuple{Float64, String, String, String}}}, Base.RefValue{Int64}})
    @ NonconvexIpopt ~/.julia/packages/NonconvexIpopt/jkxtN/src/ipopt.jl:52
 [41] #optimize#131
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/models/vec_model.jl:63 [inlined]
 [42] optimize(::Model{Vector{Any}}, ::IpoptAlg, ::Vector{Float64}; kwargs::Base.Pairs{Symbol, IpoptOptions{NamedTuple{(:tol, :hessian_approximation, :jac_c_constant, :jac_d_constant), Tuple{Float64, String, String, String}}}, Tuple{Symbol}, NamedTuple{(:options,), Tuple{IpoptOptions{NamedTuple{(:tol, :hessian_approximation, :jac_c_constant, :jac_d_constant), Tuple{Float64, String, String, String}}}}}})
    @ NonconvexCore ~/.julia/packages/NonconvexCore/uZAVo/src/common.jl:233
 [43] top-level scope
    @ ~/Nextcloud/Dauphine/OptiForm/Julia/Nonconvex/testTrapz.jl:27
in expression starting at /home/mc/Nextcloud/Dauphine/OptiForm/Julia/Nonconvex/testTrapz.jl:27
5

The error is saying (unhelpfully) that Trapz isn’t compatible with automatic differentiation.

6

Indeed :slight_smile: that’s why I posted this message.

7

How many decision variables are there? If not many, you can try to use ForwardDiff for the function. In the next release (in a day or 2), you can do ff = forwarddiffy(f, x...) to tell Nonconvex to use ForwardDiff for function ff(x...).

8

You could pretty easily define a manual differentiation rule (rrule) for trapz via ChainRules.jl. (This would make a good PR for the Trapz package too.)

Or you could use Integrals.jl instead of Trapz, since that provides differentiable integration routines.

3 Likes
9

Hi Maxime :slight_smile: To complete the responses here, the issue is that Zygote does not support mutation of arrays (which is anyway a pretty tricky thing to do in an AD engine). You can just precompute weights and replace trapz(x) with sum(weights .* x): you’ll need to compute the weights if you want to write a rrule as @stevengj suggested anyway. Or just don’t use Trapz and use a simple uniform integration scheme.

1 Like
10

Thanks everyone. I’m going to test the different solutions.

11

In Nonconvex 2.0, you can now use abstractdiffy to change the AD package used for a function or the entire model. See Using other AD packages · Nonconvex.jl for docs.

2 Likes