[SciPy-User] scipy.interpolate.rbf sensitive to input noise ?

[josef.pktd at gmail.com]

On Mon, Feb 22, 2010 at 7:35 AM, denis <denis-bz-gg at t-online.de> wrote:
> On Feb 19, 5:41 pm, josef.p... at gmail.com wrote:
>> On Fri, Feb 19, 2010 at 11:26 AM, denis <denis-bz... at t-online.de> wrote:
>
>> > 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.

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

>
> cheers
>  -- denis
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>