Heath - Scientific Computing (523150), страница 20
Текст из файла (страница 20)
Sparse linear systems with more generalnonzero patterns, on the other hand, require more sophisticated algorithms and data structures in order to avoid storing or operating on the zeros in the matrix (see Section 11.4.1).The properties just defined for real matrices have analogues for complex matrices, butin the complex case the ordinary matrix transpose is replaced by the conjugate transpose,denoted by a superscriptH. If γ = α + iβ is a complex number, where α and β are real√numbers and i = −1, then its complex conjugate is defined by γ̄ = α − iβ. The conjugatetranspose of a matrix A is then given by {AH }ij = āji .
Of course, for a real matrix A,AH = AT . A complex matrix is Hermitian if A = AH , and positive definite if xH Ax > 0for all complex vectors x 6= o.2.5.1Symmetric Positive Definite SystemsIf the matrix A is symmetric and positive definite, then an LU factorization can be arrangedso that U = LT , that is, A = LLT , where L is lower triangular and has positive diagonal64CHAPTER 2.
SYSTEMS OF LINEAR EQUATIONSentries (but not, in general, a unit diagonal). This is known as the Cholesky factorization ofA, and an algorithm for computing it can be derived simply by equating the correspondingentries of A and LLT and then generating the entries of L in the correct order. In the 2 × 2case, for example, we have a11 a21l11 0l11 l21,=a21 a220 l22l21 l22which implies thatl11 =√a11 ,l21 = a21 /l11 ,l22 =q2 .a22 − l21One way to write the resulting general algorithm, in which the Cholesky factor L overwritesthe original matrix A, is as follows:for j = 1 to n{ for each column j }for k = 1 to j − 1{ loop over all prior columns k }for i = j to n{ subtract a multiple ofaij = aij − aik · ajkcolumn k from column j }endend√ajj = ajjfor k = j + 1 to n{ scale column j by squareakj = akj /ajjroot of diagonal entry }endendA number of facts about the Cholesky factorization algorithm make it very attractive andpopular for symmetric positive definite matrices:• The n square roots required are all of positive numbers, so the algorithm is well-defined.• No pivoting is required for numerical stability.• Only the lower triangle of A is accessed, and hence the upper triangular portion neednot be stored.• Only about n3 /6 multiplications and a similar number of additions are required.Thus, Cholesky factorization requires only about half as much work and half as muchstorage as are required for LU factorization of a general matrix by Gaussian elimination.Unfortunately, taking advantage of this gain in storage usually requires that one triangle ofthe symmetric matrix be packed into a one-dimensional array, which is less convenient thanthe usual two-dimensional storage for a matrix.
For this reason, linear algebra softwarepackages commonly offer both packed storage and standard two-dimensional array storageversions for symmetric matrices so that the user can choose between convenience and storageconservation.In some circumstances it may be advantageous to express the Cholesky factorization inthe form A = LDLT , where L is unit lower triangular and D is diagonal with positivediagonal entries. Such a factorization can be computed by a simple variant of the standardCholesky algorithm, and it has the advantage of not requiring any square roots.
The2.5. SPECIAL TYPES OF LINEAR SYSTEMS65diagonal entries of D in the LDLT factorization are simply the squares of the diagonalentries of L in the LLT factorization.Example 2.12 Cholesky Factorization. To illustrate the algorithm, we compute theCholesky factorization of the symmetric positive definite matrix5.0 0 2.5A = 0 2.50 .2.5 0 2.125The successive transformations of the lower triangle of the matrix will be shown, as thealgorithm touches only this portion of the matrix. The first columnhas no prior columns,√so it is merely scaled by the square root of the diagonal entry, 5, to give2.236 0.2.51.118 0 2.125The second column now requires updating by subtracting a multiple of the first column.But in this case the multiplier in the second row of the first column is zero, so that thesecond column is unaffected by the first column.Thus, the second column is simply scaled√by the square root of its diagonal entry, 2.5, to give2.236 0.1.5811.11802.125Finally, the third column must be updated by subtracting multiples of the previous twocolumns.
The multipliers for the first two columns, found in the third row, are 1.118 andzero, respectively. Updating the third column accordingly gives2.236 0.1.5811.11800.875Taking the square root of the third diagonal entry then yields the final result2.236.L= 01.5811.11800.9352.5.2Symmetric Indefinite SystemsIf the matrix A is symmetric but indefinite (i.e., xT Ax can take on both positive andnegative values), then Cholesky factorization is not applicable, and some form of pivoting isgenerally required for numerical stability.
Obviously, any pivoting must be symmetric—of66CHAPTER 2. SYSTEMS OF LINEAR EQUATIONSthe form P AP T , where P is a permutation matrix—if the symmetry of the matrix is tobe preserved.We would like to obtain a factorization of the form P AP T = LDLT , where L is unitlower triangular and D is diagonal.
Unfortunately, such a factorization, with diagonal D,may not exist, and in any case it generally cannot be computed stably using only symmetricpivoting. The best we can do is to take D to be either tridiagonal or block diagonal with1 × 1 and 2 × 2 diagonal blocks. (A block matrix is a matrix whose entries are partitionedinto submatrices, or “blocks,” of compatible dimensions. In a block diagonal matrix, all ofthese submatrices are zero except those on the main block diagonal.)Efficient algorithms have been developed by Aasen for the tridiagonal factorization, andby Bunch and Parlett (with subsequent improvements in the pivoting procedure by Bunchand Kaufman) for the block diagonal factorization (see [104]).
In either case, the pivotingprocedure yields a stable factorization that requires only about n3 /6 multiplications and asimilar number of additions. Also, in either case, the subsequent solution phase requiresonly O(n2 ) work. Thus, the cost of solving symmetric indefinite systems is similar to thatfor positive definite systems using Cholesky factorization, and only about half the cost fornonsymmetric systems using Gaussian elimination.2.5.3Band SystemsGaussian elimination for band matrices differs little from the general case—the only algorithmic changes are in the ranges of the loops. Of course, one should also use a datastructure for a band matrix that avoids storing zero entries. A common choice when theband is dense is to store the matrix in a two-dimensional array by diagonals.
If pivoting isrequired for numerical stability, then the algorithm becomes slightly more complicated inthat the bandwidth can grow (but no more than double). Thus, a general-purpose bandsolver for arbitrary bandwidth is very similar to a code for Gaussian elimination for generalmatrices.For a fixed small bandwidth, however, a band solver can be extremely simple, especiallyif pivoting is not required for stability. Consider, for example, the tridiagonal matrixb1a 2A= 0 . ..0c10b2c2......... an−1···0···......bn−1an0.. . .0 cn−1bnIf pivoting is not required for stability, which is often the case for tridiagonal systems arisingin practice (e.g., if the matrix is diagonally dominant or positive definite), then Gaussianelimination reduces to the following simple algorithm:d 1 = b1for i = 2 to nmi = ai /di−1di = bi − mi ci−1end2.6.
ITERATIVE METHODS FOR LINEAR SYSTEMS67and the resulting LU factorization of A is given by1 m2L= 0 . ..001......······......···......mn−101mn0.. ... ,.01U =d1c10......d2......···00···......c2...... dn−1···00.. . .0 cn−1dnIn general, a band system of bandwidth β requires only O(βn) storage, and the factorization requires only O(β 2 n) work, both of which represent substantial savings over fullsystems if β n.2.6Iterative Methods for Linear SystemsGaussian elimination is an example of a direct method for solving linear systems, i.e., onethat produces the exact solution (assuming exact arithmetic) to a linear system in a finitenumber of steps.
Iterative methods for solving linear systems begin with an initial estimatefor the solution and successively improve it until the solution is as accurate as desired. Intheory, an infinite number of iterations might be required to converge to the exact solution,but in practice the iterations terminate after the residual kb − Axk, or some other measureof error, is as small as desired. For some types of problems, iterative methods may havesignificant advantages over direct methods.
Iterative methods for solving linear systemswill be postponed until Chapter 11, where we consider the numerical solution of partialdifferential equations, which leads to sparse linear systems that are often best solved byiterative methods.2.7Software for Linear SystemsAlmost any software library for scientific computing contains routines for solving linearsystems of various types. Table 2.1 is a list of appropriate routines for solving real, general,dense linear systems, and also for estimating the condition number, in some widely availablesoftware collections. Some packages use different prefixes or suffixes in the routine namesto indicate the data type, typically s for single-precision real, d for double-precision real, cfor single-precision complex, and z for double-precision complex; only the single-precisionreal versions are listed here.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
