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

https://rcoh.me/posts/linear-time-median-finding/
  • Перевод
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 разбивает список на две приблизительно равных части. Поэтому каждая последующая рекурсия оперирует с 12 данных предыдущего шага.

$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. Верхний правый прямоугольник содержит 35 половины строк — 3/5*1/2=3/10. Поэтому на каждом этапе мы избавляемся по крайней мере от 30% строк.

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

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

В результате мы получаем следующее уравнение полного времени выполнения 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.
Поделиться публикацией
Ой, у вас баннер убежал!

Ну. И что?
Реклама
Комментарии 34
  • +3
    Рекуррентное соотношение должно быть 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 по индукции.
    • 0

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

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

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

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

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

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

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

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

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

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

                        Поясните, что имеете в виду под случаем и что под вариантом?
                        • 0
                          можно сделать быстрее O(1)

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

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

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

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

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

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


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

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

                • +1
                  Да, это тоже верно. Однако, как отмечено и в статье, это не говорит о более быстром выполнении алгоритма. Всегда или почти всегда на практике быстрее работает случайный выбор.
            • 0
              А почему pivot выбирается именно как элемент массива с некоторым индексом? В случае quicksort это понятно, но для quickselect можно выбрать любое значение (не обязательно из массива). Например, брать среднее (O(n)), тогда массив будет гарантированно хорошо делиться.
              • 0
                Совсем не обязательно.
                • 0
                  del (уже выше написали)
                  • +3
                    Вот контрпример: массив, состоящий из одного числа 100 и 99 нулей. Среднее значение равно 1. В итоге в одной части окажутся все нули, а во второй всего одно число — наша сотня. Не самой удачное разделение, не так ли?
                    • 0
                      По приколу спрошу: ваша версия деления такого массива примерно пополам? :)
                    • 0
                      Это только на один цикл вглубь. На следующем окажется, что все элементы одинаковы.

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

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

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

                        Куда вы это обобщать собрались?
                      • 0
                        Ваши алгоритмы. Все на 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)
                        • 0
                          Дык список слишком короткий, чтобы найти различие между 1 и 2. А 3 и должна работать дольше, она только сравнений меньше делает (ну и имеет теоретическую гарантию того, что не будет тупить — для варианта 2 можно подобрать контрпример).
                        • –6

                          В 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
                          • 0

                            Вот здесь же:


                            # Для простоты мы можем отбросить все группы, которые не являются полными. 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. Правда, это легко исправляется дописыванием единичек спереди к последним числам в подсписках, но всё равно странно приводить заведомо неправильную картинку в качестве визуализации алгоритма.

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

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

                              • +1
                                Привет, я автор. Извините за автоматический перевод. Спасибо, что заметили эту ошибку! Это был ужасный надзор, и я исправил его.
                              • 0
                                # Затем мы сортируем каждый фрагмент. Каждая группа имеет фиксированную длину, поэтому каждая сортировка
                                # занимает постоянное время. Поскольку у нас есть n/5 фрагментов, эта операция
                                # тоже O(n)

                                Мне не понятно почему это время линейное. Ну, да — у нас n/5 фрагментов, как из этого следует что сортировка n/5 фрагментов по 5 элементов это O(n)?
                                • +1
                                  Потому что время сортировки одного фрагмента ограничено сверху. Ну пусть отсортировать 5 элементов требует до 100500 операций (очень плохая, очень медленная сортировка), тогда отсортировать n/5 фрагментов — это не более 20100 * n операций. 20100 * n — это O(n)

                                Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                                Самое читаемое