** Wilkinson's polynomial
:PROPERTIES:
:categories: polynomial
:date: 2014/02/21 09:54:47
:updated: 2014/02/21 09:55:18
:END:
[[http://en.wikipedia.org/wiki/Wilkinson%27s_polynomial][Wilkinson's polynomial]] is defined as
\( w(x) = \prod_{i=1}^{20} (x - i) = (x-1)(x-2) \ldots (x-20) \).
This innocent looking function has 20 roots, which are 1,2,3,...,19,20. Here is a plot of the function.
#+BEGIN_SRC python
import matplotlib.pyplot as plt
import numpy as np
@np.vectorize
def wilkinson(x):
p = np.prod(np.array([x - i for i in range(1, 21)]))
return p
x = np.linspace(0, 21, 1000)
plt.plot(x, wilkinson(x))
plt.ylim([-5e13, 5e13])
plt.savefig('./images/wilkinson-1.png')
#+END_SRC
#+RESULTS:
[[./images/wilkinson-1.png]]
Let us consider the expanded version of the polynomial. We will use sympy to expand the polynomial.
#+BEGIN_SRC python
from sympy import Symbol, Poly
from sympy.polys.polytools import poly_from_expr
x = Symbol('x')
W = 1
for i in range(1, 21):
W = W * (x-i)
print W.expand()
P,d = poly_from_expr(W.expand())
print P
#+END_SRC
#+RESULTS:
: x**20 - 210*x**19 + 20615*x**18 - 1256850*x**17 + 53327946*x**16 - 1672280820*x**15 + 40171771630*x**14 - 756111184500*x**13 + 11310276995381*x**12 - 135585182899530*x**11 + 1307535010540395*x**10 - 10142299865511450*x**9 + 63030812099294896*x**8 - 311333643161390640*x**7 + 1206647803780373360*x**6 - 3599979517947607200*x**5 + 8037811822645051776*x**4 - 12870931245150988800*x**3 + 13803759753640704000*x**2 - 8752948036761600000*x + 2432902008176640000
: Poly(x**20 - 210*x**19 + 20615*x**18 - 1256850*x**17 + 53327946*x**16 - 1672280820*x**15 + 40171771630*x**14 - 756111184500*x**13 + 11310276995381*x**12 - 135585182899530*x**11 + 1307535010540395*x**10 - 10142299865511450*x**9 + 63030812099294896*x**8 - 311333643161390640*x**7 + 1206647803780373360*x**6 - 3599979517947607200*x**5 + 8037811822645051776*x**4 - 12870931245150988800*x**3 + 13803759753640704000*x**2 - 8752948036761600000*x + 2432902008176640000, x, domain='ZZ')
The coefficients are orders of magnitude apart in size. This should make you nervous, because the roots of this equation are between 1-20, but there are numbers here that are O(19). This is likely to make any rounding errors in the number representations very significant, and may lead to issues with accuracy of the solution. Let us explore that.
We will get the roots using numpy.roots.
#+BEGIN_SRC python
import numpy as np
from sympy import Symbol
from sympy.polys.polytools import poly_from_expr
x = Symbol('x')
W = 1
for i in range(1, 21):
W = W * (x-i)
P,d = poly_from_expr(W.expand())
p = P.all_coeffs()
x = np.arange(1, 21)
print '\nThese are the known roots\n',x
# evaluate the polynomial at the known roots
print '\nThe polynomial evaluates to {0} at the known roots'.format(np.polyval(p, x))
# find the roots ourselves
roots = np.roots(p)
print '\nHere are the roots from numpy:\n', roots
# evaluate solution at roots
print '\nHere is the polynomial evaluated at the calculated roots:\n', np.polyval(p, roots)
#+END_SRC
#+RESULTS:
#+begin_example
These are the known roots
[ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20]
The polynomial evaluates to [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] at the known roots
Here are the roots from numpy:
[ 20.00032488 18.99715999 18.01122169 16.97113219 16.04827464
14.9353556 14.06527291 12.94905558 12.03344921 10.98404125
10.00605969 8.99839449 8.00028434 6.99997348 5.99999976
5.00000034 3.99999997 3. 2. 1. ]
Here is the polynomial evaluated at the calculated roots:
[40711209714176.0 15404160985600.0 8634610242048.00 3479686769152.00
1780604828160.00 694313602048.000 321293542400.000 150174387712.000
56110411264.0000 21911624192.0000 8370015744.00000 3104464384.00000
695443968.000000 125754368.000000 -947200.000000000 -9128960.00000000
-4393984.00000000 -712192.000000000 -31744.0000000000 17408.0000000000]
#+end_example
The roots are not exact. Even more to the point, the polynomial does not evaluate to zero at the calculated roots! Something is clearly wrong here. The polynomial function is fine, and it does evaluate to zero at the known roots which are integers. It is subtle, but up to that point, we are using only integers, which can be represented exactly. The roots function is evidently using some float math, and the floats are not the same as the integers.
If we simply change the roots to floats, and reevaluate our polynomial, we get dramatically different results.
#+BEGIN_SRC python
import numpy as np
from sympy import Symbol
from sympy.polys.polytools import poly_from_expr
x = Symbol('x')
W = 1
for i in range(1, 21):
W = W * (x - i)
P,d = poly_from_expr(W.expand())
p = P.all_coeffs()
x = np.arange(1, 21, dtype=np.float)
print '\nThese are the known roots\n',x
# evaluate the polynomial at the known roots
print '\nThe polynomial evaluates to {0} at the known roots'.format(np.polyval(p, x))
#+END_SRC
#+RESULTS:
#+begin_example
These are the known roots
[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.
16. 17. 18. 19. 20.]
The polynomial evaluates to [0 -8192.00000000000 -73728.0000000000 262144.000000000 716800.000000000
4055040.00000000 -200704.000000000 5767168.00000000 -13768704.0000000
152166400.000000 89210880.0000000 -146866176.000000 -91027456.0000000
-111190016.000000 405964800.000000 301989888.000000 -354531328.000000
-10256523264.0000 1316743168.00000 5308416000.00000] at the known roots
#+end_example
This also happens if we make the polynomial coefficients floats. That happens because in Python whenever one element is a float the results of math operations with that element are floats.
#+BEGIN_SRC python
import numpy as np
from sympy import Symbol
from sympy.polys.polytools import poly_from_expr
x = Symbol('x')
W = 1
for i in range(1, 21):
W = W * (x - i)
P,d = poly_from_expr(W.expand())
p = [float(x) for x in P.all_coeffs()]
x = np.arange(1, 21)
print '\nThese are the known roots\n',x
# evaluate the polynomial at the known roots
print '\nThe polynomial evaluates to {0} at the known roots'.format(np.polyval(p, x))
#+END_SRC
#+RESULTS:
:
: These are the known roots
: [ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20]
:
: The polynomial evaluates to [ 0.00000000e+00 -8.19200000e+03 -1.84320000e+04 -6.22592000e+05
: -2.04800000e+06 -1.08380160e+07 -2.31813120e+07 -5.89824000e+07
: -1.31383296e+08 -9.93280000e+07 -5.61532928e+08 -8.75003904e+08
: -1.38583245e+09 -1.97532877e+09 -3.80851200e+09 -6.02931200e+09
: -9.61910374e+09 -2.36191334e+10 -1.62105057e+10 -2.71933440e+10] at the known roots
Let us try to understand what is happening here. It turns out that the integer and float representations of the numbers are different! It is known that you cannot exactly represent numbers as floats.
#+BEGIN_SRC python
import numpy as np
from sympy import Symbol
from sympy.polys.polytools import poly_from_expr
x = Symbol('x')
W = 1
for i in range(1, 21):
W = W * (x - i)
P,d = poly_from_expr(W.expand())
p = P.all_coeffs()
print '{0:<30s}{1:<30s}{2}'.format('Integer','Float','\delta')
for pj in p:
print '{0:<30s}{1:<30f}{2}'.format(pj, float(pj), pj - float(pj))
#+END_SRC
#+RESULTS:
#+begin_example
Integer Float \delta
1 1.000000 0
-210 -210.000000 0
20615 20615.000000 0
-1256850 -1256850.000000 0
53327946 53327946.000000 0
-1672280820 -1672280820.000000 0
40171771630 40171771630.000000 0
-756111184500 -756111184500.000000 0
11310276995381 11310276995381.000000 0
-135585182899530 -135585182899530.000000 0
1307535010540395 1307535010540395.000000 0
-10142299865511450 -10142299865511450.000000 0
63030812099294896 63030812099294896.000000 0
-311333643161390640 -311333643161390656.000000 16.0000000000000
1206647803780373360 1206647803780373248.000000 112.000000000000
-3599979517947607200 -3599979517947607040.000000 -160.000000000000
8037811822645051776 8037811822645051392.000000 384.000000000000
-12870931245150988800 -12870931245150988288.000000 -512.000000000000
13803759753640704000 13803759753640704000.000000 0
-8752948036761600000 -8752948036761600000.000000 0
2432902008176640000 2432902008176640000.000000 0
#+end_example
Now you can see the issue. Many of these numbers are identical in integer and float form, but five of them are not. The integer /cannot/ be exactly represented as a float, and there is a difference in the representations. It is a small difference compared to the magnitude, but these kinds of differences get raised to high powers, and become larger. You may wonder why I used "0:<30s>" to print the integer? That is because =pj= in that loop is an object from sympy, which prints as a string.
This is a famous, and well known problem that is especially bad for this case. This illustrates that you cannot simply rely on what a computer tells you the answer is, without doing some critical thinking about the problem and the solution. Especially in problems where there are coefficients that vary by many orders of magnitude you should be cautious.
There are a few interesting webpages on this topic, which inspired me to work this out in python. These webpages go into more detail on this problem, and provide additional insight into the sensitivity of the solutions to the polynomial coefficients.
1. http://blogs.mathworks.com/cleve/2013/03/04/wilkinsons-polynomials/
2. http://www.numericalexpert.com/blog/wilkinson_polynomial/
3. http://en.wikipedia.org/wiki/Wilkinson%27s_polynomial