Участник:Grachovamc/Интегрирование методом Монте - Карло: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показано 50 промежуточных версий этого же участника)
Строка 1: Строка 1:
 +
Основной автор описания: А.С.Грачева.
 +
 
== Свойства и структура алгоритма ==
 
== Свойства и структура алгоритма ==
  
Строка 14: Строка 16:
 
Предположим, требуется взять определенный интеграл от некоторой функции. Его значение будет численно равно площади под графиком данной функции.
 
Предположим, требуется взять определенный интеграл от некоторой функции. Его значение будет численно равно площади под графиком данной функции.
  
−
Для определения этой площади можно разбить отрезок, по которому производится интегрирование, на подотрезки, подсчитать площадь под графиком функции на каждом из них и сложить. И если предположить, что в случае функции одной переменной достаточно разбиения на 50 подотрезков и, соответственно, 50  вычислений значений функции, то в случае функции <math>n</math> переменных количество подотрезков и вычислений значений функции возрастает до <math>50^n</math>. При размерности функции больше 10 задача становится очень громоздкой. Поскольку пространства большой размерности встречаются во многих задачах, необходимо иметь метод решения, вычислительная сложность которого бы не столь сильно зависела от размерности. Именно таким свойством обладает метод Монте-Карло.
+
Для определения этой площади можно разбить отрезок, по которому производится интегрирование, на подотрезки, подсчитать площадь под графиком функции на каждом из них и сложить. И если предположить, что в случае функции одной переменной достаточно разбиения на 50 подотрезков и, соответственно, 50  вычислений значений функции, то в случае функции <math>n</math> переменных количество подотрезков и вычислений значений функции возрастает до <math>50^n</math>. Таким образом, при большой размерности функции задача становится очень громоздкой. Поскольку пространства большой размерности встречаются во многих задачах, необходимо иметь метод решения, вычислительная сложность которого бы не столь сильно зависела от размерности. Именно таким свойством обладает метод Монте-Карло.
  
 
=== Математическое описание алгоритма ===
 
=== Математическое описание алгоритма ===
Строка 24: Строка 26:
 
<math>I = \int_a^b \phi(x) dx</math>  
 
<math>I = \int_a^b \phi(x) dx</math>  
  
−
Пусть <math>X</math>  - случайная величина, равномерно распределенная на отрезке интегрирования <math>[a, b]</math> с плотностью <math>f(x) = \frac{1}{b - a}</math>. Тогда <math>\phi(X)</math> также является случайной величиной, и ее математическое ожидание выражается как
+
Пусть <math>X</math>  - случайная величина, равномерно распределенная на отрезке интегрирования <math>[a, b]</math> с плотностью <math>f(x) = \frac{1}{b - a}</math>.
 +
 
 +
Тогда <math>\phi(X)</math> также является случайной величиной, и ее математическое ожидание выражается как
  
 
<math>\mathbb{E}\phi(X)=\int\limits_a^b \phi(x)f(x) dx = \frac{1}{b - a}\int\limits_a^b \phi(x) dx</math>
 
<math>\mathbb{E}\phi(X)=\int\limits_a^b \phi(x)f(x) dx = \frac{1}{b - a}\int\limits_a^b \phi(x) dx</math>
Строка 43: Строка 47:
 
<math>I = \int_a^b \phi(x) dx</math>, где подынтегральная функция неотрицательна и ограничена <math>0 \le \phi(x) \le c</math> для любого <math>x \in [a, b]</math>.
 
<math>I = \int_a^b \phi(x) dx</math>, где подынтегральная функция неотрицательна и ограничена <math>0 \le \phi(x) \le c</math> для любого <math>x \in [a, b]</math>.
  
−
Согласно геометрическому определению интеграла для вычисления значения интеграла определим площадь под графиком интегрируемой функции. Для этого применяется следующий стохастический алгоритм:
+
Согласно геометрическому определению интеграла для вычисления значения интеграла определим площадь под графиком интегрируемой функции.  
 +
 
 +
Для этого применяется следующий стохастический алгоритм:
  
 
* функция ограничивается прямоугольником площадью <math>S</math>;
 
* функция ограничивается прямоугольником площадью <math>S</math>;
 
* в данном прямоугольнике случайным образом выбираются координаты <math>N</math> точек;
 
* в данном прямоугольнике случайным образом выбираются координаты <math>N</math> точек;
−
* определяется число <math>n</math> точек, координаты которых попали под график функции; тогда выполняется соотношение <math>\frac{\int\limits_a^b \phi(x) dx}{S} \approx \frac{n}{N}</math>
+
* определяется число точек <math>n</math>, координаты которых попали под график функции; тогда выполняется соотношение <math>\frac{\int\limits_a^b \phi(x) dx}{S} \approx \frac{n}{N}</math>
 
* искомая оценка интеграла равна <math>S\frac{n}{N}</math>.
 
* искомая оценка интеграла равна <math>S\frac{n}{N}</math>.
  
Строка 54: Строка 60:
 
===== Метод усреднения подынтегральной функции =====
 
===== Метод усреднения подынтегральной функции =====
  
−
Пусть функция <math>y = f(x_1, x_2,...,x_n)</math> непрерывна в ограниченной замкнутой области D и требуется вычислить n-кратный интеграл I по области D:
+
Пусть функция <math>y = f(x_1, x_2,...,x_n)</math> непрерывна в ограниченной замкнутой области <math>D</math>, и требуется вычислить <math>n</math>-кратный интеграл <math>I</math> по области <math>D</math>:
  
 
<math>I = {\int\limits_{D}{\ldots \int{f\left( {{x}_{1}},\ldots ,{{x}_{n}} \right)d{{x}_{1}}\ldots d{{x}_{n}}}}}</math>
 
<math>I = {\int\limits_{D}{\ldots \int{f\left( {{x}_{1}},\ldots ,{{x}_{n}} \right)d{{x}_{1}}\ldots d{{x}_{n}}}}}</math>
  
−
Геометрически число <math>I</math> представляет собой <math>(m + 1)</math> - мерный объем вертикального цилиндрического тела в пространстве <math>Ox_1x_2...x_ny</math>, построенного на основании D и ограниченного сверху соответствующим куском поверхности <math>y = f(\vec x)</math>,где <math>\vec x = (x_1, ..., x_n)</math>.
+
Геометрически число <math>I</math> представляет собой <math>(m + 1)</math> - мерный объем вертикального цилиндрического тела в пространстве <math>Ox_1x_2...x_ny</math>, построенного на основании <math>D</math> и ограниченного сверху соответствующим куском поверхности <math>y = f(\vec x)</math>,где <math>\vec x = (x_1, ..., x_n)</math>.
  
−
Пусть область интегрирования D расположена в <math>n</math>-мерном параллелепипеде <math>a_k \le x_k \le b_k</math>, <math>k = 1,...,n</math>.
+
Пусть область интегрирования <math>D</math> расположена в <math>n</math>-мерном параллелепипеде <math>a_k \le x_k \le b_k</math>, <math>k = 1,...,n</math>.
−
Преобразуем интеграл таким образом, чтобы новая область интегрирования целиком находилась внутри единичного <math>n</math> - мерного куба. Для этого сделаем замену переменных <math>x_k = a_k + (b_k - a_k) y_k</math>, <math>k = 1,...,n</math>. Тогда <math>0 \le y_k \le 1</math>, <math>k = 1,...,n</math>. Следовательно новая область интегрирования G целиком расположена в единичном <math>n</math>-мерном кубе.
+
Преобразуем интеграл таким образом, чтобы новая область интегрирования целиком находилась внутри единичного <math>n</math> - мерного куба. Для этого сделаем замену переменных <math>x_k = a_k + (b_k - a_k) y_k</math>, <math>k = 1,...,n</math>. Тогда <math>0 \le y_k \le 1</math>, <math>k = 1,...,n</math>. Следовательно новая область интегрирования <math>G</math> целиком расположена в единичном <math>n</math>-мерном кубе.
  
 
Таким образом, получаем выражение <math>I = \frac{D(x_1,...x_n)}{D(y_1,...y_n)}I'</math>, где якобиан преобразования имеет вид  
 
Таким образом, получаем выражение <math>I = \frac{D(x_1,...x_n)}{D(y_1,...y_n)}I'</math>, где якобиан преобразования имеет вид  
Строка 69: Строка 75:
 
<math>I' = {\int\limits_{G}{\ldots \int{g\left( {{y}_{1}},\ldots ,{{y}_{n}} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}} = {\int\limits_{G}{\ldots \int{f\left( {a_1 + (b_1 - a_1)y_1},\ldots ,{a_n + (b_n - a_n)y_n} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}}</math>,
 
<math>I' = {\int\limits_{G}{\ldots \int{g\left( {{y}_{1}},\ldots ,{{y}_{n}} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}} = {\int\limits_{G}{\ldots \int{f\left( {a_1 + (b_1 - a_1)y_1},\ldots ,{a_n + (b_n - a_n)y_n} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}}</math>,
  
−
а область интегрирования G расположена внутри <math>n</math> - мерного единичного куба <math>0 \le y_k \le 1</math>, <math>k = 1,...,n</math>.
+
а область интегрирования <math>G</math> расположена внутри <math>n</math> - мерного единичного куба <math>0 \le y_k \le 1</math>, <math>k = 1,...,n</math>.
 +
 
 +
Выберем <math>N</math> равномерно распределенных на отрезке <math>[0, 1]</math> последовательностей случайных величин <math>\vec y^{(i)} = ({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>, <math>i = 1,...,N</math>. Эти последовательности являются координатами случайных точек <math>K_i({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>.
 +
 
 +
Далее определим, какие из этих точек принадлежат области интегрирования <math>G</math>.
 +
 
 +
Пусть таких точек всего <math>m</math>. Тогда при достаточно большом <math>m</math> можно считать, что <math>g_{mean} = \mathbb{E}g(K) \approx \tilde{g}_{mean} = \frac{1}{m}\sum^{m}_{i=1}g(K_i) </math>. Следовательно оценка интеграла <math>I'</math> имеет вид
 +
<math>I' = \tilde{g}_{mean} * V_G</math>, где <math>V_G</math> - <math>n</math> - мерный объем области интегрирования <math>G</math>.
 +
 
 +
Из определения геометрической вероятности примем <math>V_G \approx \tilde{V}_G = \frac{m}{N}</math>.
  
−
Выберем N равномерно распределенных на отрезке <math>[0, 1]</math> последовательностей случайных величин <math>\vec y^{(i)} = ({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>, <math>i = 1,...,N</math>. Эти последовательности являются координатами случайных точек <math>K_i({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>. Далее определим, какие из этих точек принадлежат области интегрирования G. Пусть таких точек всего m. Тогда при достаточно большом m можно считать, что <math>\mathbb{E} g(K) = \frac{1}{m}\sum^{m}_{i=1}g(K_i) </math>. Следовательно оценка интеграла <math>I'</math> имеет вид
+
Таким образом, оценка искомого интеграла имеет вид <math>I = V_{par} * \tilde{g}_{mean} * \tilde{V}_G = \frac{V_{par}}{N}\sum^{m}_{i=1}g(K_i)</math>, где <math>V_{par} = (b_1 - a_1)...(b_n - a_n)</math> - объем параллепипеда, ограничивающего исходную область интегрирования <math>D</math>.
−
<math>I' = \mathbb{E} g(K) * V_G</math>, где <math>V_G</math> - <math>n</math> - мерный объем области интегрирования G.
 
−
Таким образом, оценка искомого интеграла имеет вид <math>I = V_{par} * \mathbb{E} g(K) * V_G = \frac{V_{par}}{N}\sum^{m}_{i=1}g(K_i)</math>, где <math>V_{par} = (b_1 - a_1)...(b_n - a_n)</math> - объем параллепипеда, ограничивающего исходную область интегрирования D.
 
  
 
===== Геометрический метод =====
 
===== Геометрический метод =====
 
В случае использования геометрического метода осуществляется та же последовательность действий, что и при использовании данного метода в случае функции одной переменной с учетом того, что функция ограничивается <math>n</math> - мерным параллелепипедом, и, соответственно, в результирующей формуле вместо площади прямоугольника берется объем ограничивающего функцию параллелепипеда.
 
В случае использования геометрического метода осуществляется та же последовательность действий, что и при использовании данного метода в случае функции одной переменной с учетом того, что функция ограничивается <math>n</math> - мерным параллелепипедом, и, соответственно, в результирующей формуле вместо площади прямоугольника берется объем ограничивающего функцию параллелепипеда.
 +
 
=== Вычислительное ядро алгоритма ===
 
=== Вычислительное ядро алгоритма ===
 
Таким образом, вычислительным ядром последовательной реализации любого из представленных вариантов метода Монте - Карло интегрирования является многократное вычисление значения функции в случайной точке.
 
Таким образом, вычислительным ядром последовательной реализации любого из представленных вариантов метода Монте - Карло интегрирования является многократное вычисление значения функции в случайной точке.
 
=== Макроструктура алгоритма ===
 
=== Макроструктура алгоритма ===
−
Согласно описанию вычислительного ядра алгоритма в качестве макровершины можно взять вычисление подынтегральной функции в случайной точке. Данное действие выполняется в алгоритме N раз.
+
В случае использования метода усреднения подынтегральной функции начальным этапом работы алгоритма является генерация случайных точек, соответствующих начальным условиям, то есть подынтегральной функции, пределам интегрирования и количеству генерируемых точек. Время на выполнение первого шага оценивается константой. Затем производится вычисление подынтегральной функции в полученных точках, которое согласно описанию вычислительного ядра алгоритма можно принять в качестве макровершины и не задумываться далее о виде и особенностях этой функции, которые будут определять время выполнения этого этапа. Далее производится суммирование результатов предыдущего шага и умножение полученной суммы на соответствующую константу. Также в качестве макровершины можно взять в целом суммирование значений функции в случайных точках и умножение этой суммы на константу, но так как в большинстве случаев время, которое уходит на вычисление подынтегральной функции гораздо больше, чем время, требуемое на осуществление операций суммирования и умножения, то основное время работы все равно приходится на вычисление подынтегральной функции.
 +
 
 +
В случае использования геометрического метода первый этап аналогичен начальному этапу метода усреднения подынтегральной функции с тем лишь отличием, что кроме случайных точек, в которых в дальнейшем будет вычисляться функция, генерируются и значения самой функции таким образом, чтобы полученная точка принадлежала параллелепипеду, ограничивающему функцию. Затем также производится вычисление подынтегральной функции в полученных точках, которое согласно описанию вычислительного ядра алгоритма можно принять в качестве макровершины. Далее пересчитывается количество точек, попавших под график исходной функции, и полученный результат умножается на соответствующую константу.
 
=== Схема реализации последовательного алгоритма ===
 
=== Схема реализации последовательного алгоритма ===
 
==== Интегрирование функции одной переменной ====
 
==== Интегрирование функции одной переменной ====
Строка 88: Строка 104:
 
1. Генерация последовательности <math>{{\{x_i\}}_{i = 1}}^{N}</math> значений случайной величины, равномерно распределенной на отрезке интегрирования <math>[a, b]</math>. Для экономии памяти лучше генерировать не всю случайную последовательность целиком, а по одному значению в цикле.  
 
1. Генерация последовательности <math>{{\{x_i\}}_{i = 1}}^{N}</math> значений случайной величины, равномерно распределенной на отрезке интегрирования <math>[a, b]</math>. Для экономии памяти лучше генерировать не всю случайную последовательность целиком, а по одному значению в цикле.  
  
−
2. Вычисление значений функции <math>f(x_i)</math> в каждой из точек последовательности <math>{{\{x_i\}}_{i = 1}}^{N}</math>.
+
2. Вычисление значений функции <math>\phi(x_i)</math> в каждой из точек последовательности <math>{{\{x_i\}}_{i = 1}}^{N}</math>.
  
−
3. Вычисление выборочного среднего <math>\frac{1}{N}\sum^{N}_{i=1}f(x_i)</math>.
+
3. Вычисление выборочного среднего <math>\frac{1}{N}\sum^{N}_{i=1}\phi(x_i)</math>.
  
 
4. Получение искомого значения интеграла согласно формуле <math>\int\limits_{a}^{b}\phi(x)\,dx\approx\frac{b-a}{N}\sum^{N}_{i=1}\phi(x_i)</math>.
 
4. Получение искомого значения интеграла согласно формуле <math>\int\limits_{a}^{b}\phi(x)\,dx\approx\frac{b-a}{N}\sum^{N}_{i=1}\phi(x_i)</math>.
Строка 119: Строка 135:
 
Последовательность выполнения метода следующая:
 
Последовательность выполнения метода следующая:
  
−
1. Преобразование интеграла таким образом, чтобы новая область интегрирования G целиком находилась внутри единичного <math>n</math> - мерного куба с использованием замены переменных <math>x_k = a_k + (b_k - a_k) y_k</math>, <math>k = 1,...,n</math>.
+
1. Преобразование интеграла таким образом, чтобы новая область интегрирования <math>G</math> целиком находилась внутри единичного <math>n</math> - мерного куба с использованием замены переменных <math>x_k = a_k + (b_k - a_k) y_k</math>, <math>k = 1,...,n</math>.
  
 
2. Вычисление якобиана преобразования и получение следующего выражения получаем выражение <math>I = \frac{D(x_1,...x_n)}{D(y_1,...y_n)}I'</math>, где якобиан преобразования имеет вид  
 
2. Вычисление якобиана преобразования и получение следующего выражения получаем выражение <math>I = \frac{D(x_1,...x_n)}{D(y_1,...y_n)}I'</math>, где якобиан преобразования имеет вид  
Строка 125: Строка 141:
 
<math>\frac{D(x_1,...x_n)}{D(y_1,...y_n)} = (b_1 - a_1)...(b_n - a_n)</math>, <math>I' = {\int\limits_{G}{\ldots \int{g\left( {{y}_{1}},\ldots ,{{y}_{n}} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}} = {\int\limits_{G}{\ldots \int{f\left( {a_1 + (b_1 - a_1)y_1},\ldots ,{a_n + (b_n - a_n)y_n} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}}</math>.
 
<math>\frac{D(x_1,...x_n)}{D(y_1,...y_n)} = (b_1 - a_1)...(b_n - a_n)</math>, <math>I' = {\int\limits_{G}{\ldots \int{g\left( {{y}_{1}},\ldots ,{{y}_{n}} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}} = {\int\limits_{G}{\ldots \int{f\left( {a_1 + (b_1 - a_1)y_1},\ldots ,{a_n + (b_n - a_n)y_n} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}}</math>.
  
−
3.Выбор N равномерно распределенных на отрезке <math>[0, 1]</math> последовательностей случайных величин <math>\vec y^{(i)} = ({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>, <math>i = 1,...,N</math>, которые являются координатами случайных точек <math>K_i({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>.  
+
3.Выбор <math>N</math> равномерно распределенных на отрезке <math>[0, 1]</math> последовательностей случайных величин <math>\vec y^{(i)} = ({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>, <math>i = 1,...,N</math>, которые являются координатами случайных точек <math>K_i({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})</math>.
 +
 
 +
4. Определение того, какие из этих точек принадлежат области интегрирования <math>G</math>. Пусть таких точек всего <math>m</math>.
  
−
4. Определение того, какие из этих точек принадлежат области интегрирования G. Пусть таких точек всего m.
+
5. Вычисление <math>g_{mean} = \mathbb{E}g(K) \approx \tilde{g}_{mean} = \frac{1}{m}\sum^{m}_{i=1}g(K_i) </math>.
  
−
5. Вычисление <math>\mathbb{E} g(K) = \frac{1}{m}\sum^{m}_{i=1}g(K_i) </math>.
+
6. Получение оценки интеграла <math>I' = \tilde{g}_{mean} * V_G</math>, где <math>V_G</math> - <math>n</math> - мерный объем области интегрирования <math>G</math>. Считаем, что <math>V_G \approx \tilde{V}_G = \frac{m}{N}</math>.
  
−
6. Получение оценки интеграла <math>I'</math>: <math>I' = \mathbb{E} g(K) * V_G</math>, где <math>V_G</math> - <math>n</math> - мерный объем области интегрирования G.
+
7. Получение оценки искомого интеграла <math>I = V_{par} * \tilde{g}_{mean} * \tilde{V}_G = \frac{V_{par}}{N}\sum^{m}_{i=1}g(K_i)</math>, где <math>V_{par} = (b_1 - a_1)...(b_n - a_n)</math> - объем параллелепипеда, ограничивающего исходную область интегрирования <math>D</math>.
  
−
7. Получение оценки искомого интеграла <math>I = V_{par} * \mathbb{E} g(K) * V_G = \frac{V_{par}}{N}\sum^{m}_{i=1}g(K_i)</math>, где <math>V_{par} = (b_1 - a_1)...(b_n - a_n)</math> - объем параллепипеда, ограничивающего исходную область интегрирования D.
 
 
==== Геометрический метод ====
 
==== Геометрический метод ====
 
Схема реализации последовательного алгоритма аналогична одномерному случаю.
 
Схема реализации последовательного алгоритма аналогична одномерному случаю.
Строка 142: Строка 159:
 
1) при использовании метода усреднения подынтегральной функции
 
1) при использовании метода усреднения подынтегральной функции
  
−
* N генераций случайных величин;
+
* <math>N</math> генераций случайных величин;
−
* N вычислений значений функции;
+
* <math>N</math> вычислений значений функции;
−
* N - 1 сложение.
+
* <math>N - 1</math> сложение.
  
 
2) при использовании геометрического метода
 
2) при использовании геометрического метода
  
−
* N генераций случайных величин;
+
* <math>N</math> генераций случайных величин;
−
* N вычислений значений функции;
+
* <math>N</math> вычислений значений функции;
−
* N сравнений.
+
* <math>N</math> сравнений.
  
 
Таким образом, алгоритм является линейным относительно количества случайных точек.
 
Таким образом, алгоритм является линейным относительно количества случайных точек.
 +
 
=== Информационный граф ===
 
=== Информационный граф ===
−
[[Файл:My perfect graph.png|thumb|center|723px|Рисунок 1. Граф алгоритма при использовании метода усреднения подынтегральной функции]]
+
==== Метод усреднения подынтегральной функции ====
−
 
+
Граф алгоритма состоит из трех групп вершин. Вершины на верхнем ярусе используют в качестве входных данных сгенерированные случайные значения и вычисляют в этих точках функцию <math>f</math>. Вершины второго яруса суммируют полученные значения. На третьем ярусе происходит умножение полученной на предыдущем ярусе суммы на соответствующую константу, формула для вычисления которой представлена ранее при описании алгоритма.
 +
[[Файл:Graph mean.jpg|thumb|center|723px|Рисунок 1. Граф алгоритма при использовании метода усреднения подынтегральной функции]]
 +
==== Геометрический метод ====
 +
Граф алгоритма состоит из четырех групп вершин. При данном методе сначала производится ограничение функции <math>f</math> параллелепипедом и генерация соответствующих случайных точек. Далее полученные значения подаются на вход вершинам верхнего яруса, которые вычисляют значение функции <math>f</math> в них. Вершины второго яруса производят сравнивание, в зависимости от которого на последующем ярусе вычисляется количество точек, попавших под график функции. Затем производится умножение полученного результата на соответствующую константу, описанную ранее.
 
[[Файл:Graph 2.png|thumb|center|810px|Рисунок 2. Граф алгоритма при использовании геометрического метода]]
 
[[Файл:Graph 2.png|thumb|center|810px|Рисунок 2. Граф алгоритма при использовании геометрического метода]]
−
 
  
 
=== Ресурс параллелизма алгоритма ===
 
=== Ресурс параллелизма алгоритма ===
  
 
Для вычисления интеграла методом Монте - Карло при использовании усреднения подынтегральной функции требуется последовательно выполнить следующие ярусы:
 
Для вычисления интеграла методом Монте - Карло при использовании усреднения подынтегральной функции требуется последовательно выполнить следующие ярусы:
−
* 1 ярус генераций случайных величин (N * p генераций, где N - изначально выбранное количество точек, p - кратность интеграла);
 
 
* 1 ярус вычислений значений функции в полученных ранее случайных точках (N вычислений);
 
* 1 ярус вычислений значений функции в полученных ранее случайных точках (N вычислений);
 
* 1 ярус сложений (N - 1 сложение);
 
* 1 ярус сложений (N - 1 сложение);
Строка 168: Строка 187:
  
 
Для вычисления интеграла методом Монте - Карло при использовании геометрической интерпретации требуется последовательно выполнить следующие ярусы:
 
Для вычисления интеграла методом Монте - Карло при использовании геометрической интерпретации требуется последовательно выполнить следующие ярусы:
−
* 1 ярус генераций случайных величин (N * p генераций, где N - изначально выбранное количество точек, p - кратность интеграла);
 
 
* 1 ярус вычислений значений функции в полученных ранее случайных точках (N вычислений);
 
* 1 ярус вычислений значений функции в полученных ранее случайных точках (N вычислений);
 
* 1 ярус сравнений (N сравнений);
 
* 1 ярус сравнений (N сравнений);
Строка 176: Строка 194:
 
Таким образом, алгоритм отлично поддается распараллеливанию, так как все вычисления функции f могут производиться независимо в общем случае.
 
Таким образом, алгоритм отлично поддается распараллеливанию, так как все вычисления функции f могут производиться независимо в общем случае.
  
−
При классификации по высоте ярусно-параллельной формы алгоритм метода Монте - Карло имеет сложность 4 или 5 в зависимости от конкретного выбора одного из двух предлагаемых вариантов. При классификации по ширине его сложность будет <math>O(N)</math>.
+
При классификации по высоте ярусно-параллельной формы алгоритм метода Монте - Карло имеет сложность 3 или 4 в зависимости от конкретного выбора одного из двух предлагаемых вариантов. При классификации по ширине его сложность будет <math>O(N)</math>.
  
 
=== Входные и выходные данные алгоритма ===
 
=== Входные и выходные данные алгоритма ===
Строка 193: Строка 211:
 
Алгоритм интегрирования методом Монте - Карло не является детерминированным, так он главным образом основывается на генерации случайных величин, в связи с чем невозможно получить точный и предопределенный заранее результат. Поэтому для достижения необходимой точности выбирается допустимая погрешность.
 
Алгоритм интегрирования методом Монте - Карло не является детерминированным, так он главным образом основывается на генерации случайных величин, в связи с чем невозможно получить точный и предопределенный заранее результат. Поэтому для достижения необходимой точности выбирается допустимая погрешность.
  
−
Согласно теории погрешностей, используя формулу оценки погрешности произведения двух величин, получим следующее соотношение:
+
Выведем формулу для вычисления значения погрешности:
 +
 
 +
<math>\Delta I = \Delta(V_{par} * \tilde{g}_{mean} * \tilde{V_G}) = V_{par} * (\Delta g_{mean} * \tilde{V_G} + \Delta V_G * |\tilde{g}_{mean}|)</math>.
 +
 
 +
По центральной предельной теореме
 +
 
 +
<math> P \left\{ |g_{mean} - \tilde{g}_{mean}| \le t_{\beta}\sqrt{\frac{\sigma_{1}^2}{m}} \right\} \approx \beta </math>, <math> P \left\{ |V_G - \tilde{V}_G| \le t_{\beta}\sqrt{\frac{\sigma_{2}^2}{N}} \right\} \approx 2\Phi(t_{\beta}) = \beta </math>,
 +
 
 +
где <math>\beta</math> - коэффициент доверия, <math>\Phi(t)</math> - функция стандартного нормального распределения, <math>\sigma_{1}^2 = D[g(K)], \sigma_{2}^2 = D[V_G]</math>.
 +
 
 +
То есть с вероятностью <math>\beta</math>
 +
 
 +
<math>\Delta g_{mean} = t_{\beta} \sqrt{\frac{\sigma_{1}^2}{m}}, \Delta V_G = t_{\beta} \sqrt{\frac{\sigma_{2}^2}{N}}</math>.
 +
 
 +
При этом <math>\tilde{g}_{mean} = \frac{1}{m}\sum^{m}_{i=1}g(K_i), \tilde{V}_{G} = \frac{m}{N} </math>.
 +
 
 +
Таким образом,
 +
 
 +
<math>\Delta I = V_{par} * t_{\beta} (\sigma_1 * \frac{\tilde{V}_G}{m} + \sigma_2 * \frac{|\tilde{g}_{mean}|}{N})</math>.
 +
 
 +
Обычно величины <math>\sigma_1^2</math> и <math>\sigma_2^2</math> оказываются неизвестными, тогда их приближенные значения вычисляются по формулам
 +
 
 +
<math>\tilde{\sigma}_1^2 \approx (\tilde{g}^2)_{mean} - (|\tilde{g}_{mean}|)^2</math>, где <math>(\tilde{g}^2)_{mean} = \frac{1}{m}\sum^{m}_{i=1}g^2(K_i)</math>;
 +
 
 +
<math>\tilde{\sigma}_2^2 \approx \tilde{V}_G (1 - \tilde{V}_G)</math>.
 +
 
 +
Таким образом, получаем следующее соотношение для погрешности интегрирования методом Монте - Карло:
 +
 
 +
<math>\Delta I = V_{par} * t_{\beta} (\tilde{\sigma}_1 * \frac{\sqrt{\tilde{V}_G}}{\sqrt{N}} + \tilde{\sigma}_2 * \frac{|\tilde{g}_{mean}|}{\sqrt{N}})</math>.
 +
 
 +
Если полученное значение погрешности велико, то можно назначить желаемую погрешность, на ее основании вычислить необходимое количество испытаний и затем определить значение интеграла с требуемым количеством испытаний.
  
−
<math>\Delta I = \Delta(</math>
+
Дуги информационного графа являются полностью локальными.
  
 
== Программная реализация алгоритма ==
 
== Программная реализация алгоритма ==
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
 +
==== Метод усреднения подынтегральной функции ====
 +
 +
Код последовательной реализации для случая использования метода усреднения подынтегральной функции при простейшем предположении, что предел интегрирования по каждой из переменных  расположен внутри отрезка [0,1], на языке C может выглядеть так:
 +
 +
<pre>
 +
double mean_method(long int number, int order, double *limits)
 +
{
 +
    double coord[order], res;
 +
    int n, i, ind;
 +
    res = 0.0;
 +
    for (n = 0; n < number; n++) {
 +
        ind = 0;
 +
        for (i = 0; i < order; i++) {
 +
            coord[i] = (double)(rand())/RAND_MAX;
 +
            if (coord[i] > limits[2 * i + 1] || coord[i] < limits[2 * i]) {
 +
                ind = 1;
 +
                break;
 +
            }
 +
        }
 +
        if (ind == 0) {
 +
            res = res + f(coord, order);
 +
        }
 +
    }                                 
 +
    return res;
 +
}
 +
</pre>
 +
 +
==== Геометрический метод ====
 +
 +
Код последовательной реализации для случая использования геометрического метода при простейшем предположении, что функцию можно ограничить параллелепипедом, каждая из сторон которого является отрезком [0,1], на языке C может выглядеть так:
 +
 +
<pre>
 +
double geometric_method(long int number, int order)
 +
{
 +
    double coord[order + 1];
 +
    int score, n, i;
 +
    score = 0;
 +
    for (n = 0; n < number; n++) {
 +
        for (i = 0; i <= order; i++) {
 +
            coord[i] = (double)(rand())/RAND_MAX;       
 +
        }
 +
       
 +
        if (coord[order] <= f(coord, order))
 +
            score++;
 +
    }
 +
    return score;
 +
}
 +
</pre>
 +
 
=== Локальность данных и вычислений ===
 
=== Локальность данных и вычислений ===
 +
 +
Локальность вычислений существенно определяется реализацией вычисления функции <math>f</math>. В предположении, что данная реализация уже имеется, будем считать, что в ней учтены требования эффективного использования памяти.
 +
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
Самым очевидным вариантом реализации параллельного алгоритма является независимое вычислении функции f в генерируемых случайных точках различными вычислительными узлами. Это может быть реализовано как при помощи технологии threads или OpenMP при запуске на одном узле, так и при использовании межпроцессорных обменов при помощи MPI. Головной узел передает остальным вычислительным узлам реализацию функции f и точки или наборы точек, в которых должна вычисляться функция.
 +
 +
Однако при такой реализации нужно уделять особое внимание генерации случайных чисел, которая должна быть независимой на каждом вычислительном узле.
 +
Поскольку при вычислениях методом Монте-Карло существенное количество операций расходуется для оперирования над случайными числами, то наличие простых и экономных способов формирования последовательности случайных чисел во многом определяет возможность практического использования этого метода. Обычно для подобных задач используется генерация равномерно распределенных на интервале (0, 1) случайных величин, так как такой способ требует мало машинного времени и обеспечивает простоту дальнейших преобразований. Но в связи с тем, что в языках программирования используются генераторы псевдослучайных чисел, важно, чтобы они имели достаточно длинный период, гарантирующий отсутствие зацикливания последовательности в пределах решаемой задачи.
 +
 +
Для более эффективной и быстрой работы алгоритма возможно использование более сложных реализаций, которые будут требовать уточнения класса подынтегральных функций и следовательно соответствующей модификации алгоритма их вычисления.
 +
 
=== Масштабируемость алгоритма и его реализации ===
 
=== Масштабируемость алгоритма и его реализации ===
 +
Согласно теоретическим оценкам алгоритмы должны быть отлично масштабируемы. Докажем это предположение на практике.
 +
 +
Для этого были созданы реализации описанных методов. В качестве тестового примера рассматривалось интегрирование функции трех переменных <math>\sqrt{x_1} + \sqrt{x_2} - \sqrt{x_3}</math>. Использование метода Монте - Карло в данном случае нерационально, поскольку оптимальные выбрать другой детерминированный метод, но так как в данной задаче вычисление подынтегральной функции рассматривается как "черный ящик", следовательно сложность вычисления функции не имеет значения.
 +
 +
Выбранный алгоритм использует в качестве метода распараллеливания технологию MPI. Для проведения экспериментов использовалась программа, написанная на языке C, с использованием библиотеки Intel MPI 5.0.1. Тесты проводились на суперкомпьютере Ломоносов<ref>http://users.parallel.ru/wiki/pages/22-config</ref> в узле test. Для компиляции программы использовалась команда mpicc, которая является надстройкой над имеющимися компиляторами и автоматически подключает заголовочные файлы и библиотеки MPI. Для демонстрации оказалось достаточным провести тесты с количеством процессоров до 128 с масштабированием по степени числа 2. Для количества точек, в которых вычислялась функция, была введена логарифмическая шкала: от 1000 до 1000000000 с шагом умножения на 10. Такая шкала была выбрана для того, чтобы полностью продемонстрировать постепенный рост производительности при увеличении числа процессоров и падение эффективности при небольших количествах точек интегрирования.
 +
 +
В случае рассматриваемого метода не представлялось возможным измерять производительность в гигафлопсах, поэтому для построения более наглядного графика в качестве показателя производительности использовалось количество вычислений функций в единицу времени, умноженное на 10000.
 +
 +
В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:
 +
 +
* минимальная эффективность реализации 1,08%;
 +
* максимальная эффективность реализации 99,8%.
 +
 +
 +
Ниже представлены графики производительности и эффективности реализации алгоритма интегрирования методом усреднения подынтегральной функции в зависимости от изменяемых параметров запуска. Так два рассматриваемых алгоритма интегрирования методом Монте - Карло имеют похожие реализации, то для геометрического метода графики будут аналогичными.
 +
 +
[[Файл:Perfomance 1.png|thumb|center|723px|Рисунок 3. Параллельная реализация интегрирования методом Монте - Карло. Изменение производительности в зависимости от числа процессоров и количества генераций случайных точек.]]
 +
 +
[[Файл:Efficiency 1.png|thumb|center|723px|Рисунок 4. Параллельная реализация интегрирования методом Монте - Карло. Изменение эффективности в зависимости от числа процессоров и количества генераций случайных точек.]]
 +
 +
Таким образом, очевидно, что эффективность для данной задачи существенно зависит от того, превышает ли количество генерируемых случайных точек количество процессоров хотя бы в 100 раз. Тогда эффективность реализации стремится к единице, и получаемое ускорение практически линейно. Из того, что в реальных задачах количество используемых процессоров даже близко не приближается к количеству генерируемых точек, следует вывод, что алгоритм действительно прекрасно масштабируем.
 +
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
=== Выводы для классов архитектур ===
 
=== Выводы для классов архитектур ===
 
=== Существующие реализации алгоритма ===
 
=== Существующие реализации алгоритма ===
 +
В рамках библиотеки GNU Scientific Library (или GSL)<ref>https://www.gnu.org/software/gsl/</ref>, которая была написана для численных вычислений в прикладной математике и науке, реализован алгоритм интегрирования функций многих переменных методом Монте - Карло<ref>https://www.gnu.org/software/gsl/doc/html/montecarlo.html</ref>. Данная библиотека является свободным программным обеспечением. Распространяется под лицензией GNU GPL 3.
 +
 +
Других реализаций обнаружено не было.
 +
 
== Литература ==
 
== Литература ==

Текущая версия на 15:22, 5 декабря 2017

Основной автор описания: А.С.Грачева.

Содержание

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Метод Монте-Карло — общее название группы численных методов, основанных на получении большого числа реализаций стохастического процесса, который формируется таким образом, чтобы его вероятностные характеристики совпадали с аналогичными величинами решаемой задачи. Метод Монте-Карло имеет множество различных приложений. Он применяется для решения задач в различных областях математики, физики, химии, биологии, экономики, оптимизации, теории управления и др.

Годом рождения метода Монте-Карло считается 1949 год, когда в свет выходит статья Метрополиса и Улама «Метод Монте-Карло»[1]. Название метода происходит от названия города в княжестве Монако, широко известного своими многочисленными казино, поскольку именно рулетка является одним из самых широко известных генераторов случайных чисел.

Особенностью данного метода является простая структура вычислительного алгоритма. Как правило, составляется программа для осуществления одного случайного испытания, затем это испытание повторяется N раз, причем каждый опыт не зависит от остальных. Для получения итогового значения усредняются результаты всех опытов. Поэтому метод Монте - Карло также иногда называют методом статистических испытаний[2].

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

Рассмотрим использование метода Монте — Карло для задачи вычисления интегралов.

Предположим, требуется взять определенный интеграл от некоторой функции. Его значение будет численно равно площади под графиком данной функции.

Для определения этой площади можно разбить отрезок, по которому производится интегрирование, на подотрезки, подсчитать площадь под графиком функции на каждом из них и сложить. И если предположить, что в случае функции одной переменной достаточно разбиения на 50 подотрезков и, соответственно, 50 вычислений значений функции, то в случае функции [math]n[/math] переменных количество подотрезков и вычислений значений функции возрастает до [math]50^n[/math]. Таким образом, при большой размерности функции задача становится очень громоздкой. Поскольку пространства большой размерности встречаются во многих задачах, необходимо иметь метод решения, вычислительная сложность которого бы не столь сильно зависела от размерности. Именно таким свойством обладает метод Монте-Карло.

1.2 Математическое описание алгоритма

1.2.1 Интегрирование функции одной переменной

Рассмотрим два наиболее распространенных варианта использования метода Монте - Карло для вычисления интегралов.

1.2.1.1 Метод усреднения подынтегральной функции

Рассмотрим интеграл вида

[math]I = \int_a^b \phi(x) dx[/math]

Пусть [math]X[/math] - случайная величина, равномерно распределенная на отрезке интегрирования [math][a, b][/math] с плотностью [math]f(x) = \frac{1}{b - a}[/math].

Тогда [math]\phi(X)[/math] также является случайной величиной, и ее математическое ожидание выражается как

[math]\mathbb{E}\phi(X)=\int\limits_a^b \phi(x)f(x) dx = \frac{1}{b - a}\int\limits_a^b \phi(x) dx[/math]

Отсюда

[math]\int_a^b \phi(x) dx = (b - a)\mathbb{E}\phi(X)[/math]

Заменив математическое ожидание [math]\mathbb{E}\phi(X)[/math] его оценкой – выборочным средним, получим оценку интеграла

[math]\int\limits_{a}^{b}\phi(x)\,dx\approx\frac{b-a}{N}\sum^{N}_{i=1}\phi(x_i)[/math],

где [math]x_i[/math] – возможные значения случайной величины [math]X[/math], [math]N[/math] – число испытаний.

1.2.1.2 Геометрический метод

Рассмотрим интеграл вида

[math]I = \int_a^b \phi(x) dx[/math], где подынтегральная функция неотрицательна и ограничена [math]0 \le \phi(x) \le c[/math] для любого [math]x \in [a, b][/math].

Согласно геометрическому определению интеграла для вычисления значения интеграла определим площадь под графиком интегрируемой функции.

Для этого применяется следующий стохастический алгоритм:

  • функция ограничивается прямоугольником площадью [math]S[/math];
  • в данном прямоугольнике случайным образом выбираются координаты [math]N[/math] точек;
  • определяется число точек [math]n[/math], координаты которых попали под график функции; тогда выполняется соотношение [math]\frac{\int\limits_a^b \phi(x) dx}{S} \approx \frac{n}{N}[/math]
  • искомая оценка интеграла равна [math]S\frac{n}{N}[/math].

1.2.2 Вычисление кратного интеграла

Случай вычисления [math]n[/math] - мерного интеграла аналогичен одномерному случаю с учетом некоторых корректировок.

1.2.2.1 Метод усреднения подынтегральной функции

Пусть функция [math]y = f(x_1, x_2,...,x_n)[/math] непрерывна в ограниченной замкнутой области [math]D[/math], и требуется вычислить [math]n[/math]-кратный интеграл [math]I[/math] по области [math]D[/math]:

[math]I = {\int\limits_{D}{\ldots \int{f\left( {{x}_{1}},\ldots ,{{x}_{n}} \right)d{{x}_{1}}\ldots d{{x}_{n}}}}}[/math]

Геометрически число [math]I[/math] представляет собой [math](m + 1)[/math] - мерный объем вертикального цилиндрического тела в пространстве [math]Ox_1x_2...x_ny[/math], построенного на основании [math]D[/math] и ограниченного сверху соответствующим куском поверхности [math]y = f(\vec x)[/math],где [math]\vec x = (x_1, ..., x_n)[/math].

Пусть область интегрирования [math]D[/math] расположена в [math]n[/math]-мерном параллелепипеде [math]a_k \le x_k \le b_k[/math], [math]k = 1,...,n[/math]. Преобразуем интеграл таким образом, чтобы новая область интегрирования целиком находилась внутри единичного [math]n[/math] - мерного куба. Для этого сделаем замену переменных [math]x_k = a_k + (b_k - a_k) y_k[/math], [math]k = 1,...,n[/math]. Тогда [math]0 \le y_k \le 1[/math], [math]k = 1,...,n[/math]. Следовательно новая область интегрирования [math]G[/math] целиком расположена в единичном [math]n[/math]-мерном кубе.

Таким образом, получаем выражение [math]I = \frac{D(x_1,...x_n)}{D(y_1,...y_n)}I'[/math], где якобиан преобразования имеет вид

[math]\frac{D(x_1,...x_n)}{D(y_1,...y_n)} = (b_1 - a_1)...(b_n - a_n)[/math],

[math]I' = {\int\limits_{G}{\ldots \int{g\left( {{y}_{1}},\ldots ,{{y}_{n}} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}} = {\int\limits_{G}{\ldots \int{f\left( {a_1 + (b_1 - a_1)y_1},\ldots ,{a_n + (b_n - a_n)y_n} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}}[/math],

а область интегрирования [math]G[/math] расположена внутри [math]n[/math] - мерного единичного куба [math]0 \le y_k \le 1[/math], [math]k = 1,...,n[/math].

Выберем [math]N[/math] равномерно распределенных на отрезке [math][0, 1][/math] последовательностей случайных величин [math]\vec y^{(i)} = ({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})[/math], [math]i = 1,...,N[/math]. Эти последовательности являются координатами случайных точек [math]K_i({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})[/math].

Далее определим, какие из этих точек принадлежат области интегрирования [math]G[/math].

Пусть таких точек всего [math]m[/math]. Тогда при достаточно большом [math]m[/math] можно считать, что [math]g_{mean} = \mathbb{E}g(K) \approx \tilde{g}_{mean} = \frac{1}{m}\sum^{m}_{i=1}g(K_i) [/math]. Следовательно оценка интеграла [math]I'[/math] имеет вид [math]I' = \tilde{g}_{mean} * V_G[/math], где [math]V_G[/math] - [math]n[/math] - мерный объем области интегрирования [math]G[/math].

Из определения геометрической вероятности примем [math]V_G \approx \tilde{V}_G = \frac{m}{N}[/math].

Таким образом, оценка искомого интеграла имеет вид [math]I = V_{par} * \tilde{g}_{mean} * \tilde{V}_G = \frac{V_{par}}{N}\sum^{m}_{i=1}g(K_i)[/math], где [math]V_{par} = (b_1 - a_1)...(b_n - a_n)[/math] - объем параллепипеда, ограничивающего исходную область интегрирования [math]D[/math].

1.2.2.2 Геометрический метод

В случае использования геометрического метода осуществляется та же последовательность действий, что и при использовании данного метода в случае функции одной переменной с учетом того, что функция ограничивается [math]n[/math] - мерным параллелепипедом, и, соответственно, в результирующей формуле вместо площади прямоугольника берется объем ограничивающего функцию параллелепипеда.

1.3 Вычислительное ядро алгоритма

Таким образом, вычислительным ядром последовательной реализации любого из представленных вариантов метода Монте - Карло интегрирования является многократное вычисление значения функции в случайной точке.

1.4 Макроструктура алгоритма

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

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

1.5 Схема реализации последовательного алгоритма

1.5.1 Интегрирование функции одной переменной

1.5.1.1 Метод усреднения подынтегральной функции

Последовательность исполнения метода следующая:

1. Генерация последовательности [math]{{\{x_i\}}_{i = 1}}^{N}[/math] значений случайной величины, равномерно распределенной на отрезке интегрирования [math][a, b][/math]. Для экономии памяти лучше генерировать не всю случайную последовательность целиком, а по одному значению в цикле.

2. Вычисление значений функции [math]\phi(x_i)[/math] в каждой из точек последовательности [math]{{\{x_i\}}_{i = 1}}^{N}[/math].

3. Вычисление выборочного среднего [math]\frac{1}{N}\sum^{N}_{i=1}\phi(x_i)[/math].

4. Получение искомого значения интеграла согласно формуле [math]\int\limits_{a}^{b}\phi(x)\,dx\approx\frac{b-a}{N}\sum^{N}_{i=1}\phi(x_i)[/math].

1.5.1.2 Геометрический метод

Пусть интегрируемая функция имеет следующий вид

0 step.png

Последовательность исполнения метода следующая:

1. Функция ограничивается прямоугольником площадью [math]S[/math]

1 step.png

2. В прямоугольнике случайным образом выбираются координаты [math]N[/math] точек

3 step.png

3. Определяется число [math]n[/math] точек, координаты которых попали под график функции

4 step.png

4. Искомая площадь вычисляется согласно формуле [math]S\frac{n}{N}[/math]

1.5.2 Вычисление кратного интеграла

1.5.2.1 Метод усреднения подынтегральной функции

Последовательность выполнения метода следующая:

1. Преобразование интеграла таким образом, чтобы новая область интегрирования [math]G[/math] целиком находилась внутри единичного [math]n[/math] - мерного куба с использованием замены переменных [math]x_k = a_k + (b_k - a_k) y_k[/math], [math]k = 1,...,n[/math].

2. Вычисление якобиана преобразования и получение следующего выражения получаем выражение [math]I = \frac{D(x_1,...x_n)}{D(y_1,...y_n)}I'[/math], где якобиан преобразования имеет вид

[math]\frac{D(x_1,...x_n)}{D(y_1,...y_n)} = (b_1 - a_1)...(b_n - a_n)[/math], [math]I' = {\int\limits_{G}{\ldots \int{g\left( {{y}_{1}},\ldots ,{{y}_{n}} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}} = {\int\limits_{G}{\ldots \int{f\left( {a_1 + (b_1 - a_1)y_1},\ldots ,{a_n + (b_n - a_n)y_n} \right)d{{y}_{1}}\ldots d{{y}_{n}}}}}[/math].

3.Выбор [math]N[/math] равномерно распределенных на отрезке [math][0, 1][/math] последовательностей случайных величин [math]\vec y^{(i)} = ({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})[/math], [math]i = 1,...,N[/math], которые являются координатами случайных точек [math]K_i({y^{(i)}}_{1}, ..., {y^{(i)}}_{n})[/math].

4. Определение того, какие из этих точек принадлежат области интегрирования [math]G[/math]. Пусть таких точек всего [math]m[/math].

5. Вычисление [math]g_{mean} = \mathbb{E}g(K) \approx \tilde{g}_{mean} = \frac{1}{m}\sum^{m}_{i=1}g(K_i) [/math].

6. Получение оценки интеграла [math]I' = \tilde{g}_{mean} * V_G[/math], где [math]V_G[/math] - [math]n[/math] - мерный объем области интегрирования [math]G[/math]. Считаем, что [math]V_G \approx \tilde{V}_G = \frac{m}{N}[/math].

7. Получение оценки искомого интеграла [math]I = V_{par} * \tilde{g}_{mean} * \tilde{V}_G = \frac{V_{par}}{N}\sum^{m}_{i=1}g(K_i)[/math], где [math]V_{par} = (b_1 - a_1)...(b_n - a_n)[/math] - объем параллелепипеда, ограничивающего исходную область интегрирования [math]D[/math].

1.5.3 Геометрический метод

Схема реализации последовательного алгоритма аналогична одномерному случаю.

1.6 Последовательная сложность алгоритма

Для вычисления определенного интеграла методом Монте - Карло требуется осуществить:

1) при использовании метода усреднения подынтегральной функции

  • [math]N[/math] генераций случайных величин;
  • [math]N[/math] вычислений значений функции;
  • [math]N - 1[/math] сложение.

2) при использовании геометрического метода

  • [math]N[/math] генераций случайных величин;
  • [math]N[/math] вычислений значений функции;
  • [math]N[/math] сравнений.

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

1.7 Информационный граф

1.7.1 Метод усреднения подынтегральной функции

Граф алгоритма состоит из трех групп вершин. Вершины на верхнем ярусе используют в качестве входных данных сгенерированные случайные значения и вычисляют в этих точках функцию [math]f[/math]. Вершины второго яруса суммируют полученные значения. На третьем ярусе происходит умножение полученной на предыдущем ярусе суммы на соответствующую константу, формула для вычисления которой представлена ранее при описании алгоритма.

Рисунок 1. Граф алгоритма при использовании метода усреднения подынтегральной функции

1.7.2 Геометрический метод

Граф алгоритма состоит из четырех групп вершин. При данном методе сначала производится ограничение функции [math]f[/math] параллелепипедом и генерация соответствующих случайных точек. Далее полученные значения подаются на вход вершинам верхнего яруса, которые вычисляют значение функции [math]f[/math] в них. Вершины второго яруса производят сравнивание, в зависимости от которого на последующем ярусе вычисляется количество точек, попавших под график функции. Затем производится умножение полученного результата на соответствующую константу, описанную ранее.

Рисунок 2. Граф алгоритма при использовании геометрического метода

1.8 Ресурс параллелизма алгоритма

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

  • 1 ярус вычислений значений функции в полученных ранее случайных точках (N вычислений);
  • 1 ярус сложений (N - 1 сложение);
  • 1 ярус умножений (1 умножение на константу).

Для вычисления интеграла методом Монте - Карло при использовании геометрической интерпретации требуется последовательно выполнить следующие ярусы:

  • 1 ярус вычислений значений функции в полученных ранее случайных точках (N вычислений);
  • 1 ярус сравнений (N сравнений);
  • 1 ярус сложений (N - 1 сложений);
  • 1 ярус умножений (1 умножение на константу).

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

При классификации по высоте ярусно-параллельной формы алгоритм метода Монте - Карло имеет сложность 3 или 4 в зависимости от конкретного выбора одного из двух предлагаемых вариантов. При классификации по ширине его сложность будет [math]O(N)[/math].

1.9 Входные и выходные данные алгоритма

Входные данные: подынтегральная функция [math]f[/math], кратность интеграла [math]p[/math], пределы интегрирования и количество генерируемых случайных точек [math]N[/math].

Выходные данные: число - искомое значение интеграла.

1.10 Свойства алгоритма

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов является линейным (отношение линейной сложности к константной).

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

При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, линейна.

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

Выведем формулу для вычисления значения погрешности:

[math]\Delta I = \Delta(V_{par} * \tilde{g}_{mean} * \tilde{V_G}) = V_{par} * (\Delta g_{mean} * \tilde{V_G} + \Delta V_G * |\tilde{g}_{mean}|)[/math].

По центральной предельной теореме

[math] P \left\{ |g_{mean} - \tilde{g}_{mean}| \le t_{\beta}\sqrt{\frac{\sigma_{1}^2}{m}} \right\} \approx \beta [/math], [math] P \left\{ |V_G - \tilde{V}_G| \le t_{\beta}\sqrt{\frac{\sigma_{2}^2}{N}} \right\} \approx 2\Phi(t_{\beta}) = \beta [/math],

где [math]\beta[/math] - коэффициент доверия, [math]\Phi(t)[/math] - функция стандартного нормального распределения, [math]\sigma_{1}^2 = D[g(K)], \sigma_{2}^2 = D[V_G][/math].

То есть с вероятностью [math]\beta[/math]

[math]\Delta g_{mean} = t_{\beta} \sqrt{\frac{\sigma_{1}^2}{m}}, \Delta V_G = t_{\beta} \sqrt{\frac{\sigma_{2}^2}{N}}[/math].

При этом [math]\tilde{g}_{mean} = \frac{1}{m}\sum^{m}_{i=1}g(K_i), \tilde{V}_{G} = \frac{m}{N} [/math].

Таким образом,

[math]\Delta I = V_{par} * t_{\beta} (\sigma_1 * \frac{\tilde{V}_G}{m} + \sigma_2 * \frac{|\tilde{g}_{mean}|}{N})[/math].

Обычно величины [math]\sigma_1^2[/math] и [math]\sigma_2^2[/math] оказываются неизвестными, тогда их приближенные значения вычисляются по формулам

[math]\tilde{\sigma}_1^2 \approx (\tilde{g}^2)_{mean} - (|\tilde{g}_{mean}|)^2[/math], где [math](\tilde{g}^2)_{mean} = \frac{1}{m}\sum^{m}_{i=1}g^2(K_i)[/math];

[math]\tilde{\sigma}_2^2 \approx \tilde{V}_G (1 - \tilde{V}_G)[/math].

Таким образом, получаем следующее соотношение для погрешности интегрирования методом Монте - Карло:

[math]\Delta I = V_{par} * t_{\beta} (\tilde{\sigma}_1 * \frac{\sqrt{\tilde{V}_G}}{\sqrt{N}} + \tilde{\sigma}_2 * \frac{|\tilde{g}_{mean}|}{\sqrt{N}})[/math].

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

Дуги информационного графа являются полностью локальными.

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

2.1.1 Метод усреднения подынтегральной функции

Код последовательной реализации для случая использования метода усреднения подынтегральной функции при простейшем предположении, что предел интегрирования по каждой из переменных расположен внутри отрезка [0,1], на языке C может выглядеть так: 

double mean_method(long int number, int order, double *limits)
{
    double coord[order], res;
    int n, i, ind;
    res = 0.0;
    for (n = 0; n < number; n++) {
        ind = 0;
        for (i = 0; i < order; i++) {
            coord[i] = (double)(rand())/RAND_MAX;
            if (coord[i] > limits[2 * i + 1] || coord[i] < limits[2 * i]) {
                ind = 1;
                break;
            }
        }
        if (ind == 0) {
            res = res + f(coord, order);
        }
    }                                   
    return res;
}

2.1.2 Геометрический метод

Код последовательной реализации для случая использования геометрического метода при простейшем предположении, что функцию можно ограничить параллелепипедом, каждая из сторон которого является отрезком [0,1], на языке C может выглядеть так: 

double geometric_method(long int number, int order)
{
    double coord[order + 1];
    int score, n, i;
    score = 0;
    for (n = 0; n < number; n++) {
        for (i = 0; i <= order; i++) {
            coord[i] = (double)(rand())/RAND_MAX;        
        }
        
        if (coord[order] <= f(coord, order))
            score++;
    }
    return score;
}

2.2 Локальность данных и вычислений

Локальность вычислений существенно определяется реализацией вычисления функции [math]f[/math]. В предположении, что данная реализация уже имеется, будем считать, что в ней учтены требования эффективного использования памяти.

2.3 Возможные способы и особенности параллельной реализации алгоритма

Самым очевидным вариантом реализации параллельного алгоритма является независимое вычислении функции f в генерируемых случайных точках различными вычислительными узлами. Это может быть реализовано как при помощи технологии threads или OpenMP при запуске на одном узле, так и при использовании межпроцессорных обменов при помощи MPI. Головной узел передает остальным вычислительным узлам реализацию функции f и точки или наборы точек, в которых должна вычисляться функция.

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

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

2.4 Масштабируемость алгоритма и его реализации

Согласно теоретическим оценкам алгоритмы должны быть отлично масштабируемы. Докажем это предположение на практике.

Для этого были созданы реализации описанных методов. В качестве тестового примера рассматривалось интегрирование функции трех переменных [math]\sqrt{x_1} + \sqrt{x_2} - \sqrt{x_3}[/math]. Использование метода Монте - Карло в данном случае нерационально, поскольку оптимальные выбрать другой детерминированный метод, но так как в данной задаче вычисление подынтегральной функции рассматривается как "черный ящик", следовательно сложность вычисления функции не имеет значения.

Выбранный алгоритм использует в качестве метода распараллеливания технологию MPI. Для проведения экспериментов использовалась программа, написанная на языке C, с использованием библиотеки Intel MPI 5.0.1. Тесты проводились на суперкомпьютере Ломоносов[3] в узле test. Для компиляции программы использовалась команда mpicc, которая является надстройкой над имеющимися компиляторами и автоматически подключает заголовочные файлы и библиотеки MPI. Для демонстрации оказалось достаточным провести тесты с количеством процессоров до 128 с масштабированием по степени числа 2. Для количества точек, в которых вычислялась функция, была введена логарифмическая шкала: от 1000 до 1000000000 с шагом умножения на 10. Такая шкала была выбрана для того, чтобы полностью продемонстрировать постепенный рост производительности при увеличении числа процессоров и падение эффективности при небольших количествах точек интегрирования.

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

В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:

  • минимальная эффективность реализации 1,08%;
  • максимальная эффективность реализации 99,8%.


Ниже представлены графики производительности и эффективности реализации алгоритма интегрирования методом усреднения подынтегральной функции в зависимости от изменяемых параметров запуска. Так два рассматриваемых алгоритма интегрирования методом Монте - Карло имеют похожие реализации, то для геометрического метода графики будут аналогичными.

Рисунок 3. Параллельная реализация интегрирования методом Монте - Карло. Изменение производительности в зависимости от числа процессоров и количества генераций случайных точек.
Рисунок 4. Параллельная реализация интегрирования методом Монте - Карло. Изменение эффективности в зависимости от числа процессоров и количества генераций случайных точек.

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

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

2.7 Существующие реализации алгоритма

В рамках библиотеки GNU Scientific Library (или GSL)[4], которая была написана для численных вычислений в прикладной математике и науке, реализован алгоритм интегрирования функций многих переменных методом Монте - Карло[5]. Данная библиотека является свободным программным обеспечением. Распространяется под лицензией GNU GPL 3.

Других реализаций обнаружено не было.

3 Литература

  1. ↑ N. Metropolis, S. Ulam The Monte Carlo Method, — J. Amer. statistical assoc. 1949 44 № 247 335—341.
  2. ↑ Соболь И. М. Метод Монте-Карло. — М.: Наука, 1968. — 64 с. — (Популярные лекции по математике). — 79 000 экз.
  3. ↑ http://users.parallel.ru/wiki/pages/22-config
  4. ↑ https://www.gnu.org/software/gsl/
  5. ↑ https://www.gnu.org/software/gsl/doc/html/montecarlo.html
Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Grachovamc/Интегрирование_методом_Монте_-_Карло&oldid=25006