RND Случайность Приём
Лас-Вегас и Монте-Карло
Также: Las Vegas and Monte Carlo algorithms · алгоритм Лас-Вегаса · алгоритм Монте-Карло · проверка Фрейвалдса · Freivalds' algorithm · лемма Шварца — Зиппеля · Schwartz–Zippel lemma · перезапуски · ZPP, RP, BPP
Случайный алгоритм может платить за скорость двумя способами. Лас-Вегас всегда отвечает верно, но время его работы случайно; Монте-Карло работает быстро и предсказуемо, но иногда ошибается — с вероятностью, которую мы выбираем. Один превращается в другой остановкой по таймеру или проверкой ответа, а ошибку можно сделать сколь угодно малой повторами.
Сервер прислал размытую фотографию: каждую точку он заменил суммой девяти соседей по строке — умножил матрицу снимка на матрицу размытия. Честно ли он посчитал? Пересчитать заново — 16,8 миллиона умножений. Есть способ проверить за двести тысяч, но он может ошибиться: сказать «всё верно», когда ошибка есть. Ошибается он не чаще чем в половине случаев, а двадцать проверок подряд пропустят ошибку с вероятностью меньше одной миллионной. Испортите точку на снимке справа и проверьте сами.
Одна испорченная точка ловится примерно в половине проверок, и это не случайное совпадение, а точная вероятность: проверка видит ошибку в столбце j, только если случайный вектор в этом столбце равен единице. Нажмите «хитрая порча»: две ошибки в одной строке с разными знаками, которые взаимно гасятся в сумме строки. Обычная контрольная сумма их не заметит, а проверка по-прежнему ловит их в половине случаев.
Идея
Случайность в алгоритме можно потратить двумя способами. Алгоритм Лас-Вегаса всегда отвечает верно, а случайно только время его работы. Алгоритм Монте-Карло работает предсказуемое время, но может ошибиться — с вероятностью, которую мы контролируем. Эти два вида переходят друг в друга: Лас-Вегас, остановленный по таймеру, становится Монте-Карло, а Монте-Карло, чей ответ можно быстро проверить, становится Лас-Вегасом. Ошибку Монте-Карло повторы делают сколь угодно малой, а долгие хвосты Лас-Вегаса срезают перезапуски.
Названия придуманы в пару: «Монте-Карло» появилось в 1940-х у физиков, моделировавших нейтроны (об этом — статья о методе Монте-Карло), а «Лас-Вегас» в 1979 году предложил Ласло Бабаи — одно казино в пару к другому. Основы — случайная быстрая сортировка как Лас-Вегас и стягивание рёбер как Монте-Карло — разобраны в главе курса «Подбросим монетку»; здесь речь о том, что с этими алгоритмами можно делать.
У ошибки Монте-Карло бывает направление. Проверка Фрейвалдса ошибается в одну сторону: если она сказала «ошибка есть», ошибка точно есть, а «всё верно» бывает ложным. Так же устроены тест простоты Миллера — Рабина (его «составное» всегда правда) и фильтр Блума (его «нет» всегда правда). Бывают и алгоритмы с двусторонней ошибкой — они могут ошибиться в любую сторону, и усиливать их приходится иначе.
Проверка Фрейвалдса
Даны матрицы A, B и C размера n × n; верно ли, что AB = C? Перемножить A и B по определению стоит n³ умножений. Русиньш Фрейвалдс в 1977 году предложил взять случайный вектор r из нулей и единиц и сравнить A(Br) с Cr. Умножение матрицы на вектор стоит n² операций, и таких умножений три.
Если AB = C, проверка всегда отвечает «верно». Если AB ≠ C, а каждая координата r независимо равна 0 или 1 с вероятностью ½, то проверка ошибочно отвечает «верно» с вероятностью не больше ½.
Первая часть очевидна: A(Br) = (AB)r = Cr. Для второй обозначим D = AB − C; проверка ошибается ровно тогда, когда Dr = 0. Раз D ≠ 0, в ней есть ненулевой элемент dij. Посмотрим на i-ю координату:
$$(Dr)_i = d_{ij}\,r_j + \sum_{k \ne j} d_{ik}\,r_k.$$Представим, что все rk с k ≠ j уже выбраны, а rj — ещё нет; это законно, потому что координаты независимы. Сумма справа стала каким-то числом s, и (Dr)i = dijrj + s. Из двух значений rj ∈ {0, 1} обнулить это выражение может самое большее одно: если бы оба, то s = 0 и dij + s = 0, откуда dij = 0. Значит, (Dr)i = 0 с вероятностью не больше ½ при любых остальных координатах, а Dr = 0 — тем более.
Из доказательства видно и то, что показывает сцена: одна ошибка в клетке (i, j) ловится ровно тогда, когда rj = 1, то есть ровно в половине случаев, и проявляется в строке i. В листинге ниже матрица A — граф ссылок трёх учебников сайта (202 страницы, 1624 ссылки), C = A·A — число путей из двух щелчков между страницами. Верное произведение проходит все 20 проверок; в произведении, где одна клетка увеличена на единицу, 64 отдельные проверки ловят ошибку 34 раза.
Если алгоритм с односторонней ошибкой ошибается с вероятностью не больше q, то k независимых запусков, где ответ «ошибка есть» засчитывается, если его дал хоть один, ошибаются с вероятностью не больше qk.
Ложный ответ «всё верно» возможен, только если ошибку пропустили все k запусков. Запуски независимы, поэтому вероятность этого — произведение k вероятностей, каждая не больше q. А при AB = C ни один запуск не скажет «ошибка есть», так что в эту сторону ошибиться нельзя вовсе.
Двадцать проверок Фрейвалдса стоят 60n² операций и ошибаются с вероятностью меньше 2−20 ≈ 10−6. При n = 256, как на фотографии, это четыре миллиона операций против 16,8 миллиона на полный пересчёт, а при n = 10 000 — 6·109 против 1012.
Лас-Вегас: верно, но неизвестно когда
Теперь обратный случай. Решим судоку из банка задач простым локальным поиском: заполним каждый квадрат 3×3 недостающими цифрами наугад и будем менять местами две свободные клетки одного квадрата, если от этого повторов в строках и столбцах не становится больше. Когда повторов не осталось, судоку решено, и ответ проверен — он верен всегда. Но сколько шагов это займёт?
Гистограмма выглядит необычно. Многие прогоны решают задачу за несколько тысяч шагов, но у первого судоку примерно половина не решила её и за сто тысяч: такой прогон застрял — в нём осталось два-три повтора, которые нельзя убрать ни одним обменом, не создав новых. Чтобы выбраться, спуску нужна редкая удача, и время его работы описывается распределением с «тяжёлым хвостом»: среднее время огромно, хотя медианное невелико. Карла Гомес, Барт Селман, Нуно Крато и Генри Каутц в 2000 году показали, что так ведут себя многие переборы для задач выполнимости и расписаний.
Что делать с таким алгоритмом? Его можно превратить в Монте-Карло: остановить по таймеру и сказать «не знаю». А можно перезапустить: оборвать прогон, который затянулся, и начать заново со случайного заполнения. В листинге ниже из восьми прогонов без перезапусков четыре не решили задачу за 20 000 шагов, а остальные справились за 1303–4031 шаг. С перезапуском каждые 2000 шагов задача решается в среднем за 5001 шаг.
Пусть время работы алгоритма Лас-Вегаса — случайная величина T ≥ 0 со средним E[T]. Если остановить его через t шагов и ответить «не знаю», это случится с вероятностью P(T > t) ≤ E[T] / t. В частности, при t = 2E[T] получаем алгоритм Монте-Карло с вероятностью отказа не больше ½.
Величина T не меньше величины, которая равна t при T > t и нулю иначе. Средние сохраняют неравенство: E[T] ≥ t·P(T > t). Остаётся разделить на t.
На судоку неравенство Маркова бесполезно: половина прогонов не уложилась даже в сто тысяч шагов, и среднее время там больше пятидесяти тысяч — оценка для t = 2000 получается больше 25, а вероятность не бывает больше единицы. Но настоящее распределение можно измерить, и вот что из него следует.
Пусть прогоны независимы и одинаково распределены, и каждый обрывается через t шагов, если не успел, после чего начинается новый. Если P(T ≤ t) > 0, то среднее число шагов до ответа равно
$$\frac{\mathbb E\,[\min(T,\ t)]}{P(T \le t)}.$$Обозначим p = P(T ≤ t) и через X — искомое среднее. Первый прогон стоит min(T, t) шагов. С вероятностью p он закончился ответом; с вероятностью 1 − p его оборвали, и всё начинается заново, а прошлое ни на что не влияет. Поэтому X = E[min(T, t)] + (1 − p)·X, откуда p·X = E[min(T, t)]. Чтобы законно вычесть (1 − p)X из обеих частей, X должно быть конечно; это так, потому что число прогонов распределено геометрически со средним 1/p, а каждый стоит не больше t шагов, так что X ≤ t/p.
Эта формула — кривая под гистограммой в сцене, а её минимум — лучшая граница. Замер ниже: 400 прогонов на первом судоку, и по ним цена решения с перезапуском каждые t шагов. При t = 2000 решение стоит в среднем 5635 шагов; при t = 500 — 33 257, потому что за 500 шагов прогон почти никогда не успевает; при t = 100 000 — 103 101, потому что застрявшие прогоны тянутся до конца.
Лучшая граница зависит от задачи, а распределение заранее неизвестно. Майкл Лаби, Алистер Синклер и Дэвид Цукерман в 1993 году нашли универсальное правило: перезапускать через 1, 1, 2, 1, 1, 2, 4, 1, 1, 2, 1, 1, 2, 4, 8, … единиц времени. Такая стратегия в среднем работает не дольше лучшей фиксированной границы, умноженной на логарифмический множитель, и они доказали, что этот множитель неизбежен: стратегия, не знающая распределения, на каком-то распределении проигрывает лучшей границе в логарифм раз. По этой последовательности перезапускает поиск, например, решатель выполнимости MiniSat.
Обратный переход — от Монте-Карло к Лас-Вегасу — возможен, когда ответ можно проверить. Если одна попытка стоит a шагов и даёт верный ответ с вероятностью p, а проверка стоит v, то «пробовать и проверять, пока не получится» — это Лас-Вегас со средним временем (a + v)/p: число попыток распределено геометрически со средним 1/p, как в доказательстве выше. Наш решатель судоку так и устроен: оборванный прогон — неудачная попытка, а «ноль повторов» — проверка.
Двусторонняя ошибка и голосование
Если алгоритм может ошибиться в обе стороны, «хоть один сказал „нет“» уже не годится: в любую сторону бывают ложные ответы. Тогда запускают его k раз и берут ответ большинства. Если каждый запуск верен с вероятностью больше ½, большинство ошибается редко — и экспоненциально реже с ростом k.
Пусть каждый из k независимых запусков ошибается с вероятностью p < ½. Тогда большинство ошибается — то есть неверных ответов не меньше половины — с вероятностью не больше (4p(1 − p))k/2.
Вероятность, что ошиблись ровно j запусков, равна C(k, j)·pj(1 − p)k−j. Перепишем множитель: pj(1 − p)k−j = (p(1 − p))k/2·(p/(1 − p))j − k/2. При j ≥ k/2 показатель последнего множителя неотрицателен, а его основание p/(1 − p) меньше единицы, так что он не больше 1. Поэтому сумма по всем j ≥ k/2 не больше (p(1 − p))k/2, умноженного на сумму биномиальных коэффициентов, а та не больше 2k. Получаем 2k(p(1 − p))k/2 = (4p(1 − p))k/2.
При p = ⅓ основание 4·⅓·⅔ = 8/9, и сто запусков ошибаются с вероятностью меньше (8/9)50 ≈ 0,003, а тысяча — меньше 10−25. Оценка грубая (точнее её даёт неравенство Чернова), но главное видно и из неё: каждый новый запуск умножает вероятность ошибки на постоянную меньше единицы. Поэтому для алгоритмов Монте-Карло не важно, ошибаются они в трети случаев или в одной сотой: несколько десятков повторов делают ошибку меньше вероятности аппаратного сбоя.
Отпечатки: многочлен в случайной точке
Самый мощный источник алгоритмов Монте-Карло — простое свойство многочленов: ненулевой многочлен не может обращаться в ноль слишком часто.
Пусть f(x1, …, xm) — ненулевой многочлен степени d с коэффициентами из поля, а S — конечное множество элементов этого поля. Если каждое xi независимо выбрать равновероятно из S, то f(x1, …, xm) = 0 с вероятностью не больше d / |S|.
Индукция по числу переменных. При m = 1 у ненулевого многочлена степени d не больше d корней: если r — корень, то f(x) делится на x − r без остатка (остаток от деления на x − r равен f(r) = 0), частное имеет степень d − 1 и остальные корни f — его корни, потому что в поле произведение равно нулю, только если ноль один из множителей. Значит, из |S| равновероятных значений не больше d обнуляют f.
Шаг. Выделим наибольшую степень k переменной x1: f = x1k·g(x2, …, xm) + (слагаемые, где степень x1 меньше k), где g ≠ 0 и степень g не больше d − k. Сначала выбираем x2, …, xm. По предположению индукции g обнулится с вероятностью не больше (d − k)/|S|. Если же g ≠ 0, то при выбранных остальных переменных f — ненулевой многочлен от одной x1 степени k, и он обнулится с вероятностью не больше k/|S|. Вероятность обнуления f не больше суммы двух этих вероятностей: d/|S|.
Применение, ради которого не нужен даже многомерный случай, — сравнение больших наборов. Два журнала веб-сервера содержат одни и те же запросы, возможно в разном порядке? Запишем каждый журнал как многочлен (z − x1)(z − x2)…(z − xn), где xi — номера запросов. От порядка он не зависит, а у разных наборов многочлены разные — их разность ненулевая и степени не больше n. Посчитаем оба многочлена в одной случайной точке z по модулю простого P = 261 − 1 — это и есть отпечаток журнала, одно число вместо сотни тысяч, и считается он за один проход без сортировки. По лемме отпечатки разных журналов совпадут с вероятностью не больше n / |S|.
В листинге — журнал сервера Космического центра Кеннеди за 13 июля 1995 года, 128 932 запроса. Отсортированная копия даёт тот же отпечаток. Копия, где один запрос потерялся, а другой записан дважды, — другой. И копия, где два номера сдвинуты навстречу друг другу на единицу, — тоже другой, хотя простая сумма номеров у неё та же. Точка z в листинге берётся из 253 значений, поэтому вероятность ложного совпадения не больше 128 932 / 253 ≈ 1,4·10−11. Та же лемма лежит в основе многих других проверок: равенства многочленов, заданных формулами, существования совершенного паросочетания (через определитель матрицы Татта) и полиномиального хеширования строк.
Что из чего получается
Теорию сложности эти превращения интересуют в общем виде. Класс ZPP — задачи, решаемые алгоритмом Лас-Вегаса за полиномиальное в среднем время. RP — алгоритмом Монте-Карло за полиномиальное время с односторонней ошибкой: при ответе «нет» он всегда прав, а при «да» ошибается с вероятностью не больше ½. coRP — то же с ошибкой в другую сторону. Классы ввёл Джон Гилл в 1970-х.
ZPP = RP ∩ coRP: задача решается Лас-Вегасом за полиномиальное в среднем время тогда и только тогда, когда для неё есть алгоритмы Монте-Карло с односторонней ошибкой в обе стороны.
Из ZPP в RP. Пусть Лас-Вегас работает в среднем t(n) шагов. Остановим его через 2t(n) шагов; если не успел, ответим «нет». По неравенству Маркова он не успевает с вероятностью не больше ½, так что на верный ответ «да» мы ошибочно скажем «нет» с вероятностью не больше ½, а «да» говорим только тогда, когда его сказал Лас-Вегас, — значит, верно. Это RP. Ответив при отказе «да», получим coRP.
Из RP ∩ coRP в ZPP. Есть алгоритм R, чьё «да» всегда верно, и алгоритм R′, чьё «нет» всегда верно. Запускаем их по очереди, пока R не скажет «да» или R′ не скажет «нет», и выдаём этот ответ — он заведомо верный. Если правильный ответ «да», то R говорит «да» с вероятностью не меньше ½ в каждом раунде; если «нет» — так же R′. Значит, число раундов распределено не хуже геометрического с вероятностью ½, и в среднем их не больше двух, а каждый полиномиален.
С двусторонней ошибкой получается класс BPP, и голосование большинством делает его ошибку сколь угодно малой. Равен ли BPP классу P, то есть можно ли любую случайность заменить детерминированным вычислением без большой потери во времени, — открытый вопрос; большинство специалистов считают, что можно.
Код
Проверка Фрейвалдса на графе ссылок сайта. A хранится списками ссылок, поэтому A(Ar) стоит два прохода по ссылкам, а заявленное произведение C — обычная плотная таблица n × n, как если бы его прислали со стороны.
Отпечатки журналов. Произведение двух 61-битных чисел в Python — обычное целое, а в Go оно не помещается в 64 бита, и его считают функцией bits.Mul64, которая отдаёт старшую и младшую половины, а остаток берут через bits.Rem64.
И Лас-Вегас с перезапусками на судоку из банка задач. Функция descent — одна попытка с потолком, with_restarts — перезапуски, пока попытка не удастся; ответ проверяется сам собой: ноль повторов.
Сколько это стоит
| Что | Время | Ошибка |
|---|---|---|
| проверить AB = C пересчётом | O(n³) (O(n2,37…) быстрыми алгоритмами) | нет |
| Фрейвалдс, k раз | O(kn²) | не больше 2−k, только «верно» при ошибке |
| отпечатки двух наборов по n элементов | O(n) умножений по модулю | не больше n / |S| |
| сравнение сортировкой | O(n log n) и память на копии | нет |
| Лас-Вегас, остановленный на c·E[T] | c·E[T] | отказ не чаще 1/c |
| Лас-Вегас с перезапуском каждые t | E[min(T, t)] / P(T ≤ t) в среднем | нет |
| Монте-Карло с проверкой, вероятность успеха p | (попытка + проверка) / p в среднем | нет |
| голосование из k запусков | k запусков | не больше (4p(1 − p))k/2 |
Общее правило читается по таблице: за повтор платят линейно, а ошибка падает экспоненциально. Поэтому алгоритм Монте-Карло с вероятностью ошибки ½ практически так же хорош, как с вероятностью 10−6: разница — двадцать повторов.
Где работает
- Проверка чужих вычислений. Мануэль Блюм и Сампат Каннан в 1989 году предложили «проверяющие программы»: вместо того чтобы доказывать правильность программы, проверять каждый её ответ быстрым случайным тестом. Проверка Фрейвалдса — классический пример.
- Сравнение и синхронизация данных. Отпечатки по модулю простого числа позволяют сравнивать файлы и наборы, не пересылая их; на той же идее держатся полиномиальные хеши строк и алгоритм Рабина — Карпа.
- Перезапуски в решателях. Решатели выполнимости и программы перебора для расписаний перезапускают поиск, чтобы срезать тяжёлые хвосты; MiniSat по умолчанию делает это по последовательности Лаби.
- Классика каталога: быстрая сортировка и быстрый выбор со случайным опорным — Лас-Вегас; тест Миллера — Рабина, минимальный разрез Каргера и фильтр Блума — Монте-Карло с односторонней ошибкой.
История
Случайные алгоритмы для «чистых» задач, а не для моделирования, появились во второй половине 1970-х. В 1976–1977 годах Майкл Рабин, Роберт Соловей и Фолькер Штрассен предложили вероятностные тесты простоты, в 1977 году Русиньш Фрейвалдс — проверку умножения матриц, в 1978–1980 годах Ричард Демилло и Ричард Липтон, Ричард Зиппель и Джейкоб Шварц независимо доказали лемму о корнях многочлена, которая теперь носит имена двух последних. Джон Гилл тогда же описал классы задач, решаемых случайными алгоритмами, а в 1979 году Ласло Бабаи, работая над проверкой изоморфизма графов, назвал алгоритмы, которые не ошибаются, но могут долго работать, алгоритмами Лас-Вегаса.
Задачи
В первой задаче нужно проверить произведение матриц, которое слишком велико, чтобы его пересчитывать. Во второй — по выборке времён работы алгоритма Лас-Вегаса выбрать лучшую границу для перезапусков; один из тестов — настоящие прогоны решателя судоку из статьи.
Проверка за сервером Когда начинать зановоДанные: Фотографии автора сайта Граф ссылок учебников legost.in: базовый курс, «Царица наук», каталог алгоритмов Журнал веб-сервера Космического центра Кеннеди, 13 июля 1995: NASA-HTTP, Internet Traffic Archive, свободно распространяется Головоломки: банк Sudoku Exchange, общественное достояние; оценки трудности — Sukaku Explainer