[SciPy-User] scipy.interpolate.rbf: how is "smooth" defined?
Mischa Schirmer
I'm doing a 2D fit using scipy.interpolate.rbf.
I have no problem with the fit itself, it works fine.
Data points are randomly scattered in a x-y plane,
and have a z-value each. The data is fairly well
behaved, in the sense that variations across-x-y
plane are slow. The data is somewhat noisy, and
thus I want to smooth it.
rbf has a 'smooth' parameter which is working well,
but it is exteremly poorly documented. Basically,
for smooth=0 no smoothing takes place, and for
smooth>0 some smoothing takes place.
Does anybody know how large the smoothing length
is? Is it internally modified depending on the
number of nodes entering the fit?
I'm asking because the number of nodes in my data
can vary from a few 10s to more than 10000.
Clearly, in the latter case I can do with a much
smaller smoothing scale, but I'd like to determine
it in some automatic fashion.
It appears that the smoothing length is
somehow normalised. At the moment I set 'smooth'
to 10% of the average distance between nodes
(in data units), resulting in a (numerically)
much larger effective smoothing scale.
Any insight is much appreciated!
Mischa Schirmer
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
josef...@gmail.com
<mis...@astro.uni-bonn.de> wrote:
> Hello,
>
> I'm doing a 2D fit using scipy.interpolate.rbf.
> I have no problem with the fit itself, it works fine.
>
> Data points are randomly scattered in a x-y plane,
> and have a z-value each. The data is fairly well
> behaved, in the sense that variations across-x-y
> plane are slow. The data is somewhat noisy, and
> thus I want to smooth it.
>
> rbf has a 'smooth' parameter which is working well,
> but it is exteremly poorly documented. Basically,
> for smooth=0 no smoothing takes place, and for
> smooth>0 some smoothing takes place.
smooth<0 also works, and, I think, is the correct sign for gaussian
>
> Does anybody know how large the smoothing length
> is? Is it internally modified depending on the
> number of nodes entering the fit?
No, smoothing factor is directly used.
My intuition mainly for gaussian case: the estimate depends on the
data and a prior that biases towards zero. The matrix in the normal
equations are a weighted average (except negative sign) between data
(kernel matrix) with weight 1, and an identity matrix with weight
smooth.
It's a bit similar to Ridge Regression, where relative weights can be
selected depending on the data (but often just given as parameter that
depends on the "prior beliefs").
I think, it's the usual bandwidth selection problem and scipy offers
no support for this (until someone adds it).
Also, I don't know the rbf literature, so I have no idea if there are
any easy or standard solutions for choosing the smoothing parameter,
but there should be bandwidth selection criteria analogous to other
smoothers.
>
> I'm asking because the number of nodes in my data
> can vary from a few 10s to more than 10000.
> Clearly, in the latter case I can do with a much
> smaller smoothing scale, but I'd like to determine
> it in some automatic fashion.
Is this true for all radial basis functions ?
If I remember correctly, I had to smooth more with many points, maybe
it depends on the density of points and the amount of noise.
Josef
denis
> rbf has a 'smooth' parameter which is working well,
> but it is exteremly poorly documented. Basically,
> for smooth=0 no smoothing takes place, and for
> smooth>0 some smoothing takes place
nope, it's just an offset for the matrix A that RBF solves:
interpolate/rbf.py line 191 (scipy 0.8.0rc3)
self.A = self._init_function(r) - eye(self.N)*self.smooth
self.nodes = linalg.solve(self.A, self.di)
Josef, we went back and forth on this back in February, and I thought
agreed
-- the -= should be a += (but no ticket)
-- for gaussian, use smooth = -1e-6 to nudge the eigenvalues of A >
0;
others have eigenvalues < 0 and > 0, so leave smooth=0 .
Mischa, are you using Gaussian ?
You could then try varying epsilon aka r0 in exp( - (r / r0)^2 )
which has a big affect on the fit.
Its "default to approximate average distance between [all] nodes" is
not so hot,
better av distance to nearby nodes, i.e. smaller Gaussian peaks.
cheers
-- denis
josef...@gmail.com
> On Aug 30, 1:10 pm, Mischa Schirmer <mis...@astro.uni-bonn.de> wrote:
>> Hello,
> ...
>> rbf has a 'smooth' parameter which is working well,
>> but it is exteremly poorly documented. Basically,
>> for smooth=0 no smoothing takes place, and for
>> smooth>0 some smoothing takes place
>
> nope, it's just an offset for the matrix A that RBF solves:
> interpolate/rbf.py line 191 (scipy 0.8.0rc3)
> self.A = self._init_function(r) - eye(self.N)*self.smooth
> self.nodes = linalg.solve(self.A, self.di)
>
> Josef, we went back and forth on this back in February, and I thought
> agreed
> -- the -= should be a += (but no ticket)
Yes, do you want to open a ticket?
> -- for gaussian, use smooth = -1e-6 to nudge the eigenvalues of A >
> 0;
1e-6 might be more like a minimum, in many cases it might need to be
much larger. But it would be useful to have better heuristics than
just trial and error.
Josef
Gael Varoquaux
> 1e-6 might be more like a minimum, in many cases it might need to be
> much larger. But it would be useful to have better heuristics than
> just trial and error.
I have not really been following the conversation (sorry), but are you
not talking about Tikhonov regularization of a covariance matrix? In
which case, I have had good experience using the Ledoit-Wolf approach to
setting the regularization factor.
I can provide numpy code for that, if needed.
Gaël
josef...@gmail.com
<gael.va...@normalesup.org> wrote:
> On Mon, Aug 30, 2010 at 11:57:47AM -0400, josef...@gmail.com wrote:
>> 1e-6 might be more like a minimum, in many cases it might need to be
>> much larger. But it would be useful to have better heuristics than
>> just trial and error.
>
> I have not really been following the conversation (sorry), but are you
> not talking about Tikhonov regularization of a covariance matrix? In
> which case, I have had good experience using the Ledoit-Wolf approach to
> setting the regularization factor.
It's the same idea, but instead of the covariance matrix, RBF uses the
kernel matrix. It will be worth a try, but I don't know whether the
choice of the regularization factor will apply in the same way here.
>
> I can provide numpy code for that, if needed.
Please do, I can also use it for other things. (I never read the small
print in Ledoit-Wolf.)
Josef
Gael Varoquaux
> Please do, I can also use it for other things. (I never read the small
> print in Ledoit-Wolf.)
There you go (BSD licensed). This will land in scikit learn at some point
because we are going to need it in different places, but I haven't found
the time to polish it.
Gael
def ledoit_wolf(x, return_factor=False):
""" Estimates the shrunk Ledoit-Wolf covariance matrix.
Parameters
----------
x: 2D ndarray, shape (n, p)
The data matrix, with p features and n samples.
return_factor: boolean, optional
If return_factor is True, the regularisation_factor is
returned.
Returns
-------
regularised_cov: 2D ndarray
Regularized covariance
regularisation_factor: float
Regularisation factor
Notes
-----
The regularised covariance is::
(1 - regularisation_factor)*cov
+ regularisation_factor*np.identity(n_features)
"""
n_samples, n_features = x.shape
if n_features == 1:
if return_factor:
return np.atleast_2d(x.std()), 0
return np.atleast_2d(x.std())
cov = np.dot(x.T, x)/n_samples
i = np.identity(n_features)
mu = np.trace(cov)/n_features
delta = ((cov - mu*i)**2).sum()/n_features
x2 = x**2
beta_ = 1./(n_features*n_samples) * np.sum(
np.dot(x2.T, x2)/n_samples - cov**2
)
beta = min(beta_, delta)
alpha = delta - beta
if not return_factor:
return beta/delta * mu * i + alpha/delta * cov
else:
return beta/delta * mu * i + alpha/delta * cov, beta/delta
josef...@gmail.com
Thanks,
If you use this with moment/normal equations for penalized
least-squares, do you have to multiply also `dot(x,y)/n_samples` by
`alpha/delta` ?
In the normalization of RBF, smooth would be `beta/delta * mu /
(alpha/delta)` or maybe additionally also divided by n_samples or
something like this.
Josef
Gael Varoquaux
> If you use this with moment/normal equations for penalized
> least-squares, do you have to multiply also `dot(x,y)/n_samples` by
> `alpha/delta` ?
I do not know, and I do not want to reply with something stupid. If you
find out the answer, I would be interest.
> In the normalization of RBF, smooth would be `beta/delta * mu /
> (alpha/delta)` or maybe additionally also divided by n_samples or
> something like this.
I believe it would be 'beta/delta * mu / (alpha/delta)', but I am not
sure whether or not the multiplication with n_samples should be applied,
as I haven't really looked at the exact formulation of the problem in the
context of RBFs (sorry, no time). All I know is that the trace of the
covariance matrix estimate should be kept constant as much as possible.
Gaël
josef...@gmail.com
<gael.va...@normalesup.org> wrote:
> On Mon, Aug 30, 2010 at 12:45:27PM -0400, josef...@gmail.com wrote:
>> If you use this with moment/normal equations for penalized
>> least-squares, do you have to multiply also `dot(x,y)/n_samples` by
>> `alpha/delta` ?
>
> I do not know, and I do not want to reply with something stupid. If you
> find out the answer, I would be interest.
Maybe you find out before me, with all the other things I don't find
time for, Ridge, informative Bayesian priors and Tychonov will have to
wait.
Josef
denis
>
> I'm doing a 2D fit using scipy.interpolate.rbf.
> I have no problem with the fit itself, it works fine.
>
> Data points are randomly scattered in a x-y plane,
> and have a z-value each. The data is fairly well
> behaved, in the sense that variations across-x-y
> plane are slow. The data is somewhat noisy, and
> thus I want to smooth it.
>
> rbf has a 'smooth' parameter which is working well,
> but it is exteremly poorly documented. Basically,
> for smooth=0 no smoothing takes place, and for
> smooth>0 some smoothing takes place.
Mischa,
are you seeing smoothing with smooth > 0 ? measured how ?
For varying epsilon aka r0 in exp( - (r / r0)^2 ),
r0 = side * sqrt( 2 / N )
looks to be of the right scale:
on a checkerboard of N points in a side^2 square,
a circle of that radius will enclose its 8 neighbours.
Well, the points are scattered, but.
Can you try epsilon = this r0, *2, *3
keeping smooth=0 ?
Of course any interpolation depends heavily on what you're fitting.
(Opinion:
You really want smaller Gaussian hills / smaller r0
where points are clustered, bigger where sparse.
RBF is not at all adaptive -- all hills are of size r0.
griddata (Delaunay triangulation + Natural Neighbour) is adaptive,
Invdisttree (kd tree + inverse-distance weighting) is adaptive.
)
cheers
-- denis
r0 = side * np.sqrt( 2 / N )
for smooth in ( 0, -1e-6 ):
print "smooth %.2g" % smooth
for c in range( 1, 5+1 ):
rbf = Rbf( x,y,z, function=func, epsilon=c*r0, smooth=smooth )
zi = rbf( xm.flatten(), ym.flatten() ) .reshape((ngrid,ngrid))