How does the "epsilon" parameter behave in scipy.interpolate.RBFInterpolator?
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!
However I have the impression that the
epsilonparameter 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:
epsilonscales the RBF input asr*epsilonrather thanr/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()