Vectorized Partial Derivative of Multivariate Function in Python

Question

There was a phenomenal answer posted by alko for computing a partial derivative of a multivariate function numerically in this thread.

I have a follow-up question now about enhancing this function to accept an array of input values. I have some code where I'm looping through a big long list of n-dimensional points, calculating the partial derivatives with respect to each variable, and this becomes quite computationally expensive.

It's easy enough to vectorize the function in question with np.vectorize, but it causes issues with the partial_derivative wrapper:

from scipy.misc import derivative
import numpy as np

def foo(x, y):
    return(x**2 + y**3)

def partial_derivative(func, var=0, point=[]):
    args = point[:]
    def wraps(x):
        args[var] = x
        return func(*args)
    return derivative(wraps, point[var], dx=1e-6)

vfoo = np.vectorize(foo)

>>>foo(3,1)
>>>10
>>>vfoo([3,3], [1,1])
>>>array([10,10])

>>>partial_derivative(foo,0,[3,1])
>>>6.0
>>>partial_derivative(vfoo,0,[[3,3], [1,1]])
>>>TypeError: can only concatenate list (not "float") to list

The last line should ideally return [6.0, 6.0]. In this case the two arrays supplied to the vectorized function vfoo are essentially zipped up pairwise, so ([3,3], [1,1]) gets transformed into two points, [3,1] and [3,1]. This seems to get mangled when it gets passed to the function wraps. The point that it ends up passing to the function derivative is [3,3]. In addition, there's obviously the TypeError thrown up.

Does anyone have any recommendations or suggestions? Has anyone ever had a need to do something similar?

Edit

Sometimes I think posting on SO is just what it takes to break a mental block. I think I've got it working for anyone who might be interested:

vfoo = np.vectorize(foo)

foo(3,1)
X = np.array([3,3])
Y = np.array([1,1])
vfoo(X, Y)
partial_derivative(foo,0,[3,1])
partial_derivative(vfoo,0,[X, Y])

And the last line now returns array([ 6., 6.])

Bear in mind that `np.vectorize` is not actually any faster than doing the same thing using a normal Python `for` loop - it's just syntactic sugar. — ali_m, Aug 20 '14 at 14:04

ntg · Answer 1 · 2015-08-05T09:21:57.283

I have a small problem with args[var] = x : this might forever change args[var] , and all values have been passed by reference however small your change is. So you might not get the exact answer you are looking for. Here is an example:

In[67]: a = np.arange(9).reshape(3,3)
In[68]: b = a[:]
In[69]: b[0,0]=42
In[70]: a
Out[70]: 
array([[42,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8]])

you need to fix it by e.g.:

def wraps(x):
    tmp = args[var]
    args[var] = x
    ret= func(*args)
    args[var] = tmp
    return ret

Also, you can use numdifftools. They seem to know what they are doing. This will do all the partial derivatives:

import numpy as np
import numdifftools as nd

def partial_function(f___,input,pos,value):
    tmp  = input[pos]
    input[pos] = value
    ret = f___(*input)
    input[pos] = tmp
    return ret

def partial_derivative(f,input):
    ret = np.empty(len(input))
    for i in range(len(input)):
        fg = lambda x:partial_function(f,input,i,x)
        ret[i] = nd.Derivative(fg)(input[i])
    return ret
if __name__ == "__main__":
    f     = lambda x,y: x*x*x+y*y
    input = np.array([1.0,1.0])
    print ('partial_derivative of f() at: '+str(input))
    print (partial_derivative(f,input))

Finally: if you want your function to take an array of the parameters, e.g.:

f     = lambda x: x[0]*x[0]*x[0]+x[1]*x[1]

then replace the respective line with (removed the '*')

ret = f___(input)

Vectorized Partial Derivative of Multivariate Function in Python

Edit

1 Answers1