[SciPy-user] Polyfit may be poorly conditioned
Bruce Southey
bsouthey@gmail....
Thu Mar 27 08:38:03 CDT 2008
Hi,
To further Anne's reply, you really should at least rescale your
explanatory term say x dividing by some 'large' constant [instead of a1
+a2*x + a3x**2, first divide x by say 100 so you effectively fit
a1+a2*x/100 + a3*(x/100)**2]. This 'trick' also works well for nonlinear
models that would not or tend not to converge otherwise and really does
not change anything but the scale. The next step is using orthogonal
polynomials as these rather nice properties and usually work well as
Anne indicated.
Splines and loess transformations are extremely flexible and you really
should try these especially given the polynomial degree you are trying
to fit. Also these are semi-parametric so the model is not as important
as with polynomials and provide a 'local fit'. The latter means that you
can have a different fit depending on where you are in the data.
However, after saying that, much depends on the data and the information
content. I would also suggest somewhat strongly that you use a better
model comparison procedure than some mean square error like BIC
(Bayesian information criterion) or Akaike information criterion. That
way you can address model complexity - higher polynomial terms should
fit the data better than models with lower polynomial terms.
Regards
Bruce
massimo sandal wrote:
> 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.
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-user
mailing list