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
