Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 59
Текст из файла (страница 59)
Any banded matrix has an inversewith a special “low rank” structure; the earliest reference on the inverse of ageneral band matrix is Asplund [32, 1959]. For a recent survey on the inversesof symmetric tridiagonal and block tridiagonal matrices see Meurant [751,1992 ].For symmetric positive definite tridiagonal A the standard way to solve306C ONDITION N UMBER E STIMATIONAX = b is by using a Cholesky or LDLT factorization, rather than an LUfactorization.
The LINPACK routine SPTSL uses a nonstandard “LUB” factorization resulting from the BABE (“burn at both ends”) algorithm, whicheliminates from the middle of the matrix to the top and bottom simultaneously (see the LINPACK Users’ Guide [307, 1979, Chap. 7] and Higham [531,19 86] ). The results of §14.5 are applicable to all these factorization, withminor modifications.14.6.1. LAPACKAlgorithm 14.4 is implemented in routine xLACON , which has a reverse communication interface. The LAPACK routines xPTCON and xPTRFS for symmetric positive definite tridiagonal matrices compute condition numbers using(14.10); the LAPACK routines xGTCON and xGTRFS for general tridiagonalmatrices use Algorithm 14.4. LINPACK’S tridiagonal matrix routines do notincorporate condition estimation.Problems14.1.
(Higham [551, 1992]) The purpose of this problem is to develop a noniterative method for choosing a starting vector for Algorithm 14.1. The ideais to choose the components of x in the order x1, x2, . . . . xn in an attemptto maximize ||Ax||p/||x||p. Suppose x1, . . . . xk-1 satisfying ||x(1:k-1)||p = 1have been determined and let γk -1 = ||A( : , 1 :k-1)x(1:k-1)||p. We now tryto choose xk, and at the same time revise x(1:k –1), to give the next partialproduct a larger norm. Definingwe setwhereThen ||x(1:k)||p = 1 andDevelop this outline into a practical algorithm. What can you prove aboutthe quality of the estimate ||Ax||p /||x||p that it produces?307PROBLEMS14.2. (Higham [543,(tij ) be defined by1990])Let the n × n symmetric tridiagonal matrix Tn ( a ) =For example, T6 (a) is given byNote that, for all a, ||Tn (a)e n- 1||1 = ||Tn (a)||1.
Show that if Algorithm 14.3is applied to Tn (a) with 0 < a < 1 then x = ei -1 on the ith iteration, fori = 2, . . ., n, with convergence on the nth iteration. Algorithm 14.4, however,terminates after five iterations with y5 = Tn (a)e4, andShow that the extra estimate saves the day, so that Algorithm 14.4 returns afinal estimate that is within a factor 3 of the true norm, for any a < 1.14.3. Let PA = LU be an LU factorization with partial pivoting of A 1Show that14.4.
(Higham [537, 1988]) Investigate the behaviour of Algorithms 14.3 and14.4 for the Pei matrix, A = aI + eeT ( a > 0), and for the upper bidiagonalmatrix with 1s on the diagonal and the first superdiagonal.14.5. (Ikebe [600, 1979], Higham [531, 1986]) Letbe nonsingular,tridiagonal, and irreducible. By equating the last columns in AA -1 = Iand the first rows in A –1 A = I, show how to compute the vectors x andy in Theorem 14.9 in O(n) flops. Hence obtain an O(n) flops algorithm forcomputingwhere d > 0.308CONDITION NUMBER ESTIMATION14.6.
The representation of Theorem 14.9 for the inverse of nonsingular, tridiagonal, and irreducibleinvolves 472 parameters, yet A depends onlyon 3n – 2 parameters. Obtain an alternative representation that involves only3n – 2 parameters. (Hint: symmetrize the matrix.)14.7. (R ESEARCH P ROBLEM ) (Demmel [286, 1992 ]) Show that estimating||A - 1|| to within a factor depending only on the dimension of A is at least asexpensive as computing A– 1 .14.8. (RESEARCH PROBLEM) Letbe diagonally dominant by rows,let A = LU be an LU factorization, and let y > 0.
What is the maximumsize ofThis is an open problem raised in [541,1990]. In a small number of numerical experiments with full random matricesthe ratio has been found to be less than 2 [541, 1990], [790, 1986].PreviousHomeChapter 15The Sylvester EquationWe must commence, not with a square,but with an oblong arrangement of terms consisting, suppose,of m lines and n columns.This will not in itself represent a determinant,but is, as it were, a Matrix out of which we may formvarious systems of determinants by fixing upon a number p,and selecting at will p lines and p columns,the squares corresponding to which may be termeddeterminants of the pth order.— J.
J. SYLVESTER, Additions to the Articles, “On a New Classof Theorems, ” and “On Pasta/’s Theorem” (1850)I have in previous papers defined a “Matrix” as a rectangular array of terms,out of which different systems of determinants may be engendered,as from the womb of a common parent;these cognate determinants beingby no means isolated in their relations to one another,but subject to certain simple laws ofmutual dependence and simultaneous deperition.— J.
J. SYLVESTER, On the Relation Between the MinorDeterminants of Linear/y Equivalent Quadratic Functions (1851)309NextTHE SYLVESTER EQUATION310The linear matrix equationAX – XB = C,(15.1)whereandare given andis to be determined, is called the Sylvester equation. It is of pedagogicalinterest because it includes as special cases several important linear equationproblems:1 . linear system: Ax = c,2 . multiple right-hand side linear system: AX = C,3 . matrix inversion: AX = I,4 . eigenvector corresponding to given eigenvalue b: (A – bI)x = 0,5 .
commuting matrices: AX – XA = 0.The Sylvester equation arises in its full generality in various applications. Forexample, the equationsshow that block-diagonalizing a block triangular matrix is equivalent to solving a Sylvester equation. The Sylvester equation can also be produced fromfinite difference discretization of a separable elliptic boundary value problemon a rectangular domain, where A and B represent application of a differenceoperator in the “y” and “x” directions, respectively [935, 1991].That (15.1) is merely a linear system is emphasized by writing it in theform(15.2)where AB := (aij B) is a Kronecker product and the vec operator stacksthe columns of a matrix into one long vector. For future reference, we notethe useful relation(See Horn and Johnson [581, 1991, Chap.
4] for a detailed presentation ofproperties of the Kronecker product and the vec operator). The mn × mncoefficient matrix in (15.2) has a very special structure, illustrated for n = 3by15.1 S OLVINGTHES YLVESTER E QUATION311In dealing with the Sylvester equation it is vital to consider this structure andnot treat (15.2) as a general linear system.Since the mn eigenvalues of InA – BTIm are given by(15.3)the Sylvester equation is nonsingular precisely when A and B have no eigenvalues in common.In this chapter we briefly discuss the Schur method for solving the Sylvesterequation and summarize its rounding error analysis. Then we determine thebackward error for the Sylvester equation, investigate its relationship withthe residual, and derive a condition number.
All these results respect thestructure of the Sylvester equation and are relevant to any solution method.We also consider the special case of the Lyapunov equation and mention howthe results extend to generalizations of the Sylvester equation.15.1. Solving the Sylvester EquationOne way to solve the Sylvester equation is to apply Gaussian elimination withpartial pivoting (GEPP) to the “big” system (15.2), but the structure of thecoefficient matrix cannot be exploited and the cost is a prohibitive O( m 3 n 3 )flops. A more efficient method, requiring O( m 3 + n3) flops, is obtained withthe aid of Schur decompositions of A and B. Let A and B have the real SchurdecompositionsA = URU T ,(15.4)B = VSVT,where U and V are orthogonal and R and S are quasi-triangular, that is, blocktriangular with 1 × 1 or 2 × 2 diagonal blocks, and with any 2 × 2 diagonalblocks having complex conjugate eigenvalues.
(See Golub and Van Loan [470,1989, §7.4.1] for more details of the real Schur decomposition.)With the decompositions (15.4), the Sylvester equation transforms to(15.5)or, equivalently, Pz = d, where P = I nR – STI m , z = vec(Z) andd = v e t(D). If R and S are both triangular then P is block triangularwith triangular diagonal blocks, so Pz = d can be solved by substitution.Expressed in the notation of (15.5), the solution process take the form of nsubstitutions: if S is upper triangular thenSuppose now that R and S are quasi-triangular, and for definiteness assume that they are both upper quasi-triangular. Partitioning Z = (Zij ) con-312T HE S YLVESTER E QUATIONformally with R = (Rij ) and S = (Sij ) we have(15.6)These equations can be used to determine the blocks of Z working up theblock columns from first to last.
Since Rii and Sjj are both of order 1 or 2,each system (15.6) is a linear system of order 1, 2, or 4 for Zij ; in the lattertwo cases it is usually solved by GEPP (or even Gaussian elimination withcomplete pivoting—see Problem 15.4).This Schur decomposition method for solving the Sylvester equation isdue to Bartels and Stewart [74, 1972].
What can be said about its numericalstability? In the case where R and S are both triangular, Theorem 8.5 showsthat(15.7)where cm,n denotes a constant depending on the dimensions m and n (in fact,we can take c m,nmn ). Thuswhich implies theweaker inequality(15.8)If R or S is quasi-triangular then the error analysis depends on how thesystems of dimension 2 or 4 are solved. If GEPP followed by fixed precisioniterative refinement is used for each of these systemsand if foreach systemis not too ill conditioned and the vectoris not too badlyscaled, then (15.7) and (15.8) remain valid (see §11.2). Otherwise, we haveonly a normwise boundBecause the transformation of a matrix to Schur form (by the QR algorithm)is a backward stable process15 it is true overall that(15.9)Thus the relative residual is guaranteed to be bounded by a modest multipleof the unit roundoff u.Golub, Nash, and Van Loan [464, 1979] suggested a modification of theBartels-Stewart algorithm in which A is reduced only to upper Hessenbergform: A = UHU T .
The reduced system HZ – ZS = D can be solvedby solving n systems that are either upper Hessenberg or differ from upper15See Golub and Van Loan [470, 19 8 9 , §7.5.6]. A proof is outside the scope of this book,but the necessary tools are Lemmas 18.3 and 18.8 about Householder and Givens rotations.15.2 B ACKWARD E RROR313Hessenberg form by the addition of an extra subdiagonal. As shown in [464,1979 ], the Hessenberg–Schur algorithm has a smaller flop count than theBartels-Stewart algorithm, with the improvement depending on the relativesizes of m and n.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
