The Square Root Function of a Matrix

Georgia State University Digital Archive @ GSU Mathematics Theses Department of Mathematics and Statistics 4-24-2007 The Square Root Function of a Matrix Crystal Monterz Gordon [email protected] This Thesis is brought to you for free and open access by the Department of Mathematics and Statistics at Digital Archive @ GSU. It has been accepted for inclusion in Mathematics Theses by an authorized administrator of Digital Archive @ GSU. For more information, please contact [email protected] .Recommended Citation Gordon, Crystal Monterz, "The Square Root Function of a Matrix" (2007). Mathematics Theses. Paper 24. http://digitalarchive.gsu.edu/math_theses/24 THE SQUARE ROOT FUNCTION OF A MATRIX by Crystal Monterz Gordon Under the Direction of Marina Arav and Frank Hall ABSTRACT Having origins in the increasingly popular Matrix Theory, the square root func- tion of a matrix has received notable attention in recent years. In this thesis, we discuss some of the more common matrix functions and their general properties, but we specifically explore the square root function of a matrix and the most effi- cient method (Schur decomposition) of computing it. Calculating the square root o fa2×2 matrix by the Cayley-Hamilton Theorem is highlighted, along with square roots of positive semidefinite matrices and general square roots using the Jordan Canonical Form. Keywords : Cayley-Hamilton Theorem, Interpolatory Polynomials, Jordan Canon- ical Form, Matrix Theory, Functions of Matrices, Positive Semidefinite Matrices, Schur’s Theorem, Square Roots of Matrices THE SQUARE ROOT FUNCTION OF A MATRIX by Crystal Monterz Gordon A Thesis Presented in Partial Fulfillment of the Requirements for the Degree of Master of Science in College of Arts and Sciences Georgia State University 2007 Copyright by Crystal Monterz Gordon 2007 THE SQUARE ROOT FUNCTION OF A MATRIX by Crystal Monterz Gordon Major Professors: Marina Arav and Frank Hall Committee: Rachel Belinsky Zhongshan LiMichael Stewart Electronic Version Approved: Office of Graduate Studies College of Arts and Sciences Georgia State University May 2007 ACKNOWLEDGEMENTS The author wishes to gratefully acknowledge the assistance of Dr. Marina Arav, Dr. Frank J. Hall, Dr. Rachel Belinsky, Dr. Zhongshan Li, and Dr. Michael Stew- art without whose guidance this thesis would not have been possible. She would also like to thank Drs. George Davis, Mihaly Bakonyi, Johannes Hattingh, Lifeng Ding, Alexandra Smirnova, Guantao Chen, and Yichuan Zhao for their support and encouragement in her course work and in her research towards this thesis. iv TABLE OF CONTENTS ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . iv 1 . I n t r od u c t i o n ................................ 1 2. Functions of Matrices . . . . . . . . . . . . . . . . . . . . . . . 5 3. The Square Root of a 2 ×2 Matrix . . . . . . . . . 12 4. Positive Semidefinite Matrices . . . . . . . . . . . . . . 175. General Square Roots . . . . . . . . . . . . . . . . . . . . . . 206. Computational Method . . . . . . . . . . . . . . . . . . . . . 28R e f e r e n c e s .................................... 3 7 v 1 1. Introduction As stated in [1] and [19], the introduction and development of the notion of a matrix and the subject of linear algebra followed the development of determinants. Gottfried Leibnitz, one of the two founders of calculus, used determinants in 1693 arising from the study of coefficients of systems of linear equations. Additionally, Cramer presented his determinant-based formula, known as Cramer’s Rule, for solving systems of linear equations in 1750. However, the first implicit use of matrices occurred in Lagrange’s work on bilinear forms in the late 1700’s in his method now known as Lagrange’s multipliers. Some research indicates that the concept of a determinant first appeared between 300 BC and AD 200, almost 2000 years before its invention by Leibnitz, in the Nine Chapters of the Mathematical Art by Chiu Chang Suan Shu. There is no debate that in 1848 J.J. Sylvester coined the term, “matrix”, which is the Latin word for womb, as a name for an array of numbers. Matrix algebra was nurtured by the work of Arthur Cayley in 1855. He studied compositions of linear transformations and was led to define matrix multiplication, so that the matrix of coefficients for the composite transformation ABis the product of the matrix Atimes the matrix B. Cayley went on to study the algebra of these compositions including matrix inverses and is famous for the Cayley-Hamilton theorem, which is presented later in this thesis. In mathematics, a matrix is a rectangular table of numbers, or more generally, a table consisting of abstract quantities. Matrices are used to describe linear equa- tions, keep track of coefficients of linear transformations, and to record data that depend on two parameters. Matrices can be added, multiplied, and decomposed in various ways, which makes them a key concept in linear algebra and matrix theory, two of the fundamental tools in mathematical disciplines. This makes intermediate facts about matrices necessary to understand nearly every area of mathematical 2 science, including but not limited to differential equations, probability, statistics, and optimization. Additionally, continuous research and interest in applied mathe- matics created the need for the development of courses devoted entirely to another key concept, the functions of matrices. In this thesis, we provide a detailed overview of the basic functions of matrices while focusing on the square root function of a matrix and a few of the most common computational methods. We discuss the specific case of a square root of a 2 ×2 matrix before outlining results on square roots of positive semidefinite matrices and general square roots. Although the theory of matrix square roots is rather complicated, simplifica- tion occurs for certain classes of matrices. Consider, for example, symmetric pos- itive semi(definite) matrices. Any such matrix has a unique symmetric positive semi(definite) square root, and this root finds use in the theory of the generalized eigenproblem [16] (section 15-10), and preconditioned methods [4, 10]. More gener- ally, any matrix Ahaving no nonpositive real eigenvalues has a unique square root, for which every eigenvalue has a positive real part, and it is this square root, de- notedA1 2and sometimes called the principal square root , that is usually of interest (e.g. the application in boundary value problems, [17]). There is a vast amount of references available focusing on the square root func- tion of a matrix, many of which are listed in the References section. While some of the references were used explicitly, all provided insight and assistance in the completion of this thesis. We begin now by defining key terms used throughout this thesis for clarity and cohesiveness. Definitions As in [8] and [9], we let Mndenote the set of all n×ncomplex matrices. We note that some authors use the notation Cn×n. Now letA∈Mn. Then a nonzero 3 vectorx∈Cnis said to be an eigenvector ofAcorresponding to the eigenvalueλ, if Ax=λx. The set containing all of the eigenvalues of Ais called the spectrum of Aand is denoted,σ(A). LetA,B∈Mn. ThenBis asquare root ofA,i fB2=A. A matrixD=[dij]∈Mnis called a diagonal matrix ,i fdij= 0 whenever i/negationslash=j. LetA,B∈Mn. ThenAissimilar toB, denotedA∼B, if there is a nonsingular matrixSsuch thatS−1AS=B.I fA∼B, then they have the same characteristic polynomial and therefore the same eigenvalues with the same multiplicities. LetA∈Mn. ThenAisdiagonalizable ,i fAis similar to a diagonal matrix. A matrixU∈Mnis said to be unitary ,i fU∗U=I. A matrixA∈Mnis said to be unitarily equivalent orunitarily similar to B∈Mn, if there is an unitary matrix U∈Mnsuch thatU∗AU=B.I fUmay be taken to be real (and therefore real orthogonal), then Ais said to be (real) orthogonally equivalent toB. If a matrix A∈Mnis unitarily equivalent to a diagonal matrix, Ais said to be unitarily diagonalizable . LetA∈Mn. ThenAisHermitian ,i fA∗=A, whereA∗=¯AT=[ ¯aji]. If A∈Mnis Hermitian, then the following statements hold: (a) All eigenvalues of Aare real; and (b)Ais unitarily diagonalizable. 4 The minimal polynomial of A , denotedm(t), is the monic annihilating polyno- mial of the least possible degree. Ann×nmatrixAis called upper triangular ,i faij= 0 fori>j , i.e. all of the entries below the main diagonal are zero. Ann×nmatrixAis called lower triangular ,i faij= 0 fori<j , i.e. all of the entries above the main diagonal are zero. The functions of matrices appear widely in many areas of linear algebra and are linked to numerous applications in both science and engineering. While the most common matrix function is the matrix inverse (usually mentioned with terms: in- vertible or nonsingular), other general matrix functions are the matrix square root, the trigonometric, the exponential and the logarithmic functions. The following are the definitions of the matrix functions mentioned above. Examples of General Matrix Functions A matrixAisinvertible ornonsingular , if there exists a unique inverse denoted byA−1, whereA−1A=IandAA−1=IandIis the identity matrix. Letp(t)=aktk+···+a1t+a0be a polynomial. Then, by definition, p(A)=akAk+···+a1A+a0I. The exponential of A∈Mn, denotedeAorexp(A), is defined by eA=I+A+A2 2!+···+Ak k!+···. LetA∈Mn.A n yXsuch thateX=Ais alogarithm ofA. The sineandcosine ofA∈Mnare defined by cos(A)=I−A2 2!+···+(−1)k (2k)!A2k+···, sin(A)=A−A3 3!+···+(−1)k (2k+ 1)!A2k+1+···. 5 2. Functions of Matrices We provide a detailed overview of the basic ideas of functions of matrices to aid the reader in the understanding of the “connectivity” of the fundamental principles (many of which are defined in the introduction) of matrix theory. One can easily show that if Ax=λxandp(t) is a polynomial, then p(A)x= p(λ)x, so that ifxis an eigenvector of Acorresponding to λ, thenxis an eigenvector ofp(A) corresponding to the eigenvalue p(λ). We will shortly obtain an even stronger result. Perhaps the most fundamentally useful fact of elementary matrix theory is that any matrix A∈Mnis unitarily equivalent to an upper triangular (also to a lower triangular) matrix T. Representing the simplest form achievable under unitary equivalence, we now recall one of the most useful theorems in all of matrix theory, Schur’s Theorem. Schur’s Theorem: IfA∈Mn, thenAis unitarily triangularizable, that is, there exists a unitary matrix Uand an upper-triangular matrix Tsuch thatU∗AU=T. Through the use of Schur’s Theorem, one can prove that if A∈Mnwithσ(A)= {λ1,...,λ n}andp(t) is a polynomial, then σ(p(A)) ={p(λ1),...,p (λn)}. The proof goes as follows: U∗p(A)U=p(U∗AU)=p(T), which is upper- triangular with p(λ1),...,p (λn) on the diagonal. The proof follows from the simi- larity ofp(A) andp(T). We now shift our focus from polynomials to general functions.LetA∈M nand suppose that λ1,λ2,...,λsare the distinct eigenvalues of A,s o that m(t)=(t−λ1)m1(t−λ2)m2···(t−λs)ms 6 is the minimal polynomial of Awith degree m=m1+m2+ ...+ms. Thenmk is the index of the eigenvalue λk,i.e. it is the size of the largest Jordan block associated with λkand is equal to the maximal degree of the elementary divisors associated with λk(1≤k≤s). Now, a function f(t)i sdefined on the spectrum of A, if the numbers f(λk),f/prime(λk),...,f(mk−1)(λk),k =1,2,...,s, are defined (exist). These numbers are called the values of f (t)on the spectrum of A, where ifmk=1 ,f(mk−1)isf(0)or simplyf. Many of the succeeding results can be found in [12], but we will provide more details here. Proposition 2.1: Every polynomial is defined on the spectrum of any matrix inMn. For the polynomial m(t), the values of m(λk),m/prime(λk),...,m(mk−1)(λk),k =1,2,...,s, are all zero. Proof: The first statement is clear. Next, each m(λk)=0.So, m/prime(t)=(t−λ1)m1d dt[(t−λ2)m2···(t−λs)ms]+[(t−λ2)m2···(t−λs)ms]·m1(t−λ1)m1−1. Therefore, m/prime(λ1)=0·d dt[(t−λ2)m2···(t−λs)ms]+[(t−λ2)m2···(t−λs)ms]·0=0,ifm1>1. Similarly, for the other λkand the higher order derivatives. Proposition 2.2: For the two polynomials p1(t) andp2(t),p1(A)=p2(A)i f and only if p1(t) andp2(t) have the same values on the spectrum of A. 7 Proof: ⇒Supposep1(A)=p2(A).Letp0(t)=p1(t)-p2(t). Then,p0(A)=0 . So,m(t) is a factor of p0(t), i.e.p0(t)=q(t)m(t) for some polynomial q(t). Now, each term of p(j) 0(t) is a product, which involves one of the terms: m(t),m/prime(t),...,m(j)(t). Hence, by Proposition 2.1, p(j) 1(λk)−p(j) 2(λk)=p(j) 0(λk)=0, forj=0,1,...,m k−1,and 1 ≤k≤s.So,p(j) 1(λk)=p(j) 2(λk) for the values of j andk. ⇐We assume that p1(t) andp2(t) have the same values on the spectrum of A. Letp0(t)=p1(t)−p2(t), then p(j) 0(λk) = 0 forj=0,1,2,...,m k−1. So,λkis a zero of p0(t) with multiplicity of at least mk, i.e. (t−λk)mkis a factor ofp0(t). Hence,m(t) is a factor of p0(t), wherep0(t)=q(t)m(t) and therefore, p0(A)=0.Thus,p1(A)=p2(A). Proposition 2.3 (Interpolatory Polynomial): Given distinct numbers λ1,λ2,...,λ s, positive integers m1,m2,...,m swithm=s/summationdisplay k=1mk, and a set of numbers fk,0,f k,1,..., f k,m k−1,k=1,2,...,s, there exists a polynomial p(t) of degree less than msuch that p(λk)=fk,0,p(1)(λk)=fk,1,..., p(mk−1)(λk)=fk,m k−1fork=1,2,...,s. (1) Proof: It is easily seen that the polynomial pk(t)=αk(t)ψk(t) (note: ifs=1 , then by definition ψ1(t)≡1), where 1 ≤k≤sand αk(t)=αk,0+αk,1(t−λk)+···+αk,m k−1(t−λk)mk−1, 8 ψk(t)=s/productdisplay j=1,j/negationslash=k(t−λj)mj, has degree less than mand satisfies the conditions pk(λi)=p(1) k(λi)=···=p(mi−1) k (λi)=0 fori/negationslash=kand arbitrary αk,0,αk,1,···,αk,m k−1.Hence, the polynomial p(t)=p1(t)+p2(t)+···+ps(t) (2) satisfies conditions (1) if and only if pk(λk)=fk,0,p(1) k(λk)=fk,1,..., p(mk−1) k (λk)=fk,m k−1for each 1 ≤k≤s.(3) By differentiation, p(j) k(λk)=j/summationdisplay i=0/parenleftbiggj i/parenrightbigg α(i) k(λk)ψ(j−i) k(λk) for 1≤k≤s,0≤j≤mk−1.Using Eqs.(3) and recalling the definition of αk(λ), we have for k=1,2,...,s ,j=0,1,...,m k−1, fk,j=j/summationdisplay i=0/parenleftbiggj i/parenrightbigg i!αk,iψ(j−i) k(λk). (4) Sinceψk(λk)/negationslash= 0 for each fixed k, Eqs. (4) can now be solved successively (beginning withj= 0) to find the coefficients αk,0,...,α k,m k−1for which (3) holds. Thus, a polynomial p(t) of the form given in (2) satisfies the required conditions. The interpolatory polynomial referred to in Proposition 2.3 is known as the Her- mite interpolating polynomial . It is in fact unique, but the proof of the uniqueness is omitted, since it is quite cumbersome. If f(t) is defined on the spectrum of A, we definef(A)t ob ep(A), wherep(t) is the interpolating polynomial for f(t)o n the spectrum of A. Theorem 2.4: IfA∈Mnis a block-diagonal matrix, A= diag[A1,A2,...,A t], 9 and the function f(t) is defined on the spectrum of A, then f(A) = diag[f(A1),f(A2),...,f (At)]. (5) Proof : It is clear that for any polynomial q(t), q(A) = diag[q(A1),q(A2),...,q (At)]. Hence, ifp(t) is the interpolatory polynomial for f(t) on the spectrum of A, we have f(A)=p(A) = diag[p(A1),p(A2),...,p (At)]. Since the spectrum of Aj(1≤j≤t) is obviously a subset of the spectrum of A, the function f(t) is defined on the spectrum of Ajfor eachj=1,2,...,t . (Note also that the index of an eigenvalue of Ajcannot exceed the index of the same eigenvalue of A.) Furthermore, since f(t) andp(t) assume the same values on the spectrum of A, they must also have the same values on the spectrum of Aj(j=1,2,...,t ). Hence, f(Aj)=p(Aj) and we obtain Eq. (5). Theorem 2.5: IfA,B,S ∈Mn, whereB=SAS−1, andf(t) is defined on the spectrum of A, then f(B)=Sf(A)S−1. (6) Proof : Since A and B are similar, they have the same minimal polynomial. Thus, ifp(t) is the interpolatory polynomial for f(t) on the spectrum of A, then it is also the interpolatory polynomial for f(t) on the spectrum of B. Thus, we have f(A)=p(A), f(B)=p(B), 10 p(B)=Sp(A)S−1, so the relation (6) follows. Theorem 2.6: LetA∈Mnand letJ= diag[J j]t j=1be the Jordan canonical form of A, where A=SJS−1andJjis thejthJordan block of J. Then f(A)=Sdiag[f(J1),f(J2),...,f (Jt)]S−1. (7) The last step in computing f(A) by use of the Jordan form of Aconsists of the following formula. Theorem 2.7: LetJ0be a Jordan block of size lassociated with λ0: J0= λ 01 λ0... ...1 λ0 . Iff(t) is an (l−1)-times differentiable function in a neighborhood of λ 0, then f(J0)= f(λ 0)1 1!f/prime(λ0)...1 (l−1)!f(l−1)(λ0) 0f(λ0)...... .........1 1!f/prime(λ0) 0... 0f(λ0) . (8) Proof: The minimal polynomial of J0is (t−λ0)land the values of f(t)o n the spectrum of J0are therefore f(λ0),f/prime(λ0),. . . ,f(l−1)(λ0). The interpolatory polynomial p(t), defined by the values of f(t) on the spectrum {λ0}ofJ0, is found by putting s=1,mk=1,λ1=λ0, andψ1(t)≡1. One obtains p(t)=l−1/summationdisplay i=01 i!f(i)(λ0)(t−λ0)i. The fact that the polynomial p(t) solves the interpolation problem p(j)(λ0)= f(j)(λ0),1≤j≤l−1, can also be easily checked by a straightforward calcula- tion. 11 We then have f(J0)=p(J0) and hence f(J0)=l−1/summationdisplay i=01 i!f(i)(λ0)(J0−λ0I)i. Computing the powers of J0−λ0I, we obtain (J0−λ0I)i= 01 0... ...1 0 i = 0...01 0...0 ...... ... 0 ... 1 0 ... 0  with 1’s in the i-th super-diagonal positions, and zeros elsewhere, and Eq.(8) follows. Thus, given a Jordan decomposition of the matrix A, the matrix f(A) is easily found by combining Theorems 2.6 and 2.7. From Theorems 2.6 and 2.7, we have the following results. Theorem 2.8: Using the notation of Theorem 2.6, f(A)=Sdiag[f(J1),f(J2),...,f (Jt)]S−1, wheref(Ji)(i=1,2,...,t ) are upper triangular matrices of the form given in Eq.(8). Theorem 2.9: Ifλ1,λ2,...,λ nare the eigenvalues of the matrix A∈Mnandf(t) is defined on the spectrum of A, then the eigenvalues of f(A) aref(λ1),f(λ2),...,f (λn). This follows from the fact that the eigenvalues of an upper triangular matrix are its diagonal entries. 12 3. The Square Root of a 2×2Matrix IfA,B∈MnandAis similar to B, thenAhas a square root if and only if B has a square root. The standard method for computing a square root of an n×n diagonalizable matrix Ais easily stated. Suppose S−1AS=D for some nonsingular matrix Sand diagonal matrix D. Then A=SDS−1, and by substitution we have A=(SˆDS−1)(SˆDS−1)=SDS−1, where ˆDis a square root of D. In general, the matrix Dwill have 2ndistinct square roots (when Ahasnnonzero eigenvalues, which are obtained by taking the square roots of the diagonal elements of Dwith all possible sign choices). If D1/2 is any square root of D, it follows that B=SD1/2S−1is a square root of A, that isB2=A.However, even in some 2 ×2 cases, the computations can become quite messy. Not every 2 ×2 matrix has a square root. For example, by direct calculation, we can show that/bracketleftbigg01 00/bracketrightbigg has no square root. On the other hand, if b∈C, /bracketleftbigg bb −b−b/bracketrightbigg gives an infinite number of square roots of /bracketleftbigg00 00/bracketrightbigg 13 and ifb/negationslash=1 , 1 b−1/bracketleftbigg b1 11/bracketrightbigg/bracketleftbigg −10 01/bracketrightbigg/bracketleftbigg 1−1 −1b/bracketrightbigg =1 b−1/bracketleftbigg −b−12b −2b+1/bracketrightbigg yields an infinite number of square roots of /bracketleftbigg10 01/bracketrightbigg . We next recall another useful theorem in matrix analysis, the Cayley-Hamilton Theorem. Cayley-Hamilton Theorem: IfA∈MnandpA(t) = det(tI−A) is the char- acteristic polynomial of A, thenpA(A)=0. In [15] and [18], the authors show how the Cayley-Hamilton Theorem may be used to determine explicit formulae for all the square roots of 2 ×2 matrices. These formulae indicate exactly when a 2 ×2 matrix has square roots, and the number of such roots. Suppose Ais 2×2 and X2=A. (9) However, for each 2 ×2 matrtixX, the Cayley-Hamilton Theorem states that X2−(trX)X + (det X)I = 0 . (10) Thus, if a 2 ×2 matrixAhas a square root X, then we may use (10) to eliminate X2from (9) to obtain tr(X)X = A + (det X)I . Now, since (det X)2= detX2= detA, then detX=/epsilon11√ detA, /epsilon1 1=±1, that is det√ A=/epsilon11√ detA, so that the above result simplifies to the identity: (trX)X = A + /epsilon11√ det AI,/epsilon11=±1. (11) 14 Case 1:Ais a scalar matrix. IfAis a scalar matrix, A=aI, then (11) gives (trX)X = (1 + /epsilon11)aI,/epsilon11=±1. Hence, either (trX)X = 0, or (trX)X = 2aI. The first of these possibilities deter- mines the general solution of (9) as X=/parenleftbigg αβ γ−α/parenrightbigg ,α2+βγ=a, (12) and it covers the second possibility, if a= 0. On the other hand, if a/negationslash= 0, then the second possibility, (trX)X = 2aI, implies Xis scalar and has only one pair of solutions X=±√ aI. (13) For this case, we conclude that if Ais a zero matrix, then it has a double-infinity of square roots as given by (12) with a= 0, whereas if Ais a nonzero, scalar matrix, then it has a double-infinity of square roots plus two scalar square roots as given by (12) and (13). Case 2:Ais not a scalar matrix. IfAis not a scalar matrix, then trX /negationslash=0i n (11). Consequently, every square root Xhas the form: X=τ−1(A+/epsilon11√ detAI),τ/negationslash=0. Substituting this expression for Xinto (9) and using the Cayley-Hamilton theorem forA,w efi n d A2+( 2/epsilon11√ detA−τ2)A+ (detA)I=0 ((trA)A −(det A)I) + (2 /epsilon11√ det A−τ2)A + (det A)I = 0 (trA+2/epsilon11√ detA−τ2)A=0. SinceAis not a scalar matrix, then Ais not a zero matrix, so τ2=trA+2/epsilon11√ detA,(τ/negationslash=0,/epsilon11=±1). (14) 15 If (trA)2/negationslash= 4 det A, then both values of /epsilon11may be used in (14) without reducing τ to zero. Consequently, it follows from (11) that we may write X, the square root ofA,a s X=/epsilon12A+/epsilon11√ detAI /radicalbig trA+2/epsilon11√ detA. (15) Here each/epsilon1i=±1, and if det A/negationslash= 0, the result determines exactly four square roots forA. However, if det A= 0, then the result (15) determines two square roots for Aas given by X=±1 √ trAA. (16) Alternatively, if ( trA)2= 4 detA/negationslash= 0, then one value of /epsilon11in (14) reduces τto zero, whereas the other value yields the results 2 /epsilon11√ detA= trA andτ2= 2 trA.In this case,Ahas exactly two square roots given by X=±1 √ 2trA(A+1 2(trA)I). (17) Finally, if (trA)2= 4 det A = 0, then both values of /epsilon11reduceτto zero in (14). Hence it follows by contradiction that Ahas no square roots. For this case, we conclude that a nonscalar matrix, A, has square roots, if and only if, at least one of the numbers, trA and det A, is nonzero. Then the matrix has four square roots given by (15), if (trA)2/negationslash= 4 det A,det A/negationslash=0 and two square roots given by (16) or (17), if (trA)2/negationslash= 4 det A,det A = 0 or (trA)2= 4 det A,det A/negationslash=0. It is worth noting from (15) that tr X = tr√ A=/epsilon12/radicalBig trA + 2/epsilon11√ det A. 16 Hence using the identity, det√ A=/epsilon11√ detAas applied in (11), result (15) may be rewritten as √ A=1 tr√ A(A+ det√ AI), which is equivalent to the Cayley-Hamilton Theorem for the matrix√ A. This same deduction can be made, of course, for all other cases under which√ Aexists. In [2], the author is concerned with the determination of algebraic formulas yielding all of the solutions of the matrix equation Bn=A, wherenis a positive integer greater than 2 and Ais a 2×2 matrix with real or complex elements. If Ais a 2 ×2 scalar matrix, the equation Bn=Ahas infinitely many solutions, and one can obtain the explicit formulas giving all of the solutions. If A is a non- scalar matrix, the equation Bn=Ahas only a finite number of solutions. While the author’s concern is beyond the scope of this thesis, it outlines a process for obtaining other roots with expressed preciseness. 17 4. Positive Semidefinite Matrices By definition, an n×nHermitian matrix A is called positive semidefinite ,i f x∗Ax≥0 for all nonzero x∈Cn. Theorem 4.1: A Hermitian matrix A∈Mnis positive semidefinite if and only if all of its eigenvalues are nonnegative. Proof: SinceAis Hermitian, there exists a unitary matrix Uand a diagonal matrixDsuch thatU∗AU=D.Then x∗Ax=x∗UDU∗x=y∗Dy=n/summationdisplay i=1di¯yiyi=n/summationdisplay i=1di|yi|2, whereU∗x=y. ⇒Let the Hermitian matrix A∈Mnbe positive semidefinite. Then, from the above, n/summationdisplay i=1di|yi|2≥0,for ally∈Cn. Lettingy=ei, theny∗Dy=di. Hence, all of the eigenvalues of Aare nonnegative. ⇐LetA∈Mnbe a Hermitian matrix and suppose that all λi(A)≥0. Then n/summationdisplay i=1di|yi|2≥0 and hence, x∗Ax≥0 for allx∈Cn. Corollary 4.2: IfA∈Mnis positive semidefinite, then so are all the powers Ak,k=1,2,3,.... Proof: If the eigenvalues of Aareλ1,λ2, ...,λn, then the eigenvalues of Akare λ1k,λ2k, ...,λnk. 18 A positive semidefinite matrix can have more than one square root. However, it can only have one positive semidefinite matrix square root. The proof of the following result is adapted from [8]. Theorem 4.3: LetA∈Mnbe positive semidefinite and let k≥1 be a given integer. Then there exists a unique positive semidefinite Hermitian matrix Bsuch thatBk=A. We also have (a)BA=ABand there is a polynomial p(t) such that B=p(A); (b) rankB= rankA,s oBis a positive definite if Ais; and (c)Bis real ifAis real. Proof: We know that the Hermitian matrix A can be unitarily diagonalized as A=UDU∗withD= diag(λ1,λ2,....,λ n) and allλi≥0. We define B=UD1 kU∗, whereD1/k≡diag(λ1/k 1,λ1/k 2,...,λ1/k n), and the unique nonnegative kth root is taken in each case. Clearly, Bk=AandBis Hermitian and positive semidefinite. Also, AB=UDU∗UD1 kU∗=UDD1 kU∗=UD1 kDU∗=UD1 kU∗UDU∗=BA, andB is positive semidefinite because all λi(and hence their kthroots) are nonnegative. The rank of Bis just the number of nonzero λiterms, which is also the rank of A. IfAis real and positive semidefinite, then we know that Umay be chosen to be a real orthogonal matrix, so it is clear that Bcan be chosen to be real in this case. It remains only to consider the question of uniqueness. Notice first that there is a polynomial p(t) such that p(A)=B; we need only choosep(t) the Lagrange interpolating poynomial for the set {(λ1,λ1 k 1),..., (λn,λ1 kn)} to getp(D)=D1 kandp(A)=p(UDU∗)=Up(D)U∗=UD1 kU∗=B. But then ifCis any positive semidefinite Hermitian matrix such that Ck=A,w e haveB=p(A)=p(Ck) so thatCB =Cp(Ck)=p(Ck)C=BC. SinceB andCare commuting Hermitian matrices, they may be simultaneously unitarily 19 diagonalized; that is, there is some unitary matrix Vand diagonal matrices D1and D2with nonnegative diagonal entries such that B=VD 1V∗andC=VD 2V∗. Then from the fact that Bk=A=Ckwe deduce that Dk 1=Dk 2. But since the nonnegative kthroot of a nonnegative number is unique, we conclude that (Dk 1)1 k=D1=D2=(Dk 2)1 kandB=C. The most useful case of the preceding theorem is for k= 2. The unique positive (semi)definite square root of the positive (semi)definite matrix Ais usually denoted byA1 2. Similarly, A1 kdenotes the unique positive (semi)definite kthroot ofAfor eachk=1,2,.... Ann×nHermitian matrix Ais called positive definite ,i f x∗Ax > 0 for all nonzero x∈Cn. IfAis a real symmetric positive definite matrix, then Acan be factored into a productLDLT, whereLis a lower triangular matrix with 1’s along the diagonal andDis a diagonal matrix, whose diagonal entries are all positive. Then A= (LD1 2)(D1 2LT), which is the Cholesky Decomposition of A. It is worth noting that many applications can be recast to use a Cholesky Decomposition instead of the square root. The Cholesky Decomposition becomes referred to as a square root decomposition and has many applications such as in multiwavelet representations, predictive control, and square-root filters. 20 5. General Square Roots Our main tool in this chapter is the Jordan Canonical Form. Theorem 5.1 (Jordan Form Theorem): LetA∈Mn. Then there is a non- singular matrix Ssuch thatS−1AS=Jis a direct sum of Jordan blocks. Further- more,Jis unique up to a permutation of the Jordan blocks. LetA=SJS−1be the Jordan canonical form of the given matrix A, so that if X2=A=SJS−1, thenS−1X2S=(S−1XS)2=J. It suffices, therefore, to solve the equation X2=J. But ifXis such that the Jordan canonical form of X2is equal toJ, then there is some nonsingular Tsuch thatJ=TX2T−1=(TXT−1)2. Thus, it suffices to find an Xsuch that the Jordan canonical form of X2is equal toJ. If the Jordan canonical form of Xitself isJX, then the Jordan canonical form of X2is the same as that of ( JX)2, so it suffices to find a Jordan matrix JXsuch that the Jordan canonical form of ( JX)2is equal toJ. Finally, if JX=Jm1(µ1)⊕···⊕Jmr(µr), then the Jordan canonical form of ( JX)2is the same as the direct sum of the Jordan canonical forms of Jmi(µi)2,i=1,...,r . Thus, to solve X2=A, it suffices to consider only whether there are choices of scalars µand positive intergers msuch that the given Jordan canonical form Jis the direct sum of the Jordan canonical forms of matrices of the form Jm(λ)2.I fµ/negationslash= 0, we know that the Jordan canonical form ofJm(µ)2isJm(µ2), so every nonsingular Jordan block Jk(λ) has a square root; in fact, it has square roots that lie in two distinct similarity classes with Jordan canonical forms Jk(±√ λ). If necessary, these square roots of Jordan blocks can be computed explicitly. Thus, every nonsingular matrix A∈Mnhas a square root, and it has square roots in at least 2µdifferent similarity classes, if Ahasµdistinct eigenvalues. It has square roots in at most 2νdifferent similarity classes, if the Jordan canonical form 21 ofAis the direct sum of νJordan blocks; there are exactly 2νsimilarity classes, if all the Jordan blocks with the same eigenvalue have different sizes, but there are fewer than 2ν, if two or more blocks of the same size have the same eigenvalue, since permuting blocks does not change the similarity class of a Jordan canonical form. Some of these nonsimilar square roots may not be “functions” of A, however. If each of theνblocks has a different eigenvalue (i.e., Ais nonderogatory), a Lagrange- Hermite interpolation polynomial can always be used to express any square root ofAas a polynomial in A. If the same eigenvalue λoccurs in two or more blocks, however, polynomial interpolation is possible only if the same choice is made for λ1 2for all of them; if different choices are made in this case, one obtains a square root ofAthat is not a “function” of Ain the sense that it cannot be obtained as a polynomial in A, and therefore cannot be a primary matrix function f(A) with a single-valued function f(·). What happens if Ais singular? Since each nonsingular Jordan block of Ahas a square root, it suffices to consider the direct sum of all the singular Jordan blocks ofA.I fAhas a square root, then this direct sum is the Jordan canonical form of the square of a direct sum of singular Jordan blocks. Which direct sums of singular Jordan blocks can arise in this way? Letk>1. We know that the Jordan canonical form of Jk(0)2consists of exactly two Jordan blocks Jk/2(0)⊕Jk/2(0) ifk>1 is even, and it consists of exactly two Jordan blocks J(k+1)/2(0)⊕J(k−1)/2(0) ifk>1 is odd. Ifk=1,J1(0)2= [0] is a 1-by-1 block, and this is the only Jordan block that is similar to the square of a singular Jordan block. Putting together this information, we can determine whether or not a given singular Jordan matrix Jhas a square root as follows: Arrange the diagonal blocks 22 ofJby decreasing size, so J=Jk1(0)⊕Jk2(0)⊕···⊕Jkp(0) withk1≥k2≥ k3≥ ··· ≥kp≥1.Consider the differences in sizes of successive pairs of blocks: ∆1=k1−k2,∆3=k3−k4,∆5=k5−k6,etc., and suppose Jis the Jordan canonical form of the square of a singular Jordan matrix ˜J. We have seen that ∆ 1=0o r1 because either k1= 1 [in which case J1(0)⊕J1(0) corresponds to ( J1(0)⊕J1(0))2 or toJ2(0)2]o rk1>1 andJk1(0)⊕Jk2(0) corresponds to the square of the largest Jordan block in ˜J, which has size k1+k2. The same reasoning shows that ∆ 3,∆5,... must all have the value 0 or 1 and an acceptable square root corresponding to the pairJki(0)⊕Jki+1(0) isJki+ki+1(0),i=1,3,5,....I fp(the total number of singular Jordan blocks in J) is odd, then the last block Jkp(0) is left unpaired in this process and must therefore have size 1 since it must be the square of a singular Jordan block. Conversely, if the successive differences (and kp,i fpis odd) satisfy these conditions, then the pairing process described constructs a square root for J. SupposeA∈Mnis singular and suppose there is a polynomial r(t) such that B=r(A) is a square root of A. Thenr(0) = 0,r(t)=tg(t) for some polynomial g(t), andA=B2=A2g(A)2, which is clearly impossible if rank A2<rankA. Thus, rank A= rankA2in this case, which means that every singular Jordan block ofAis 1-by-1. Conversely, if Ais singular and has minimal polynomial qA(t)=t(t−λ1)r1···(t−λµ)rµwith distinct nonzero λ1,...,λ µand allri≥1, let g(t) be a polynomial that interpolates the function f(t)=1/√ tand its derivatives at the (necessarily nonzero) roots of the polynomial qA(t)/t= 0, and let r(t)≡tg(t). For each nonsingular Jordan block Jni(λi)o fA,g(Jni(λi)) = [Jni(λi)]−1 2and hence, we haver(Jni(λi)) =Jni(λi)[Jni(λi)]−1 2=Jni(λi)1 2. Since all the singular Jordan blocks ofAare 1-by-1 and r(0) = 0, we conclude that r(A) is a square root of A. Thus, a given singular A∈Mnhas a square root that is a polynomial in Aif and only if rank A= rankA2. Since this latter condition is trivially satisfied if A is nonsingular (in which case we already know that Ahas a square root that is a 23 polynomial in A), we conclude that a given A∈Mnhas a square root that is a polynomial in Aif and only if rank A= rankA2. If we agree that a “square root” of a matrix A∈Mnis any matrix B∈Mn such thatB2=A, we can summarize what we have learned about the solutions of the equation X2−A= 0 in Theorem 5.2. Theorem 5.2: LetA∈Mnbe given. (a) IfAis nonsingular and has µdistinct eigenvalues and νJordan blocks in its Jordan canonical form, it has at least 2µand at most 2νnonsimilar square roots. Furthermore, at least one of its square roots can be expressed as a polynomial in A. (b) IfAis singular and has Jordan canonical form A=SJS−1, letJk1(0)⊕ Jk2(0)⊕···⊕Jkp(0) be the singular part of Jwith the blocks arranged in decreasing order of size: k1≥k2≥···≥kp≥1. Define ∆ 1=k1−k2,∆3=k3−k4,.... Then Ahas a square root if and only if ∆ i= 0 or 1 for i=1,3,5,...and, ifpis odd, kp=1. Furthermore, Ahas a square root that is a polynomial in Aif and only if k1= 1, a condition that is equivalent to requiring that rank A= rankA2. (c) IfAhas a square root, its set of square roots lies in finitely many different similarity classes. Since the sizes and numbers of the Jordan blocks Jk(λ) of a matrix Acan be inferred from the sequence of ranks of the powers ( A−λI)k,k=1,2,..., the necessary and sufficient condition on the sizes of the singular Jordan blocks of A in part (b) of the preceding theorem can be restated in terms of ranks of powers. LetA∈Mnbe a given singular matrix, and let r0=n, r k=rank Akfork= 1,2,.... The sequence r0,r1,r2,...is decreasing and eventually becomes constant. Ifrk1−1>r k1=rk1+1=..., then the largest singular Jordan block in Ahas sizek1, which is the index of the matrix with respect to the eigenvalue λ= 0. The difference 24 rk1−1−rk1is the number of singular Jordan blocks of size k1. If this number is even, the blocks of size k1can all be paired together in forming a square root. If this number is odd, then one block is left over after the blocks are paired and A can have a square root only if either k1= 1 (so that no further pairing is required), or there is at least one singular Jordan block of size k1−1 available to be paired with it; this is the case only if rk1−2−rk1−1>r k1−1−rk1, sincerk1−2−rk1−1equals the total number of singular Jordan blocks of sizes k1andk1−1. This reasoning is easily continued backward through the sequence of ranks rk. If all the differences ri−ri+1are even,i=k1−1,k1−3,..., thenAhas a square root. If any difference ri−ri+1is odd, then ri−1−rimust have a larger value, if Ais to have a square root. Since r0−r1is the total number of singular Jordan blocks of all sizes, if r0−r1is odd, we must also require that there be at least one block of size 1, that is, 1≤(# of singular blocks of all sizes ≥1)−(# of singular blocks of all sizes ≥ 2 )=(r0−r1)−(r1−r2)=r0−2r1+r2.Notice that rk≡n, ifAis nonsingular, so all the successive differences ri−ri+1are zero and Atrivially satisfies the criteria for a square root in this case. This theorem largely results from [3] but is presented more clearly in [9].Corollary 5.3: LetA∈Mnand letr0=n,r k= rankAkfork=1,2,.... Then Ahas a square root if and only if the sequence {rk−rk+1},k=0,1,2,... does not contain two successive occurrences of the same odd integer and, if r0−r1 is odd,r0−2r1+r2≥1. As mentioned in chapter 3, we can use direct calculation to show that there is no matrix B=/bracketleftbiggab cd/bracketrightbigg 25 such that B2=A=/bracketleftbigg01 00/bracketrightbigg . We can also use the criteria in Theorem 5.2 and Corollary 5.3 given by [9] to show that no such matrix Bexist. Although we now know exactly when a given complex matrix has a complex square root, one sometimes needs to answer a slightly different question: When does a given real matrix A∈Mn(R) have a real square root? The equivalent criteria in Theorems 5.2 and Corollary 5.3 are still necessary, of course, but they do not guarantee that any of the possible square roots are real. The crucial observation needed here is that if one looks at the Jordan canonical form of a real matrix, then the Jordan blocks with nonreal eigenvalues occur only in conjugate pairs, i.e., there is an even number of Jordan blocks of each size for each nonreal eigenvalue. Moreover, a given complex matrix is similar to a real matrix if and only if the nonreal blocks in its Jordan canonical form occur in conjugate pairs. If there is some B∈Mnsuch thatB2=A, then any Jordan block of Awith a negative eigenvalue corresponds to a Jordan block of Bof the same size with a purely imaginary eigenvalue. If Bis real, such blocks must occur in conjugate pairs, which means that the Jordan blocks of Awith negative eigenvalues must also occur in pairs, just like the nonreal blocks of A. Conversely, let Jbe the Jordan canonical form of the real matrix A∈Mn(R), suppose all of the Jordan blocks in Jwith negative eigenvalues occur in pairs, and suppose Asatisfies the rank conditions in Corollary 5.3. Form a square root forJusing the process leading to Theorem 5.2 for the singular blocks, and using the primary-function method (found in [9]) for each individual nonsingular Jordan block, but be careful to choose conjugate values for the square root for the two members of each pair of blocks with nonreal or negative eigenvalues; blocks or 26 groups of blocks with nonnegative eigenvalues necessarily have real square roots. Denote the resulting, possibly complex, block diagonal upper triangular matrix by C,s oC2=J. Each diagonal block of Cis similar to a Jordan block of the same size with the same eigenvalue, so Cis similar to a real matrix Rbecause of the conjugate pairing of its nonreal Jordan blocks. Thus, the real matrix R2is similar toC2=J, andJis similar to the real matrix A,s oR2is similar to A. Recall that two real matrices are similar if and only if they are similar via a real similarity, since they must have the same real Jordan canonical form, which can always be attained via a real similarity. Thus, there is a real nonsingular S∈Mn(R) such thatA=SR2S−1=SRS−1SRS−1=(SRS−1)2and the real matrix SRS−1is therefore a real square root of A. In the above argument, notice that if Ahas any pairs of negative eigenvalues, the necessity of choosing conjugate purely imaginary values for the square roots of the two members of each pair precludes any possibility that a real square root of Acould be a polynomial in Aor a primary matrix function of A. The following theorem summarizes these observations. Theorem 5.4: LetA∈Mn(R) be a given real matrix. There exists a real B∈Mn(R) withB2=Aif and only if Asatisfies the rank condition given in Corollary 5.3 and has an even number of Jordan blocks of each size for every negative eigenvalue. If Ahas any negative eigenvalues, no real square root of Acan be a polynomial in Aor a primary matrix function of A. The same reasoning used before to analyze the equation X2=Acan be used to analyzeXm=Aform=3,4,.... Every nonsingular A∈Mnhas anmth root, in fact, a great many of them, and the existence of an mth root of a singular matrix is determined entirely by the sequence of sizes of its singular Jordan blocks. 27 From [11], it is important to note that if σ(A)∩(−∞,0] =∅, thenAhas a unique square root B∈Mnwithσ(B) in the open right (complex) half plane. 28 6. Computational Method To evaluate matrix square root functions, it is suggested that the most stable way is to use the Schur decomposition. Recalling Schur’s Theorem, we know that for each complex matrix Athere exist a unitary matrix qand upper triangular matrixt, such that A=qtq−1. A square root bof the upper triangular factor t could be computed by directly solving the equation b2=t.The choices of signs on the diagonal of b, b mm=√ tmmdetermine which square root is obtained. {tmm}are eigenvalues of A,{bmm}are eigenvalues of b, and the principal√ Ahas nonnegative eigenvalues or eigenvalues with a nonnegative real part. Using Schur decomposition, we have A=qtq−1, b2=t, c=qbq−1, c2=qbq−1qbq−1=qb2q−1=qtq−1=A. The most time consuming if done by hand is Schur decomposition. Fortunately, MatLab does it for us. Computing bthat satisfies b2=tis the next time consuming (if done by hand) portion of this process. Consider the following 4 ×4 matrix.  b11b12b13b14 0b22b23b24 00b33b34 000 b44  b11b12b13b14 0b22b23b24 00b33b34 000 b44 = t11t12t13t14 0t22t23t24 00t33t34 000 t44 . Since the eigenvalues of the upper triangular matrix tare the main diagonal entries, we can use them to calculate the eigenvalues of the matrix b. First we 29 compute entries of the main diagonal: bmm=√ tmm, then we compute entries on the first diagonal parellel to the main one. t12=b11b12+b12b22=b12(b11+b22) b12=t12 b11+b22 t23=b22b23+b23b33=b23(b22+b33) b23=t23 b22+b23 t34=b33b34+b34b44=b34(b33+b44) b34=t34 b33+b44 After that, we compute elements of the second diagonal parellel to the main one. t13=b11b13+b12b23+b13b33 t13−b12b23=b13(b11+b33) b13=t13−b12b23 b11+b33 t24=b22b24+b23b34+b24b44 t24−b23b34=b24(b22+b44) b24=t24−b23b34 b22+b44 Finally, we compute elements of the third diagonal parallel to the main one. In case of the 4-th order, it consists of only one element. t14=b11b14+b12b24+b13(b34+b14b44) 30 t14−b12b24−b13b34=b14(b11+b44) b14=t14−b12b24−b13b34 b11+b44. Therefore, we have b12=t12 b11+b22, b23=t23 b22+b23, b34=t34 b33+b44, b13=t13−b12b23 b11+b33, b24=t24−b23b34 b22+b44, b14=t14−b12b24−b13b34 b11+b44. Now we can derive the following algorithm for a matrix of n-th order. Form=1,2,...,n : bmm=√ tmm Form=1,2,...,n −1: bm,m +1=tm,m +1 bm,m+bm+1,m+1 31 Forr=2,3,...,n −1 andm=1,2,...,n −r: b(m,m +r)=t(m,m +r)−m+r+1/summationdisplay k=m+1b(m,k)·b(k,m +r) bm,m+bm+r,m+r Below is a script file (as it should be typed for use in MatLab) to find the square root of a matrix using Schur decomposition: n=input (‘Enter size of a matrix: ’) A=input (‘Enter n x n matrix:’) a=A+0.000001i*norm(A)*eye(n,n); eigvala=eig(a) [q,t]=schur(a); b=zeros(n,n); for m = 1:n b(m,m)= sqrt(t(m,m)); end; for m=1:n-1 b(m,m+1)=t(m,m+1)/(b(m,m)+b(m+1,m+1)); end; for r=2:n-1 for m=1:n-r B=0; for k=(m+1):(m+r-1) B=B+b(m,k)*b(k,m+r); end; b(m,m+r)=(t(m,m+r)-B)/(b(m,m)+b(m+r,m+r)); end; 32 end; b c=q*b*q’ eigvalc=eig(c) csquare=c*c A er=norm(A-c*c)/norm(A) Explanation of the Program n=input (‘Enter size of a matrix: ’): Indicate the size of the matrixA=input (‘Enter n x n matrix:’): Enter the entries of the matrix enclosed in brackets. Separate each entry with a comma and each row with a semicolon. a=A+0.000001i*norm(A)*eye(n,n);: If a matrix Ais real with complex con- jugate eigenvalues, MatLab will automatically return a real matrix with Jordan blocks instead of an upper triangular complex matrix unless we indicate that weare interested in the complex triangular matrix. It is more difficult to compute the square root of a real matrix with Jordan blocks than it is of a triangular one. The term ‘eye’ in our command refers to the identity matrix, I. The addition, A+/epsilon1I=qtq /prime+/epsilon1qIq/prime=q(t+/epsilon1I)q/prime,will change the matrix tby/epsilon1I, which is very small for small number /epsilon1. Our computation will show that our results has error less than/epsilon1. eigvala=eig(a): This command yields the eigenvalues of the matrix A. [q,t]=schur(a);: This command shows the breakdown of the matrix Ausing Schur decomposition. Here qis unitary and tis upper triangular. 33 b=zeros(n,n);: This command sets all entries of the triangular matrix bto zero to begin the cycle. for m = 1:n b(m,m)=sqrt(t(m,m)); end;: This cycle yields the main diagonal entries for the matrix b. for m=1:n-1 b(m,m+1)=t(m,m+1)/(b(m,m)+b(m+1,m+1)); end;: This cycle yields the other diagonal entries {b12,b23,b34}for the matrix b. The following section of the program yields the general formula for finding the matrix entries of b: bm,m +r=tm,m +r−m+r+1/summationdisplay k=m+1bm,k·bk,m+r bm,m+bm+r,m+r,r≥2. This formula has two nonnegative terms (or terms with nonnegative real parts) in the denominator. They are square roots of eigenvalues of the given matrix. Even if the given matrix has 2 zero eigenvalues, the denominator will not be zero because of the added matrix 0.000001i*norm(A)*eye(n,n). But in this case the result is not reliable, it may have a larger error. for r=2:n-1 for m=1:n-r B=0; for k=(m+1):(m+r-1) B=B+b(m,k)*b(k,m+r); end; b(m,m+r)=(t(m,m+r)-B)/(b(m,m)+b(m+r,m+r)); end; 34 end; b: Prints the matrix b c=q*b*q’: Gives the final result c=√ A. Sinceqis a unitary matrix, we can use transpose, q/prime, instead of the inverse, q−1in our program. eigvalc=eig(c): This command prints the eigenvalues of the matrix c. csquare=c*c: This command calculates c2. A: Prints the matrix Afor comparison with matrix c2. er=norm(A-c*c)/norm(A): This commands computes the relative error. Here is an example of computation by the program above: /greatermuchSchurSqrt(name given to program) Enter the size of a matrix: 4 Output: n = 4 Enter the matrix n×n: [2, 3, 1, 5; 0, 2, 5, -3; -1, 2, 3, 0; 2, 4, -2, 1] Output: A= 23 1 5 02 5 −3 −12 3 024−21  Output: eigvala= 5.1193 + 2.3938i 5.1193−2.3936i −1.1193 + 2.3938i −1.1193−2.3936i 35 Output: q= −0.3131 + 0.6646i−0.3443−0.3407i0.3651 + 0.0808i−0.1874−0.2253i −0.4184−0.1222i−0.3389 + 0.5834i0.0488−0.3534i0.2551−0.4030i −0.3214−0.0658i−0.2637 + 0.4194i−0.0033 + 0.4943i−0.3791 + 0.5088i −0.1982 + 0.3512i−0.2096−0.1442i−0.6745−0.1832i0.3723 + 0.3815i  Output: t= 5.1193 + 2.3938i2.3321−0.4383i−2.1173 + 1.5117i−0.3298−2.4575i 05.1193−2.3936i2.6075−1.9311i0.8210 + 2.3928i 00 −1.1193 + 2.3938i0.1782 + 1.6295i 00 0 −1.1193−2.3936i  Output: b= 2.3206 + 0.5158i0.5025−0.0944i−0.2775 + 0.7764i0.3460−0.6983i 02.3206−0.5157i0.6106−0.7684i−0.2512 + 0.4469i 00 0 .8727 + 1.3715i0.1021 + 0.9336i 00 0 0 .8727−1.3714i  Output: c= 1.2751 + 0.0000i0.0942 + 0.0000i0.9015−0.0000i1.7204−0.0000i 0.2889−0.0000i1.7738 + 0.0000i0.9570 + 0.0000i−1.1722 + 0.0000i −0.4103 + 0.0000i0.4295 + 0.0000i1.8323 + 0.0000i0.3623−0.0000i 0.4166 + 0.0000i1.3518−0.0000i−1.0993 + 0.0000i1.5054 + 0.0000i  Output: eigvalc= 0.8727 + 1.3715i 0.8727−1.3714i 2.3206 + 0.5158i 2.3206−0.5157i 36 Output: (csquare) c2= 23 1 5 02 5 −3 −12 3 0 24−21  Output: A= 23 1 5 02 5 −3 −12 3 024−21  Output: er = 1 .0000e−006 Once again, to find a matrix c 2=Ausing Schur decomposition, we begin with A=qtq−1, whereqis a unitary matrix and tis an upper triagular matrix. From here we have, b2=t, b=√ t, where√ tis the principal square root and c=qbq−1. 37 References [1] S. Athloen and R. McLaughlin, Gauss-Jordan Reduction: A Brief History , American Mathematical Monthly 94, 130-142, 1987. [2] A. Choudhry, Extraction of nth Roots of 2×2Matrices , Linear Algebra and its Applications, 387:183-192, 2004. [3] G.W. Cross, P. Lancaster, Square roots of Complex Matrices, Linear and Mul- tilinear Algebra , 289-293, 1974. [4] L.A. Hageman, D. M. Young, Applied Iterative Methods , Academic Press, New York, 1981. [5] N.J. Higham, Computing Real Square Roots of a Real Matrix , Linear Algebra and its Applications, 88/89:405-430, 1987. [6] N.J. Higham, Functions of Matrices , Chapter 11 in Handbook of Linear Alge- bra, Chapman/CRC Press, Boca Raton, 2007. [7] N.J. Higham, Stable Iterations for the Matrix Square Root , Numerical Algo- rithms 15, 227-242, 1997. [8] R.A. Horn, C.R. Johnson, Matrix Analysis , Cambridge University Press, 1985. [9] R.A. Horn, C.R. Johnson, Topics in Matrix Analysis , Cambridge University Press, 1991. [10] T.R. Hughes, I.Levit, and J.Winget, Element-by-Element Implicit Algorithms for Heat Conduction , J. Eng. Mech., 109(2), 576-585, 1983. [11] C.R. Johnson, K. Okubo, and R. Reams, Uniqueness of Matrix Square Roots and an Application , Linear Algebra and its Applications, 323:51-60, 2001. 38 [12] P. Lancaster, M. Tismenetsky, The Theory of Matrices Second Edition with Applications , Academic Press, 1985. [13] D.C. Lay, Linear Algebra and Its Applications-3rd edition , Pearson Education, Inc., 2003. [14] S.J. Leon, Linear Algebra with Applications-6th edition , Prentice-Hall,Inc., 2002. [15] B.W. Levinger, The Square Root of a 2×2Matrix , Mathematics Magazine, Vol 53, No.4, 1980. [16] P.N. Parlett, The Symmetric Eigenvalue Problem , Prentice-Hall, Englewood Cliffs, NJ, 1980. [17] B.A. Schmitt, An Algebraic Approximation for the Matrix Exponential in Sin- gularly Perturbed Boundary Value Problems , SIAM J. Numer. Anal. 27(1), 51-66, 1990. [18] D. Sullivan, The Square Root of 2×2Matrices , Mathematics Magazine, Vol 66, No. 5, 1993. [19] A. Tucker, The Growing Importance of Linear Algebra in Undergraduate Math- ematics , The College Mathematics Journal, 24, 3-9, 1993.