# IPNewton and ForwardDiff with preallocated Arrays

Hi !

I am new here so please let me know if I post at the right place. I could not find that issue on other threads.

I need to solve a constrained optimization problem, and I’d like to use Optim.optimize with the IPNewton method. The cost function makes use of many arrays so I would like to preallocate all of them. I also want to use automatic differentiation.

The preallocated arrays need to have the right type (which seems to be

``````ForwardDiff.Dual{ForwardDiff.Tag{var"#J#409"{Array{Float64,1}},Float64},Float64,5}
``````

in my case), but the optimize function fails.

Here is a toy script that gives exactly the same error as mine. It is just looking for square roots of an array of Floats.

``````using Optim, ForwardDiff

function test()
data = [1.0, 2.0, 3.0, 4.0, 5.0]
U = fill(1.0, length(data))
function J(U)
#makes use of a preallocated array called aux
for i in 1:length(U)
aux[i] = U[i]^2
end
sum((aux .- data).^2)
end
h! = (hess, x) -> ForwardDiff.hessian!(hess, J, x)
type = eltype(cfg)
aux = ones(type, length(data)) #here it is
df = TwiceDifferentiable(J, g!, h!, U)
#lower = fill(zero(type), length(data))
lower = fill(zero(eltype(U)), length(data))
upper = fill(Inf, length(data))
dfc = TwiceDifferentiableConstraints(lower, upper)
opt = optimize(df, dfc, U, IPNewton())
return opt
end

test()
``````

And here is the error I get :

``````ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#J#419"{Array{Float64,1}},Float64},Float64,5})
Closest candidates are:
Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
Float64(::T) where T<:Number at boot.jl:716
Float64(::Float32) at float.jl:255
...
Stacktrace:
[1] convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#J#419"{Array{Float64,1}},Float64},Float64,5}) at ./number.jl:7
[2] setproperty!(::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Symbol, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#J#419"{Array{Float64,1}},Float64},Float64,5}) at ./Base.jl:34
[5] initial_state(::IPNewton{typeof(Optim.backtrack_constrained_grad),Symbol}, ::Optim.Options{Float64,Nothing}, ::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::TwiceDifferentiableConstraints{NLSolversBase.var"#131#134",NLSolversBase.var"#132#135",NLSolversBase.var"#133#136",Float64}, ::Array{Float64,1}) at /Users/sylvainmoi/.julia/packages/Optim/Yd5Zq/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl:125
[6] optimize(::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::TwiceDifferentiableConstraints{NLSolversBase.var"#131#134",NLSolversBase.var"#132#135",NLSolversBase.var"#133#136",Float64}, ::Array{Float64,1}, ::IPNewton{typeof(Optim.backtrack_constrained_grad),Symbol}, ::Optim.Options{Float64,Nothing}) at /Users/sylvainmoi/.julia/packages/Optim/Yd5Zq/src/multivariate/solvers/constrained/ipnewton/interior.jl:228 (repeats 2 times)
[7] test() at /Users/sylvainmoi/Documents/These/Optim_Julia/cas_ecole.jl:24
[8] top-level scope at none:1
``````

It has to do with the function Optim.value_gradient!. If I change

``````aux = ones(type, length(data))
``````

to

``````aux = ones(Real, length(data))
``````

it works but it’s useless, because each

``````aux[i] = U[i]^2
``````

in the for-loop needs then to allocate every time.

Thanks a lot if you can have a look at that.

Cheers,

You have the right idea of using `eltype` when you want to pre-allocate things. But I think the problem is that `aux` needs to be a different type for different calls. For example, calling `J` with the initial vector to get the objective value means that you want `aux` to be a `Vector{Float64}`, but for the AD calls you’ll want it to be a `Vector{Dual{T,V,N}} where{T,V,N}`, for example. The reason it works for `Real`, an abstract type, is because both `Float64` and `Dual` are subtypes of `Real`.

I could be incorrect, but it may be easiest to refactor your problem a bit. Presumably your real application actually does need `aux`, but in this case you don’t actually need it, and a similar idea of this refactor and moving of `J` outside `test` might also work for your real problem.

``````using Optim, ForwardDiff

# Probably cleaner to make this a separate function that takes more arguments
# and create a closure in test.
_J(U, data) = sum(z->abs2(z[1]^2-z[2]), zip(U, data))

function test()
data = [1.0, 2.0, 3.0, 4.0, 5.0]
U    = ones(length(data))
obj  = x -> _J(x, data)
h! = (hess, x) -> ForwardDiff.hessian!(hess, obj, x)
df = TwiceDifferentiable(obj, g!, h!, U)
lower = fill(0.0, length(data))
upper = fill(Inf, length(data))
dfc = TwiceDifferentiableConstraints(lower, upper)
opt = optimize(df, dfc, U, IPNewton())
return opt
end

result = test()
``````
1 Like

I see indeed a type problem : inside the function test() above, I tried to add

``````  ...
aux = ones(type, length(data))

println(J(U)) #works
println(J(aux)) #works
I’m a bit confused because I can read in the ForwardDiff doc that the `x` argument of `ForwardDiff.gradient` must be of type `x::AbstractArray` and I assumed that `Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#J#354"{Array{Float64,1}},Float64},Float64,5},1}` (ie `aux` type) would match with `AbstractArray`.