Алгоритм Ванга — Ландау

Алгоритм Ванга-Ландау, предложенный Фугао Вангом и Дэвидом Ландау, это метод Монте-Карло, предназначенный для расчета плотности состояний системы. Метод выполняет немарковские случайные переходы для построения плотности состояний, посещая все возможные состояния. Алгоритм Ванга и Ландау, это важный для получения плотности состояний метод, требуемый для выполнения мультиканонического моделирования

Алгоритм Ванга-Ландау может быть применен к любой системе, которая характеризуется некоторым параметром (например, энергией, объемом и др.). К примеру, он может быть использован для численного интегрирования и моделирования белков.

ОписаниеПравить

Алгоритм Ванга-Ландау является реализацией метода энтропического моделирования, в котором изучается плотность состояний с помощью блуждания в пространстве энергий с равновероятным посещением всех энергетических состояний. Алгоритм решает проблему подбора подходящих вероятностей перехода для получения требуемого при энтропическом моделировании равномерного посещения энергетических состояний и, следовательно, позволяет получить плотность состояний $\mathrm{\Omega <!--\; \Omega \; -->}(E)$$\Omega (E)$ .

АлгоритмПравить

Рассмотрим систему в фазовом пространстве $\mathrm{\Omega <!--\; \Omega \; -->}$$\Omega$  и энергию $E$$E$ , изменение которой ограничено диапазоном $E_{min}\le <!--\; \le \; -->E\le <!--\; \le \; -->E_{max}$$E_{min}\leq E\leq E_{max}$ . Пусть рассматриваемая система имеет плотность вероятности $\mathrm{\rho <!--\; \rho \; -->}(E)\sim <!--\; \sim \; -->\exp <!--\; \; -->(S(E))$$\rho (E)\sim \exp(S(E))$ , которую нам требуется посчитать. Поскольку алгоритм Ванга и Ландау работает с дискретным спектром энергии, диапазон $E_{min}\le <!--\; \le \; -->E\le <!--\; \le \; -->E_{max}$$E_{min}\leq E\leq E_{max}$  разбивается на конечное число ($M$$M$ ) равных отрезков («ящиков»), размер которых равен $\mathrm{\Delta <!--\; \Delta \; -->}$$\Delta$ . Таким образом:

$M=\frac{E_{max}-<!--\; -\; -->E_{min}}{\mathrm{\Delta <!--\; \Delta \; -->}}$$M={\frac {E_{max}-E_{min}}{\Delta }}$ .

С учетом этого дискретного спектра, алгоритм имеет следующие начальные условия:

• $S(E_i)=0,i=1,2,...,M$$S(E_{i})=0,i=1,2,...,M$ 
• $f=1$$f=1$  (используется далее, как добавка к энтропии $S(E_i)$$S(E_{i})$  после каждого принятого шага)
• Выбирается случайная конфигурация системы $\stackrel{\to <!--\; \to \; -->}{r}\in <!--\; \in \; -->\mathrm{\Omega <!--\; \Omega \; -->}$${\vec {r}}\in \Omega$ 

Алгоритм выполняет моделирование в мультиканоническом ансамбле: случайное блуждание Метрополиса-Гастингса по фазовому пространству $\mathrm{\Omega <!--\; \Omega \; -->}$$\Omega$  системы с распределением вероятности $P(\stackrel{\to <!--\; \to \; -->}{r})=1/\mathrm{\rho <!--\; \rho \; -->}(E(\stackrel{\to <!--\; \to \; -->}{r}))=\exp <!--\; \; -->(-<!--\; -\; -->S(E(\stackrel{\to <!--\; \to \; -->}{r})))$$P({\vec {r}})=1/\rho (E({\vec {r}}))=\exp(-S(E({\vec {r}})))$  и вероятностью генерации нового состояния, данной распределением вероятности $g(\stackrel{\to <!--\; \to \; -->}{r}\to <!--\; \to \; -->{\stackrel{\to <!--\; \to \; -->}{r}}^{\prime })$$g({\vec {r}}\rightarrow {\vec {r}}')$ , которое выбирается произвольно (обычно любое состояние может быть сгенерировано с равной вероятностью). В процессе моделирования, посещение каждого «ящика» записывается в гистограмму $H(E)$$H(E)$  (то есть, значение $H(E)$$H(E)$  увеличивается на единицу). Как и в алгоритме Метрополиса-Гастингса, генерация и принятие нового состояния выполняется следующим образом:

1. генерация нового состояния ${\stackrel{\to <!--\; \to \; -->}{r}}^{\prime }\in <!--\; \in \; -->\mathrm{\Omega <!--\; \Omega \; -->}$${\vec {r}}'\in \Omega$  согласно распределению вероятности $g(\stackrel{\to <!--\; \to \; -->}{r}\to <!--\; \to \; -->{\stackrel{\to <!--\; \to \; -->}{r}}^{\prime })$$g({\vec {r}}\rightarrow {\vec {r}}')$ 
2. принятие/отклонение нового состояния, производится следующим образом:

Если энтропия нового состояния меньше текущего, то оно сразу принимается. Если же энтропия увеличилась, то новое состояние принимается с вероятностью: $A=e^{S-<!--\; -\; -->S^{\prime }}\frac{g({\stackrel{\to <!--\; \to \; -->}{r}}^{\prime }\to <!--\; \to \; -->\stackrel{\to <!--\; \to \; -->}{r})}{g(\stackrel{\to <!--\; \to \; -->}{r}\to <!--\; \to \; -->{\stackrel{\to <!--\; \to \; -->}{r}}^{\prime })}$$A=e^{S-S'}{\frac {g({\vec {r}}'\rightarrow {\vec {r}})}{g({\vec {r}}\rightarrow {\vec {r}}')}}$ , где $S=S(E(\stackrel{\to <!--\; \to \; -->}{r}))$$S=S(E({\vec {r}}))$  и $S^{\prime }=S(E({\stackrel{\to <!--\; \to \; -->}{r}}^{\prime }))$$S'=S(E({\vec {r}}'))$ .

То есть, общая формула выглядит следующим образом:

$A(\stackrel{\to <!--\; \to \; -->}{r}\to <!--\; \to \; -->{\stackrel{\to <!--\; \to \; -->}{r}}^{\prime })=min\left(1,e^{S-<!--\; -\; -->S^{\prime }}\frac{g({\stackrel{\to <!--\; \to \; -->}{r}}^{\prime }\to <!--\; \to \; -->\stackrel{\to <!--\; \to \; -->}{r})}{g(\stackrel{\to <!--\; \to \; -->}{r}\to <!--\; \to \; -->{\stackrel{\to <!--\; \to \; -->}{r}}^{\prime })}\right)$$A({\vec {r}}\rightarrow {\vec {r}}')=\min \left(1,e^{S-S'}{\frac {g({\vec {r}}'\rightarrow {\vec {r}})}{g({\vec {r}}\rightarrow {\vec {r}}')}}\right)$ .

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

После каждого шага генерации-принятия система переходит в некоторое состояние $E_i$$E_{i}$ , значение $H(E_i)$$H(E_{i})$  увеличивается на единицу, а также выполняется следующее изменение:

$S(E_i)\leftarrow <!--\; \leftarrow \; -->S(E_i)+f$$S(E_{i})\leftarrow S(E_{i})+f$ 

Это важный шаг алгоритма, и это то, что делает алгоритм Ванга-Ландау немарковским: случайный процесс теперь зависит от истории процесса. Таким образом, когда в следующий раз будет предложено состояние с энергией $E_i$$E_{i}$ , это состояние будет отклонено с большей вероятностью; в этом смысле, алгоритм принуждает систему посещать все состояния с одинаковой частотой. Как следствие, гистограмма $H(E)$$H(E)$  становится все более и более плоской. Хотя, эта равномерность зависит от того, насколько посчитанная энтропия близка к точной энтропии, что зависит от $f$$f$ . Для улучшения приближения точной энтропии (и, таким образом, равномерности гистограммы), $f$$f$  уменьшается после $M$$M$  шагов генерации-принятия:

$f\leftarrow <!--\; \leftarrow \; -->f/2$$f\leftarrow f/2$ 

Через некоторое время было показано, что при изменении $f$$f$ постоянным делением на два алгоритм может не сходиться. Небольшая модификация метода Ванга-Ландау позволяет избежать этого: производится деление не на два, а на $t$$t$ , при чем $t$$t$  пропорционально шагу моделирования.

В результате использования этого алгоритма происходит автоматическая настройка весов вероятности перехода, которые одновременно определяют плотности состояний. По окончании расчета вычисляется массив $\mathrm{\Omega <!--\; \Omega \; -->}(E)=\exp <!--\; \; -->S(E)$$\Omega (E)=\exp S(E)$ и нормируется на единицу.

Пример кодаПравить

Ниже показан пример кода на Python, в котором предполагается симметричность функции распределения $g$$g$ :

$g({\mathit{x}}^{\prime }\to <!--\; \to \; -->\mathit{x})=g(\mathit{x}\to <!--\; \to \; -->{\mathit{x}}^{\prime })$$g({\boldsymbol {x}}'\rightarrow {\boldsymbol {x}})=g({\boldsymbol {x}}\rightarrow {\boldsymbol {x}}')$ 

currentEnergy = system.randomConfiguration() # случайная генерация начального состояния системы
while (f > epsilon):
    system.proposeConfiguration() # генерация новой конфигурации
    proposedEnergy = system.proposedEnergy() # вычисление энергии нового состояния

    if (random() < exp(entropy[currentEnergy]-entropy[proposedEnergy])):
        # если принято, обновляем энергию и систему
        currentEnergy = proposedEnergy
        system.acceptProposedConfiguration()
    else:
        # если отклонено
        system.rejectProposedConfiguration()
    
    H[currentEnergy] += 1
    entropy[currentEnergy] += f
    
    if (isFlat(H)): # isFlat проверяет достаточно ли гладкая гистограмма (например, 95%)
        H[:] = 0
        f *= 0.5 # refine the f parameter