Interpreting Python code for partial derivatives of multi-variable function

Question

I'm trying to understand how the following code works to numerically approximate the partial derivative of a function of multiple variables. To quote the book:

When f is a function of many variables, it has multiple partial derivatives, each indicating how f changes when we make small changes in just one of the input variables. We calculate its ith partial derivative by treating it as a function of just its ith variable, holding the other variables fixed:

def partial_difference_quotient(f, v, i, h):
    """compute the ith partial difference quotient of f at v"""
    w = [v_j + (h if j == i else 0) for j, v_j in enumerate(v)] # add h to just the ith element of v

    return (f(w) - f(v)) / h

def estimate_gradient(f, v, h=0.00001):
    return [partial_difference_quotient(f, v, i, h) for i, _ in enumerate(v)]

Unfortunately, the books doesn't give any example of how to use the code; it just states that this will give you the partial derivative and then moves on.

The first thing I did is to try to rewrite the code so that it is more readable for someone with my background (I'm more used to seeing C++ code), but I'm not sure if I've done it right.

def partial_difference_quotient1(f, v, i, h):   # i is the variable that this function is being differentiated with respect to
    w = []
    for j, v_j in enumerate(v):
        if i == j:      # add h to just the ith element of v
            w.append(v_j + h)
        else:
            w.append(v_j)

    return (f(w) - f(v)) / h

def estimate_gradient1(f, v, h=0.00001):
    list1 = []
    for i, _ in enumerate(v):
        list1.append(partial_difference_quotient1(f, v, i, h))

    return list1

If the code is correct, then I'm still not sure how to use it. For instance, let's say I would like to find the partial derivative of f = x^2 + x*y^2 at y = 3, then I tried

def f1(x, y):
    return (x*x) + (x*y*y)

v = range(-10,10)
f = partial(f1, y=3)
f_x = estimate_gradient1(f, v, h=0.00001)

I can see that this wouldn't work because, for isntance, you can't multiply a list (i.e. v) with itself (x*x in the function f1). But it is clear from the code that v should be a list, so I've how to use the code.

Answer 1

Generally, you can work out the partial of a vector, you can do this:

For d/dx x^2+2x+1, the easiest way is to simply calculate the gradient for y in terms of an array of x:

from numpy import arange, gradient

x = arange(0, 50, 0.01)  # Range

y = ((x**2) + (2*x)) + 1  # Function

ddx = gradient(y)  # Gradient

which yields the following:

Another arbitrary example of dx/dt, d2/dy2, and phase plot would be as follows:

from scipy.integrate import odeint
from numpy import arange, sin, cos, gradient


def sdot(y, t, *args):
    val = sin(t*2)+(sin(y**(t/4)))
    return val

init = 0
x = arange(0, 30, 0.01)

args = (None)

dydt = odeint(sdot, y0=init, t=x, args=(args,))
d2dy2 = gradient(dydt[:, 0])

For additional information of gradient, see the documentations here, and for odeint see here. Also, you may find some additional information on vectorisation (aka. array programming) useful too.

For actual partial derivatives (i.e. through integration), you may also do this:

from scipy.misc import derivative

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)

Although this would be rather more challenging for someone who is not a hardcore (Python) programmer. But if you're keen to learn it, you can find the documentations here; and learn more about wrappers and decorators in this nice tutorial from the Code Ship in here.

Hope this gives you an idea of how you can handle derivatives in Python.

Answer 2

Here is how to use it:

def f(v):
    x, y = v
    return (x*x) + (x*y*y)

x = 3
y = 3
v = [x, y]
grad = estimate_gradient1(f, v, h=0.00001)
print(grad)

This outputs:

[15.000009999965867, 18.000030000564493]

which is approximately correct.