On Sat, Jan 28, 2012 at 11:14 PM, Charles R Harris < charlesr.har...@gmail.com> wrote:
> > > On Sat, Jan 28, 2012 at 11:15 AM, eat <e.antero.ta...@gmail.com> wrote: > >> Hi, >> >> Short demonstration of the issue: >> In []: sys.version >> Out[]: '2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit >> (Intel)]' >> In []: np.version.version >> Out[]: '1.6.0' >> >> In []: from numpy.polynomial import Polynomial as Poly >> In []: def p_tst(c): >> ..: p= Poly(c) >> ..: r= p.roots() >> ..: return sort(abs(p(r))) >> ..: >> >> Now I would expect a result more like: >> In []: p_tst(randn(123))[-3:] >> Out[]: array([ 3.41987203e-07, 2.82123675e-03, 2.82123675e-03]) >> >> be the case, but actually most result seems to be more like: >> In []: p_tst(randn(123))[-3:] >> Out[]: array([ 9.09325898e+13, 9.09325898e+13, 1.29387029e+72]) >> In []: p_tst(randn(123))[-3:] >> Out[]: array([ 8.60862087e-11, 8.60862087e-11, 6.58784520e+32]) >> In []: p_tst(randn(123))[-3:] >> Out[]: array([ 2.00545673e-09, 3.25537709e+32, 3.25537709e+32]) >> In []: p_tst(randn(123))[-3:] >> Out[]: array([ 3.22753481e-04, 1.87056454e+00, 1.87056454e+00]) >> In []: p_tst(randn(123))[-3:] >> Out[]: array([ 2.98556327e+08, 2.98556327e+08, 8.23588003e+12]) >> >> So, does this phenomena imply that >> - I'm testing with too high order polynomials (if so, does there exists a >> definite upper limit of polynomial order I'll not face this issue) >> or >> - it's just the 'nature' of computations with float values (if so, >> probably I should be able to tackle this regardless of the polynomial order) >> or >> - it's a nasty bug in class Polynomial >> >> > It's a defect. You will get all the roots and the number will equal the > degree. I haven't decided what the best way to deal with this is, but my > thoughts have trended towards specifying an interval with the default being > the domain. If you have other thoughts I'd be glad for the feedback. > > For the problem at hand, note first that you are specifying the > coefficients, not the roots as was the case with poly1d. Second, as a rule > of thumb, plain old polynomials will generally only be good for degree < 22 > due to being numerically ill conditioned. If you are really looking to use > high degrees, Chebyshev or Legendre will work better, although you will > probably need to explicitly specify the domain. If you want to specify the > polynomial using roots, do Poly.fromroots(...). Third, for the high degrees > you are probably screwed anyway for degree 123, since the accuracy of the > root finding will be limited, especially for roots that can cluster, and > any root that falls even a little bit outside the interval [-1,1] (the > default domain) is going to evaluate to a big number simply because the > polynomial is going to h*ll at a rate you wouldn't believe ;) > > For evenly spaced roots in [-1, 1] and using Chebyshev polynomials, things > look good for degree 50, get a bit loose at degree 75 but can be fixed up > with one iteration of Newton, and blow up at degree 100. I think that's > pretty good, actually, doing better would require a lot more work. There > are some zero finding algorithms out there that might do better if someone > wants to give it a shot. > > In [20]: p = Cheb.fromroots(linspace(-1, 1, 50)) > > In [21]: sort(abs(p(p.roots()))) > Out[21]: > array([ 6.20385459e-25, 1.65436123e-24, 2.06795153e-24, > 5.79026429e-24, 5.89366186e-24, 6.44916482e-24, > 6.44916482e-24, 6.77254127e-24, 6.97933642e-24, > 7.25459208e-24, 1.00295649e-23, 1.37391414e-23, > 1.37391414e-23, 1.63368171e-23, 2.39882378e-23, > 3.30872245e-23, 4.38405725e-23, 4.49502653e-23, > 4.49502653e-23, 5.58346913e-23, 8.35452419e-23, > 9.38407760e-23, 9.38407760e-23, 1.03703218e-22, > 1.03703218e-22, 1.23249911e-22, 1.75197880e-22, > 1.75197880e-22, 3.07711188e-22, 3.09821786e-22, > 3.09821786e-22, 4.56625520e-22, 4.56625520e-22, > 4.69638303e-22, 4.69638303e-22, 5.96448724e-22, > 5.96448724e-22, 1.24076485e-21, 1.24076485e-21, > 1.59972624e-21, 1.59972624e-21, 1.62930347e-21, > 1.62930347e-21, 1.73773328e-21, 1.73773328e-21, > 1.87935435e-21, 2.30287083e-21, 2.48815928e-21, > 2.85411753e-21, 2.85411753e-21]) > Thanks, for a very informative feedback. I'll study those orthogonal polynomials more detail. Regards, - eat > > > Chuck > > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion > >
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
