rcond in polyfit

Charles R Harris charlesr.harris at gmail.com
Sat Oct 14 18:25:48 CDT 2006

On 10/14/06, A. M. Archibald <peridot.faceted at gmail.com> wrote:
> On 14/10/06, Charles R Harris <charlesr.harris at gmail.com> wrote:
> >
> >
> > On 10/13/06, A. M. Archibald <peridot.faceted at gmail.com> wrote:
> > > On 13/10/06, Tim Hochberg <tim.hochberg at ieee.org> wrote:
> > > > Charles R Harris wrote:
> >


Numerical Recipes (http://www.nrbook.com/a/bookcpdf/c15-4.pdf )
> recommend setting rcond to the number of data points times machine
> epsilon (which of course is different for single/double). We should
> definitely warn the user if any singular value is below s[0]*rcond (as
> that means that there is effectively a useless basis function, up to
> roundoff).

Well, that would work. On the other hand, it seems overly pessimistic from
my experiments here. What seems to be a better guide is rank reduction. For
instance, I can do a perfectly decent fit to Greg's data with single
precision, degree 10,  and ~1000 data points with rcond = ~5e-7 (effectively
single precision precision). Degree 11 blows up entirely even though it
loses rank. It is also true that rcond=1e-3 fails for degree 11 even though
rank is strongly reduced. What looks to be taking place is roundoff error in
evaluating the polynomial in single precision in the presence of higher
order terms that still belong to the reduced basis functions and rank
reduction is a good indicator of this. Bear in mind that I am now
normalizing x by dividing it by its largest element, which futzes with the
condition number. The condition number of the unscaled fit doesn't bear
thinking about.

Hmmm, I wonder if we have a dictionary of precisions indexed by dtype

> I don't get the impression that the warnings module is much tested; I
> had similar headaches.

I see some folks bitchin' at it when I google. There are remarkably few
hits, though. For the moment I am going with the smaller rcond numbers and
raising an error on rank reduction. I suspect something similar should be
done for pinv.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/numpy-discussion/attachments/20061014/5f4cd8d9/attachment.html 
-------------- next part --------------
Using Tomcat but need to do more? Need to support web services, security?
Get stuff done quickly with pre-integrated technology to make your job easier
Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo
-------------- next part --------------
Numpy-discussion mailing list
Numpy-discussion at lists.sourceforge.net

More information about the Numpy-discussion mailing list