# Is there a way to compute the AD Hessian of a function using purely Zygote functions and not ForwardDiff?

To solve the issue raised in this post, I’ve been trying to write a Jacobian function using Zygote functions only so I can compute the AD Hessian of my function without using Dual type from ForwardDiff, which doesn’t work with FFTW. Here is the Jacobian function inspired by this.

``````function jacobian(f,x)
y,back  = Zygote.pullback(f,x)
k  = length(y)
n  = length(x)
J  = Matrix{eltype(y)}(undef,k,n)
e_mat = Matrix(I,k,k)
@inbounds for i = 1:k
J[i,:] = back(e_mat[:,i])[1]
end
(J,)
end

hessian(f, x) = jacobian(x -> gradient(f, x)[1], x)

f(x) = x[1]*x[2]^2
x0=rand(Float64,(2,1))
hessian(f,x0)
``````

However, this gives the following error:

``````Mutating arrays is not supported

error(::String)@error.jl:33
(::Zygote.var"#364#365")(::Nothing)@array.jl:58
(::Zygote.var"#150#151"{Zygote.var"#2245#back#366"{Zygote.var"#364#365"},Tuple{Tuple{Nothing,Nothing},Tuple{Nothing}}})(::Nothing)@lib.jl:191
#356@array.jl:38[inlined]
(::typeof(∂(λ)))(::Tuple{Array{Bool,2},Nothing})@interface2.jl:0
(::typeof(∂(λ)))(::Tuple{Nothing,Array{Bool,2},Nothing})@interface2.jl:0
literal_getindex@builtins.jl:15[inlined]
(::typeof(∂(λ)))(::Tuple{Nothing,Array{Bool,2},Nothing})@interface2.jl:0
f@Other: 1[inlined]
(::typeof(∂(λ)))(::Tuple{Nothing,Array{Bool,1}})@interface2.jl:0
#41@interface.jl:40[inlined]
(::typeof(∂(λ)))(::Tuple{Array{Bool,1}})@interface2.jl:0
I cannot think of a way to avoid mutating array, since we need it to construct the Jacobian matrix. I’m just trying to figure out if trying to write a Jacobian function using only Zygote is doomed? `f` is a very simple function (and another example here) yet it’s giving me errors every time. I can’t imagine how hard it would be to make it compatible to a function involving complex-numbers and FFT. Thank you so much for any guidance, and happy new year!