Polynôme de Wilkinson

Tracé du polynôme de Wilkinson

Tracé de ${\displaystyle \operatorname {sgn}(w(x))\ln(1+|w(x)|)}$

En analyse numérique, le polynôme de Wilkinson est un polynôme réel qui a été utilisé par James H. Wilkinson en 1963 pour illustrer un problème mal conditionné lors de la recherche de zéros d'un polynôme : la position des racines peut être très sensible aux perturbations des coefficients du polynôme.

Le polynôme est

${\displaystyle w(x)=\prod _{i=1}^{20}(x-i)=(x-1)(x-2)\cdots (x-20).}$

Parfois, le nom polynôme de Wilkinson est également utilisé pour désigner d'autres polynômes apparaissant dans l'étude de Wilkinson.

Préliminaires

Le polynôme de Wilkinson est apparu dans l'étude des algorithmes permettant de trouver les racines d'un polynôme

${\displaystyle p(x)=\sum _{i=0}^{n}c_{i}x^{i}.}$

C'est une question courante en analyse numérique de se demander si le problème de trouver les racines de p à partir des coefficients c_i est bien conditionné. Autrement dit, on espère qu’une petite variation dans les coefficients entraînera une petite variation dans les racines. Le polynôme de Wilkinson est un exemple de cas pathologique.

Le problème est mal conditionné lorsque le polynôme a une racine multiple. Par exemple, le polynôme x² a une racine double en 0. Cependant, le polynôme x ² − ε (avec une perturbation de taille ε ) a deux racines simples en ±√ ε, qui est beaucoup plus grand que ε lorsque ε est petit.

On en déduit qu'il faut s'attendre à ce qu’un mauvais conditionnement apparaisse aussi lorsque le polynôme a des zéros très proches. Cependant, le problème peut également être extrêmement mal conditionné pour les polynômes dont les zéros sont bien séparés. Wilkinson a utilisé le polynôme w(x) pour illustrer ce point (Wilkinson 1963).

En 1984, il décrit l’impact personnel de cette découverte :

« Pour ma part, je considère cela comme l’expérience la plus traumatisante de ma carrière en analyse numérique ^[1]. »

Le polynôme de Wilkinson est souvent utilisé pour illustrer le caractère indésirable du calcul naïf des valeurs propres d'une matrice en calculant d'abord les coefficients du polynôme caractéristique de la matrice, puis en trouvant ses racines, car l'utilisation des coefficients comme étape intermédiaire peut introduire un très mauvais conditionnement même si le problème d'origine était bien conditionné^[2].

Conditionnement du polynôme de Wilkinson

Le polynôme de Wilkinson

${\displaystyle w(x)=\prod _{i=1}^{20}(x-i)=(x-1)(x-2)\cdots (x-20)}$

a clairement 20 racines simples, situées en x = 1, 2, ..., 20. Ces racines sont très éloignées. Cependant, le polynôme est encore très mal conditionné.

En développant le polynôme, on trouve

${\displaystyle {\begin{aligned}w(x)={}&x^{20}-210x^{19}+20615x^{18}-1256850x^{17}[phone]x^{16}\\&{}-1672280820x^{15}[phone]x^{14}-756111184500x^{13}\\&{}[phone]x^{12}-135585182899530x^{11}\\&{}+1307535010540395x^{10}-10142299865511450x^{9}\\&{}+63030812099294896x^{8}-311333643161390640x^{7}\\&{}+1206647803780373360x^{6}-3599979517947607200x^{5}\\&{}+8037811822645051776x^{4}-12870931245150988800x^{3}\\&{}+13803759753640704000x^{2}-8752948036761600000x\\&{}+2432902008176640000.\end{aligned}}}$

Trace du logarithme du polynôme de Wilkinson et de sa version perturbée

Si le coefficient de x¹⁹ est réduit de 2⁻²³, soit de −210 à −210,0000001192, alors la valeur w(20) diminue de 0 à −2⁻²³ × 20¹⁹ =−6,25 ×10¹⁷, et la racine en x = 20 devient x ≈ 20,8. Les racines en x = 18 et x = 19 se rapprochent jusqu'à devenir une racine double en x ≈ 18,62 qui se transforme en deux racines conjuguées complexes en x ≈ 19,5 ± 1,9 i à mesure que la perturbation augmente. Les 20 racines deviennent (à 5 décimales)

${\displaystyle {\begin{array}{rrrrr}1,00000&2,00000&3,00000&4,00000&5,00000\\[8pt]6,00001&6,99970&8.00727&8.91725&20.84691\\[8pt]10,09527\pm {}&11,79363\pm {}&13,99236\pm {}&16,73074\pm {}&19,50244\pm {}\\[-3pt]0,64350\mathrm {i} &1,65233\mathrm {i} &2,51883\mathrm {i} &2,81262\mathrm {i} &1,94033\mathrm {i} \end{array}}}$

Certaines racines sont fortement déplacées, même si la modification du coefficient est minime et que les racines originales semblent très espacées. Wilkinson a montré par l'analyse de stabilité discutée dans la section suivante que ce comportement est lié au fait que certaines racines α (telles que α = 15) ont de nombreuses racines β qui sont « proches » au sens où |α − β| est plus petit que |α|.

Wilkinson a choisi la perturbation de 2⁻²³ parce que son ordinateur Pilot ACE avait des mantisses à virgule flottante de 30 bits, donc pour les nombres autour de 210, 2⁻²³ était une erreur dans la position du premier bit non représentée dans l'ordinateur. Les deux nombres réels, −210 et −210 − 2⁻²³, sont représentés par le même nombre à virgule flottante, ce qui signifie que 2⁻²³ est l'erreur inévitable en représentant un coefficient réel proche de −210 par un nombre à virgule flottante sur cet ordinateur. L'analyse des perturbations montre que la précision du coefficient de 30 bits est insuffisante pour séparer les racines du polynôme de Wilkinson.

Analyse de stabilité

On suppose que l'on perturbe un polynôme p(x) = Π(x − α_j) avec les racines α _j simples, en ajoutant un petit multiple t·c(x) d'un polynôme c(x), et on souhaite étudier le comportement des racines α_j . Au premier ordre, la variation des racines sera contrôlée par la dérivée

${\displaystyle {\mathrm {d} \alpha _{j} \over \mathrm {d} t}=-{c(\alpha _{j}) \over p^{\prime }(\alpha _{j})}.}$

Lorsque la dérivée est grande, les racines seront plus stables sous les variations de t, et inversement si cette dérivée est petite les racines seront instables. En particulier, si α _j est une racine multiple, alors le dénominateur s'annule. Dans ce cas, α _j n'est généralement pas dérivable par rapport à t (à moins que c ne s'y annule), et les racines seront extrêmement instables.

Pour de petites valeurs de t, la racine perturbée est donnée par le développement de Taylor t :

${\displaystyle \alpha _{j}+{\mathrm {d} \alpha _{j} \over \mathrm {d} t}t+{\mathrm {d} ^{2}\alpha _{j} \over \mathrm {d} t^{2}}{t^{2} \over 2!}+\cdots }$

et on s'attend à des problèmes quand |t| est plus grand que le rayon de convergence de cette série, qui est donné par la plus petite valeur de |t| tel que la racine α _j devienne multiple. Une estimation très grossière de ce rayon consiste à prendre la moitié de la distance entre α _j et la racine la plus proche et à la diviser par la dérivée ci-dessus.

Dans l'exemple du polynôme de Wilkinson de degré 20, les racines sont données par α_j = j pour j = 1, ..., 20, et c(x) est égal à x¹⁹. La dérivée est donc donnée par

${\displaystyle {\mathrm {d} \alpha _{j} \over \mathrm {d} t}=-{\alpha _{j}^{19} \over \prod _{k\neq j}(\alpha _{j}-\alpha _{k})}=-\prod _{k\neq j}{\alpha _{j} \over \alpha _{j}-\alpha _{k}}.\,\!}$

Cela montre que la racine α_j sera moins stable s'il existe de nombreuses racines α _k proches de α _j, dans le sens où la distance |α_j − α_k| entre eux est plus petite que |α _j| .

Exemple de bon comportement. Pour la racine α₁ = 1, la dérivée est égale à 1/19! ce qui est très petit ; cette racine est stable même pour de grands changements de t. En effet, toutes les autres racines β en sont loin, dans le sens où |α₁ − β | = 1, 2, 3, ..., 19 est supérieur à |α₁| = 1. Par exemple, même si t est aussi grand que –10000000000, la racine α₁ ne change que de 1 à environ 0,99999991779380 (ce qui est très proche de l'approximation du premier ordre 1 + t /19! ≈ 0,99999991779365). De même, les autres petites racines du polynôme de Wilkinson sont insensibles aux changements de t.

Exemple de mauvais comportement . Pour la racine α₂₀ = 20, la dérivée est égale à −20¹⁹/19! ce qui est très grand (environ 4 300 000), donc cette racine est très sensible aux petits changements de t. Les autres racines β sont proches de α₂₀, dans le sens où |β − α₂₀| = 1, 2, 3, ..., 19 est inférieur à |α₂₀| = 20. Pour t = −2⁻²³, l'approximation du premier ordre 20 − t × 20¹⁹/19! = 25,137... à la racine perturbée 20,84... est très mauvaise ; ceci est encore plus évident pour la racine α₁₉ où la racine perturbée a une grande partie imaginaire mais l'approximation du premier ordre (et d'ailleurs toutes les approximations d'ordre supérieur) sont réelles. La raison de cet écart est que |t| ≈ 0,000000119 est supérieur au rayon de convergence de la série entière mentionnée ci-dessus (qui est d'environ 0,0000000029, légèrement inférieur à la valeur 0,00000001 donnée par l'estimation brute), donc l'approximation linéaire ne s'applique pas. Pour une valeur telle que t = 0,000000001 qui est nettement inférieure à ce rayon de convergence, l'approximation du premier ordre 19,9569... est raisonnablement proche de la racine 19,9509...

À première vue, les racines α₁ = 1 et α₂₀ = 20 du polynôme de Wilkinson semblent similaires, car elles se trouvent aux extrémités opposées d'un ensemble symétrique de racines et ont le même ensemble de distances 1, 2, 3, ..., 19 par rapport aux autres racines. Cependant, l'analyse ci-dessus montre que cela est extrêmement trompeur : la racine α₂₀ = 20 est moins stable que α₁ = 1 (à de petites perturbations du coefficient de x¹⁹) d'un facteur de 20¹⁹ = 5,2488 × 10²⁴.

Deuxième exemple de Wilkinson

Le deuxième exemple considéré par Wilkinson est

${\displaystyle w_{2}(x)=\prod _{i=1}^{20}\left(x-{\frac {1}{2^{i}}}\right)=\left(x-{\frac {1}{2}}\right)\left(x-{\frac {1}{4}}\right)\cdots \left(x-{\frac {1}{2^{20}}}\right).}$

Les vingt zéros de ce polynôme sont dans une progression géométrique de raison 2, et donc le quotient

${\displaystyle \alpha _{j} \over \alpha _{j}-\alpha _{k}}$

ne peut pas être grand. En effet, les zéros de w₂ sont assez stables aux changements relatifs importants des coefficients.

Influence de la base de polynômes

Le développement

${\displaystyle p(x)=\sum _{i=0}^{n}c_{i}x^{i}}$

exprime le polynôme dans une base particulière, à savoir celle des monômes. Si le polynôme est exprimé dans une autre base, alors le problème de trouver ses racines peut cesser d'être mal conditionné. Par exemple, dans une forme de Lagrange, un petit changement dans un (ou plusieurs) coefficients ne modifie pas nécessairement trop les racines. En effet, les polynômes de base pour l'interpolation aux points 0, 1, 2, ..., 20 sont

${\displaystyle \ell _{k}(x)=\prod _{i\in \{0,\ldots ,20\}\setminus \{k\}}{\frac {x-i}{k-i}},\qquad {\text{pour}}\quad k=0,\ldots ,20.}$

Tout polynôme (de degré 20 ou moins) peut être exprimé sur cette base :

${\displaystyle p(x)=\sum _{i=0}^{20}d_{i}\ell _{i}(x).}$

Pour le polynôme de Wilkinson, on trouve

${\displaystyle w(x)=(20!)\ell _{0}(x)=\sum _{i=0}^{20}d_{i}\ell _{i}(x)\quad {\text{avec}}\quad d_{0}=(20!),\,d_{1}=d_{2}=\cdots =d_{20}=0.}$

Compte tenu de la définition du polynôme de base de Lagrange ℓ₀(x), un changement du coefficient d₀ ne produira aucun changement dans les racines de w. Cependant, une perturbation des autres coefficients (tous égaux à zéro) modifiera légèrement les racines. Le polynôme de Wilkinson est donc bien conditionné sur cette base.

Remarques

1. (en) James H. Wilkinson, Studies in Numerical Analysis, Mathematical Association of America, 1984, 3 p. (ISBN 978-0-88385-126-5), « The perfidious polynomial »
2. (en) Lloyd N. Trefethen et David Bau, Numerical Linear Algebra, SIAM, 1997

Références

Wilkinson a discuté de « son » polynôme dans

• (en) J. H. Wilkinson, « The evaluation of the zeros of ill-conditioned polynomials. Part I. », Numerische Mathematik, vol. 1,‎ 1959, p. 150–166 (DOI 10.1007/BF01386381).
• (en) J. H. Wilkinson, Rounding Errors in Algebraic Processes, Englewood Cliffs, New Jersey, Prentice Hall, 1963.

Il est mentionné dans les manuels standards d'analyse numérique, comme

• (en) F. S. Acton, Numerical methods that work, 201 p. (ISBN 978-0-88385-450-1).

Autres références :

• (en) Ronald G. Mosier, « Root neighborhoods of a polynomial », Mathematics of Computation, vol. 47, n^o 175,‎ juillet 1986, p. 265–273 (lire en ligne).
• (en) J. H. Wilkinson, « The perfidious polynomial », Studies in Numerical Analysis,‎ 1984 (lire en ligne)

Un calcul numérique de haute précision est présenté dans :

• (en) Ray Buvel, « Polynomials And Rational Functions », partie du RPN Calculator User Manual (en Python), vérifié le 27 septembre 2013.