Алгоритмы вычисления БПФ параллельно

Я пытаюсь распараллелить вычисление FFT на файлах сигналов терабайтового размера. Прямо сейчас такое БПФ, использующее библиотеку с открытым исходным кодом, занимает много часов, даже через CUDA на самом быстром из имеющихся у меня графических процессоров. Основой, которую я пытаюсь адаптировать к этому процессу, является Hadoop. В самых основных терминах Hadoop распределяет проблему по любому количеству узлов сервера следующим образом:

• Вы разбили входной файл на пары (ключ, значение).
• Эти пары подаются в алгоритм «Карта», который преобразует ваши пары (ключ, значение) в некоторые другие пары (ключ, значение) в зависимости от того, что вы поместили в карту.
• Затем платформа собирает все выходные данные (ключ, значение) из карт и сортирует их по ключу, а также объединяет значения с одним и тем же ключом в одну пару, поэтому в итоге вы получаете (ключ, список (значение1, значение2, ..)) пары
• Эти пары затем передаются в алгоритм «Reduce», который, в свою очередь, выводит больше пар (ключ, значение) в качестве вашего конечного результата (записанного в файл).

Существует множество приложений для этой модели в таких практических вещах, как обработка журналов сервера, но мне трудно применить среду для разбиения FFT на «карту» и «сокращения» задач, тем более что я не очень знаком с DSP.

Я не буду беспокоить вас программированием Mumbo Jumbo, так как это DSP Q & A. Однако я не совсем понимаю, какие существуют алгоритмы для параллельного вычисления БПФ; Задачи Map и Reduce не могут (технически) взаимодействовать друг с другом, поэтому БПФ должно быть разделено на независимые задачи, из которых результаты могут быть каким-то образом объединены в конце.

Я запрограммировал простую реализацию Cooley-Tukey Radix 2 DIT, которая работает на небольших примерах, но использовать ее для рекурсивного вычисления DFT с нечетным / четным индексом для миллиарда байтов не будет. Я потратил несколько недель на чтение многих работ, в том числе одного по алгоритму MapReduce FFT (написанного Tsz-Wo Sze как часть его статьи о умножении SSA, я не могу связать более двух гиперссылок) и «четырехступенчатого БПФ» ( здесь и здесь), которые кажутся похожими друг на друга и на то, что я пытаюсь достичь. Однако я безнадежно плохо разбираюсь в математике, и применение любого из этих методов вручную к простому набору чего-то вроде {1,2, 3, 4, 5, 6, 7, 8} (со всеми мнимыми компонентами равными 0) дает мне дико неверные результаты. Может ли кто-нибудь объяснить мне эффективный параллельный алгоритм FFT на простом английском языке (тот, который я связал, или любой другой), чтобы я мог попробовать его запрограммировать?

Редактировать: Джим Клэй и любой другой, кто может быть смущен моим объяснением, я пытаюсь сделать один FFT файла терабайта. Но я хочу сделать это одновременно на нескольких серверах, чтобы ускорить процесс.

Я думаю, что ваша главная проблема не в том, как распараллелить алгоритм (что на самом деле можно сделать), а в его числовой точности. БПФ такого размера в числовом отношении довольно сложны. Коэффициенты БПФ имеют вид и если N очень велико, вычисление коэффициента становится шумным. Допустим, у вас и вы используете 64-битную арифметику с двойной точностью. Первые 1000 коэффициентов имеют действительную часть, которая в точности равна единице (хотя так не должно быть), поэтому вам потребуется математика с более высокой точностью, которая очень неэффективна и громоздка в использовании. N=2$e^{-j\cdot 2\cdot \pi \cdot \frac{k }{N}}$$N = 2^{40}$

Вы также будете накапливать множество ошибок округления и усечения, поскольку само число операций, которые входят в одно выходное число, также очень велико. Из-за того, что БПФ «каждый выход зависит от каждого входа», распространение ошибок является безудержным.

Я не знаю простого способа обойти это. Ваш запрос необычный. Большинство приложений, которые выполняют спектральный анализ больших наборов данных, выполняют текущий анализ, где у вас нет такой проблемы. Возможно, если вы сможете описать свое приложение и его ограничения, но мы можем указать вам более подходящее решение.

Вместо того , чтобы пытаться переписать FFT вы можете попробовать использовать существующую реализацию FFT (например, в FFTW , например) и применить его повторно по всей длине вашего сигнала (независимо от того , насколько велико она есть) либо через перекрытие или перекрывающийся сохранить методы. Это возможно, выражая БПФ как свертку .

Эти БПФ более короткой длины не должны связываться друг с другом, и вся схема соответствует шагам сокращения карты.

В целом, то, что вы хотели бы сделать, это разделить сигнал X на более мелкие сегменты, которые также могут перекрываться (например, X [0:10], X [5:15], X [10:20] ... .). Выполните БПФ на этих маленьких сегментах и ​​объедините их в конце, чтобы получить последний. Это очень хорошо подходит для операторов уменьшения карты.

Во время «карты» вы можете генерировать пары (ключ, значение), где «ключ» - это некоторый последовательный идентификатор каждого сегмента (0,1,2,3,4,5, ....), а «значение» - это INDEX (или позиция файла) первого значения сегмента в файле вашего сигнала. Так, например, если ваш файл полон INT32, тогда индекс второго сегмента (выше) равен 5 * sizeof (INT32). (Или, если это в любом другом формате, у вас может быть библиотека для этого)

Теперь каждый работник получает (ключ, значение), открывает файл, ищет нужную точку, читает из него M выборок (где M на 10 выше), выполняет БПФ и сохраняет его в файл с некоторым именем, например " RES_ [INKEY] .dat "и возвращает пару (ключ, значение). В этом случае «ключ» будет INDEX («значение» входящего кортежа (ключ, значение)), а «значение» будет именем файла, содержащего результаты FFT. (мы вернемся к этому)

В рамках функции «Reduce» теперь вы можете реализовать либо перекрытие-добавление, либо сохранение-перекрытие, приняв (ключ, значение) из шага «карта», открыв этот файл, загрузив результаты FFT, выполнив oa или os, а затем сохранив их в правильный ИНДЕКС в вашем выходном файле. (См. Псевдокод в этом (или в этом ) шаге «map» обрабатывает «yt = ...» параллельно, а шаг «Reduce» обрабатывает часть «y (i, k) = ...».)

Здесь может потребоваться некоторая подтасовка файлов, чтобы уменьшить трафик в сети или нагрузку на сервер, который может содержать ваш фактический файл данных.

Предположим , что размер данных . В противном случае дополняется нулями. В вашем случае, поскольку вы упоминаете размеры в терабайтном масштабе, мы возьмем N = 40.$2^N$

Поскольку - это большой - но абсолютно разумный для одного компьютера - размер БПФ, я предлагаю вам сделать всего одну итерацию Кули-Тьюки с основанием , а затем дать подходящую библиотеку БПФ (как FFTW) выполняет работу на каждой машине для меньшего размера .$2^{N / 2}$$N / 2$$2 ^{N / 2}$

Чтобы быть более точным, нет необходимости использовать MR на протяжении всей рекурсии, это будет действительно неэффективно. Ваша проблема может быть разбита на внутренние и внешние БПФ размером в миллионы мегабайт, и эти БПФ в мегабайтах могут быть идеально рассчитаны с использованием БПФ или тому подобного. MR будет просто отвечать за контроль перетасовки и рекомбинации данных, а не за фактические вычисления FFT ...

Моей самой первой идеей было бы следующее, но я подозреваю, что это можно сделать в одном MR с более умным представлением данных.

Пусть будет вашим входным сигналом,$s$$R = 2^{N / 2}$

Первый МР: внутренний БПФ

Карта: выполнять прореживание во времени, группировать выборки по блокам для внутреннего БПФ

вход: где - индекс выборки в ; значение, принятое$(k, v)$$k$$0 .. 2^N - 1$$v$$s[k]$

emit: - где% представляет деление по модулю и / целое число.$(k \% R, (k / R, v))$

Уменьшить: вычислить внутреннее БПФ

вход: где - индекс блока; и это список пар$(k, vs)$$k$$vs$$(i, v)$

заполнить вектор размере таким образом, что для всех значений в списке.$in$$R$$in[i] = v$

выполнить размер FFT на , чтобы получить вектор размера$R$$in$$out$$R$

для в , выброс$i$$0 .. R - 1$$(k, (i, out[i]))$

Второе МР: внешнее БПФ

Карта: сгруппируйте образцы для внешнего БПФ и примените факторы твида

input: где - индекс блока, выборка внутреннего БПФ для этого блока.$(k, (i, v))$$k$$(i, v)$

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

Уменьшить: выполнить внешнее БПФ

вход: где - индекс блока; и это список пар$(k, vs)$$k$$vs$$(i, v)$

заполнить вектор размере таким образом, что для всех значений в списке.$in$$R$$in[i] = v$

выполнить размер FFT на , чтобы получить вектор размера$R$$in$$out$$R$

для в испускаем$i$$0 .. R - 1$$(i \times R + k, out[i]))$

Доказательство концепции Python-кода здесь.

Как видите, Mappers только перетасовывают порядок данных, поэтому при следующих допущениях:

• прореживание во времени (Mapper 1) может быть выполнено на предыдущем этапе (например, программой, которая преобразует данные в правильный формат ввода).
• Ваша среда MR поддерживает запись редукторов в ключ, отличный от их ключа ввода (в реализации Google редукторы могут выводить данные только на тот же ключ, который они получили, я думаю, это связано с использованием SSTable в качестве формата вывода).

Все это может быть сделано в одном MR, внутреннем БПФ в преобразователе, внешнем БПФ в редукторе. Доказательство концепции здесь .

Если ваш сигнал многомерен, то распараллеливание БПФ может быть выполнено довольно легко; сохранить одно измерение непрерывным в процессе MPI, выполнить БПФ и транспонировать (altoall) для работы в следующем измерении. FFTW делает это.

Если данные 1D, проблема гораздо сложнее. FFTW, например, не записал 1D FFT с использованием MPI. Если кто-то использует алгоритм прореживания по частоте радиуса-2, то первые несколько этапов могут быть выполнены как наивные ДПФ, позволяя использовать 2 или 4 узла без потери точности (это потому, что корни единства для первые этапы - либо -1, либо я, с которыми приятно работать).

Кстати, что вы планируете делать с данными после их преобразования? Это может быть что-то сделать, если вы знаете, что происходит с выводом (например, свертка, фильтр нижних частот и т. Д.).