Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 57
Текст из файла (страница 57)
The algorithm is limited to a minimum of two and amaximum of five iterations. Further, convergence is declared after computing ξ if the new ξ is the same as the previous one; this event signals thatconvergence will be obtained on the current iteration and that the next (andfinal) multiplication ATξ is unnecessary. Convergence is also declared if thenew ||y||1 is no larger than the previous one. This nonincrease of the normcan happen only in finite precision arithmetic and signals the possibility of avertex ej being revisited—the onset of “cycling.”The improved algorithm is as follows. This algorithm is the basis of allthe condition number estimation in LAPACK.Algorithm 14.4 (LAPACK norm estimator). Giventhis algorithm computes γ and υ = Aw such that γ < ||A|| 1 with ||υ||1 /||w|| 1 = γ (wis not returned).υ = A(n – le)if n = 1, quit with γ = |υ 1|, endγ = ||υ| |1ξ = sign(υ)296CONDITION NUMBER ESTIMATIONx = A Tξk=2repeatj = min{i: |xi | =υ = Aejγ = ||υ| |1if sign(υ) = ξ or γ <ξ = sign(υ)x = A Tξk = k+ 1}goto (*), enduntili+ 1(*) xi = (–1)x = Axi f 2||x||1 /(3n) > γ thenυ = xγ = 2||x||1 /(3n)endAlgorithm 14.4 can still be “defeated”: it returns an estimate 1 for matricesA(θ ) of the formA(θ) = I + θP, where P = PT, Pe = 0, Pe1 = 0, Pb = 0.(14.6)(P can be constructed as I – Q where Q is the orthogonal projection ontospan{e, e1, b}.) Indeed, the existence of counterexamples is intuitively obvious since Algorithm 14.4 samples the behaviour of A on less than n vectorsinNumerical counterexamples (not parametrized) can be constructedautomatically by direct search, as described in §24.3.1.
Despite these weaknesses, practical experience with Algorithm 14.4 shows that it is very rarefor the estimate to be more than three times smaller than the actual norm,independent of the dimension n. Therefore Algorithm 14.4 is, in practice, avery reliable norm estimator. The number of matrix-vector products requiredis at least 4 and at most 11, and averages between 4 and 5.There is an analogue of Algorithm 14.3 for complex matrices, in which ξiis defined as yi /|yi | if yi0 and 1 otherwise. In the corresponding version ofAlgorithm 14.4 the test for repeated ξ vectors is removed, because ξ now hasnoninteger, complex components and so is unlikely to repeat.It is interesting to look at a performance profile of Algorithm 14.4.
A performance profile is a plot of some measure of the performance of an algorithmversus a problem parameter. In this case, the natural measure of performanceis the underestimation ratio, γ/||A|| 1. Figure 14.1 shows the performanceprofile for a 5 × 5 matrix A( θ) of the form (14.6), with P constructed as14.4 O THER C ONDITION E STIMATORS297Figure 14.1. Underestimation ratio for Algorithm 14.4 for 5 × 5 matrix A(O) of (14.6)with 150 equally spaced values of[0,10].described above (because of rounding errors in constructing A(θ ) and withinthe algorithm, the computed norm estimates differ from those that would beproduced inexact arithmetic).
The jagged nature of the performance curve istypical for algorithms that contain logical tests and branches. Small changesin the parameter θ, which themselves result in different rounding errors, cancause the algorithm to visit different vertices in this example.14.4. Other Condition EstimatorsThe first condition estimator to be widely used is the one employed in LINPACK.
It was developed by Cline, Moler, Stewart, and Wilkinson [216, 1979].The idea behind this condition estimator originates with Gragg and Stewart [476, 1976], who were interested in computing an approximate null vectorrather than estimating the condition number itself.We will describe the algorithm as it applies to a triangular matrix TThere are three steps:1. Choose a vector d such that ||y|| is as large as possible relative to ||d||,where TTy = d.2. Solve Tx = y.CONDITION NUMBER ESTIMATION2983. EstimateIn LINPACK the norm is the l-norm, but the algorithm can also be usedfor the 2-norm or the cm-norm. The motivation for step 2 is based on a singularvalue decomposition analysis. Roughly, if ||y||/||d||is large thenwill almost certainly be at least as large, and it could bea much better estimate.
Notice that TTTx = d, so the algorithm is related tothe power method on the matrix (TTT) –1 with the specially chosen startingvector d.To examine step 1 more closely, suppose that T = UT is lower triangularand note that the equation Uy = d can be solved by the following columnoriented (saxpy) form of substitution:endThe idea is to choose the elements of the right-hand side vector d adaptivelyas the solution proceeds, with d j = ±1.
At the jth stage of the algorithmd . . . , dj+1 have been chosen and yn , . . . , yj+1 are known. The next elementdj{+1, –1} is chosen so as to maximize a weighted sum of dj – pj and thepartial sums p1, . . . , pj, which would be computed during the next executionof statement (*) above. Hence the algorithm looks ahead, trying to gaugethe effect of the choice of d j on future solution components. This heuristicalgorithm for choosing d is expressed in detail as follows.Algorithm 14.5 (LINPACK condition estimator). Given a nonsingular upper triangular matrixand a set of nonnegative weights {w i }, thisalgorithm computes a vector y such that Uy = d, where the elements dj = ±1are chosen to make ||y|| large.14.4 O THER C ONDITION E STIMATORS299p(1:j– 1 ) = p -(1:j–1)endendcost: 4n 2 flops.LINPACK takes the weights wj 1, though another possible (but moreexpensive) choice would be wj = 1/|ujj|, which corresponds to how pj isweighted in the expression yj = (dj – pj)/ujj.To estimate ||A- 1 || for a full matrix A, the LINPACK estimator makesuse of an LU factorization of A.
Given PA = LU, the equations solved areUTz = d, LTy = z, and AX = PTy, where for the first system d is constructedby the analogue of Algorithm 14.5 for lower triangular matrices; the estimateis ||x||1 /||y||1||A - 1 ||1. Since d is chosen without reference to L, there isan underlying assumption that any ill condition in A is reflected in U. Thisassumption may not be true; see Problem 14.3.In contrast to the LAPACK norm estimator, the LINPACK estimator requires explicit access to the elements of the matrix.
Hence the estimatorcannot be used to estimate componentwise condition numbers. Furthermore,separate code has to be written for each different type of matrix and factorization. Consequently, while LAPACK has just a single norm estimation routine,which is called by many other routines, LINPACK has multiple versions of itsalgorithm, each tailored to the specific matrix or factorization.Several years after the LINPACK condition estimator was developed, several parametrized counterexamples were found by Cline and Rew [217, 1983].Numerical counterexamples can also be constructed by direct search, as shownin §24.3.1. Despite the existence of these counterexamples the LINPACK estimator has been widely used and is regarded as being almost certain to producean estimate correct to within a factor 10 in practice.A 2-norm condition estimator was developed by Cline, Corm, and VanLoan [218, 1982, Algorithm 1]; see also Van Loan [1043, 1987] for anotherexplanation.
The algorithm builds on the ideas underlying the LINPACK estimator by using “look-behind” as well as look-ahead. It estimates σ min (R) =or σm a x (R) = ||R||2 for a triangular matrix R, where σ min and σ m a xdenote the smallest and largest singular values, respectively. Full matricescan be treated if a factorization A = QR is available ( Q orthogonal, R upper triangular), since R and A have the same singular values. The estimatorperforms extremely well in numerical tests, often producing an estimate thathas some correct digits [218, 1982], [534, 1987]. No counterexamples to theestimator were known until Bischof [103, 1990] obtained counterexamples asa by-product of the analysis of a different but related method, mentioned atthe end of this section.All the methods described so far have the property that when appliedrepeatedly to a given matrix they always produce the same estimate. Another300C ONDITION N UMBER E STIMATIONapproach is to introduce some randomness, so that the output of the methoddepends on the particular random numbers chosen.
A natural idea along theselines is to apply the power method to the matrix (AAT) -1 with a randomlychosen starting vector. If a factorization of A is available, the power methodvectors can be computed inexpensively by solving linear systems with A andAT. Analysis based on the singular value decomposition suggests that thereis a high probability that a good estimate of ||A - 1 ||2 will be obtained. Thisnotion is made precise by Dixon [306, 1983], who proves the following result.Theorem 14.6 (Dixon).
Letbe nonsingular and let θ > 1 be aconstant. Ifis a random vector from the uniform distribution on theunit sphere Sn =then the inequality(14.7)holds with probability at least 1 – 0.8θ - k / 2 n 1/2(k > 1).Note that the left-hand inequality in (14.7) always holds; it is only theright-hand inequality that is in question.For k = 1, (14.7) can be written aswhich suggests the simple estimatewhere x is chosen randomly from the uniform distribution on Sn. Such vectors x can be generatedfrom the formulawhere z1, .
. . , zn are independent random variables from the normal N(0, 1)distribution [668, 1981, p. 130]. If, for example, n = 100 and θ has the ratherlarge value 6400 then inequality (14.7) holds with probability at least 0.9.In order to take a smaller constant θ, for fixed n and a desired probability,we can use larger values of k. If k = 2j is even then we can simplify (14.7),obtaining(14.8)and the minimum probability stated by the theorem is 1 - 0.8θ - j n 1/2.
Takingj = 3, for the same value n = 100 as before, we find that (14.8) holds withprobability at least 0.9 for the considerably smaller value θ = 4.31.Probabilistic condition estimation has not yet been adopted in any major software packages, perhaps because the other techniques work so well.For more on the probabilistic power method approach see Dixon [306, 1983],Higham [534, 1987], and Kuczynski and Wozniakowski [676, 1992] (who alsoanalyse the more powerful Lanczos method with a random starting vector).For a probabilistic condition estimation method of very general applicability14.5 C ONDITION N UMBERSOFT RIDIAGONAL M ATRICES301see Kenney and Laub [652, 1994] and Gudmundsson, Kenney, and Laub [485,1995 ].The condition estimators described above assume that a single estimateis required for a matrix given in its entirety. Condition estimators have alsobeen developed for more specialized situations.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
