Algoritmos para calcular FFT em paralelo

Estou tentando paralelizar o cálculo de uma FFT em arquivos de sinal do tamanho de terabytes. No momento, essa FFT usando uma biblioteca de código aberto leva muitas horas, mesmo sendo executada no CUDA na GPU mais rápida que eu tenho. A estrutura que estou tentando adaptar a esse processo é o Hadoop. Em termos muito básicos, o Hadoop distribui um problema por qualquer número de nós do servidor da seguinte maneira:

• Você divide seu arquivo de entrada em pares (chave, valor).
• Esses pares são inseridos no algoritmo "Mapa", que transforma seus pares (chave, valor) em outros pares (chave, valor) com base no que você coloca dentro do Mapa.
• A estrutura coleta todas as saídas (chave, valor) dos Mapas e as classifica por chave, além de agregar valores com a mesma chave a um único par, para que você acabe com (chave, lista (valor1, valor2, ..)) pares
• Esses pares são inseridos no algoritmo "Reduzir", que, por sua vez, gera mais pares (chave, valor) como resultado final (gravado em um arquivo).

Existem muitos aplicativos para esse modelo em itens práticos, como o processamento de logs do servidor, mas estou tendo dificuldade em aplicar a estrutura para dividir uma FFT em tarefas de "mapa" e "redução", especialmente porque eu não estou realmente familiarizado com o DSP.

Não vou incomodá-lo com o mumbo jumbo de programação, pois este é um Q & A de DSP. No entanto, estou confuso sobre quais algoritmos existem para calcular FFTs em paralelo; As tarefas Mapear e Reduzir não podem (tecnicamente) se comunicar, portanto a FFT deve ser dividida em problemas independentes, dos quais os resultados podem de alguma forma ser recombinados no final.

Programei uma implementação simples do Cooley-Tukey Radix 2 DIT que funciona em pequenos exemplos, mas usá-lo para calcular recursivamente DFTs de índices ímpares / pares para um bilhão de bytes não funcionará. Passei algumas semanas lendo muitos artigos, incluindo um sobre o algoritmo MapReduce FFT (escrito por Tsz-Wo Sze como parte de seu artigo sobre multiplicação de SSA, não consigo vincular mais de 2 hiperlinks) e o “FFT em quatro etapas” ( aqui e aqui), que parecem semelhantes entre si e com o que estou tentando realizar. No entanto, sou irremediavelmente ruim em matemática, e aplicar qualquer um desses métodos manualmente a um conjunto simples de algo como {1,2, 3, 4, 5, 6, 7, 8} (com todos os componentes imaginários sendo 0) me resultados extremamente incorretos. Alguém pode me explicar um algoritmo FFT paralelo eficiente em inglês simples (um que vinculei ou outro) para que eu possa tentar programá-lo?

Edit: Jim Clay e qualquer outra pessoa que possa estar confusa com a minha explicação, estou tentando fazer uma única FFT do arquivo terabyte. Mas quero fazê-lo simultaneamente em vários servidores para acelerar o processo.

Eu acho que o seu principal problema não é como paralelizar o algoritmo (o que realmente pode ser feito), mas é a precisão numérica. FFTs de tamanho grande são numericamente bastante complicadas. Os coeficientes de FFT têm a forma e se N for muito grande, o cálculo do coeficiente fica barulhento. Digamos que você tenhaN=240e use aritmética de precisão dupla de 64 bits. Os primeiros 1000 coeficientes têm uma parte real que é exatamente a unidade (embora não deva ser assim), portanto você precisará de uma matemática de maior precisão, que é muito ineficiente e difícil de usar.$e^{-j\cdot 2\cdot \pi \cdot \frac{k }{N}}$$N = 2^{40}$

Você também acumulará muitos erros de arredondamento e truncamento, pois o grande número de operações que entram em um único número de saída também é muito grande. Devido à natureza "cada saída depende de cada entrada" da FFT, a propagação de erros é desenfreada.

Não conheço uma maneira fácil de solucionar isso. Seu pedido é incomum. A maioria dos aplicativos que fazem análise espectral de grandes conjuntos de dados faz uma análise em execução na qual você não tem esse problema. Talvez se você puder descrever sua aplicação e suas restrições, mas ainda mais, podemos apontar uma solução mais adequada.

Em vez de tentar reescrever a FFT, você pode tentar usar uma implementação existente da FFT (como a FFTW, por exemplo) e aplicá-la repetidamente ao longo do comprimento do seu sinal (não importa o tamanho), através da adição ou sobreposição de sobreposição salvar métodos. Isso é possível expressando a FFT como uma convolução .

Esses FFTs de menor comprimento não precisam se comunicar e todo o esquema corresponde às etapas de redução de mapa.

Em geral, o que você gostaria de fazer é dividir seu sinal X em segmentos menores que também podem se sobrepor (por exemplo, X [0:10], X [5:15], X [10:20] ... .). Execute a FFT nesses pequenos segmentos e recombine-os no final para produzir o final. Isso se encaixa muito bem com os operadores de redução de mapa.

Durante o "mapa", você pode gerar pares (chave, valor) com a "chave" sendo algum ID seqüencial de cada segmento (0,1,2,3,4,5, ....) e o "valor" sendo o INDEX (ou posição do arquivo) do primeiro valor de um segmento no arquivo do seu sinal. Portanto, por exemplo, se seu arquivo estiver cheio de INT32s, o índice do segundo segmento (acima) estará em 5 * sizeof (INT32). (Ou, se estiver em qualquer outro formato, você pode ter uma lib para isso)

Agora, cada trabalhador recebe uma (chave, valor) abre um arquivo, procura o ponto certo, lê M amostras dele (onde M está 10 acima), executa a FFT e a salva em um arquivo com algum nome, por exemplo " RES_ [INKEY] .dat "e retorna um par (chave, valor). Nesse caso, "chave" seria o ÍNDICE (o "valor" da tupla de entrada (chave, valor)) e "valor" seria o nome do arquivo que contém os resultados da FFT. (voltaremos a isso)

Em "reduzir", agora você pode implementar sobreposição-adição ou sobreposição-economia aceitando um (chave, valor) da etapa "mapa", abrindo esse arquivo, carregando os resultados da FFT, executando oa ou os e salvando-os em o índice correto no seu arquivo de saída. (Veja o pseudocódigo neste (ou neste ), a etapa "map" lida com o "yt = ..." em paralelo e a etapa "reduzir" lida com a parte "y (i, k) = ...".)

Pode ser necessário algum malabarismo de arquivo aqui para reduzir o tráfego na rede ou a carga de um servidor que pode conter seu arquivo de dados real.

Vamos supor que o tamanho dos seus dados seja $2^N$ . Almofada com zeros caso contrário. No seu caso, como você menciona tamanhos em "escala de terabytes", usaremos N = 40.

Como é um tamanho de FFT grande - mas absolutamente razoável para uma única máquina -, sugiro que você faça apenas uma única iteração Cooley-Tukey do radix N / 2 e, em seguida, deixe uma biblioteca FFT adequada (como FFTW ) faça o trabalho em cada máquina para o tamanho menor 2 N /$2^{N / 2}$$N / 2$$2 ^{N / 2}$ .

Para ser mais explícito, não há necessidade de usar MR ao longo de toda a recursão, isso será realmente ineficiente. Seu problema pode ser dividido em FFTs internas e externas de um milhão de megabytes, e essas FFTs de megabytes podem ser perfeitamente calculadas usando FFTW ou algo semelhante. O MR será apenas responsável por supervisionar a mistura e recombinação de dados, e não o cálculo real da FFT ...

Minha primeira idéia seria a seguinte, mas suspeito que isso possa ser feito em um único MR com representação de dados mais inteligente.

$s$$R = 2^{N / 2}$

Primeiro MR: FFT interna

Mapa: realize a dizimação no tempo, agrupe amostras em blocos para FFT interna

$(k, v)$$k$$0 .. 2^N - 1$$v$$s[k]$

$(k \% R, (k / R, v))$

Reduzir: calcular a FFT interna

$(k, vs)$$k$é o índice de bloco; e$vs$ é uma lista de $(i, v)$ pares

preencher um vetor $in$ de tamanho $R$ de tal modo que $in[i] = v$ para todos os valores na lista.

executar um tamanho $R$ FFT ativado $in$ para pegar um vetor $out$ de tamanho $R$

para $i$ dentro $0 .. R - 1$, emitir $(k, (i, out[i]))$

Segundo MR: FFT externa

Mapa: agrupe amostras para fft externo e aplique fatores de variação

entrada: $(k, (i, v))$ Onde $k$ é um índice de bloco, $(i, v)$ uma amostra da FFT interna para este bloco.

emitir $(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.

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).