View Question
Q: Fourier transforms for reducing iterations ( Answered ,   0 Comments )
 Question
 Subject: Fourier transforms for reducing iterations Category: Science > Math Asked by: cwd-ga List Price: \$15.00 Posted: 28 May 2004 18:38 PDT Expires: 27 Jun 2004 18:38 PDT Question ID: 353406
 ```Fourier transforms are used for DEs and signal analysis, of course. But I've also heard they're used for reducing the number of iterations for a lengthy algorithm outside the realm of optics and signal processing. How does this work? Thanks.``` Request for Question Clarification by mathtalk-ga on 08 Jun 2004 12:31 PDT ```Hi, cwd-ga: I can think of several applications that might fit your description. Perhaps the most elementary one would be extended precision multiplication (finding the product of two extremely large integers, for example). When you ask "How does this work?", are you seeking a mathematical/theoretical explication (hey, this is the Science > Math department!), or more of a programming/practical "How do they do that?" account? regards, mathtalk-ga``` Clarification of Question by cwd-ga on 09 Jun 2004 20:50 PDT ```Either formulas or a computational example would suit me great. I'd enjoy/appreciate both theoretical and practical. Arfken has 30 pages on Fourier transforms, but the applications are on how this transform is a special case of that transform and the like. I'd heard that Fourier transforms can be used to bring an O(N) program down to O(lnN), or something of that order, and I'm curious how this is possible and what the calculations/code could look like. Sorry this is so vague. I'm an actuary and it occurred to me earlier today that I might find more what I'm looking for by Googling on "insurance" along with "Fourier transform". I've found a UConn publication that doesn't provide public access, so I might look into subscribing. Anyway, hope that helps.```
 ```Hi, cwd-ga: In 1965 Cooley and Tukey published an algorithm for computing the discrete Fourier transform on 2^n nodes using O(n log n) operations: [Fast Fourier Transform - Cooley and Tukey] http://en.wikipedia.org/wiki/Fast_Fourier_transform It was later found that the idea had been used by C.F. Gauss in 1805 and independently re-discovered by several others in the intervening 160 years. Nonetheless the paper by Cooley and Tukey gave rise to considerable interest in using these kinds of numerical computations to increase algorithmic efficiency. It should be mentioned that the fast Fourier transform, even when done with real valued input, will inherently involving complex arithmetic in intermediate steps. See the article above for descriptions of and links to other fast Fourier transform algorithms, applicable in varying degree to numbers of nodes other than a power of 2. A nice explanation with specifics was presented by my colleague rbnn-ga in thie earlier Answer: [Q: Fast Fourier Transform] http://answers.google.com/answers/threadview?id=92858 for how a computation for general values of n might be consumated using a power of 2 computation. In any case for the sake of illustration of the utility of fast Fourier transforms, I'll choose an application in which the number of nodes can be chosen with some flexibility, multiplying two extended precision integers. The key idea here is that in passing from the original space of array values to their discrete Fourier transform values, we have a relationship between the convolution of two arrays on one side and their pointwise multiplication on the other. Now a convolution of two arrays is just what happens when two polynomials are multiplied modulo X^N - 1, where N is the number of nodes. Polynomial arithmetic modulo X^N -1 causes the coefficients to "wrap around" after power X^(N-1) and "go back" to constants, etc. So let's assume we have two extended precision integers U and V, and let's regard these integers as polynomials: u(X) = SUM u_i X^i v(X) = SUM v_i X^i which are to be evaluated at X = radix (for whatever number base the integers are written in, e.g. binary if desired). Then the product of these polynomials will have coefficients which are a convolution of their coefficients: w(X) = SUM (SUM u_k * v_{i-k} FOR k=0,..,i) X^i or to write this in another way: w_i = SUM u_k * v_{i-k} FOR k=0,..,i Now if N is chosen to be sufficiently large, say N more than twice the maximum degree of u(X) or v(X), then there is no practical distinction between their product w(X) and their product modulo X^N - 1 (at least so far as the coefficients are concerned). So let N be a sufficiently large power of two, and let: F(u) = discrete Fourier transform of (u_i: i = 0,..,N-1) F(v) = discrete Fourier transform of (v_i: i = 0,..,N-1) in which the higher order coefficients are padded with zeros. Now F(u), F(v) are themselves arrays of length N, albeit with complex entries in the general case. And the relationship between convolution and pointwise multiplication gives this relation: F(w) = F(u)*F(v) where here we mean by "*" the entry-by-entry product of two arrays. Thus to recover the polynomial w(X) it's needed at this point to "undo" the discrete Fourier transform. But the inverse of this operation is closely related to itself: F(F(w)) = N w where here we mean on the right-hand side that each entry in array w has been scaled by a factor of N. [In fact the Cooley-Tukey fast Fourier transform involves n steps where N = 2^n, and on the way "back" if exact arithmetic were required, we could extract a factor of 2 at each such step, so that the factor N disappears along the way.] In other words undoing the discrete Fourier transform is really not more difficult than doing it in the first place. Having obtained the polynomial product w(X), the actual extended precision integer product W that obtains by evaluating X = radix, will involve "carries" at potentially all of the coefficients. However this additional computation, involving only additions and radix shifts, will for practical purposes be O(N) complexity. For further details of this see here: [FFT based multiplication of large numbers] http://numbers.computation.free.fr/Constants/Algorithms/fft.html Note that the "naive" (pencil and paper) multiplication of u(X) and v(X) requires something like the product of their degrees number of operations, for each coefficient ("digit") of one must be multiplied by each coefficient ("digit") in the other. With the FFT approach we have only N multiplications to perform (on the discrete Fourier transforms, which are complex values to be sure), and N is no more than twice the size of the maximum of these degrees (so within a constant of being linear in their degrees). This would be useless, though, without the FFT to save our effort. With that computation included the work becomes O(N log N). [Recall that the FFT had n steps where N = 2^n.] This is a big savings compared to the O(N^2) complexity of the naive algorithm. Many extended precision numerical packages will include an FFT implementation as one of several options. See here for description of one such package: [A Tour of NTL: NTL Implementation and Portability] http://www.wschnei.de/number-theory-software/ntl/tour-impl.html To summarize the fast Fourier transform allows us to multiply two N/2 "digit" integers, obtaining an exact N "digit" product, with O(N log N) operations instead of the O(N^2) operations needed by the pencil-and-paper approach we learn in grade school. Please let know by way of a Request for Clarification if additional information would be helpful. regards, mathtalk-ga``` Clarification of Answer by mathtalk-ga on 27 Jun 2004 23:15 PDT ```The above is an example of a "direct" method as opposed to an "iterative" one, so it's a bit easier to state the savings in complexity (since an exact answer is obtained in some finite number of operations). However I'd like to give some illustrations of FFT in iterative methods as well. First, however, let me work through a simple exercise to illustrate the ideas above. This will also give me an opportunity to clear up a couple of points of confusion in my presentation. I said that F(F(w)) = N w. Actually one needs to make a complex conjugation in how the second "Fourier transform" is done to get this result. To be clear, let's do a multiplication of two 2-digit numbers. For simplicity I've decided to square the number 14, which will also illustrate how "carries" come into the picture at the end. So we think of the (decimal number) 14 as being a polynomial, and we'll pad it out to four coefficients (because a 2-digit number times a 2-digit number will give at most a 4-digit product): 0X^3 + 0X^2 + 1X + 4 or simply, as the array of coefficients: u = [4, 1, 0, 0] putting the least significant term first in the array. The discrete Fourier transform is a linear transformation, which means it can be represented by a matrix multiplication (with possibly complex values). In the case of 4 nodes, the matrix is: [FFT - Some Simple Transforms] http://www.eptools.com/tn/T0001/PT01.HTM#Head687 +1 +1 +1 +1 +1 -i -1 +i +1 -1 +1 -1 +1 +i -1 -i If we multiply this matrix times its complex conjugate, we get 4I (verify!). Now multiply our vector u = [4, 1, 0, 0] by it: F(u) = [5, 4-i, 3, 4+i] Okay, so since we wanted to square 14, here we just multiply each entry of F(u) by itself: F(w) = [25, 15-8i, 9, 15+8i] Now multiply F(w) by the complex conjugate of our matrix above: 4w = [64, 25 + i(15-8i) - 9 - i(15+8i), 4, 25 - i(15-8i) - 9 + i(15+8i)] = [64, 32, 4, 0] w = [16, 8, 1, 0] Okay, so the polynomial w(X) is: 0X^3 + 1X^2 + 8X + 16 When we "evaluate" X = 10, then (proceeding from lowest order terms to highest), we remove any "excess" greater than or equal to 10 and promote it to the next higher coefficient. Thus we find: 14 squared = 1 6 8 1 -------------- 1 9 6 = 196 (base ten) That's how the carries work out, and in general it suffices to make a single sweep from lowest order to highest order (hence the estimate of O(N) work). Now this illustrates the discrete Fourier transform. What makes the method "fast". That would be a matter of cleverly arranging the arithmetic (in this particular case) into two steps (because 4 = 2^2), essentially decomposing the 4 node case into two 2 node cases. The cleverness involves some shuffling of entries, some adding and subtracting, and some multiplying by "roots of unity". The case at issue here would be computed something like this: 0 0 1 4 (original coefficients) 0 1 0 4 (shuffling based on bit reversals of indexes) 1 1 4 4 (sum and difference of adjacent pairs) -i 1 4 4 (multiplying by roots of unity) 4+i 3 4-i 5 (sum and difference of doubly spaced pairs) Next I'd like to discuss the appearance of the fast Fourier transform in the numerical solution of partial differential equations, an area in which iterative methods will find their place. regards, mathtalk-ga```