Gradiens of gradients using ReverseDiff


This is a duplicate of the issue I created here:

I am trying to get Hessian-vector product, without having to explicitly compute the Hessian matrix.
I have the following example:

using ReverseDiff

function f(x)
    m = length(x)
    return 100.0 * sum((x[i] - x[i - 1]^2)^2 for i=2:m) + (1.0 - x[1])^2

n = 2
x = [0.150369, 0.8463333]
u = [0.284309, 0.927797]

F = x -> ReverseDiff.gradient(f, x)
ϕᵤ(x) = F(x)'*u

Which gives the exact same thing as if I did: ReverseDiff.hessian(f, x)*u.

I know that calling ReverseDiff.gradient(f, x) isn’t the most efficient way to use ReverseDiff, so I tried to improve my code by doing: (everything before the definition of F remains the same)

tape = ReverseDiff.compile(ReverseDiff.GradientTape(f, rand(n)))
F = x -> ReverseDiff.gradient!(g, tape, x)
ϕᵤ(x) = F(x)'*u

But if I do that, the final answer is a Vector of zeros, which is not the expected result. What am I doing wrong?
Is there a way to do what I want efficiently and correctly?

Thanks for your help!