[SciPy-user] calculating numerical jacobian

[Rob Falck]

When I added scipy.optimize.fmin_slsqp for Scipy 0.70.0 I included a
function approx_jacobian that just does a basic forward finite difference
like approx_fprime.

http://www.scipy.org/scipy/scipy/attachment/ticket/570/slsqp.py

If you don't yet have 0.7.0, heres the code:

_epsilon = sqrt(finfo(float).eps)

def approx_jacobian(x,func,epsilon,*args):
    """Approximate the Jacobian matrix of callable function func

       * Parameters
         x       - The state vector at which the Jacobian matrix is
desired
         func    - A vector-valued function of the form f(x,*args)
         epsilon - The peturbation used to determine the partial derivatives
         *args   - Additional arguments passed to func

       * Returns
         An array of dimensions (lenf, lenx) where lenf is the length
         of the outputs of func, and lenx is the number of

       * Notes
         The approximation is done using forward differences

    """
    x0 = asfarray(x)
    f0 = func(*((x0,)+args))
    jac = zeros([len(x0),len(f0)])
    dx = zeros(len(x0))
    for i in range(len(x0)):
       dx[i] = epsilon
       jac[i] = (func(*((x0+dx,)+args)) - f0)/epsilon
       dx[i] = 0.0
    return jac.transpose()




On Sat, Nov 15, 2008 at 8:17 AM, Pauli Virtanen <pav at iki.fi> wrote:

> 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
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>



-- 
- Rob Falck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20081115/05ece5c0/attachment.html>