✬ ✩ Structured matrix methods for the computation of multiple roots of univariate polynomials John D. Allan and Joab R. Winkler Department of Computer Science The University of Sheffield, United Kingdom ✫ ✪
✬ ✩ Overview • Introduction • Structured and unstructured condition numbers • A root finder that computes root multiplicities • The Sylvester resultant matrix • The method of structured total least norm • The principle of minimum description length • A polynomial root solver • Conclusions ✫ ✪ 1
✬ ✩ 1. Introduction Classical methods yield unsatisfactory results for the computation of multiple roots of a polynomial. Example Consider the fourth order polynomial x 4 − 4 x 3 + 6 x 2 − 4 x + 1 = ( x − 1) 4 whose root is x = 1 with multiplicity 4 . The roots function in MATLAB returns 1 . 0002 , 1 . 0000 ± 0 . 0002 i, 0 . 9998 which shows that roundoff errors are sufficient to cause a relative error in the solution of 2 × 10 − 4 . � ✫ ✪ 2
✬ ✩ Example The roots of the polynomial ( x − 1) 100 were computed by the roots function in MATLAB. 4 3 2 1 Imag 0 −1 −2 −3 −4 0 1 2 3 4 5 6 Real Figure 1: The computed roots of ( x − 1) 100 . ✫ ✪ 3
✬ ✩ Example Consider the computed roots of ( x − 1) n , n = 1 , 6 , 11 , . . . , when the constant term is perturbed by 2 − 10 . 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0 0.5 1 1.5 2 Figure 2: The computed roots of ( x − 1) n , n = 1 , 6 , 11 , . . . , when the constant term is perturbed. ✫ ✪ 4
✬ ✩ 2. Structured and unstructured condition numbers The componentwise and normwise condition numbers of a root x 0 of the polynomial m � f ( x ) = a i φ i ( x ) i =0 are � 1 � n m 1 1 m ! � κ c ( x 0 ) = | a i φ i ( x 0 ) | � � � f ( m ) ( x 0 ) 1 − 1 | x 0 | � ε m i =0 c and � 1 � m 1 1 m ! κ n ( x 0 ) = � � a � 2 � φ ( x 0 ) � 2 � � � f ( m ) ( x 0 ) 1 − 1 | x 0 | ε m n ✫ ✪ 5
✬ ✩ Compare with the structured condition number of ( x − x 0 ) n that preserves the multiplicity of the root, but not the value of the root � 1 � � n � n � 2 ( x 0 ) 2 i 2 � ( x − x 0 ) n � 2 ρ ( x 0 ) = ∆ x 0 1 1 i =0 i ∆ f = = � 2 ( x 0 ) 2 i � n − 1 � ( x − x 0 ) n − 1 � 2 � n − 1 n | x 0 | n | x 0 | i =0 i where ∆ f = � δf � 2 ∆ x 0 = | δx 0 | and � f � 2 | x 0 | ✫ ✪ 6
✬ ✩ Example The condition number ρ (1) of the root x 0 = 1 of the polynomial ( x − 1) n is � � 1 � � n � � n � 2 � 2 n � � 2 ρ (1) = 1 = 1 � = 1 2(2 n − 1) ≈ 2 � i =0 i n � � 2( n − 1) � 2 � n − 1 � n − 1 n n n n n n − 1 i =0 i if n is large. Compare with the componentwise and normwise condition numbers for random (un - structured) perturbations κ c (1) ≈ | δx 0 | κ n (1) ≈ | δx 0 | and ε c ε n � ✫ ✪ 7
✬ ✩ Conclusion • A multiple root of a polynomial is ill-conditioned, and the ill-conditioning increases as its multiplicity increases • If the multiplicity is sufficiently large, its condition numbers are proportional to the signal-to-noise ratio • A multiple root of a polynomial is well-conditioned with respect to perturbations that preserve its multiplicity ✫ ✪ 8
✬ ✩ 3. A root finder that computes root multiplicities If there is prior evidence that the theoretically exact form of the polynomial has multi- ple roots, then it is desirable to compute these multiple roots • Given an inexact polynomial, what is the nearest polynomial that has multiple roots? • How can multiple roots be computed because roundoff error is sufficient to cause them to break up into simple roots? ✫ ✪ 9
✬ ✩ Example Let w i ( x ) , i = 1 , . . . , m max , be the product of all factors of degree i of f ( x ) f ( x ) = w 1 ( x ) w 2 2 ( x ) w 3 3 ( x ) · · · w m max m max ( x ) and thus � � f ( x ) , f (1) ( x ) = w 2 ( x ) w 2 3 ( x ) w 3 4 ( x ) · · · w m max − 1 r 0 ( x ) = GCD ( x ) m max Similarly � � r 0 ( x ) , r (1) w 3 ( x ) w 2 4 ( x ) w 3 5 ( x ) · · · w m max − 2 r 1 ( x ) = 0 ( x ) = ( x ) GCD m max � � r 1 ( x ) , r (1) w 4 ( x ) w 2 5 ( x ) w 3 6 ( x ) · · · w m max − 3 r 2 ( x ) = 1 ( x ) = ( x ) GCD m max � � r 2 ( x ) , r (1) w 5 ( x ) w 2 6 ( x ) · · · w m max − 4 r 3 ( x ) = 2 ( x ) = ( x ) GCD m max . . . ✫ and the sequence terminates at r m max − 1 ( x ) , which is a constant. ✪ 10
✬ ✩ A sequence of polynomials h i ( x ) , i = 1 , . . . , m max , is defined such that f ( x ) h 1 ( x ) = = w 1 ( x ) w 2 ( x ) w 3 ( x ) · · · r 0 ( x ) r 0 ( x ) h 2 ( x ) = = w 2 ( x ) w 3 ( x ) · · · r 1 ( x ) r 1 ( x ) h 3 ( x ) = = w 3 ( x ) · · · r 2 ( x ) . . . r m max − 2 h m max ( x ) = = w m max ( x ) r m max − 1 and thus all the functions, w 1 ( x ) , w 2 ( x ) , · · · , w m max ( x ) , are determined from w 1 ( x ) = h 1 ( x ) w 2 ( x ) = h 2 ( x ) , w m max − 1 ( x ) = h m max − 1 ( x ) h 2 ( x ) , h 3 ( x ) , · · · h m max ( x ) until w m max ( x ) = h m max ( x ) ✫ ✪ 11
✬ ✩ The equations w 1 ( x ) = 0 , w 2 ( x ) = 0 , · · · , w m max ( x ) = 0 contain only simple roots: In particular • All the roots of w i are roots of multiplicity i of f ( x ) � The algorithm requires repeated GCD and polynomial division computations: • The GCD computations are achieved by using the Sylvester matrix and its sub- resultant matrices • The polynomial division computations are implemented by a deconvolution oper- ation ✫ ✪ 12
✬ ✩ 4. The Sylvester resultant matrix The Sylvester resultant matrix S ( f, g ) of the polynomials m n � � a i x m − i b j x n − j f ( x ) = g ( x ) = and i =0 j =0 has two important properties: • The determinant of S ( f, g ) is equal to zero if and only if f ( x ) and g ( x ) have a non-constant common divisor • The rank of S ( f, g ) is equal to ( m + n − d ) , where d is the degree of the greatest common divisor (GCD) of f ( x ) and g ( x ) ✫ ✪ 13
✬ ✩ The matrix S ( f, g ) is � � a 0 b 0 � � � a 1 a 0 b 1 b 0 � � . . ... ... . . � a 1 b 1 . . � � . . ... ... � . . a m − 1 a 0 b n − 1 b 0 . . � � � ... ... � a m a m − 1 a 1 b n b n − 1 b 1 � � . . ... ... . . � a m b n . . � � ... ... � a m − 1 b n − 1 � � � a m b n � ✫ ✪ 14
✬ ✩ The k ’th Sylvester matrix, or subresultant, S k ∈ R ( m + n − k +1) × ( m + n − 2 k +2) is a submatrix of S ( f, g ) that is formed by: • Deleting the last ( k − 1) rows of S ( f, g ) • Deleting the last ( k − 1) columns of the coefficients of f ( x ) • Deleting the last ( k − 1) columns of the coefficients of g ( x ) The integer k satisfies 1 ≤ k ≤ min ( m, n ) , and a subresultant matrix is defined for each value of k . • Start with k = k 0 = min ( m, n ) and decrease by one until a solution is obtained. ✫ ✪ 15
✬ ✩ Each matrix S k is partitioned into: • A vector c k ∈ R m + n − k +1 , where c k is the first column of S k • A matrix A k ∈ R ( m + n − k +1) × ( m + n − 2 k +1) , where A k is the matrix formed from the remaining columns of S k � � � � = S k � A k c k � � � � � � = c k coeffs. of f ( x ) coeffs. of g ( x ) � � ���� � �� � � �� � 1 n − k m − k + 1 ✫ ✪ 16
✬ ✩ Theorem 1 Consider the polynomials f ( x ) and g ( x ) , and let k be a positive integer, where 1 ≤ k ≤ min ( m, n ) . Then 1. The dimension of the null space of S k is greater than or equal to one if and only if the over determined equation A k y = c k possesses a solution. 2. A necessary and sufficient condition for the polynomials f ( x ) and g ( x ) to have a common divisor of degree greater than or equal to k is that the dimension of the null space of S k is greater than or equal to one. ✫ ✪ 17
✬ ✩ 5. The method of structured total least norm The problem to be solved is: For a given value of k , compute the smallest perturbations to f ( x ) and g ( x ) such that ( A k + E k ) y = c k + h k has a solution, where • E k has the same structure as A k • h k has the same structure as c k • Start with k = min ( m, n ) and decrease by one until a solution is obtained. ✫ ✪ 18
✬ ✩ • z i be the perturbation of a i , i = 0 , . . . , m , of f ( x ) • z m +1+ j be the perturbation of b j , j = 0 , . . . , n, of g ( x ) The perturbation of the Sylvester resultant matrix is: z 0 z m +1 z 1 z 0 z m +2 . . ... ... . . z 1 . . . ... ... . z m − 1 z 0 z m + n z m +1 � . � � � = � E k h k ... ... z m z m − 1 z 1 z m + n +1 z m +2 . . ... ... . . z m . . ... ... z m − 1 z m + n z m z m + n +1 ✫ ✪ 19
Recommend
More recommend
