Algorithmes de calcul de la FFT en parallèle

J'essaie de paralléliser le calcul d'une FFT sur des fichiers de signaux de taille téraoctet. À l'heure actuelle, une telle FFT utilisant une bibliothèque open-source prend de nombreuses heures, même en passant par CUDA sur le GPU le plus rapide que j'ai. Le cadre que j'essaie d'adapter à ce processus est Hadoop. En termes très basiques, Hadoop distribue un problème sur n'importe quel nombre de nœuds de serveur de la manière suivante:

• Vous divisez votre fichier d'entrée en paires (clé, valeur).
• Ces paires sont introduites dans un algorithme de «carte», qui transforme vos paires (clé, valeur) en d'autres paires (clé, valeur) en fonction de ce que vous mettez à l'intérieur de la carte.
• Le cadre collecte ensuite toutes les sorties (clé, valeur) des cartes et les trie par clé, ainsi que l'agrégation des valeurs avec la même clé en une seule paire, de sorte que vous vous retrouvez avec (clé, liste (valeur1, valeur2, ..)) paires
• Ces paires sont ensuite introduites dans un algorithme «Réduire», qui à son tour génère plus de paires (clé, valeur) comme résultat final (écrit dans un fichier).

Il existe de nombreuses applications pour ce modèle dans des trucs pratiques comme le traitement des journaux de serveur, mais j'ai du mal à appliquer le cadre pour découper une FFT en «cartographier» et «réduire» les tâches, d'autant plus que je ne suis pas vraiment familier avec DSP.

Je ne vais pas vous déranger avec le jumbo de programmation mumbo, car il s'agit d'un Q&A DSP. Je suis cependant confus sur les algorithmes qui existent pour calculer les FFT en parallèle; Les tâches de cartographie et de réduction ne peuvent pas (techniquement) communiquer entre elles, la FFT doit donc être divisée en problèmes indépendants à partir desquels les résultats peuvent en quelque sorte être recombinés à la fin.

J'ai programmé une implémentation simple de Cooley-Tukey Radix 2 DIT qui fonctionne sur de petits exemples, mais en l'utilisant pour calculer récursivement les DFT d'indice impair / pair pour un milliard d'octets ne fonctionnera pas. J'ai passé quelques semaines à lire de nombreux articles, dont un sur un algorithme MapReduce FFT (écrit par Tsz-Wo Sze dans le cadre de son article sur la multiplication SSA, je ne peux pas lier plus de 2 hyperliens) et la «FFT en quatre étapes» ( ici et ici), qui semblent similaires les uns aux autres et à ce que j'essaie d'accomplir. Cependant, je suis désespérément mauvais en mathématiques, et appliquer l'une de ces méthodes à la main à un simple ensemble de quelque chose comme {1,2, 3, 4, 5, 6, 7, 8} (avec tous les composants imaginaires étant 0) donne moi des résultats extrêmement incorrects. Quelqu'un peut-il m'expliquer un algorithme FFT parallèle efficace en anglais simple (celui que j'ai lié ou tout autre) afin que je puisse essayer de le programmer?

Edit: Jim Clay et toute autre personne qui pourrait être confuse par mon explication, j'essaie de faire une seule FFT du fichier téraoctet. Mais je veux pouvoir le faire simultanément sur plusieurs serveurs afin d'accélérer le processus.

Je pense que votre problème principal n'est pas de savoir comment paralléliser l'algorithme (ce qui peut en fait être fait) mais sa précision numérique. Les FFT d'une taille aussi grande sont numériquement assez délicates. Les coefficients FFT sont de la forme et si N est très grand le calcul du coefficient devient bruyant. Disons que vous avez et que vous utilisez une arithmétique double précision 64 bits. Les 1000 premiers coefficients ont une partie réelle qui est exactement l'unité (bien que cela ne devrait pas être le cas), vous aurez donc besoin de mathématiques de plus grande précision, ce qui est très inefficace et lourd à utiliser. N=2$e^{-j\cdot 2\cdot \pi \cdot \frac{k }{N}}$$N = 2^{40}$

Vous accumulerez également de nombreuses erreurs d'arrondi et de troncature, car le nombre d'opérations qui entrent dans un seul numéro de sortie est également très important. En raison de la nature «chaque sortie dépend de chaque entrée» de la FFT, la propagation des erreurs est endémique.

Je ne suis pas au courant d'un moyen facile de contourner cela. Votre demande est inhabituelle. La plupart des applications qui effectuent une analyse spectrale de grands ensembles de données effectuent une analyse en cours là où vous n'avez pas ce problème. Peut-être que si vous pouvez décrire votre application et ses contraintes, mais nous pouvons vous indiquer une solution plus appropriée.

Au lieu d'essayer de réécrire la FFT vous pouvez essayer d' utiliser une implémentation FFT existante (comme le FFTW par exemple) et l' appliquer répétitivement le long de la longueur de votre signal (peu importe comment elle est grande) soit par le chevauchement ajouter ou overlap- enregistrer les méthodes. Ceci est possible en exprimant la FFT sous forme de convolution .

Ces FFT de longueur plus courte n'ont pas besoin de communiquer entre elles et l'ensemble du schéma correspond aux étapes de réduction de carte.

En général, ce que vous viseriez à faire, c'est de diviser votre signal X en segments plus petits qui pourraient également se chevaucher (par exemple X [0:10], X [5:15], X [10:20] ... .). Effectuer la FFT sur ces petits segments et les recombiner à la fin pour produire le dernier. Cela correspond très bien aux opérateurs de réduction de carte.

Pendant la "carte", vous pouvez générer des paires (clé, valeur) avec la "clé" étant un ID séquentiel de chaque segment (0,1,2,3,4,5, ....) et la "valeur" étant le INDEX (ou position du fichier) de la première valeur d'un segment dans le fichier de votre signal. Ainsi, par exemple, si votre fichier est plein d'INT32, alors l'index du deuxième segment (ci-dessus) est à 5 * sizeof (INT32). (Ou si c'est dans un autre format, vous pourriez avoir une bibliothèque pour ça)

Maintenant, chaque travailleur reçoit une (clé, valeur) ouvre un fichier, cherche au bon point, en lit M échantillons (où M est 10 ci-dessus), effectue la FFT et l'enregistre dans un fichier avec un nom, par exemple " RES_ [INKEY] .dat "et renvoie une paire (clé, valeur). Dans ce cas, "clé" serait l'INDEX (la "valeur" du tuple entrant (clé, valeur)) et "valeur" serait le nom du fichier qui contient les résultats de la FFT. (nous y reviendrons)

Dans "réduire", vous pouvez maintenant implémenter soit chevauchement-ajout, soit chevauchement-sauvegarde en acceptant une (clé, valeur) à partir de l'étape "carte", en ouvrant ce fichier, en chargeant les résultats FFT, en faisant oa ou os puis en les enregistrant le bon INDEX dans votre fichier de sortie. (Voir le pseudocode dans ceci (ou cela ), l'étape "map" gère le "yt = ..." en parallèle et l'étape "réduire" gère la partie "y (i, k) = ...".)

Certains jonglages de fichiers peuvent être nécessaires ici pour réduire le trafic sur le réseau ou la charge d'un serveur qui pourrait contenir votre fichier de données réel.

Supposons que votre taille de données est . Pad avec des zéros sinon. Dans votre cas, puisque vous mentionnez des tailles "à l'échelle du téraoctet", nous prendrons N = 40.$2^N$

Étant donné que est une grande taille FFT, mais tout à fait raisonnable pour une seule machine, je vous suggère de faire une seule itération Cooley-Tukey de Radix , puis de laisser une bibliothèque FFT appropriée (comme FFTW) faire le travail sur chaque machine pour la plus petite taille .$2^{N / 2}$$N / 2$$2 ^{N / 2}$

Pour être plus explicite, il n'est pas nécessaire d'utiliser MR tout au long de la récursivité, ce sera en effet assez inefficace. Votre problème peut être décomposé en millions de FFT internes et externes de taille mégaoctet, et ces FFT mégaoctets peuvent parfaitement être calculés à l'aide de FFTW ou similaire. MR sera juste responsable de superviser le brassage et la recombinaison des données, pas le calcul FFT réel ...

Ma toute première idée serait la suivante, mais je pense que cela peut être fait dans un seul MR avec une représentation des données plus intelligente.

Soit votre signal d'entrée,$s$$R = 2^{N / 2}$

Premier MR: FFT intérieur

Carte: effectuer la décimation dans le temps, grouper les échantillons en blocs pour la FFT intérieure

entrée: où est l'indice de l'échantillon dans ; la valeur prise par$(k, v)$$k$$0 .. 2^N - 1$$v$$s[k]$

émettre: - où% représente la division modulo et / entier.$(k \% R, (k / R, v))$

Réduire: calculer la FFT intérieure

entrée: où est l'indice de bloc; et est une liste de paires$(k, vs)$$k$$vs$$(i, v)$

peupler un vecteur la taille de telle sorte que pour toutes les valeurs dans la liste.$in$$R$$in[i] = v$

effectuer une taille FFT sur pour obtenir un vecteur la taille$R$$in$$out$$R$

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

Deuxième MR: FFT externe

Carte: regrouper les échantillons pour la zone extérieure et appliquer les facteurs de torsion

entrée: où est un index de bloc, un échantillon de la FFT interne pour ce bloc.$(k, (i, v))$$k$$(i, v)$

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

Réduire: effectuer une FFT externe

entrée: où est l'indice de bloc; et est une liste de paires$(k, vs)$$k$$vs$$(i, v)$

peupler un vecteur la taille de telle sorte que pour toutes les valeurs dans la liste.$in$$R$$in[i] = v$

effectuer une taille FFT sur pour obtenir un vecteur la taille$R$$in$$out$$R$

pour in , emit0 . . R - 1 ( i × R + k , o u t [ i ] ) $i$$0 .. R - 1$$(i \times R + k, out[i]))$

Preuve de concept du code python ici.

Comme vous pouvez le voir, les mappeurs ne font que mélanger l'ordre des données, donc sous les hypothèses suivantes:

• la décimation dans le temps (Mapper 1) peut être effectuée dans une étape précédente (par exemple par le programme qui convertit les données dans le bon format d'entrée).
• votre infrastructure MR prend en charge les réducteurs écrivant sur une clé différente de leur clé d'entrée (dans l'implémentation de Google, les réducteurs ne peuvent sortir les données que sur la même clé qu'ils l'ont reçue, je pense que cela est dû au fait que SSTable est utilisé comme format de sortie).

Tout cela peut être fait dans un seul MR, la FFT intérieure dans le mappeur, la FFT extérieure dans le réducteur. Preuve de concept ici .

Si votre signal est multidimensionnel, la parallélisation de la FFT peut être effectuée assez facilement; conserver une dimension contiguë dans un processus MPI, effectuer la FFT et transposer (altoall) pour travailler sur la dimension suivante. La FFTW le fait.

Si les données sont 1D, le problème est beaucoup plus difficile. FFTW, par exemple, n'a pas écrit de FFT 1D à l'aide de MPI. Si l'on utilise un algorithme de décimation en fréquence radix-2, alors les premières étapes peuvent être effectuées comme un DFT naïf, permettant d'utiliser 2 ou 4 nœuds sans aucune perte de précision (c'est parce que les racines de l'unité pour le les premières étapes sont -1 ou i, ce qui est agréable avec lequel travailler).

Par ailleurs, que comptez-vous faire avec les données une fois que vous les avez transformées? Cela pourrait être de faire quelque chose si l'on sait ce qui arrive à la sortie (c'est-à-dire une convolution, un filtre passe-bas, etc.).