Algorithms for computing FFT in parallel

Question

I am trying to parallelize the computation of an FFT on terabyte-sized signal files. Right now such an FFT using an open-source library takes many hours, even running through CUDA on the fastest GPU I have. The framework I am trying to adapt to this process is Hadoop. In very basic terms, Hadoop distributes a problem over any number of server nodes in the following manner:

• You split your input file into (key, value) pairs.
• These pairs are fed into a “Map” algorithm, which transforms your (key, value) pairs into some other (key, value) pairs based on what you put inside the Map.
• The framework then collects all the (key, value) outputs from the Maps and sorts them by key, as well as aggregating values with the same key to a single pair, so you end up with (key, list(value1, value2, ..)) pairs
• These pairs are then fed into a “Reduce” algorithm, which in turn outputs more (key, value) pairs as your final result (written to a file).

There are many applications for this model in practical stuff like processing server logs, but I am having a hard time applying the framework to chopping up an FFT into “map” and “reduce” tasks, especially since I am not really familiar with DSP.

I will not bother you with the programming mumbo jumbo, as this is a DSP Q&A. I am, however, confused on what algorithms exist for computing FFTs in parallel; Map and Reduce tasks can’t (technically) talk to each other, so the FFT must be split into independent problems from which the results can somehow be recombined at the end.

I have programmed a simple implementation of Cooley-Tukey Radix 2 DIT that works on small examples, but using it for recursively calculating odd/even index DFTs for a billion bytes will not work. I have spent a few weeks on reading many papers, including one on a MapReduce FFT algorithm (written by Tsz-Wo Sze as part of his paper on SSA multiplication, I cannot link more than 2 hyperlinks) and the “four-step FFT” (here and here), which seem similar to each other and to what I am trying to accomplish. However, I am hopelessly bad at mathematics, and applying any of those methods by hand to a simple set of something like {1,2, 3, 4, 5, 6, 7, 8} (with all imaginary components being 0) gives me wildly incorrect results. Can anyone explain an efficient parallel FFT algorithm to me in plain English (one that I linked or any other) so that I may try and program it?

Edit: Jim Clay and anyone else who may be confused by my explanation, I am trying to do a single FFT of the terabyte file. But I want to able to do it concurrently on multiple servers in order to speed up the process.

Answer 1

I think your main problem isn't how to parallel the algorithm (which actually can be done) but it's numeric precision. FFTs of a size that large are numerically quite tricky. The FFT coefficients are of the form $e^{-j\cdot 2\cdot \pi \cdot \frac{k }{N}}$ and if N is very large the coefficient calculation gets noisy. Lets say you have $N = 2^{40}$ and you use 64 bit double precision arithmetic. The first 1000 coefficients have a real part that's exactly unity (although it shouldn't be that way) so you will need higher precision math, which is very inefficient and cumbersome to use.

You will also amass lots of rounding and truncation errors as the sheer number of operations that go into a single output number is very large too. Due to the "every output depends on every input" nature of the FFT, error propagation is rampant.

I'm not aware of an easy way to work around that. Your request is an unusual one. Most applications that do spectral analysis of large data sets do a running analysis where you don't have that problem. Perhaps if you can describe your application and it's constraints a but more, we can point you to a more suitable solution.

Answer 2

Instead of trying to re-write the FFT you could try using an existing FFT implementation (such as the FFTW for example) and apply it repetitively along your signal's length (no matter how big it is) through either the overlap-add or overlap-save methods. This is possible by expressing the FFT as a convolution.

These shorter length FFTs do not need to communicate with each other and the whole scheme matches the map-reduce steps.

In general, what you would aim to do is to have your signal X split into smaller segments that also could be overlapping (for example X[0:10],X[5:15],X[10:20]....). Perform the FFT on these small segments and recombine them in the end to produce the final one. This fits very well with the map-reduce operators.

During "map" you can generate (key,value) pairs with the "key" being some sequential ID of each segment (0,1,2,3,4,5,....) and the "value" being the INDEX (or file position) of the first value of a segment in your signal's file. So, for instance, if your file is full of INT32s then the index of the second segment (above) is at 5*sizeof(INT32). (Or if it's in any other format you might have a lib for it)

Now, each worker receives a (key,value) opens a file, seeks to the right point, reads M samples from it (where M is 10 above), performs the FFT and saves it to a file with some name, for example "RES_[INKEY].dat" and returns a (key,value) pair. In this case, "key" would be the INDEX (the "value" of the incoming (key,value) tuple) and "value" would be the name of the file that contains the FFT results. (we'll get back to this)

Within "reduce" you can now implement either overlap-add or overlap-save by accepting a (key,value) from the "map" step, openning that file, loading the FFT results, doing either o-a or o-s and then saving them to the right INDEX in your output file. (See pseudocode in this (or this), the "map" step handles the "yt=..." in parallel and the "reduce" step handles the "y(i,k)=..." part.)

Some file juggling might be needed here to reduce the traffic on the network or the load of a server that might contain your actual data file.

Answer 3

Let us assume that your data size is $2^N$. Pad with zeros otherwise. In your case, since you mention "Terabyte-scale" sizes, we'll take N = 40.

Since $2^{N / 2}$ is a big -- but absolutely reasonable for a single machine -- FFT size, I suggest you to do just one single Cooley-Tukey iteration of radix $N / 2$, and then let a proper FFT library (like FFTW) do the job on each machine for the smaller size $2 ^{N / 2}$.

To be more explicit, there's no need to use MR along the whole recursion, this will be quite inefficient indeed. Your problem can be broken down into a million of megabyte size inner and outer FFTs, and those megabyte FFTs can perfectly be computed using FFTW or the like. MR will just be responsible for supervising the data shuffling and recombination, not the actual FFT computation...

My very first idea would be the following, but I suspect this can be done in one single MR with smarter data representation.

Let $s$ be your input signal, $R = 2^{N / 2}$

First MR: inner FFT

Map: perform decimation in time, group samples in blocks for inner FFT

input: $(k, v)$ where $k$ is the sample index in $0 .. 2^N - 1$ ; $v$ the value taken by $s[k]$

emit: $(k \% R, (k / R, v))$ - where % represents modulo and / integer division.

Reduce: compute inner FFT

input: $(k, vs)$ where $k$ is the block index ; and $vs$ is a list of $(i, v)$ pairs

populate a vector $in$ of size $R$ such that $in[i] = v$ for all values in the list.

perform a size $R$ FFT on $in$ to get a vector $out$ of size $R$

for $i$ in $0 .. R - 1$, emit $(k, (i, out[i]))$

Second MR: outer FFT

Map: group samples for outer fft and apply twiddle factors

input: $(k, (i, v))$ where $k$ is a block index, $(i, v)$ a sample of the inner FFT for this block.

emit $(i, (k, v \times \exp \frac{-2 \pi j i k}{2 ^ N}))$

Reduce: perform outer FFT

input: $(k, vs)$ where $k$ is the block index ; and $vs$ is a list of $(i, v)$ pairs

populate a vector $in$ of size $R$ such that $in[i] = v$ for all values in the list.

perform a size $R$ FFT on $in$ to get a vector $out$ of size $R$

for $i$ in $0 .. R - 1$, emit $(i \times R + k, out[i]))$

Proof of concept python code here.

As you can see, the Mappers are only shuffling the order of data, so under the following assumptions:

• decimation in time (Mapper 1) can be done in a previous step (for example by the program that converts the data into the right input format).
• your MR framework supports Reducers writing to a key different from their input key (In Google's implementation reducers can only output data to the same key as they received it, I think it's due to SSTable being used as an output format).

All this can be done in one single MR, the inner FFT in the mapper, the outer FFT in the reducer. Proof of concept here.

Answer 4

If your signal is multi-dimensional, then parallelizing the FFT can be accomplished fairly easily; keep one dimension contiguous in an MPI process, perform the FFT, and transpose (altoall) to work on the next dimension. FFTW does this.

If the data is 1D, the problem is much more difficult. FFTW, for example, did not write a 1D FFT using MPI. If one uses a radix-2 decimation-in-frequency algorithm, then the first few stages can be performed as a naive DFT, allowing one to use 2 or 4 nodes without any loss of precision (this is because the roots of unity for the first stages are either -1 or i, which are nice with which to work).

Incidentally, what are you planning on doing with the data once you've transformed it? It might be do something if one knows what happens to the output (i.e. a convolution, low-pass filter, etc).