Algoritmos para calcular FFT en paralelo

Estoy tratando de paralelizar el cálculo de una FFT en archivos de señal de tamaño terabyte. En este momento, una FFT que utiliza una biblioteca de código abierto lleva muchas horas, incluso a través de CUDA en la GPU más rápida que tengo. El marco que estoy tratando de adaptar a este proceso es Hadoop. En términos muy básicos, Hadoop distribuye un problema en cualquier número de nodos del servidor de la siguiente manera:

• Divide su archivo de entrada en pares (clave, valor).
• Estos pares se introducen en un algoritmo de "Mapa", que transforma sus pares (clave, valor) en otros pares (clave, valor) en función de lo que coloque dentro del Mapa.
• Luego, el marco recopila todos los resultados (clave, valor) de los mapas y los ordena por clave, además de agregar valores con la misma clave a un solo par, por lo que termina con (clave, lista (valor1, valor2, ..)) pares
• Estos pares se introducen en un algoritmo "Reducir", que a su vez genera más pares (clave, valor) como resultado final (escrito en un archivo).

Hay muchas aplicaciones para este modelo en cuestiones prácticas como el procesamiento de registros del servidor, pero me resulta difícil aplicar el marco para dividir un FFT en tareas de "mapa" y "reducir", especialmente porque no estoy realmente familiarizado con DSP.

No te molestaré con la programación mumbo jumbo, ya que este es un DSP Q&A. Sin embargo, estoy confundido sobre qué algoritmos existen para calcular FFT en paralelo; Las tareas de mapeo y reducción no pueden (técnicamente) comunicarse entre sí, por lo que la FFT debe dividirse en problemas independientes a partir de los cuales los resultados se pueden recombinar de alguna manera al final.

He programado una implementación simple de Cooley-Tukey Radix 2 DIT que funciona en pequeños ejemplos, pero usarlo para calcular recursivamente los DFT de índice impar / par por mil millones de bytes no funcionará. He pasado algunas semanas leyendo muchos documentos, incluido uno en un algoritmo MapReduce FFT (escrito por Tsz-Wo Sze como parte de su artículo sobre la multiplicación SSA, no puedo vincular más de 2 hipervínculos) y el "FFT de cuatro pasos" ( aquí y aquí), que parecen similares entre sí y con lo que estoy tratando de lograr. Sin embargo, soy irremediablemente malo en matemáticas, y aplicar cualquiera de esos métodos a mano a un conjunto simple de algo como {1,2, 3, 4, 5, 6, 7, 8} (con todos los componentes imaginarios siendo 0) da Me resultados muy incorrectos. ¿Alguien puede explicarme un algoritmo FFT paralelo eficiente en inglés simple (uno que vinculé u otro) para que pueda intentar programarlo?

Editar: Jim Clay y cualquier otra persona que pueda estar confundida por mi explicación, estoy tratando de hacer una sola FFT del archivo de terabytes. Pero quiero poder hacerlo simultáneamente en varios servidores para acelerar el proceso.

Creo que su principal problema no es cómo paralelizar el algoritmo (que en realidad se puede hacer) sino la precisión numérica. Los FFT de un tamaño tan grande son numéricamente bastante complicados. Los coeficientes FFT tienen la forma y si N es muy grande, el cálculo del coeficiente se vuelve ruidoso. Digamos que tiene y utiliza aritmética de doble precisión de 64 bits. Los primeros 1000 coeficientes tienen una parte real que es exactamente la unidad (aunque no debería ser así), por lo que necesitará una matemática de mayor precisión, que es muy ineficiente y engorrosa de usar.$e^{-j\cdot 2\cdot \pi \cdot \frac{k }{N}}$$N = 2^{40}$

También acumulará muchos errores de redondeo y truncamiento, ya que la gran cantidad de operaciones que van a un solo número de salida también es muy grande. Debido a la naturaleza de "cada salida depende de cada entrada" de la FFT, la propagación de errores es rampante.

No conozco una manera fácil de evitarlo. Su solicitud es inusual. La mayoría de las aplicaciones que realizan análisis espectrales de grandes conjuntos de datos hacen un análisis en ejecución donde no tiene ese problema. Quizás si puede describir su aplicación y sus limitaciones, pero más, podemos indicarle una solución más adecuada.

En lugar de tratar de re-escribir la FFT se podría tratar de usar una implementación FFT existente (como la FFTW por ejemplo) y aplicarlo de forma repetitiva a lo largo de la longitud de su señal (no importa lo grande que es) a través de ya sea el solapesuma o se traslapan guardar métodos Esto es posible expresando la FFT como una convolución .

Estas FFT de menor longitud no necesitan comunicarse entre sí y todo el esquema coincide con los pasos de reducción de mapas.

En general, lo que querría hacer es dividir su señal X en segmentos más pequeños que también podrían superponerse (por ejemplo, X [0:10], X [5:15], X [10:20] ... .). Realice la FFT en estos pequeños segmentos y recombínelos al final para producir el último. Esto encaja muy bien con los operadores de reducción de mapas.

Durante el "mapa" puede generar pares (clave, valor) con la "clave" como una identificación secuencial de cada segmento (0,1,2,3,4,5, ....) y el "valor" es el ÍNDICE (o posición del archivo) del primer valor de un segmento en el archivo de su señal. Entonces, por ejemplo, si su archivo está lleno de INT32s, entonces el índice del segundo segmento (arriba) está en 5 * sizeof (INT32). (O si está en cualquier otro formato, puede que tenga una lib para ello)

Ahora, cada trabajador recibe una (clave, valor), abre un archivo, busca el punto correcto, lee M muestras de él (donde M es 10 arriba), realiza el FFT y lo guarda en un archivo con algún nombre, por ejemplo " RES_ [INKEY] .dat "y devuelve un par (clave, valor). En este caso, "clave" sería el ÍNDICE (el "valor" de la tupla entrante (clave, valor)) y "valor" sería el nombre del archivo que contiene los resultados FFT. (volveremos a esto)

Dentro de "reducir" ahora puede implementar overlap-add o overlap-save aceptando una (clave, valor) del paso "map", abriendo ese archivo, cargando los resultados de FFT, haciendo ya sea o os y luego guardándolos en el ÍNDICE correcto en su archivo de salida. (Ver pseudocódigo en este (o este ), el paso "mapa" maneja el "yt = ..." en paralelo y el paso "reducir" maneja la parte "y (i, k) = ...".)

Es posible que se necesiten algunos malabarismos de archivos aquí para reducir el tráfico en la red o la carga de un servidor que pueda contener su archivo de datos real.

Supongamos que el tamaño de los datos es . Almohadilla con ceros de lo contrario. En su caso, dado que menciona los tamaños de "escala de terabyte", tomaremos N = 40.$2^N$

Como es un tamaño FFT grande, pero absolutamente razonable para una sola máquina, le sugiero que haga una sola iteración Cooley-Tukey de radix , y luego deje una biblioteca FFT adecuada (como FFTW) haga el trabajo en cada máquina para el tamaño más pequeño .$2^{N / 2}$$N / 2$$2 ^{N / 2}$

Para ser más explícito, no hay necesidad de usar MR a lo largo de toda la recursión, esto será bastante ineficiente. Su problema puede dividirse en un millón de FFT de tamaño interno y externo de megabytes, y esos FFT de megabyte se pueden calcular perfectamente utilizando FFTW o similar. MR solo será responsable de supervisar la combinación de datos y la recombinación, no el cálculo real de FFT ...

Mi primera idea sería la siguiente, pero sospecho que esto se puede hacer en un solo MR con una representación de datos más inteligente.

Sea su señal de entrada,$s$$R = 2^{N / 2}$

Primer MR: FFT interno

Mapa: realizar diezmado en el tiempo, agrupar muestras en bloques para FFT interna

entrada: donde es el índice de muestra en ; el valor tomado por$(k, v)$$k$$0 .. 2^N - 1$$v$$s[k]$

emitir: - donde% representa la división de módulo / entero.$(k \% R, (k / R, v))$

Reducir: calcular FFT interno

entrada: donde es el índice de bloque; y es una lista de pares$(k, vs)$$k$$vs$$(i, v)$

poblar un vector de tamaño de tal manera que para todos los valores en la lista.$in$$R$$in[i] = v$

realice un tamaño FFT para obtener un vector del tamaño$R$$in$$out$$R$

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

Segundo MR: FFT externo

Mapa: agrupe las muestras para fft externo y aplique factores de giro

entrada: donde es un índice de bloque, una muestra de la FFT interna para este bloque.$(k, (i, v))$$k$$(i, v)$

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

Reducir: realizar FFT externa

entrada: donde es el índice de bloque; y es una lista de pares$(k, vs)$$k$$vs$$(i, v)$

poblar un vector de tamaño de tal manera que para todos los valores en la lista.$in$$R$$in[i] = v$

realice un tamaño FFT para obtener un vector del tamaño$R$$in$$out$$R$

para en , emitir0 . . R - 1 ( i × R + k , o u t [ i ] ) $i$$0 .. R - 1$$(i \times R + k, out[i]))$

Prueba de concepto de código python aquí.

Como puede ver, los mapeadores solo están barajando el orden de los datos, por lo que bajo los siguientes supuestos:

• la decimación en el tiempo (Mapper 1) se puede realizar en un paso anterior (por ejemplo, mediante el programa que convierte los datos al formato de entrada correcto).
• su marco MR admite reductores que escriben en una clave diferente de su clave de entrada (en la implementación de Google, los reductores solo pueden enviar datos a la misma clave que la recibieron, creo que se debe a que SSTable se utiliza como formato de salida).

Todo esto se puede hacer en una sola MR, la FFT interna en el mapeador, la FFT externa en el reductor. Prueba de concepto aquí .

Si su señal es multidimensional, entonces la paralelización de la FFT se puede lograr con bastante facilidad; mantenga una dimensión contigua en un proceso MPI, realice la FFT y transponga (todo) para trabajar en la siguiente dimensión. FFTW hace esto.

Si los datos son 1D, el problema es mucho más difícil. FFTW, por ejemplo, no escribió un 1D FFT usando MPI. Si se usa un algoritmo de decimación en frecuencia radix-2, entonces las primeras etapas se pueden realizar como un DFT ingenuo, lo que permite usar 2 o 4 nodos sin pérdida de precisión (esto se debe a las raíces de la unidad para el las primeras etapas son -1 o i, que son buenas para trabajar).

Por cierto, ¿qué planeas hacer con los datos una vez que los hayas transformado? Podría hacer algo si uno sabe lo que le sucede a la salida (es decir, una convolución, filtro de paso bajo, etc.).