rbf with gaussian kernel seems to produce more noise than original data

[eyurtsev]

This is probably user error, but I can't get rbf interpolation working in scipy. (I tried reading through the code, but it's a bit difficult to figure out what's going on...)

Basically, the more x data points I add the noisier the interpolation gets:

Below is a short script showing the problem.

I used sklearn's KneighborsRegressor to show what I expected to see as my output. I also re-implemented the same thing myself which gives identical output to that of sklearn (not included below.)

I am on scipy 0.15.1.

%pylab inline
import scipy
from scipy.interpolate import Rbf
from sklearn.neighbors import KNeighborsRegressor

# Generate some data
x = np.linspace(0, 10, 50)
z = x + 4.0 * randn(len(x))

# plot data
scatter(x, z, color='black', label='data')

# RBF regression using sklearn. Use all neighbors and gaussian kernel
def weights_func(dist_array):
    return exp(-dist_array**2)

est = KNeighborsRegressor(n_neighbors=50, weights=weights_func)
est.fit(x.reshape(-1, 1), z)
ym = est.predict(x.reshape(-1, 1))
plot(x, ym, color='purple', linestyle='-', linewidth=2, label='rbf using KneighborsRegressor')


# Rbf from scipy
rbf = Rbf(x, z, function='gaussian', epsilon=1)
plot(x, rbf(x), label='scipy-rbf')
legend(loc='center left', bbox_to_anchor=(1.05, 0.5), frameon=False)

[ev-br]

FWIW, the result seems to depend on the smooth parameter, eg rbf(x, y, function='gaussian', smooth=5.) definitely has less noise than the original data :-). Whether this is just bad set of defaults or there's a deeper problem with RBF as implemented in scipy, I've no idea.

[eyurtsev]

I am not familiar with the smoothness parameter, so the logic may be wrong, but based on the name of the parameter I expected to only see additional smoothness. I also expected the default setting to correspond to no additional smoothness (smoothness = 0?), such that the default interpolation would simply implemented weighted mean.

RBF interpolation using a Gaussian Kernel:

Given a data set composed of 2-tuples (xi, yi). The y value at location x is Y[x], given by:

Y[x] = sum_i(w(x, xi) * y[xi]) / (sum_i (w(x, xi))

And

w(x, xi) = exp(-(x-xi)^2 / r

[josef-pkt]

I haven't looked at this in a long time, and didn't check whether there are changes since then.

Try a negative smooth value for gaussian
see old thread here http://mail.scipy.org/pipermail/scipy-user/2010-March/024514.html
(no changes in scipy's RBF where made following that and similar discussions.)

[argriffing]

The 1/epsilon*r in this code is different than described in http://en.wikipedia.org/wiki/Radial_basis_function.

    def _h_multiquadric(self, r):
        return sqrt((1.0/self.epsilon*r)**2 + 1)

    def _h_inverse_multiquadric(self, r):
        return 1.0/sqrt((1.0/self.epsilon*r)**2 + 1)

    def _h_gaussian(self, r):
        return exp(-(1.0/self.epsilon*r)**2)

In Python, 1/e*r means (1/e)*r not 1/(e*r).

[josef-pkt]

@argriffing comparing your copied code with wikipedia, I think the only difference is that epsilon in wikipedia is 1/self.epsilon r is on the same side of the fraction as in wikipedia

e.g. exp(-(r * 1/eps)**2)

AFAICS, and whatever self.epsilon is

[argriffing]

@josef-pkt Yes thanks, I was just looking at that. I don't know what's going on with the scipy code. Either it is using the wrong formula or it's using a weird notation.

[argriffing]

The epsilon parameter is described as having the default

- defaults to approximate average distance between nodes
(which is a good start).

but I think this is implemented in scipy as the average distance over all pairs of nodes. For @eyurtsev's example problem this is like the average distance between points in the interval [0, 10] which is between 3 and 4. If I re-interpret the docstring as the average distance from a node to its nearest neighboring node, this is something more like 0.2. When I use epsilon=0.2 as input to the rbf function, I do not see the weirdness shown in the plot.

[argriffing]

I guess that because the example uses a weight function that has epsilon effectively 1 and uses all nodes as neighbors, this should be the same as the rbf function with epsilon=1 which has been specified explicitly in the function call. So my points about epsilon were pretty useless.

The smoothness docstring says

smooth : float, optional
        Values greater than zero increase the smoothness of the
        approximation.  0 is for interpolation (default), the function will
        always go through the nodal points in this case.

which suggests that without specifying a smoothness value, the scipy rbf will attempt to go exactly through the input points. I think this is not what @eyurtsev was expecting. Now that I read this, the scipy rbf sounds more like kriging...

[pv]

Afaik, determining the coefficients by solving the interpolation conditions is what "RBF interpolation" usually means. (Kriging is apparently this + using extra information on the covariances etc.)

Nevertheless, whether the regularization in Scipy's Rbf implementation makes sense would need to be checked, and ditto for the default epsilon values. Also, often 1st-order polynomials are thrown in the basis function set, as this guarantees some solvability conditions at least for the thin-plate splines.

[argriffing]

FWIW here's a simpler implementation for the 1d case with epsilon=1, no smoothing. There are stability issues as mentioned in http://duhctaep.blogspot.com/2010/05/kriging-with-numpy.html.

from __future__ import print_function, division

from functools import partial

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

def distances_1d(a, b):
    return np.absolute(a[np.newaxis, :] - b[:, np.newaxis])

def _solve_krige(x, w, u):
    r = distances_1d(x, u)
    B = np.exp(-r*r)
    v = np.dot(B, w)
    return v

def krige(x, z):
    r = distances_1d(x, x)
    A = np.exp(-r*r)
    w = scipy.linalg.solve(A, z)
    return partial(_solve_krige, x, w)

def main():
    x = np.linspace(0, 10, 15)
    z = x + 4.0 * np.random.randn(len(x))

    # plot data
    plt.scatter(x, z, color='black', label='data')

    u = np.linspace(0, 10, 50)

    # Rbf from scipy
    rbf = scipy.interpolate.Rbf(x, z, function='gaussian', epsilon=1)
    plt.plot(u, rbf(u), label='scipy-rbf')

    # custom function that is the same as scipy Rbf for 1d
    f = krige(x, z)
    plt.plot(u, f(u), color='purple', linestyle='-', linewidth=2, label='krige')

    plt.legend(loc='center left', bbox_to_anchor=(1.05, 0.5), frameon=False)
    plt.show()

main()

Edit: the point of this is to show @eyurtsev what Rbf is actually doing for his example.

[argriffing]

If these kernel functions are positive definite it may make more sense to use cholesky decomposition/solution. If you try that with this example (with 50 points) you can see that the decomposition raises an exception because the decomposition algorithm thinks the matrix is not positive definite, presumably because it is badly conditioned. This poor conditioning is probably also responsible for the weird plots.

[josef-pkt]

That's why we have a smoothing or penalization parameter.

[argriffing]

I think it's weird that some of the rbf kernel functions are not positive definite. Is this deliberate? I wanted to write a patch changing solve to cholesky for the kernel solver, but this is causing test failures because of the indefiniteness.

[josef-pkt]

don't use cholesky. I don't think we can assume positive definiteness (stronger than numerical noise)
extreme case, IIUC, if you have two identical observations, which I guess might be allowed.

maybe you should ask on the scikit-learn mailing list iwhether someone has more information about scipy's RBF. At least they have gaussian processes which is similar to the gaussian kernel RBF, IIUC or IIRC.

[pv]

The spectrum of the matrix is p.d. only for the inverse multiquadric and gaussian kernels (just check eigenvalues of rbf_obj.A), the other functions are not p.d. although they are commonly used in RBF.

To fix it up:

• Google scholar up some literature, as it's a mostly solved problem
• Use lstsq instead of solve (LAPACK probably has specialized functions for symmetric), or fix up the spectral shift approach to depend on the choice of the basis function. The shift probably needs to be negative for multiquadrics (only one positive eigenvalue, which is large, the rest are negative), and positive for p.d. basis functions. Properties of the spectrum of the popular kernels are known...
• The smoothing parameter should be given relative to the largest eigenvalue, similarly as condition number, so the largest eigenvalue should be estimated.
• The matrix A needs to be discarded in __init__ and not stored in self as it's not needed later (and may be big)
• For large data, the p.d. kernels are likely amenable for solution via iterative methods. The scaling of the required smoothing needs knowledge about the largest eigenvalue, but this is also easily found.
• There seem to be known ways to estimate the epsilon parameter. Mean distance is probably not a good choice.

Fixing this will break strict backward compat, but it seems this is warranted.