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).

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

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).

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.

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:

• Mathematica notebook can be viewed or downloaded
• Download MATLAB m-file
• The IPython HTML notebook can be viewed or downloaded.

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.