A loss depends on many parameters at once. Its gradient is the vector of partial derivatives, one per parameter: how the loss changes when only that one parameter moves. The gradient points uphill, the direction of steepest increase, so gradient descent steps against it.

A numerical gradient moves each coordinate in turn by h in both directions and takes the central difference. It costs two evaluations per parameter, far too slow for training, and exactly right for checking an analytic gradient on a small case.

import numpy as np

def f(x):
    return x[0] ** 2 + 3 * x[0] * x[1]

x = np.array([1.0, 2.0])
grad = np.zeros_like(x)
for i in range(len(x)):
    e = np.zeros_like(x)
    e[i] = 1e-5
    grad[i] = (f(x + e) - f(x - e)) / 2e-5
print(grad)  # [8. 3.]

Compare the two gradients with a relative error, the size of their difference over the size of both. Below about 1e-7 they agree; above 1e-3 something is wrong.

Write numerical_gradient(f, x, h=1e-5) returning the central-difference gradient of f at the float array x, leaving x unchanged. Write relative_error(a, b) returning |a - b| / max(|a| + |b|, 1e-12) with L2 norms, as a float.

The explanation opens here once every check passes.