Случайный порядок вставки: Делоне в ожидании и наименьшая окружность

Если ты знаешь, что такое рандомизированное инкрементное построение, как работает обратный анализ и почему наименьшая окружность строится за линейное ожидаемое время, а Делоне - за n · log n, листай до «Теперь сам».

Как думаешь, что решает порядок

Нужна наименьшая окружность, внутри которой все точки, - например, где поставить вышку, чтобы покрыть сорок домов наименьшим радиусом. Алгоритм простой: держим окружность для уже обработанных точек; новая точка внутри - ничего не делаем; снаружи - строим окружность заново, зная, что новая точка лежит на её границе.

Ответь себе до чтения: что будет, если подать точки от центра наружу - тогда новая точка почти всегда снаружи? А если в случайном порядке? Ответ: на сорока точках 3942 проверки против 203 - в двадцать раз больше, и при удвоении n второе число растёт почти впятеро - быстрее квадрата.

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

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

Механизм: починка, обратный анализ и цена

Инкрементное построение. Ответ A(i) для первых i объектов. Объект i + 1 либо не нарушает ответ - шаг почти бесплатный, - либо нарушает, и тогда ответ строят заново, но для задачи поменьше: известно, что новый объект войдёт в ответ особым образом. Для окружности это «новая точка на границе», и задача с одной точкой на границе решается тем же приёмом с двумя точками на границе, а с двумя - перебором третьей.

Обратный анализ. Сколько стоит шаг i в случайном порядке? Удобно смотреть назад: первые i точек - это случайное множество, а последняя вставленная - случайная точка из него. Шаг дорогой, только если эта точка была нужна ответу - если без неё окружность меньше. Окружность опирается на две или три точки, значит, вероятность ≤ 3 / i. Дорогой шаг - это починка на i точках. Её цена в ожидании O(i) - по тому же доводу на уровень ниже: внутри починки одна опорная точка уже известна, и её собственные дорогие шаги так же редки. В худшем случае починка дороже, но в ожидании цена шага 3/i · O(i) = O(1), а всех n шагов - O(n).

Делоне. Сначала добавляют три далёкие вершины - огромный треугольник, накрывающий все точки, - чтобы каждая новая точка оказалась внутри, а не за оболочкой; в конце их убирают вместе с рёбрами. Вставка точки: найти треугольник, где она лежит, разрезать его на три (точка ровно на ребре режет два треугольника на четыре) и переворачивать рёбра, ставшие нелегальными, - все они вокруг новой точки. Сколько переворотов? Столько, сколько рёбер у новой точки после вставки, минус три. Обратный анализ: последняя точка - случайная вершина Делоне из i точек, а сумма степеней всех вершин равна удвоенному числу рёбер, ≤ 2(3i − 6). Средняя степень меньше шести, ожидаемое число переворотов на вставку - O(1). Поиск треугольника идёт по истории: все треугольники, когда-либо существовавшие, хранятся со ссылками на те, что их заменили, и точку ведут от первого огромного треугольника вниз по ссылкам. Путь в ожидании O(log n). Итого O(n · log n) (Гибас, Кнут, Шарир, 1992).

Линейное программирование на плоскости - та же схема (Зайдель, 1991): ограничения вставляются в случайном порядке, нарушенное ограничение проходит через новый оптимум, и задача становится одномерной. Ожидаемое время O(n) при фиксированной размерности.

Правило

Рандомизированное инкрементное построение: перемешай вход, вставляй по одному, чини ответ при нарушении. Цену считает обратный анализ: вероятность, что последний объект был нужен ответу, умножить на цену починки. Наименьшая окружность - O(n) в ожидании (на неё опираются ≤ 3 точки), Делоне - O(n · log n) в ожидании (средняя степень < 6), линейное программирование на плоскости - O(n).

Разбор: от центра наружу против случайного

Окружность хранится своими опорными точками - двумя (диаметр) или тремя, - и проверка «точка внутри» точная: скалярное произведение для диаметра, определитель из прошлых уроков для трёх точек. Корней и центров не считаем вовсе.

import random
schet = [0]

def povorot(a, b, c):
    return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])

def v_kruge(a, b, c, d):
    m = []
    for p in [a, b, c]:
        dx = p[0] - d[0]
        dy = p[1] - d[1]
        m.append([dx, dy, dx * dx + dy * dy])
    return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
            - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
            + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]))

def vnutri(krug, p):
    # круг задан двумя точками диаметра или тремя точками на окружности;
    # граница считается внутренней; всё в целых
    schet[0] = schet[0] + 1
    if len(krug) == 2:
        a, b = krug
        dx = (p[0] - a[0]) * (p[0] - b[0])
        return dx + (p[1] - a[1]) * (p[1] - b[1]) <= 0
    a, b, c = krug
    if povorot(a, b, c) < 0:
        b, c = c, b
    return v_kruge(a, b, c, p) >= 0

def s_dvumya(P, q1, q2):
    krug = [q1, q2]
    for p in P:
        if not vnutri(krug, p):
            krug = [q1, q2, p]
    return krug

def s_odnoj(P, q):
    krug = [q, P[0]]
    for j in range(1, len(P)):
        if not vnutri(krug, P[j]):
            krug = s_dvumya(P[:j], P[j], q)
    return krug

pereschetov = [0]

def minimalnyj(P):
    krug = [P[0], P[1]]
    for i in range(2, len(P)):
        if not vnutri(krug, P[i]):
            pereschetov[0] = pereschetov[0] + 1
            krug = s_odnoj(P[:i], P[i])
    return krug
def nabor(n):
    tochki = []
    for i in range(n):
        x = random.randint(-1000, 1000)
        tochki.append([x, random.randint(-1000, 1000)])
    return tochki

def ot_centra(tochki):
    # неудачный порядок: от центра наружу, точка часто вне круга
    pary = []
    for p in tochki:
        pary.append([p[0] * p[0] + p[1] * p[1], p])
    pary.sort()
    otvet = []
    for para in pary:
        otvet.append(para[1])
    return otvet

random.seed(12)
for n in [20, 40]:
    tochki = nabor(n)
    random.shuffle(tochki)
    schet[0] = 0
    krug = minimalnyj(tochki)
    s1 = schet[0]
    schet[0] = 0
    drugoj = minimalnyj(ot_centra(tochki))
    print(n, s1, schet[0], sorted(krug) == sorted(drugoj))
print(krug)

20 103 816 True
40 203 3942 True
[[-831, 981], [486, -894]]

На двадцати точках случайный порядок потратил 103 проверки, порядок «от центра наружу» - 816; на сорока - 203 против 3942. Удвоение n удвоило первое число и увеличило второе почти в пять раз. Окружность в обоих случаях одна и та же - алгоритм верен при любом порядке, порядок влияет только на цену. Для сорока точек она опирается на две точки: это диаметр от (−831, 981) до (486, −894). Точка p внутри круга с диаметром ab, когда угол a p b не острый, - это и проверяет скалярное произведение (a − p) · (b − p) ≤ 0. Точек на входе нужно хотя бы две - minimalnyj начинает с P[0] и P[1]; для одной точки ответ - круг нулевого радиуса.

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

Разбор: обратный анализ в числах

Шесть случайных перестановок для 50 и 100 точек. Считаем проверки на точку и число дорогих шагов верхнего уровня - сравниваем со суммой 3/3 + 3/4 + … + 3/n из обратного анализа.

import random
schet = [0]

def povorot(a, b, c):
    return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])

def v_kruge(a, b, c, d):
    m = []
    for p in [a, b, c]:
        dx = p[0] - d[0]
        dy = p[1] - d[1]
        m.append([dx, dy, dx * dx + dy * dy])
    return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
            - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
            + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]))

def vnutri(krug, p):
    # круг задан двумя точками диаметра или тремя точками на окружности;
    # граница считается внутренней; всё в целых
    schet[0] = schet[0] + 1
    if len(krug) == 2:
        a, b = krug
        dx = (p[0] - a[0]) * (p[0] - b[0])
        return dx + (p[1] - a[1]) * (p[1] - b[1]) <= 0
    a, b, c = krug
    if povorot(a, b, c) < 0:
        b, c = c, b
    return v_kruge(a, b, c, p) >= 0

def s_dvumya(P, q1, q2):
    krug = [q1, q2]
    for p in P:
        if not vnutri(krug, p):
            krug = [q1, q2, p]
    return krug

def s_odnoj(P, q):
    krug = [q, P[0]]
    for j in range(1, len(P)):
        if not vnutri(krug, P[j]):
            krug = s_dvumya(P[:j], P[j], q)
    return krug

pereschetov = [0]

def minimalnyj(P):
    krug = [P[0], P[1]]
    for i in range(2, len(P)):
        if not vnutri(krug, P[i]):
            pereschetov[0] = pereschetov[0] + 1
            krug = s_odnoj(P[:i], P[i])
    return krug
random.seed(13)
for n in [50, 100]:
    tochki = []
    for i in range(n):
        x = random.randint(-1000, 1000)
        tochki.append([x, random.randint(-1000, 1000)])
    harm = 0
    for i in range(3, n + 1):
        harm = harm + 3 / i
    schet[0] = 0
    pereschetov[0] = 0
    for t in range(6):
        random.shuffle(tochki)
        minimalnyj(tochki)
    print(n, round(schet[0] / 6 / n, 2), round(pereschetov[0] / 6, 2),
          round(harm, 2))

50 5.33 7.17 9.0
100 5.76 8.33 11.06

Проверок на точку 5,3 и 5,8 - число почти не меняется при удвоении n: время линейно. Дорогих шагов в среднем 7,2 и 8,3, а оценка сверху - 9,0 и 11,1. Оценка выше среднего, потому что «не больше трёх опорных точек» - иногда их две, и тогда вероятность 2 / i. Оценка растёт как логарифм, 3 · ln n плюс константа, а среднее держится под ней.

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

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

Вход «как пришёл». Данные часто приходят упорядоченными - по времени, по координате, по обходу контура. Это как раз худшие порядки для инкрементных алгоритмов. Перемешивание стоит O(n) и снимает проблему.

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

Окружность через центр и радиус во float. Точка на границе - то внутри, то снаружи, и починки запускаются лишний раз или пропускаются. Опорные точки и определитель решают это точно.

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

Ограничивающие окружности и сферы в играх и физических движках - алгоритм Вельцля (1991) в той же схеме. Где поставить вышку или склад, чтобы покрыть всех наименьшим радиусом. Построение Делоне в CGAL и других библиотеках - вставка с перемешиванием; на практике порядок ещё и упорядочивают по кривой Гильберта внутри случайных пачек (так называемый смещённый случайный порядок, BRIO), чтобы соседние вставки трогали соседние места в памяти, - случайность сохраняет гарантию, упорядоченность - кеш.

Теперь сам

1. Точек 1000, окружность всегда опирается на три точки. Сколько в среднем дорогих шагов верхнего уровня?
2. Почему в Делоне ожидаемое число переворотов на вставку не больше примерно трёх?
3. Можно ли в алгоритме окружности вместо перемешивания один раз отсортировать точки по x?