Hi,
I am not sure you are generating your eps correctly: rand is a uniform in [0, 1], not a normal.
Also the value in the formula is for a variance close to 0.
We actually use this to check our gradients here pytorch/gradcheck.py at master · pytorch/pytorch · GitHub 27
The main things are:

Always use double precision otherwise the precision is just really bad
Use delta=1e-6 for double precision. No need to actually do sampling
You can do centered difference if you need more precision: f(z + e_t * eps/2) - f(z - e_t * eps/2). Even even more precise with an additional point at the center. You can check the wikipedia page for the exact formula you need to use in this case.
Don’t create the full I matrix. Just one vector and write a 1 at the proper index into it.
Disable autograd if you’re using nn.Module to avoid extra allocation with @torch.no_grad()