# IPNewton and ForwardDiff with preallocated Arrays

**URL:** <https://discourse.julialang.org/t/ipnewton-and-forwarddiff-with-preallocated-arrays/55013>\
**Category:** Optimization (Mathematical)\
**Created:** [February 10, 2021, 5:56pm UTC](https://discourse.julialang.org/t/ipnewton-and-forwarddiff-with-preallocated-arrays/55013 "2021-02-10T17:56:27Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sylvain](https://avatars.discourse-cdn.com/v4/letter/s/7ba0ec/32.png) [@Sylvain](https://discourse.julialang.org/u/Sylvain)\
**Post date:** [February 10, 2021, 5:56pm UTC](https://discourse.julialang.org/t/ipnewton-and-forwarddiff-with-preallocated-arrays/55013/1 "2021-02-10T17:56:27Z")

</div>

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

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

```julia
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
  cfg = ForwardDiff.GradientConfig(J, U)
  g! = (grad, x) -> ForwardDiff.gradient!(grad, J, x, cfg)
  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 :

```julia
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
 [3] value_gradient!!(::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}) at /Users/sylvainmoi/.julia/packages/NLSolversBase/QPnui/src/interface.jl:82
 [4] value_gradient! at /Users/sylvainmoi/.julia/packages/NLSolversBase/QPnui/src/interface.jl:69 [inlined]
 [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

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

```

to

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

```

it works but it’s useless, because each

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

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [February 10, 2021, 9:29pm UTC](https://discourse.julialang.org/t/ipnewton-and-forwarddiff-with-preallocated-arrays/55013/2 "2021-02-10T21:29:43Z")

</div>

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.

```julia
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)
  cfg = ForwardDiff.GradientConfig(obj, U)
  g! = (grad, x) -> ForwardDiff.gradient!(grad, obj, x, cfg)
  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()

```

---

<div class="post-metadata">

**Author:** ![Sylvain](https://avatars.discourse-cdn.com/v4/letter/s/7ba0ec/32.png) [@Sylvain](https://discourse.julialang.org/u/Sylvain)\
**Post date:** [February 11, 2021, 11:36am UTC](https://discourse.julialang.org/t/ipnewton-and-forwarddiff-with-preallocated-arrays/55013/3 "2021-02-11T11:36:39Z")

</div>

Thank you for your answer !

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

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

  println(J(U)) #works
  println(J(aux)) #works
  println(ForwardDiff.gradient(J, U)) #works
  println(ForwardDiff.gradient(J, aux)) #reports MethodError
  ...

```

I’m a bit confused because I can read in the [ForwardDiff doc](https://juliadiff.org/ForwardDiff.jl/stable/user/api/#Gradients-of-f(x::AbstractArray)::Real) 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`.

> [@cgeoga](#):
>
> 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.

I guess you are right. I don’t absolutely need to use preallocated arrays in my problem but I hoped they could speed up my code. I take note of your advice to write cleaner code.