Цепь под заказ: алгоритм Метрополиса и время смешивания

Если ты запишешь правило Метрополиса, проверишь для него равновесие между двумя состояниями и объяснишь, от чего зависит время смешивания, - листай до «Теперь сам».

Как думаешь, можно ли выбирать по весам, не зная их суммы

Нужно выбирать случайный объект так, чтобы вероятность каждого была пропорциональна его весу. Объектов больше, чем атомов во Вселенной: это, например, все раскраски графа или все расписания. Вес любого одного объекта посчитать легко.

Ответь себе до чтения: возможно ли это? Чтобы превратить вес в вероятность, его делят на сумму всех весов, а сумму по такому множеству не посчитать никогда. Возможно. Достаточно уметь сравнивать веса двух соседних объектов - и долго блуждать.

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

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

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

Механизм: равновесие между каждой парой

Дано множество состояний и вес w(x) > 0 у каждого. Цель - цепь, стационарное распределение которой π(x) = w(x) / Z, где Z - сумма всех весов, которую мы не знаем.

Алгоритм Метрополиса. У каждого состояния есть соседи; способ выбрать случайного соседа симметричен: y предлагается из x с той же вероятностью, что x из y. Если у состояний разное число соседей, «случайный сосед» несимметричен; тогда выбирают одно из d направлений, где d - наибольшее число соседей, и остаются на месте, если в выбранном направлении соседа нет.

текущее состояние x:
    предложить случайного соседа y
    с вероятностью min(1, w(y) / w(x)):  перейти в y
    иначе:                               остаться в x

В правило входит только отношение w(y) / w(x). Сумма Z сократилась бы в нём, если бы мы подставили вероятности, - поэтому её и не нужно знать.

Почему стационарное распределение то самое. Достаточно проверить условие сильнее обычного равновесия: для каждой пары состояний поток вероятности из x в y равен потоку из y в x. Это детальное равновесие: π(x) · P(x → y) = π(y) · P(y → x).

Пусть w(y) ≤ w(x), и пусть q - вероятность предложить один другому. Из x в y переход принимается с вероятностью w(y) / w(x), из y в x - всегда:

из x в y:   π(x) · q · w(y)/w(x)  =  w(x)/Z · q · w(y)/w(x)  =  q · w(y) / Z
из y в x:   π(y) · q · 1          =  q · w(y) / Z

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

Сколько ждать. Цепь сходится, но не мгновенно. Расстояние по вариации между двумя распределениями - половина суммы модулей разностей вероятностей по всем состояниям; это наибольшая разница, с которой два распределения оценивают вероятность какого-либо события. Время смешивания - число шагов, после которого расстояние от распределения цепи до стационарного становится меньше заданного порога при любом начале.

Правило

Алгоритм Метрополиса: предложить соседа симметричным способом и принять переход с вероятностью min(1, w(y) / w(x)). Полученная цепь удовлетворяет детальному равновесию π(x) · P(x → y) = π(y) · P(y → x) для π(x), пропорционального w(x), и потому имеет это распределение стационарным. Расстояние по вариации - 1/2 · сумма |p(x) − π(x)|; время смешивания - число шагов до малого расстояния. Состояние цепи после времени смешивания - почти честная выборка; соседние состояния цепи зависимы.

Разбор: пять состояний с весами 1, 2, 3, 4 и 10

Состояния стоят по кругу, сосед - следующий или предыдущий, выбирается монеткой. Сорок тысяч шагов; печатается доля времени в каждом состоянии и целевые доли.

import random

vesa = [1, 2, 3, 4, 10]
n = len(vesa)

random.seed(77)
shagov = 40000
byl = [0] * n
x = 0
for t in range(shagov):
    if random.getrandbits(1) == 1:
        y = (x + 1) % n
    else:
        y = (x - 1) % n
    if random.random() < vesa[y] / vesa[x]:
        x = y
    byl[x] = byl[x] + 1

zamer = []
cel = []
for i in range(n):
    zamer.append(round(byl[i] / shagov, 3))
    cel.append(vesa[i] / sum(vesa))
print(zamer)
print(cel)

[0.048, 0.099, 0.152, 0.204, 0.497]
[0.05, 0.1, 0.15, 0.2, 0.5]

Доли времени отличаются от целевых не больше чем на 0,004. Сумма весов в цикле не используется: sum(vesa) стоит только в строке, которая печатает ответ для сравнения. Условие random.random() < vesa[y] / vesa[x] - это и есть «принять с вероятностью min(1, отношение)»: когда отношение больше единицы, случайное число от 0 до 1 меньше него всегда.

Из состояния с весом 10 цепь уходит неохотно: сосед с весом 4 принимается с вероятностью 0,4, сосед с весом 1 - с вероятностью 0,1. Обратно в него цепь переходит всегда, когда предложено. Так и получается, что в тяжёлом состоянии проводится половина времени.

Разбор: кольцо против полного графа

Две цепи на 16 состояниях, обе с равномерным стационарным распределением и обе ленивые: с вероятностью 1/2 стоят. Первая ходит по кольцу - к левому или правому соседу. Вторая - по полному графу: к любому другому состоянию. Обе начинают из состояния 0. Здесь нет случайных чисел: считается точное распределение после t шагов и его расстояние по вариации до равномерного.

def rasstoyanie(r, n):
    s = 0
    for x in r:
        s = s + abs(x - 1 / n)
    return s / 2

def shag_po_kolcu(r, n):
    novoe = []
    for v in range(n):
        novoe.append(r[v] / 2 + r[(v - 1) % n] / 4 + r[(v + 1) % n] / 4)
    return novoe

def shag_po_polnomu(r, n):
    novoe = []
    for v in range(n):
        novoe.append(r[v] / 2 + (1 - r[v]) / 2 / (n - 1))
    return novoe

n = 16
kolco = [0] * n
polnyj = [0] * n
kolco[0] = 1
polnyj[0] = 1
t = 0
for otmetka in [1, 4, 16, 64, 256]:
    while t < otmetka:
        kolco = shag_po_kolcu(kolco, n)
        polnyj = shag_po_polnomu(polnyj, n)
        t = t + 1
    print(t, round(rasstoyanie(kolco, n), 4), round(rasstoyanie(polnyj, n), 4))

1 0.8125 0.4375
4 0.6172 0.0445
16 0.3472 0.0
64 0.0525 0.0
256 0.0 0.0

Столбцы: число шагов, расстояние для кольца, расстояние для полного графа. Полный граф смешался за четыре шага: расстояние 0,04. Кольцу после шестнадцати шагов до равновесия ещё далеко - 0,35, и только к шестьдесят четвёртому шагу расстояние падает до 0,05. Это диффузия из урока про блуждания: чтобы вероятность растеклась по кольцу длины n, нужно порядка n² шагов.

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

Тот же вопрос в житейской форме - сколько раз тасовать колоду. Тасование - это шаг цепи на множестве перестановок. Если каждый раз перекладывать верхнюю карту в случайное место, для 52 карт нужно около 52 · ln 52 ≈ 205 перекладываний. Для обычного тасования врезкой известен результат: около семи раз.

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

Несимметричное предложение без поправки. Если соседа y из x предлагают чаще, чем x из y, правило с одним отношением весов даёт другое стационарное распределение. В вероятность принятия нужно добавить отношение вероятностей предложений.

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

Начальный участок включают в выборку. Первые шаги помнят, откуда цепь стартовала. Их отбрасывают. Сколько именно, зависит от времени смешивания, которое редко известно точно; на практике запускают несколько цепей из далёких начал и отбрасывают участок, пока их оценки не сойдутся.

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

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

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

Оптимизация. Если взять вес, быстро убывающий с ростом «стоимости» решения, цепь будет проводить почти всё время в дешёвых решениях. Постепенно делая вес всё более резким, получают метод имитации отжигаРазбирается в разделе 22, глава «Локальный поиск, отжиг, эвристики»: Отжиг: иногда соглашаться на худшее: поиск, который соглашается на ухудшения сначала часто, потом всё реже и поэтому не застревает в первом же локальном минимуме.

Приближённый счёт. Число раскрасок графа или число паросочетаний точно посчитать трудно, но если уметь выбирать случайный объект почти равномерно, количество оценивают цепочкой долей: какая часть раскрасок остаётся раскрасками, если добавить в граф ещё одно ребро, - и перемножают эти доли, начав с графа без рёбер, где ответ известен. Так выборка из огромного множества превращается в оценку его размера.

Теперь сам

1. Веса w(x) = 4, w(y) = 1. С какой вероятностью принимается переход из x в y и из y в x?
2. Два состояния с весами 2 и 6, каждое предлагает другое с вероятностью 1/2. Проверь детальное равновесие для π = (0,25, 0,75).
3. Найди расстояние по вариации между распределениями (0,5, 0,5, 0) и (1/3, 1/3, 1/3).
4. По таблице из разбора: во сколько раз кольцо из 16 состояний смешивается медленнее полного графа, если мерить числом шагов до расстояния около 0,05?
5. Цепь Метрополиса на множестве, состоящем из двух плотных областей с единственным переходом между ними, запустили на миллион шагов. Доли времени стабилизировались. Можно ли верить выборке?