scipy的RegularGridInterpolator可以通過一次調用返回值和漸變嗎？

我正在使用scipy.interpolate.RegularGridInterpolator和method='linear'。獲取內插值非常簡單（請參閱https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html上的示例）。獲取漸變和插值的好方法是什麼？scipy的RegularGridInterpolator可以通過一次調用返回值和漸變嗎？

一種可能性是多次調用插值器並使用有限差分「手動」計算梯度。這讓人覺得很浪費，因爲每次調用內插器都可能已經在計算引擎蓋下的漸變。那是對的嗎？如果是這樣，我該如何修改RegularGridInterpolator以返回插值函數值及其漸變？很明顯，我對我要插入的函數的「真實」梯度不感興趣 - 只是線性逼近的梯度，例如，在https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html的示例中的my_interpolating_function的梯度。

這裏是一個例子。我有一個函數f，我建立了一個線性內插器f_interp，我對f_interp（與f的梯度相反）的梯度感興趣。我可以使用有限差分來計算它，但有沒有更好的方法？我假設RegularGridInterpolator已經在計算引擎蓋下的漸變 - 並且快速完成。我如何修改它以返回漸變以及插值值？

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]

這裏的f和f_interp圖片。我感興趣的f_interp的梯度（實線）：

來源

2016-08-20 Adrian

號

這裏是scipy.interpolate.RegularGridInterpolator引擎蓋::

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)

下做使用scipy.ndimage.map_coordinates進行線性插值。座標包含像素座標中的位置。您應該能夠使用這些權重，並且可以使用每個維度的較低值和較高值來確定插值的陡峭程度。

但是，梯度還取決於角點的值。

你可以在這裏找到數學：https://en.wikipedia.org/wiki/Trilinear_interpolation

來源

2016-11-28 02:54:53 j13r

scipy的RegularGridInterpolator可以通過一次調用返回值和漸變嗎？

回答

相關問題