Approach: compute analytical gradients via autograd, compute numerical Jacobian via central finite differences for each scalar input element, handle vector outputs by flattening outputs and comparing Jacobian rows, and compare using configurable absolute and relative tolerances.python
import torch
def numerical_gradcheck(fn, inputs, eps=1e-6, atol=1e-5, rtol=1e-3, device=None):
"""
fn: function mapping tuple(inputs) -> tensor (possibly multi-dim). fn should be in eval mode.
inputs: tuple of torch.Tensor with requires_grad=True
Returns: (ok:bool, details:str)
"""
# Ensure inputs are tensors requiring grad
inputs = tuple(inp.clone().detach().requires_grad_(True) for inp in inputs)
# Evaluate output and analytic gradients
out = fn(*inputs)
out_flat = out.detach().reshape(-1)
# Compute analytic Jacobian: for each output scalar, backprop a unit vector
analytic_jacobians = []
for i in range(out_flat.numel()):
grads = torch.autograd.grad(out_flat[i], inputs, retain_graph=True, create_graph=False)
analytic_jacobians.append(torch.cat([g.reshape(-1) for g in grads]).detach())
analytic = torch.stack(analytic_jacobians) # shape [M, N] where M outputs, N total input elems
# Numerical Jacobian via central differences
num_elems = sum(inp.numel() for inp in inputs)
numerical = torch.zeros_like(analytic)
# Precompute flattened views
flat_inputs = [inp.reshape(-1) for inp in inputs]
for idx in range(num_elems):
# map idx to (tensor_i, offset)
cum = 0
for t_idx, flat in enumerate(flat_inputs):
if idx < cum + flat.numel():
offset = idx - cum
break
cum += flat.numel()
# perturb
orig = flat_inputs[t_idx][offset].item()
# +eps
flat_inputs[t_idx][offset] = orig + eps
out_p = fn(*[x.reshape(inp.shape) for x, inp in zip(flat_inputs, inputs)]).detach().reshape(-1)
# -eps
flat_inputs[t_idx][offset] = orig - eps
out_m = fn(*[x.reshape(inp.shape) for x, inp in zip(flat_inputs, inputs)]).detach().reshape(-1)
# restore
flat_inputs[t_idx][offset] = orig
# central difference
numerical[:, idx] = (out_p - out_m) / (2 * eps)
# Compare
abs_err = (analytic - numerical).abs()
rel_err = abs_err / (numerical.abs().clamp(min=1e-12) + analytic.abs().clamp(min=1e-12))
failed = (abs_err > atol) & (rel_err > rtol)
if failed.any():
first = torch.nonzero(failed)[0]
i_out, j_in = int(first[0]), int(first[1])
return False, f"Mismatch at output {i_out}, input idx {j_in}: analytic={analytic[i_out,j_in]}, numerical={numerical[i_out,j_in]}, abs_err={abs_err[i_out,j_in]}, rel_err={rel_err[i_out,j_in]}"
return True, "OK"
Key concepts:- Use central differences for O(eps^2) accuracy.- Flatten inputs and outputs to build full Jacobian for vector outputs.- Computationally expensive: O(M*N) forward evaluations (M outputs, N input elements).Numerical pitfalls:- eps too large → truncation error; too small → catastrophic cancellation due to floating-point; typical eps ~1e-6 for float32, 1e-8-1e-6 for float64.- Use central difference rather than forward/backward.- Compare using both absolute and relative tolerances because small true gradients need absolute checks while large values need relative checks.- Ensure deterministic fn (no randomness) and avoid in-place ops that break autograd.