Ранг в списке за линейную работу

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

Как думаешь, где тратится работа

В уроке о применениях скана прыжки ранжировали список из тысячи элементов за десять раундов и 11 000 операций. Последовательный обход справился бы за тысячу. Восемь ядер против одного: при работе n · log₂ n на миллионе элементов прыжки делают двадцать миллионов обновлений, по два с половиной миллиона на ядро - в два с половиной раза больше, чем один обход на одном ядре.

Прежде чем читать дальше, ответь себе: какие элементы в прыжках трудятся зря? И можно ли их убрать из работы, не потеряв логарифмической глубины целиком?

На что это похоже

Перекличка в длинной колонне, где каждый знает только того, кто стоит за ним. Можно всем разом перекрикиваться через одного, через двух, через четырёх - так быстро, но кричат все и на каждом круге. А можно иначе: часть людей ненадолго выходит из колонны, сказав соседу впереди «прибавь мою единицу». Колонна стала короче, её перекличку провести дешевле, а вышедшие потом встают на место и узнают свой номер у соседа сзади. Главное условие: из колонны нельзя выходить двоим соседям сразу, иначе некому будет передать единицу.

Механизм: сжать, решить, собрать

Список, как и раньше, лежит в массиве вперемешку: sled[i] - следующий, у хвоста sled указывает на него самого, начальный rang у всех 1, у хвоста 0. Нужен для каждого элемента ранг - сумма rang от него до хвоста. Это обобщение (у Сандерса и соавторов, раздел 3.3) удобно тем, что после сжатия начальные ранги уже не единицы, а алгоритм тот же.

Выбор. Каждый элемент бросает монету, хвост всегда выкидывает решку. Уходит элемент, у которого орёл, а у следующего - решка. Двое соседей уйти не могут: у ушедшего следующий выбросил решку, значит сам он не ушёл. Предыдущего проверять не нужно - это следует из его собственного правила.

Сжатие. Оставшиеся получают новые номера исключающим сканом флагов «остаюсь» - как при упаковке. Если следующий за оставшимся ушёл, оставшийся прибавляет себе ранг ушедшего и берёт его указатель: rang += rang[sled], sled = sled[sled]. Это один прыжок, но только там, где он нужен.

Рекурсия и сборка. Короткий список решается тем же способом, пока длина не станет меньше порога B, а там - прыжками по указателям. Оставшиеся забирают ранг из решения, ушедшие делают rang += rang[sled]: их следующий остался и ранг уже знает.

Почему четверть. Монеты независимы, поэтому вероятность «у меня орёл и у следующего решка» - произведение: ½ · ½ = ¼. Ожидаемое число ушедших по линейности матожидания из главы о вероятностном анализе - сумма этих четвертей по всем элементам, при n ≥ 2 точно n/4: у n−2 элементов вероятность ухода 1/4, у предпоследнего 1/2, у хвоста 0. При одном хвосте никто не уходит.

Почему работа линейна. Каждый уровень проходит свои элементы дважды, а элементов на уровне в среднем ¾ от прошлого: n + ¾n + (¾)²n + … не больше n / (1 − ¾) = 4n. Уровней - сколько раз n делится на 4/3, пока не станет меньше B, то есть log по основанию 4/3 от n/B. На каждом уровне скан стоит O(log n) шагов, отсюда и глубина. Сандерс берёт B = n / log n: прыжки на остатке тогда стоят B · log B = O(n), а уровней - O(log log n).

Правило

Работоэффективный ранг в списке: на каждом уровне удалить независимое множество (уходит тот, у кого орёл, а у следующего решка, - в среднем n/4), сшить оставшихся сканом, решить короткий список рекурсивно, ушедшим добавить ранг следующего. При пороге B = n / log n с прыжками в основании - ожидаемые работа O(n) и глубина O(log n · log log n). Неэффективный алгоритм допустим, когда задача уже сжата настолько, что его лишняя работа не видна.

Разбор: удвоение против сжатия на 65 536 элементах

Программа ранжирует случайно перемешанный список обоими способами, сверяет ответ с настоящими рангами и считает работу (обновлений элемента) и глубину (шагов, которые нельзя сделать одновременно: раунды прыжков плюс шаги сканов на каждом уровне).

import math
import random


def udvoenie(sled, rang, uchet):
    # прыжки по указателям; у хвоста sled = сам он и rang = 0
    sled, rang = sled[:], rang[:]
    n = len(sled)
    while any(sled[i] != sled[sled[i]] for i in range(n)):
        rang = [rang[i] + rang[sled[i]] for i in range(n)]
        sled = [sled[sled[i]] for i in range(n)]
        uchet["работа"] += n
        uchet["раундов"] += 1
    return rang


def isr(sled, rang, B, rng, uchet):
    n = len(sled)
    if n <= B:
        return udvoenie(sled, rang, uchet)
    uchet["уровней"] += 1
    orel = [rng.random() < 0.5 and sled[i] != i for i in range(n)]
    # уходит тот, у кого орёл, а у следующего решка
    ujti = [orel[i] and not orel[sled[i]] for i in range(n)]
    nomer, m = [], 0                  # исключающий скан «остаётся»
    for i in range(n):
        nomer.append(m)
        m += not ujti[i]
    sled2, rang2 = [0] * m, [0] * m
    for i in range(n):
        if not ujti[i]:
            r = sled[i]
            j = nomer[i]
            rang2[j] = rang[i] + (rang[r] if ujti[r] else 0)
            sled2[j] = nomer[sled[r]] if ujti[r] else nomer[r]
    uchet["работа"] += 2 * n          # сжатие и обратная сборка
    uchet["шагов скана"] += math.ceil(math.log2(n))
    r2 = isr(sled2, rang2, B, rng, uchet)
    otvet = [0] * n
    for i in range(n):
        if not ujti[i]:
            otvet[i] = r2[nomer[i]]
    for i in range(n):
        if ujti[i]:
            otvet[i] = rang[i] + otvet[sled[i]]
    return otvet


rng = random.Random(5)
print("      n | удвоение: работа/n  раундов | ISR: работа/n"
      "  уровней  глубина")
for n in [1024, 16384, 65536]:
    poryadok = list(range(n))
    rng.shuffle(poryadok)
    sled = list(range(n))
    for k in range(n - 1):
        sled[poryadok[k]] = poryadok[k + 1]
    rang = [0 if sled[i] == i else 1 for i in range(n)]
    verno = [0] * n
    for k, v in enumerate(poryadok):
        verno[v] = n - 1 - k
    u1 = {"работа": 0, "раундов": 0}
    assert udvoenie(sled, rang, u1) == verno
    u2 = {"работа": 0, "раундов": 0, "уровней": 0, "шагов скана": 0}
    B = n // int(math.log2(n))
    assert isr(sled, rang, B, rng, u2) == verno
    gl = u2["шагов скана"] + u2["раундов"]
    print(f"{n:7} | {u1['работа'] / n:16.1f}  {u1['раундов']:7} |"
          f" {u2['работа'] / n:10.1f}  {u2['уровней']:7}  {gl:7}")

      n | удвоение: работа/n  раундов | ISR: работа/n  уровней  глубина
   1024 |             10.0       10 |        7.9        8       79
  16384 |             14.0       14 |        8.1       10      136
  65536 |             16.0       16 |        8.2       10      158

Работа на элемент у прыжков растёт как log₂ n: 10, 14, 16. У сжатия она стоит на месте около 8 - это и есть линейная работа: все уровни вместе содержат около 4n элементов, каждый проходится дважды, плюс немного на прыжки в основании. Уровней 10 на 65 536 элементах при пороге 4096: (¾)¹⁰ ≈ 0,056 < 1/16, а девяти мало - ровно как предсказывает доля ушедших. Цена видна в последней колонке: глубина 158 против 16. Сжатие не бесплатно - каждый уровень платит за свой скан, и при малом числе ядер выигрывает тот, кто меньше работает, а при тысячах ядер - тот, у кого меньше глубина.

Разбор: нельзя ли выбрасывать больше

Наивная мысль: монета нечестная, орёл с вероятностью 0,9 - уйдёт больше. Вторая наивная мысль: не проверять соседа, пусть уходят все с орлом - уйдёт половина. Проверим обе на списке из ста тысяч.

import random

rng = random.Random(8)
n = 100000
poryadok = list(range(n))
rng.shuffle(poryadok)
sled = list(range(n))
for k in range(n - 1):
    sled[poryadok[k]] = poryadok[k + 1]

print("p орла  ушло, доля   p(1 - p)")
for p in [0.1, 0.3, 0.5, 0.7, 0.9]:
    orel = [rng.random() < p and sled[i] != i for i in range(n)]
    ujti = sum(orel[i] and not orel[sled[i]] for i in range(n))
    print(f"{p:6}  {ujti / n:11.4f}   {p * (1 - p):8.4f}")

# без проверки соседа: уходит каждый, у кого орёл
orel = [rng.random() < 0.5 and sled[i] != i for i in range(n)]
sosedi = sum(orel[i] and orel[sled[i]] for i in range(n))
print("ушло бы:", sum(orel), " из них пар соседей:", sosedi)

p орла  ушло, доля   p(1 - p)
   0.1       0.0909     0.0900
   0.3       0.2099     0.2100
   0.5       0.2503     0.2500
   0.7       0.2107     0.2100
   0.9       0.0900     0.0900
ушло бы: 49769  из них пар соседей: 24754

Доля ушедших - это p · (1 − p): орёл у себя и решка у следующего. Произведение двух чисел с суммой 1 больше всего, когда они равны, поэтому честная монета даёт максимум ¼, а 0,9 даёт столько же, сколько 0,1, - девять процентов. Без проверки соседа уходит половина, но половина ушедших стоит парами: у элемента, следующий за которым ушёл, указатель ведёт к тоже ушедшему, и прибавить его ранг одним шагом уже нельзя - сшивка требует цепочки, то есть того же последовательного обхода.

Что ломается без этого

Прыжки на всём списке. Шестнадцать проходов по 65 536 элементам вместо восьми: на четырёх ядрах это медленнее одного ядра с простым циклом, хотя глубина и красивая.

Порог B слишком мал. Рекурсия до списков из пары элементов даёт log по основанию 4/3 от n уровней, каждый со сканом за log n шагов: глубина log² n. Работа при этом остаётся линейной - и никто не замечает, что выигрыш съела глубина.

Ушедшие соседи. Выбор «все с орлом» без проверки следующего ломает сшивку: ранги выходят неверными, а ошибка видна только на списке, где выпала пара орлов подряд.

Зачем это в работе

Связный список и дерево на видеокарте приходят в виде указателей, а не массива: рёбра эйлерова обхода дерева, цепочки слияния в компонентах связности, списки смежности после фильтрации. Прежде чем считать по ним сканы, их ранжируют, то есть превращают в массив. Приём «сжать независимое множество, решить меньшую задачу, достроить ответ» работает и шире: им же сжимают деревья (выбрасывают листья и вершины с одним ребёнком), а по той же многоуровневой схеме разбивает графы METIS - только сжимает он паросочетанием.

Теперь сам

1. Монеты списка из восьми элементов по порядку от головы к хвосту: О Р О О Р О Р Р (последняя - хвост). Какие элементы уйдут?
2. Сколько уровней сжатия нужно списку из миллиона элементов, чтобы дойти до порога B = 50 000, если каждый уровень оставляет ровно ¾?