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!

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!

• python
• scipy

2 Answers 2

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/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()

Basically, what you're seeing isn’t a glitch—it's just a change in how things are done. The old version of Rbf used epsilon to directly scale the distance, so even small tweaks could really shake things up. But the new RBFInterpolator has been reworked with extra stabilization and scaling tricks that essentially "normalize" epsilon’s role. This means that even if you change epsilon a lot, its impact is much more subtle, especially if your data is already nicely scaled.