Question: I've been trying to port some code from using scipy.interpolate.Rbf to scipy.interpolate.RBFInterpolator. However I have the impression that the epsilon parameter has a different behavior in the latter -- in fact in my tests it seems like at least with the multiquadric kernel I can vary this parameter by multiple orders of magnitude with no appreciable change in output -- and from the existing scipy documentation it is not clear to me how this should be working (it describes the smoothing for RBFInterpolator very nicely, but only the Rbf documentation seems to explicitly show how epsilon enters the kernel functions as a scale parameter). To demonstrate this phenomenon, I have the following test code -- admittedly not a MWE since I made use of ROOT for visualizing the output. import sys import numpy as np import scipy import ROOT as rt def true_func(x,y,sigma,zscale): # first Gaussian, at center g1 = zscale * np.exp(-0.5 * (np.square(x) + np.square(y)) / np.square(sigma)) # second Gaussian, offset height_scale = 0.5 xp = x - 3. * sigma g2 = height_scale * zscale * np.exp(-0.5 * (np.square(xp) + np.square(y)) / np.square(sigma)) # add a couple sharper peaks positions = [(0,-2 * sigma), (0, -1 * sigma), (0, 2 * sigma)] spikes = 0 pow = 1.1 sig = sigma / 10 height_scale = 2. for pos in positions: xp = x - pos[0] yp = y - pos[1] spikes += height_scale * zscale * np.exp(-0.5 * (np.power(np.abs(xp),pow) + np.power(np.abs(yp),pow)) / np.power(sig,pow)) return g1 + g2 + spikes def test(new=False): N = 15 xscale = 100 xlin = np.linspace(-xscale*N,xscale*N,2 * N) ylin = np.linspace(-xscale*N,xscale*N,2 * N) x,y = np.meshgrid(xlin,ylin) # generate our z values rng = np.random.default_rng() zscale = 10 sigma = xscale * N / 4 z = true_func(x,y,sigma,zscale) z += 0.1 * zscale * rng.uniform(size=z.shape) xf = x.flatten() yf = y.flatten() zf = z.flatten() # Create two interpolators with different values of epsilon, # keep everything else the same between them. basis = 'multiquadric' rbf_dict = {} epsilon_vals = [0.1, 1000] if(new): for epsilon in epsilon_vals: rbf_dict[epsilon] = scipy.interpolate.RBFInterpolator( np.vstack((xf,yf)).T, zf, kernel=basis, epsilon=epsilon ) else: for epsilon in epsilon_vals: rbf_dict[epsilon] = scipy.interpolate.Rbf( xf,yf, zf, kernel=basis, epsilon=epsilon ) # now evaluate the two interpolators on the grid points = np.stack((x.ravel(), y.ravel()), axis=-1) if(new): evals = {key:val(points) for key,val in rbf_dict.items()} else: evals = {key:val(x,y) for key,val in rbf_dict.items()} diffs = {} for i,(key,val) in enumerate(evals.items()): if(i == 0): continue diffs[key] = (val - evals[epsilon_vals[0]]) print(np.max(diffs[key])) # now plot things dims = (1600,1200) c = rt.TCanvas('c1','c1',*dims) c.Divide(2,2) c.cd(1) true_graph = rt.TGraph2D(len(zf),xf,yf,zf) true_graph.SetName('true_graph') true_graph.Draw('SURF2Z') if(new): true_graph.SetTitle("scipy.interpolate.RBFInterpolator Test") else: true_graph.SetTitle("scipy.interpolate.Rbf Test") true_graph.GetXaxis().SetTitle("x") true_graph.GetYaxis().SetTitle("y") true_graph.GetZaxis().SetTitle("z") true_graph.SetNpx(80) true_graph.SetNpy(80) # now draw the two interpolations with the largest difference in epsilon interp1 = rt.TGraph2D(len(zf),xf,yf,evals[epsilon_vals[0]]) interp1.SetName('interp1') interp2 = rt.TGraph2D(len(zf),xf,yf,evals[epsilon_vals[-1]]) interp2.SetName('interp2') interp1.SetLineColor(rt.kRed) interp2.SetLineColor(rt.kGreen) # interp2.SetLineWidth(2) interp2.SetLineStyle(rt.kDotted) for g in (interp1, interp2): g.SetNpx(80) g.SetNpy(80) interp1.Draw('SAME SURF1') interp2.Draw('SAME SURF1') c.cd(2) diff_graph = rt.TGraph2D(len(zf),xf,yf,diffs[epsilon_vals[-1]]) diff_graph.SetName('diff_graph') diff_graph.SetTitle('Difference between interpolations, epsilon #in [{}, {}]'.format(epsilon_vals[0],epsilon_vals[-1])) diff_graph.Draw('SURF1') rt.gPad.SetLogz() c.cd(3) interp1.Draw('CONTZ') interp1.SetTitle('Interpolation with epsilon = {}'.format(epsilon_vals[0])) c.cd(4) interp2.Draw('CONTZ') interp2.SetTitle('Interpolation with epsilon = {}'.format(epsilon_vals[-1])) c.Draw() if(new): c.SaveAs('c_new.pdf') else: c.SaveAs('c_old.pdf') return def main(args): test(new=False) test(new=True) if(__name__=='__main__'): main(sys.argv)  This produces the following output plots: Output when using scipy.interpolate.Rbf Output when using scipy.interpolate.RBFInterpolator Maybe I am running into the consequences of some other changes to the RBF method, but it seems odd to me that the results are so different -- and that with using the scipy.interpolate.RBFInterpolator, I am seeing basically no difference between two interpolations using very different values of epsilon. I've taken a glance at what I think is the relevant scipy source code but so far have not figured out what is going on. Any help or advice would be much-appreciated. Thank you!

Answer:  However I have the impression that the epsilon parameter has a different behavior in the latter  It is different. Specifically, Rbf's epsilon divides distance, and RBFInterpolator multiplies it. Here is a source explaining why RBFInterpolator changes this from Rbf:  epsilon scales the RBF input as r*epsilon rather than r/epsilon, which is consistent with RBF literature (for example: https://www.sciencedirect.com/science/article/pii/S0898122107002210)  (Source.) If I change your RBFInterpolator call like this:  rbf_dict[epsilon] = scipy.interpolate.RBFInterpolator( np.vstack((xf,yf)).T, zf, kernel=basis, epsilon=1/epsilon )  I now get much more similar results between the two methods. Here is a comparison between large/small epsilon for old and new methods.   and that with using the scipy.interpolate.RBFInterpolator, I am seeing basically no difference between two interpolations using very different values of epsilon.  The issue here is connected with the first issue. You're not trying an epsilon value which is small enough. If the epsilon value is too large, then varying it will have no effect: the multiquadratic kernel is not scale-free, but it is approximately scale-free for large values of r. (Here is a plot showing that multiquadratic is approximately -r for large r. The function -r is one of the scale-free functions.) Here's the full test code. I changed epsilon, and also replaced the plotting code that used ROOT with pyplot. import sys import numpy as np import scipy from mpl_toolkits.mplot3d import Axes3D # Axes3D import has side effects, it enables using projection='3d' in add_subplot import matplotlib.pyplot as plt import random import numpy as np # import ROOT as rt def true_func(x,y,sigma,zscale): # first Gaussian, at center g1 = zscale * np.exp(-0.5 * (np.square(x) + np.square(y)) / np.square(sigma)) # second Gaussian, offset height_scale = 0.5 xp = x - 3. * sigma g2 = height_scale * zscale * np.exp(-0.5 * (np.square(xp) + np.square(y)) / np.square(sigma)) # add a couple sharper peaks positions = [(0,-2 * sigma), (0, -1 * sigma), (0, 2 * sigma)] spikes = 0 pow = 1.1 sig = sigma / 10 height_scale = 2. for pos in positions: xp = x - pos[0] yp = y - pos[1] spikes += height_scale * zscale * np.exp(-0.5 * (np.power(np.abs(xp),pow) + np.power(np.abs(yp),pow)) / np.power(sig,pow)) return g1 + g2 + spikes def test(new=False): N = 15 xscale = 100 xlin = np.linspace(-xscale*N,xscale*N,2 * N) ylin = np.linspace(-xscale*N,xscale*N,2 * N) x,y = np.meshgrid(xlin,ylin) # generate our z values rng = np.random.default_rng() zscale = 30 sigma = xscale * N / 4 z = true_func(x,y,sigma,zscale) z += 0.1 * zscale * rng.uniform(size=z.shape) xf = x.flatten() yf = y.flatten() zf = z.flatten() # Create two interpolators with different values of epsilon, # keep everything else the same between them. basis = 'multiquadric' rbf_dict = {} epsilon_vals = [0.1, 1000] if(new): for epsilon in epsilon_vals: rbf_dict[epsilon] = scipy.interpolate.RBFInterpolator( np.vstack((xf,yf)).T, zf, kernel=basis, epsilon=1/epsilon ) else: for epsilon in epsilon_vals: rbf_dict[epsilon] = scipy.interpolate.Rbf( xf,yf, zf, kernel=basis, epsilon=epsilon ) # now evaluate the two interpolators on the grid points = np.stack((x.ravel(), y.ravel()), axis=-1) if(new): evals = {key:val(points) for key,val in rbf_dict.items()} else: evals = {key:val(x,y) for key,val in rbf_dict.items()} diffs = {} for i,(key,val) in enumerate(evals.items()): if(i == 0): continue diffs[key] = (val - evals[epsilon_vals[0]]) print(np.max(diffs[key])) for epsilon in epsilon_vals: fig = plt.figure() plt.title(f"new {new}, epsilon {epsilon}") ax = fig.add_subplot(111, projection='3d') interp = rbf_dict[epsilon] if new: z = interp(points).reshape(x.shape) else: z = interp(x, y) ax.plot_surface(x, y, z) plt.show() def main(): test(new=False) test(new=True) if(__name__=='__main__'): main()