Can scipy's RegularGridInterpolator return both values and gradients with a single call?

Question

I'm using scipy.interpolate.RegularGridInterpolator with method='linear'. Getting interpolated values is straightforward (see example at https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html). What is a good way of getting gradients along with interpolated values?

One possibility is to call the interpolator multiple times and compute the gradient "manually" using finite differences. This feels wasteful given that each call to the interpolator is probably already computing the gradient under the hood. Is that correct? If so, how can I modify RegularGridInterpolator to return both the interpolated function value and its gradient?

To be clear, I'm not interested in the "true" gradient of the function which I'm interpolating -- just the gradient of the linear approximation, e.g. the gradient of my_interpolating_function in the example at https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html.

Here is an example. I have a function f, I build a linear interpolator f_interp, and I am interested in the gradient of f_interp (as opposed to the gradient of f). I can compute it using finite differences, but is there a better way? I assume RegularGridInterpolator is already computing the gradient under the hood -- and doing it fast. How can I modify it to return the gradient along with the interpolated values?

import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as interp

def f(x, y, z):
    return 0.01*x**2 + 0.05*x**3 + 5*x*y + 3*x*y*z + 0.1*x*(y**2)*z + 9*x*z**2 - 2*y

x_grid = np.linspace(0.0, 10.0, 20)
y_grid = np.linspace(-10.0, 10.0, 20)
z_grid = np.linspace(0.0, 20.0, 20)

mesh = np.meshgrid(x_grid, y_grid, z_grid, indexing="ij")

f_on_grid = f(mesh[0], mesh[1], mesh[2])
assert np.isclose(f_on_grid[0, 0, 0], f(x_grid[0], y_grid[0], z_grid[0]))  # Sanity check

grid = (x_grid, y_grid, z_grid)
f_interp = interp.RegularGridInterpolator(grid, f_on_grid, method="linear",
                                          bounds_error=False, fill_value=None)

dense_x = np.linspace(0.0, 20.0, 400)
plt.plot(dense_x, [f_interp([x, 1.0, 2.0])[0] for x in dense_x], label="y=1.0, z=2.0")
plt.plot(dense_x, [f_interp([x, 1.0, 4.0])[0] for x in dense_x], label="y=1.0, z=4.0")
plt.legend()
plt.show()

f_interp([0.05, 1.0, 2.0])  # Linearly interpolated value, distinct from f(0.05, 1.0, 2.0)

## Suppose I want to compute both f_interp and its gradient at point_of_interest
point_of_interest = np.array([0.23, 1.67, 5.88])
f_interp(point_of_interest)  # Function value -- how can I get f_interp to also return gradient?

## First gradient calculation using np.gradient and a mesh around point_of_interest +- delta
delta = 0.10
delta_mesh = np.meshgrid(*([-delta, 0.0, delta], ) * 3, indexing="ij")
delta_mesh_long = np.column_stack((delta_mesh[0].flatten(),
                                   delta_mesh[1].flatten(),
                                   delta_mesh[2].flatten()))
assert delta_mesh_long.shape[1] == 3
point_plus_delta_mesh = delta_mesh_long + point_of_interest.reshape((1, 3))
values_for_gradient = f_interp(point_plus_delta_mesh).reshape(delta_mesh[0].shape)
gradient = [x[1, 1, 1] for x in np.gradient(values_for_gradient, delta)]
gradient  # Roughly [353.1, 3.8, 25.2]

## Second gradient calculation using finite differences, should give similar result
gradient = np.zeros((3, ))
for idx in [0, 1, 2]:
    point_right = np.copy(point_of_interest)
    point_right[idx] += delta
    point_left = np.copy(point_of_interest)
    point_left[idx] -= delta
    gradient[idx] = (f_interp(point_right)[0] - f_interp(point_left)[0]) / (2*delta)

gradient  # Roughly [353.1, 3.8, 25.2]

Here's a picture of f and f_interp. I'm interested in the gradient of f_interp (solid lines):

score 1 · Accepted Answer · edited Jul 06 '18 at 13:48

No.

Here is what scipy.interpolate.RegularGridInterpolator does under the hood:

class CartesianGrid(object):
    """
    Linear Multivariate Cartesian Grid interpolation in arbitrary dimensions
    This is a regular grid with equal spacing.
    """
    def __init__(self, limits, values):
        self.values = values
        self.limits = limits

    def __call__(self, *coords):
        # transform coords into pixel values
        coords = numpy.asarray(coords)
        coords = [(c - lo) * (n - 1) / (hi - lo) for (lo, hi), c, n in zip(self.limits, coords, self.values.shape)]

        return scipy.ndimage.map_coordinates(self.values, coords, 
            cval=numpy.nan, order=1)

https://github.com/JohannesBuchner/regulargrid/blob/master/regulargrid/cartesiangrid.py

It uses scipy.ndimage.map_coordinates to do the linear interpolation. coords contains the location in pixel coordinates. You should be able to use these weights, and the lower and upper values at each dimension to figure out how steep the interpolation rises.

However, the gradient also depends on the values of the corner points.

You can find the math here: https://en.wikipedia.org/wiki/Trilinear_interpolation

Can scipy's RegularGridInterpolator return both values and gradients with a single call?

1 Answers1