Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 41
Текст из файла (страница 41)
In tests during the development of LINPACK, the largest valueobserved was= 23, occurring for a random matrix of 1s, 0s, and -1s [307,1979, p. 1.21]. Macleod [720, 1989] recorded a value= 35.1, which occurred for a symmetric matrix with elements from the uniform distributionon [-1, 1]. In one MATLAB test of 1000 matrices of dimension 100 from thenormal N(0, 1) distribution, I found the largest growth factor to be= 9.59.Gould [474, 1991] used the optimization LANCELOT [236, 1992] to maximize the nth pivot for complete pivoting as a function of about n 3/3 variablescomprising the intermediate elementsof the elimination; constraints wereincluded that normalize the matrix A, describe the elimination equations, andimpose the complete pivoting conditions.
Gould’s package found many localmaxima, and many different starting values had to be tried in order to locate the matrix for which> 13. In an earlier attempt at maximizing thegrowth factor, Day and Peterson [271, 1988] used a problem formulation inwhich the variables are the n2 elements of A, which makes the constraints andobjective function substantially more nonlinear than in Gould’s formulation.Using the package NPSOL [444, 1986], they obtained “largest known” growthfactors for 5 < n < 7.Theoretical progress into understanding the behaviour of the growth factor for complete pivoting has been made by Day and Peterson [271, 1988],Puschmann and Cortés [849, 1983], Puschmann and Nordio [850, 1985], andEdelman and Mascarenhas [345, 1995].A novel alternative to partial pivoting for stabilizing GE is proposed byStewart [942, 1974].
The idea is to modify the pivot element to make it suitably large, and undo this rank one change later using the Sherman-Morrisonformula. Stewart gives error analysis that bounds the backward error for thismodified form of GE.Theorem 9.8 is proved for matrices diagonally dominant by columns byWilkinson [1085, 1961, pp. 288-289]. Theorem 9.9 is proved in the same paper.That pn < 2 for matrices diagonally dominant by rows does not appear tobe well known, but it is proved by Wendroff [1074, 1966, pp. 122-123].
forexample.The results in §9.5 for tridiagonal matrices are taken from Higham [541,1990 ]. Another method for solving tridiagonal systems is cyclic reduction,which was developed in the 1960s [171, 1970]. Error analysis given by Amodioand Mazzia [15, 1994] shows that cyclic reduction is normwise backward stablefor diagonally dominant tridiagonal matrices.The chapter “Scaling Equations and Unknowns” of Forsythe and Moler[396, 1967] is a perceptive, easy to understand treatment that is still well worthreading. Early efforts at matrix scaling for GE were directed to equilibratingeither just the rows or the rows and columns simultaneously (so that all therows and columns have similar norms).
An algorithm with the latter aimis described by Curtis and Reid [259, 1972]. Other important references on198LU FACTORIZATIONANDLINEAR E QUATIONSscaling are the papers by van der Sluis [1040, 1970] and Stewart [945, 1977].which employ normwise analysis, and those by Skeel [919, 1979], [921, 1981],which use componentwise analysis.Much is known about the existence and stability of LU factorizations ofM-matrices and related matrices. A is an H-matrix if the comparison matrixM(A) (defined in (8.6)) is an M-matrix. Funderlic, Neumann, and Plemmons [410, 1982] prove the existence of an LU factorization for an H-matrixA that is generalized diagonally dominant, that is, DA is diagonally dominant by columns for some nonsingular diagonal matrix D: they show thatthe growth factor satisfies pn < 2 maxi |dii |/min i |dii |.
Neumann and Plemmons [791, 1984] obtain a growth factor bound for an inverse of an H -matrix.Ahac, Buoni, and Olesky [7, 1988] describe a novel column-pivoting schemefor which the growth factor can be bounded by n then A is an H-matrix.The normwise bounds in Theorem 9.14 are due to Barrlund [71, 1991]and the componentwise ones to Sun [972, 1992]. Similar bounds are givenby Stewart [951, 1993] and Sun [973, 1992].
Barrlund [72, 1992] describes ageneral technique for deriving matrix perturbation bounds using integrals.Interval arithmetic techniques (see §24.4) are worth considering if high accuracy or guaranteed accuracy is required when solving a linear system. Wemention just one paper, that by Demmel and Krückeberg [297, 1985], whichprovides a very readable introduction to the subject and contains further references.For several years Edelman has been collecting information on the solutionof large, dense linear algebra problems. His papers [337, 1991], [341, 1993],[342, 1994] present statistics and details of the applications in which largedense problems arise. Edelman also discusses relevant issues such as whatusers expect of the computed solutions and how best to make use of parallelcomputers. Table 9.2 contains “world records” for linear systems from Edelman’s surveys.
For all the records shown the matrix was complex and thesystem was solved in double precision arithmetic by some version of LU factorization. Most of the very large systems currently being solved come fromthe solution of boundary integral equations, a major application being theanalysis of radar cross sections; the resulting systems have coefficient mat rices that are complex symmetric (but not Hermitian). A recent reference isWang [1064, 1991].9.10.1.
LAPACKDriver routines xGESV (simple) and xGESVX (expert) use LU factorization withpartial pivoting to solve a general system of linear equations with multipleright-hand sides. The expert driver incorporates iterative refinement, condition estimation, and backward and forward error estimation and has an optionto scale the system AX = B to=before solution,199P ROBLEMSTable 9.2. Records for largest dense linear systems solved (dimension n).Year19911992/319941995n55,29675,26476,800128,600ComputerConnection Machine CM-2Intel iPSC/860Connection Machine CM-5Intel ParagonTime4.4 days22/3days4.1 days1 hourwhere D R = diag(ri ) = diag(maxj |aij|) and DC = diag(maxi ri |aij|); thescaling is done by the routine xGEEQU.
The LU factorization is computed bythe routine xGETRF, which uses a partitioned outer product algorithm. Theexpert driver also returns the quantity ||A||/||U||, where ||A|| := maxi,j |aij|,which is an estimate of the reciprocal of the growth factor,A valuemuch less than 1 signals that the condition estimate and forward error boundcould be unreliable.For band matrices, the driver routines are xGBSV and xGBSVX, and fortridiagonal matrices, xGTSV and xGTSVX; again, these use LU factorizationwith partial pivoting.Problems9.1. (Completion of proof of Theorem 9.1.) Show that if a singular matrixhas a unique LU factorization then Ak is nonsingular for k =1 :n - 1.9.2.
Define A(σ) = A - σI, whereandFor how manyvalues of σ, at most, does A(σ) fail to have an LU factorization withoutpivoting?9.3. Show thatto the field of valueshas a unique LU factorization if 0 does not belong9.4. State analogues of Theorems 9.3 and 9.4 for LU factorization with rowand column interchanges: PAQ = LU.9.5. Give a 2 × 2 matrix A having an LU factorization A = LU such that|L||U| < c|A| does not hold for any c, yetis of order 1.9.6. Show that ifis nonsingular and totally nonnegative it has anLU factorization A = LU with L > 0 and U > 0.
(Hint: use the inequalitywhich holds for any totally nonnegative A [414, 1959, p. 100].) Deduce thatthe growth factor for GE without pivoting p n1.200LU FACTORIZATIONANDLINEAR E QUATIONS9.7. Show that ifis nonsingular and its inverse is totally nonnegative then it has an LU factorization A = LU with |A| = |L||U| . (Use the factthat if C is totally nonnegative and nonsingular then JC - 1 J is totally nonnegative, where J = diag((-1) i +1) (this can be proved using determinantalidentities; see [21, 1987, Thm. 3.3]).)9.8. Show that Theorem 9.5 is valid for GE without pivoting, with a differentconstant.9.9. Suppose that GE without pivoting is applied to a linear system Ax = b ,whereis nonsingular, and that all operations are performed exactlyexcept for the division determining a single multiplier lij (where i > j andA = LU ), which is computed with relative errorl i j = l i j(1 + ).
Evaluatethe difference x between the exact and computed solutions. (The answerallows us to see clearly the effect of a computational blunder, which could, forexample. be the result of the malfunction of a computer’s divide operation.)9.10. Show that θ in Theorem 9.7 satisfiesHence, for g(n) defined in (9.14) and Sn in (9.11). deduce a larger lower boundthan g (2n) >9.11. Explain the errors in the following criticism of GE with complete pivoting.Gaussian elimination with complete pivoting maximizes the pivotat each stage of the elimination. Since the product of the pivots isthe determinant (up to sign), which is fixed, making early pivotslarge forces later ones to be small.
These small pivots will have largerelative errors due to the accumulation of rounding errors during thealgorithm, and dividing by them therefore introduces larger errors.9.12. In sparse matrix computations the choice of pivot in GE has to bemade with the aim of preserving sparsity as well as maintaining stability. Inthreshold pivoting, a pivot element is chosen from among those elements incolumn k that satisfywhereis a parameter(see, for example, Duff, Erisman, and Reid [325, 1986, §5.4]).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
