How to compute numerical derivatives with SciPy

Numerical derivatives let a Python workflow estimate slopes or curvature when a function is available only as callable code. SciPy exposes finite-difference tools in scipy.differentiate for checking scalar derivatives, Jacobian matrices, and Hessian matrices without deriving formulas by hand.

The derivative() function handles elementwise real scalar functions, while jacobian() and hessian() handle multivariable functions with vectorized call signatures. Each call returns a result object with the estimated derivative data plus convergence fields such as success, status, and error.

Finite differences sample the function near the requested point, so smooth functions and step sizes that are resolvable at the scale of the input matter. Use analytic derivatives or automatic differentiation when exact gradients are required, and check the SciPy status and error fields before feeding numerical derivatives into optimization or modeling code.

1. Create a Python script that checks a scalar derivative, a Jacobian, and a Hessian against known values.

   derivative_demo.py

   import numpy as np
   from scipy.differentiate import derivative, hessian, jacobian
    
    
   x0 = np.pi / 4
   scalar = derivative(np.sin, x0)
    
   print(f"derivative: {scalar.df:.8f}")
   print(f"expected: {np.cos(x0):.8f}")
   print(f"status: {int(scalar.status)}")
   print(f"error_estimate: {scalar.error:.2e}")
    
    
   def vector_function(x):
       return np.array([
           x[0]**2 + x[1],
           x[0] * x[1],
       ])
    
    
   point = np.array([2.0, 3.0])
   expected_jacobian = np.array([
       [4.0, 1.0],
       [3.0, 2.0],
   ])
   jac = jacobian(vector_function, point)
    
   print("jacobian:")
   for index, row in enumerate(jac.df):
       print(f"row {index}: {row}")
   print("jacobian success:", bool(jac.success.all()))
   print(f"jacobian max_abs_diff: {np.max(np.abs(jac.df - expected_jacobian)):.2e}")
    
    
   def quadratic(x):
       return x[0]**2 + 3 * x[0] * x[1] + 2 * x[1]**2
    
    
   expected_hessian = np.array([
       [2.0, 3.0],
       [3.0, 4.0],
   ])
   hes = hessian(quadratic, point)
    
   print("hessian:")
   for index, row in enumerate(hes.ddf):
       print(f"row {index}: {row}")
   print("hessian success:", bool(hes.success.all()))
   print(f"hessian max_abs_diff: {np.max(np.abs(hes.ddf - expected_hessian)):.2e}")

   jacobian() and hessian() expect callables that accept vectorized input. Indexing x[0] and x[1] lets SciPy evaluate nearby finite-difference points in one call.

2. Run the derivative script.

   $ python derivative_demo.py
   derivative: 0.70710678
   expected: 0.70710678
   status: 0
   error_estimate: 1.84e-12
   jacobian:
   row 0: [4. 1.]
   row 1: [3. 2.]
   jacobian success: True
   jacobian max_abs_diff: 1.87e-14
   hessian:
   row 0: [2. 3.]
   row 1: [3. 4.]
   hessian success: True
   hessian max_abs_diff: 5.74e-13

3. Confirm the scalar derivative result fields.

   status 0 means the estimate met the requested tolerance. The df value matches cos(pi/4) to eight printed decimals, and the reported error estimate is far below the displayed precision.

4. Read the Jacobian matrix from jac.df.

   For the vector function in the script, the two rows are output components and the two columns are input variables. The maximum absolute difference from the known matrix is 1.87e-14.

5. Read the Hessian matrix from hes.ddf.

   hessian() returns second-derivative values in ddf. The quadratic function has a constant Hessian, so the rows should match [2, 3] and [3, 4] within floating-point error.

6. Adjust finite-difference options when the function needs one-sided steps or a looser convergence target.

   Use step_direction near a domain boundary. If the true derivative may be exactly zero, pass a tolerances dictionary with an absolute tolerance such as {"atol": 1e-12}.

7. Remove the demo script when the check is complete.

   $ rm derivative_demo.py