finite_difference_hessian is mutating arrays

[YichengDWu]

julia> function f(x,p)
       hess = FiniteDiff.finite_difference_hessian(y -> sum(y.^3), x)
       return hess * p
       end
f (generic function with 1 method)

julia> x=rand(3)
3-element Vector{Float64}:
 0.45298015977579165
 0.6696731824795704
 0.8825460798816239

julia> p=rand(3)
3-element Vector{Float64}:
 0.8610782341642329
 0.35801727566906194
 0.8214679367160949

julia> f(x,p)
3-element Vector{Float64}:
 2.3403081413900946
 1.438527409400547
 4.34989985490966

julia> Zygote.gradient(p->sum(f(x,p)), p)[1]
ERROR: Mutating arrays is not supported -- called setindex!(Matrix{Float64}, ...)
This error occurs when you ask Zygote to differentiate operations that change
the elements of arrays in place (e.g. setting values with x .= ...)

Possible fixes:
- avoid mutating operations (preferred)
- or read the documentation and solutions for this error
  https://fluxml.ai/Zygote.jl/dev/limitations.html#Array-mutation-1

Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:33
  [2] _throw_mutation_error(f::Function, args::Matrix{Float64})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\lib\array.jl:70
  [3] (::Zygote.var"#444#445"{Matrix{Float64}})(#unused#::Nothing)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\lib\array.jl:82
  [4] (::Zygote.var"#2496#back#446"{Zygote.var"#444#445"{Matrix{Float64}}})(Δ::Nothing)
    @ Zygote C:\Users\Luffy\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:67
  [5] Pullback
    @ C:\Users\Luffy\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\matmul.jl:514 [inlined]
  [6] Pullback
    @ C:\Users\Luffy\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\matmul.jl:510 [inlined]
  [7] (::typeof(∂(copytri!)))(Δ::Nothing)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
  [8] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:139 [inlined]
  [9] (::typeof(∂(#finite_difference_hessian!#43)))(Δ::Nothing)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [10] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:61 [inlined]
 [11] (::typeof(∂(finite_difference_hessian!##kw)))(Δ::Nothing)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [12] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:41 [inlined]
 [13] (::typeof(∂(#finite_difference_hessian#41)))(Δ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [14] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:39 [inlined]
 [15] (::typeof(∂(finite_difference_hessian##kw)))(Δ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [16] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:31 [inlined]
 [17] (::typeof(∂(#finite_difference_hessian#40)))(Δ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [18] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:30 [inlined]
 [19] (::typeof(∂(finite_difference_hessian)))(Δ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [20] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:30 [inlined]
 [21] (::typeof(∂(finite_difference_hessian)))(Δ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [22] Pullback
    @ .\REPL[19]:2 [inlined]
 [23] (::typeof(∂(f)))(Δ::FillArrays.Fill{Float64, 1, Tuple{Base.OneTo{Int64}}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [24] Pullback
    @ .\REPL[23]:1 [inlined]
 [25] (::typeof(∂(#15)))(Δ::Float64)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [26] (::Zygote.var"#60#61"{typeof(∂(#15))})(Δ::Float64)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface.jl:41
 [27] gradient(f::Function, args::Vector{Float64})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface.jl:76
 [28] top-level scope
    @ REPL[23]:1

[urbainvaes]

This issue seems to persist. Is it being worked on?

[ChrisRackauckas]

It is not. It's a bit tricky to do but if someone wants to put together a fully out of place version,

using FiniteDiff, Zygote

function f(x,p)
    hess = FiniteDiff.finite_difference_hessian(y -> sum(y.^3), x,
    FiniteDiff.HessianCache(x, Val(:hcentral), false))
    return hess * p
end
x=rand(3)
p=rand(3)
f(x,p)
Zygote.gradient(p->sum(f(x,p)), p)[1]

in theory that call should work. Right now it fails at

ArrayInterface.allowed_setindex!(H,(f(_xpp) - 2*fx + f(_xmm)) / epsilon^2,i,i)

ArrayInterface.allowed_setindex!(H,(f(_xpp) - f(_xpm) - f(_xmp) + f(_xmm))/(4*epsiloni*epsilonj),i,j)

i.e. the allowed_setindex calls which are used for a reason. If you're going to do it without this, then you need to build up all of the values in a map operation. It would probably be best to just write that as a separate dispatch for the inplace=false case.