Gradients of gradients

Hi!
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
end

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

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

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
ReverseDiff.gradient(ϕᵤ,x)

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!

Start by comparing the ReverseDiff.gradient and gradient! calls, the compiled GradientTape, and the hessian entry point shown in the example. Trace why nested differentiation produces zeros and determine the supported efficient approach for a Hessian-vector product. Done means the example returns the expected nonzero result without explicitly forming the Hessian.