Numerical experiments, Tips, Tricks and Gotchas

## Experiments with Wilkinson's polynomial

### Introduction

Wilkinson's polynomial is defined as [1, 2] $$W_{n}(x)=\prod_{m=1}^{n}(x-m)\label{eq:wp}$$ where $n=20$ . The expanded form of (\ref{eq:wp}) $$W_{n}(x)=\sum_{m=0}^{n}p_{m}x^{m}$$ is used for the analysis of numerical stability. Obviously the exact roots are $$\{r_{k}\}=\{1,2,3,\cdots,19,20\}\label{eq:rk}$$ 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).
 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)$.

 # 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 $$p_{19}=-210\rightarrow-210-2^{-23} = -\frac{1761607681}{8388608}\label{eq:modif}$$ 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).

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.

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.