Мой любимый алгоритм: нахождение медианы за линейное время +44



image

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

Нахождение медианы за O(n log n)


Самым прямолинейным способом нахождения медианы является сортировка списка и выбор медианы по её индексу. Самая быстрая сортировка сравнением выполняется за O(n log n), поэтому от неё зависит время выполнения1, 2.

def nlogn_median(l):
    l = sorted(l)
    if len(l) % 2 == 1:
        return l[len(l) / 2]
    else:
        return 0.5 * (l[len(l) / 2 - 1] + l[len(l) / 2])

У этого способа самый простой код, но он определённо не самый быстрый.

Нахождение медианы за среднее время O(n)


Следующим нашим шагом будет нахождение медианы в среднем за линейное время, если нам будет везти. Этот алгоритм, называемый «quickselect», разработан Тони Хоаром, который также изобрёл алгоритм сортировки с похожим названием — quicksort. Это рекурсивный алгоритм, и он может находить любой элемент (не только медиану).

  1. Выберем индекс списка. Способ выбора не важен, на практике вполне подходит и случайный. Элемент с этим индексом называется опорным элементом (pivot).
  2. Разделим список на две группы:
    1. Элементы меньше или равные pivot, lesser_els
    2. Элементы строго большие, чем pivot, great_els
  3. Мы знаем, что одна из этих групп содержит медиану. Предположим, что мы ищем k-тый элемент:
    • Если в lesser_els есть k или больше элементов, рекурсивно обходим список lesser_els в поисках k-того элемента.
    • Если в lesser_els меньше, чем k элементтов, рекурсивно обходим список greater_els. Вместо поиска k мы ищем k-len(lesser_els).

Вот пример алгоритма, выполняемого для 11 элементов:

Возьмём представленный ниже список. Мы хотим найти медиану.
l = [9,1,0,2,3,4,6,8,7,10,5]
len(l) == 11, поэтому мы ищем шестой наименьший элемент 
Сначала нам нужно выбрать опорный элемент (pivot). Мы случайным образом выбираем индекс 3. 
Значение элемента с этим индексом равно 2.

Разбиваем список на группы согласно pivot:
[1,0,2], [9,3,4,6,8,7,10,5]
Нам нужен шестой элемент. 6-len(left) = 3, поэтому нам нужен
третий наименьший элемент в правом массиве 

Теперь мы ищем третий наименьший элемент в следующем массиве:
[9,3,4,6,8,7,10,5]
Мы случайным образом выбираем индекс, который будет нашим pivot. 
Мы выбрали индекс 2, значение в котором равно l[2]=6

Разбиваем на группы согласно pivot:
[3,4,5,6] [9,7,10]
Нам нужен третий наименьший элемент, поэтому мы знаем, что это
третий наименьший элемент в левом массиве

Теперь мы ищем третий наименьший в следующем массиве:
[3,4,5,6]
Мы случайным образом выбираем индекс, который будет нашим pivot.
Мы выбрали индекс 1, значение в котором равно l[1]=4
Разбиваем на группы согласно pivot:
[3,4] [5,6]
Нам нужен третий наименьший элемент, поэтому мы знаем, что это
наименьший элемент в правом массиве.

Теперь мы ищем наименьший элемент в следующем массиве:
[5,6]

На этом этапе у нас есть базовый вариант, выбирающий наибольший
или наименьший элемент на основании индекса.
Нам нужен наименьший элемент, то есть 5.
return 5

Чтобы найти с помощью quickselect медиану, мы выделим quickselect в отдельную функцию. Наша функция quickselect_median будет вызывать quickselect с нужными индексами.

import random
def quickselect_median(l, pivot_fn=random.choice):
    if len(l) % 2 == 1:
        return quickselect(l, len(l) / 2, pivot_fn)
    else:
        return 0.5 * (quickselect(l, len(l) / 2 - 1, pivot_fn) +
                      quickselect(l, len(l) / 2, pivot_fn))


def quickselect(l, k, pivot_fn):
    """
    Выбираем k-тый элемент в списке l (с нулевой базой)
    :param l: список числовых данных
    :param k: индекс
    :param pivot_fn: функция выбора pivot, по умолчанию выбирает случайно
    :return: k-тый элемент l
    """
    if len(l) == 1:
        assert k == 0
        return l[0]

    pivot = pivot_fn(l)

    lows = [el for el in l if el < pivot]
    highs = [el for el in l if el > pivot]
    pivots = [el for el in l if el == pivot]

    if k < len(lows):
        return quickselect(lows, k, pivot_fn)
    elif k < len(lows) + len(pivots):
        # Нам повезло и мы угадали медиану
        return pivots[0]
    else:
        return quickselect(highs, k - len(lows) - len(pivots), pivot_fn)

В реальном мире Quickselect отлично себя проявляет: он почти не потребляет лишних ресурсов и выполняется в среднем за O(n). Давайте докажем это.

Доказательство среднего времени O(n)


В среднем pivot разбивает список на две приблизительно равных части. Поэтому каждая последующая рекурсия оперирует с 1?2 данных предыдущего шага.

$C=n+\frac{n}{2}+\frac{n}{4}+\frac{n}{8}+…=2n=O(n)$


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

Quickselect даёт нам линейную скорость, но только в среднем случае. Что, если нас не устраивает среднее, и мы хотим гарантированного выполнения алгоритма за линейное время?

Детерминированное O(n)


В предыдущем разделе я описал quickselect, алгоритм со средней скоростью O(n). «Среднее» в этом контексте означает, что в среднем алгоритм будет выполняться за O(n). С технической точки зрения, нам может очень не повезти: на каждом шаге мы можем выбирать в качестве pivot наибольший элемент. На каждом этапе мы сможем избавляться от одного элемента из списка, и в результате получим скорость O(n^2), а не O(n).

С учётом этого, нам нужен алгоритм для подбора опорных элементов. Нашей целью будет выбор за линейное время pivot, который в худшем случае удаляет достаточное количество элементов для обеспечения скорости O(n) при использовании его вместе с quickselect. Этот алгоритм был разработан в 1973 году Блумом (Blum), Флойдом (Floyd), Праттом (Pratt), Ривестом (Rivest) и Тарьяном (Tarjan). Если моего объяснения вам не хватит, то можете изучить их статью 1973 года. Вместо того, чтобы описывать алгоритм, я подробно прокомментирую мою реализацию на Python:

def pick_pivot(l):
    """
    Выбираем хорошй pivot в списке чисел l
    Этот алгоритм выполняется за время O(n).
    """
    assert len(l) > 0

    # Если элементов < 5, просто возвращаем медиану
    if len(l) < 5:
        # В этом случае мы возвращаемся к первой написанной нами функции медианы. 
        # Поскольку мы выполняем её только для списка из пяти или менее элементов, она не 
        # зависит от длины входных данных и может считаться постоянным
        # временем.
        return nlogn_median(l)

    # Сначала разделим l на группы по 5 элементов. O(n)
    chunks = chunked(l, 5)

    # Для простоты мы можем отбросить все группы, которые не являются полными. O(n)
    full_chunks = [chunk for chunk in chunks if len(chunk) == 5]


    # Затем мы сортируем каждый фрагмент. Каждая группа имеет фиксированную длину, поэтому каждая сортировка
    # занимает постоянное время. Поскольку у нас есть n/5 фрагментов, эта операция 
    # тоже O(n)
    sorted_groups = [sorted(chunk) for chunk in full_chunks]

    # Медиана каждого фрагмента имеет индекс 2
    medians = [chunk[2] for chunk in sorted_groups]

    # Возможно, я немного повторюсь, но я собираюсь доказать, что нахождение
    # медианы списка можно произвести за доказуемое O(n).
    # Мы находим медиану списка длиной n/5, поэтому эта операция также O(n)
    # Мы передаём нашу текущую функцию pick_pivot в качестве создателя pivot алгоритму
    # quickselect. O(n)
    median_of_medians = quickselect_median(medians, pick_pivot)
    return median_of_medians

def chunked(l, chunk_size):
    """Разделяем список `l` на фрагменты размером `chunk_size`."""
    return [l[i:i + chunk_size] for i in range(0, len(l), chunk_size)]

Давайте докажем, что медиана медиан является хорошим pivot. Нам поможет, если мы представим визуализацию нашего алгоритма выбора опорных элементов:


Красным овалом обозначены медианы фрагментов, а центральным кругом — медиана медиан. Не забывайте, мы хотим, чтобы pivot разделял список как можно ровнее. В худшем возможном случае каждый элемент в синем прямоугольнике (слева вверху) будет меньше или равен pivot. Верхний правый прямоугольник содержит 3?5 половины строк — 3/5*1/2=3/10. Поэтому на каждом этапе мы избавляемся по крайней мере от 30% строк.

Но достаточно ли нам отбрасывать 30% элементов на каждом этапе? На каждом этапе наш алгоритм должен выполнять следующее:

  • Выполнять работу O(n) по разбиению элементов
  • Для рекурсии решать одну подзадачу размером в 7?10 от исходной
  • Для вычисления медианы медиан решать одну подзадачу размером с 1?5 от исходной

В результате мы получаем следующее уравнение полного времени выполнения T(n):

$T(n)=T(\frac{n}{5})+7T(\frac{n}{10})+n$


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

Подводим итог


У нас есть quickselect, алгоритм, который находит медиану за линейное время при условии наличия достаточно хорошей опорного элемента. У нас есть алгоритм медианы медиан, алгоритм O(n) для выбора опорного элемента (который достаточно хорош для quickselect). Соединив их, мы получили алгоритм нахождения медианы (или n-ного элемента в списка) за линейное время!

Медианы за линейное время на практике


В реальном мире почти всегда достаточно случайного выбора медианы. Хотя подход с медианой медиан всё равно выполняется за линейное время, на практике его вычисление длится слишком долго. В стандартной библиотеке C++ используется алгоритм под названием introselect, в котором применено сочетание heapselect и quickselect; предел его выполнения O(n log n). Introselect позволяет использовать обычно быстрый алгоритм с плохим верхним пределом в сочетании с алгоритмом, который медленнее на практике, но имеет хороший верхний предел. Реализации начинают с быстрого алгоритма, но возвращаются к более медленному, если не могут выбрать эффективные опорные элементы.

В завершение приведу сравнение элементов, используемых в каждой из реализаций. Это не скорость выполнения, а общее количество элементов, которые рассматривает функция quickselect. Здесь не учитывается работа по вычислению медианы медиан.

image

Именно этого мы и ожидали! Детерминированный опорный элемент почти всегда рассматривает при quickselect меньшее количество элементов, чем случайный. Иногда нам везёт и мы угадываем pivot с первой попытки, что проявляется как впадины на зелёной линии. Математика работает!



  1. Это может стать интересным применением поразрядной сортировки (radix sort), если вам нужно найти медиану в списке целых чисел, каждое из которых меньше 232.
  2. На самом деле в Python используется Timsort, впечатляющее сочетание теоретических пределов и практической скорости. Заметки о списках в Python.

Вы можете помочь и перевести немного средств на развитие сайта



Комментарии (32):

  1. mikhaelkh
    /#10626472 / +1

    Рекуррентное соотношение должно быть T(n) <= T(0.2 * n) + T(0.7 * n) + O(n), нам важно, что (0.2 * n + 0.7 * n) / n = 0.9 < 1. Использовать то, что T(n) = O(n) при доказательстве этого же факта — это facepalm.

    Не так уж просто доказать, почему это равно O(n)
    На самом деле очень просто — подберём константу c такую, что для всех n выполняется T(n) <= T(0.2 * n) + T(0.7 * n) + c * n или T(n) <= 10 * c * n , тогда T(n) <= 10 * c * n по индукции.

  2. yorko
    /#10626552 / -1

    Кажется, что такой выбор опорного элемента (медиана из медиан) должен помогать и quicksort-у, поскольку при каждом разделении получаются более равные по длине 2 списка. Однако для quicksort доказано, что оптимальная стратегия выбора опорного элемента – случайная.

    • splav_asv
      /#10626660

      Можно ссылочку на доказательство?

      • yorko
        /#10626714

        Я это помню по курсу Тима Рафгардена. В видео "Choosing a good pivot" лектор про это говорит и доказывает, что так получится среднее время O(n logn), но формальное доказательство, что случайная стратегия – наилучшая, наверное в его книге надо искать.

        • splav_asv
          /#10626726

          Обычно везде говорится, что этастратегия, не имеет детерминированной «kill sequence» (данные при которых сваливаемся в O(n^2)). На разных данных оптимальный выбор будет отличаться. Если сортировать последовательность равномерно распределённых случайных величин(или близкую к ней), то вообще можно брать первый попавшийся pivot: среднее не пострадает, а константа лучше.

          Это я всё к чему — оптимальной стратегии нет. Есть размен ухудшения константы в среднем на отсутствие детерминированного худшего случая. А дальше зависит от данных и требований.

          • yorko
            /#10626742

            Пожалуй, про "оптимальную стратегию" – это громко сказано, запомнилось просто, как акцентировано это Рафгарден отмечал. В книге Сэджвика "Algorithms" про quicksort нашел только утверждение (с. 295, Proposition L), что в худшем случае quicksort делает ~ N^2/2 сравнений (очевидно), но случайное перемешивание на каждом шаге снижает число сравнений обратно до ~ n log n. Но да, это далеко не то же самое, о чем я говорил в начале.

          • chersanya
            /#10626866

            Нужно как минимум понимать, что при рассмотрении quicksort со случайным выбором элемента среднее берётся не по входным данным, а по реализациям случайных решений в ходе алгоритма. Это и значит, что не существует набора данных, дающих n^2 — среднее количество операций для любого массива будет O(n log n).

            Можно даже посчитать вероятность того, что алгоритм будет выполняться например в 5 раз медленнее, чем в среднем. Эта вероятность получается весьма малой, особенно при больших n — поэтому на практике имея адекватный генератор случайных чисел нет смысла что-то оптимизировать по сравнению со случайным выбором элемента.

            • splav_asv
              /#10626880

              Не вижу противоречия.

              • chersanya
                /#10626934

                Например:

                Есть размен ухудшения константы в среднем на отсутствие детерминированного худшего случая.

                «детерминированного худшего случая» нет ни в алгоритме со случайным выбором pivot, ни в описанном в статье.

                • splav_asv
                  /#10626950

                  Ну так я и не говорю, что в этих вариантах есть худший случай. Случайный выбор можно сделать быстрее O(1), данный алгоритм O(n). Первый не влияет на константу по идее, второй влияет. Для первого в каждом конкретном случае есть плохой вариант. Для второго, исключая массив повторяющихся величин, нет.
                  В данном случае такой размен.

                  • chersanya
                    /#10626982

                    Что-то я вас вообще не понял.

                    Случайный выбор можно сделать быстрее O(1)

                    Алгоритмов быстрее О(1) не бывает.
                    Для первого в каждом конкретном случае есть плохой вариант.

                    Поясните, что имеете в виду под случаем и что под вариантом?

                    • splav_asv
                      /#10627016

                      можно сделать быстрее O(1)

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

                      В данном контексте случай — реализация генератора псевдослучайных чисел. Вариант — набор данных, который нужно подать на вход quicksort (зная состояние ГПСЧ), чтобы свалиться в O(n^2).

    • MikailBag
      /#10626678

      Насколько я помню, quicksort используется исключительно из-за низкой константы, а при использовании "безопасного" поиска медианы он работает медленнее других алгоритмов сортировки.

    • chersanya
      /#10626690

      Причём тут случайная стратегия, если речь о детерминированных алгоритмах? Детерминированный quicksort стандартно пишется как раз с таким же выбором опорного элемента, как и в алгоритме поиска медианы. Единственное отличие quicksort от quickselect — то, что он обрабатывает обе части массива, а не только одну из них.

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

      • yorko
        /#10626718

        Вывод по графику в статье


        Детерминированный опорный элемент почти всегда рассматривает при quickselect меньшее количество элементов, чем случайный

        Мой вопрос (или мысли вслух): работает ли то же самое для quicksort. Будет ли детерминированный quicksort with median pivots рассматривать меньше элементов, чем недетерминированный quicksort with random pivot.

        • chersanya
          /#10626864

          Да, это тоже верно. Однако, как отмечено и в статье, это не говорит о более быстром выполнении алгоритма. Всегда или почти всегда на практике быстрее работает случайный выбор.

  3. nikitadanilov
    /#10626554

    А почему pivot выбирается именно как элемент массива с некоторым индексом? В случае quicksort это понятно, но для quickselect можно выбрать любое значение (не обязательно из массива). Например, брать среднее (O(n)), тогда массив будет гарантированно хорошо делиться.

    • andreybotanic
      /#10626614

      Вот контрпример: массив, состоящий из одного числа 100 и 99 нулей. Среднее значение равно 1. В итоге в одной части окажутся все нули, а во второй всего одно число — наша сотня. Не самой удачное разделение, не так ли?

      • dmih
        /#10626788

        По приколу спрошу: ваша версия деления такого массива примерно пополам? :)

      • netch80
        /#10629588

        Это только на один цикл вглубь. На следующем окажется, что все элементы одинаковы.

        (Хотя, чтобы корректно выйти из этой ситуации, исходный алгоритм должен делить вход не на две части, а на три, или же пополамить набор элементов, равных опорному. Для quicksort некоторые методы деления отрезка специально заточены на оптимизацию такого варианта. Это можно перенести и на quickselect, но явно.)

  4. nikitazvonarev
    /#10626876

    На самом деле, нас могут интересовать не только медианы, но и квартили и вообще произвольные квантили, короче говоря, k-е порядковые статистики (т.е. k-ый элемент в отсортированной последовательности). quickselect (который обрубленный quicksort) спокойно обобщается на этот случай с сохранением асимптотики среднего времени, а детерминированный подход можно ли обобщить?

    В статье это никак не раскрыто.

    • chersanya
      /#10626930

      Описанный в статье алгоритм это и есть детерминированный quickselect, и он абсолютно так же обобщается.

    • khim
      /#10627030

      Вы наверное какую-ту другую статью читали. Потому что то, что описано в статье сначала рассматривает нахождение медианы как нахождение k-й порядковой статистики с номером N/2, а потом, рекурсивно решает задачу нахождения k-й порядковой статистики.

      Обратите внимание что один шаг алгоритма: он заключается в том, что мы отбрасываем либо 30% «слишком больших» чисел, либо 30% слишком маленьких, после чего ищем… медиану? Фиг вам: так как «перекосили» массив, то нам теперь нужна не медиана, нам нужен элемент с другим (впрочем легко высчитываемым) номером.

      Куда вы это обобщать собрались?

  5. snowSTAFF
    /#10626878

    Ваши алгоритмы. Все на Python 3.6.3
    Что-то мне не везло. Это все на списке из 10 000 рандомных от 0 до 1000
    Причем, 3-я функция работает с ошибкой.
    Для 1-ой:
    2.47 ms ± 39.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
    Для 2-ой:
    7.92 ms ± 159 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
    Для 3-ей:
    11.2 ms ± 73.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

    • alexeykuzmin0
      /#10626966

      Дык список слишком короткий, чтобы найти различие между 1 и 2. А 3 и должна работать дольше, она только сравнений меньше делает (ну и имеет теоретическую гарантию того, что не будет тупить — для варианта 2 можно подобрать контрпример).

  6. sergio_nsk
    /#10626920

    В C++ уже всё сделано для вас:


    #include <iostream>
    #include <vector>
    #include <algorithm>
    #include <functional>
    
    int main()
    {
        std::vector<int> v{5, 6, 4, 3, 2, 6, 7, 9, 3};
    
        std::nth_element(v.begin(), v.begin() + v.size()/2, v.end());
        std::cout << "The median is " << v[v.size()/2] << '\n';
    }

    Output:


    The median is 5

  7. Mingun
    /#10627272

    Вот здесь же:


    # Для простоты мы можем отбросить все группы, которые не являются полными. O(n)
    full_chunks = [chunk for chunk in chunks if len(chunk) == 5]

    можно упростить до:


    last_chunk = chunks[len(chunks) - 1]
    full_chunks = chunks
    if len(last_chunk) < 5:
        full_chunks.pop()

    и будет О(1). Разбивка ведь оригинального списка ведётся так, что меньше 5-ти элементов может быть только в последнем подсписке.


    Кроме того, визуализация алгоритма неправильная же. В первом и последнем подсписках выбраны вовсе не медианы! В первом подсписке медиана 41, а не 99, а во втором — 65 вместо 116. Правда, это легко исправляется дописыванием единичек спереди к последним числам в подсписках, но всё равно странно приводить заведомо неправильную картинку в качестве визуализации алгоритма.

    • amakhrov
      /#10627766

      можно упростить
      и будет О(1)

      Хорошее замечание. Стоит уточнить, что хотя асимптотика всего алгоритма от этого не меняется, мы уменьшаем константу.

    • rcoh
      /#10630124 / +1

      Привет, я автор. Извините за автоматический перевод. Спасибо, что заметили эту ошибку! Это был ужасный надзор, и я исправил его.