Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 60
Текст из файла (страница 60)
The computed solutionagain satisfies (15.9).The use of iterative methods to solve the Sylvester equation has attractedattention recently for applications where A and B are large and sparse [572,1995 ], [588, 1992 ], [935, 1991 ], [1059, 19 88]. The iterations are usually terminated when an inequality of the form (15.9) holds, so here the size of therelative residual is known a priori (assuming the method converges).15.2. Backward ErrorWe saw in the last section that standard methods for solving the Sylvesterequation are guaranteed to produce a small relative residual.
Does a smallrelative residual imply a small backward error? The answer to this questionfor a general linear system is yes (Theorem 7. 1). But for the highly structuredSylvester equation the answer must be no, because for the special case of matrix inversion we know that a small residual does not imply a small backwarderror (§13.1). In this section we investigate the relationship between residualand backward error for the Sylvester equation.The normwise backward error of an approximate solution Y to (15.1) isdefined by(15.10)The tolerances a, β, and γ provide some freedom in how we measure theperturbations. Of most interest is the choice a = ||A||F, β = ||B||F, γ =||C||F , for which we will call q the normwise relative backward error.
Theequation (A + ∆A)Y – Y(B + ∆B) = C + ∆C may be written(15.11)where the residual R = C – (AY – YB). A small backward error implies asmall relative residual since, using the optimal perturbations from (15. 10) in(15.11), we have(15.12)To explore the converse question of what the residual implies about thebackward error we begin by transforming (15.11) using the SVD of Y, Y =UΣ VT, whereandare orthogonal and Σ = diag(σi ) 1The numbers σ1 > σ 2 > . . . > σ m i n (m , n) > 0 are the singular values314T HE S YLVESTER E QUATIONof Y and we define, in addition, σ m i n (m , n)+l = . .
. = σ m a x (m , n) = 0. Equation(15.1 1) transforms to(15.13)whereThis is an underdetermined system, with mn equations in m 2 + n 2 + mnunknowns. We can write it in the uncoupled form16(15.14)For each i and j it is straightforward to show that the minimum ofsubject to (15.14) is attained forThese matrices minimizeSince η(Y ) is the minimum value of max{||a - 1 ∆ A||F,||β - 1∆ B||F,||γ - 1 ∆C||F} ,it follows that.(15.15)where(15.16)16For notational convenience we extend(if m < n) or(if m > n) to dimensionm × n; the “fictitious” elements will be set to zero by the minimization.15.2 B ACKWARD E RROR315This expression shows that the backward error is approximately equal not tothe normwise relative residualbut to a componentwise residual corresponding to the diagonalized equation (15.13).From (15.15) and (15.16) we deduce that(15.17)where(15.18)The scalar µ > 1 is an amplification factor that measures by how much, atworst, the backward error can exceed the normwise relative residual.
We nowexamine µ more closely, concentrating on the normwise relative backwarderror, for which a = ||A||F, β = ||B||F, and γ = ||C||F.First, note that if n = 1 and B = 0, so that the Sylvester equationreduces to a linear system Ay = c, then σ 1 = ||y||2 and σ k = 0 for k > 1, soand sowe recover Theorem 7.1 (for the 2-norm) from (15.12) and (15.17), to withina constant factor.If m = n then(15.19)We see that µ is large only when(15.20)that is, when Y is ill conditioned and Y is a large-normed solution to theSylvester equation. In the general case, with mn, one ofandisalways zero and hence µ can be large for a third reason: A (if m < n ) or B(if m > n) greatly exceeds the rest of the data in norm; in these cases theSylvester equation is badly scaled.
However, if we set a = β = ||A||F + ||B||F,which corresponds to regarding A and B as comprising a single set of data,then bad scaling does not affect µ.If we allow only A and B to be perturbed in (15.10) (as may be desirableif the right-hand side C is known exactly), then γ = 0 and (15.19) and (15.20)remain valid with ||C||F replaced by zero. In this case µ > ||Y||F||Y+|| 2κ 2 (Y ) (for any m and n), so µ is large whenever Y is ill conditioned (andincluded in this case is matrix inversion).
Conditions involving controllabilitywhich guarantee that the solution to the Sylvester equation with m = n isnonsingular are given by Hearon [508, 1977], while Datta [265, 1988] gives adeterminantal condition for nonsingularity. It is an open problem to deriveT HE S YLVESTER E QUATION316conditions for the Sylvester equation to have a well-conditioned solution (seeProblem 15.5).The following numerical example illustrates the above analysis. This particular example was carefully chosen so that the entries of A and B are of asimple form, but equally effective examples are easily generated using random,ill-conditioned A and B of dimension m, n > 2. LetDefine C by the property that vet(C ) is the singular vector corresponding tothe smallest singular value of InA – BT Im.
With a = 10-6, we solved theSylvester equation in MATLAB by the Bartels-Stewart algorithm and foundthat the computed solution X satisfiesAlthoughhas a very acceptable residual (as it must in view of (15.9)), itsbackward error is eight orders of magnitude larger than is necessary to achievebackward stability. We solved the same Sylvester equation using GEPP on thesystem (15.2). The relative residual was again less than u, but the backwarderror was appreciably larger:One conclusion we can draw from the analysis is that standard methods forsolving the Sylvester equation are at best conditionally backward stable, sincethere exist rounding errors such thatis the only nonzero element ofand then (15.17) is an approximate equality, with µ possibly large.15.2.1. The Lyapunov EquationIf we put B = –AT in the Sylvester equation we obtainAX+ XA T = C,which is called the Lyapunov equation.
This equation plays a major role incontrol and systems theory and it can be solved using the same techniques asfor the Sylvester equation.If C = C T then C = AX + XA T = X T A T + AX T = C T , so X andTX are both solutions to the Lyapunov equation. If the Lyapunov equationis nonsingular (equivalently,for all i and j, by (15.3)) ittherefore has a unique symmetric solution.31715.2 B ACKWARD E RRORWe assume that C is symmetric and that Y is a symmetric approximatesolution. The definition of backward error is nowThe analogue of (15.11) isLet Y = UΛUT be a spectral decomposition, with A = diagresidual equation transforms towherewritten in uncoupled form asandThen theThis system can be(15.21)We can obtain the minimum value ofby minimizingsubject to (15.21), for i, j = 1:n.
The solution is(Note thatwhereis symmetric sinceis.) It follows thatwhere the last inequality is usually a good approximation. Comparing with(15.16) we see that respecting the extra structure of the Lyapunov equationhas essentially no effect on the backward error.Finally, the analogue of (15.17) and (15.18) iswhereT HE S YLVESTER E QUATION31815.3. Perturbation ResultTo derive a perturbation result we consider the perturbed Sylvester equationwhich, on dropping second-order terms, becomesThis system may be written in the form(15.22)where P = InA – BTI m . If we measure the perturbations normwise bywhere a, β, and γ are tolerances as in (15.10), then(15.23)is a sharp bound (to first order in), where(15.24)is the corresponding condition number for the Sylvester equation.
The bound(15.23) can be weakened to(15.25)whereIfthen twice the upper bound in (15.25) can be shownto be a strict bound for the error. The perturbation bound (15.25) witha ||A||F, β = ||B||F, and γ = ||C||F is the one that is usually quoted inthe literature for the Sylvester equation (see [464, 1979] and [522, 1988], forexample), and corresponds to applying standard perturbation theory for Ax =b to (15.2). Note that ||P- 1||2 = sep(A, B)-1, where sep is the separation ofA and B,(15.26)15.3 PERTURBATION RESULT319The sep function is an important tool for measuring invariant subspace sensitivity [470, 1989, §7.2.5], [940, 1973], [1050, 1979].For the Lyapunov equation, a similar derivation to the one above showsthat the condition number is(15.27)where Π is the vet-permutation matrix, which is defined by the property thatvec(A T) = Π vet(A).How much can the bounds (15.23) and (15.25) differ? The answer is byan arbitrary factor.
To show this we consider the case where B is normal(or equivalently, A is normal if we transpose the Sylvester equation). Wecan assume B is in Schur form, thus B = diag(µ j ) (with the µ j possiblycomplex). Then P = diag(A – ujjIm), and it is straightforward to show thatif X = [x1, . . . , xn ], and if we approximate the 2-norms in the definitions ofand Φ by Frobenius norms, thenwhileThese formulae show that in general Ψ and Φ will be of similar magnitude,and we know that Ψ < Φ from the definitions.
However, Ψ can be muchsmaller than Φ. For example, suppose that γ = 0 andThen ifwe have Ψ << Φ. Such examples are easily constructed. To illustrate, letA = diag(2, 2, . . . , 2, 1) and B = diag(1/2, 1/2, . . . , 1/2, 1 –with> 0, soand let X = (A – µ n n I m )Y,that A – µ nnIm = diagwhere Y = [y, y, . . . , y, 0] with ||(A–µ n n I m )y||2 = ||A–µ n n I m || 2 and ||y||2 = 1.Then, if γ =320T HE S YLVESTER E QUATIONTo summarize, the “traditional” perturbation bound (15.25) for the Sylvester equation can severely overestimate the effect of a perturbation on thedata when only A and B are perturbed, because it does not take accountof the special structure of the problem.
In contrast, the perturbation bound(15.23) does respect the Kronecker structure, and consequently is attainablefor any given A, B, and C.To obtain an a posteriori error bound for a computed solution∆ X we can set ∆A = 0, ∆B = 0, and ∆C = AX – XB – C = R in (15.22),which leads to(15.28)A similar but potentially much smaller bound is described in the next section.15.4.
Practical Error BoundsFor the Sylvester equation we can obtain an analogue of the practical errorbound (7.27) by identifying Ax = b with (15.2). For the computed residual ofa computed solution X we haveTherefore the bound is(15.29)where ||X|| := maxi,j |xij|. After transformation by the technique illustratedin (14.1), this bound can be estimated by the LAPACK norm estimator (Algorithm 14.4) at the cost of solving a few linear systems with coefficient matricesInA – B T Im and its transpose—in other words, solving a few Sylvesterequations AX – XB = C and A T X – XB T = D. If the Bartels-Stewartalgorithm is used, these solutions can be computed with the aid of the previously computed Schur decompositions of A and B. The condition numberΨ in (15.24) and sep(A, B) =can both be estimated in much thesame way; alternatively, the power method can be used (see Ghavimi andLaub [440, 1995]).
Other algorithms for efficiently estimating sep(A, B) givenSchur decompositions of A and B are given by Byers [172, 1984] and Kågströmand Poromaa [621, 1992].The attraction of (15.29) is that large elements in the j th column of P- 1may be countered by a small jth element of vec+ vec(R u ), making thebound much smaller than (15.28). In this sense (15.29) has better scaling32115.5 E XTENSIONSproperties than (15.28), although (15.29) is not actually invariant under diagonal scalings of the Sylvester equation.We give a numerical example to illustrate the advantage of (15.29) over(15.28). Letwheredenotes a Jordan block of size n with eigenvalueSolvingthe Sylvester equation by the Bartels–Stewart algorithm we found that thebounds are(where in evaluating (15.28) we replaced R by+ Ru, as in (15.29)).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
