Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 80
Текст из файла (страница 80)
No erroranalysis is given in [135, 1994], [178, 1992], or [865, 1991].The O(n2) complexity of the algorithms mentioned above for solving Vandermonde-like systems is not optimal. Lu [714, 1994], [715, 1995], [716, 1996]derives O(n log n log p) flops algorithms, where p is the number of distinctpoints. The numerical stability and practical efficiency of the algorithms remain to be determined. Bini and Pan [98, 1994] give an O ( n log2 n) algorithmfor solving a dual Vandermonde system that involves solving related Toeplitzand Hankel systems.Since Vandermonde systems can be solved in less than O(n3) flops it isnatural to ask whether the O(mn 2 ) complexity of QR factorization of anm x n matrix can be bettered for a Vandermonde matrix.
QR factorizationalgorithms with O(mn) flop counts have been derived by Demeure [277, 1989],[278, 1990], and for Vandermonde-like matrices where the polynomials satisfya three-term recurrence by Reichel [864, 1991]. No error analysis has beenpublished for these algorithms. Demeure’s algorithms are likely to be unstable,because they form the normal equations matrix VTV.Problems21.1. Derive a modified version of Algorithm 21.1 in which the scale factor442V ANDERMONDE S YSTEMSs = q(αi ) is computed directly asWhat is the flop count for this version?21.2. (Calvetti and Reichel [179, 1993]) Generalize Algorithm 21.1 to theinversion of a Vandermonde-like matrix for polynomials that satisfy a threeterm recurrence relation.21.3.
Investigate the stability of Algorithm 21.1 and the modified versionof Problem 21.1. (a) Evaluate the left and right residuals of the computedinverses; compare the results with those for GEPP. (b) Show that Algorithm 21.1 always performs subtractions of like-signed numbers and explainthe implications for stability. (Does Algorithm 21.1 share the high accuracyproperties discussed at the end of 521.3.1?) (c) (RESEARCH PROBLEM) Deriveand explore forward error bounds and residual bounds for both algorithms.Extend the analysis to the algorithms of Calvetti and Reichel [179, 1993].21.4. By summing (21.1) for i = 0:n, show thatWhat doesthis imply about the sign pattern of V–1 ? What is the sum of all the elementso f V –l ?be a Chebyshev–Vandermonde21.5.
Let T = T ( α0, α1, . . . , α n ) =matrix (Tj is the Chebyshev polynomial of degree j), with T–1 = (uij).Analogously to (21.1) we haveHencewhere V-l = V(α0, α1, . . . ,αn )-l = (wij). Show thatand hence that(Hint: show that T = LV for a lower triangular matrix L and use the factthat21.6. Show that for a nonconfluent Vandermonde-like matrix P =where the pi satisfy (21.6),(Hint: study the structure of (21.22 ).)PROBLEMS44321.7. Show that for the Chebyshev-Vandermonde-like matrix T = T(α 0 ,α 1 , . .
. , α n ),1.K 2 (T)2.K2 (T)=for αi = cos((i + ½)π/(n + 1)) (zeros of Tn+1 ).< 2, for αi = cos(iπ/n ) (extrema of Tn ).(Hint: use the discrete orthogonality properties of the Chebyshev polynomials;see, e.g., Hamming [501, 1973, pp.472–473].)21.8. Derive an O(n2 ) flops algorithm that reorders the distinct points α0 ,α1, . . . .
αnaccording to the same permutation that would be produced byGEPP applied to PT(α0, α1,. . . . , αn ). (Cf. Problem 5.4. ) Can your algorithmever produce the increasing ordering?21.9. Derive Algorithm 21.8 by differentiating repeatedly the original Clenshaw recurrence (which is Algorithm 21.8 with k = 0) and resealing so as toconsign factorial terms to a cleanup loop at the end.
Derive an algorithm forcomputing the residual for the primal system in a similar way, using recurrences obtained by differentiating (21.6).21.10. (Higham [533, 1987]) A structured condition number for the primalVandermonde system Vx = b, where V = V(α0, α1, . . . .
,α n ), can be definedbyShow thatcond(V) =and derive a corresponding condition number for the dual system VTa = f.Previous HomeChapter 22Fast Matrix MultiplicationA simple but extremely valuable bit of equipment in matrix multiplicationconsists of two plain cards,with a re-entrant right angle cut out of one or both of themif symmetric matrices are to be multiplied.In getting the element of the ith row and jth column of the product,the ith row of the first factor and the jth column of the secondshould be marked by a card beside, above, or below it.— HAROLD HOTELLING, Some New Methods in Matrix Calculation (1943)It was found that multiplication of matrices using punched card storagecould be a highly efficient process on the Pilot ACE,due to the relative speeds of the Hollerith card reader used for input(one number per 16 ins.) and the automatic multiplier (2 ins.).While a few rows of one matrix were held in the machinethe matrix to be multiplied by it was passed through the card reader.The actual computing and selection of numbers from storeoccupied most of the time between the passage ofsuccessive rows of the cards through the reader,so that the overall time was but little longerthan it would have been if the machinehad been able to accommodate both matrices.— MICHAEL WOODGER, The History and Present Use ofDigital Computers at the National Physical Laboratory (1958)445Next446FAST M ATRIX M ULTIPLICATION22.1.
MethodsA fast matrix multiplication method forms the product of two n x n matricesinarithmetic operations, where ω < 3. Such a method is more efficientasymptotically than direct use of the definition(22.1)which requires O(n3) operations. For over a century after the development ofmatrix algebra in the 1850s by Cayley, Sylvester and others, this definitionprovided the only known method for multiplying matrices. In 1967, however,to the surprise of many, Winograd found a way to exchange half the multiplications for additions in the basic formula [1105, 1968].
The method restson the identity, for vectors of even dimension n,(22.2)When this identity is applied to a matrix product AB, with x a row of A andy a column of B, the second and third summations are found to be commonto the other inner products involving that row or column, so they can becomputed once and reused. Winograd’s paper generated immediate practicalinterest because on the computers of the 1960s floating point multiplicationwas typically two or three times slower than floating point addition. (Ontodays’ machines these two operations are usually similar in cost).Shortly after Winograd’s discovery, Strassen astounded the computer science community by finding anoperations method for matrix multiplication (log2 7 2.807).
A variant of this technique can be used to computeA-l (see Problem 22.8) and thereby to solve AX = b, both inoperations. Hence the title of Strassen’s 1969 paper [962, 1969], which refers tothe question of whether Gaussian elimination is asymptotically optimal forsolving linear systems.Strassen’s method is based on a circuitous way to form the product of apair of 2 x 2 matrices in 7 multiplications and 18 additions, instead of the usual8 multiplications and 4 additions. As a means of multiplying 2 x 2 matrices theformulae have nothing to recommend them, but they are valid more generallyfor block 2 x 2 matrices. Let A and B be matrices of dimensions m x n andn x p respectively, where all the dimensions are even, and partition each of A,B, and C = AB into four equally sized blocks:(22.3)22.1 M ETHODS447Strassen’s formulae areP1 = (A ll + A 22 )(B II + B 22 )P2 = (A21 + A22 )B11 ,P3 = A11 (B12 – B22 ),P4P5P6P7====A22 (B21 – B11 ),(A11 + A12 )B22 ,(A 21 – A II )(B II + B 12 )(A12 – A22 )(B21 + B22 )(22.4)C11 = P1 + P4 – P5 + P7,C12 = P3 + P5,C21 = P2 + P4,C22 = P1 + P3 – P2 + P6.Counting the additions (A) and multiplications (M) we find that while conventional multiplication requiresmnpM + m(n – 1)pA,Strassen’s algorithm, using conventional multiplication at the block level, requiresThus, if m, n, and p are large, Strassen’s algorithm reduces the arithmeticby a factor of about 7/8.
The same idea can be used recursively on themultiplications associated with the Pa. In practice, recursion is only performeddown to the “crossover” level at which any savings in floating point operationsare outweighed by the overheads of a computer implementation.To state a complete operation count, we suppose that m = n = p = 2 kand that recursion is terminated when the matrices are of dimension no = 2 r ,at which point conventional multiplication is used. The number of multiplications and additions can be shown to beM(k) = 7 k - r8r ,(22.5)The sum M(k) + A(k) is minimized over all integers r by r = 3; interestingly,this value is independent of k. The total operation count for the “optimal”no = 8 is less thanHence, in addition to having a lower exponent, Strassen’s method has a reasonable constant.448FAST M ATRIX M ULTIPLICATIONWinograd found a variant of Strassen’s formulae that requires the samenumber of multiplications but only 15 additions (instead of 18).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
