# Is there an efficient way to compute the Hessian of a NN?

**URL:** <https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971>\
**Category:** Machine Learning\
**Tags:** flux\
**Created:** [July 30, 2019, 10:08am UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971 "2019-07-30T10:08:33Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![HenriDeh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrideh/32/8316_2.png) [@HenriDeh](https://discourse.julialang.org/u/HenriDeh)\
**Post date:** [July 30, 2019, 10:08am UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/1 "2019-07-30T10:08:33Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![alandion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alandion/32/5591_2.png) [@alandion](https://discourse.julialang.org/u/alandion)\
**Post date:** [July 30, 2019, 4:12pm UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/2 "2019-07-30T16:12:20Z")

</div>

One way to compute the Hessian that usually works pretty efficiently is to apply the reverse-on-forward trick, as outlined in

> **[1502.05767.pdf](https://arxiv.org/pdf/1502.05767.pdf)**
>
> 573.54 KB

If one replaces the floating point type with a dual number type, such as  
[http://www.juliadiff.org/ForwardDiff.jl/latest/dev/how\_it\_works.html#Dual-Number-Implementation-1](http://www.juliadiff.org/ForwardDiff.jl/latest/dev/how_it_works.html#Dual-Number-Implementation-1)

Then the usual backpropagation will, in addition to giving the gradient in the value of the dual number output, encode the Hessian-vector product in the dual number partial – the “vector” in Hessian-vector product being the list of partials with which your variables were initialized.

to compute the first diagonal element of the Hessian, you would initialize the partial of the first element of `x` to be 1.0 and all other partials to be 0.0

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 30, 2019, 4:24pm UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/3 "2019-07-30T16:24:48Z")

</div>

Reverse mode is going to give columns, and I don’t think you need that. Using double forward mode will be the fastest here. You’ll need to `mapchildren` to remove the tracker information (or use the `Flux#zygote` branch) and then just forward diff (or use a hyperdual)

If you do want to Forward-over-Reverse for Hess-vec products though, it is implemented in SparseDiffTools.jl

[https://github.com/JuliaDiffEq/SparseDiffTools.jl#jacobian-vector-and-hessian-vector-products](https://github.com/JuliaDiffEq/SparseDiffTools.jl#jacobian-vector-and-hessian-vector-products)

but note that our tests don’t show that using Zygote here is the fastest yet 🤷‍♂️

---

<div class="post-metadata">

**Author:** ![HenriDeh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrideh/32/8316_2.png) [@HenriDeh](https://discourse.julialang.org/u/HenriDeh)\
**Post date:** [August 2, 2019, 2:44pm UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/4 "2019-08-02T14:44:52Z")

</div>

Thank you both for your answers. And thank you for this nice paper, I understand AD much better now. I wanted to attempt both RoF and double forward methods. But I encounter a stack overflow issue when I use ForwardDiff methods on a neural network. Only the jacobian method takes an `AbstractArray` as input and accepts `AbstractArrays` as output. The other assume one of the two is Real. However:

```julia
f = Chain(Dense(10,5,elu),Dense(5,3))
x = rand(10)
ForwardDiff.jacobian(f, x) 

```

Throws a stack overflow. I assume this method does n forward passes, one for each input. But the underlying `ForwardDiff.gradient` cannot take multi-dimensional output functions as argument. Isn’t it weird for a forward accumulation mode to not be able to do that ? According to the paper from Baydin et al. it is quite straightforward to do. Same for the `ForwardDiff.derivative` function, it only handle functions with `Real ` inputs.

Am I missing something ?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 2, 2019, 3:04pm UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/5 "2019-08-02T15:04:53Z")

</div>

That will stack overflow 🤷‍♂️ an issue in flux (pun intended)

---

<div class="post-metadata">

**Author:** ![HenriDeh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrideh/32/8316_2.png) [@HenriDeh](https://discourse.julialang.org/u/HenriDeh)\
**Post date:** [August 5, 2019, 2:16pm UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/6 "2019-08-05T14:16:39Z")

</div>

I opened an issue then. By the way, indeed double forward mode is exactly what I need. The first pass will compute [dy\_1/dx\_1 … dy\_m/dx\_1] and the derivative with respect to x\_1 will give the dx\_1^2.

So I’ll go with this mode once I figure out how to use that with a `::Chain`. Reverse-on-forward would be best If you need the full hessian.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 5, 2019, 2:18pm UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/7 "2019-08-05T14:18:57Z")

</div>

It’s hardcoded in Flux to not work (and stackoverflow… it should throw a better error), but kind of for a good reason. We are working for a solution to this right now for other reasons, so stay tuned. We may start chatting about it in #autodiff

---

<div class="post-metadata">

**Author:** ![HenriDeh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrideh/32/8316_2.png) [@HenriDeh](https://discourse.julialang.org/u/HenriDeh)\
**Post date:** [August 5, 2019, 2:26pm UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/8 "2019-08-05T14:26:49Z")

</div>

Ha okay ! Great, thank you I’ll keep track of the progress then !

---

<div class="post-metadata">

**Author:** ![stochastic\_guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stochastic_guy/32/9877_2.png) [@stochastic\_guy](https://discourse.julialang.org/u/stochastic_guy)\
**Post date:** [August 27, 2019, 7:40pm UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/9 "2019-08-27T19:40:08Z")

</div>

Did you have any luck @HenriDeh ? I’m facing similar issues.

> [@\[Flux\] Issue while computing Hessian MethodError: objects of type Tracker.Grads are not callable](https://discourse.julialang.org/t/flux-issue-while-computing-hessian-methoderror-objects-of-type-tracker-grads-are-not-callable/27980):
>
> I’m new to Flux and I’m trying to compute the hessian. Here’s the code: grad = Flux.gradient(() -\> loss(x, y), Flux.params(model)) hess = Flux.gradient(grad, Flux.params(model)) But when running the line to compute hessian, it’s throwing the following error MethodError: objects of type Tracker.Grads are not callable

---

<div class="post-metadata">

**Author:** ![HenriDeh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrideh/32/8316_2.png) [@HenriDeh](https://discourse.julialang.org/u/HenriDeh)\
**Post date:** [August 28, 2019, 10:28am UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/10 "2019-08-28T10:28:23Z")

</div>

Hi, no unfortunately I haven’t heard of any changes regarding this. I know that Flux is getting a complete revamp with zygote and other stuff so this might not be their priority at the moment. It could also fix this in the process, I don’t know. I have put this aspect of my project on hold in the meantime.

---

<div class="post-metadata">

**Author:** ![elperkerson62](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elperkerson62/32/597_2.png) [@elperkerson62](https://discourse.julialang.org/u/elperkerson62)\
**Post date:** [November 11, 2019, 11:10pm UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/11 "2019-11-11T23:10:34Z")

</div>

> [@HenriDeh](#):
>
> `Flux.jacobian(::Chain, ::AbstractArray)` gives the full jacobian of a neural network. `Flux.hessian` however throws the error “output is not scalar”.

I had this same issue and resolved it using the `sum` function:

```julia
using Flux

m = Chain(Dense(2, 1))
x = zeros(2)
# Flux.hessian(m, x) throws "ERROR: Function output is not scalar"
# Flux.hessian(v->m(v)[1], x) throws "ERROR: Nested AD not defined for getindex"
Flux.hessian(v->sum(m(v)), x) # works

Tracked 2×2 Array{Float64,2}:
 0.0 0.0
 0.0 0.0

```

---

<div class="post-metadata">

**Author:** ![HenriDeh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrideh/32/8316_2.png) [@HenriDeh](https://discourse.julialang.org/u/HenriDeh)\
**Post date:** [November 12, 2019, 8:25am UTC](https://discourse.julialang.org/t/is-there-an-efficient-way-to-compute-the-hessian-of-a-nn/26971/12 "2019-11-12T08:25:43Z")

</div>

This only works if your network has one output. With multiple outputs nodes you cannot use a sum.