Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 61
Текст из файла (страница 61)
Here,sep(A, B) = 1.67 × 10–16, and the bound (15.29) is small because relativelylarge columns of P-1 are nullified by relatively small elements of |vec+vec(R u ). For this example, with a = ||A||F, β = ||B||F, γ = ||C||F, we haveΨ = 7.00 × 109,Φ = 1.70 × 1016,confirming that the usual perturbation bound (15.25) for the Sylvester equation can be very pessimistic. Furthermore,so we have an example where the backward error is small despite a large µ .15.5. ExtensionsThe Sylvester equation can be generalized in two main ways.
One retains thelinearity but adds extra coefficient matrices, yielding the generalized Sylvesterequations(15.30)AXB + CXD = EandAX – YB = C, DX – YE = F.(15.31)These two forms are equivalent, under conditions on the coefficient matrices[210, 1987]; for example, defining Z := XB and W := –CX, (15.30) becomesAZ – WD = E, CZ + WB = 0. Applications of generalized Sylvester equations include the computation of stable eigendecompositions of matrix pencils[294, 198 7], [295, 1988], [622, 1993], [623, 1994] and the implementation ofT HE S YLVESTER E QUATION322numerical methods for solving implicit ordinary differential equations [353,19 8 0 ].The second generalization incorporates a quadratic term, yielding the algebraic Riccati equation(15.32)This general Riccati equation and its symmetric counterpart with B = ATand F and G symmetric are widely used in control theory.The backward error results and perturbation theory of this chapter can begeneralized in a straightforward way to (15.31) and (15.32).
See Kågström [620,1994] for (15.31) and Ghavimi and Laub [440, 1995] for (15.32). The backward error derivations do not extend to (15.30), because in this equation thecoefficient matrices appear nonlinearly.A variation of the Lyapunov equation called the discrete-time Lyapunovequation has the formAs in (15.30), the data appears nonlinearly. Ghavimiwhereand Laub [441, 1995] show how to derive an approximation to the backwarderror by linearizing an equation characterizing the optimal perturbations.Another generalization of the Sylvester equation, mainly of theoreticalinterest, iswhereassociated theory.andSee Lancaster [685,1970]for15.6.
Notes and ReferencesThis chapter is based on Higham [556, 1993]. The backward error derivationsmake use of ideas of Ghavimi and Laub [440, 1995].The Sylvester equation is so named because Sylvester considered the homogeneous version of the equation [985, 1884].Bhatia and Rosenthal [96, 1996] give a survey of theoretical results for theSylvester equation in both finite- and infinite-dimensional spaces.For details of the role of the Sylvester equation in the eigenproblem seeBai, Demmel, and McKenney [38, 1993], [40, 1993] and the references therein.Iterative methods that make use of matrix inversion to solve the Sylvesterequation are described by Miller [755, 1988] and Roberts [875, 1980].Hammarling [496, 198 2] gives a method for solving the Lyapunov equation AX + XAT = –C in the case where A has eigenvalues with negative15.6 N OTESANDR EFERENCES323real parts and C is positive semidefinite; his method directly computes theCholesky factor of the solution (which is indeed symmetric positive definite—see Problem 15.2).A survey of the vec operator, the Kronecker product, and the vec-permutation matrix is given together with historical comments by Henderson andSearle [514, 1981].
Historical research by Henderson, Pukelsheim, and Searle[513, 1983] indicates that the Kronecker product should be called the Zehfussproduct, in recognition of an 1858 paper by Zehfuss that gives a determinantalresult involving the product.The vet-permutation matrix Π (which appears in (15.27)) is given explicitly byand has the property thatApplications of the Lyapunov equation in control theory, including specialsituations where an approximate solution of low rank is required, are discussedby Hodel [574, 1992].
A much older reference to applications is Barnett andStorey [68, 1968].Algorithms and software for solving (15.30) are developed by Gardiner,Wette, Laub, Amato, and Moler [417, 1982], [418, 1982].Perturbation theory for Lyapunov and Riccati equations can be found inthe work of Byers [173, 1985], Hewer and Kenney [522, 1988], [650, 1990], andGahinet, Laub, Kenney, and Hewer [411, 1990].Chu [210, 1987] determines conditions for the existence of unique solutionsto the generalized Sylvester equations (15.30) and (15.31).
The appropriateconditions for (15.30) are that the pencils A +and D +are regularand the spectra of the pencils have an empty intersection, which neatly generalizes the conditions for the Sylvester equation to have a unique solution;the conditions for (15.31) are analogous.There is much work on algorithms and software for solving the algebraicRiccati equation. For a sampling, see Laub [691, 1979], Arnold and Laub [30,1984], Byers [174, 1987], Gardiner, and Laub [416, 1991], and Kenney, Laub,and Papadopoulos [653, 1992].An algorithm for estimating a generalization of sep that occurs in perturbation theory for the generalized Sylvester equation (15.31) is developed byKågström and Westin [624, 1989].Another generalization of the Sylvester equation is to take just one equation from (15.31), AX – YB = C ((15.13) is of this form).
This equation canbe underdetermined or overdetermined, depending on the dimensions of thecoefficient matrices. Conditions involving generalized inverses that are bothnecessary and sufficient for the existence of a solution are given by Baksalaryand Kala [49, 1979]. Zietak examines the inconsistent case [1130, 1985] for oneT HE S YLVESTER E QUATION324choice of dimensions giving an overdetermined system. Stewart [949, 1992]shows how to compute a minimum Frobenius norm least squares solution.The even more general equation AXB + CYD = E has also been analysedby Baksalary and Kala [50, 198 0], who again give necessary and sufficientconditions for the existence of a solution.15.6.1. LAPACKThe computations discussed in this chapter can all be done using LAPACK.The Bartels-Stewart algorithm can be implemented by calling xGEES to compute the Schur decomposition, using the level-3 BLAS routine xGEMM to transform the right-hand side C, calling xTRSYL to solve the (quasi-) triangularSylvester equation, and using xGEMM to transform back to the solution X.The error bound (15.29) can be estimated using xLACON in conjunction withthe above routines.
A Fortran 77 code dggsvx [556, 1993] of Higham followsthis outline and may appear in a future release of LAPACK.Routine xLASY2 solves a real Sylvester equation AX ± XB = σC in whichA and B have dimension 1 or 2 and σ is a scale factor. It is called by xTRSYL .Kågström and Poromaa have developed codes for solving (15.31), whichare intended for a future release of LAPACK [622, 1993], [623, 1994].Problems15.1.
Show that the Sylvester equation AX – XA = I has no solution.15.2. (Bellman [89,1970,§10.18]) Show that if the expressionexists for all C it represents the unique solution of the Sylvester equationAX + XB = C. (Hint: consider the matrix differential equation dZ/dt =AZ(t) + Z(t)B, Z(0) = C.) Deduce that the Lyapunov equation AX +XAT = –C has a symmetric positive definite solution if A has eigenvalueswith negative real parts and C is symmetric positive definite.15.3. (Byers and Nash [176,1987])Letand considerShow that there exists a minimizer X that is either symmetric or skewsymmetric.15.4. How would you solve a Sylvester equation AX – XB = C in which Aand B are of dimension 1 or 2? Compare your method with the one used inthe LAPACK routine XLASY2 .15.5.
(RESEARCH PROBLEM) Derive conditions for the Sylvester equation tohave a well-conditioned solution.Previous Home NextChapter 16Stationary Iterative MethodsI recommend this method to you for imitation.You will hardly ever again eliminate directly,at least not when you have more than 2 unknowns.The indirect [iterative] procedure can be done while half asleep,or while thinking about other things.17— CARL FRIEDRICH GAUSS, Letter to C. L.
Gerling (1823)The iterative method is commonly called the “Seidel process, ”or the “Gauss–Seidel process. ”But, as Ostrowski (1952) points out,Seidel (1874) mentions the process but advocates not using it.Gauss nowhere mentions it.— GEORGE E. FORSYTHE,Solving Linear Algebraic Equations Can Be Interesting (1953)The spurious contributions in null(A)grow at worst linearly andif the rounding errors are small the scheme can be quite effective.— HERBERT B. KELLER,On the Solution of Singular and SemidefiniteLinear Systems by Iteration (1965)17Gauss refers here to his relaxation method for solving the normal equations. Thetranslation is taken from Forsythe [387, 1951 ].325326S TATIONARY I TERATIVE M ETHODSTable 16.1.
Dates of publication of selected iterative methods. Based on Young [1123,19 8 9 ].1845187419101938–19391940s1952Jacobi methodJacobiGauss-Seidel methodSeidelRichardsonRichardson’s methodMethod of steepest descentTempleVarious (analysis by Successive overrelaxationYoung and Frankel) (SOR) methodHestenesand Stiefel Conjugate gradient methodIterative methods for solving linear systems have along history, going backat least to Gauss.
Table 16.1 shows the dates of publication of selected methods. It is perhaps surprising, then, that rounding error analysis for iterativemethods is not well developed. There are two main reasons for the paucityof error analysis. One is that in many applications accuracy requirementsare modest and are satisfied without difficulty, resulting in little demand forerror analysis.
Certainly there is no point in computing an answer to greateraccuracy than that determined by the data, and in scientific and engineeringapplications the data often has only a few correct digits. The second reasonis that rounding error analysis for iterative methods is inherently more difficult than for direct methods, and the bounds that are obtained are harder tointerpret.In this chapter we consider a simple but important class of iterative methods, stationary iterative methods, for which a reasonably comprehensive erroranalysis can be given.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
