[SciPy-user] calculating numerical jacobian

[Pauli Virtanen]

Fri, 14 Nov 2008 13:33:45 -0500, Tony S Yu wrote:
> Does scipy provide any functions to calculate the Jacobian of a
> function? I noticed that minpack (used by scipy.optimize) approximates a
> Jacobian numerically, but scipy doesn't seem to provide a wrapper for
> these minpack functions. Any ideas?

If you're happy with naive differentiation (as you may well be --
optimization and root finding typically isn't too picky), something
like this may work:

{{{
import numpy as np

def jacobian(f, u, eps=1e-6):
    """Evaluate partial derivatives of f(u) numerically.

    :note:
        This routine is currently naive and could be improved.

    :returns:
        (*f.shape, *u.shape) array ``df``, where df[i,j] ~= (d f_i / u_j)(u)
    """
    f0 = _N.asarray(f(u)) # asarray: because of matrices

    u_shape = u.shape
    nu = np.prod(u_shape)

    f_shape = f0.shape
    nf = np.prod(f_shape)

    df = np.empty([nf, nu])
    
    for k in range(nu):
        du = np.zeros(nu)
        du[k] = max(eps*abs(u.flat[k]), eps)
        f1 = np.asarray(f(u + _N.reshape(du, u_shape)))
        df[:,k] = np.reshape((f1 - f0) / eps, [nf])

    df.shape = f_shape + u_shape
    return df
}}}

Requires the function be vectorized.

-- 
Pauli Virtanen