Hi,
You can find in this 126 gist a function that shows how to compute the hessian.
But you can compute the dot product with f-y much more efficiently using Rop from this 40 gist. In particular it will take Rop(loss, theta, f(theta) - y).