Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 82
Текст из файла (страница 82)
Assuming that these four matrices are added inthe order indicated, we haveClearly, the same bound holds for the other three ||∆Cij|| terms. Thus, overall,= AB + ∆C ,||∆C|| < (46 · 2 k-1 + 12ck- l )u||A|| ||B|| + O(u2).A comparison with (22.15) shows that we need to define the ck byc k = 12ck-1 +46 · 2 k-1, k > r, cr = 4 r ,where cr =(22.16)Solving this recurrence we obtainwhich gives (22.14).The forward error bound for Strassen’s method is not of the componentwiseform (22.10) that holds for conventional multiplication, which we know itcannot be by Miller’s result.
One unfortunate consequence is that while thescaling AB(AD)(D-1 B), where D is diagonal, leaves (22.10) unchanged,it can alter (22.14) by an arbitrary amount.The reason for the scale dependence is that Strassen’s method adds together elements of A matrix-wide (and similarly for B); for example, in (22.4)A11 is added to A22 , A12 , and A21 . This intermingling of elements is particularly undesirable when A or B has elements of widely differing magnitudesbecause it allows large errors to contaminate small components of the product.This phenomenon is well illustrated by the examplewhich is evaluated exactly in floating point arithmetic if we use conventionalmultiplication.
However, Strassen’s method computes22.2 E RROR ANALYSIS455Because c22 involves subterms of order unity, the error c22 –will be oforder u. Thus the relative error |c22 –|/|c22| =which is muchlarger than u if ε is small. This is an example where Strassen’s method doesnot satisfy the bound (22.10). For another example, consider the productX = P32 E, where Pn is the n x n Pascal matrix (see 26.4) and eij = 1/3.With just one level of recursion in Strassen’s method we find in MATLAB thatis of order 10-5, so that, again, some elements of thecomputed product have high relative error.It is instructive to compare the bound (22.14) for Strassen’s method withthe weakened, normwise version of (22.10) for conventional multiplication:(22.17)The bounds (22.
14) and (22. 17) differ only in the constant term. For Strassen’smethod, the greater the depth of recursion the bigger the constant in (22.14):if we use just one level of recursion (n0 = n/2) then the constant is 3n2 +25n, whereas with full recursion (n0 = 1) the constant is6n3.585 – 5n. It is also interesting to note that the bound for Strassen’s method(minimal for n 0 = n) is not correlated with the operation count (minimal forn0 = 8).Our conclusion is that Strassen’s method has less favorable stability properties than conventional multiplication in two respects: it satisfies a weakererror bound (normwise rather than componentwise) and it has a larger constant in the bound (how much larger depending on no).Another interesting property of Strassen’s method is that it always involvessome genuine subtractions (assuming that all additions are of nonzero terms).This is easily deduced from the formulae (22.4).
This makes Strassen’s methodunattractive in applications where all the elements of A and B are nonnegative(for example, in Markov processes). Here, conventional multiplication yieldslow componentwise relative error because, in (22.10), |A||B| = |AB| = |C|,yet comparable accuracy cannot be guaranteed for Strassen’s method.An analogue of Theorem 22.2 holds for Winograd’s variant of Strassen’smethod.Theorem 22.3. Let A, Bwhere n = 2 k. Suppose that C = AB iscomputed by Winograd’s variant (22.6) of Strassen’s method and that n0 = 2 ris the threshold at which conventional multiplication is used. The computedproductsatisfies(22.18)Proof. The proof is analogous to that of Theorem 22.2, but more tedious.It suffices to analyse the computation of C12, and the recurrence correspondingFAST M ATRIX M ULTIPLICATION456to (22.16) isck= 1 8c k - l + 8 9 · 2 k – 1 , k > r, c r = 4 r .Note that the bound for the Winograd–Strassen method has exponentlog2 18 4.170 in place of log2 123.585 for Strassen’s method, suggestingthat the price to be paid for a reduction in the number of additions is anincreased rate of error growth.
All the comments above about Strassen’smethod apply to the Winograd variant.Two further questions are suggested by the error analysis:llHow do the actual errors compare with the bounds?Which formulae are the more accurate in practice, Strassen’s or Winograd’s variant?To give some insight we quote results obtained with a single precision Fortran 90 implementation of the two methods (the code is easy to write if weexploit the language’s dynamic arrays and recursive procedures). We takerandom n x n matrices A and B and look at ||AB – fl(AB)||/(||A|| ||B||) forn0 = l, 2, .
. . , 2k = n (note that this is not the relative error, since the denominator is ||A|| ||B|| instead of ||AB||, and note that no = n correspondsto conventional multiplication). Figure 22.2 plots the results for one randommatrix of order 1024 from the uniform [0, 1] distribution and another matrix ofthe same size from the uniform [– 1, 1] distribution. The error bound (22.14)for Strassen’s method is also plotted. Two observations can be made.lWinograd’s variant can be more accurate than Strassen’s formulae, forall no, despite its larger error bound.lThe error bound overestimates the actual error by a factor up to 1.8 x 106for n = 1024, but the variation of the errors with no is roughly aspredicted by the bound.22.2.3.
Bilinear Noncommutative AlgorithmsBini and Lotti [97, 1980] have analysed the stability of bilinear noncommutative algorithms in general. They prove the following result.Theorem 22.4 (Bini and Lotti). Let A, B(n = h k ) and let theproduct C = AB be formed by a recursive application of the bilinear noncommutative algorithm (22.7), which multiplies h x h matrices using t nonscalarmultiplications. The computed productsatisfies(22.19)22.2 E RROR A NALYSIS457Figure 22.2.
Errors for Strassen’s method with two random matrices of dimensionn = 1024. Strassen’s formulae: “x”, Winograd’s variant: "o". X-axis contains log2of recursion threshold n0, 1 < n0 < n. Dot-dash line is error bound for Strassen’sformulae.where α and β are constants that depend on the number of nonzero terms inthe matrices U, V and W that define the algorithm.The precise definition of α and β is given in [97, 1980]. If we take k = 1,so that h = n, and if the basic algorithm (22.7) is chosen to be conventionalmultiplication, then it turns out that α = n – 1 and β = n, so the boundof the theorem becomes (n – 1)nu||A|| ||B|| + O(u2), which is essentially thesame as (22.17). For Strassen’s method, h = 2 and t = 7, and α = 5, β = 12,whichso the theorem produces the boundis a factor log2 n larger than (22.14) (with n0 = 1).
This extra weakness ofthe bound is not surprising given the generality of Theorem 22.4.Bini and Lotti consider the set of all bilinear noncommutative algorithmsthat form 2 x 2 products in 7 multiplications and that employ integer constantsof the form ±2i , where i is an integer (this set breaks into 26 equivalenceclasses). They show that Strassen’s method has the minimum exponent βin its error bound in this class (namely, β = 12). In particular, Winograd'svariant of Strassen’s method has β = 18, so Bini and Lotti’s bound has thesame exponent log2 18 as in Theorem 22.3.FAST M ATRIX M ULTIPLICATION45822.2.4.
The 3M MethodA simple example reveals a fundamental weakness of the 3M method. Consider the computation of the scalarIn floating point arithmetic, if y is computed in the usual way, as y = θ(1/θ ) +( 1 /θ) θ, then no cancellation occurs and the computed has high relativeThe 3M method computesaccuracy:If |θ| is large this formula expresses a number of order 1 as the differencewill almost certainly be contaminatedof large numbers. The computedby rounding errors of order uθ2, in which case the relative error is large:However, if We measure the error inrelative to z,then it is acceptably small:This example suggests thatthe 3M method may be stable, but in a weaker sense than for conventionalmultiplication.To analyse the general case, consider the product C 1 + iC 2 = (A1 +iA2 )( B1 + iB2 ), where Ak, Bk, Ckk = 1:2.
Using (22.10) we findthat the computed product from conventional multiplication,satisfies(22.20)(22.21)For the 3M method C1 is computed in the conventional way, and so (22.20)holds. It is straightforward to show thatsatisfies(22.22)Two notable features of the bound (22.22) are as follows. First, it is ofa different and weaker form than (22.21); in fact, it exceeds the sum of thebounds (22.20) and (22.21). Second and more pleasing, it retains the propertyof (22.20) and (22.21) of being invariant under diagonal scalingsC = ABD 1 AD 2 · D 2 -1 BD 3 = D 1 CD 3 ,D j diagonal,in the sense that the upper bound ∆C2 in (22.22) scales also according toD 1 ∆C2D3.
(The “hidden” second-order terms in (22.20)–(22.22) are invariantunder these diagonal scalings. )22.3 N OTESANDR EFERENCES459The disparity between (22.21) and (22.22) is, in part, a consequence ofthe differing numerical cancellation properties of the two methods. It is easyto show that there are always subtractions of like-signed numbers in the 3Mmethod, whereas if A1, A2, Bl, and B2 have nonnegative elements (for example) then no numerical cancellation takes place in conventional multiplication.We can define a measure of stability with respect to which the 3M methodmatches conventional multiplication by taking norms in (22.21) and (22.22).We obtain the weaker bounds(22.23)(22.24)(having used || |A1| + |A2| || <||Al + iA2 ||).
Combining these with an analogous weakening of (22.20), we find that for both conventional multiplicationand the 3M method the computed complex matrix satisfieswhere cn = O(n).The conclusion is that the 3M method produces a computed productwhose imaginary part may be contaminated by relative errors much largerthan those for conventional multiplication (or, equivalently, much larger thancan be accounted for by small componentwise perturbations in the data Aand B).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
