Second order gradient in numpy

There's no universal right answer for numerical gradient calculation. Before you can calculate the gradient about sample data, you have to make some assumption about the underlying function that generated that data. You can technically use np.diff for gradient calculation. Using np.gradient is a reasonable approach. I don't see anything fundamentally wrong with what you are doing---it's one particular approximation of the 2nd derivative of a 1-D function.

I'll second @jrennie's first sentence - it can all depend. The numpy.gradient function requires that the data be evenly spaced (although allows for different distances in each direction if multi-dimensional). If your data does not adhere to this, than numpy.gradient isn't going to be much use. Experimental data may have (OK, will have) noise on it, in addition to not necessarily being all evenly spaced. In this case it might be better to use one of the scipy.interpolate spline functions (or objects). These can take unevenly spaced data, allow for smoothing, and can return derivatives up to k-1 where k is the order of the spline fit requested. The default value for k is 3, so a second derivative is just fine. Example:

spl = scipy.interpolate.splrep(x,y,k=3) # no smoothing, 3rd order spline
ddy = scipy.interpolate.splev(x,spl,der=2) # use those knots to get second derivative 

The object oriented splines like scipy.interpolate.UnivariateSpline have methods for the derivatives. Note that the derivative methods are implemented in Scipy 0.13 and are not present in 0.12.

Note that, as pointed out by @JosephCottham in comments in 2018, this answer (good for Numpy 1.08 at least), is no longer applicable since (at least) Numpy 1.14. Check your version number and the available options for the call.

As I keep stepping over this problem in one form or the other again and again, I decided to write a function gradient_n, which adds an differentiation oder functionality to np.gradient. Not all functionalities of np.gradient are supported, like differentiation of mutiple axis.

Like np.gradient, gradient_n returns the differentiated result in the same shape as the input. Also a pixel distance argument (d) is supported.

import numpy as np

def gradient_n(arr, n, d=1, axis=0):
    """Differentiate np.ndarray n times.

    Similar to np.diff, but additional support of pixel distance d
    and padding of the result to the same shape as arr.

    If n is even: np.diff is applied and the result is zero-padded
    If n is odd: 
        np.diff is applied n-1 times and zero-padded.
        Then gradient is applied. This ensures the right output shape.
    """
    n2 = int((n // 2) * 2)
    diff = arr

    if n2 > 0:
        a0 = max(0, axis)
        a1 = max(0, arr.ndim-axis-1)
        diff = np.diff(arr, n2, axis=axis) / d**n2
        diff = np.pad(diff, tuple([(0,0)]*a0 + [(1,1)] +[(0,0)]*a1),
                    'constant', constant_values=0)

    if n > n2:
        assert n-n2 == 1, 'n={:f}, n2={:f}'.format(n, n2)
        diff = np.gradient(diff, d, axis=axis)

    return diff

def test_gradient_n():
    import matplotlib.pyplot as plt

    x = np.linspace(-4, 4, 17)
    y = np.linspace(-2, 2, 9)
    X, Y = np.meshgrid(x, y)
    arr = np.abs(X)
    arr_x = np.gradient(arr, .5, axis=1)
    arr_x2 = gradient_n(arr, 1, .5, axis=1)
    arr_xx = np.diff(arr, 2, axis=1) / .5**2
    arr_xx = np.pad(arr_xx, ((0, 0), (1, 1)), 'constant', constant_values=0)
    arr_xx2 = gradient_n(arr, 2, .5, axis=1)

    assert np.sum(arr_x - arr_x2) == 0
    assert np.sum(arr_xx - arr_xx2) == 0

    fig, axs = plt.subplots(2, 2, figsize=(29, 21))
    axs = np.array(axs).flatten()

    ax = axs[0]
    ax.set_title('x-cut')
    ax.plot(x, arr[0, :], marker='o', label='arr')
    ax.plot(x, arr_x[0, :], marker='o', label='arr_x')
    ax.plot(x, arr_x2[0, :], marker='x', label='arr_x2', ls='--')
    ax.plot(x, arr_xx[0, :], marker='o', label='arr_xx')
    ax.plot(x, arr_xx2[0, :], marker='x', label='arr_xx2', ls='--')
    ax.legend()

    ax = axs[1]
    ax.set_title('arr')
    im = ax.imshow(arr, cmap='bwr')
    cbar = ax.figure.colorbar(im, ax=ax, pad=.05)

    ax = axs[2]
    ax.set_title('arr_x')
    im = ax.imshow(arr_x, cmap='bwr')
    cbar = ax.figure.colorbar(im, ax=ax, pad=.05)

    ax = axs[3]
    ax.set_title('arr_xx')
    im = ax.imshow(arr_xx, cmap='bwr')
    cbar = ax.figure.colorbar(im, ax=ax, pad=.05)

test_gradient_n()

The double gradient approach fails for discontinuities in the first derivative. As the gradient function takes one data point to the left and to the right into account, this continues/spreads when applying it multiple times.

On the other hand side, the second derivative can be calculated by the formula

d^2 f(x[i]) / dx^2 = (f(x[i-1]) - 2*f(x[i]) + f(x[i+1])) / h^2

compare here. This has the advantage to just take the two neighboring pixels into account.

In the picture the double np.gradient approach (left) and the above mentioned formula (right), as implemented by np.diff are compared. As f(x) has only one kink at zero, the second derivative (green) should only there have a peak. As the double gradient solution takes 2 neighboring points in each direction into account, this leads to finite second derivative values at +/- 1.

In some cases, however, you may want to prefer the double gradient solution, as this is more robust to noise.

I am not sure why there is np.gradient and np.diff, but a reason might be, that the second argument of np.gradient defines the pixel distance (for each dimension) and for images it can be applied for both dimensions simultaneously gy, gx = np.gradient(a).

Code

import numpy as np
import matplotlib.pyplot as plt

xs = np.arange(-5,6,1)
f = np.abs(xs)
f_x = np.gradient(f)
f_xx_bad = np.gradient(f_x)
f_xx_good = np.diff(f, 2)
test = f[:-2] - 2* f[1:-1] + f[2:]





# lets plot all this

fig, axs = plt.subplots(1, 2, figsize=(9, 3), sharey=True)

ax = axs[0]
ax.set_title('bad: double gradient')
ax.plot(xs, f, marker='o', label='f(x)')
ax.plot(xs, f_x, marker='o', label='d f(x) / dx')
ax.plot(xs, f_xx_bad, marker='o', label='d^2 f(x) / dx^2')
ax.legend()

ax = axs[1]
ax.set_title('good: diff with n=2')
ax.plot(xs, f, marker='o', label='f(x)')
ax.plot(xs, f_x, marker='o', label='d f(x) / dx')
ax.plot(xs[1:-1], f_xx_good, marker='o', label='d^2 f(x) / dx^2')
ax.plot(xs[1:-1], test, marker='o', label='test', markersize=1)
ax.legend()