Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 58
Текст из файла (страница 58)
Bischof [103, 1990] develops amethod for estimating the smallest singular value of a triangular matrix whichprocesses the matrix a row or a column at a time. This “incremental conditionestimation” method can be used to monitor the condition of a triangular matrix as it is generated, and so is useful in the context of matrix factorizationsuch as the QR factorization with column pivoting. The estimator is generalized to sparse matrices by Bischof, Lewis, and Pierce [104, 1990]. Barlow andVemulapati [67, 1992] develop a l-norm incremental condition estimator withlook-ahead for sparse matrices.Condition estimates are also required in applications where a matrix factorization is repeatedly updated as a matrix undergoes low rank changes.Algorithms designed for a recursive least squares problem and employingthe Lanczos method are described by Ferng, Golub, and Plemmons [372,1991 ].
Pierce and Plemmons [831, 1992 ] describe an algorithm for use withthe Cholesky factorization as the factorization is updated, while Shroff andBischof [918, 1992] treat the QR factorization.14.5. Condition Numbers of Tridiagonal MatricesFor a bidiagonal matrix B, |B–1| = M(B)–1 (see §8.3), so the condition numbers κ E,f and condE,f can be computed exactly with an order of magnitudeless work than is required to compute B-1 explicitly. This property holdsmore generally for several types of tridiagonal matrix, as a consequence ofthe following result. Recall that the LU factors of a tridiagonal matrix arebidiagonal and may be computed using the formulae (9.16).Theorem 14.7.
If the nonsingular tridiagonal matrixhas the LUfactorization A = LU and |L||U| = |A|, then |U - 1 |L- 1 | = |A- 1|.Proof. Using the notation of (9.15), |L||U| = |A| = |LU| if and only if,for all 2,that is, if(14.9)Using the formulae302CONDITION NUMBER ESTIMATIONwe haveThus, in view of (14.9), it is clear that |U - 1 L - 1|ij = (|U- 1 ||L - 1 |) z j, as required.Since L and U are bidiagonal, |U- 1 | = M(U)-1 and |L- 1 | = M(L)- 1 .Hence, if |A| = |L||U|, then, from Theorem 14.7,(14.10)It follows that we can compute any of the condition numbers or forward errorbounds of interest exactly by solving two bidiagonal systems.
The cost isO(n) flops, as opposed to the O(n2) flops needed to compute the inverse of atridiagonal matrix.When does the condition |A| = |L||U| hold? Theorem 9.11 shows that itholds if the tridiagonal matrix A is symmetric positive definite, totally positive, or an M-matrix. So for these types of matrix we have a very satisfactoryway to compute the condition number.If A is tridiagonal and diagonally dominant by rows, then we can computein O(n) flops an upper bound for the condition number that is not more thana factor 2n – 1 too large.14.5 C ONDITION N UMBERSOFT RIDIAGONAL M ATRICES303Theorem 14.8.
Suppose the nonsingular, row diagonally dominant tridiagonal matrixhas the LU factorization A = LU. Then, if y > 0,Proof. We have L-1 = UA– 1 , sowhere the bidiagonal matrix V = diag(u i i )-1 U has υi i|ei /ui | < 1 (see the proof of Theorem 9.12). Thus1 and |υi , i +1| =and the result follows on taking norms.In fact, it is possible to compute |A- 1|y exactly in O(n) operations for anytridiagonal matrix. This is a consequence of the special form of the inverse ofa tridiagonal matrix.Theorem 14.9 (Ikebe).
Letbe tridiagonal and irreducible (thatis, ai+ 1 ,i and ai,i+1 are nonzero for all i). Then there are vectors x, y, p, andq such thatThis result says that the inverse of an irreducible tridiagonal matrix isthe upper triangular part of a rank-1 matrix joined along the diagonal to thelower triangular part of another rank-1 matrix. If A is reducible then it has(or its transpose), and this blocking can be appliedthe block formrecursively until the diagonal blocks are all irreducible, at which point thetheorem can be applied to the diagonal blocks.The vectors x, y, p, and q in Theorem 14.9 can all be computed in O(n)flops, and this enables the condition numbers and forward error bounds to becomputed also in O(n) flops (see Problem 14.5). Unfortunately, the vectorsx, y, p, and q can have a huge dynamic range, causing the computation tobreak down because of overflow or underflow.
For example, for the diagonallydominant tridiagonal matrix with a ii4, ai+ 1 ,i = ai,i+ 1 1, we have (x1 =1) |x n | θ n -1, |y1|θ -1, and |yn|θ - n , where θ = 2 +3.73. Thesenumerical problems can be overcome, but only at a nontrivial increase incost. Therefore, we do not recommend the use of Theorem 14.9 for computingcondition numbers. For a general tridiagonal matrix it is probably better toestimate the condition number using Algorithm 14.4.304CONDITION NUMBER ESTIMATION14.6.
Notes and ReferencesThe clever trick (14.1) for converting the norminto the norm ofa matrix with which products are easily formed is due to Arioli, Demmel, andDuff [24, 1989].The p-norm power method was first derived and analysed by Boyd [139,1974] and was later investigated by Tao [993, 1984]. Tao applies the methodto an arbitrary mixed subordinate norm ||A||a,β (see (6.6)), while Boyd takesthe a and β-norms to be p-norms (possibly different). Algorithm 14.1 can beconverted to estimate ||A||a,β by making straightforward modifications to thenorm-dependent terms. An algorithm that estimates ||A||p using the powermethod with a specially chosen starting vector is developed by Higham [551,1992]; the method for obtaining the starting vector is outlined in Problem 14.1.The estimate produced by this algorithm is always within a factor n 1 – 1 /p of||A||p and the cost is about 70n 2 flops.
A MATLAB M-file pnorm implementingthis method is part of the Test Matrix Toolbox (see Appendix E).The finite convergence of the power method for p = 1 and p = 00 holdsmore generally: if the power method is applied to the norm ||·|| a,β andone of the a and β norms is polyhedral (that is, its unit ball has a finitenumber of extreme points), then the iteration converges in a finite numberof steps. Moreover, under a reasonable assumption, this number of steps canbe bounded in terms of the number of extreme points of the unit balls in thea-norm and the dual of the β-norm. See Bartels [75, 1991] and Tao [993,1984] for further details.Hager [492, 1984] gave a derivation of the 1-norm estimator based on subgradients and used the method to estimate κ 1 (A).
That the method is ofwide applicability y because it accesses A only through matrix-vector productswas recognized by Higham, who developed Algorithm 14.4 and its complexanalogue and wrote Fortran 77 implementations, which use a reverse communication interface [537, 1988], [543, 1990]. These codes are used in LAPACK,the NAG library, and various other program libraries. A version of Algorithm 14.4 dedicated to estimating κ 1 (A) is supplied with MATLAB as M-filecondest.
Algorithm 14.4 is also implemented in ROM on the Hewlett-PackardHP 48G and HP 48GX calculators (along with several other LAPACK routines), in a form that estimates κ 1 (A). The Hewlett-Packard implementationis instructive because it shows that condition estimation can be efficient evenfor small dimensions: on a standard HP 48G, inverting A and estimating itscondition number (without being given a factorization of A in either case)both take about 5 seconds for n = 10, while for n = 20 inversion takes 30seconds and condition estimation only 20 seconds.Moler [769, 1978] describes an early version of the LINPACK conditionestimator and raises the question of the existence of counterexamples.
Anearly version without look-ahead was incorporated in the Fortran code decomp14.6 N OTESANDR EFERENCES305in the book of Forsythe, Malcolm, and Moler [395, 1977].Matrices for which condition estimators perform poorly can be very hardto find theoretically or with random testing, but for all the estimators described in this chapter they can be found quite easily by applying directsearch optimization to the under- or overestimation ratio; see §24.3.1.Both LINPACK and LAPACK return estimates of the reciprocal of thecondition number, in a variable rcond < 1. Overflow for a very ill conditioned matrix is thereby avoided, and rcond is simply set to zero when singularity is detected.
MATLAB has a built-in function rcond that implementsthe LINPACK condition estimation algorithm.A simple modification to the LINPACK estimator that can produce alarger estimate is suggested by O’Leary [804, 198 0]. For sparse matrices,Grimes and Lewis [483, 1981] suggest a way to reduce the cost of the scalingstrategy used in LINPACK to avoid overflow in the condition estimation.Zlatev, Wasniewski, and Schaumburg [1134, 1986] describe their experiencein implementing the LINPACK condition estimation algorithm in a softwarepackage for sparse matrices.Stewart [946, 1980] describes an efficient way to generate random matricesof a given condition number and singular value distribution (see §26.3) andtests the LINPACK estimator on such random matrices.Condition estimators specialized to the (generalized) Sylvester equationhave been developed by Byers [172, 1984], Kågström and Westin [624, 1989],and Kågström and Poromaa [621, 1992].A survey of condition estimators up to 1987, which includes counterexamples and the results of extensive numerical tests, is given by Higham [534,19 8 7 ].Theorems 14.7 and 14.8 are from Higham [541, 1990].
Thatcanbe computed in O(n) flops for symmetric positive definite tridiagonal A wasfirst, shown in Higham [531, 1986].Theorem 14.9 has a long history, having been discovered independentlyin various forms by different authors. The earliest reference we know forthe result as stated is Ikebe [600, 1979], where a more general result forHessenberg matrices is proved. A version of Theorem 14.9 for symmetrictridiagonal matrices was proved by Bukhberger and Emel’yanenko [157, 1973].The culmination of the many papers on inverses of tridiagonal and Hessenbergmatrices is a result of Cao and Stewart on the form of the inverse of a blockmatrix (Aij) with Aij = 0 for i > j + s [184, 1986]; despite the generality ofthis result, the proof is short and elegant.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
