Wilkinson's polynomial
Posted February 21, 2014 at 09:54 AM | categories: polynomial | tags:
Updated February 21, 2014 at 09:55 AM
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.
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')
Let us consider the expanded version of the polynomial. We will use sympy to expand the polynomial.
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
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.
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)
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]
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.
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))
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
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.
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))
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.
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))
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
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.
- http://blogs.mathworks.com/cleve/2013/03/04/wilkinsons-polynomials/
- http://www.numericalexpert.com/blog/wilkinson_polynomial/
- http://en.wikipedia.org/wiki/Wilkinson%27s_polynomial
Copyright (C) 2014 by John Kitchin. See the License for information about copying.
Org-mode version = 8.2.5h