`Flux.jacobian(::Chain, ::AbstractArray)`

gives the full jacobian of a neural network. `Flux.hessian`

however throws the error “output is not scalar”.

Let `m = Chain(Dense(5,3,elu), Dense(3,3))`

. What I can do is compute the gradient of each value in the jacobian matrix w.r. to each input using `Flux.gradient(x->Flux.jacobian(m, x)[i], input)`

where `i`

the index of the jacobian matrix we are interested in differentiating, and `input`

is the input vector fed to the neural network (i.e. the point where we want to differentiate). The jacobian has 3x5 dimensions, that means that to compute the hessian for each `i`

I must compute 15 times the same jacobian. That’s not efficient.

Moreover, I’m in fact only interested in what I would call the “second order jacobian” and not the full hessian. That is, the diagonals of the hessian (dy/dx^2). So computing the full hessian is wasteful anyway. In other words, what I would like is to compute the second order derivative of each output of a neural network with respect to each input.