Download
# olynomial Roots from Companion Matrix Eigen alues AlanEdelman H PDF document - DocSlides

stefany-barnette | 2014-12-12 | General

### Presentations text content in olynomial Roots from Companion Matrix Eigen alues AlanEdelman H

Show

Page 1

olynomial Roots from Companion Matrix Eigen alues AlanEdelman H.Murak ami Jan uary 1,1994 Abstract In classical linear algebra, the eigen alues of a matrix are sometimes de ned as the ro ots of the c har- acteristic p olynomial. An algorithm to compute the ro ots of a p olynomial b y computing the eigen alues of the corresp onding companion matrix turns the tables on the usual de nition. e deriv e a rst order error analysis of this algorithm that sheds ligh t on b oth the underlying geometry of the problem as w ell as the n umerical error of the algorithm. Our error analysis expands on w ork b y V an Do oren and Dewilde in that it states that the algorithm is bac kw ards norm wise stable in a sense that m ust b e de ned carefully Regarding the stronger concept of a small comp onen wise bac kw ards error, our analysis predicts a small suc h error on a test suite of eigh t random p olynomial s suggested b y T oh and T refethen. Ho ev er, w construct examples for whic h a small comp onen wise relativ bac kw ards error is neither predicted nor obtained in practice. e extend our results to p olynomial matrices, where the result is essen tially the same, but the geometry b ecomes more complicated. In troduction Computing ro ots of p olynomials ma y b e p osed as an eigen alue problem b y forming the companion matrix. The eigen alues of this matrix ma y b e found computing the eigen alues of this nonsymmetric matrix using standard ersions of balancing (v ery imp ortan t for accuracy!!) [7] follo ed b y the QR algorithm as ma y b e found in LAP CK or its precursor EISP CK. This is ho w the MA TLAB command roots p erforms its computation. As Clev e Moler has p oin ted out in [6 ], this metho d ma y not b e the b est p ossible b ecause it uses order storage and order time. An algorithm designed sp eci cally for p olynomial ro ots migh t use order storage and time. And the roundo errors in tro duced corresp ond to p erturbations in the elemen ts of the companion matrix, not the p olynomial co ecien ts. Moler con tin ues b y p oin ting out that e don't kno w of an y cases where the computed ro ots are not the exact ro ots of a sligh tly p erturb ed p olynomial , but suc h cases migh t exist. This pap er in estigates whether suc h cases migh t exist. Let = 1 ; : : : ; n denote the ro ots of ) = : : : that are computed b y this metho d. F urther assume that the are the exact ro ots of ) = ^ + ^ : : : + ^ Departmen of Mathematics Ro om 2-380, Massac usetts Institute of ec hnology Cam bridge, MA 02139, edelman@ma th. mit .ed Supp orted b y NSF gran t DMS-9120852 and the Applied Mathematical Sciences subprogram of the Oce of Energy Researc h, U.S. Departmen t of Energy under Con tract DE-A C03-76SF00098. Quan tum Chemistry Lab., Departmen t of Chemistry , Hokk aido Univ ersit , Sapp oro 060, Japan, hiroshi@ch em2 .ho kud ai .ac .jp

Page 2

What do es it mean to sa y that is a sligh t p erturbation of e giv e four answ ers, the b est one is sa ed for last. In the de nitions, ) is mean t to signify a small though unsp eci ed ultiple of mac hine precision. 1. The \calculus" de nition w ould require that rst order p erturbations of the matrix lead to rst order p erturbations of the co ecien ts. 2. normwise answ er that is compatible with standard eigen alue bac kw ard error analyses is to require that max where denotes the companion matrix corresp onding to , and denotes the mac hine precision. 3. A second norm wise answ er that w ould b e ev en b etter is to require that max balance( Standard algorithms for eigen alue computations balance a matrix y nding a diagonal matrix suc h that C T has a smaller norm than 4. The strongest requiremen t should b e that max This is what w e will mean b y a small omp onentwise p erturb ation If = 0, then one often w an ts ^ to b e 0 to o, i.e., ideally one preserv es the sparsit y structure of the problem. It w as already sho wn b y V an Do oren and Dewilde [11 , p.576] b y a blo c k Gaussian elimination argumen that the calculus de nition holds. Their argumen is alid in the more complicated case of p olynomial matrices. It imm ediately follo ws that go o d norm wise answ ers are ailable, though it is not clear what exactly are the constan ts in olv ed. W e found that in tegrating w ork b y Arnold [1 ] with Newton-lik e iden tities allo ws for an illumi nating geometrical deriv ation of a bac kw ard error b ound that is precise to rst order. Our ork impro es on [11 in that deriv the exact rst order p erturbation expression whic test against umerical exp erimen ts. Numerical exp erimen ts oh and refethen [9 compare this algorithm with the Jenkins-T raub or the Madsen-Reid algorithm. These exp erimen ts indicate that all three algorithms ha roughly similar stabilit prop erties, and further that there app ears to b e close link b et een the pseudosp ectra of the balanced companion matrix and the pseudozero sets of the corresp onding p olynomial. Error Analysis 2.1 ProblemStatemen and Solution Let ) = : : : b e an y monic p olynomial. The companion matrix of this p olynomial (1) has c haracteristic p olynomial det ( z I ) = ).

Page 3

If is dense p erturbation matrix with \small" en tries, the natural error analysis question is the computation of a a : : : a In MA TLAB st yle notation, are studying poly poly ) to rst order. Our result is Theorem 2.1 o rst or der, the c ecient of in is =0 +1 i;i =1 i;i (2) wher is de ne d to b In particular, w e see that a small p erturbation in tro duces errors in the co ecien ts that are linear in the ij Since it is w ell kno wn that standard eigen alue pro cedures compute eigen alues of matrices with a \small" bac kw ard error, w e ha e a precise sense in whic h w e claim that there is a p olynomial near whose exact ro ots are computed in this manner. or con enience, w e state this result in a matrix times v ector format. Let k;d =1 i;i and k;d i;i These are the forw ard and bac kw ard \plus-scans" or \pre xes" of the th diagonal of Our result is that the matrix-v ector pro duct a a a a : : : ;n ;n ;n : : : ;n ;n ;n n; 1) n; 2) n; 3) : : : n; : : : n; = 1 (3) is correct to rst order. The + 1) matrix ab o e con tains bac kw ard pre xes in the lo er triangular part and forw ard pre xes in the upp er triangular part. The last ro w of the matrix equation states that p erturbing the trace of a matrix p erturbs the (negativ e of the) co ecien t of y the same amoun t. If further assume that the is the bac kw ards error matrix computed standard eigen alue routine, then e migh t as ell assume that is nearly upp er triangular. There will also b e elemen ts on the sub diagonal, and p ossibly on the next lo er diagonal, dep ending on exactly whic h eigen alue metho d is used. A result analagous to Theorem 2.1 for matrix p olynomials ma y b e found in Section 4. 2.2 Geometryof angen Spaces and ransv ersalit Figure 1 illustrates matrix space ( ) with the usual F rob enius inner pro duct: A; B tr( AB The basic pla ers are a p olynomial (not sho wn) corresp onding companion matrix Orb manifold of matrices similar to an tangen t space to this manifol d = C X X C Normal normal space to this manifold = ) : is a p olynomial Syl \Sylv ester" family of companion matrices

Page 4

Tan Normal Syl = Companion Matrices n x n Orb Figure 1 Matrix space near a companion matrix The curv ed surface represen ts the orbit of matrices similar to These are all the non-derogatory matrices that ha e eigen alues equal to the ro ots of with the same m ultiplicities. (The dimension of this space turns out to b e An informal argumen t is that only parameters are sp eci ed, i.e., the eigen alues. The tangen space to this manifold, famili ar to an one who has studied elemen tary Lie algebra (or p erformed a rst order calculation), is the set of ommutators C X X C (Dimension .) It is easy to c hec that if is a p olynomial, then an y matrix of the form ) is p erp endicular to ev ery matrix of the form C X X C Only sligh tly more dicult [1] is to v erify that all matrices in the normal space are of the form ). (Dimension = .) The set of companion matrices (also kno wn as the \Sylv ester" family) is ob viously an ane space through (Dimension = .) Prop osition 2.1 The Sylvester sp ac of c omp anion matric es is transv ersal to the tangent sp ac i.e. every matrix may b e expr esse d as a line ar c ombination of a c omp anion matrix and a matrix in the tangent sp ac e. This fact ma y b e found in [1 ]. When w e resolv in to comp onen ts in these t o directions, the former displa ys the c hange in the co ecien ts (to rst order) and the latter do es not a ect the co ecien ts (to rst order). Actually , a stronger statemen t holds: the resolution is unique. In the language of singularit y theory not only do w e ha e transv ersalit , but also universality unique + transv ersal. e explicitly p erform the resolution, and thereb y pro e the prop osition, in the next subsection. In n umerical analysis, there are man y examples where p erturbations in \non-ph ysical" directions cause umerical instabilit In the companion matrix problem, only p erturbations to the p olynomial co ecien ts are relev an t to the problem. Other p erturbations are b y-pro ducts of our n umerical metho d. It is fortunate in this case that the error pro duced b y an en tire dimensional space of extraneous directions is absorb ed y the tangen t space.

Page 5

2.3 Algebra The cen tralizer Let us tak a close lo ok at the cen tralizer of This is the dimensional linear space of p olynomial s in , b ecause is non-derogatory aking the as in (1), for = 0 ; : : : ; n 1 and letting =1, = 0 for [0 ; n ], w e de ne matrices Clearly = 0 ; M and the = 1 ; : : : n span the cen tralizer of An in teresting relationship (whic h is closely related to syn thetic division) is that )( C I =1 [2 , p. 85, Eq. (32)]. The trace of the ab o e equation is the Newton iden tities in generating function form; a ariation on the ab o e equation giv es the exact in erse of a V andermonde matrix [10 ]. A more imp ortan t relationship for our purp oses is that C M +1 I ; = 0 ; : : : from whic h w e can inductiv ely pro e (the easiest w y is bac kw ards from ) that +1 +1 = 1 = 1 +1 = 1 This is almost the T o eplitz matrix with ( i; j ) en try equal to except that the left side of the matrix is lo er triangular and the righ t side is strictly upp er triangular. e are no w ready to resolv e an in to tan syl (4) All w e need is the relationship that expresses ho tan is p erp endicular to the normal space: tr ( tan ) = 0 for = 1 ; : : : ; n: e therefore conclude from (4) that tr( ) = tr ( syl ) = syl k;n (5) Because of the almost T o eplitz nature of the , the trace of in olv es partial sums of elemen ts of along certain diagonals. riting out this trace w e ha e that syl k;n =0 +1 i;i =1 i;i The ab o e expression giv es the co ecien t of the p erturb ed haracteristic p olynomial correct to rst order. (Since the co ecien ts are negated in (1), w e are in terested in syl k;n .) Since syl is zero in ev ery column other than the last, w e ma y also use (4) to calculate tan , should w ho ose.

Page 6

Numerical Experimen ts 3.1 pair of 2x2 examples The companion matrices and for not to o small illustrate man y subtleties that o ccur in oating p oin t arithmetic. or con enience, assume our arithmetic follo ws the IEEE double precision standard for whic 27 is large enough to illustrate our p oin t. o mac hine accuracy , the eigen alues of are 2 and 2 The eigen alues of are 1 and 1. LAP CK and MA TLAB compute 2 and 0 for the eigen alues of , while the eigen alues of are computed to b e 52 and 1. Neither of these matrices are a ected b y balancing. Both of these answ ers are consisten with the error estimate in (2). Neither of these matrices giv es answ ers with a small comp onen wise bac kw ard error. In the rst case, the giv en pro duct of the ro ots is 1 while the exact pro duct of the computed ro ots is 0. In the second case, the giv en sum of the ro ots is 2 81 , while the exact sum of the computed ro ots is 52 Ho ev er, MA TLAB and LAP CK could b e more accurate! Both pac ages compute the Sc ur form of a 2 2 matrix using an algorithm that is more unstable for the smaller eigen alue than is necessary prop ose that suc h high qualit y pac ages should compute the eigen alues of a general 2 2 matrix b y solving the quadratic equation as accurately as p ossible giv en the rounded v alues of the trace and the determinan t. If w e ha e a 2 2 companion matrix, then there will b e no roundo error in the trace and the determinan t. The lesson of these examples is that the ro ots could b e computed far more accurately than w ould b e pre- dicted b y our error b ound (2), but curren tly LAP CK and MA TLAB return eigen alues that are consisten with our b ound. The other lesson is that without further assumptions it is imp ossible to require a small comp onen wise bac kw ard error. ortunately , these examples are rather pathological. As w e will see in the next subsection, in practice w e do exp ect to compute ro ots with a small comp onen wise bac kw ards error. 3.2 more \realistic"set of tests In this subsection use Theorem 2.1 to predict the comp onen wise bac kw ard error. also p erform umerical exp erimen ts to measure the comp onen wise bac kw ard error. Our results sho w that Theorem 2.1 alw ys predicts a small bac kw ard error and is most only p essimistic b y one, t o, or ma yb e three digits. o predict the error, w e m ust mo del the bac kw ards error in tro duced b y the QR algorithm. W e decided to ho ose a mo del that is designed to comp ensate for p essimistic factors often found in rounding error analyses. Therefore, the mo del do es not pro vide a rigorous b ound, but at least in our test cases it seems to w ork w ell in practice. What w e do is consider an error matrix with en tries = 2 52 in all en tries ( i; j ) with 2. or example, when = 6, This structure allo ws (with some o erkill) for the p ossibilit y of double shifting in the eigen alue algorithm. A dense p erturbation matrix did not mak e a substan tial di erence in our test cases.

Page 7

The standard algorithms balance the matrix b y nding a diagonal matrix suc h that AT has a smaller norm than Our mo del will b e that the eigen alue algorithm computes the exact eigen alues of a matrix , where j , i.e. eac h elemen t of has absolute v alue at most ab o e the second lo er diagonal and is otherwise. Therefore are computing the exact eigen alues of T E o rst order, then, the error in the co ecien ts is b ounded b y the absolute v alue of the matrix times the absolute alue of the v ector in the pro duct (3), where the scans are computing using T E T This is ho w w e predict the urther details app ear in our MA TLAB co de in the App endix. ollo wing [9] exactly , w e explore the follo wing degree 20 monic p olynomials: (1) \Wilkinson p olynomial" : zeros 1,2,3,..., 20. (2) the monic p olynomial with zeros equally spaced in the in terv al [ 9]. (3) ) = (20!) 20 =0 =k !. (4) the Bernoulli p olynomial of degree 20. (5) the p olynomial 20 19 18 + 1. (6) the monic p olynomial with zeros 2 10 ; : : : ; (7) the Cheb yshev p olynomial of degree 20. (8) the monic p olynomial with zeros equally spaced on a sine curv e, viz., (2 = 19( + 0 5)) + sin(2 = 19( + 0 5)) ; k 10 ; : : : ; 9. Our exp erimen tal results consist of three columns for eac h p olynomial. T o b e precise, w e rst computed the co ecien ts either exactly or with 30 decimal precision using Mathematica. e then read these n um b ers in to MA TLAB whic h ma y not ha e rounded the last bit or t o correctly Though w e could ha e rounded the result correctly in Mathematica, w e c hose not to do so. Rather w e to ok the rounded p olynomial s stored in MA TLAB to b e our \ocial" test cases. Column 1: Log Predicted Error: In the rst column w e mo del the predicted error from (2) in the manner e describ ed ab o e. First w e compute the mo deled a , and then displa y the rounded v alue of log 10 a =a Column 2: Log Computed Error: e computed the eigen alues using MA TLAB, and then translated the IEEE double precision n um b ers to Mathematica without an y error. e then computed the exact p olynomial with these computed ro ots and compared the relativ e bac kw ard error in the co ecien ts. Column 3: essimism Index: By taking the ratio of the computed error in column 2 to the predicted error describ ed in column and then taking the logarithm and rounding, obtain p essimism index. Indices suc h as 0,-1, and -2 indicate that w e are p essimistic b y at most t o orders of magnitude; and index of -19 indicates a p essimism of 19 orders of magnitude. (Since w e are using negativ e n um b ers, p erhaps w should more prop erly call this an optimism index.) There are no en tries where the p olynomial co ecien t is zero. (The Bernoulli p olynomial is a little funn in that it has a 19 term, but no other o dd degree term.) The computed relativ e error w ould b e in nite in man y of these cases. The top co ecien is the log relativ e error in the determinan t, i.e. the co ecien of the constan t term; the b ottom co ecien t refers to the trace, i.e., the co ecien t of 19 ry en tering 1.000000000 000 000 18 1.00000000 000 000 019 , and 1.00000000 000 00 001 899 999 999 in to MA TLAB ( fteen zeros in eac h n um b er). The results that w e obtained on a SUN Sparc Station 10 w ere 1, 1 + 2 52 , and 1 52 resp ectiv ely , though the correctly rounded result should b e 1 + 2 52 in all instances. Clev e Moler has resp onded that a b etter string parser is needed.

Page 8

Relativ e errors in the Co ecien ts (log base 10 for eigh t di eren t degree 20 p olynomials Key to eac h panel in the table b elo Col 1 Col 2 Col 3 Predicted Observ ed essimism Rel Error Rel Error Index (1) (2) (3) (4) (5) (6) (7) (8) det -12 -13 -1 -12 -12 -13 -14 -1 -13 -14 -1 -14 -14 -1 -8 -14 -6 -12 -14 -2 -13 -14 -2 -12 -14 -2 -13 -14 -1 -13 -14 -1 -13 -14 -1 -8 -14 -5 -12 -13 -2 -12 -12 -13 -14 -1 -13 -15 -2 -13 -14 -1 -9 -14 -5 -12 -14 -2 -13 -15 -2 -12 -13 -2 -13 -14 -1 -13 -14 -1 -13 -16 -3 -9 -14 -5 -12 -15 -3 -12 -14 -2 -13 -14 -2 -13 -14 -1 -13 -14 -1 -9 -14 -5 -12 -14 -2 -13 -15 -2 -12 -16 -4 -13 -14 -1 -13 -14 -1 -13 -14 -1 -10 -14 -4 -12 -14 -2 -13 -14 -1 -13 -14 -1 -13 -14 -1 -13 -14 -1 -10 -15 -5 -13 -14 -1 -13 -15 -2 -12 -15 -2 -13 -15 -2 -13 -14 -1 -13 -14 -1 -10 -14 -4 -12 -15 -2 -13 -15 -2 -13 -15 -2 -13 -14 -1 -13 -14 -1 -11 -15 -5 -13 -14 -1 -13 -15 -1 -13 -15 -3 -13 -15 -2 -13 -16 -3 -13 -14 -1 -11 -14 -3 -13 -15 -2 -13 -14 -1 -13 -15 -2 -13 -14 -1 -13 -14 -1 -11 -14 -3 -13 -15 -1 -13 -15 -1 10 -13 -15 -2 -13 -15 -1 -13 -15 -1 -13 -14 -1 -11 -14 -3 11 -13 -15 -3 -13 -14 -1 -13 -15 -1 -13 -14 -1 -13 -14 -1 -12 -15 -4 -14 -15 -1 -14 -15 -1 12 -13 -14 -1 -13 -14 -1 -13 -14 -1 -13 -14 -1 -12 -15 -3 13 -13 -14 -1 -13 -14 -1 -13 -14 -1 -14 -14 -1 -13 -14 -1 -12 -15 -2 -14 -15 -1 -14 -14 -1 14 -13 -15 -2 -13 -14 -1 -13 -14 -1 -13 -14 -1 -13 -14 -2 15 -13 -15 -2 -14 -14 -1 -13 -14 -1 -13 -15 -2 -13 -15 -1 -13 -16 -3 -14 -15 -1 -14 -15 -1 16 -13 -15 -2 -14 -14 -1 -13 -14 -1 -13 -15 -2 -13 -15 -2 17 -14 -16 -2 -14 -15 -1 -14 -14 -1 -14 -16 -3 -13 -14 -1 -14 -15 -1 -14 -15 -1 -14 -15 -1 18 trace -14 -16 -2 -14 -14 -14 -14 -1 -14 -16 -2 -14 -14 -14 -16 -2 19 In all cases, see that the computed bac kw ard relativ error as excellen t. With the exception of column (6), this is fairly w ell predicted usually with a p essimism index of one to three orders of magnitude. Column (6) is an exception, where the bac kw ard error is far more fa orable than w e predict. Wh this migh t b e p ossible as explained in the previous subsection. kno the determinan to full precision (v ery rare when p erforming matrix computations!!). So long as our QR shifts and de ation criteria manage not to destro y the determinan t, it will remain in tact. In column (6) w e see a case where this o ccurred. By con trast, in the previous subsection w e illustrated a case where this did not o ccur. e susp ect that for an y companion matrix, it is often if not alw ys p ossible to c ho ose shifts and de ation criteria to guaran tee high (bac kw ard relativ e) accuracy ev en with the smallest of determinan ts, but w e ha not pro en this. Generalization to Matrix olynomials A monic matrix p olynomial has the form ) = : : : where assume that the co ecien ts (and are matrices, and is scalar. It is of in terest [3, 4, 5 ] to nd the for whic h det ( )) = 0 and the corresp onding eigen ectors suc that

Page 9

) = 0. Suc h information ma y b e obtained from the pn pn blo c k companion matrix (6) The Sylv ester space no w has dimension np while the tangen space to the orbit of generically has dimension np , though it can b e smaller. It seems that if p > 1, w e ha e to o man y dimensions! will no w sho w that w e ma y pro ceed in a manner that is analogous to that of Section 2.3 to obtain what is roughly the same answ er, but to do so w e m ust carefully pic k a natural subspace of the tangen t space that will giv e a univ ersal decomp osition. This is not necessary when = 1. The natural subspace of the tangen space consists of all matrices of the form C X X C where the last ro ws of are 0. Lemma 4.1 De ne wher denotes the Kr one cker pr duct and is the identity matrix of or der Then C M +1 and +1 +1 +1 Pro of These statemen ts are readily v eri ed b y induction. e no w in tro duce the blo c k trace of a matrx: De nitio n 4.1 If is a pn pn matrix whose blo cks ar e denote ij , then we de ne tr =1 ii Notice that tr ) is a matrix, not a scalar. Theorem 4.1 Given the rst blo ck c olumns of a pn pn matrix , the c ondition that 0 = tr Z M = 0 ; : : : ; n (7) is e quivalent to the c ondition that C X X C for some with ottom blo ck r ow (8) Mor over, either c ondition determines the nal blo ck c olumn of uniquely given the r emaining c olumns.

Page 10

Pro of The n; k blo c en try of is and this determines kn uniquely from (7). If has as its b ottom blo c k ro w, it is easy to v erify that the map from to the rst 1 blo c k columns of C X X C has a trivial n ullspace. Th us is uniquely determined b y (8). What remains is to sho w (8) implies (7). Supp ose that has a zero b ottom blo c k ro w. Then tr C Y ) = tr Y C ) = =1 i;i +1 Therefore, if has a zero b ottom blo c k ro w, then tr C X M ) = tr X M ) b ho osing X M Finally tr X C X C )) 0 for an matrix b ecause ) is 0 except for the last blo c k column. Therefore X M X C M , from whic h w conclude that tr C X M ) = tr X M ) = tr X C M ). e no w summarize the geometry Corollary 4.1 In the dimensional sp ac e of np np matric es, the np subsp ac e of the tangent sp ac e of the orbit of de ne d either by (7) or (8) is tr ansver al at to the np dimensional Sylvester sp ac onsisting of blo ck c omp anion matric es. e ma y no w resolv e an y p erturbation matrix in to tan syl (9) as in (4), except no tan ust b e in this np dimensional subspace, and syl is 0 except in the nal blo ck column. Because tr Z M = tr ), the correct result is that syl k;n =0 +1 i;i =1 i;i The ab o e expression giv es the blo c k co ecien t of a p erturb ed matrix p olynomial correct to rst order. 10

Page 11

Ac kno wledgemen ts e w ould lik e to thank Jim Demmel for man y in teresting discussions and exp erimen ts on transv ersalit y and p olynomial ro ot nding. References [1] V.I. Arnold, On matrices dep ending on parameters, RussianMath.Surv eys26 (1971), 29{43. [2] F.R. Gan tmac her, TheTheoryofMatrices , V ols. 1, 2, Chelsea, New Y ork, 1959. [3] I. Goh b erg, P . Lancaster, and L. Ro dman, MatrixP olynomial , Academic Press, 1982. [4] D. Mano c ha and J. Demmel, Algorithms for in tersecting parametric and algebraic curv es I: simple in ter- sections, CMT ransactionsonGraphics,13 , (1994), 73{100. [5] D. Mano c ha and J. Demmel, Algorithms for in tersecting parametric and algebraic curv es I I: ultiple in tersections, Computergraphics,visionandimageprocessing:Graphicalmodelsandimageprocessing to app ear. [6] C. Moler, Clev e's Corner: OOTS Of p olynomials, that is, The MathW orks Newsletter, v. n. (Spring 1991), 8{9. [7] B. P arlett and C. Reinsc h, Balancing a matrix for calculation of eigen alues and eigen ectors, Numer. Math.,13 , (1969), 293{304. [8] W.H. Press, et al., NumericalRecipes , Cam bridge Univ ersit y Press, 1986. [9] K. T oh and L.N. T refethen, Pseudozeros of p olynomials and pseudosp ectra of companion matrices, Cornell Univ ersit y Departmen t of Computer Science preprin t, TR 93-1360, June 1993 (to app ear in Numerisc he Mathematik). [10] J.F. raub, Asso ciated p olynomials and uniform metho ds for the solution of linear problems, SIAM Review,8 , (1966), 277-301. [11] . V an Do oren and P . Dewilde, The eigenstructure of an arbitrary p olynomial matrix: computational asp ects, LinearAlg.Appl.50 , (1983), 545{579. 11

Page 12

MA TLAB program used in experimen ts The umerical part of our exp erimen ts as p erformed in MA TLAB, while the exact comp onen of the exp erimen ts w as p erformed with Mathematica. W e repro duce only the main MA TLAB co de b elo w for the purp ose of sp ecifying precisely the crucial comp onen ts of our exp erimen ts. The subroutines scan and scanr compute the plus-scan and the rev ersed plus-scan of a v ector resp ectiv ely Predict and compute the componen twi sebackward error in the eight polynomialtest cases suggestedby Toh and Trefethen --- Alan Edelman, October 1993 Step ... Run MathematicaProgram to ComputeCoefficient s. Output will be read into matlab as the array d1 consistingof eight columns and 21 rows from the constant term (det) on top to the x^19 term (trace), then on the bottom (Code not shown.) Step ... Form the eight companionmatrices for i=1:8, eval(['m'num2str(i)'=zeros(20 ); ']) eval(['m'num2str(i)'(2:21:380 )= one s(1 ,1 9); '] ); eval(['m'num2str(i)'(:,20)=-d 1( 1:2 0,'num2str(i ),');']); end Step ... Obtain an "unbalance d"error matrix. e=eps*ones (2 0); e= tri u(e ,- 2); for i=1:8, eval(['[t, b]= ba lan ce( m'num2str(i ');']); ti=diag(1. /di ag (t) ); eval(['e'num2str(i)'=(t*e*ti) *n orm (b) ;' ]); end Step ... Compute the first order perturbat io matrix: er. for j=1:8, eval(['e=e num2str(j )';']); forw=zeros (20 ); bac k=z er os( 20 );e r= zer os( 20 ,21 ); for i=1:19 forw forw diag(scan( di ag( e, i-1 )), i- 1); back back diag(scanr (d iag (e ,-i )), -i ); end er(1:19,1: 19) =b ack (2: 20 ,1: 19 );e r( :,2 :21 )= er( :, 2:2 1)- fo rw; eval(['er'num2str(j)'=er;']); end Step ... Compute the predictedrelative errors in the coefficie nt predicted= zeros(20,8) for i=1:8, eval(['pr ed ict ed (:, num2str( i)')=abs(er'num2str(i)')*abs(d1 (:, i) );' ]) end predicted= ab s(p re dic ted ./ d1( 1: 20, :) ); Step ... Compute the exact relative errors in the coefficien ts using Mathemati caand finally display all relevant quantitiesby taking the base 10 logarithm. (Code not shown.) 12

Murak ami Jan uary 11994 Abstract In classical linear algebra the eigen alues of a matrix are sometimes de ned as the ro ots of the c har acteristic p olynomial An algorithm to compute the ro ots of a p olynomial b y computing the eigen alues of the ID: 23002

- Views :
**160**

**Direct Link:**- Link:https://www.docslides.com/stefany-barnette/olynomial-roots-from-companion
**Embed code:**

Download this pdf

DownloadNote - The PPT/PDF document "olynomial Roots from Companion Matrix Ei..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.

Page 1

olynomial Roots from Companion Matrix Eigen alues AlanEdelman H.Murak ami Jan uary 1,1994 Abstract In classical linear algebra, the eigen alues of a matrix are sometimes de ned as the ro ots of the c har- acteristic p olynomial. An algorithm to compute the ro ots of a p olynomial b y computing the eigen alues of the corresp onding companion matrix turns the tables on the usual de nition. e deriv e a rst order error analysis of this algorithm that sheds ligh t on b oth the underlying geometry of the problem as w ell as the n umerical error of the algorithm. Our error analysis expands on w ork b y V an Do oren and Dewilde in that it states that the algorithm is bac kw ards norm wise stable in a sense that m ust b e de ned carefully Regarding the stronger concept of a small comp onen wise bac kw ards error, our analysis predicts a small suc h error on a test suite of eigh t random p olynomial s suggested b y T oh and T refethen. Ho ev er, w construct examples for whic h a small comp onen wise relativ bac kw ards error is neither predicted nor obtained in practice. e extend our results to p olynomial matrices, where the result is essen tially the same, but the geometry b ecomes more complicated. In troduction Computing ro ots of p olynomials ma y b e p osed as an eigen alue problem b y forming the companion matrix. The eigen alues of this matrix ma y b e found computing the eigen alues of this nonsymmetric matrix using standard ersions of balancing (v ery imp ortan t for accuracy!!) [7] follo ed b y the QR algorithm as ma y b e found in LAP CK or its precursor EISP CK. This is ho w the MA TLAB command roots p erforms its computation. As Clev e Moler has p oin ted out in [6 ], this metho d ma y not b e the b est p ossible b ecause it uses order storage and order time. An algorithm designed sp eci cally for p olynomial ro ots migh t use order storage and time. And the roundo errors in tro duced corresp ond to p erturbations in the elemen ts of the companion matrix, not the p olynomial co ecien ts. Moler con tin ues b y p oin ting out that e don't kno w of an y cases where the computed ro ots are not the exact ro ots of a sligh tly p erturb ed p olynomial , but suc h cases migh t exist. This pap er in estigates whether suc h cases migh t exist. Let = 1 ; : : : ; n denote the ro ots of ) = : : : that are computed b y this metho d. F urther assume that the are the exact ro ots of ) = ^ + ^ : : : + ^ Departmen of Mathematics Ro om 2-380, Massac usetts Institute of ec hnology Cam bridge, MA 02139, edelman@ma th. mit .ed Supp orted b y NSF gran t DMS-9120852 and the Applied Mathematical Sciences subprogram of the Oce of Energy Researc h, U.S. Departmen t of Energy under Con tract DE-A C03-76SF00098. Quan tum Chemistry Lab., Departmen t of Chemistry , Hokk aido Univ ersit , Sapp oro 060, Japan, hiroshi@ch em2 .ho kud ai .ac .jp

Page 2

What do es it mean to sa y that is a sligh t p erturbation of e giv e four answ ers, the b est one is sa ed for last. In the de nitions, ) is mean t to signify a small though unsp eci ed ultiple of mac hine precision. 1. The \calculus" de nition w ould require that rst order p erturbations of the matrix lead to rst order p erturbations of the co ecien ts. 2. normwise answ er that is compatible with standard eigen alue bac kw ard error analyses is to require that max where denotes the companion matrix corresp onding to , and denotes the mac hine precision. 3. A second norm wise answ er that w ould b e ev en b etter is to require that max balance( Standard algorithms for eigen alue computations balance a matrix y nding a diagonal matrix suc h that C T has a smaller norm than 4. The strongest requiremen t should b e that max This is what w e will mean b y a small omp onentwise p erturb ation If = 0, then one often w an ts ^ to b e 0 to o, i.e., ideally one preserv es the sparsit y structure of the problem. It w as already sho wn b y V an Do oren and Dewilde [11 , p.576] b y a blo c k Gaussian elimination argumen that the calculus de nition holds. Their argumen is alid in the more complicated case of p olynomial matrices. It imm ediately follo ws that go o d norm wise answ ers are ailable, though it is not clear what exactly are the constan ts in olv ed. W e found that in tegrating w ork b y Arnold [1 ] with Newton-lik e iden tities allo ws for an illumi nating geometrical deriv ation of a bac kw ard error b ound that is precise to rst order. Our ork impro es on [11 in that deriv the exact rst order p erturbation expression whic test against umerical exp erimen ts. Numerical exp erimen ts oh and refethen [9 compare this algorithm with the Jenkins-T raub or the Madsen-Reid algorithm. These exp erimen ts indicate that all three algorithms ha roughly similar stabilit prop erties, and further that there app ears to b e close link b et een the pseudosp ectra of the balanced companion matrix and the pseudozero sets of the corresp onding p olynomial. Error Analysis 2.1 ProblemStatemen and Solution Let ) = : : : b e an y monic p olynomial. The companion matrix of this p olynomial (1) has c haracteristic p olynomial det ( z I ) = ).

Page 3

If is dense p erturbation matrix with \small" en tries, the natural error analysis question is the computation of a a : : : a In MA TLAB st yle notation, are studying poly poly ) to rst order. Our result is Theorem 2.1 o rst or der, the c ecient of in is =0 +1 i;i =1 i;i (2) wher is de ne d to b In particular, w e see that a small p erturbation in tro duces errors in the co ecien ts that are linear in the ij Since it is w ell kno wn that standard eigen alue pro cedures compute eigen alues of matrices with a \small" bac kw ard error, w e ha e a precise sense in whic h w e claim that there is a p olynomial near whose exact ro ots are computed in this manner. or con enience, w e state this result in a matrix times v ector format. Let k;d =1 i;i and k;d i;i These are the forw ard and bac kw ard \plus-scans" or \pre xes" of the th diagonal of Our result is that the matrix-v ector pro duct a a a a : : : ;n ;n ;n : : : ;n ;n ;n n; 1) n; 2) n; 3) : : : n; : : : n; = 1 (3) is correct to rst order. The + 1) matrix ab o e con tains bac kw ard pre xes in the lo er triangular part and forw ard pre xes in the upp er triangular part. The last ro w of the matrix equation states that p erturbing the trace of a matrix p erturbs the (negativ e of the) co ecien t of y the same amoun t. If further assume that the is the bac kw ards error matrix computed standard eigen alue routine, then e migh t as ell assume that is nearly upp er triangular. There will also b e elemen ts on the sub diagonal, and p ossibly on the next lo er diagonal, dep ending on exactly whic h eigen alue metho d is used. A result analagous to Theorem 2.1 for matrix p olynomials ma y b e found in Section 4. 2.2 Geometryof angen Spaces and ransv ersalit Figure 1 illustrates matrix space ( ) with the usual F rob enius inner pro duct: A; B tr( AB The basic pla ers are a p olynomial (not sho wn) corresp onding companion matrix Orb manifold of matrices similar to an tangen t space to this manifol d = C X X C Normal normal space to this manifold = ) : is a p olynomial Syl \Sylv ester" family of companion matrices

Page 4

Tan Normal Syl = Companion Matrices n x n Orb Figure 1 Matrix space near a companion matrix The curv ed surface represen ts the orbit of matrices similar to These are all the non-derogatory matrices that ha e eigen alues equal to the ro ots of with the same m ultiplicities. (The dimension of this space turns out to b e An informal argumen t is that only parameters are sp eci ed, i.e., the eigen alues. The tangen space to this manifold, famili ar to an one who has studied elemen tary Lie algebra (or p erformed a rst order calculation), is the set of ommutators C X X C (Dimension .) It is easy to c hec that if is a p olynomial, then an y matrix of the form ) is p erp endicular to ev ery matrix of the form C X X C Only sligh tly more dicult [1] is to v erify that all matrices in the normal space are of the form ). (Dimension = .) The set of companion matrices (also kno wn as the \Sylv ester" family) is ob viously an ane space through (Dimension = .) Prop osition 2.1 The Sylvester sp ac of c omp anion matric es is transv ersal to the tangent sp ac i.e. every matrix may b e expr esse d as a line ar c ombination of a c omp anion matrix and a matrix in the tangent sp ac e. This fact ma y b e found in [1 ]. When w e resolv in to comp onen ts in these t o directions, the former displa ys the c hange in the co ecien ts (to rst order) and the latter do es not a ect the co ecien ts (to rst order). Actually , a stronger statemen t holds: the resolution is unique. In the language of singularit y theory not only do w e ha e transv ersalit , but also universality unique + transv ersal. e explicitly p erform the resolution, and thereb y pro e the prop osition, in the next subsection. In n umerical analysis, there are man y examples where p erturbations in \non-ph ysical" directions cause umerical instabilit In the companion matrix problem, only p erturbations to the p olynomial co ecien ts are relev an t to the problem. Other p erturbations are b y-pro ducts of our n umerical metho d. It is fortunate in this case that the error pro duced b y an en tire dimensional space of extraneous directions is absorb ed y the tangen t space.

Page 5

2.3 Algebra The cen tralizer Let us tak a close lo ok at the cen tralizer of This is the dimensional linear space of p olynomial s in , b ecause is non-derogatory aking the as in (1), for = 0 ; : : : ; n 1 and letting =1, = 0 for [0 ; n ], w e de ne matrices Clearly = 0 ; M and the = 1 ; : : : n span the cen tralizer of An in teresting relationship (whic h is closely related to syn thetic division) is that )( C I =1 [2 , p. 85, Eq. (32)]. The trace of the ab o e equation is the Newton iden tities in generating function form; a ariation on the ab o e equation giv es the exact in erse of a V andermonde matrix [10 ]. A more imp ortan t relationship for our purp oses is that C M +1 I ; = 0 ; : : : from whic h w e can inductiv ely pro e (the easiest w y is bac kw ards from ) that +1 +1 = 1 = 1 +1 = 1 This is almost the T o eplitz matrix with ( i; j ) en try equal to except that the left side of the matrix is lo er triangular and the righ t side is strictly upp er triangular. e are no w ready to resolv e an in to tan syl (4) All w e need is the relationship that expresses ho tan is p erp endicular to the normal space: tr ( tan ) = 0 for = 1 ; : : : ; n: e therefore conclude from (4) that tr( ) = tr ( syl ) = syl k;n (5) Because of the almost T o eplitz nature of the , the trace of in olv es partial sums of elemen ts of along certain diagonals. riting out this trace w e ha e that syl k;n =0 +1 i;i =1 i;i The ab o e expression giv es the co ecien t of the p erturb ed haracteristic p olynomial correct to rst order. (Since the co ecien ts are negated in (1), w e are in terested in syl k;n .) Since syl is zero in ev ery column other than the last, w e ma y also use (4) to calculate tan , should w ho ose.

Page 6

Numerical Experimen ts 3.1 pair of 2x2 examples The companion matrices and for not to o small illustrate man y subtleties that o ccur in oating p oin t arithmetic. or con enience, assume our arithmetic follo ws the IEEE double precision standard for whic 27 is large enough to illustrate our p oin t. o mac hine accuracy , the eigen alues of are 2 and 2 The eigen alues of are 1 and 1. LAP CK and MA TLAB compute 2 and 0 for the eigen alues of , while the eigen alues of are computed to b e 52 and 1. Neither of these matrices are a ected b y balancing. Both of these answ ers are consisten with the error estimate in (2). Neither of these matrices giv es answ ers with a small comp onen wise bac kw ard error. In the rst case, the giv en pro duct of the ro ots is 1 while the exact pro duct of the computed ro ots is 0. In the second case, the giv en sum of the ro ots is 2 81 , while the exact sum of the computed ro ots is 52 Ho ev er, MA TLAB and LAP CK could b e more accurate! Both pac ages compute the Sc ur form of a 2 2 matrix using an algorithm that is more unstable for the smaller eigen alue than is necessary prop ose that suc h high qualit y pac ages should compute the eigen alues of a general 2 2 matrix b y solving the quadratic equation as accurately as p ossible giv en the rounded v alues of the trace and the determinan t. If w e ha e a 2 2 companion matrix, then there will b e no roundo error in the trace and the determinan t. The lesson of these examples is that the ro ots could b e computed far more accurately than w ould b e pre- dicted b y our error b ound (2), but curren tly LAP CK and MA TLAB return eigen alues that are consisten with our b ound. The other lesson is that without further assumptions it is imp ossible to require a small comp onen wise bac kw ard error. ortunately , these examples are rather pathological. As w e will see in the next subsection, in practice w e do exp ect to compute ro ots with a small comp onen wise bac kw ards error. 3.2 more \realistic"set of tests In this subsection use Theorem 2.1 to predict the comp onen wise bac kw ard error. also p erform umerical exp erimen ts to measure the comp onen wise bac kw ard error. Our results sho w that Theorem 2.1 alw ys predicts a small bac kw ard error and is most only p essimistic b y one, t o, or ma yb e three digits. o predict the error, w e m ust mo del the bac kw ards error in tro duced b y the QR algorithm. W e decided to ho ose a mo del that is designed to comp ensate for p essimistic factors often found in rounding error analyses. Therefore, the mo del do es not pro vide a rigorous b ound, but at least in our test cases it seems to w ork w ell in practice. What w e do is consider an error matrix with en tries = 2 52 in all en tries ( i; j ) with 2. or example, when = 6, This structure allo ws (with some o erkill) for the p ossibilit y of double shifting in the eigen alue algorithm. A dense p erturbation matrix did not mak e a substan tial di erence in our test cases.

Page 7

The standard algorithms balance the matrix b y nding a diagonal matrix suc h that AT has a smaller norm than Our mo del will b e that the eigen alue algorithm computes the exact eigen alues of a matrix , where j , i.e. eac h elemen t of has absolute v alue at most ab o e the second lo er diagonal and is otherwise. Therefore are computing the exact eigen alues of T E o rst order, then, the error in the co ecien ts is b ounded b y the absolute v alue of the matrix times the absolute alue of the v ector in the pro duct (3), where the scans are computing using T E T This is ho w w e predict the urther details app ear in our MA TLAB co de in the App endix. ollo wing [9] exactly , w e explore the follo wing degree 20 monic p olynomials: (1) \Wilkinson p olynomial" : zeros 1,2,3,..., 20. (2) the monic p olynomial with zeros equally spaced in the in terv al [ 9]. (3) ) = (20!) 20 =0 =k !. (4) the Bernoulli p olynomial of degree 20. (5) the p olynomial 20 19 18 + 1. (6) the monic p olynomial with zeros 2 10 ; : : : ; (7) the Cheb yshev p olynomial of degree 20. (8) the monic p olynomial with zeros equally spaced on a sine curv e, viz., (2 = 19( + 0 5)) + sin(2 = 19( + 0 5)) ; k 10 ; : : : ; 9. Our exp erimen tal results consist of three columns for eac h p olynomial. T o b e precise, w e rst computed the co ecien ts either exactly or with 30 decimal precision using Mathematica. e then read these n um b ers in to MA TLAB whic h ma y not ha e rounded the last bit or t o correctly Though w e could ha e rounded the result correctly in Mathematica, w e c hose not to do so. Rather w e to ok the rounded p olynomial s stored in MA TLAB to b e our \ocial" test cases. Column 1: Log Predicted Error: In the rst column w e mo del the predicted error from (2) in the manner e describ ed ab o e. First w e compute the mo deled a , and then displa y the rounded v alue of log 10 a =a Column 2: Log Computed Error: e computed the eigen alues using MA TLAB, and then translated the IEEE double precision n um b ers to Mathematica without an y error. e then computed the exact p olynomial with these computed ro ots and compared the relativ e bac kw ard error in the co ecien ts. Column 3: essimism Index: By taking the ratio of the computed error in column 2 to the predicted error describ ed in column and then taking the logarithm and rounding, obtain p essimism index. Indices suc h as 0,-1, and -2 indicate that w e are p essimistic b y at most t o orders of magnitude; and index of -19 indicates a p essimism of 19 orders of magnitude. (Since w e are using negativ e n um b ers, p erhaps w should more prop erly call this an optimism index.) There are no en tries where the p olynomial co ecien t is zero. (The Bernoulli p olynomial is a little funn in that it has a 19 term, but no other o dd degree term.) The computed relativ e error w ould b e in nite in man y of these cases. The top co ecien is the log relativ e error in the determinan t, i.e. the co ecien of the constan t term; the b ottom co ecien t refers to the trace, i.e., the co ecien t of 19 ry en tering 1.000000000 000 000 18 1.00000000 000 000 019 , and 1.00000000 000 00 001 899 999 999 in to MA TLAB ( fteen zeros in eac h n um b er). The results that w e obtained on a SUN Sparc Station 10 w ere 1, 1 + 2 52 , and 1 52 resp ectiv ely , though the correctly rounded result should b e 1 + 2 52 in all instances. Clev e Moler has resp onded that a b etter string parser is needed.

Page 8

Relativ e errors in the Co ecien ts (log base 10 for eigh t di eren t degree 20 p olynomials Key to eac h panel in the table b elo Col 1 Col 2 Col 3 Predicted Observ ed essimism Rel Error Rel Error Index (1) (2) (3) (4) (5) (6) (7) (8) det -12 -13 -1 -12 -12 -13 -14 -1 -13 -14 -1 -14 -14 -1 -8 -14 -6 -12 -14 -2 -13 -14 -2 -12 -14 -2 -13 -14 -1 -13 -14 -1 -13 -14 -1 -8 -14 -5 -12 -13 -2 -12 -12 -13 -14 -1 -13 -15 -2 -13 -14 -1 -9 -14 -5 -12 -14 -2 -13 -15 -2 -12 -13 -2 -13 -14 -1 -13 -14 -1 -13 -16 -3 -9 -14 -5 -12 -15 -3 -12 -14 -2 -13 -14 -2 -13 -14 -1 -13 -14 -1 -9 -14 -5 -12 -14 -2 -13 -15 -2 -12 -16 -4 -13 -14 -1 -13 -14 -1 -13 -14 -1 -10 -14 -4 -12 -14 -2 -13 -14 -1 -13 -14 -1 -13 -14 -1 -13 -14 -1 -10 -15 -5 -13 -14 -1 -13 -15 -2 -12 -15 -2 -13 -15 -2 -13 -14 -1 -13 -14 -1 -10 -14 -4 -12 -15 -2 -13 -15 -2 -13 -15 -2 -13 -14 -1 -13 -14 -1 -11 -15 -5 -13 -14 -1 -13 -15 -1 -13 -15 -3 -13 -15 -2 -13 -16 -3 -13 -14 -1 -11 -14 -3 -13 -15 -2 -13 -14 -1 -13 -15 -2 -13 -14 -1 -13 -14 -1 -11 -14 -3 -13 -15 -1 -13 -15 -1 10 -13 -15 -2 -13 -15 -1 -13 -15 -1 -13 -14 -1 -11 -14 -3 11 -13 -15 -3 -13 -14 -1 -13 -15 -1 -13 -14 -1 -13 -14 -1 -12 -15 -4 -14 -15 -1 -14 -15 -1 12 -13 -14 -1 -13 -14 -1 -13 -14 -1 -13 -14 -1 -12 -15 -3 13 -13 -14 -1 -13 -14 -1 -13 -14 -1 -14 -14 -1 -13 -14 -1 -12 -15 -2 -14 -15 -1 -14 -14 -1 14 -13 -15 -2 -13 -14 -1 -13 -14 -1 -13 -14 -1 -13 -14 -2 15 -13 -15 -2 -14 -14 -1 -13 -14 -1 -13 -15 -2 -13 -15 -1 -13 -16 -3 -14 -15 -1 -14 -15 -1 16 -13 -15 -2 -14 -14 -1 -13 -14 -1 -13 -15 -2 -13 -15 -2 17 -14 -16 -2 -14 -15 -1 -14 -14 -1 -14 -16 -3 -13 -14 -1 -14 -15 -1 -14 -15 -1 -14 -15 -1 18 trace -14 -16 -2 -14 -14 -14 -14 -1 -14 -16 -2 -14 -14 -14 -16 -2 19 In all cases, see that the computed bac kw ard relativ error as excellen t. With the exception of column (6), this is fairly w ell predicted usually with a p essimism index of one to three orders of magnitude. Column (6) is an exception, where the bac kw ard error is far more fa orable than w e predict. Wh this migh t b e p ossible as explained in the previous subsection. kno the determinan to full precision (v ery rare when p erforming matrix computations!!). So long as our QR shifts and de ation criteria manage not to destro y the determinan t, it will remain in tact. In column (6) w e see a case where this o ccurred. By con trast, in the previous subsection w e illustrated a case where this did not o ccur. e susp ect that for an y companion matrix, it is often if not alw ys p ossible to c ho ose shifts and de ation criteria to guaran tee high (bac kw ard relativ e) accuracy ev en with the smallest of determinan ts, but w e ha not pro en this. Generalization to Matrix olynomials A monic matrix p olynomial has the form ) = : : : where assume that the co ecien ts (and are matrices, and is scalar. It is of in terest [3, 4, 5 ] to nd the for whic h det ( )) = 0 and the corresp onding eigen ectors suc that

Page 9

) = 0. Suc h information ma y b e obtained from the pn pn blo c k companion matrix (6) The Sylv ester space no w has dimension np while the tangen space to the orbit of generically has dimension np , though it can b e smaller. It seems that if p > 1, w e ha e to o man y dimensions! will no w sho w that w e ma y pro ceed in a manner that is analogous to that of Section 2.3 to obtain what is roughly the same answ er, but to do so w e m ust carefully pic k a natural subspace of the tangen t space that will giv e a univ ersal decomp osition. This is not necessary when = 1. The natural subspace of the tangen space consists of all matrices of the form C X X C where the last ro ws of are 0. Lemma 4.1 De ne wher denotes the Kr one cker pr duct and is the identity matrix of or der Then C M +1 and +1 +1 +1 Pro of These statemen ts are readily v eri ed b y induction. e no w in tro duce the blo c k trace of a matrx: De nitio n 4.1 If is a pn pn matrix whose blo cks ar e denote ij , then we de ne tr =1 ii Notice that tr ) is a matrix, not a scalar. Theorem 4.1 Given the rst blo ck c olumns of a pn pn matrix , the c ondition that 0 = tr Z M = 0 ; : : : ; n (7) is e quivalent to the c ondition that C X X C for some with ottom blo ck r ow (8) Mor over, either c ondition determines the nal blo ck c olumn of uniquely given the r emaining c olumns.

Page 10

Pro of The n; k blo c en try of is and this determines kn uniquely from (7). If has as its b ottom blo c k ro w, it is easy to v erify that the map from to the rst 1 blo c k columns of C X X C has a trivial n ullspace. Th us is uniquely determined b y (8). What remains is to sho w (8) implies (7). Supp ose that has a zero b ottom blo c k ro w. Then tr C Y ) = tr Y C ) = =1 i;i +1 Therefore, if has a zero b ottom blo c k ro w, then tr C X M ) = tr X M ) b ho osing X M Finally tr X C X C )) 0 for an matrix b ecause ) is 0 except for the last blo c k column. Therefore X M X C M , from whic h w conclude that tr C X M ) = tr X M ) = tr X C M ). e no w summarize the geometry Corollary 4.1 In the dimensional sp ac e of np np matric es, the np subsp ac e of the tangent sp ac e of the orbit of de ne d either by (7) or (8) is tr ansver al at to the np dimensional Sylvester sp ac onsisting of blo ck c omp anion matric es. e ma y no w resolv e an y p erturbation matrix in to tan syl (9) as in (4), except no tan ust b e in this np dimensional subspace, and syl is 0 except in the nal blo ck column. Because tr Z M = tr ), the correct result is that syl k;n =0 +1 i;i =1 i;i The ab o e expression giv es the blo c k co ecien t of a p erturb ed matrix p olynomial correct to rst order. 10

Page 11

Ac kno wledgemen ts e w ould lik e to thank Jim Demmel for man y in teresting discussions and exp erimen ts on transv ersalit y and p olynomial ro ot nding. References [1] V.I. Arnold, On matrices dep ending on parameters, RussianMath.Surv eys26 (1971), 29{43. [2] F.R. Gan tmac her, TheTheoryofMatrices , V ols. 1, 2, Chelsea, New Y ork, 1959. [3] I. Goh b erg, P . Lancaster, and L. Ro dman, MatrixP olynomial , Academic Press, 1982. [4] D. Mano c ha and J. Demmel, Algorithms for in tersecting parametric and algebraic curv es I: simple in ter- sections, CMT ransactionsonGraphics,13 , (1994), 73{100. [5] D. Mano c ha and J. Demmel, Algorithms for in tersecting parametric and algebraic curv es I I: ultiple in tersections, Computergraphics,visionandimageprocessing:Graphicalmodelsandimageprocessing to app ear. [6] C. Moler, Clev e's Corner: OOTS Of p olynomials, that is, The MathW orks Newsletter, v. n. (Spring 1991), 8{9. [7] B. P arlett and C. Reinsc h, Balancing a matrix for calculation of eigen alues and eigen ectors, Numer. Math.,13 , (1969), 293{304. [8] W.H. Press, et al., NumericalRecipes , Cam bridge Univ ersit y Press, 1986. [9] K. T oh and L.N. T refethen, Pseudozeros of p olynomials and pseudosp ectra of companion matrices, Cornell Univ ersit y Departmen t of Computer Science preprin t, TR 93-1360, June 1993 (to app ear in Numerisc he Mathematik). [10] J.F. raub, Asso ciated p olynomials and uniform metho ds for the solution of linear problems, SIAM Review,8 , (1966), 277-301. [11] . V an Do oren and P . Dewilde, The eigenstructure of an arbitrary p olynomial matrix: computational asp ects, LinearAlg.Appl.50 , (1983), 545{579. 11

Page 12

MA TLAB program used in experimen ts The umerical part of our exp erimen ts as p erformed in MA TLAB, while the exact comp onen of the exp erimen ts w as p erformed with Mathematica. W e repro duce only the main MA TLAB co de b elo w for the purp ose of sp ecifying precisely the crucial comp onen ts of our exp erimen ts. The subroutines scan and scanr compute the plus-scan and the rev ersed plus-scan of a v ector resp ectiv ely Predict and compute the componen twi sebackward error in the eight polynomialtest cases suggestedby Toh and Trefethen --- Alan Edelman, October 1993 Step ... Run MathematicaProgram to ComputeCoefficient s. Output will be read into matlab as the array d1 consistingof eight columns and 21 rows from the constant term (det) on top to the x^19 term (trace), then on the bottom (Code not shown.) Step ... Form the eight companionmatrices for i=1:8, eval(['m'num2str(i)'=zeros(20 ); ']) eval(['m'num2str(i)'(2:21:380 )= one s(1 ,1 9); '] ); eval(['m'num2str(i)'(:,20)=-d 1( 1:2 0,'num2str(i ),');']); end Step ... Obtain an "unbalance d"error matrix. e=eps*ones (2 0); e= tri u(e ,- 2); for i=1:8, eval(['[t, b]= ba lan ce( m'num2str(i ');']); ti=diag(1. /di ag (t) ); eval(['e'num2str(i)'=(t*e*ti) *n orm (b) ;' ]); end Step ... Compute the first order perturbat io matrix: er. for j=1:8, eval(['e=e num2str(j )';']); forw=zeros (20 ); bac k=z er os( 20 );e r= zer os( 20 ,21 ); for i=1:19 forw forw diag(scan( di ag( e, i-1 )), i- 1); back back diag(scanr (d iag (e ,-i )), -i ); end er(1:19,1: 19) =b ack (2: 20 ,1: 19 );e r( :,2 :21 )= er( :, 2:2 1)- fo rw; eval(['er'num2str(j)'=er;']); end Step ... Compute the predictedrelative errors in the coefficie nt predicted= zeros(20,8) for i=1:8, eval(['pr ed ict ed (:, num2str( i)')=abs(er'num2str(i)')*abs(d1 (:, i) );' ]) end predicted= ab s(p re dic ted ./ d1( 1: 20, :) ); Step ... Compute the exact relative errors in the coefficien ts using Mathemati caand finally display all relevant quantitiesby taking the base 10 logarithm. (Code not shown.) 12

Today's Top Docs

Related Slides