Jacobian Matrix Calculation

Problem

Compute the Jacobian matrix of a vector-valued function f: R^n -> R^m. The Jacobian is the m x n matrix of all first-order partial derivatives, where entry (i, j) is df_i/dx_j.

Solution

import numpy as np

def jacobian(f, x: np.ndarray, epsilon: float = 1e-5) -> np.ndarray:
    x = np.array(x, dtype=float)
    f0 = np.array(f(x), dtype=float)
    n = len(x)
    m = len(f0)
    J = np.zeros((m, n))

    for j in range(n):
        x_plus = x.copy()
        x_minus = x.copy()
        x_plus[j] += epsilon
        x_minus[j] -= epsilon
        f_plus = np.array(f(x_plus), dtype=float)
        f_minus = np.array(f(x_minus), dtype=float)
        J[:, j] = (f_plus - f_minus) / (2 * epsilon)

    return J

Explanation

1. Evaluate the function at the given point x to determine the output dimension m.
2. For each input dimension j, apply the central difference formula: perturb x_j by +epsilon and -epsilon, evaluate f at both points, and compute (f(x+e_j*eps) - f(x-e_j*eps)) / (2*eps).
3. This gives the j-th column of the Jacobian, containing all partial derivatives with respect to x_j.
4. Central differences give O(epsilon^2) accuracy, better than forward differences.

Complexity

• Time: O(n * cost(f)) where n is the input dimension and cost(f) is the time to evaluate f
• Space: O(m * n) for the Jacobian matrix