[SciPy-user] Polyfit may be poorly conditioned

[massimo sandal]

Anne Archibald ha scritto:
> You should worry. That message indicates that the least-squares
> fitting is trying to solve a linear problem that is "ill-conditioned",
> that is, for which the solution depends very sensitively on the input.
> In combination with numerical errors, this can be a serious problem
> (though the fact that you're smoothing should reduce the issues).
> There are two things going on here.
> 
> First, high-order polynomials like to thrash wildly. This is basically
> because polynomials are very "stiff": they are very smooth, and every
> coefficient depends strongly on all the data. So if there's a kink at
> one place in your data (say), all your coefficients get bent out of
> shape trying to accomodate it. This can give polynomials extremely
> erratic behaviour even before the fit becomes ill-conditioned; it also
> means that the polynomials are going to become ill-conditioned quite
> soon.

I know and understand that polynomial behaviour.

Just to let you understand better what I'm trying to do: The data I'm 
fitting are basically a straight line, with noise, with a gentle 
oscillation superimposed. What I want is to remove the oscillation (a 
well known instrumental artefact). The fit of this line of practically 
pure noise acts as a "reference" for the true data set, that contains 
the oscillation AND the signal superimposed. So what I do is fitting the 
reference data, and subtracting the fit to the actual data.

(If you know about AFM force spectroscopy, the "reference" is the 
approaching curve, the data is the retracting curve, and the oscillation 
is optical interference)

Since the noise amplitude is very small compared to the large 
oscillation, spikes in the noise practically seem not to affect the fit. 
Of course, rare occurrences may happen where the fit goes wild, but I 
still haven't seen that -and a rare occurrence of that, say, 1%, is OK now.

> The second issue is that representing polynomials as a0 + a1*x +
> a2*x**2 + ... + an*x**n is not necessarily the best approach. If, say,
> your x value lies between -1 and +1, and your data lies between -1 and
> +1, then getting a polynomial to fit will require very large
> coefficients which lead to delicate cancellation. Numerical errors
> begin to be a problem surprisingly soon. A partial solution is to
> write your polynomials differently: instead of representing them as a
> linear combination of 1, x, x**2, ..., x**n, it sometimes helps to
> represent them as a linear combination of a set of orthogonal
> polynomials (Chebyshev polynomials, for example). This won't help if,
> like scipy, the implementation of orthogonal polynomials breaks down,
> but in fact the Chebyshev polynomials have a nice clean representation
> as cos(n*arccos(x)) (IIRC; look it up before using!). I have found
> that using this sort of representation can drastically increase the
> degree of polynomial I can fit before numerics become a problem. You
> can't use polyfit any more, unfortunately (though at some point
> someone might write a version of the polynomial class that used a
> different basis), but it's not too hard to write the least-squares
> fitting using scipy.linalg.lstsq (plus or minus a few vowels).

Thanks for the tip! I'll look into it if I meet problems with raw 
polynomials.

> Since you are computing polynomial coefficients (which is unstable)
> only to go back and compute the polynomial values at the same points
> you fit, it should be possible to set things up so the numerical
> instabilities are not really a problem - if some linear combination of
> polynomials gives nearly zero values, then you might as well just take
> it to be zero. If you want to look into this, you might try
> implementing your own polynomial fitting using the scipy.linalg.svd,
> which will give you control over how you deal with this situation.

Thanks for this tip, again.

>>  [1]this is needed to "flatten" the data set from systematic, irregular
>>  low-frequency interferences.
> 
> You might also think about whether you need to be subtracting
> polynomials. scipy's smoothing splines are extremely useful, and much
> more stable, numerically, than fitting polynomials. (This is primarily
> because they're less affected by distant data points.) You might also
> consider fitting a few sinusoids of prescribed frequencies (this may
> have easier-to-understand spectral properties).

Fitting sinusoids is something I've thought about, but having the easy 
polynomial fit at hand, I tried that and it seems it works, apart from 
the annoying warning.

As for splines, I've not tried that. I asked people before about splines 
and they told me that they are even more unstable than polynomials, 
however! so I didn't even try. Again, however, thank you to telling me that.

> I should mention that low-frequency "red" noise can pose serious
> problems for spectral analysis - if you're planning to do a Fourier
> transform later it's possible that strong low-frequency components may
> "leak" into your high frequency bins, contaminating them. Deeter and
> Boynton 1982 talk about how to analyze data with this kind of
> contamination, in the context of pulsar timing.

Spectral analysis should not be needed, however this is an interesting 
read, surely.

> In all I'd say your easiest solution is to use scipy's smoothing
> splines with a lot of smoothing. If you have to use polynomials, try
> using (scaled) Chebyshev polynomials as your basis. Only if these both
> fail is it worth going into some of my more exotic solutions.

Thanks a lot! I'll look more deeply if I effectively have troubles with 
polynomials, and I will surely think about your advices. I learned a lot!

m.

-- 
Massimo Sandal
University of Bologna
Department of Biochemistry "G.Moruzzi"

snail mail:
Via Irnerio 48, 40126 Bologna, Italy

email:
massimo.sandal at unibo.it

tel: +39-051-2094388
fax: +39-051-2094387
-------------- next part --------------
A non-text attachment was scrubbed...
Name: massimo.sandal.vcf
Type: text/x-vcard
Size: 274 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20080327/ef1b1873/attachment.vcf>