N: 100 ngrid: 50 sigma: 0.1
# min max max |delta| Rbf
-1.0 1.0 0.5 gaussian
-1.0 1.2 0.3 linear
-1.4 2.1 2.5 thin-plate
-1.0 2.2 3.0 inverse multiquadric
-1.0 2.4 3.2 multiquadric
-2.1 12.8 7.4 quintic
-10.6 7.0 11.7 cubic
Should Rbf be moved off to noisymethods/... until experts can look it
over,
or have I done something stupid ?
(A back-of-the-envelope method for choosing epsilon aka r0 would be
nice,
and a rough noise amplification model would be nice too.)
Also, the doc should have a link to
http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data
and should say that the initial rbf() is N^2.
cheers
-- denis
<pre><code>
# from http://scipy.org/Cookbook/RadialBasisFunctions
import sys
import numpy as np
from rbf import Rbf # lstsq
# from scipy.interpolate import Rbf
Rbfuncs = ( 'gaussian', 'linear', 'thin-plate', 'inverse
multiquadric',
'multiquadric', 'quintic', 'cubic', )
N = 100
ngrid = 50
sigma = .1
exec "\n".join( sys.argv[1:] ) # N=
np.random.seed(1)
np.set_printoptions( 2, threshold=50, edgeitems=5,
suppress=True ) # .2f
x = np.sort( np.random.uniform( 0, 10, N+1 ))
y = np.sin(x) + np.random.normal( 0, sigma, x.shape )
xi = np.linspace( 0, 10, ngrid+1 )
print "N: %d ngrid: %d sigma: %.2g" % (N, ngrid, sigma)
print "# min max max |delta| Rbf"
for func in Rbfuncs:
rbf = Rbf( x, y, function=func )
fi = rbf(xi)
maxdiff = np.amax( np.abs( np.diff( fi )))
print "%5.1f %4.1f %4.1f %s" % (
fi.min(), fi.max(), maxdiff, func )
</pre></code>
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
Why would you do this? This does not magically turn Rbf interpolation
into a smoothing approximation. Use the "smooth" keyword to specify a
smoothing parameter.
> shows that all but linear and gaussian are noisy,
> giving interpolants way outside the input range:
>
> N: 100 ngrid: 50 sigma: 0.1
> # min max max |delta| Rbf
> -1.0 1.0 0.5 gaussian
> -1.0 1.2 0.3 linear
> -1.4 2.1 2.5 thin-plate
> -1.0 2.2 3.0 inverse multiquadric
> -1.0 2.4 3.2 multiquadric
> -2.1 12.8 7.4 quintic
> -10.6 7.0 11.7 cubic
>
> Should Rbf be moved off to noisymethods/... until experts can look it
> over,
> or have I done something stupid ?
I'm not sure what else you would expect from interpolating uniform
random noise. Many interpolation methods tend to give such
oscillations on such data.
--
Robert Kern
"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
-- Umberto Eco
> > Running rbf on 100 1d random.uniform input points
> > (after changing the linalg.solve in rbf.py to lstsq)
>
> Why would you do this? This does not magically turn Rbf interpolation
> into a smoothing approximation. Use the "smooth" keyword to specify a
> smoothing parameter.
Robert,
i'm interpolating y = np.sin(x) + np.random.normal( 0, .1 ), not
random noise.
The idea was to look at gaussian vs thin-plate;
true, that 1d snippet doesn't say much,
but my 2d plots were so noisy that I went down to 1d.
Use "smooth" ? rbf.py just does
self.A = self._function(r) - eye(self.N)*self.smooth
and you don't know A .
Bytheway (googling), http://www.farfieldtechnology.com/products/toolbox
...
have an O(N lnN) FastRBF TM in matlab, $
cheers
-- denis
Okay. You said otherwise. It would help if you simply attached the
code that you are using.
> The idea was to look at gaussian vs thin-plate;
> true, that 1d snippet doesn't say much,
> but my 2d plots were so noisy that I went down to 1d.
>
> Use "smooth" ? rbf.py just does
> self.A = self._function(r) - eye(self.N)*self.smooth
> and you don't know A .
I have no idea what you mean by the last part of your sentence. Did
you actually try using the smooth keyword? What were your results?
> Bytheway (googling), http://www.farfieldtechnology.com/products/toolbox
> ...
> have an O(N lnN) FastRBF TM in matlab, $
Okay. And?
--
Robert Kern
"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
-- Umberto Eco
rbf is a global interpolator, and as Robert said for some cases you
get very unstable results.
with a large number of points self._function(r) is not very well
behaved and you need the penalization to make it smoother. In a
variation of your example, I printed out the eigenvalues for self.A
and they don't look nice without penalization.
There are still some strange things that I don't understand with the
behavior rbf, but when I looked at it in the past, I got better
results by using only local information, i.e. fit rbf only to a
neighborhood of the points, although I think thinplate did pretty well
also with a larger number of points.
Josef
> > Use "smooth" ? rbf.py just does
> > self.A = self._function(r) - eye(self.N)*self.smooth
> > and you don't know A .
That's a line from scipy/interpolate/rbf.py: it solves
(A - smooth*I)x = b instead of
Ax = b
Looks to me like a hack for A singular, plus the caller doesn't know A
anyway.
I had looked at the eigenvalues too (100 points, 2d, like
test_rbf.py):
gauss : 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1.2e-10 ... 44
linear : -1.8e+02 ... -0.027 5.6e+02
thinplate : -4.7e+03 ... -5.4e+02 0.0032
So gauss is singular => don't do that then.
(Odd that linalg.solve didn't LinAlgError though.)
> rbf is a global interpolator, and as Robert said for some cases you
> get very unstable results.
>
> with a large number of points self._function(r) is not very well
> behaved and you need the penalization to make it smoother. In a
> variation of your example, I printed out the eigenvalues for self.A
> and they don't look nice without penalization.
> There are still some strange things that I don't understand with the
> behavior rbf, but when I looked at it in the past, I got better
> results by using only local information, i.e. fit rbf only to a
> neighborhood of the points, although I think thinplate did pretty well
> also with a larger number of points.
Yes -- but you don't know which points are local
without some k-nearest-neighbor algorithm ? in 2d, might as well
triangulate.
Gaussian weights nearer points automatically BUT rbf gauss looks to me
singular / unusable.
It's not a hack it's a requirement, ill-posed inverse problems need
penalization, this is just Ridge or Tychonov with a kernel matrix. A
is (nobs,nobs) and the number of features is always the same as the
number of observations that are used. (I was looking at "Kernel Ridge
Regression" and "Gaussian Process" before I realized that rbf is
essentially the same, at least for 'gauss')
I don't know anything about thinplate.
I still don't understand what you mean with "the caller doesn't know
A". A is the internally calculated kernel matrix (if I remember
correctly.)
>
> I had looked at the eigenvalues too (100 points, 2d, like
> test_rbf.py):
> gauss : 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 1.2e-10 ... 44
> linear : -1.8e+02 ... -0.027 5.6e+02
> thinplate : -4.7e+03 ... -5.4e+02 0.0032
>
> So gauss is singular => don't do that then.
> (Odd that linalg.solve didn't LinAlgError though.)
>
>> rbf is a global interpolator, and as Robert said for some cases you
>> get very unstable results.
>>
>> with a large number of points self._function(r) is not very well
>> behaved and you need the penalization to make it smoother. In a
>> variation of your example, I printed out the eigenvalues for self.A
>> and they don't look nice without penalization.
>> There are still some strange things that I don't understand with the
>> behavior rbf, but when I looked at it in the past, I got better
>> results by using only local information, i.e. fit rbf only to a
>> neighborhood of the points, although I think thinplate did pretty well
>> also with a larger number of points.
>
> Yes -- but you don't know which points are local
> without some k-nearest-neighbor algorithm ? in 2d, might as well
> triangulate.
I used scipy.spatial for this, or if there are not too many points use
the full dense distance matrix setting far away points to zero (kind
of). My example script should be on the mailing list.
rbf works for nd not just 1d or 2d, I think only the number of
observations is important, not the number of features.
> Gaussian weights nearer points automatically BUT rbf gauss looks to me
> singular / unusable.
unusable without large enough smooth>0 , e.g. I tried smooth=0.4 in bad cases.
Cheers,
Josef
> It's not a hack it's a requirement, ill-posed inverse problems need
OK, I must be wrong; but (sorry, I'm ignorant) how can (A - smooth)
penalize ?
For gauss the eigenvalues are >= 0, many 0, so we're shifting them
negative ??
Or is it a simple sign error, A + smooth ?
> penalization, this is just Ridge or Tychonov with a kernel matrix. A
> is (nobs,nobs) and the number of features is always the same as the
> number of observations that are used. (I was looking at "Kernel Ridge
> Regression" and "Gaussian Process" before I realized that rbf is
> essentially the same, at least for 'gauss')
> I don't know anything about thinplate.
>
> I still don't understand what you mean with "the caller doesn't know
> A". A is the internally calculated kernel matrix (if I remember
> correctly.)
Yes that's right; how can the caller of Rbf() give a reasonable value
of "smooth"
to solve (A - smoothI) inside Rbf, without knowing A ? A is wildly
different for gauss, linear ... too.
Or do you just shut your eyes and try 1e-6 ?
Thanks Josef,
>> I still don't understand what you mean with "the caller doesn't know
>> A". A is the internally calculated kernel matrix (if I remember
>> correctly.)
>
> Yes that's right; how can the caller of Rbf() give a reasonable value
> of "smooth"
> to solve (A - smoothI) inside Rbf, without knowing A ? A is wildly
> different for gauss, linear ... too.
Ah! *That's* what you meant.
> Or do you just shut your eyes and try 1e-6 ?
Basically, yes. You can try different values while watching the
residuals if you want a more rigorous approach. But it really does
work. Try using 'multiquadric' and smooth=0.1.
--
Robert Kern
"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
-- Umberto Eco
ouch, standard Ridge is A + smooth * identity_matrix to make it
positive definite.
I don't know why there is a minus. When I checked the eigenvalues, I
found it strange that there were some large *negative* eigenvalues of
A, but I didn't have time to figure this out.
Generically, A - smooth*eye would still make it invertible, although
not positive definite
I haven't looked at the sign convention in rbf, but if you figure out
what's going on, I'm very interested in an answer.
I just briefly ran an example with a negative smooth (-0.5 versus
0.5), rbf with gauss seems better, but multiquadric seems worse. If
smooth is small 1e-6, then there is not much difference.
Even with negative smooth, all except for gauss still have negative
eigenvalues. I have no idea, I only looked at the theory for gaussian
process and don't know how the other ones differ.
>
>> penalization, this is just Ridge or Tychonov with a kernel matrix. A
>> is (nobs,nobs) and the number of features is always the same as the
>> number of observations that are used. (I was looking at "Kernel Ridge
>> Regression" and "Gaussian Process" before I realized that rbf is
>> essentially the same, at least for 'gauss')
>> I don't know anything about thinplate.
>>
>> I still don't understand what you mean with "the caller doesn't know
>> A". A is the internally calculated kernel matrix (if I remember
>> correctly.)
>
> Yes that's right; how can the caller of Rbf() give a reasonable value
> of "smooth"
> to solve (A - smoothI) inside Rbf, without knowing A ? A is wildly
> different for gauss, linear ... too.
> Or do you just shut your eyes and try 1e-6 ?
That's the usual problem of bandwidth selection for non-parametric
estimation, visual inspection, cross-validation, plug-in, ... I don't
know what's recommended for rbf.
Cheers,
Josef
My intuition for gaussan might not be correct for the other rbfs,
gauss is decreasing in distance, all others are increasing in distance
>>> for func in Rbfuncs: print func, Rbf( x, y, function=func ,smooth=smooth)._function(np.arange(5))
gaussian [ 1. 0.9272 0.739 0.5063 0.2982]
linear [0 1 2 3 4]
thin-plate [ 0. 0. 2.7726 9.8875 22.1807]
multiquadric [ 1. 1.0371 1.1413 1.2964 1.4866]
quintic [ 0 1 32 243 1024]
cubic [ 0 1 8 27 64]
and the distance matrix itself has negative eigenvalues
Josef
Looking at Wikipedia Tikhonov (thanks Josef) reminded me
of lstsq_near on advice.mechanicalkern:
minimize |Ax-b| and w|x| together, i.e. Tikhonov with Gammamatrix = I.
But for RBF, why minimize |x| ? don't we really want to minimize
|Ax-b|2 + w2 xRx
where R is a roughness penalty ?
On a regular grid Rij ~ Laplacian smoother, not much like I.
Here I'm over my head;
any ideas for a q+d roughness matrix for scattered data ?
(just for fun -- for RBF we've reached the point of diminishing
returns).
Not necessarily. It seems to work well in at least some cases. Find a
reference that says otherwise if you want it changed.
> for gauss, start with A + 1e-6*I to move eigenvalues away from 0
> others have pos/neg eigenvalues, don't need smooth.
No. If you need a smoothed approximation to noisy data rather than
exact interpolation than you use the smooth keyword. Otherwise, you
don't.
--
Robert Kern
"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
-- Umberto Eco
chapter 2 page 16, for gaussian process. As I said I don't know about
the other methods
www.gaussianprocess.org/gpml/chapters/RW2.pdf
Josef
google links are very short, missed a part
I've found a couple of things on RBFs specifically that agree.
However, the original source on which our implementation is based does
subtract. He may have a source that uses a negative sign.
Playing around, it seems to me that some of the radial functions
create smoother approximations with large positive values (with the
current implementation) while others create smoother approximations
with large negative values. It's possible that it's just a convention
as to which sign you use.
>From the examples I checked, it looks like all other than gaussian
have large positive and negative eigenvalues, so adding or subtracting
doesn't change whether the matrix is definite.
For the gaussian case it is different, with smooth=0 the smallest
eigenvalues are zero and all others are positive. Only by *adding* a
penalization, the matrix becomes positive definite. (I've never seen a
case where eigenvalues are made negative to be able to invert a matrix
in the Ridge regression type of penalization.)
So it might be a convention if the kernel is increasing in distance
but not for the gaussian kernel that is decreasing.
If it's a convention than we could flip the sign and make it
correspond to the gaussian case, and the other cases would not be
strongly affected.
I will look at the reference. I haven't found much in terms of
references directly for general rbf.
Josef