# [SciPy-user] calculating numerical jacobian

Pauli Virtanen pav@iki...
Sat Nov 15 07:17:25 CST 2008

```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

```