Heath - Scientific Computing (523150), страница 16

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 16)

But then 1 01 1LU ==6= A.1/ 10 −1/1 0Using a small pivot, and a correspondingly large multiplier, has caused an unrecoverableloss of information in the transformed matrix. If we interchange rows, on the other hand,then the pivot is 1 and the resulting multiplier is −, so that we get the elimination matrix1 0M=,− 12.2. SOLVING LINEAR SYSTEMSand henceL=1 0 147and U = 111=0 1−0in floating-point arithmetic. We therefore have 1 111 0=LU = 10 1111,1which is the correct result after permutation.Although the foregoing example is rather extreme, the principle holds in general thatlarger pivots produce smaller multipliers and hence smaller errors. In particular, if thelargest entry on or below the diagonal in each column is used as pivot (partial pivoting),then the multipliers are bounded in magnitude by 1.

In Example 2.6, we did not use rowinterchanges, and some of the multipliers were greater than 1. For illustration, we nowrepeat that example, this time using partial pivoting. The system in Example 2.6 is   24 −2x12 49 −3x2 =8.−2 −37x310The largest entry in the first column is 4, sopermutation matrix0P1 =  10obtaining the permuted system4P1 Ax =  2−2we interchange the first two rows using the10000,1   9 −3x184 −2   x2  =  2  = P1 b.−37x310To annihilate the subdiagonal entries of the first column, we use the elimination matrix1 0 0M1 =  − 21 1 0  ,10 12obtaining the transformed system49M1 P1 Ax = 0 − 12302  −3x18− 12   x2  =  −2  = M1 P1 b.11x3142The largest entry in the second column on or below the diagonal isthe last two rows using the permutation matrix1 0 0P2 =  0 0 1  ,0 1 032,so we interchange48CHAPTER 2.

SYSTEMS OF LINEAR EQUATIONSobtaining the permuted system x1811  x2  =  14  = P2 M1 P1 b.21−2x3−2−3493P2 M1 P1 Ax =  020 − 12To annihilate the subdiagonal entry of the second1 0M2 = 0 10 31obtaining the transformed system4M2 P2 M1 P1 Ax =  00932column, we use the elimination matrix00,1  x1811   =  14  = M2 P2 M1 P1 b.x2284x333−30We have therefore reduced the original system to an equivalent upper triangular system,which can now be solved by back-substitution to obtain the same answer as before.To write out the LU factorization explicitly, we have010100L = M −1 = (M2 P2 M1 P1 )−1 = P1T L1 P2T L2 =  11 0 010 001 0 021 0 =  10   21 1 0   0 0 1   0− 12 0 10 1 00 − 13 1− 121− 13 10 0,1 0and hence2A= 4−2  14 −229 −3  =  1−37− 12− 13 140 001 009320−311243 = LU .Note that L is not lower triangular, but it is triangular in the more general sense (it is apermutation of a lower triangular matrix).

Alternatively, we can take 1 0 00 1 00 1 0P = P2 P1 =  0 0 1   1 0 0  =  0 0 1  ,0 1 00 0 11 0 0and1L =  − 12120 01 0,− 13 1so that0PA = 0110002140−2 4 −219 −3 = − 121−37240 001 0− 13 109320−311243 = LU ,2.2. SOLVING LINEAR SYSTEMS49where L now really is lower triangular but A is permuted.As we have just seen, pivoting is generally required for Gaussian elimination to be stable.There are some classes of matrices, however, for which Gaussian elimination is stable withoutpivoting.

For example, if the matrix A is diagonally dominant by columns, which meansthat each diagonal entry is larger in magnitude than the sum of the magnitudes of the otherentries in its column,nX|aij | < |ajj |, j = 1, . . . , n,i=1, i6=jthen pivoting is not required in computing its LU factorization by Gaussian elimination.If partial pivoting is used on such a matrix, then no row interchanges will actually occur.Another important class for which pivoting is not required is matrices that are symmetricand positive definite, which will be defined in Section 2.5.

Avoiding an unnecessary pivotsearch can save a significant amount of time in computing the factorization.2.2.5Implementation of Gaussian EliminationGaussian elimination, or LU factorization, has the general form of a triple-nested loop,forforforaij = aij − (aik /akk )akjendendendwhere the indices i, j, and k of the for loops can be taken in any order, for a total of 3! = 6different ways of arranging the loops. Some of the indicated arithmetic operations can bemoved outside the innermost loop for greater efficiency, depending on the specific indicesinvolved, and additional reorderings of the operations that do not have strictly nestedloops are also possible.

These variations of the basic algorithm have different memoryaccess patterns (e.g., accessing memory row-wise or column-wise), and also differ in theirability to take advantage of the architectural features of a given computer (e.g., cache,paging, vectorization, multiple processors). Thus, their performance may vary widely on agiven computer or across different computers, and no single arrangement may be uniformlysuperior.Numerous implementation details of the algorithm are subject to variation in this way.For example, the partial pivoting procedure we described searches along columns and interchanges rows, but alternatively, one could search along rows and interchange columns.

Wehave also taken L to have unit diagonal, but one could instead arrange for U to have unitdiagonal. Some of these variations of Gaussian elimination are of sufficient importance tohave been given names, such as the Crout and Doolittle methods.Although the many possible variations on Gaussian elimination may have a dramaticeffect on performance, they all produce essentially the same factorization. Provided therow pivot sequence is the same, if we have two LU factorizations P A = LU = L̂Û , then50CHAPTER 2.

SYSTEMS OF LINEAR EQUATIONSthis expression implies that L̂−1 L = Û U −1 = D is both lower and upper triangular, andhence diagonal. If both L and L̂ are assumed to be unit lower triangular, then D mustin fact be the identity matrix I, and hence L = L̂ and U = Û , so that the factorizationis unique. Even without this assumption, however, we may still conclude that that theLU factorization is unique up to diagonal scaling of the factors.

This uniqueness is madeexplicit in the LDU factorization P A = LDU , where L is unit lower triangular, U is unitupper triangular, and D is diagonal.Storage management is another important implementation issue. The numerous matrices we considered—the elementary elimination matrices Mk , their inverses Lk , and thepermutation matrices Pk —merely describe the factorization process formally.

They are notformed explicitly in an actual implementation. To conserve storage, the L and U factorsoverwrite the initial storage for the input matrix A, with the transformed matrix U occupying the upper triangle of A (including the diagonal), and the multipliers that make upthe strict lower triangle of L occupying the (now zero) strict lower triangle of A. The unitdiagonal of L need not be stored.To minimize data movement, the row interchanges required by pivoting are not usuallycarried out explicitly.

Instead, the rows remain in their original locations, and an auxiliaryinteger vector is used to keep track of the new row order. Note that a single such vectorsuffices, because the net effect of all of the interchanges is still just a permutation of theintegers 1, . . . , n.2.2.6Complexity of Solving Linear SystemsThe Gaussian elimination process for computing the LU factorization requires about n3 /3floating-point multiplications and a similar number of additions.

Solving the resultingtriangular system for a single right-hand-side vector by forward- and back-substitutionrequires about n2 multiplications and a similar number of additions. Thus, as the order nof the matrix grows, the LU factorization phase becomes increasingly dominant in the costof solving linear systems.We can also solve a linear system by explicitly inverting the matrix so that the solutionis given by x = A−1 b.

But computing A−1 is tantamount to solving n linear systems: itrequires an LU factorization of A followed by n forward- and back-substitutions, one foreach column of the identity matrix. The total operation count is about n3 multiplicationsand a similar number of additions (taking advantage of the zeros in the right-hand-sidevectors for the forward-substitution).Thus, explicit inversion is three times as expensive as LU factorization. The subsequent matrix-vector multiplication x = A−1 b to solve a linear system requires about n2multiplications and a similar number of additions, which is similar to the total cost offorward- and back-substitution.

Hence, even for multiple right-hand-side vectors, matrixinversion is more costly than LU factorization for solving linear systems. In addition, explicit inversion gives a less accurate answer. As a simple example, if we solve the 1 × 1linear system 3x = 18 by division, we get x = 18/3 = 6, but explicit inversion would givex = 3−1 × 18 = 0.333 × 18 = 5.99 using three-digit arithmetic. In this small example, inversion requires an additional arithmetic operation and obtains a less accurate result. Thedisadvantages of inversion become worse as the size of the system grows.2.2. SOLVING LINEAR SYSTEMS51Explicit matrix inverses often occur as a convenient notation in various formulas, butthis practice does not mean that an explicit inverse is required to implement such a formula.One merely need solve a linear system with an appropriate right-hand side, which mightitself be a matrix. Thus, for example, a product of the form A−1 B should be computedby LU factorization of A, followed by forward- and back-substitutions using each columnof B.

It is extremely rare in practice that an explicit matrix inverse is actually needed, sowhenever you see a matrix inverse in a formula, you should think “solve a system” ratherthan “invert a matrix.”Another method for solving linear systems that should be avoided is Cramer’s rule, inwhich each component of the solution is computed as a ratio of determinants. Though oftentaught in elementary linear algebra courses, this method is astronomically expensive for fullmatrices of nontrivial size. Cramer’s rule is useful mostly as a theoretical tool.2.2.7Gauss-Jordan EliminationThe motivation for Gaussian elimination is to reduce a general matrix to triangular form,because the resulting linear system is easy to solve.

Diagonal linear systems are even easierto solve, however, so diagonal form would appear to be an even more desirable target. GaussJordan elimination is a variation of standard Gaussian elimination in which the matrix isreduced to diagonal form rather than merely to triangular form. The same type of rowcombinations are used to eliminate matrix entries as in standard Gaussian elimination,but they are applied to annihilate entries above as well as below the diagonal.

Характеристики

Список файлов книги