chapter7

Chapter 7 Factorization Theorems This chapter highlights a few of the many factorization theorems for ma- trices. While some factorization resul ts are relatively direct, others are it- erative. While some factorization results serve to simplify the solution tolinear systems, others are concerned with revealing the matrix eigenvalues.We consider both types of results here. 7.1 The PLU Decomposition The PLU decomposition (or factorization) To achieve LU factorization werequire a modi fied notion of the row reduced echelon form. Definition 7.1.1. The modified row echelon form of a matrix is that form which satis fies all the conditions of the modi fied row reduced echelon form except that we do not require zeros to be above leading ones, and moreoverwe do not require leading ones, just nonzero entries. For example the matrices below are in row echelon form. A= 123 001000 B= 1 230 04−76 0 001  Most of the factorizations A∈M n(C) studied so far require one essential ingredient, namely the eigenvectors of A. While it was not emphasized when we studied Gaussian elimination, there is a LU-type factorization there. Assume for the moment that the only operations needed to carry Ato its 201 202 CHAPTER 7. FACTORIZATION THEOREMS modi fied row echelon form are those that add a multiple of one row to another. The modified row echelon form of a matrix is that form which satisfies all the conditions of the modi fied row reduced echelon form except that we do not require zeros to be above leading ones, and moreover wedo not require leading ones, just nonzero entries. Naturally it is easy to make the leading nonzero entries into leading ones by the multiplication by an appropriate identity matrix. That is not the point here. What wewant to observe is that in this case the reduction is accomplished by the leftmultiplication of Aby a sequence of lower triangular matrices of the form. L= 1 01 0 ...01 c... 0··· 1  Since we pivot at the (1 ,1)-entry first, we eliminate all the entries in the first column below the first row. The product of all the matrices Lto accomplish this has the form L 1= 1 c 2110 c3101 ...... cn10··· 1  where c k1=−ak1 a11.Thus, with the notation that A=A1has entries a(1) ijthis first phase of the reduction renders the matrix A2with entries a(2) ij A2=L1A1= a (2) 11 ··· a(2) 1n 0a(2) 22 ···... 0a(2) 32a(2)33 ......... 0a(2) n2··· a(2) nn  Since we have assumed that no row interchanges are necessary to carry out the reduction we know that a (2) 226=0.The next part of the reduction process is the elimination of the elements in the second column below the second 7.1. THE PLU DECOMPOSITION 203 row, i.e. a(2) 32→0, ...a(2) n2→0.Correspondingly, this can be achieved by a matrix of the form L2= 1 01 0 0c 22 1 ......... 0cn2··· 1  (What are the values c k2?) The result is the matrix A3given by A3=L2A2=L2L1A1= a (3) 11 ··· a(3) 1n 0a(3) 22 ···... 00 a(3) 33............... 00 a (3) 3n a(3) nn  Proceeding in this way through all the rows (columns) there results A n=Ln−1An−1=Ln−1···L2L1A1= a (3) 11 ··· a(3) 1n 0a(3) 22 ···... 00 a(3) 33............... 000 a(3) nn  The right side of the equation above is an upper triangular matrix. Denote it by U.Since each of the matrices L i,i=1,...n−1i s i n v e r t i b l e w e c a n write A=L−1 1···L−1 n−1U The lemma below is useful in this. Lemma 7.1.1. Suppose the lower triangular matrix L∈Mn(C)has the 204 CHAPTER 7. FACTORIZATION THEOREMS form L= 1 0... 0 1 01 ......c k+1,k... ...... 0··· 0cnk 1 ←−k throw Then Lis invertible with inverse given by L−1= 1 0... 0 101 ......−c k+1,k... ...... 0··· 0−cnk 1 ←−k throw Proof. Trivial Lemma 7.1.2. Suppose L1,L2,···,Ln−1are the matrices given above. Then the matrix L=L−1 1···L−1 n−1has the form L= 1 −c 21 10 −c31−c32 1 1 .........−ck+1,k... ...... −cn1−cn2···−cnk ··· 1  Proof. Trivial. Applying these lemmas to the present situation we can say that when no row interchanges are needed we can factor and matrix A∈M n(C)a s A=LU,where Lis lower triangular and Uis upper triangular. When row 7.1. THE PLU DECOMPOSITION 205 interchanges are needed and we let Pbe the permutation matrix that creates these row interchanges then the LU-factorization above can be carried outfor the matrix PA. Thus PA=LU, where Lis lower triangular and Uis upper triangular. We call this the PLU factorization. Let us summarize this in the following theorem. Theorem 7.1.1. LetA∈M n(C). Then there is a permutation matrix P∈Mn(C)and lower Land upper Utriangular matrices ( ∈Mn(C)), such thatPA=LU. Moreover, Lcan be taken to have ones on its diagonal. That is,`ii=1,i=1,...n . By applying the result above to ATit is easy to see that the matrix U can be taken to have the ones in its diagonal. The result is stated as a corollary. Corollary 7.1.1. LetA∈Mn(C). Then there is a permutation matrix P∈Mn(C)and lower and upper triangular matrices ( ∈Mn(C)) respec- tively, such that PA=LU. Moreover, Ucan be taken to have ones on its diagonal ( uii=1,i=1,...n ). The PLU decomposition can be put in service to solving the system Ax=bas follows. Assume that A∈Mn(C) is invertible. Determine the permutation matrix Pin order that PA=LU, where Lis lower triangular andUis upper triangular. Thus, we have Ax =b PAx =Pb LUx =Pb Solve the systems Ly =Pb Ux =y Then LUx =Ly=Pb.Hence xis a solution to the system. The advantages of this formulation over the direct Gaussian elimination is that the systemsLy=PbandUx=yare triangular and hence are easy to solve. For example for the first of the systems, Ly=Pb,let the vector Pb=h ˆb 1,..., ˆbniT . Then it is easy to see that “back substitution” (aka “forward substitution”) 206 CHAPTER 7. FACTORIZATION THEOREMS can be used to determine y. That is, we have the recursive relations y1=ˆb1 l11 y2=ˆb2−l21y1 l22 ... yn=à ˆbn−n−1X m=1lnmym! l−1 nn A similar formula applies to solve Ux=y. I nt h i sc a s ew es o l v e first for xn=yn/unn.The general formula is recursive with xkbeing determined after xk+1,...,x n.are determined using the formula xk=à yk−nX m=k+1ukmym! u−1 kk In practice the step of determining and then multiplying by the per- mutation matrix is not actually carried out. Rather, an index array is generated, while the elimination step is accomplished that e ffectively inter- changes a “pointer” to the row interchanges. This saves considerable timein solving potentially very large systems. More general and instructive methods are available for accomplishing this LU factorization. Also, conditions are available for when no (nontrivial)permutation is required. We need the following lemma. Lemma 7.1.3. LetA∈M n(C)have the LU factorization A=LU,w h e r e Lis lower triangular and Uis upper triangular. For any partition of the matrix of the form A=·A11A12 A21A22¸ there are corresponding decompositions of the matrices LandU L=·L11 0 L21L22¸ and U=·U11U12 0U22¸ 7.1. THE PLU DECOMPOSITION 207 where the Liiand the Uii.are lower and upper triangular respectively. More- over, we have A11=L11U11 A21=L21U11 A12=L12U22 A22=L21U12+L22U22 Thus L11U11is a LU factorization of A11. With this lemma we can establish that almost every matrix can have a LU factorization. Definition 7.1.2. LetA∈Mn(C) and suppose that 1 ≤j≤n.T h e expression det( A{1,...,j }) means the determininant of the upper left j×j submatrix of A. These quaditities for j=1,...,n are called the principal determinants of A. Theorem 7.1.2. LetA∈Mn(C)and suppose that Ahas rank k.If det(A{1,...,j })6=0 forj=1,...,k (1) thenAhas a LU factorization A=LU,w h e r e Lis lower triangular and U is upper triangular. Moreover, the factorization may be taken so that either LorUis nonsingular. In the case k=nbothLandUwill be nonsingular. Proof. We carry out this LU factorization as a direct calculation in compar- ison to the Gaussian elimination method above. Let us propose to solve the equation LU=Aexpressed as  l 11 l21l22 0 l31l32l33 ............ ... ln1ln2··· ··· lnn  u 11u12u13··· u1n u22u23··· u2n u33 0...... ... unn  = a 11a12a13 a1n a21a22a23 a2n a31a32a33 ............... ... an1an2··· ··· ann  208 CHAPTER 7. FACTORIZATION THEOREMS It is easy to see that l11u11=a11.We can take, for example l11=1a n d solve for u11.The detminant condition assures us that u116=0.Next solve for the (2 ,1)-entry. We have l21u11=a21.Since u116=0,solve for l21. For the (1 ,2)-entry we have l11u12=a12,w h i c hc a nb es o l v e df o r u12since l116= 0. Finally, for the (2 ,2)-entry, l12u12+l22u22=a22is an equation with two unknowns. Assign l22= 1 and solve for u22.What is important to note is that the process carried out this way gives the factorization of theupper left 2 ×2 submatrix of A.Thus ·l 110 l21l22¸·u11u12 0u22¸ =·a11a12 a21a22¸ Since detµ·a11a12 a21a22¸¶ 6=0,it follows that detµ·u11u12 0u22¸¶ 6=0a n d we know that·l110 l21l22¸ is nonsingular as the diagonal elements are ones. Continue the factorization process through the k×kupper left submatrix ofA. Now consider the blocked matrix form form A A=·A11A12 A21A22¸ where A11isk×kand has rank k. Thus we know that the rows of the lower (n−k)×nmatrix above, that is£ A21A22¤ c a nb ew r i t t e na sau n i q u e linear combination of the rows of the upper k×nmatrix£ A11A12¤ .Thus £ A21A22¤ =C£ A11A12¤ for some ( n−k)×kmatrix C.Of course this means: A21=CA 11and A22=CA 12. We consider the factorization A=·A11A12 A21A22¸ =·L11 0 L21L22¸·U11U12 0U22¸ where the blocks L11andU11have just been determined. From the equations in the lemma above we solve to get U12=L−1 11A12andL21= 7.2. LRLRLRFACTORIZATION 209 A12U−1 11.T h e n A22=L21U12+L22U22 =A12U−1 11L−1 11A12+L22U22 =A12A−1 11A12+L22U22 =CA 11A−1 11A12+L22U22 =CA 12+L22U22 =A22+L22U22 Thus we solve L22U22=0.Obviously, we can take for L22any nonsingular matrix we wish and solve for U22or conversely. 7.2LRLRLRfactorization While the PLU factorization is useful for solving systems, the LR factoriza- tion can be used to determine eigenvalues. . LetA∈Mnbe given. Then A=A1=L1R1. Then L−1 1A1L1=R1L1≡A2 A2=L2R2 L−1 2A2L2=R2L2≡A3. Continue in this fashion to obtain L−1 kAkLk=RkLk≡Ak+1 (?) We de fine Pk=L1L2...L k Qk=Rk...R 2R1. Then PkAk+1=A1Pk 210 CHAPTER 7. FACTORIZATION THEOREMS for Ak+1=L−1 kAkLk =L−1 kL−1 k−1Ak−1Lk−1Lk ... =P−1 kA1Pk or PkAk+1=A1Pk. Hence PkQk=Pk−1AkQk−1 =A1Pk−1Qk−1 =A1Pk−2Ak−1Qk−2 =A2 1Pk−2Qk−2 ... =Ak 1. Theorem 7.2.1 (Rutishauser). LetA∈Mnbe given. Assume the eigen- values of Asatisfy |λ1|>|λ2|>···>|λn|>0. Then A∼Λ=diag(λ1...λn). Assume A=SΛS−1,a n d Y≡S−1=LyRy X=S=LxRx where LyandLxare lower unit triangular matrices and RyandRxare upper triangular. Then Akdefined by (?)satisfy the result limAkis upper triangular. Proof. (Wilkinson) We have Ak 1=XΛkY =XΛkLyRy =XΛkLyΛ−kΛkRy. 7.3. THE QRALGORITHM 211 By the strict inequalities between the eigenvalues we have (ΛkLyΛ−k)ij=  1 i=j µλi λj¶k `iji>j 0 i<j . HenceΛkLyΛ−k→I(because|λi| |λj|<1i fi>j ). Hence with Ak 1=LxRx(ΛkLyΛ−k)ΛkRy and Ak 1=PkQk we conclude that lim k→∞Pk=Lx. Therefore Lk=P−1 k−1Pk→I. Finally we have that Akmust be upper triangular because L−1 kAk=Rk is upper triangular. This exposes all the eigenvalues of A.T h e r e f o r e t h e e i g e n v e c t o r s c a n b e determined. 7.3 The QRQRQRalgorithm Certain numerical problems with the LUalgorithm have led to the QR algorithm, which is based on the decomposition of the matrix Aas A=QR where Qis unitary and Ris upper triangular. Theorem 7.3.1 (QR-factorization). (i) Suppose Ais inMn,mandn≥ m. Then there is a matrix Q∈Mn,mwith orthogonal columns and an u p p e rt r i a n g u l a rm a t r i x R∈Mmsuch that A=QR. 212 CHAPTER 7. FACTORIZATION THEOREMS (ii) If n=m,t h e n Qis unitary. If Ais nonsingular the diagonal entries ofRcan be chosen to be positive. (iii) If Ais real; then QandRm a yb ec h o s e nt ob er e a l . Proof. (i) We proceed inductively. Let a1,... , a ndenote the columns ofAandq1,q2,... ,q mdenote the columns of Q. The basic idea of the QR-factorization is to orthogonalize the columns of Afrom left to right. Then the columns can be expressed by the formulas ak=Pk i=1ckqk,k =1,...,n .T h e c o e fficients of the expansion become, respectively, the entries of the kthcolumn of R,c o m p l e t e db y n−k zeros. (Of course, if the rank of Ais less than m,w e fill in arbitrary orthogonal vectors which we know exist as m≤n.) For the details, first de fineq1=a1/ka1k. To compute q2we use the Gram—Schmidt procedure. ˆq2=a2−hq1,a1iq1 q2=ˆq2/kˆq2k. Tracing backwards note that a2=ˆq2+hq1,a1iq1 =kˆq2kq2+hq1,a1iq1. So we have ·a1a2a3 ↓↓↓ ...¸ =·q1q2q3 ↓↓↓ ...¸ ka 1khq1,a1i... 0 kˆq2k ...0 00 . Instead of the full inductive step we compute q 3andfinish at that point ˆq3=a3−hq1,a3iq1−hq2,a3iq2 q3=ˆq3/kˆq3k. Hence a3=kˆq3kq3+hq1,a3iq1+hq2,a3iq2. 7.3. THE QRALGORITHM 213 The third column of Ris thus given by r3=[hq1,a3i,hq2,a3i,kˆq3k,0,0,... , 0]T. In this way we see that the columns of Qare orthogonal and the matrix Ris upper triangular, with an exception. That is the possibility that ˆqk=0f o rs o m e k. In this degenerate case we take qkto be any vector orthogonal to the span of a1,a2,... ,a m,a n dw et a k e rkj=0 , j=k,k+1...m . A l s ow en o t et h a ti fˆ qk=0 ,t h e n akis linearly dependent on a1,a2,... ,a k−1, and hence on q1,q2,...q k−1. Select the coefficients r1k,... ,r k−1kto reflect this dependence. (ii) If m=n, the process above yields a unitary matrix. If Ais nonsingu- lar, the process above yields a matrix Rwith a positive diagonal. (iii) If Ais a real, all operators above can be carried out in real arithmetic. Now what about the uniqueness of the decomposition? Essentially the uniqueness is true up to a multiplication by a diagonal matrix, except inthe case when the matrix has rank is less than m, when there is no form of uniqueness. Suppose that the rank of Aism. Then application of the Gram-Schmidt procedure yields a matrix Rwith positive diagonal. Suppose that Ahas two QR factorizations, QRandPS with upper triangular factors having positive diagonals. Then P ∗Q=SR−1 We have that SR−1is upper triangular and moreover has a positive diagonal. Also, P∗Qis unitary. We know that the only upper triangular unitary matrices are diagonal matrices, and finally the only unitary matrix with a positive diagonal is the identity matrix. Therefore P∗Q=I,w h i c hi st o say that P=Q.We summarize as Corollary 7.3.1. Suppose Ais in Mn,mandn≥m.I f r a n k (A)=m then the QR factorization of A=QRwith upper triangular matrix Rhaving a positive diagonal is unique. 214 CHAPTER 7. FACTORIZATION THEOREMS TheQRalgorithm TheQRalgorithm parallels the LRalgorithm almost identically. Suppose Ais inMnDefine A1=Q1R1 A2≡R1Q1. Also Q∗ 1A1Q=A2. Then decompose A2into a QRdecomposition A2=Q2R2 and Q∗ 2A2Q2=R2Q2≡A3. Also Q∗ 2Q∗1A1Q1Q2=R2Q2=A3. Proceed sequentially Ak=QkRk Ak+1=RkQk Q∗ kAkQk=Ak+1. Let Pk=Q1Q2...Q k Tk=RkRk−1...R 1. Then P∗ kA1Pk=Ak+1. whence PkAk+1=A1Pk. 7.3. THE QRALGORITHM 215 Also we have PkTk=Pk−1QkRkTk−1 =Pk−1AkTk−1 =A1Pk−1Tk−1 =... =Ak 1. Theorem 7.3.2. LetA∈Mnbe given, and assume the eigenvalues of A satisfy |λ1|>|λ2|>···>|λn|>0. Then the iterations Akconverge to a triangular matrix. Proof. Our hypothesis gives that Ais diagonalizable, and we write A∼Λ= diag(λ1...λn). That is, A1=SΛS−1 whereΛ= diag(λ1...λn). Let X=S=QxRx hereQR Y=S−1=LyUyhereLU. Then Ak 1=QxRxΛkLyUy =QxRxΛkLyΛ−kΛkUy =Qx(I+RxEkR−1 x)RxΛkUY where Ek=ΛkLyΛ−k−I (Ek)ij=  0 i=j (λi/λj)k`iji>j 0 i<j . It follows that I+RxEkR−1 x→I,a n d RxΛ−kUyis upper triangular. Thus Qx(I+RxEkR−1 x)RxΛkUy=PkTk. 216 CHAPTER 7. FACTORIZATION THEOREMS The matrix I+RxEkR−1 xcan be QR factored as ˜Uk˜Rk, and since I+ RxEkR−1 x→I, it follows that we can assume both ˜Uk→Iand ˜Rk→I. Hence Ak 1=Qx˜Uk[˜Rk(I+RxEkR−1 x)RxΛkUy]=PkTk. with the first factor unitary and the second factor upper triangular. Since we have assumed (by the eigenvalue condition) that Ais nonsingular, this factorization is essentially unique, where possibly a multiplication by a di- agonal matrix must be applied to give the upper triangular factor on theright a positive diagonal. Just what is the form of the diagonal matrix canbe seen from the following. Let Λ=|Λ|Λ 1,w h e r e |Λ|is the diagonal matrix of moduli of the elements of Λand where Λ1is the unitary matrix of the signs of each eigenvalue respectively. We also take Uy=Λ2(Λ∗ 2Uy)w h e r e Λ2is a unitary matrix chosen so that Λ∗ 2Uhas a positive diagonal. Then Ak 1=Qx˜UkΛ2Λk 1[³ Λ2Λk 1´−1˜Rk(I+RxEkR−1 x)Rx³ Λ2Λk 1´ |Λ|k(Λ∗ 2Uy)] =PkTk. From this we obtain Pkis essentially asymptotic to Qx˜UkΛ2Λk 1and from this we obtain that Qk=P−1 k−1Pk→Λ1 which is diagonal. Finally, it follows that Akis upper triangular since Q−1 kAk=Rk In the limit therefore Ais similar to an upper triangular matrix. Example 7.3.1. Apply the QR method to the matrix A:= 2.31 2 22 2 .1 320  The matrix Ahas eigvenvalues 5 .45,0.723,−1.87. The successive iterations are 7.4. LEAST SQUARES 217 A2= 5.10−0.511 2 .13 0.631 0 .662 0 .136 1.42−0.0202−1.44  A3= 5.51−1.02−0.36 −0.0146 0 .666 0 .482 0.513 0 .240−1.84  A4= 5.46−1.41 0 .482 −0.0372 0 .495 0 .672 0.169 0 .815−1.62  A5= 5.47−0.366−1.26 −0.0404−0.462 1 .39 0.0430 1 .21−0.677  A6= 5.46−1.13−0.687 −0.0184−1.52 0 .813 0.00826 0 .983 0 .381  A7= 5.45 0 .529−1.18 −0.00682−1.78 0 .585 0.00115 0 .414 0 .638  A8= 5.43 0 .684−1.09 −0.000822 −1.87 0 .229 0.0000215 0 .0659 0 .729  Note the gradual appearance of the eigenvalues on the diagonal. Remark. These iterations were carried out in precision 3 arithmetic, whichaffects the rate of convergence to triangular form. 7.4 Least Squares As we know, if A∈Mn,mwith m<n it is generally not possible to solve the overdetermined system Ax=b. For example, suppose we have the data {(xi,yi)}n i=1,w i t ht h e x-coordinates distinct. We may wish to “ fit” a straight to this data. This means we want tofind coefficients mandbso that b+mxi=yi,i =1,... ,n . (?) Taking the matrix and data vector A= 1x 1 1x2 ... 1xn b= y 1 y2 ... yn  andz=[b, m]T, the system ( ?) becomes Az=b.U s u a l l y nÀ2. Hence there is virtually no hope to determine a unique solution to system. However, there are numerous ways to determine constants mandbso that the resulting line represents the data. For example, owing to the dis- tinctness of the x-coordinates, it is possible to solve any 2 ×2 subsystem of 218 CHAPTER 7. FACTORIZATION THEOREMS Az=b. Other variations exist. A new 2 ×2 system could be created by creating two averages of the data, say left and right, and solving. Assume the sequence {xj}is ordered from least to greatest. De finex`=1 kkP j=1xjand xr=1 n−knP j=k+1xj.L e t y`andyrdenote the corresponding averages for the ordinates. Then de fine the intercept band slope mby solving the system ·1x` 1xr¸·b m¸ =·y` yr¸ While this will normally give a reasonable approximating line, its value has little utility beyond its naive simplicity and visual appearance. What is desired is to establish a criteria for choosing the line. Define the residual of the approximation r=b−Az.I tm a k e sp e r f e c t sense to consider finding z=[b, m]Tfor which the residual is minimized in some norm. Any norm can be selected here, but on practical grounds thebest norm to use is the Euclidean norm k·k 2.T h e v e c t o r Azthat yields the minimal norm residual is the one for which ( b−Az)⊥Aw, for we are seeking the nearest value in the Awto the vector b. It can be found by select the one for the solution, Az,f o rw h i c h b−Ax⊥Aw allw. This means hb−Ax, Ay i=0 a l l y or hAT(b−Ay),yi=0 a l l y or AT(b−Ay)=0 ATAy=ATb.Normal Equations Theleast squares solution to Ax=bis given by the solution to the normal equation ATAy=ATb. 7.5. EXERCISES 219 Suppose we have the QRdecomposition for A.T h e ni f Ais real ATA=RTQTQR=RTR ATy=RTQy. Hence the normal equations become RTRx=RTQy. Assuming that the rank of Aism,w em u s th a v et h a t Rand hence RT is invertible. Therefore we have the least squares solution is given by the triangular system Rx=Qy. 7.5 Exercises 1. If A∈M(C)h a sr a n k k, show that there is a permutation matrix P such that PAhas its firstkprincipal determinants nonzero. 2. For the least squares fit of a straight line determine RandQ. 3. In the case of data ATA=·nΣxi ΣxiΣx2 i¸ ATb=·Σyi Σxiyi¸ . 4. In attempting to solve a quadratic fitw eh a v et h em o d e l c+bxi+ax2 i=yi i=1,... ,n . The system is A= 1x1x2 1......... 1xnx2 n  b= y1 y2 ... yn . The normal equations have the matrix and data given by ATA= nΣxiΣx2 i ΣxiΣx2 iΣx3 i Σx2 iΣx2 iΣx4 i  ATb= Σyi Σxiyi Σx2 iyi . 5. Find the normal equations for the least squares fito fd a t at oap o l y - nomial of degree k.