Интегрирование методом Монте-Карло

Общий подход

Давайте теперь рассмотрим, как можно обобщить этот метод. Пусть нам необходимо вычислить какой-то многомерный интеграл по какой-то области \(\Omega \subset \mathbb{R}^D\) следующего вида:

\[ I=\int_{\boldsymbol{r} \in \Omega} d^D \boldsymbol{r} f(\boldsymbol{r}), \quad \boldsymbol{r}=\left(x_1, \ldots, x_D\right) \]

Сразу напрашивается тривиальное обобщение: давайте просто случайным образом выбирать точку из области \(\boldsymbol{r}_i \in \Omega\) (предположим, у нас имеется такой алгоритм); тогда мы могли бы записать:

\[ I=V_{\Omega} \int_{\boldsymbol{r} \in \Omega} \frac{d^D \boldsymbol{r}}{V_{\Omega}} f(\boldsymbol{r})=V_{\Omega}\langle f(\boldsymbol{r})\rangle_{\Omega} \underset{N \gg 1}{\approx} V_{\Omega} \times \frac{\sum_{i=1}^N f\left(\boldsymbol{r}_i\right)}{N} . \]

(где \(V_{\Omega}-\) объём области \(\Omega\), который, как мы предполагаем, мы тоже умеем вычислять).

По сути, это обобщение той самой предыдущей формулы для вычисления числа \(\pi\) при \(V_{\Omega}-\) квадрат, \(f(\boldsymbol{r})\) - индикаторная функция единичного круга.

Какие в этом метода есть подводные камни?

Проблема некомпактного носителя

Первый подводный камень заключается в том, что будет если функция \(f(\boldsymbol{r})\) отлична от нуля во всём пространстве? Формально , это соответствует \(\Omega=\mathbb{R}^D\) и \(V_{\Omega}=\infty\).

Первый «наивный» подход - обрезка функции.

Вместо области \(\Omega\) посчитаем интеграл на некоторой меньшую компактную область \(\Omega^{\prime} \subset \Omega\), где функция существенно отлична от нуля. Некомпактный остаток на \(\Omega \backslash \Omega^{\prime}\) оценим аналитически.

У этого подхода есть много проблем:

• Сложность аналитической оценки остатка

• Увеличение области \(\Omega\) повышает погрешность метода Монте-Карло, т.к. \(\sigma_I^2 \propto V_{\Omega}^2\)

• В случае непрямоугольной \(\Omega^{\prime}\) (например, какая-нибудь кривая), то выбирать случайным образом из неё точки тоже является сложной задачей

Второй «крутой» подход - выбор функции распределения.

Вообще говоря, нам необязательно выбирать точку из равномерного распределения. Давайте модифицируем наш интеграл следуюшим образом:

\[ I=\int d^D \boldsymbol{r} f(\boldsymbol{r}) \equiv \int P(\boldsymbol{r}) d^D \boldsymbol{r} \cdot \frac{f(\boldsymbol{r})}{P(\boldsymbol{r})} \]

Если функция \(P(\boldsymbol{r}) \geq 0\) и \(\int d^D \boldsymbol{r} P(\boldsymbol{r})=1\), то эту формулу можно интерпретировать следующим образом:

\[ I \equiv\left\langle\frac{f(\boldsymbol{r})}{P(\boldsymbol{r})}\right\rangle \approx I_N \equiv \frac{1}{N} \sum_{i=1}^N \frac{f\left(\boldsymbol{r}_i\right)}{P\left(\boldsymbol{r}_i\right)} \]

где \(\boldsymbol{r}_i\) выбираются из известного распределения \(P(\boldsymbol{r})\).

Наша задача заключается в выборе из этого набора распределения \(P(\boldsymbol{r})\), который максимально похож на функцию \(f(\boldsymbol{r})\) - тем самым мы минимизируем дисперсию \(\left\langle\left\langle(f(\boldsymbol{r}) / P(\boldsymbol{r}))^2\right\rangle\right\rangle\).

Возникает вопрос, а как нам генерировать \(\boldsymbol{r}_i\) из определенного распределения \(P(\boldsymbol{r})\)? Конечно, в python уже встроены алгоритмы генерации многих распределений и для многих простых случаев они подойдут, но в случае практических задач придётся научиться генерировать точки из нужного нам распределения.

Пример

Давайте разберём тривиальный пример. Пусть мы хотим сосчитать такой интеграл:

\[ I=\int_{-\infty}^{\infty} d x \exp \left(-|x|^3\right)=2 \Gamma\left(\frac{4}{3}\right) \approx 1.786 \]

(тут \(\Gamma(x)\) - так называемая Гамма-функция Эйлера). Пусть мы хотим вычислить интеграл с точностью \(\epsilon=10^{-3}\). Действовать будем «на глаз».

Уменьшение области интегрирования

Построим подынтегральную функцию, и пристально на неё посмотрим. Судя по графику, кажется, что существенно отлична от нуля примерно на отрезке \((-2,2)\). Это будет нашей отправной точкой.

Сама функция на краях этого отрезка уже достаточно мала: \(e^{-8} \approx 3 \cdot 10^{-4}\), это даёт нам наивную оценку сверху для погрешности от такого приближения (точная оценка даёт \(2 \int_2^{\infty} \exp \left(-x^3\right) d x \approx 5 \cdot 10^{-5}\), но чтобы её получить, нужно уметь этот интеграл считать…). Такая погрешность нас вполне удовлетворяет, поэтому мы реализуем следующий алгоритм:

1. Выбираем \(N\) случайных точек \(x_i\) из отрезка \((-2,2)\)

2. Вычисляем \(I_N=4 \frac{1}{N} \sum_{i=1}^N f\left(x_i\right)\).

Погрешность этого алгоритма будет равна

\[ \left\langle\left\langle I_N^2\right\rangle\right\rangle=\frac{16}{N}\left\langle\left\langle f^2(x)\right\rangle\right\rangle=\frac{16}{N}\left(\int_{-2}^2 \frac{d x}{4} f^2(x)-\left(\int_{-2}^2 \frac{d x}{4} f(x)\right)^2\right) \approx \frac{2.48}{N} \]

Выбор функции распределения

А давайте теперь вместо этого попробуем приблизить эту функцию при помощи нормального распределения \(P(x)=e^{-x^2 / 2 \sigma^2} / \sqrt{2 \pi \sigma^2}\).

Дисперсию этого распределения \(\sigma\) мы опять выберем «на глаз», чтобы функция распределения была максимально похожей на подынтегральное выражение. Автору хорошим приближением показалось значение \(\sigma=0.5\). Так и поступим:

1. Выбираем \(N\) случайных точек \(x_i\) из нормального распределения с нулевым средним и дисперсией \(\sigma=0.5\) ( распределение \(P(x)\) ).

2. Вычисляем \(I_N=\frac{1}{N} \sum_{i=1}^N \frac{f\left(x_i\right)}{P\left(x_i\right)}\).

Погрешность этого алгоритма будет равна

\[ \left\langle\left\langle I_N^2\right\rangle\right\rangle=\frac{1}{N}\left\langle\left\langle\left(\frac{f(x)}{P(x)}\right)^2\right\rangle\right\rangle=\frac{1}{N}\left(\int_{-\infty}^{\infty} d x P(x) \frac{f^2(x)}{P^2(x)}-\left[\int_{-\infty}^{\infty} d x P(x) \frac{f(x)}{P(x)}\right]^2\right) \approx \frac{0.43}{N} \]

Ключевой вывод, который мы можем сделать из этого примера: выбирая правильным образом распределение \(P(\boldsymbol{r})\), можно заметным образом (в нашем случае - почти в шесть раз!) ускорить алгоритм.

Лирика: когда имеет смысл использовать Монте-Карло алгоритмы?

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

\[ \delta I / I \propto \epsilon \sim 1 / N \]

Тогда как в методе Монте-Карло:

\[ \delta I / I \propto \epsilon \sim 1 / \sqrt{N} \]

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

Метод Монте-Карло нужен, когда интегрирование происходит в многомерном пространстве достаточно большой размерности \(D\).

Ведь если мы хотим использовать квадратуру, то нам будет необходимо дискретизовать вдоль каждой пространственной размерности; и если мы апроксимируем интеграл, используя \(N\) точек, то каждое пространственное измерение мы будем вынуждены дискретизовать, вводя гораздо меньше - порядка \(N^{1 / D}\) - точек.

Как следствие, общая погрешность квадратурного метода будет вести себя уже как \(\delta I / I \propto \frac{1}{N^{m / D}}\), и при достаточно больших \(D>2 m\) Монте-Карло метод будет более эффективным (на него \(D\) не влияет никак в ассимптотике)

Поэтому ключевой вывод, который отсюда стоит вынести: метод Монте-Карло эффективен при интегрировании в пространстве достаточно большой размерности. И как раз основное применение методов Монте-Карло - случай, когда размерность пространства оказывается астрономически большой (например, фазовое пространство системы большого количество частиц).