Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 84
Текст из файла (страница 84)
VAN LOAN, Computation/Frameworks for the Fast Fourier Transform (1992)465466T HE FAST FOURIER T RANSFORMANDA PPLICATIONS23.1. The Fast Fourier TransformThe matrix–vector product y = Fn x, whereis the key computation in the numerical evaluation of Fourier transforms.
Ifthe product is formed in the obvious way then O(n2) operations are required.The fast Fourier transform (FFT) is a way to compute y in just O ( n log n)operations. This represents a dramatic reduction in complexity.The FFT is best understood (at least by a numerical analyst!) by interpreting it as the application of a clever factorization of the discrete Fouriertransform (DFT) matrix Fn .Theorem 23.1 (Cooley–Tukey radix 2 factorization).
If n = 2 t then the DFTmatrix Fn may be factorized asFn = At . . . A1Pn,(23.1)where Pn is a permutation matrix andProof. See Van Loan [1044, 1992, Thm. 1.3.3].The theorem shows that we can write y = Fnx asy = At . . . A1Pnx.which is formed as a sequence of matrix-vector products.
It is the sparsity ofthe Ak (two nonzeros per row) that yields the O ( n log n) operation count.We will not consider the implementation of the FFT, and therefore wedo not need to define the “bit reversing: permutation matrix Pn in (23.1).However, the way in which the weightsare computed does affect the accuracy. We will assume that computed weightsare used that satisfy, for allj and k,(23.2)Among the many methods for computing the weights are ones for which wecan take µ = cu, µ = cu log j, and µ = cuj, where c is a constant that dependson the method; see Van Loan [1044, 1992, 1.4].We are now ready to prove an error bound.46723.1 T HE F AST F OURIER T RANSFORMTheorem 23.2.
Let = fl(F n x), where n = 2t, be computed using theCooley-Tukey radix 2 FFT, and assume that (23.2) holds. ThenProof. Denote byThenthe matrix defined in terms of the computedusing the fact that each Ak has only two nonzeros per row, and recalling thatwe are using complex arithmetic. In view of (23.2),Hence, overall,Invoking Lemma 3.7 we find thatusing Lemma 3.1 for the second inequality.
NowHence(23.3)Finally, because Fn istimes a unitary matrixTheorem 23.2 says that the Cooley–Tukey radix 2 FFT yields a tiny forward error, provided that the weights are computed stably. It follows immediately that the computation is backward stable, since= y + ∆y = Fnx + ∆yimplies= F n ( x + ∆x) with ||∆x||2/||x||2 = ||∆y||2/||y||2. If we formy = Fn x by conventional multiplication using the exact Fn , then (ProbHence whenlem 3.7)soµ is of order u, the FFT has an error bound smaller than that for conventional multiplication by the same factor as the reduction in complexity of themethod.
To sum up, the FFT is perfectly stable.468T HE FAST FOURIER T RANSFORMANDA PPLICATIONSFigure 23.1. Error in FFT followed by inverse FFT (“o”). Dotted line is errorbound.The inverse transform x = F n - 1 y = n – 1 F n *y can again be computedin O(n log n) operations using the Cooley–Tukey radix 2 factorization, andsatisfies the same bound as in Theorem 23.2. (Strictly, weshould replace t by t + 1 in the bound to account for the rounding error individing by n)For other variations of the FFT, based on different radices or differentfactorization of Fn , results analogous to Theorem 23.2 hold.A simple way to test the error bounds is to compute the FFT followedby the inverse transform,and to evaluate e n = ||x –Our analysis gives the bound en < 2n1/2 log2 nη + O(η 2).
Figure 23.1 plots en and the approximate error bound n1/2 log2 nu for n = 2k ,k = 0:16, with random x from the normal N(0,1) distribution (the FFTswere computed using MATLAB’s fft and ifft functions). The errors grow atroughly the same rate as the bound and are on average about a factor of 10smaller than the bound.23.2. Circulant Linear SystemsA circulant matrix (or circulant, for short) is a special Toeplitz matrix inwhich the diagonals wrap around:23.2 C IRCULANT LINEAR S YSTEMS469Circulant matrices have the important property that they are diagonalized bythe DFT matrix Fn :Fn CFn -l = D = diag(d i).Moreover, the eigenvalues are given by d = Fnc. Hence a linear system Cx = bcan be solved in O(n log n) operations with the aid of the FFT as follows:(1)(2)(3)(4)d = Fnc,g = Fnb,h = D -l g,x = Fn-1h.The computation involves two FFTs, a diagonal scaling, and an inverse FFT.We now analyse the effect of rounding errors.
It is convenient to write theresult of Theorem 23.2 in the equivalent form (from (23.3))(23.4)Steps (1) and (2) yield= (Fn + ∆1 )c,= (Fn + ∆2 )b,||∆1||2 < f(n, u),||∆2||2 < f(n, u).(23.5)Equation (23.5) implies that= D + ∆ 0 , ||∆ 0 || 2 < f(n, u)||C||2.For steps (3) and (4) we have, using Lemma 3.5,Putting these equations together we have(23.6)or, working to first order,We obtain a backward stability result.470T HE FAST FOURIER T RANSFORMANDA PPLICATIONSTheorem 23.3. Let Cbe a circulant and suppose the system Cx = bis solved by the FFT process described above, where the FFT satisfies (23.4).Then the computed solutionsatisfies (C + ∆C ) = b, whereThe conclusion is that the FFT method for solving a circulant system isnormwise backward stable provided that the FFT itself is computed stably.Although we have shown that solves a slightly perturbed system, theperturbed matrix C + ∆C is not a circulant.
In fact, does not, in general,solve a nearby circulant system, as can be verified experimentally by computing the "(circulant backward error” using techniques from [527, 1992]. Thebasic reason for the lack of this stronger form of stability is that there are toofew parameters in the matrix onto which to throw the backward error.A forward error bound can be obtained by applying standard perturbation theory to Theorem 23.3: the forward error is bounded by a multiple ofk2 (C)u. That the forward error can be as large as k2 (C)u is clear from theanalysis above, because (23.5) shows that the computed eigenvalueis contaminated by an error of order u23.3. Notes and ReferencesFor a unified treatment of the many different versions of the FFT, includingimplementation details, see Van Loan [1044, 1992].For a comprehensive survey of the discrete Fourier transform see Briggsand Henson [148, 1995].The Cooley–Tukey radix 2 FFT was presented in [240, 1965], which is oneof the most cited mathematics papers of all time [554, 1993, p.
171].The history of the FFT is discussed by Cooley [238, 1990], [239, 1994]and Cooley and Tukey [241, 1993]. Cooley [238, 1990] states that the earliestknown reference to the FFT is an obscure 1866 paper of Gauss in neoclassicLatin, and he recommends that researchers not publish papers in neoclassicLatin!Theorem 23.2 is not new, but the proof is more concise than most inthe literature. Early references that derive error bounds using the matrixfactorization formulation of the FFT are Gentleman and Sande [434, 1966]and Ramos [859, 1971]. A full list of references for error analysis of the FFTis given by Van Loan [1044, 1992, 1.4].Linzer [708, 1992] shows that the FFT-based circulant solver is forwardstable and poses the question of whether or not the solver is backward stable.Our backward error analysis answers this question positively and thereforealso proves forward stability.PROBLEMS471One application of circulant linear systems is in the preconditioned conjugate gradient method for solving Toeplitz systems.
The idea of using acirculant preconditioned was suggested by Strang [960, 1986], and the theoryand practice of this technique is now well developed. For more details seeChan, Nagy, and Plemmons [191, 1994] and the references therein. A goodsource of results about circulant matrices is the book by Davis [266, 1979].Problems23.1. (Bailey [44, 1993]) In high-precision multiplication we have two integern -vectors x and y representing high-precision numbers and we wish to formthe termsBy padding the vectors with n zeros, these products can be expressed in the formwherek + 1 – j is interpreted as k + 1 – j + n if k + 1 – j is negative.
These products represent a convolution: a matrix–vector product involving a circulantmatrix. Analogously to the linear system solution in 23.2, this product canbe evaluated in terms of discrete Fourier transforms as z =where the dot denotes componentwise multiplication of two vectors. Since xand y are integer vectors, the zi should also be integers, but in practice theywill be contaminated by rounding errors.
Obtain a bound on z – and deducea sufficient condition for the “nearest integer vector to to be the exact z.Previous HomeChapter 24Automatic Error AnalysisGiven the pace of technology,I propose we leave math to the machines and go play outside.— CALVIN, Calvin and Hobbes by Bill Watterson (1992)To analyse a given numerical algorithm we proceed as follows.A number which measures the effect of roundoff erroris assigned to each set of data.“Hill-climbing” procedures are then applied to search forvalues large enough to signal instability.— WEBB MILLER, Software for Roundoff Analysis (1975)473Next474A UTOMATIC E RROR A NALYSISAutomatic error analysis is any process in which we use the computer to helpus analyse the accuracy or stability of a numerical computation. The ideaof automatic error analysis goes back to the dawn of scientific computing.For example, running error analysis, described in 3.3, is a form of automaticerror analysis; it was used extensively by Wilkinson on the ACE machine.Various forms of automatic error analysis have been developed.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
