Perfidious Polynomials and Elusive Roots

1. Perfidious Polynomials and Elusive Roots Zhonggang Zeng Northeastern Illinois University Nov. 2, 2001 at Northern Illinois University

2. What are the “best” functions? 1. Linear functions 2. Polynomials 3. Exponential functions 4. Trigonometric functions ... ... Weierstrass Theorem: Any continuous function can be approximated to any accuracy by a polynomial of sufficiently high degree. Why by a polynomial ?

3. Polynomial root-finding xn + a1 xn-1 + a2 xn-2 + ... + an-1x + a0 = 0 has played a key role in the history of mathematics - Rational/irrational numbers, complex numbers - Fundamental Theorem of Algebra - Abel and Galois Theorem ... ...

4. Can you solve (x-1.0)100 = 0 Can you solve x100-100 x99 +4950 x98 - 161700 x97+3921225x96 - ... - 100 x +1 = 0

5. Eigenvalues of 10 1 1 1 1 ... ... 1 1 0 1 1 A = X X-1

6. The Wilkinson polynomial p(x) = (x-1)(x-2)...(x-20) = x20 - 210 x19 + 20615 x18 + ... Wilkinson wrote in 1984: Speaking for myself I regard it as the most traumatic experience in my career as a numerical analyst.

7. Classical textook methods for multiple roots Newton’s iteration xj+1 = xj - f(xj)/f’(xj), j=0,1,2,... converges locally to a multiple root of f(x) with a linear rate. The modified Newton’s iteration xj+1 = xj - mf(xj)/f’(xj), j=0,1,2,... converges locally to a m-fold root of f(x) with a quadratic rate. Newton’s iteration applied to g(x) = f(x)/f’(x) converges locally and quadratically to a root of f(x) regardliss of its multiplicity. None of them work!

8. Example: f(x) = (x-2)7(x-3)(x-4) in expanded form. Modified Newton’s iteration with m = 7 intended for root x = 2: x1 = 1.9981 x2 = 1.7481 x3 = 1.9892 x4 = 0.4726 x5 = 1.8029 x6 = 1.9931 x7 = 4.2681 x8 = 3.3476 ... ...

9. Myths on polynomial roots: - multiple roots are ill-conditioned, or even intractable - extension of machine precision is necessary to calculate multiple roots - there is an “attainable precision” for multiple roots: machine precision attainable precision = ----------------------------- multiplicity Example: for 100-fold root, to get 5 digits right 500 digits in machine precision 5 digits precision = ----------------------------------------- 100 in multiplicity

10. How do we justify the answer?

11. The condition number [Forward error] ~ [Condition number] [Backward error] From problem From computational method A large condition number <=> The problem is sensitive or, ill-conditioned

12. The forward error: 5 -- Ouch! Who’s responsible? The backward error: 5 x 10-10 -- method is good! Conclusion: the problem is “bad”

13. If the answer is highly sensitive to perturbations, you have probably asked the wrong question. Maxims about numerical mathematics, computers, science and life, L. N. Trefethen. SIAM News A: “Customer” B: Numerical analyst Who is asking a wrong question? A: The polynomial B: The computing objective What is the wrong question?

14. The question we used to ask: (Fundamental Theorem of Algebra) Given a polynomial p(x) = xn + a1 xn-1+...+an-1 x + an find z = (z1, ..., zn)such that p(x) = ( x - z1)( x - z2) ... ( x - zn ) Right - or - Wrong ?

15. Kahan’s pejorative manifolds All n-polynomials having certain multiplicity structure form a pejorative manifold xn + a1xn-1+...+an-1 x + an<=> (a1, ..., an-1 , an) Example: ( x-t)2 = x2 + (-2t) x + t2 Pejorative manifold: a1= -2t a2= t2

16. Pejorative manifolds of 3-polynomials ( x - s)( x - t)2 = x3 + (-s-2t) x2 + (2st+t2) x + (-st2) a1= -s-2t a2= 2st+t2 a3= -st2 Pejorative manifold of multiplicity structure [1,2] ( x - s)3 =x3 + (-3s) x2 + (3s2) x + (-s3) a1 = -3s a2 = 3s2 a3 = -s3 Pejorative manifold of multiplicity structure [ 3 ]

17. Pejorative manifolds of 3-polynomials The wings: a1= -s-2t a2= 2st+t2 a3= -st2 The edge: a1 = -3s a2 = 3s2 a3 = -s3 General form of pejorative manifolds u = G(z)

18. W. Kahan, Conserving confluence curbs ill-condition, 1972 1. Ill-condition occurs when a polynomial is near a pejorative manifold. 2. A small “drift” by a polynomial on that pejorative manifold does not cause large forward error to the multiple roots, except 3. If a multiple root is sensitive to small perturbation on the pejorative manifold, then the polynomial is near a pejorative submanifold of higher multiplicity. Ill-condition is caused by solving polynomial equations on a wrong manifold

19. Pejorative manifolds of 3-polynomials The wings: a1= -s-2t a2= 2st+t2 a3= -st2 The edge: a1 = -3s a2 = 3s2 a3 = -s3

20. Given a polynomialp(x) = xn + a1 xn-1+...+an-1 x + an / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / Find( z1, ..., zn ) such that p(x) = ( x - z1)( x - z2) ... ( x - zn) The wrong question: because you are asking for simple roots! Find distinct z1, ..., zm such that p(x) = ( x - z1 ) s1( x - z2 )s2 ... ( x - zm )sm s1+...+ sm = n, m < n The right question: do it on the pejorative manifold!

21. For ill-conditioned polynomial p(x)= xn + a1 xn-1+...+an-1 x + an ~ a = (a1, ..., an-1 , an) The objective:find u*=G(z*) that is nearest to p(x)~a

22. Let ( x - z1 ) s1( x - z2)s2 ... ( x - zm)sm = xn + g1( z1, ..., zm )xn-1+...+gn-1( z1, ..., zm ) x + gn( z1, ..., zm ) Then, p(x) = ( x - z1 ) s1( x - z2)s2 ... ( x - zm)sm <==> g1( z1, ..., zm )=a1 g2( z1, ..., zm )=a2 ... ... ... gn( z1, ..., zm )=an n (m<n) m I.e. An over determined polynomial system G(z) = a

23. The polynomial a Project to tangent plane u1 = G(z0)+J(z0)(z1- z0) ~ tangent plane P0 : u = G(z0)+J(z0)(z- z0) pejorative root u*=G(z*) initial iterate u0=G(z0) new iterate u1=G(z1) Pejorative manifold u = G( z ) Solve G( z ) = a for nonlinear least squares solutionz=z* Solve G(z0)+J(z0)(z - z0 ) = a for linear least squares solutionz =z1 G(z0)+J(z0)(z - z0 ) = a J(z0)(z - z0 ) = - [G(z0) - a ] z1 = z0 - [J(z0)+][G(z0) - a]

24. zi+1=zi - J(zi )+[ G(zi )-a ], i=0,1,2 ... Theorem: If z=(z1, ..., zm) with z1, ..., zm distinct, then the Jacobian J(z) of G(z) is of full rank. Theorem: Let u*=G(z*) be nearest to p(x)~a, if 1. z*=(z*1, ..., z*m) with z*1, ..., z*m distinct; 2. z0 is sufficiently close to z*; 3. a is sufficiently close to u* then the iteration converges with a linear rate. Further assume that a = u* , then the convergence is quadratic.

25. Calculation of G(z)and J(z) ( x - z1 ) s1( x - z2)s2 ... ( x - zm)sm = xn + g1( z1, ..., zm )xn-1+...+gn-1( z1, ..., zm ) x + gn( z1, ..., zm ) Explicit polynomail formula gj( z1, ..., zm ) is neither necessary nor desirable. G(z) and every column of J(z) are calculated via simple numerical procedures, as key components of the algorithm.

26. The conventional condition number Let f(x*) = 0 f(x(e)) + e g(x(e)) =0, x(0) = x* f’(x(e)) x’(e) + g(x(e)) + e g’(x(e)) x’(e) =0 ----- x* ----- x* / / / / / / / / / / / / x’(0) = - g(x*)/ f’(x*) x(e) = x* - e g(x*)/ f’(x*) + h.o.t. | x(e) - x* | < | 1 / f’(x*)| .| e g(x*)| + h.o.t. Condition number | 1 /f’(x*)| = infinity at multiple a root

27. The “pejorative” condition number v = G(z) u = G(y) || u - v ||2 = backward error || y - z ||2 = forward error u - v = G(y) - G(z) = J(z) (y - z) + h.o.t. || u - v ||2 = || J(z) (y - z) ||2 > s || y - z ||2 || y - z ||2 < (1/s) || u - v ||2 1/s is the pejorative condition number where s is the smallest singular value of J(z) .

28. Example 1(x-1.0)100 = 0 To make the problem interesting: round the coefficients to 5 digits Step iterates ... ... ... ... 43 1.007 44 1.001 45 1.0001 46 1.0000003 47 .999999998 conventional condition: infinity pejorative condition: 0.0017 Is it ill-conditioned?

29. Example 2(x-0.9)18(x-1.0)10(x-1.1)16 = 0 Step z1 z2 z3 -------------------------------------------------------------------- 0 .92 .95 1.12 1 .87 1.05 1.10 2 .92 .95 1.11 3 .88 1.01 1.10 4 .90 .97 1.12 5 .901 .992 1.101 6 .89993 .9998 1.1002 7 .9000003 .999998 1.1000007 8 .899999999997 .999999999991 1.100000000009 9 .900000000000006 .99999999999997 1.10000000000001 forward error: 6 x 10-15 backward error: 8 x 10-16 Pejorative condition: 58 Even clustered multiple roots are pejoratively well conditioned

30. If you are on a wrong manifold ... Multiplicity pejorative roots backward error pejorative condition structure ------------------------------------------------------------------------------------------------------------- [1,1,...1] erratic .0000000000000006 1390704851032436. [18,10,16] (.9000, 1.0000, 1.1000) .0000000000000008 58.2 [17,11,16] (.8980, .9934, 1.1006) .0000005 51.9 [14,16,14] (.8890, .9892, 1.1090) .000006 29.0 [10,24,10] (.8711, .9906, 1.1316) .00001 25.8 [ 2,40, 2] (.7393, .9917, 1.3282) .0001 23.6 [ 1,43] (.6559, 1.0053 ) .004 12.9 [44] ( .9925 ) .04 .0058

31. Example 3(x-.3-.6i)100 (x-.1-.7i) 200 (x - .7-.5i) 300 (x-.3-.4i) 400 =0 Scary enough? Round coefficients to 6 digits. Z1 z2 z3 z4 .289 +.601i .100 +.702i .702 +.498i .301 +.399i .309 +.602i .097 +.698i .698 +.499i .299 +.401i .293 +.596i .101 +.7003i .7002 +.5005i .3007 +.4003i .300005 +.600006i .099998 +.6999992i .69999992+.4999993i .2999992 +.3999992i .3000002+.60000005i .09999995+.69999998i .69999997+.49999998i .29999997+.400000002i Roots are correct up to 7 digits! Pejorative condition: 0.58

32. Example 4: The Wilkinson polynomial p(x) = (x-1)(x-2)...(x-20) = x20 - 210 x19 + 20615 x18 + ... There are 605 manifolds in total. It is near some manifolds, but which ones? Multiplicity backward error condition Estimated structure number error ------------------------------------------------------------------------ [1,1,1,1,1,1,1,1,1...,1] .000000000000003 550195997640164 1.6 [1,1,1,1,2,2,2,4,2,2,2] .000000003 29269411 .09 [1,1,1,2,3,4,5,3] .0000001 33563 .003 [1,1,2,3,4,6,3] .000001 6546 .007 [1,1,2,5,7,4] .000005 812 .004 [1,2,5,7,5] .00004 198 .008 [1,3,8,8] .0002 25 .005 [2,8,10] .003 6 .02 [5,15] .04 1 .04 [20] .9 .2 .2 What are the roots of the Wilkinson polynomial? Choose your poison!

33. If the estimated error is the sole criterion... Roots multiplicity ---------------------------- .9999012 * 2.0012667 * 2.9758883 * 4.3848008 ** 6.9548807 *** 10.3585762 **** 15.0276741 ***** 19.2719863 *** The polynomial with these roots matches W. Polynomial for 6 digits The condition 33563 is manageable.

34. Summary of the simultaneous method: - Given a polynomial p(x) ~ a, decide the structure ( x - z1 ) s1( x - z2)s2 ... ( x - zm)sm implicitly define the coefficient operator G(z) and its Jacobian J(z) - Have an initial guess z0 to the “pejorative root” - Start the nonlinear least squares iteration zi+1=zi - J(zi )+[ G(zi )-a], i=0,1,2 ... - If it converges verify (i) the backward error || G(z*)-a|| (ii) the pejorative condition 1/s - Accept the pejorative root if both are tiny

35. Isolated multiple roots: The polynomial p(x)has an m-fold root p(x) = ( x-z*)m . q(x-z*) The clue: p(z*) = p’(z*) = p”(z*) = ... = p(m-1)(z*) = 0 p(m)(z*) =/= 0 t1 p(x) - t1 p(m)(x) t2 p’(x) - t2 p(m)(x) ... ... ... ... F = tm-1 p(m-1)( x) - tm-1 p(m)(x) x p(m-1)( x)

36. t1 p(x) - t1 p(m)(x) t2 p’(x) - t2 p(m)(x) ... ... ... ... F = tm-1 p(m-1)( x) - tm-1 p(m)(x) x p(m-1)( x) t1 t2 ... DF diag(-p(m)(x), -p(m)(x),..., -p(m)(x), p(m)(x) ) = tm-1 x z*is an m-fold root of p <==> (0,0,..., z* ) is a simple zero of F

37. Theorem: Let p(x) = [( x-z*)m . q(x-z*)] + [em-1(x-z*)m-1 + ... + e1(x-z*)+ e0] Then F has a simple zero t1 0 e0 t2 0 e1 ... ... ... = + [1/q(0)] tm-1 0 em-2 z z* em-1/m z = z* + [1/(mq(0))] [em-1]

38. Conclusion 1. Ill-condition is cause by a wrong “identity” 2. Multiple roots are pejoratively well conditioned, thereby tractable. 3. Extension ofmachine precision is not needed, a change in computing concept is. 4. To calculate multiple roots, one has to figure out the multiplicity structure (how?) Question: Can Jordan Canonical form be computed?