[SciPy-user] Jacobian matrix

[John Gleeson]

On Mar 17, 2005, at 9:59 AM, Nils Wagner wrote:

> Hi all,
>
> Is there a simpler way (built-in function) to compute the Jacobian of  
> a vector function ?
> How about higher order derivatives of vector-valued functions ?
>
>
> from scipy import *
>
> def f(x):
> #
>   tmp = zeros(5,Float)
>   tmp[0] = x[0] + x[1] + x[2]**2 + x[3]**2 + x[4]**2 - 2
>   tmp[1] = x[0] - x[1] + x[2]**2 + x[3]**2 + x[4]**2
>   tmp[2] = -x[2]**2 + x[3]**2 + x[4]**2
>   tmp[3] =  x[2]**2 - x[3]**2 + x[4]**2
>   tmp[4] =  x[2]**2 + x[3]**2 - x[4]**2
>   return tmp
>
> x0 = array(([1.02,1.02,0.02,0.02,0.02]))
> eps = 1.e-5
> J = zeros((5,5),Float)
> for i in arange(0,5):
>
>  ei = zeros(5,Float)
>  ei[i] = 1.0
>  J[:,i] = (f(x0+eps*ei)-f(x0))/eps
>
>

This approach doesn't uses the ScientificPython extension to scipy, so  
it doesn't satisfy your request for a built-in function, but it is  
fairly simple, and even provides higher order derivatives:

 >>> def f(x):
... #
...   tmp = zeros(5,PyObject)  # --- Notice change in typecode ---
...   tmp[0] = x[0] + x[1] + x[2]**2 + x[3]**2 + x[4]**2 - 2
...   tmp[1] = x[0] - x[1] + x[2]**2 + x[3]**2 + x[4]**2
...   tmp[2] = -x[2]**2 + x[3]**2 + x[4]**2
...   tmp[3] =  x[2]**2 - x[3]**2 + x[4]**2
...   tmp[4] =  x[2]**2 + x[3]**2 - x[4]**2
...   return tmp

 >>> x0 =  
array(([DerivVar(1.02,0,2),DerivVar(1.02,1,2),DerivVar(0.02,2,2),DerivVa 
r(0.02,3,2),DerivVar(0.02,4,2)]),PyObject)
 >>> D = f(x0)
 >>> D
array([(0.041199999999999903, [1.0, 1.0, 0.040000000000000001,  
0.040000000000000001, 0.040000000000000001], [[0.0, 0.0, 0.0, 0.0,  
0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0, 0.0], [0.0, 0.0,  
0.0, 2.0, 0.0], [0.0, 0.0, 0.0, 0.0, 2.0]]) ,
             (0.0012000000000000001, [1.0, -1.0, 0.040000000000000001,  
0.040000000000000001, 0.040000000000000001], [[0.0, 0.0, 0.0, 0.0,  
0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0, 0.0], [0.0, 0.0,  
0.0, 2.0, 0.0], [0.0, 0.0, 0.0, 0.0, 2.0]]) ,
             (0.00040000000000000002, [0.0, 0.0, -0.040000000000000001,  
0.040000000000000001, 0.040000000000000001], [[0.0, 0.0, 0.0, 0.0,  
0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, -2.0, 0.0, 0.0], [0.0, 0.0,  
0.0, 2.0, 0.0], [0.0, 0.0, 0.0, 0.0, 2.0]]) ,
             (0.00040000000000000002, [0.0, 0.0, 0.040000000000000001,  
-0.040000000000000001, 0.040000000000000001], [[0.0, 0.0, 0.0, 0.0,  
0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0, 0.0], [0.0, 0.0,  
0.0, -2.0, 0.0], [0.0, 0.0, 0.0, 0.0, 2.0]]) ,
             (0.00040000000000000002, [0.0, 0.0, 0.040000000000000001,  
0.040000000000000001, -0.040000000000000001], [[0.0, 0.0, 0.0, 0.0,  
0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0, 0.0], [0.0, 0.0,  
0.0, 2.0, 0.0], [0.0, 0.0, 0.0, 0.0, -2.0]]) ],'O')

# To get them into a more natural form:

 >>> f0 = array(([x[0] for x in D]))
 >>> f0
array([ 0.0412,  0.0012,  0.0004,  0.0004,  0.0004])
 >>> Df0 = array(([x[1] for x in D]))
 >>> Df0
array([[ 1.  ,  1.  ,  0.04,  0.04,  0.04],
        [ 1.  , -1.  ,  0.04,  0.04,  0.04],
        [ 0.  ,  0.  , -0.04,  0.04,  0.04],
        [ 0.  ,  0.  ,  0.04, -0.04,  0.04],
        [ 0.  ,  0.  ,  0.04,  0.04, -0.04]])
 >>> D2f0 = array(([x[2] for x in D]))
 >>> D2f0
array([[[ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  2.,  0.,  0.],
         [ 0.,  0.,  0.,  2.,  0.],
         [ 0.,  0.,  0.,  0.,  2.]],
        [[ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  2.,  0.,  0.],
         [ 0.,  0.,  0.,  2.,  0.],
         [ 0.,  0.,  0.,  0.,  2.]],
        [[ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0., -2.,  0.,  0.],
         [ 0.,  0.,  0.,  2.,  0.],
         [ 0.,  0.,  0.,  0.,  2.]],
        [[ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  2.,  0.,  0.],
         [ 0.,  0.,  0., -2.,  0.],
         [ 0.,  0.,  0.,  0.,  2.]],
        [[ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  0.,  0.,  0.],
         [ 0.,  0.,  2.,  0.,  0.],
         [ 0.,  0.,  0.,  2.,  0.],
         [ 0.,  0.,  0.,  0., -2.]]])