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]) Chuck
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion