Menu

Параллельное вычисление кумулятивных значений на примере поразрядного сложения

Мы коснёмся параллельного подхода в некоторых задачах, которые обычно решаются однопроходным алгоритмом, последовательно перебирая все элементы массива длиной $$ N $$. Сюда может входить вычисление кумулятивной суммы, бэктест в трейдинге и многое другое. В некоторых случаях удаётся достичь времени работы $$ O(\log N) $$, что является существенным ускорением на больших массивах, делая рассматриваемый подход применимым особенно на GPU, где можно достичь тысяч и более параллельных потоков. Ниже мы рассмотрим классический пример, в котором реализуем побитовое сложение двух длинных чисел. ## Рассматриваемая задача Нам даны два массива бит длиной $$ N $$ некоторых целых неотрицательных чисел от $$ 0 $$ до $$ 2^N - 1 $$. Необходимо построить параллельный алгоритм, который вычисляет биты суммы данных чисел. Наивно можно запустить побитовое сложение, начиная от младшего разряда, складывая соответствующие биты и не забывая про перенос "1 в уме", влияющий на значение следующего разряда. Именно из-за этого переноса нельзя сложить все биты параллельно и независимо, сразу предъявив финальный результат, как это было бы, например, для побитовых операций AND, OR, XOR и тд, имея время выполнения $$ O(1) $$ при достаточно большом количестве параллельных потоков. Тем не менее существует алгоритм, учитывающий перенос и работающий за $$ O(\log N) $$, основанный на бинарном префиксном дереве. О нём и пойдёт речь ниже. Его идея может быть обобщена до существенного класса задач. ## Алгоритм Для простоты в качестве $$ N $$ будем рассматривать степени двойки: $$ N = 2^n $$ Пусть первое число $$ a $$, а второе $$ b $$. Запишем массивы их бит от старшего разряда к младшему, начиная от 0. Пусть для числа $$ a $$ биты представляют из себя $$ a_0, a_1, ... a_{N-1} \in \{0,1\} $$. Аналогично $$ b_0, b_1, ... b_{N-1} \in \{0,1\} $$ для числа $$ b $$. Код для вычисления разрядов суммы `s` и переноса `c` (может быть распаралелен): ```python s_bits = [(ai ^ bi) for (ai, bi) in zip(a_bits, b_bits)] c_bits = [(ai & bi) for (ai, bi) in zip(a_bits, b_bits)] ``` Будет ли конкретный i-й разряд производить перенос, можно задать числом: $$ c_i \in \{0,1,2\} $$, где $$ c_i = 0 $$ означает, что переноса точно нет, $$ c_i = 1 $$ - перенос точно есть, а $$ c_i = 2 $$ - перенос происходит в том случае, если на предыдущем разряде случился перенос. Код для проставления неопределённости $$ c=2 $$ (может быть распараллелен): ```python c_bits = [ 2 if s_bits[i] == 1 else ci for i, ci in enumerate(c_bits) ] ``` ### Ассоциативность Очень важное наблюдение состоит в том, что данными числами $$ c $$ можно характеризовать блоки из последовательных разрядов, причём, их совокупное $$ c_{ij} \in \{0,1,2\} $$ от разряда i до j можно вычислять в любом порядке (не нарушая последовательность). Напомним, что в математике ассоциативностью называется следующее правило раскрытия скобок, в частности, для сложения: $$ (x+y)+z = x+(y+z) $$. Действительно, целый блок бит от i до j можно охарактеризовать $$ c=0 $$, если он точно не производит перенос, $$ c=1 $$ - если точно производит и $$ c=2 $$ - если производит при условии, что перед ним случился перенос. Очевидно правило по объединению блоков. Если для левого блока $$ c_{\text{left}} \in \{0,1\} $$, то берётся именно это значение для результата $$ c_{\text{result}} $$. Если же для левого блока $$ c_{\text{left}} = 2 $$, то в качестве результата берётся значение правого блока $$ c_{\text{right}} $$. Объединять можно только соседние блоки. $$ c_{\text{result}} = \begin{cases} c_{\text{left}}, c_{\text{left}} \in \{0,1\} \\ c_{\text{right}}, c_{\text{left}} = 2 \\ \end{cases} $$ Можно доказать, что данная операция объединения блоков является ассоциативной, что позволяет вычислять значения соседних блоков независимо (параллельно), в конце объединив их одной операцией. ### Дерево из блоков Мы можем построить один блок из всех $$ 2^n $$ элементов. Под ним два последовательных блока по $$ 2^{n-1} $$ элементов. Под каждым из них по два блока из $$ 2^{n-2} $$ элементов. И так далее, пока не дойдём до блоков, состоящих из единственного бита. Удобно вычислять $$ c $$ для каждого блока в обратном порядке, объединяя мелкие в более крупные. Учитывая параллельность объединения соседних блоков, построение такого дерева делается за время $$ O(\log N) = O(n) $$. Код для вычисления переносов блоков (внутренний цикл может быть распараллелен). Здесь для оптимизации памяти значения `c` пишутся в тот же массив, подразумевая, что `c[0]` теперь относится к корневому блоку размером $$ 2^n $$, `c[n/2]` - ко второму блоку длиной `n/2` (первый левый блок нам для вычислений на самом деле не нужен) и тд. ```python d = 1 while d < n: for i in range(0, n, 2*d): if c_bits[i] == 2: c_bits[i] = c_bits[i+d] d *= 2 ``` ### Вычисление неопределённости $$ c=2 $$ Когда дерево построено, некоторые блоки могут иметь $$ c=2 $$, что считается неопределённостью. Древовидная структура позволяет теперь уже сверху вниз разрешить эту неопределённость, установив её либо в $$ c=0 $$, либо в $$ c=1 $$ в зависимости от значения соседнего справа блока. Если соседнего блока нет, то ставится $$ c=0 $$. В остальных случаях удобно в качестве соседнего блока рассматривать только те, у которых длина вдвое больше, так как они уже были вычислены на предыдущем шаге. Таким образом разрешение неопределённости также происходит за время $$ O(\log N) = O(n) $$. Код для разрешения неопределённости (внутренний цикл может быть распараллелен): ```python d = n while d > 0: for i in range(d, n, 2*d): if c_bits[i] == 2: c_bits[i] = c_bits[i+d] if i+d < n else 0 d //= 2 ``` ### Финализация Теперь осталось эти $$ c $$ добавить к уже вычисленным разрядам побитовой суммы со сдвигом на 1. Это и будет итог. Код (может быть распараллелен): ```python s_bits = [ si ^ c_bits[i+1] if i < n-1 else si for i, si in enumerate(s_bits) ] ``` ## Промежуточный результат Ключевым пунктом оказалась возможность объединять переносы в блоки, причём, операция объединения обладала свойством ассоциативности. Это позволило сделать ряд операций параллельными, а также без дополнительного расхода памяти. Рассматриваемый подход можно обобщить для произвольных вычислений, которые точно так же можно ассоциативно объединять в блоки и разрешать уже в построенном дереве. Ссылка на полный пример кода: [carry-propagation.py](/backend/media/files/carry-propagation.py) ## Обобщение результата Теперь, проводя аналогии, обобщим результат для вычисления с произвольной ассоциативной функцией $$ f $$. Пусть у нас есть массив входных данных $$ x $$ длиной $$ N $$, индексация начинается с 0 и идёт от фактически последнего элемента к начальному. Мы хотим построить кумулятивный массив $$ y $$, который является некоторой агрегацией всех $$ x $$ от начала (последнего элемента по выбранной нумерации) до конкретного индекса включительно. Такая непривычная индексация выбрана для простоты описания алгоритмов ниже. Значение $$ y_i $$ вычисляется рекурентно: $$ y_i = f(x_i, y_{i+1}), i = 0,1,2,...,(N-2) $$ Ассоциативность функции $$ f $$ на математическом языке записывается как: $$ f(f(a,b),c) = f(a,f(b,c)) $$ Например, если мы имеем дело с кумулятивной суммой, то $$ f(a,b) = a+b $$. Оно обладает ассоциативностью $$ (a+b)+c = a+(b+c) $$. Тогда: $$ y_i = x_i + y_{i+1}, i = 0,1,2,...,(N-2) $$ $$ y_i = \sum_{j=i}^{N-1} x_j, i = 0,1,2,...,(N-2) $$ В примере выше функция $$ f $$ была связана с накоплением переноса, что также было ассоциативно. Переходим к самому алгоритму. На первом шаге мы копируем массив $$ x $$ в $$ y $$. Далее, объединяем $$ y $$ в блоки, чтобы добиться параллельности. Сначала мы берём блоки по два элемента, применяем к ним $$ f $$ и пишем результат в левый элемент, сохраняя в нем данные о всём блоке из двух элементов. Мы фактически стираем старое значение, но нам оно больше не потребуется. Далее мы попарно объединяем блоки по два элемента и результат пишем в первый элемент левого блока. Повторяем процедуру, пока весь массив не станет единым блоком. Если $$ N $$ не является степенью двойки, то всегда можно слева дополнить нулями. Пример кода (внутренний цикл может быть распараллелен): ```python d = 1 while d < n: for i in range(0, n, 2*d): y[i] = f(y[i], y[i+d]) d *= 2 ``` После этого мы делаем разрешение, которое заключается в объединении результатов каждого блока со всеми элементами справа, таким образом в самом конце получая результаты агрегации для каждого элемента. Мы также можем несколько расширить задачу, положив, что старт агрегации начинается не с нуля, а с некоторого значения $$ y_{\text{init}} $$. Пример кода (внутренний цикл может быть распараллелен): ```python d = n while d > 0: for i in range(d, n, 2*d): y[i] = f(y[i], y[i+d] if i+d < n else y_init) d //= 2 ``` Данный подход в теории алгоритмов называется префиксным деревом.

© Alexander Khlebushchev 2015 - now