'1 I m U) r7 5:5 lllllllllllllllllllllllllllllllllllllllllllllllllllllllll 31293 01823 210 This is to certify that the dissertation entitled Date/i“ m1 mi m6 elve— lTa V‘CIAM/ WWW {Wm if 5. WW. presented by 4"“? l {ow/xjva WWI/Cg has been accepted towards fulfillment of the requirements for P l’\ . 0 degree in WAW§ 1% %fi—vé7 Majoerrofessor Date % [1;] n /é ¢7 MSU is an Affirmative Action/Equal Opportuniry Institution 0-12771 LIBRARY Michigan State University PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 1/98 chlRC/DateDue.p65-p.14 DETERMINING THE JORDAN NORMAL FORM OF A MATRIX By Tianjun Wang A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics 1999 ABSTRACT DETERMINING THE JORDAN NORMAL FORM OF A MATRIX By Tianj un Wang In this dissertation, we present a new algorithm for calculating the Jordan normal form of an n X n matrix. The algorithm first determines the structure of the Jordan blocks of a matrix without solving for its eigenvalues, and then calculates eigenvalues of the matrix, as zeros of its minimal polynomials, to fill in the diagonal entries of the Jordan blocks. The approximation of eigenvalues of the matrix has no influence on the structure. There are two crucial steps for determining the structure of the Jordan blocks of a matrix in our algorithm. One is to determine the minimal polynomials of cyclic subspaces of the matrix by using a Las Vegas type algorithm with probability one for obtaining the correct answer. The other is to ascertain the structures of the minimal polynomials. The most important feature of our algorithm is that the im- plementations of these two steps can all be carried out by the symbolic computation, which consists of exact rational operations. Therefore, if the entries of a given matrix are all rational, then the structure of the Jordan blocks in its Jordan normal form can be determined exactly without approximating its eigenvalues. Copyright © by Tianjun Wang 1999 To my family iv ACKNOWLEDGMENTS I am most indebted to my dissertation advisor, University Distinguished Professor Tien—Yien Li, for his support during my graduate study at Michigan State University. I would like to thank him for suggesting the problem and giving invaluable guidance which makes this dissertation a reality. Special thanks go to Professor Zhinan Zhang at Xinjiang University for his dis- cussions, suggestions and collaborations when he visited Michigan State University. I would like to thank my dissertation committee members Professor Chichia Chiu, Professor Richard Hill, Professor Jay Kurtz and Professor Zhengfang Zhou for their valuable suggestions. My greatest thanks are reserved for my wife, Diane D. Dong, and my family for their enduring patience, constant support and encouragement. TABLE OF CONTENTS INTRODUCTION 1 PRELIMINARIES 2 THEORETICAL BACKGROUND 2.1 Main theorems ................................ 2.2 A constructive proof of Theorem 2.1.1 ................... 2.3 A constructive proof of Theorem 2.1.2 ................... 3 THE ALGORITHM 4 IMPLEMENTATION OF THE ALGORITHM 4.1 Implementation of Step A .......................... 4.2 Implementation of Step B .......................... 4.3 Implementation of Step C .......................... 4.4 Computational result ............................. APPENDIX BIBLIOGRAPHY vi INTRODUCTION The numerical computation for the Jordan normal form of a given matrix can be traced back to the 1970’s, e.g., see [3], [4], [9] and [11]. Usually, the singular value decomposition (SVD) method is used to ascertain the structure of the Jordan blocks and compute the Jordan normal form of a matrix. The difficulty of this method lies in a critical distinction between equal and unequal eigenvalues for computing the rank and determining finally the block structures. In numerical computation, it is, in general, difficult to determine whether a matrix has exactly equal eigenvalues. Since the sensitivity of computing the rank of a matrix subject to perturbations degrades the reliability in determining the rank, the method can sometimes alter the structure of the Jordan normal form radically[8]. Recently, the computation of the Jordan normal form of a matrix has been studied in several articles. The algebraic theory of matrices is used in [13] to present a parallel algorithm for ascertaining the structures of the Jordan blocks of a matrix. The Jordan normal form of a matrix consists of two parts: the structure of the Jordan blocks and the eigenvalues of the matrix. In this dissertation, we present a new constructive proof of a fundamental decomposition theorem in [5] which provides the geometric theory of elementary divisors of a matrix. Based on the new proof of the theorem, we develop an algorithm for calculating the Jordan normal form of an n x 12 matrix. The algorithm first determines the structure of the Jordan blocks of a matrix without solving for its eigenvalues, and then calculates eigenvalues of the matrix, as zeros of its minimal polynomials, to fill in the diagonal entries of the Jordan blocks. The approximation of eigenvalues of the matrix has no influence on the structure. There are two crucial steps for determining the structure of the Jordan blocks of a matrix in our algorithm. The first one is to determine the minimal polynomials of cyclic subspaces of the matrix by using a Las Vegas type algorithm with probability one for obtaining the correct answer. The other is to ascertain the structures of those minimal polynomials. The most important feature of our algorithm is that the implementations of these two steps can all be carried out by symbolic computation, which consists of exact rational operations. Therefore, if the entries of a given matrix are all rational, then the structure of the Jordan blocks in its Jordan normal form can be determined exactly without approximating its eigenvalues. The dissertation is organized as follows. In Chapter 1, we introduce some of the basic concepts, results and notations that will be utilized repeatedly in subsequent chapters. In Chapter 2, we study the theory of determining the Jordan normal form of a matrix, and give a new constructive proof of Theorem 3 on page 187 in [5]. In Chap- ter 3, we develop an algorithm of determining the Jordan normal form of a matrix based on the theoretical results obtained in Chpater 2. In Chapter 4, the imple- mentation of the algorithm presented in Chapter 3 is discussed in detail, meanwhile, computational result is provided to illustrate the effectiveness of the algorithm. CHAPTER 1 PRELIMINARIES We use R" to denote the n-dimensional vector space over the real number field R, and 7?,[A] is the polynomial ring on R. Let L be a linear operator mapping 7?," to R”, then L can be represented by an n x n matrix A with elements in R under a basis of R". Conversely, an n X n matrix A with elements in ’R can also be regarded as a linear operator from R" to R” under the standard basis 81,62, . . . ,e,, of ’R", where e, =(1,0,...,0)T, e2=(0,1,0,...,0)T, en: (0,...,0,1)T. Definition 1.0.1 A subspace S of R” is invariant under A (or is an invariant sub- space of A) if Ax E 5' whenever x E S. For x 52$ 0 in 72", consider the sequence x, Ax, Azx, ..., Alx, (1.1) Since R” is an n-dimensional vector space, these vectors defined in (1.1) cannot all be linearly independent, so there exists a positive integer 1 g k S n such that Akx is the first vector (i.e., the one with the lowest exponent) which is linearly dependent on its predecessors x, Ax, Azx, . . . ,Ak"1x. k—l Definition 1.0.2 The subspace spanned by x, Ax, . . . , A x is called a cyclic sub- space of R" generated by x. Definition 1.0.3 If the monic polynomial go(/\) = X" — (amAm’1 + - - - + a2/\ + a1) 6 ’R[/\] satisfies (,0(A)x = 0 (1-2) for certain vector x E R", then the polynomial Rn/V, {X} '——> {AX}, where {x} = {y E R" : y E x(mod V)} denotes an element of ’R" / V. The map is well defined since V is an invariant subspace of A. Since dim(R”/V) = dimmn-k) = n — k, there is an isomorphism o from R" / V to R"_k. We then have a linear operator H = o o G o 0‘1 from ’R""’° to Rn‘k and the following diagram commutes: Rn/V —G—> ran/V a] ] o H Rn—k 3 7211—]: Lemma 2.2.1 w()\) is the relative minimal polynomial of R"(m0d V) with respect to A if and only if w()\) is the minimal polynomial of R"“'° with respect to H. Proof. For a polynomial (MA) 6 RM], since H = o o G 0 0‘1, we have (NH) = 0 O ¢(G) 0 0'1- Therefore, ¢(H)0({X}) = 0 4:> ¢(G)({X}) = 0: V X E R"- By G({x}) = {Ax} for any vector x E ’R", we obtain ¢(G)({X}) = {¢(A)x}: V X E 73'" (2-5) This implies ¢(G)({x}) z 0 ¢=> ¢(A)x E 0(mod V), \7’ x E 72". 8 Hence, ¢(H)o({x}) = 0 4:) ¢(A)x E 0(mod V), V x E R". (2.6) Let z/J1(A) be the relative minimal polynomial of R"(mod V) with respect to A and w2()\) be the minimal polynomial of Rn‘k with respect to H. Then we have w1(A)x E 0(mod V), V x E R", and ¢2(H)0({X}) = 0. V x E 72". It follows from (2.6) that w2(A)x E 0(mod V), V x E R”, and ¢1(H)U({X}) = O, V x E ”R". Therefore w1()\)[w2()\) and w2()\)[t,b1()\). This implies w1()\) = 1,0200 and the proof is completed. I Lemma 2.2.2 w()\) is the minimal polynomial of a vector y E R""‘ with respect to H if and only if w()\) is the relative minimal polynomial of every x E o‘1(y)(mod V) with respect to A. Proof. For a polynomial a 6 RM], ¢(H) = o o ¢(G) o 0-1 implies ¢(H)y = 0 <=> ¢(G)(0‘1(y)) = 0- Let x e a-1(y), then a-1(y) = {x}. Thus ¢(H)y = 0 <=> ¢(G)({x}) = 0- By (2.5), ¢(G)({x}) = 0 <=> {¢(A)x} = 0 <=> ¢(A)X E 0(m0d V)- 9 Hence ¢(H)y = 0 «=5 ¢(A)x a 0(mod V). (2.7) Let w1()\) be the relative minimal polynomial of x E o‘1(y)(mod V) with respect to A, and w2()\) be the minimal polynomial of y with respect to H. Then we have w1(A)x E 0(mod V), and ¢2(H)y = 0 By (2-7), $204)): E 0(m0d V), and ¢1(H)y = 0 Therefore ammo) and ammo). This implies an) = no). I Lemma 2.2.3 If V is a k-dimensional cyclic subspace generated by a regular vector x0 of A, then w()\) is the minimal polynomial of a vector y E ’Rn‘k with respect to H if and only if w()\) is the minimal polynomial of some vector z E o_1(y) with respect to A. Proof. We first prove the sufficiency. For a given x E o_1(y), there exists a vector v E V such that x = z + v. It follows that ¢(A)X = MAW + WA)“ Since w(A)z = 0, we have w(A)x = w(A)v. Note that w(A)v E V. Hence, w(A)x E 0(mod V), V x E o"1(y). This implies w()\) is the relative minimal polynomial of every x E o_1(y)(mod V) with respect to A. By Lemma 2.2.2, w(/\) is the minimal polynomial of y. 10 Next we prove the necessity. By Lemma 2.2.2 again, w()\) is the relative minimal polynomial of every x E o‘1(y) (mod V) with respect to A. In particular, for a fixed 20 E o‘1(y), w(A)zo E V. This means that there is a polynomial {(A) of degree k — 1 such that ¢(A>zo = roux... (2.8) Let ¢()\) be the minimal polynomial of A, then gb()\) is divisible by w()\). This implies that there exists a polynomial n(/\) such that $0) = MAMA)- (2.9) By (2.8) and (2.9), we have 77(A)€(A)Xo = ¢(A)ZO = 0- Since x0 is a regular vector of A, ¢()\) is the minimal polynomial of XO. Consequently, 77(A)§(A) is divisible by 9250‘). Therefore, there exists a polynomial 6(A) such that and so, {(A) : w()\)9()\). Let z : zo — 6(A)x0. Then, WAlz = WAXZO — 6(A)x0) = ¢(A)Z0 — {(A)Xo = 0- (2-10) Let i/IO‘) be the minimal polynomial of z. It follows from (2.10), yb(A)|z/J()\). Since w()\) is the relative minimal polynomial of z, we also have w(/\)|ib()\). Therefore w()\) = JO), and w()\) is the minimal polynomial of z. I Lemma 2.2.4 In Lemma 2.2.3, let m be the degree of w()\). Then x0, Axo, ..., Ak‘lxo, z, Az, ..., Am”1z are linearly independent. Proof. Suppose x0, Axo, ..., Ak‘lxo, z, Az, ..., Am‘lz are linearly dependent. Then there are constants a0, a1, . . . ,ak_1 and b0, b1, . . ., bm_1 such that 00X0 + (11AXO + ' ' ' + ak_1Ak—1x0 + boZ + blAz + ' ' ' 'l— bm_1Am_1Z = O. (2.11) 11 Let n().) = a0 + a1). + - - - + WAN-1, {(A) = —b() — blA — ° ' ' — bm_1)\m—1, then (2.11) becomes {(A)z = "(Alxo E V- By Lemma 2.2.2, w()\) is also the relative minimal polynomial of z(mod V). Thus tb(/\)|€()\). This leads to a contradiction since w()\) is of degree m and {(A) is of degree m — 1. I Before we give a new constructive proof of Theorem 2.1.1, we first define the following notations recursively. Let n1 = n and L1 = A. Suppose for integer j Z 1, we have defined the linear operator L,- in 72"!” . Let \Ilj()\) E 7?.[A] be the minimal polynomial of 72’” with respect to L,- and w,- E 72’” be a regular vector, let W, be the cyclic subspace generated by fir,- and k, = dim(Wj). If nj+1 = n,- — k,- > 0, we define Gj I an/Wj—man/Wj, (2.12) {X} *—> {Lax}- Since dim(’R"J‘ /W,-) = dim(’R"J‘+1), there is an isomorphism e,- : 7am /W,- —> 727%. (2.13) Let Lj+1 : O'j O Gj O 0.771 Ian+1 ——> an+1. (2.14) The above procedure will terminate in finite many steps, say t, that is, when n, = kt. Remark 2.2.1 From the above definition, we have n 2 k1 + k2 + . - - + kt. Lemma 2.2.5 For the above ‘Ilj()\)(j = 1,. . . ,t), we have \I’t()\)|\Ilt_1()\)| - - « [‘1'1()\). (2.15) 12 Proof. For a givenj (2 S j g t), by Lemma 2.2.1, \Ilj()\) is the relative minimal polynomial of ’Rni-1(mod Wj_1) with respect to Lj_1 and \I’j_1()\) is the minimal polynomial of ’Rni—I with respect to L,_1. Therefore, ‘II,(A)|\Il,-_1(A). I Lemma 2.2.6 For w, E R"i(j = 1,...,t) defined above, there are w,- E R"(j = 1,. . . ,t) such that W1 = $91 and n- : mu --ae<{al<{w.-}>}> - - -}>, 5': 2, . . . ,t; (2.16) moreover, (DJ-(A) is the minimal polynomial of W,- with respect to A(= L1) for j = 1, t... Proof. Since W1 E R", we simply take WI 2 v31. By definition, \II1()\) is the minimal polynomial of WI with respect to L1. For 2 g j S t, since \I'j(>\) is the minimal polynomial of w,- with respect to L,, it follows from Lemma 2.2.3 that there (1') j_1 E 03:11 (69,-) such that \II,()\) is the minimal polynomial of u(-j_)1 with is a vector u J respect to Lj_1. By Lemma 2.2.3 again, ‘I’j()\) is the minimal polynomial of a vector (i) (j) —1 uj_2 E oj_2(u j_1) with respect to L,_2. Repeating the procedure, we eventually obtain a vector w,- = u[j) E of1(ugj)) C R" such that \I'j(/\) is the minimal polynomial of w,- with respect to L1. It is easy to see from w,- = oj_1({u§j_)l}) (j = 2, . . . ,t) that W, satisfies (2.16). I Lemma 2.2.7 For w,(i : 1,... ,t) obtained in Lemma 2.2.6, let W,- be the cyclic subspace generated by W,- for A, then R":WleBW2@-~eawt. Proof. For i = 1, . . . ,t, let n, = dim(l/V,-). By Remark 2.2.1, n = n1 + - - - + nt. Therefore, we only need to show that k-l k—l k.—1 w1,Aw1,...,A1 w1,W2,Aw2,...,A2 W2,...,Wt,AWt,...,A wt 13 are linearly independent. We will show this by means of mathematical induction. Suppose W1,AW1, . . . , Ak1_lW1, W2, AW2, . . . ,Ak2—1W2, . . . ,Wj, AWJ', . . . , Akj—IWJ' are linearly independent, we will show that W1,AW1,...,Ak1-1W1,W2,AW2,...,Ak2—1W2,...,Wj+1,AWJ-+1,...,Akj+1_1WJ-+1 are linearly independent. In fact, if W1,AW1, . . . , Ak1’1W1,w2,AW2, . . . ,Ak2—1w2, . . . ,Wj+1,AWJ-+1, . . . ,AkHI—IWjH are linearly dependent, then there are polynomials {1(x\), §2(/\), . . ., §,+1(A) of degrees k1 — 1, k2 — 1, . . ., kj+1 — 1, respectively, such that §j+1(A)Wj+1 = {1(Alwl + 52(A)W2 + ' ° ' + €j(A)Wj- (2-17) By the definitions of Lj, G,- and 0,, we have for any vector x E 72’” 05({ij}) = 05(G:({x})) = 0: 0 G: 0 05—1 0 05({x}) (2.18) = Lj+10j({X}), j: 1,. . . ,t — 1, and oj({vT/',-}) = 0, j: 1,... ,t — 1. (2.19) By L1 = A and (2.17), we have {€j+1(L1)WJ-+1} = {61(L1)W1} + - ° - + {5,-(L1)w,-} = {€2(L1)W2} + ° - - + {€5(L1)Wj}9 Therefore, 01({§,+1(L1)w,-+1}) E 01({€2(L1)W2}) + ' ° ' + 01({§,-(L1)w,}). 14 It follows from (2.18) and 01({w2}) 2 W2 (refer to (2.16)) that §j+1(L2)01({Wj+1}) = {2(L2)01({W2})+'°'+€j(L2)01({W2‘}) = 52(L2)V~V2 + ' ' ' + §j(L2)01({Wj})- So, {§j+1(L2)01({Wj+1})} E {53(L2)01({W3})} + ' ' ' + {€j(L2)01({Wj})}- Furthermore, 02({€j+1(L2)01({Wj+1})}) = 02({€3(L2)01({W3})}) + ° - - + 02({€j(L2)01({Wj})})- From (2.18) and 02(ol({w3}) = a, (refer to (2.16)), we obtain €j+1(L3)02({01({WJ-+1})}) = 63(L3WV3 + ' - ' + €j(L3)02({01({W5-})})- Repeating the same procedure, eventually we have than-aw... = 0. Obviously, the above procedure is equivalent to acting on both sides of (2.17) by 05({- ' - 02({01({'})}) - - -})- Since \Ilj+1()\) is the minimal polynomial of v~vj+1 with respect to LjH, we have \II,+,(,\)|§,-+1(,\). But the degree of €j+1()\) is less than the degree of \II,+1()\). Therefore we have {24100 = 0- This together with (2.17) contradicts the assumption of the induction. Therefore, W1,AW1,.. .,Ak1_1w1,w2,Aw2,.. . ,Ak2'1w2,. .. ,wj+1,ij+1,...,A'°J'+1—1w,-+1 are linearly independent. I 15 Remark 2.2.2 Note that the relations between ‘36,, u§j+1), . . . , ugt) in Rn)" for L,- are the same as the relations between W1, W2, . . . ,wt in ’R," for A in the proof of Lemma 2.2.6. Therefore, forj : 2,... ,t, ~ ~ k,-—1~ (1+1) 0+1) kj+l-1 (1+1) (0 (t) Wj,Ljo,...,Lj Wj,uj ,Ljuj ,oo-,Lj uj ,... u Lju ... are also linearly independent. Lemma 2.2.8 For Wj(j : 1,. . . ,t) defined in Lemma 2.2. 7, \I!,-+1(/\) is the relative minimal polynomial of R"(mod (W1 63 W2 69 - - . EB Wj)) forj = 1,. . . ,t - 1. Proof. Ifj : 1, then W1 = W1 and W1 2 W1. By Lemma 2.2.1, \II2()\) is the relative minimal polynomial of R”(mod W1) with respect to L1. Now for a given 3' (2 g j g t — 1) and any vector x E R”, there are, by Lemma 2.2.7, polynomials §1(A), (20‘), . . ., {t()\) with degree k1 — 1, k2 — 1, . . ., kt — 1, respectively, such that ‘I’j+1(A)x E {1(A)W1 + 52(A)W2 + ' ' ' + {,(A)wt. (2-20) Let y = 05({' ' - 02({01({x})}) ° ' -}), then y 6 RM“. Recall that \Ilj+1()\) is the minimal polynomial of 72"?“ for Lj+1~ Acting on both sides of (2.20) by oj({- - '02({01({'})}) - - }) and utilizing (2.18), we have 0 E \I'j+1(L,-+1)y E Ujli'"02({01({‘I’j+1(A)x})})"'}) ~ ' 2 E €j+1(Lj+1)Wj+1 + €j+2le+1)u§J++1 ) + ' ” + {Khaki-531- It follows immediately from Remark 2.2.2, £j+1(/\) = = €50) = 0- 16 Hence (2.20) implies that for all x E ”R." ‘I’j+1(A)X E 0(mod (W1 69 W2 63 ' ' ° 63 Wj». (2.21) Let {JR-HQ) be the relative minimal polynomial of 7?."(mod (W1 69 W2 EB - - - EB Wj)) with respect to A, then Tj+1(A)|\II,-+1(A). By Lemma 2.2.6, \I'j+1()\) is the minimal polynomial of {Cg-+1 with respect to Lj+1 and ‘i’le+1)V~Vj+1 E 01H“"02({01({‘T’(A)Wj+1})})°°°}) E 0- Thus \II,+1(/\)|\II,-+1(A). Consequently, \II,+1()\) = {IQ-HO). Therefore \Ilj+1()\) is the relative minimal polynomial of ’R"(mod (W1 69 W2 EB ~ - - 69 W,)) with respect to A. I Proof of Theorem 2.1.1. Combining Lemma 2.2.5, Lemma 2.2.6, Lemma 2.2.7 and Lemma 2.2.8, the proof of Theorem 2.1.1 follows immediately. I Remark 2.2.3 From the results obtained in this section, it is not necessary to modify each regular vector w,- of L,- in R": to a regular vector w,- in ’R,"( j = 2,. . . ,t) if we only need to determine the minimal polynomials \I'1()\), \Ilg()\), ..., \Ilt()\) defined in Theorem 2.1.1. 2.3 A constructive proof of Theorem 2.1.2 In this section, we give a constructive proof of Theorem 2.1.2. Lemma 2.3.1 [14] For a given \II(/\) E R[/\], via finitely many purely rational oper- ations one can always find a squarefree relatively prime basis {901, (p2, . . . , 903} C RH] of \Il()\) such that ‘1’(/\) = [cpi(/\)lm‘[902(/\)lm2 ' ' ' l¢s(>\)lm‘, (93-22) where m,- > 0(i = 1,2,...,s) are integers and m1 < m2 < < m,. 17 Proof. Group all simple factors of \II()\) in the complex number field C which have the same multiple, then \I’()\) can be represented by ‘1’”) E @10le1 [90200],”2 ' ' ' l90s(/\)lm’i (2.23) with m,- > 0 for i = 1,2,...,s and m1 < m2 < < m,, where gcd(<,o,-,%<,o,~) = 1 and gcd(\)”-905(/\)5 77ml+1()‘) : dm1(’\)/dm1+1()‘) : 902(A) ' ' ° (FAA), < 3 (2.26) 77mz()‘) : dm2—1()‘)/dm2()‘) : 902(A) ' ° ° (FAA), ”uni-10‘) : dm2(’\)/dm2+1()‘) = 903(A)""Ps()‘)1 77m.()‘) : dme—1()‘)/dma()‘) : 904A). 18 Hence, 901(/\) : ”milAVUsz‘): $020) = nmz()‘)/7lm3()‘)a < s (2.27) SOs—1(A) : ”ma—1(A)/nma(A)7 (Ps()‘) : 77mg()‘) Note that (2.24), (2.26) and (2.27) can all be accomplished by finitely many ra- tional operations over ’R[)\], hence 901, (p2, .. ., (0, must belong to 7?.[A]. I Lemma 2.3.2 [14] For given ‘I’(/\) and @(A) in 72D] with (A)|\II()\), one can always find a common squarefree relatively prime basis {901,902, . . .,go,} via finitely many purely rational operations such that ‘I’W E [¢1(/\)lm1[902(A)lm2 ”1905091“: (2-28) and @(A) E [<.01(/\)l"1[WOW2 ° ' ' [905(A)l”’ (2-29) withn, 3 mi, fori = 1,...,s. Proof. By Lemma 2.3.1, we may assume \II(A) has a squarefree relatively prime basis {171, 172, . . . ,nk} with decomposition ‘I’(/\) E [77100]“ [71209]“ ”17715091“. and ()\) can be expressed as 4’0) E [771(/\)]q1[772()~)lq2 °"[77k(/\)lq’°91(>\)a where each q, s p,- is a nonnegative integer and OI is not divisible by 77,- for i = 1,2,...,k. 19 If 7',- : gcd(17,-, OI) has positive degree, then 77,- has a nontrivial decomposition 775' E Tigi- (2'30) For those i’s with 7',- = 1, let 5,- : 77,-. Then Ol()\) can be rewritten as ('91 = 797'? ---T,:’°Og, where each l,- (i = 1,2,. . . ,k) is a nonnegative integer. Since ()\)|\II()\), not all the 73’s can have zero degree, and hence, deg 9;; < deg 91. Replacing 7', and g, for 77,, a new squarefree relatively prime basis of \II contained in {71,61, 7'2, £2, . . . ,Tk, {k} is obtained such that w = 551315521552 - 9 - are, and _ ql +11 91 ‘12 +12 92 Q]: +11: qk q)_Tl 172 2 "'71: £11: 92' Repeating the preceding procedure, a common squarefree relatively prime basis {901, . . ., 90,} of \Il()\) can be obtained, via finitely many steps, for which ‘I’(/\) E [(A)|\II(/\) implies O(A)|\II()\), there exists at least one (0,-(1 g i g s) such that nglioiiG) Té 1 I Remark 2.3.1 Even if @(A) is not divisible by (NA) in Lemma 2.3.2, a common squarefree relatively prime basis of both ‘II()\) and ()\) can still be obtained via finitely many purely rational operations. In fact, if OM) in (2.32) is not a constant, then, by Lemma 2.3.1, O(A) has a squarefree relatively prime basis {71, 7'2, . . . ,n}. Since 20 gcd(go,-(/\),O) = 1 fori = 1,2,...,s, {cp1,<,02,...,<,0,,7'1,7'2,...,Tl} is a squarefree relatively prime basis of (I). Obviously, it is also a squarefree relatively prime basis of \I'. Proof of Theorem 2.1.2[14]. The assertion of the theorem is a direct conse- quence of the above lemma by mathematical induction. Suppose the assertion of the theorem holds for 2 S j < t, that is, {IQ-Hz, has a common squarefree relatively prime basis {771, 772, . . . ,173}. Since {n1,172, . . . ,7],} is a squarefree relatively prime basis of ‘1’,- and \I/j+1|\Il,-, a new squarefree relatively prime basis of II,- and ‘11,“ can be obtained by repeating the same procedure described in Lemma 2.3.2. Moreover, each n,- is either in the new squarefree relatively prime basis or can be rewritten as a product of polynomials in the new squarefree relatively prime basis. It follows that the new squarefree relatively prime basis {901, 902, . . . ,go,} of I, . o t 0 .—'1 and ‘11,“ 1S also a common squarefree relatlvely prlme baSIS of {\II, 3:1. I Remark 2.3.2 Squarefree decompositions of the minimal polynomials \Ill,\112,. . . , It in Theorem 2.1.2 can be implemented by using symbolic computation. From Theorem 2.1.3, as long as those decompositions in (2.3) are achieved, the structure of the Jordan normal form of A can be determined exactly, without knowing the zeros of (ok’s. The approximation of those zeros of (pk, k = 1, . . . , r, has no influence on this structure. 21 CHAPTER 3 THE ALGORITHM Following what was described in the last chapter, we outline our algorithm for deter- mining the Jordan normal form of an n X n matrix A in ’R". 0 Step A. Find the minimal polynomials \Ill()\),\Ilg()\), . . .,\Ilt(/\) in Theorem 2.1.1. A0. j=1,n1=n,A1=A. A1. Find a regular vector w,- in 72’” and its minimal polynomial \Ilj(/\) with respect to A,. Denote the degree of \Ilj()\) by 19,-. A2. Let n,“ = n,- — k,. If nj+1 = 0, then Step A is finished. A3. Extend wj,A,-w,-,...,Aj’ w,- to a bas13 w,,Ajw,-,...,Aj’ wj, vkj+1, n' . vkj+2, ..., v", of ’R 1. Then {vkjfl}, {ij+2}, --~,{v,,j} form a baSlS of the quotient space R’” [VT/3, where W,- is the cyclic subspace generated by «7,. With this basis, we may define the isomorphism o,- as follows: 03' I an/Wj ___> {Ian-“'1, T yF—) (C1,...,an+1) , where y E levkj+ll + szvkj+2} + ' ' ' + an+1{Vk,-+n,-+1}- 22 A4. A5. Compute the matrix A,“ of the linear operator Lj+1 = o,- o G,- 0 0,71: 72“?“ ——> 72"?“ under the standard basis {e[j+1),e§j+1),...,elf,:i)} of 72W“. Letj=j+1andgotoA1. 0 Step B. Determine the structure of the Jordan blocks in the Jordan normal form of A. B1. B2. Find a common squarefree relatively prime basis {901, 902, ..., (0,} of \III, \Ilg, ..., \II,_1, ‘1', in ’R[A] such that r e, = prison - ”901"”, ‘1’2 = Spinal???” . . .9031”, t ...... (3.1) ‘I’t—i E 90:1"1’19031E1'2 . . . soft—1w, \ \I’t = 9011725490?” ' ”90:71”; where mu. 2 mm 2 2 mth 2 0 for k : 1,2,”Wr. To achieve this, first use Lemma 2.3.1 to obtain a squarefree relatively prime basis of ‘15; secondly, modify the squarefree relatively prime basis of ‘111 so that it becomes a common squarefree relatively prime basis for both ‘111 and ‘15; then apply Lemma 2.3.2 inductively to modify the known squarefree relatively prime basis of \I11,\Ilz,~-- ,\I'k_1 to a new squarefree relatively prime basis of \III, ‘112, - - - , \Ilk until k = t. From Theorem 2.1.3, corresponding to each 902”” of (3.1) for k = 1,. . . ,r and i = 1,... ,t; there are pk = deg 90,, Jordan blocks of order mu, of the form 23 (Ah, 1 0) Ala! , l=1,...,pk, (3.2) 1 10 1 mi,kxmi,k in the Jordan normal form of A, where the A, ,’s are zeros of (pk. Notice that the structure of the Jordan normal form of A has been determined without knowing those Auk. 0 Step C. Calculate the Jordan normal form of A. After the structure of the Jordan normal form of A has been determined in Step B, the remaining task is to approximate all zeros of the polynomials (pl, 902, - - - , (,0, obtained in Step B1 to fill in the diagonal entries of the matrices in (3.2). 24 CHAPTER 4 IMPLEMENTATION OF THE ALGORITHM In this chapter, we will discuss in detail the implementation of the algorithm outlined in Chapter 3, and present an example to demonstrate the implementation of the algorithm. 4.1 Implementation of Step A Let n1 2 n and A1 = A. For j = 1,2,..., choose a vector w E 72’” at random. Consider the n,- X (2n,- + 1) matrix H, : (w,A,-w,...,A;."w, 1,”), (4.1) where In, 2 (.29), - - - , elf?) is the n,- X n,- identity matrix in 72’” . In general, we denote the k X 16 identity matrix in R,“ by Ik. Suppose r,- is the largest integer for which w, Ajw, . . . , Agrlw are linearly inde— pendent. In general, Gaussian elimination with row pivoting can be used to reduce H ,- to the following form 25 / h1,1 ° ' ‘ h1,r_,- h1,r,-+1 * * \ O hr-r' th,7‘j+l * * * H; 2 ... I... (4.2) J 0 0 . . . 0 K 0 0 . . . 0 j W ‘Nf—I where hm E 0 for i = 1,. . . ,r,, and 1;], is the n,- x n,- matrix resulting from applying the above Gaussian elimination process on Inj. Notice that if computation is symbolic, then pivoting can be replaced by finding a nonzero element in Gaussian elimination process. In Appendix, an efficient algorithm of reducing H, to H 3'" is developed to reduce the memory storage and drastically speed up the computation. Solving the linear system (hm h1,2 hm, \ K at,” \) {h1,rj+l \ h2,2 (12¢,- 0(3) = h2,r,-+1 (4'3) \ O hrw‘i A \ “ii-1 A ) hurt“ ) yields the minimal polynomial of W \I/(A) : A” — (agllA’V-l + . - - + a[j)/\ + a8”). (4.4) Let nj+1 = n, — r,. If nj+1 = 0, then w,- : w, \I'j()\) 2 \Il()\) is the minimal polynomial of 72’” with h,- = r,, and an = Wj Z Span{vir,-, Ajo, . . . ,Arj—IWJ'}. .7 26 Step A is finished. Now let us assume nj+1 > 0. For the n,- x 77.,- matrix 1;], in (4.2), write 1;], = (ugj), - - - , ugh. It is easy to see that there are nj+1 columns of 1;), say ug), . . . ,ufiifl, for which ug) = egfil,u,(~:) : eglzvumfiifl = cg). For instance, when H; is J achieved by applying Gaussian elimination without requiring row pivoting to H ,- then Iiij : (11(1)). ' ' ,ug), Biz-L1) ' ' ' ’62)), namely, H _ H _ H _ ui: _ Uri-+1) — erg-+13) p — 13 ' ' ' inj-i-l' The same consequence takes place if the pivoting only occur within the first r,- rows. If, during the elimination process, say only the pth row (p S r,) and the (r,- + k)th row have been interchanged, then * . . - 1 - ' ' . ' . = (.5, . . . .ui’3.,e£j1..,u§..2n . . . 529:2. - - - ,es.:.z..-.u:::ae::ha - - - >9 that is, uli) __ uli) _ 8(1) for _ 1 . t — k ip _ T'j‘l'P— Tj+p p— ,...,n,+1 excep p— ’ while ufi) = rug) 2 eng. When Gaussian elimination is applied to transform H ,- to H J" , in essence, a non- singular matrix P,- is constructed making PJ-Hj 2 H31“. Accordingly, \ Pj(W,AjW,---,A1:j—1W)= = {hm hl srj where O hum and PJ' : Film“ 2 Pj(eij)a ‘ ' ' reg!) : (11(3),. ' ' with?) = 11:,- (4'5) with Pjezlj) = uf-j), i = 1,... ,n,. In particular, Pee) : um (i) .1 i1 i1 E erj+1a Pjefii’ E Uii) E 6.3.192: (4.6) It follows that w, Ajw, . . . ,A;j_lw and elf), . . . , e]?+1 are linearly independent, since .7 P,- is nonsingular and (hm h,,., \ P, eff) P,- egg,“ B< n,-+1 matrix. Since {Wj, ---, Ak’ _lw em e0) } is a basis of RT” for - 1 n- there are constants ' j) i inj+1 7 p_ 3'”) 3+1) .7 6,1, b1“), (32,19, . . ., bnj,P SUCh that ~ k,-—1 ~ ' ' Ajei: E ’31.ij + - - - + 5155.ij “'2' + bk5+1.peif) + ' ' ' + bn,,,,eff,)+l ( b1,p \ ,jA k’_1~ (j) 6(3) ) bki’p (4'14) Wj’ 61,1 1 )6 inj+1 bkj+1,p K bum / Z Tj(b1.pv ' ' ' abnj.p)T- 3O This implies {Ajeip} Z bkj+1,p{eg)} + ' ' ' T bnj,p{ei.{:+,}v P : 11° -- ,n,-+1. (4'15) Since 0, is defined by 0',‘ I an/Wj—anj+la y l——> (C1) ° ' ' acnj+1)T7 Where y = Cl{eg)} + C2{ez(.‘:)} + . . - -[— an+1{€,(-'::+1 }, we have o,({e,(-:)})=e,(,j+1), p: 1,...,n,-+1. For p = 1, . . . ,n,-+1, it follows from the definition of L,-+1 that . (1+1) _ . . _1 (3+1) L,+1ep — o, o G, o o, (e p ) E 0i 0 Gj({€ip}) E 0j({Aj€z‘p})- (4-16) By (4.15) and (4.16), (bkj+1,p, - - - ,bn,,p)T is the p—th column of A,-+1 forp : 1,... ,n,-+1. That is, K bkj+1,1 ' . ' bk,+1,n,-+1 \ Aj+1 E (4.17) K bnivl bniini+1 A In the following, we deduce a simple representation of A,+1. By (4.14), we have { blil blvnj+1 \ . 5+1) —_— bkj+1,1 biog-+1.11)“ . (4.18) \ bn,,1 ' ' ' bn,,n,-+1 / It follows from (4.12) and (4.13), 31 _____...-__ q .v—r‘u 5-- -_‘__ '. 'S‘l'flfz Aj+1 Z ( O ,Inm)T,.‘1(A,-e,-,,-~,A,e,,,j+1) M2 M4 B2 (4.19) B = (M2, M4) = M231+M432. B2 From (4.8), both matrices M2 and M4 are available after Gaussian elimination is applied in the first step. Therefore, A,-+1 can be obtained with minimal computation effort. In particular, if no pivoting is employed in the elimination process, M4 = In, +1 and Aj+1 = M2B1 + Bg. Remark 4.1.1 It follows from G,({e,-p}) = {A,-e,-p} and {4.15) that Gj({eip}) Z bkj+19pleg)} + ' ' ' + I’m-pkg“ }a P Z 1: ' ' - ,n,-+1. Thus, A,+1 is also the matrix of G,- under the basis {e,-,}, . . . , {ein,+1} of Rni/Wj. 4.2 Implementation of Step B Now we discuss the detail of determining the structure of the Jordan blocks in the Jordan normal form of A. A common squarefree relatively prime basis {$01, (02, - - -, (0,} in R[/\] for the min- imal polynomials I1, ‘112, . . . , \I’t constructed in Step A satisfying the decomposition, ‘1’: = salm’lsoém'z - - - 90?”. ‘112 Z (plnzgsognm ' ° ' 90in”: (4.20) in = wimp?” ' - - 90W, 32 where m1), 2 m2), 2 - -- 2 mth 2 O for k = 1, 2, . . . ,r, can be generated by means of symbolic computation as follows. For k1 2 deg ‘I'1(/\), let d0()\) = \I'1(/\) and dl (110‘) E 8Cd(‘1’1()\)ig)7‘1’1(>\)> 7 l E 1)"')k1° In the proof of Lemma 2.3.1, it was shown that if {61, . . . ,9,} is a squarefree relatively prime basis of \Ill()\) such that 210‘) E [01(A)lm1l92(/\)lm2 ° ' ' [9300]” with m, > 0 (i = 1,...,s) and m1 < m2 < < ms, then for 77,-(A) = d,_1(/\)/d,-()\), i = 1,. . . ,m, — 1, we have, ignoring the constant factors, l ab) = dam/din) = alone) - - no), 77m10‘) Z dm1—1()‘)/dm1(’\) Z 61006200 ' ' ' 950‘), 77ml+1()\) : dm1(’\)/dm1+1(/\) : 02(A) ' ' ° 63(A)) 7751520) E dm2-1()‘)/dmz()‘):02()‘)°”63()‘)7 77m2+1()\) Z dm2()‘)/dm2+l()‘) : 63(A) ° ' ' 68(’\)1 7755.1%) E dm.—1(>\)/dm.(/\)=9s(>\) 77m.+1(>\) E dm.()\)/dm.+1()\) E L This implies {-77—1,@,~-,—nm' }={1,...,1,91()K),1,...,1,92()K),---, 1,...,1,63()K)}. m1”- m2_m1_ ma—ma-l— 33 Therefore, to generate {61, . . . ,6,} and find m1, m2, . . . ,m,, we may compute 171 E E 772’ 03’ 774, yielding the sequence 1,...,1, ""1“ 751, 1,...,1, W¢L 1, RH 7flni+2 WW nn1+n2+3 n1 112 and m1 27214-1, 61(A) 2 EL, ”m1+1 m2 =m1+n2+1, 02(/\) 2 EL, nmz+1 "m- m,- =m,-_1 +le+1, 9,-(A) = 7—2;, "‘1' The process stops when m1 deg 61 + + m, deg 6, 2 k1 for certain s, and clearly, {61, . . . ,03} constitutes a squarefree relatively prime basis of \Ill().) for which rho) = 161r11621m2---16.(A>re. (4.21) By the polynomial “long division” in R[)\][15], \I/2()\) can be decomposed as ‘I’2(/\) E [91()\)lb1l92()\)lb2 ° ° ° [93(A)lb‘91()\)5 (4-22) where 0 S b, S m, and Ol()\) is not divisible by any 6,, for i = 1,... ,3. With this decomposition, we calculate T,()\) = gcd(6,, 91), i = 1, . . . , s by symbolic computation [18]. For those i’s with deg T,()\) > 0, 6, has a non-trivial decomposition 0:0) E Til/96:00: 2' E 1. - - ~ , s. For T,()\) = 1, we let {,(A) = 0,()K). Writing 91W E #1007520) ' - - T§‘(>\)@2(>K): 34 and replacing 6,()\) in (4.21) and (4.22) by r,()\)§,()\), for i = 1,... ,5; yields ‘1’1 E Tinléin172m§2n2°”7'?‘§;n‘i (4.23) \112 : Tf1+11€i>1T2bz+lz€gz . . . T:,+l,€f,®2()\). Obviously, ignoring those T,()\)’s with degree 0 in {7-1, {1, . . . ,7',, {3}, a new square- free relatively prime basis of \Ill()K) is constructed. By repeating the same procedure for (4.23) until the factor O,()\) of ‘1/2 which is not divisible by polynomials in the basis becomes a constant, a common squarefree relatively prime basis of both \Ill()\) and I20.) is produced. Inductively, when a common squarefree relatively prime basis of \II1()\), ..., \Il,_1(/\) has been constructed, a common squarefree relatively prime basis of both \II,_1()\) and \II,(/\) can be constructed by the above procedure to serve as a common squarefree relatively prime basis of 4110.), . . ., \II,()\), and the implementation of the part B] is accomplished. Now we consider the implementation of the part B2. Corresponding to each (pini’k in (4.20) for k = 1,... ,r and i = 1,. .. ,t; there are pk '2 deg go], Jordan blocks of order mm, of the form {AM 1 O) )‘k,l , l=1,...,pk, (4.24) 1 (0 up.) mi,kxmi,k in the Jordan normal form of A, where the Aw’s are zeros of (pk. Without knowing those zeros, the structure of the Jordan normal form of A is determined at this stage. Remark 4.2.1 The computation in Step A involves only Gaussian eliminations and divisions between polynomials, and the major computation in Step B is to generate a 35 common squarefree relatively prime basis. Those computations can all be carried out by symbolic computation. Therefore, if all the entries in the given matrix A consist of only rational numbers, then the structures of the Jordan blocks in the Jordan normal form of A can be determined exactly. 4.3 Implementation of Step C In order to fill in the diagonal entries of the matrix in (4.24), the zeros of 901, . . . ,go, in (4.20), which may not be in ’R, need to be evaluated. The approximation of those zeros can be computed by well established subroutines, such as the Jenkins-Traub routine [2] in the IMSL library. When those zeros are available, the complete Jordan normal form of A is achieved. Although those AkJ’S are eigenvalues of A, it is inappropriate to approximate those numbers by calculating the eigenvalues of A using standard codes such as EISPACK [16], or more advanced LAPACK [1], because if A“ is obtained as an eigenvalue of A, it might fail to provide the information of the location to which ’\k.l belongs. 4.4 Computational result Golub and Van Loan presented a 10 x 10 matrix [7] as shown below to illustrate the difficulty of calculating the Jordan normal form of a matrix by numerical computation: ( 1 1 1 —2 1 —1 2 —2 4 —3 ) —1 2 3 —4 2 —2 4 —4 8 —6 —1 0 5 —5 3 —3 6 —6 12 —-9 —1 0 3 —4 4 —4 8 —8 16 —12 A: —1 0 3 —6 5 —4 10 —1o 20 —15 —1 o 3 —6 2 —2 12 —12 24 —1s —1 0 3 —6 2 -5 15 —13 28 —21 —1 0 3 —6 2 —5 12 —11 32 —24 —1 o 3 —6 2 —5 12 —14 37 —26 ( —1 0 3 —6 2 —5 12 —14 36 —25 ) Here, we use the algorithm described in the above sections to obtain the Jordan normal form of A. For symbolic computation used in Step A and Step B, we take advantage of routines available in Mathematica[18]. The main procedure and result are described as follows. Let A1 = A and 110 : (e1,e2, - - - ,610) be the identity matrix of order 10. Ran— domly choose a non-zero vector w 6 R10, for example, we choose 1 l l l l 1 l A) 3’4’5’6’7’8’9’10' 1 1949 2537 2137 _ 233 21339 _ 20792 7597 _ 79013 _ 2275377 _ 1322597 \ 1230 1230 1230 123 1230 315 33 123 1230 3' 1'5 1 209 _ 901 _ 1279 _ 977 _ 30 29 _ 45334 _ 7357 _ 72293 _ 2033957 _ 2332;114 2 330 1230 315 33 315 13 33 330 315 1 _2_ _ 563 _ 2593 _ 293 _ 23729 _ 3434 _ 3317 _ 32335 _ 1300997 _ 434339 3' 35 4 o _420 _1 _420 —35 —12 42 —420 _3—5 1 _ 43 _ 2357 _ 4301 _ 1573 __ 23739 _ 39303 _ 5701 _ 115013 _ 1371077 _ 4911133 4 315 1230 330 33 315 315 9 33 315 315 1 _ 37 _ 2773 _ 10337 _ 17333 _ 107101 _ 15752 _ 131173 _ 1340957 _ 7395737 _ 5371193 5 2—52 1230 1230 "We _1230 03 "—130 _330 —i'2_60_ “'3_15 1 _ 37 _ 1g _ 992 _ 3233 _ 19729 _ 9734 _ 3219 _ 34947 _ 1513039 _ 2275274 3 105 23 105 105 210 35 10 35 210 105 1 _ 1019 __ 321 _ 24339 _ 10121 __ 35341 _ 73357 _ 2173921 _ 1312159 _ 3347211 _ 23933991 7 2529 "2520" "3—15 330 _252 To Tea W Too 1 _ 221 _ 3791 _ 51 13 _ 33129 _ 254331 _ 759197 _ 321407 _ 3332373 _ 19934037 _ 30255353 3 504' 2520 504 2' 520 2520' 2520 330 2520 "2' 5' 2'0' ' 2520 1 _ 577 _ 1147 _ 4297 _ 3977 __ 3537 _ 31347 _ 53933 _ 530719 _ 3353341 __ 2527373 9 1' 2' 3' ‘0 ‘ 4T 0 To W1 Ta —1 o 5 _0 o T271 .176— T1 ( 1 _ 197 _ 91 _ 2531 _ 1043 _ 123039 __ 191039 _ 131301 _ 341032 10030937 _ 15137233 ) 10 42' 0' 2:52 252 315 ' 12' 3'0 330 " 1'30- ""‘31‘5' 1230 330 Utilizing Gaussian elimination with row pivoting for (K [w],110), we obtain (K* [w], {0) where ( 1 1949 2537 2137 _ 233 _ 21339 _ 20792 _ 7597 _ 79013 _ 2275377 _ 1322597 1230 1260' '1230 173 1230 315 33 123 —1"23' 0' "—315 0 _ 97 _ 421 __ 4133 _ 405 _ 34009 _ 11353 _ 7277 _ 23351 _ 1953317 _ 431719 We 2'3— 70‘3 '2'3' Eli—0‘4 ‘1—05 ‘27]— ‘23— '_BTO'— T _ 337 _ 7405 _ 5713 _ 135333 _ 123923 _ 1722455 _ 322241 _ 1 12334 _5119131 0 0 10473 1423 5233' 31423 "—7357 31423 "‘1‘74T ' 31133" E 2019 0 0 0 _ 47 _ _ 2309 _ 1702 __ 57799 _ 37523 _ 323331 _ 242519 23 2 —4440 76—7442 T303 ‘3‘??? Z'T‘o4 WIT 2265 ' K * [W] — § 1 — _ 331 _ 332 _ 4035 __ 1.3519 _ 27921 _ 19253 _ 33351 0 0 0 0 19740 4935 3943 3290 1313 235 1974 0 0 0 0 0 _ 5 _ 35 _ 250 _ 495 _ 5025 _ 22945 193' 3' 93 993 331' 332 362 0 0 0 0 0 0 0 0 0 0 0 0 O 0 0 0 O 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 K 0 0 0 0 o 0 0 o o 0 0 37 and / 1 0 0 000 0 000K —§ 1 0 000 0 000 3% gggg 1 000 0 000 {+091 333; —% 100 0 000 f0: —% % —% fgfioo 1 000 (4.25) 17:31 _20743 4371 M 0 0 _flg 1 0 0 6620 6620 3310 1655 662 —§ 3 0 —§10 0 000 —§ g 3245—201 0 000 ——-———oo :2—310 \————oo 3—23—01) It follows from 4.25 that W1,A1W1, . . . , A5w1, 65,86,639, 310 is a basis of R10. 1 Solving the linear system 1 1949 2537 2137 _133 _21339 K K a K { _20792 K 1230 1230 1260 123 1230 0 315 0 __ fl _4_21 _4133 _flf; _34009 a _ 11353 230 230 340 23 340 1 105 0 O _ 337 _ 7405 _ 5713 _ 135333 a _ 123923 10473 31423 5233 31423 2 _ 7357 _ 47 _ 77 _ 2309 _1702 0 0 0 23424 4404 23424 33 3303 331 332 4035 0 0 0 0 — 19740 4935 34 3943 5 35 K 0 0 0 0 0 1933 j K 35 j K 1933 ) yields the minimal polynomial of w ‘II()\) = A6 — (a5A5 + a4A4 + a3A3 + (1.2)\2 + 31A + 30) = A6 — 13A5 + 69,\4 — 19113 + 29012 — 228A + 72. Since \II(A1)e,- = 0 for 7' = 5,6, 9, 10, \II()\) is the minimal polynomial of R10 and W is a regular vector. Therefore, k1 = deg\II = 6, LLH :: {gala/qleala'. '111i661}1 moreover, {e5}, {36}, {39}, {310} form a basis of the quotient space RIO/W1. From (4.8) and (4.25), we have 1 4 _3 —5 5 0 5 1 0 0 0 0 0 —§ g g —% 0 1 0 0 0 0 (M2,M4)= 39 392 152 632 28 48 m —m —m 375 0 0 fi ——5 1 0 164 539 78 844 42 248 575 _fi ‘55 m 0 0 75 _W 0 1 According to (4.13), we have 1 —1 4 -3 2 —2 8 -6 3 —3 12 —9 4 —4 16 —12 —4 2 —1 B1 5 0 5 (A185, A163,A189,A1€10) = = 2 —2 24 —18 32 2 —5 28 —21 2 -5 32 —24 2 —5 37 —26 2 —5 36 —25 Therefore, 142 I: (1772,1VT4) 1: A45 131 d—.A4; 132 Z: 32 1 —1 4 —3 —§ 3 0 —§ 1 0 2 —2 8 —6 2 7 2 12 —5 5 5 ‘5 0 1 3 —3 12 —9 + fl_fl _gg $300 4—416—12 125 375 75 375 5 —4 20 —15 1% _& _fl fl 0 0 375 375 125 375 39 (0 000\{2—528—21\ 0 000 2_532_24 _4400 28 48 z: 199 76 fi—fiio 2—537—23 W—fill 42 243 2 —5 36 —25) 33 39 \3—701)K \2—5—53—45} Let 14 be the identity matrix in 724, we choose a non-zero vector 3 4 ,Z,g)eR4. COIN W=(%, Using Gaussian elimination with row pivoting on (K [w], I4 ) where {1 ——;- —2 —10 —32K 2 3 4 3 —2 —10 —32 —88 K[w] = (W,A2W, A2w, A2w, A2W) = , 1 212 fl 2_02_9. 35_7 3 10 50 100 25 K1 13 2241 2.4% -1544) 4 300 300 300 300 yields the matrix (1% —2 _10 —322 1 0 00K 0 —§ —9 _27 —72 —% 1 0 0 (K*[W]JZ)= 139 2141 927 I 1219 769 0 0 ’fi _W _2—5 - “W W5 *1 0 507 2535 I 2743 829 751 K0 0 0 —m W - W3 71—7 *2—75 1) Solving {1% —2 —10K/aoK {~32K 0 —§ —9 ~27 31 ~72 : 7 139 2141 927 0 0 _W ‘W “2 *7 \o o o -%)\an —) we obtain the minimal polynomial \I'()\) of w in R4 \I’(/\) = A4 — (33A3 + 072)? + 01A + 0,0) = /\4 — 10%3 + 37/\2 — 60A + 36. 40 Since deg \I'O‘) = 4 = 10 — deg \I/1(/\), we have \I/2()\) = \II()\). For ‘I/1(/\) and @201), a common squarefree relatively prime basis { (A — 1), (A —— 2), (A —— 3) } in ’R[)\] is found (by means of symbolic computation) for which 91(3) = (3—2)3(A—3)2(A-1), \I'2(/\) = (3-2)2(A—3)2- It follows that the Jordan normal form J of A is of the form (2100000000) 0210000000 0020000000 0003100000 J_ 0000300000 — 0000010000 0000002100 0000000200 0000000031 K0000000003) 41 APPENDIX 42 APPENDIX: An elimination algorithm 1 An elimination algorithm For a nonzero vector w E R" and an n X 71 matrix A, we consider the following n X (271 + 1) matrix H = (w,Aw,...,A"W, In). In this appendix, we develop an efficient elimination algorithm to reduce H to the following form O hrr hr,r+1 * * It 11*: ' I n 0 0 0 \ 0 0 0 / V—v—-’ \_\,_/ where 7‘ > 0 and hm 75 0 for 7' = 1, . . . ,7. Recall that we are only interested in obtaining the first 7‘ + 1 columns and the I; part in the above matrix. Our algorithm is therefore designed towards this main goal. Let W = (101, - - - ,wn) E ’R" with 711,74 0 (1 S 2' g n). The vector wi+1 ’wnT gi:(01'”101 1'”) )1 wi wz' 43 is known as the Gauss vector and is the pseudo Gauss vector of g,. The matrix of the form 3:3—5£ is known as the elimination matrix, for which we have (1... 0 On-0\{wl\ {w1\ 0 1 00 w,- w,- EiW: : 0 ... _w_;:i 1 ... 0 1,032+] 0 (0... —— 1/ \w/ w Here, the elimination procedure E,- w : (In — g,- ei-r) W is simply called elimination on W by gi- An n x n permutation matrix Pm- is obtained by interchanging the i-th and j-th rows of the identity matrix In. Therefore, Pm- is of the form (1 ) 1 0 0 0 1 i 0 1 - 0 0 H,j:(el, ,e,, ,e,, ,3”): ' 0 0 1 0 1 0 . 0 0 J l 44 Obiously, Pm- = In when i = j. We now describe our elimination algorithm. Denote the vector (0) w1,1 W1 : 1 (0) wn,1 and (0) w1,i+1 Wi+1: E ZAWz', 131371- (0) wn,i+1 For W1, let 0 0 l'wl £1 = max IwEBI, 11 13:3,” 1 and for the permutation matrix [3,1,1 write (1) 791,1 P51,1W1 = : (1) wn,1 (1) (0) Where 101,1 2 wi1,1' For E1 = In — g1 e? where (0) mo) 2,] w 1) 0 _ 1,1 _ 81 — — 1 . G1 112.33 1 \ wii ) We have (1 ( ‘fl \ w1,1 E P — E — 0 1 41,1 W1 — 1 — W 45 where V1 2 (10(0) ). Denote the first column of E1 by i1,1 1 T1 2 _Gl Now, for W2, we write (0), E1P41,1W2 = :1, Bl with B. = (1133,. - 92>? If Bl = 0, then V2 : (102(3) and the elimination procedure terminates. Otherwise, let 1w“) | = max lw§,12)l, ”'2 251571 then for the permutation matrix Pig; we write (0) 2 “141,2 1052) Pi2,2 : Bl where 10323 2 112(1) Let E2 = In — g2 372’" where ’ (m O 1053 82 = “’2; = . 02 mm K 1173 / Then, K 433. \ 7 V, \ we wm E2 P122 E1 P413 W2 = E2 P422 “’2 = E2 12,2 = 0 Bl I \ «15.2 1 K 0 l 46 (0) (1) )T where V2 2 (111,13, 10,2,2 . Denote the second column of E2 by _G2 and update T1 by T1 (_— E2 32,2 T1. In general, for 1 < k < n + 1, suppose 1),-1,1, E1, ..., Pm,“ En, and T1, ..., Tn are computed, let us consider the following elimination procedure on wk“. Write / 10510;“ \ Ek PinJc ' ‘ ' E1 P113 Wk+1 = (k-l) ik,k+1 KBk} w with Bk 2 (105211“, ’ ' ' 1105511“)? If Bk 2 0, then VH1 = (10:31“, - - - ,wEfiBflT and the elimination procedure termi- nates. Otherwise, let IwEleI = kpgx 1333.111. - _ T For matrices Pik+11k+1 and Ek+1 — In — gk+1 ek+1 where (91 __ 111(k) _ gk+l _ k+2,k+1 _. wik+1,k+1 woe) n,k+1 K "" / wik+1vk+1 47 we have (k-l) g wik,k+1 0 Ek+1 IDik+1,lc+1 = (k) = w(k—l) wik+1,k+1 ik1k+1 with Vk+1 = (w§?,)k+1,- - - ,wEfifl, 1031““)? Denote the (k + 1)-th column of Ek+1 by K 0 K Tk+1 = 0 1 K ”Gk+1 } and update T1, . . . ,Tk as follows, r T1 <— Ek+1Pik+1,k+1 T1, 1 Th F— Ek+1 Pin+1,k+1 Tk- Obviously, if the elimination procedure is terminated at the (k + l)-th step, then I; can be expressed as I; = Ek+1 Pin+1,k+1"'E1 P413 = Mk+1 Pk+1 where Pk-l-l : Pi1,1 Pig,2 ' ' ' ‘Rik+1,k+1 and Mk+1 = T1, ”'1 Tk+11 In—k—l 48 Remark In implementing the above elimination algorithm, we only need to store the integer list {i1,i2, . . .}, the vector lists {V1, V2, . . .}, {G1, 02, . . .} and {T1,T2, . . .}. 2 The algorithm of determining the minimal polynomials in Theorem 2.1.1 by using the elimination algorithm Based on the above elimination algorithm, we use the following pseudo code to de- scribe our algorithm of determining the minimal polynomials \I11()\), \I12(A), . . ., \I't()\) of an n x n matrix A in Theorem 2.1.1. Input: An n X n matrix A. Output: The minimal polynomials \IJ1(A), \I!2()\), ..., \I’t(/\). t = O (the number of minimial polynomials); m = n; (Q) P: 0 (the list of permutation indices {i1, 72, . . .}); L: (I) (the list of the pseudo Gauss vectors {01, G2,. . .}); R: (l) (the list of {V1, V2,. . .} generated by elimination); T: (l) (the list of the columns {T1, T2, . . .} of the inverse matrix ); generate a random vector x = (:31, - - . ,ccm)T E 72’”; v = w = x; r = 0 (the degree of a minimal polynomial); For i = 1 to m + 1 For j = 1 to i — 1 If i, > j, then interchange 7113- and 11173., i.e., wj 1:? wt,- Endif update w by using Gj; V.- = V.- U {wj} 49 End If i i, then interchange w,- and win, i.e., w, (32? wk“ and the corresponding elements in T1, . . . ,T;_1 Endif V,=V,~ U {10,-},andR2RU {Vi}; generate 0;, and L = L U {0,}; update T1, . . . ,T,~_1 by using G1; generate T,- by using 0,, and T = T U {7}}; generate the next vector w 2 Av; v=w; r=r+1 Endif End t=t+1; (Jo) generate the minimal polynomial \I’t()\) of x; If m—r>0, then If \Ilt()\) is not the minimal polynomial of 72’", then 50 M) (4) Endif form the matrices M2, M4 and B1, Bz; generate the new matrix A = M2 31 + M4 82 in Rm_'; m = m — r; g_0tg (4) Endif stop. 3 Computational example We use the following 8 x 8 Hadamard matrix (11111111) K1—1—1 1—1 1 1—1} to demonstrate our elimination algorithm described above in computing its Jordan normal form. Since all the entries of this matrix are integers, we may proceed our computation by the symbolic computation. And, as we mentioned before, in this case the decision of row pivoting will base on whether or not the pivot element is zero rather than how big its magnitute is. Let A1 = A. Choose a nonzero vector __ 1 1 1 1 1 1 1 T 8 X-(1151514155715) ER: andletv=w=x. 51 Since the first entry of w is 1, pivoting is not necessary. Using Gauss vector 0 g1 = G] with the pesudo Gauss vector _ 1 1 1 1 1 1 1 T GI — (515121515135) on the elimination of w, i.e., w <— (18 — g1 eff)w, we have w = (1,0,0,0, 0,0,0,O)T and 1 S H A ...: v .5 ll _Gl Compute _ -313323231211217w2g3flT W _ A1V _ (280’ 840’ 840’ 280’ 840 ’ 280’ 280’ 840) and let v = w. followed by updating w by using g1 _ T _fl_1217§_27719011_9flfl7‘ we (18 glellw_(2301 133012101 1120’2100’240’980’960) ' 1217 7630’ pivoting is not needed. Elimination on w by Since the second entry of w is — the Gauss vector II o 82 6'2 where G = (_ 184 831 7604 133 5052 413) 2 1217’ 2434’ 6085 1217’ 8519’ 4868 and updating T1 by 13 — g2 6;, yields w <— (Is—Ezele‘”:(%1—%101010101010)T1 E _ T _ _1 _1493 _ 193 _5019 _ 303 _3743 _ 315 T T1 ([8 g232)T1 —( 1 21 3351’ 24347 6085’ 3651’ 8519’ 4333) 1 52 and hence 0 22 280 W = , and T2 2 1 _1217 1680 _Gz Next, we calculate and let v = w. Updating w by g1 and g2, we have W <— (I8 _ g2 8;)(I3 _— gl 8:111) W : (810101010101O1O)T1 and, therefore the elimination procedure for A1 should be terminated here. Based on V1 and V2, we solve the linear system 731 1 m Go 8 _ 1217 0 1330 31 0 to obtain the minimal polynomial \II(A) = A2 — 8 of x. Since pivoting is not used in the elimination process of the above two vectors, we may choose {33}, . . . , {38} as the basis of the quotient space R8 / W1 where W1 :— span{x,A1x}. Since ‘II(A1)e,- = O for i = 3,4, 5,6, 7,8, \I’1(A) = \II()\) is the minimal polynomial of R8. From _ 1493 134 3651 1217 0 _ 193 _ 331 4 _ _ 2434 2434 [8 _ T11 T21 "" 1 I _5019 7304 6 6085 6085 I 808 133 6 " 3351 121 _ 3743 505 3519 3516 __ 315 413 K 4333 4333' 53 and the notations in (4.8), we have { 1493 184 \ _3351 1217 __ 193 _ 831 2434 2434 _5019 7304 I 6085 6085 (Zhié,IVI4) :: 6 _303 133 3651 1217 _3743 5052 8519 8519 _315 413 4868 4868 It follows from (4.19) f 1 1 1 1 1 ll 1 —1 1 —1 1 —1 —1 —1 1 1 —1 —1 Bl 1 — 1 1 —1 —1 l (14163114164,A165,A1€3,A1€7,A1€3) = = 32 1 1 —1 —1 —1 —1 1 —1 —1 1 —1 1 —1 —1 —1 —1 1 1 \—1 1 —1 1 1 —1) Therefore, B1 A2 : (M2,M4) = M2B1+B2 32 {_4592 _5393 2710 1303 _4592 _5393 \ 3651 3651 3651 3651 3651 3651 _1729 1533 705 _ 393 _1729 1533 17 1217 1217 1217 1217 1217 1734 _3533 700 _13703 _ 700 _13703 __ 1217 6085 1217 6085 1217 6085 3242 _4353 _4030 2444 _4030 2444 3651 3651 3651 3651 3651 3651 _1030 _17314 _1030 _17314 1404 _273 1217 8519 1217 8519 1217 8519 _2335 910 _2335 910 2233 _1524 K 2434 1217 2434 1217 2434 1217 Next, choose a nonzero vector _ 1 1 1 1 1 1 T R6 ){'—' )2, §114151 6 ES 1 andletv=w=x. 54 Since the first entry of w is 1, again no pivoting is needed. Elimination on w by the Gauss vector 81 where i.e., W <— (13 — g1 ef)w, we obtain W = (1, 0, 0,0, 0,0)T, and hence 1 m = (1), and T1 = —Gi Compute _30029 __ 10403 12732 334 _ 320377 __ 13573)T 335101 121701 132551 33511 25555701 13255 W=A2V=( and let v = w. Updating w by g1 yields _ T _ _30029 17311 3337 22223 _1270341 _113347 T W i (I6 3161 )W _ ( 3351017302011095301433301 3339251 219030) Since the second entry of w is %, again no pivoting is needed. Elimination on W by 0 32 = O 02 with G _ (3374 33339 _432124 _113347)T 2 — 523331352221 530351 52333 and updating T1 by I6 — g2 6; yields T _ 30029 17311 W <— (13—g262)W—(‘fifi5’fifi5’0’0’0’m7 _ T _ _1 _ 4753 24592 _242239 _22743 T T1 E <16 g232)T1 - (1’ 21 173111352221 530351 17311) ° 55 It follows that _ 30029 V2 36510 17611 73020 Compute w = A2V : (8,4, and let v = w. Updating w by g1 and g2 yields W (— (16 — g2 6;)(I6 — g1 6?) W : (87010103010)T Solving the linear system 1 () _ 30029 36510 17611 73020 results in the minimal polynomial \II()\) = A2 — 8 of x. It is easy to see that the basis of the quotient space R6/W1 is {33}, {e4}, {e5}, {es} where W1 2 span{x,A2x}. Since ‘II(A2)e,- = 0 for z' = 3,4,5,6, 412(k) = @(A) is the minimal polynomial of R6. Again, from 1: 6 : T1) T2) we have (M2,M4) = \ f _ 1 0 1 0 __§ 1 __ 4753 __ 3374 __ 17311 52333 ‘ , 24529 __33339 35222 35222 _ 242239 432124 14 53035 53035 __22743 113347 17311 52333 4753 __ 3374 \ 17311 52333 24529 __33339 35222 35222 4 __242239 432124 53035 53035 __22743 113347 ) 17311 52333 5&3 in (4.8), and from (4.19), { 2710 1303 __4592 __5393 \ 3351 3351 3351 3351 705 __ 393 __1729 1533 1217 1217 1217 1217 131 __ 700 __13703 __ 700 __13703 __ __ 1217 3035 1217 3035 (A283, A264, A265, A263) — — 13 __4030 2444 __4030 2444 2 3351 3351 3351 3351 __1030 __17314 1404 __ 273 1217 3519 1217 3519 __2335 910 2233 __1524 K 2434 1217 2434 1217 Therefore, 81 143 Z: (1VL2,]V£4) =2 .A45.E31'+‘132 B2 { __14950 __313932 __ 2954 __247344 \ 7311 234135 52333 33055 __173745 125342 74095 __143252 __ 105333 52333 105333 52333 24272 __43124 __24313 2373223 33321 4303 4303 133105 \ ___25995 __73202 __ 39125 33423 ) 35222 52333 105333 17311 Following exactly the same procedures as described above, the minimal polynomial ‘II3()\) = A2 — 8 of ”R4 with respect to A3 and the matrix __1702102 688850416 398853 41879565 144 =2 _ 47135 _ 1702102 75972 398853 can be achieved. And the minimal polynomial \I’4()\) = A2 — 8 of R2 with respect to A4 can be determined similarly. For \Ill()\), \I/2()\), \I/3(/\) and \II4()\), a common squarefree relatively prime basis {(A — 2\/§), (A + 2J5” in 72D] can be found by means of symbolic computation such that 91(3) = 92(1) = 93(1) 2 94(1) = (A — 275m + 27'2"). Therefore, it follows from Theorem 2.1.3 that the Jordan normal form J of A is of 57 the form ol (Ni —2\/§ 0 0 —2\/§ 0 0 0 0 —2\/§ 272 0 0 —2\/§) 0 58 BIBLIOGRAPHY 59 BIBLIOGRAPHY [1] E. Anderson, et al., LAPACK user’s guide, SIAM Philadelphia, 1992. [2] A. Ralston and P. Rabinowitz, A first course in numerical analysis (Second edition), McGraw-Hill, 1978. [3] T. Beelen and P. Van Dooren, Computational aspects of the Jordan canonical form, in Reliable Numerical Computation, M.G. Cox and S. Hammerling, eds., Clarendon Press, Oxford, 1990, pp. 57—72. [4] J. Demmel, A numerical analyst’s Jordan canonical form, Ph.D. Dissertation, Computer Science Division, University of California, Berkeley, 1983. [5] FR. Gantmacher, The theory of matrices, Vol. 1, Chelsea, New York, 1959. [6] J. von zur Gathen, Parallel algorithms for algebraic problems, SIAM J. Comput., 13 (1934), pp. 302—324. [7] CH. Golub and CF. Van Loan, Matrix Computation (First edition), The Johns Hopkins University Press, 1983. [8] CH. Golub and CF. Van Loan, Matrix Computation (Third edition), The Johns Hopkins University Press, 1996. [9] CH. Golub and J .H. Wilkinson, Ill-conditioned eigensystems and the computa- tion of the Jordan canonical form, SIAM Rev., 18 (1976), pp. 578-619. [10] AS. Householder, The theory of matrices in numerical analysis, Ginn (Blaisdell), Boston, 1964. [11] B. Kégstréim and A. Ruhe, An algorithm for numerical computation of the Jordan normal form of a complex matrix, ACM Trans. Math. Soft., 6 (1980), pp. 398- 419. [12] E. Kaltofen, Sparse Hensel lifting, in Proceedings of EUROCAL’85, Vol. 2, Lec- ture Notes in Comput. Sci. 204, Springer-Verlag, pp. 4—17, 1985. 60 [13] E. Kaltofen and M.F. Krishnamoorthy and ED. Saunders, Parallel algorithms for matrix normal forms, Linear algebra appl., 136 (1990), July 15, pp. 189-208. [14] T.Y. Li, Zhinan Zhang and Tianjun Wang, Determine the structure of the Jor- dan normal form of a matrix by symbolic computation, Linear Algebra and Its Applications, 252 (1997), pp. 221-259. [15] S. MacLane and C. Birkhoff, Algebra, The Macmillan Company, New York, 1967. [16] ET. Smith, et al., Matrix eigensystem routines - EISPACK guide (Second edi- tion), Springer, New York, 1976. [17] J.H. Wilkinson, The algebraic eigenvalue problem, Clarendon Press, Oxford, 1965. [18] S. Wolfram, Mathematica: A system for doing mathematics by computer, Addison-Wesley, Redwood City, California, 1991. 61