Structured matrix methods for polynomial computations Joab R. - PowerPoint PPT Presentation

✬ ✩ Structured matrix methods for polynomial computations Joab R. Winkler Department of Computer Science The University of Sheffield United Kingdom Structured Matrix Days, Limoges, France. May 2012 ✫ ✪

✬ ✩ Contents • Introduction – Difficulties of performing reliable computations on polynomials – The condition number of a root of a polynomial • Multiple roots • The geometry of multiple roots • A simple polynomial root finder – GCD computations – Polynomial division • Examples • Summary ✫ ✪ 1

✬ ✩ 1. INTRODUCTION Problems in which it is required to compute the greatest common divisor (GCD) of two polynomials, and the roots of polynomials: • Computer-aided geometric design: The calculation of the points of intersection of curves and surfaces • Control theory • Cryptography ✫ ✪ 2

✬ ✩ 1.1 Difficulties of performing reliable computations on polynomials Commonly used methods and algorithms for computing the roots of a polynomial include: • Bairstow, Graeffe, Jenkins-Traub, Laguerre, M¨ uller, Newton, QR decomposition Commonly used software for computing the roots of a polynomial include: • M ATHEMATICA , M ATLAB , NAG, MPS OLVE Satisfactory results are obtained if: • The polynomial is of moderate degree • The roots are simple and well-separated • A good starting point in the iterative scheme is used ✫ ✪ 3

✬ ✩ Example 1.1 Consider the polynomial y 4 − 4 y 3 + 6 y 2 − 4 y + 1 = ( y − 1) 4 whose root is y = 1 with multiplicity 4 . M ATLAB returns the roots 1.0000 ± 0.0002 i , � 1.0002, 0.9998 Example 1.2 The roots of the polynomial ( y − 1) 100 were computed by M ATLAB . 4 Figure 1.1: The 3 2 computed roots of 1 Imag ( y − 1) 100 . 0 −1 −2 −3 −4 0 1 2 3 4 5 6 Real � ✫ ✪ 4

✬ ✩ Example 1.3 The computation of multiple roots in the presence of noise. f(x) = (x−0.6) 3 (x−1) 5 ε c = 1e−007 f(x) = (x−0.7) 3 (x−1) 5 ε c = 1e−007 0.08 0.1 0.06 0.04 0.05 0.02 Imag Imag 0 0 −0.02 −0.04 −0.05 −0.06 −0.08 −0.1 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 0.7 0.8 0.9 1 1.1 Real Real f(x) = (x−0.8) 3 (x−1) 5 ε c = 1e−007 f(x) = (x−0.9) 3 (x−1) 5 ε c = 1e−007 0.15 0.15 0.1 0.1 0.05 0.05 0 0 Imag Imag −0.05 −0.05 −0.1 −0.1 −0.15 −0.15 −0.2 −0.2 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 0.8 0.9 1 1.1 1.2 1.3 Real Real ✫ ✪ � Figure 1.2: The distribution of the roots of four perturbed polynomials. 5

✬ ✩ Example 1.4 Consider the polynomial ( y − 1) 100 and perturb the constant term by 2 − 100 g ( y ) = ( y − 1) 100 − 2 − 100 The roots of g ( y ) are y k = 1 + 1 � � 2 kπ � � 2 kπ �� cos + i sin , k = 0 , . . . , 99 2 100 100 A change in the constant coefficient of 2 − 100 ≈ 10 − 32 causes a change in the � magnitude of the roots from 1 to 1.5. ✫ ✪ 6

✬ ✩ 2. MULTIPLE ROOTS Consider the polynomial m � a i y m − i f ( y ) = i =0 that has a root of multiplicity r f ( k ) ( y ) = 0 , f (0) ( y ) = f ( y ) , k = 0 , . . . , r − 1 • Why are multiple roots ill -conditioned? • Are there conditions for which they are well-conditioned, and can they be utilised in a stable algorithm for their computation? ✫ ✪ 7

✬ ✩ Example 2.1 componentwise noise amplitude componentwise signal amplitude = 10 − 8 ε c = exact exact computed computed relative root mult. root mult. error -6.7547e-01 4 -6.7547000082e-01 4 1.2139725913e-09 5.7335e+00 6 5.7335000822e+00 6 1.4344694923e-08 2.1747e+00 7 2.1746999924e+00 7 3.5069931355e-09 -9.5568e+00 10 -9.5567996740e+00 10 3.4111255034e-08 -6.5553e+00 11 -6.5553001701e+00 11 2.5954947075e-08 � ✫ ✪ 8

✬ ✩ Example 2.2 componentwise noise amplitude componentwise signal amplitude = 10 − 8 ε c = exact exact computed computed relative root mult. root mult. error -1.1539e+00 4 -1.1539000316e+00 4 2.7389266674e-08 4.0809e+00 5 4.0808998890e+00 5 2.7198581384e-08 -2.1059e+00 6 -2.1058999294e+00 6 3.2521823798e-08 3.6683e+00 7 3.6683000481e+00 7 1.3110114060e-08 -9.6084e+00 13 -9.6084001121e+00 13 1.1664214687e-08 � ✫ ✪ 9

✬ ✩ 3. THE GEOMETRY OF MULTIPLE ROOTS William Kahan (UC Berkeley) showed it is necessary to distinguish between a random perturbation, and a structured perturbation, applied to the coefficients of a polynomial f ( y ) when considering the numerical condition of its roots. Let the root y i of f ( y ) have multiplicity m i : m r r a i y m − i = � � � ( y − y i ) m i , f ( y ) = m i = m i =0 i =1 i =1 Consider the perturbed polynomials g ( y ) and h ( y ) : m � ( a i + δa i ) y m − i , g ( y ) = | δa i | ≤ ε | a i | i =0 r ( y − ( y i + δy i )) m i , � h ( y ) = | δy i | ≪ | y i | i =1 ✫ ✪ 10

✬ ✩ • Consider the polynomial g ( y ) m � ( a i + δa i ) y m − i , g ( y ) = | δa i | ≤ ε | a i | i =0 If the roots of g ( y ) are y i + δy i , then the roots y i are ill -conditioned. ∆ y i ∆ y i = | δy i | ≫ 1 , | y i | , i = 1 , . . . , m ε • Consider the polynomial h ( y ) r ( y − ( y i + δy i )) m i , � h ( y ) = | δy i | ≪ | y i | i =1 Then ∆ y i ∆ y i = | δy i | ∆ f = || f − h � ∆ f = 0(1) , | y i | , i = 1 , . . . , r, � f � ✫ ✪ 11

✬ ✩ Summary: • If random perturbations are applied to the coefficients of a polynomial f ( y ) , then its multiple roots are ill -conditioned and break up into complex conjugate simple roots. • If structured perturbations are applied to the coefficients of f ( y ) , such that the multiplicities of its roots are preserved, then its roots are stable. Kahan introduced the pejorative manifold of f ( y ) : r m r ( y − y i ) m i = � � � a i y m − i , f ( y ) = m i = m i =1 i =0 i =1 • The pejorative manifold of f ( y ) defines the constraints satisfied by the coeffi- cients a i such that the multiplicities of its roots are m i . ✫ ✪ 12

✬ ✩ Example 3.1 Consider a quartic polynomial f ( y ) with real roots y 1 , y 2 , y 3 and y 4 : y 4 − ( y 1 + y 2 + y 3 + y 4 ) y 3 + f ( y ) = ( y 1 y 2 + y 1 y 3 + y 1 y 4 + y 2 y 3 + y 2 y 4 + y 3 y 4 ) y 2 − ( y 1 y 2 y 3 + y 1 y 2 y 4 + y 1 y 3 y 4 + y 2 y 3 y 4 ) y + y 1 y 2 y 3 y 4 • If f ( y ) has a quartic root, then y 1 = y 2 = y 3 = y 4 , then f ( y ) = y 4 − 4 y 1 y 3 + 6 y 2 1 y 2 − 4 y 3 1 y + y 4 1 The pejorative manifold of this quartic polynomial is � � 6 y 2 − 4 y 3 y 4 − 4 y 1 1 1 1 ✫ ✪ 13

✬ ✩ • If f ( y ) has one treble root and one simple root, then y 1 = y 2 = y 3 � = y 4 , then f ( y ) = y 4 − (3 y 1 + y 4 ) y 3 + 3( y 2 1 + y 1 y 4 ) y 2 − ( y 3 1 + 3 y 2 1 y 4 ) y + y 3 1 y 4 The pejorative manifold of this quartic polynomial is � � 3( y 2 − ( y 3 1 + 3 y 2 y 3 − (3 y 1 + y 4 ) 1 + y 1 y 4 ) 1 y 4 ) 1 y 4 • If f ( y ) has two double roots, then y 1 = y 2 and y 3 = y 4 � = y 1 , then y 4 − 2( y 1 + y 4 ) y 3 + ( y 2 1 + 4 y 1 y 4 + y 2 4 ) y 2 f ( y ) = − 2( y 2 1 y 4 + y 1 y 2 4 ) y + y 2 1 y 2 4 The pejorative manifold of this quartic polynomial is � � ( y 2 1 + 4 y 1 y 4 + y 2 − 2( y 2 1 y 4 + y 1 y 2 y 2 1 y 2 − 2( y 1 + y 4 ) 4 ) 4 ) 4 � ✫ ✪ 14

✬ ✩ 4. A SIMPLE POLYNOMIAL ROOT FINDER Gauss developed a polynomial root finder that is fundamentally different from the methods mentioned earlier. Let w 1 ( y ) be the product of all linear factors of f ( y ) Let w 2 ( y ) be the product of all quadratic factors of f ( y ) · · · Let w i ( y ) be the product of all factors of degree i of f ( y ) f ( y ) = w 1 ( y ) w 2 2 ( y ) w 3 3 ( y ) · · · w r r ( y ) ✫ ✪ 15

✬ ✩ Perform a sequence of GCD computations: f ( y ) , f (1) ( y ) w 2 ( y ) w 2 3 ( y ) w 3 4 ( y ) · · · w r − 1 � � q 1 ( y ) = = ( y ) GCD r � � q 1 ( y ) , q (1) w 3 ( y ) w 2 4 ( y ) w 3 5 ( y ) · · · w r − 2 q 2 ( y ) = 1 ( y ) = ( y ) GCD r � � q 2 ( y ) , q (1) w 4 ( y ) w 2 5 ( y ) w 3 6 ( y ) · · · w r − 3 q 3 ( y ) = 2 ( y ) = ( y ) GCD r � � q 3 ( y ) , q (1) w 5 ( y ) w 2 6 ( y ) w 3 7 ( y ) · · · w r − 4 q 4 ( y ) = 3 ( y ) = ( y ) GCD r . . . The sequence terminates at q r ( y ) , which is a constant ✫ ✪ 16

✬ ✩ A set of polynomials h i ( y ) , i = 1 , . . . , r, is defined such that f ( y ) h 1 ( y ) = = w 1 ( y ) w 2 ( y ) w 3 ( y ) · · · q 1 ( y ) q 1 ( y ) h 2 ( y ) = = w 2 ( y ) w 3 ( y ) · · · q 2 ( y ) q 2 ( y ) h 3 ( y ) = = w 3 ( y ) · · · q 3 ( y ) . . . q r − 1 ( y ) h r ( y ) = = w r ( y ) q r ( y ) The functions, w 1 ( y ) , w 2 ( y ) , · · · , w r ( y ) , are determined from w 1 ( y ) = h 1 ( y ) w 2 ( y ) = h 2 ( y ) w r − 1 ( y ) = h r − 1 ( y ) h 2 ( y ) , h 3 ( y ) , · · · , h r ( y ) until w r ( y ) = h r ( y ) ✫ ✪ 17

✬ ✩ The equations w 1 ( y ) = 0 , w 2 ( y ) = 0 , · · · , w r ( y ) = 0 contain only simple roots, and they yield the simple, double, triple, etc. roots of f ( y ) • If y 0 is a root of w i ( y ) , then it is a root of multiplicity i of f ( y ) � • This is a divide -and-conquer algorithm. • The multiplicities of the roots are defined by the values of the subscripts i of the polynomials w i ( y ) . • There are two essential operations – GCD computations – Polynomial divisions ✫ ✪ but these operations are ill-posed. 18