Home

Recent
Archive

Numerical experiments, Tips, Tricks and Gotchas

Numerically speaking

Experiments with Wilkinson's polynomial

Introduction

Wilkinson's polynomial is defined as [1, 2] \begin{equation} W_{n}(x)=\prod_{m=1}^{n}(x-m)\label{eq:wp} \end{equation} where $n=20$ . The expanded form of (\ref{eq:wp}) \begin{equation} W_{n}(x)=\sum_{m=0}^{n}p_{m}x^{m} \end{equation} is used for the analysis of numerical stability. Obviously the exact roots are \begin{equation} \{r_{k}\}=\{1,2,3,\cdots,19,20\}\label{eq:rk} \end{equation} Let's explore this polynomial in various environments. We will use Mathematica, Matlab and Python (IPython notebook). In all systems it easy to generate Wilkinson's polynomial \[ W_{20}(x)=2432902008176640000-8752948036761600000\, x+\cdots-210\, x^{19}+x^{20} \] This polynomial is used in the numerical experiments.

Numerical root finding

All systems have built-in root finding routines. Below are the results (float = double precision).
Table 1. Numerical roots
 Mathematica
 absolute precision
 Matlab  Mathematica
 float
 Python  float
1  0.999999999999841  1.000000000000009  1.00000000
2  0.999999999997379  1.999999999999323  2.00000000
3  3.000000001049189  2.999999999984618  3.00000000
4  3.999999967630562  4.000000000014978  4.00000004
5  5.000000341909659  5.00000003983424  4.99999937
6  5.999999755869895  5.999998867455453  6.00000798
7  6.999973481009383  7.000015281443138  6.99993101
8  8.000284343435283  7.999873221974908  8.00039949
9  8.998394492431256  9.00071853610688  8.99845889
10 10.006059681252784  9.99708642835949 10.00386198
11 10.984041283443625 11.00892080340981 10.99484196
12 12.033449121964885 11.98006772341798 11.99847940
13 12.949055715498519 13.03569538754464 13.02424289
14 14.065272732694492 13.95372538519138 13.94558757
15 14.935355760260174 15.04620125219991 15.07956153
16 16.048274533412592 15.96503131308764 15.92526766
17 16.971132243131219 17.0180265952154 17.04459974
18 18.011221676102622 17.99333969776423 17.98069544
19 18.997159990805148 19.00144511371462 19.00457196
20 20.000324878101402 19.99985435328136 19.99949308

Neither the roots from Table 1 nor the exact roots (\ref{eq:rk}) give zero if we put them into $W_{20}(x)$.

Table 2. Polynomial values at exact roots
#  Mathematica
 absolute precision
 Matlab  Mathematica
 float
 Python  float
1 0 0.0000001024 160.0 1.02400000×103
2 0 0.0000008192 -6144.0 8.19200000×103
· ·  · · · · ·  · · · · ·  · · · · ·  · · · · ·
19 0 -1.6070033408 1.97998×1012 -1.60700334×1010
20 0 -2.7029504000 3.84829×1012 -2.70295040×1010 

The only exception is Mathematica with absolute precision.

Sensitivity experiment

As it was shown above, the polynomial values and its roots are extremely sensitive to the round of errors even if the calculations are performed with double precision.

Let's perform another experiment. If we replace \begin{equation} p_{19}=-210\rightarrow-210-2^{-23} = -\frac{1761607681}{8388608}\label{eq:modif} \end{equation} then the roots become: \begin{eqnarray*} \{\tilde{r}_{k}\} & = & \{1.00000,2.00000,3.00000,4.00000,5.00000,6.00001,6.99971,8.00714,\\ & & 8.91778,10.0953\pm0.642703i,11.7935\pm1.6521i,13.9923\pm2.51879i,\\ & & 16.7307\pm2.81262i,19.5024\pm1.94033i,20.8469\} \end{eqnarray*}

The above values are from Mathematica. The other environments give similar results.

The relative change $2^{-23}/210\approx5.7\cdot10^{-10}$ in $p_{19}$ caused large distortion of the roots. A half of them even became complex.

The derivative of a root with respect to $p_{19}$ can be considered as a measure of sensitivity. It can be shown [1] that their values at $r=\{1,2,\ldots,20\}$ are \[ d_{k}=\left.\frac{\partial r}{\partial p_{19}}\right|_{r=k}=\frac{k^{19}}{D_{k}} \] where \[ D_{k}=\prod_{l=1,l\neq k}^{20}(k-l) \] Their numerical values are \begin{eqnarray*} \{d_{k}\} & = & \{-8.22\cdot10^{-18},8.19\cdot10^{-11},-1.64\cdot10^{-6},2.19\cdot10^{-3},-6.08\cdot10^{-1},\\ & & 5.82\cdot10^{1},-2.54\cdot10^{3},5.97\cdot10^{4},-8.39\cdot10^{5},7.59\cdot10^{6},\\ & & -4.64\cdot10^{7},1.99\cdot10^{8},-6.06\cdot10^{8},1.33\cdot10^{9},-2.12\cdot10^{9},\\ & & 2.41\cdot10^{9},-1.90\cdot10^{9},9.96\cdot10^{8},-3.09\cdot10^{8},4.31\cdot10^{7}\} \end{eqnarray*} At the half of roots the absolute values of $d_{k}$ exceed $10^{7}$. This explains such big distortions of the roots.

Visual analysis

The plot of Wilkinson's polynomial is presented in Fig 1. (From Matlab).

Wilkinson's polynomial

Fig. 1. Wilkinson's polynomial.

The ordinate is $y=sign(p)\left[log_{10}(1+p)\right]$ because of large variation of the polynomial ($\pm10^{15}$ in magnitude). However, this transformation preserve a linear scale near $p = 0$.

The polynomial function is so steep in the vicinity of roots then even small round-off errors cause large deviation from zero. A relatively small variation in the coefficient $p_{19}$ also causes a large distortions of the roots. The plot of the modified polynomial is presented in Fig. 2.

Modified polynomial

Fig. 2. Modified polynomial.

Implication to polynomial fitting

The polynomial $W_{20}(x)$ can be considered as an exact fit of the function $f(x) \equiv 0$ at the points (\ref{eq:rk}). The plot in Fig. 1 shows that this is the extreme case of overfitting.

Numerical experiments

Below are scripts in Mathematica, Matlab and Python implementing the experiments. The IPython HTML notebook was converted to HTML and can be viewed at the link below.

For more details:

References

  1. G. Forsythe, M. Malcolm, and C. Moler, Computer Methods for Mathematical Computations, Prentice{Hall, Englewood Cliffs, NJ, 1977.
  2. Wikipedia, Wilkinson's polynomial.

© Nikolai Shokhirev, 2012-2024

email: nikolai(dot)shokhirev(at)gmail(dot)com

Count: