Участник:Lonalone/Генерация гауссовского вектора методом линейных преобразований: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показаны 123 промежуточные версии этого же участника)
Строка 1: Строка 1:
 
Автор описания: ''Меньших И. М.''
 
Автор описания: ''Меньших И. М.''
−
----
 
 
== Свойства и структура алгоритма ==
 
== Свойства и структура алгоритма ==
  
 
=== Общее описание алгоритма ===
 
=== Общее описание алгоритма ===
−
В статье приведен алгоритм генерации n-мерного случайного вектора, распределенного по нормальному закону, с помощью метода линейных преобразований<ref name="VVV">Балдин К.В., Уткин В.Б. Информационные системы в экономике. - М.:ИТК Дашков и К, 2008. - 395 с.</ref>. Этот метод является одним из наиболее распространенных так называемых корреляционных методов, применяемых в случаях, когда при моделировании непрерывного n-мерного случайного вектора достаточно обеспечить лишь требуемые значения элементов корреляционной матрицы этого вектора (для случая нормального распределения выполнение названного требования означает выполнение достаточного условия полного статистического соответствия теоретического и моделируемого распределений<ref name="VVVVVV"> https://ru.wikipedia.org/wiki/Многомерное_нормальное_распределение</ref>) и вектора математических ожиданий компонент.<br><br>
+
 
−
Идея алгоритма заключается в линейном преобразовании случайного n-мерного вектора Y c независимыми, одинаково распределенными по стандартному нормальному закону компонентами в случайный вектор X с требуемыми корреляционной матрицей и вектором математических ожиданий компонент.
+
Из многомерных распределений особый интерес представляет нормальное. Этому закону подчиняются все результаты воздействия большого числа случайных факторов, среди которых нет превалирующих.
 +
 
 +
В статье приведен алгоритм генерации n-мерного гауссовского случайного вектора с помощью ''метода линейных преобразований''<ref name="V">Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование. Методы Монте-Карло — М.: Академия, 2006. — 368 с.</ref>. Известно, что случае нормально распределенного случайного вектора, ковариационная матрица вместе с математическим ожиданием этого вектора полностью определяют его распределение. Поэтому для полного статистического соответствия моделируемого и теоретического распределения гауссовского вектора достаточно обеспечить требуемые значения указанных параметров<ref name="VV"> https://ru.wikipedia.org/wiki/Многомерное_нормальное_распределение</ref>.
 +
 
 +
Идея алгоритма заключается в линейном преобразовании n-мерного случайного вектора <math>Y</math>, компоненты которого независимы и одинаково распределены по нормальному закону со стандартными параметрами, в случайный вектор <math>X</math> с требуемыми ковариационной матрицей и вектором математических ожиданий.
  
 
=== Математическое описание алгоритма ===
 
=== Математическое описание алгоритма ===
  
−
Даны корреляционная матрица Q и вектор математических ожиданий компонент вектора X:
+
==== Метод линейных преобразований ====
 +
 
 +
Даны ковариационная матрица <math>\Sigma</math> и вектор математических ожиданий <math>M</math>:
 +
:<math>
 +
\Sigma = \|\sigma_{ij}\| = \| \mathbb{E}[(X_{i} - m_{x_{i}})(X_{j} - m_{x_{j}})]\|, \\
 +
M = (m_{x_{1}}, m_{x_{2}}, ..., m_{x_{n}})^T.
 +
</math>
 +
 
 +
Требуется найти такую матрицу <math>B</math>, которая позволяла бы получить искомый вектор <math>X</math> с требуемыми характеристиками в результате линейного преобразования <math>X = BY + M</math>, где <math>Y</math> — n-мерный случайный вектор с независимыми нормально распределенными компонентами со стандартными параметрами.
 +
 
 +
Будем искать матрицу <math>B</math> в виде нижней треугольной матрицы. Перейдем от матричной записи к системе алгебраических уравнений:
 +
 
 +
:<math>
 +
 
 +
\begin{pmatrix} X_{1} \\
 +
X_{2} \\ \vdots \\ X_{n} \end{pmatrix} =
 +
 
 +
\begin{pmatrix} b_{11} & 0 & \cdots & 0
 +
\\b_{21} & b_{22} & \cdots & 0
 +
\\ \vdots & \vdots & \ddots & \vdots
 +
\\ b_{n1} & b_{n2} & \cdots & b_{nn}
 +
\end{pmatrix} \times
 +
 
 +
\begin{pmatrix} Y_{1} \\
 +
Y_{2} \\ \vdots \\ Y_{n} \end{pmatrix} +
 +
 
 +
\begin{pmatrix} m_{x_{1}} \\
 +
m_{x_{2}} \\ \vdots \\ m_{x_{n}} \end{pmatrix} \Rightarrow
 +
</math>
 +
<br>
 +
:<math>
 +
 
 +
\begin{cases}X_{1} - m_{x_{1}} = b_{11}Y_{1}
 +
\\X_{2} - m_{x_{2}} = b_{21}Y_{1} + b_{22}Y_{2}
 +
\\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots
 +
\\X_{n} - m_{x_{n}} = b_{n1}Y_{1} + b_{n2}Y_{2} + \cdots + b_{nn}Y_{n}
 +
\end{cases}
 +
 
 +
</math>
 +
 
 +
 
 +
Поскольку компоненты вектора <math>Y</math> независимы и имеют стандартные параметры, справедливо выражение:
 +
 
 
:<math>
 
:<math>
−
\begin{align}
+
 
−
Q = \|q_{ij}\| = \| M[(X_{i} - m_{X_{i}})(X_{j} - m_{X_{j}})]\|, \\
+
\mathbb{E}[Y_{i}Y_{j}] =
−
M = (m_{X_{1}}, m_{X_{2}}, ..., m_{X_{n}})^T.
+
\left\{\begin{matrix}
−
\end{align}
+
1, &i = j, \\
 +
0, &i \not= j.
 +
\end{matrix}\right.
 +
 
 
</math>
 
</math>
 +
 +
Почленно перемножив сами на себя и между собой соответственно левые и правые части уравнений системы и взяв от результатов перемножения математическое ожидание, получаем систему уравнений вида:
 +
 +
:<math>
 +
 +
\begin{cases}
 +
\mathbb{E}[(X_{1} - m_{x_{1}})(X_{1} - m_{x_{1}})] = \mathbb{E}[b_{11}Y_{1}b_{11}Y_{1}],
 +
\\ \mathbb{E}[(X_{1} - m_{x_{1}})(X_{2} - m_{x_{2}})] = \mathbb{E}[(b_{21}Y_{1} + b_{22}Y_{2})b_{11}Y_{1}],
 +
\\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots
 +
\end{cases}
 +
 +
</math>
 +
 +
Таким образом, в левых частях полученной системы уравнений располагаются элементы заданной ковариационной матрицы <math>\Sigma</math>, а в правых — элементы искомой матрицы <math>B</math>. Последовательно решая эту систему, получаем формулы для расчета элементов <math>\sigma_{ij}</math>:
 +
 +
:<math>
 +
b_{11}=\sqrt{\sigma_{11}}; b_{21}=\frac{\sigma_{12}}{\sqrt{\sigma_{11}}}; b_{22} = \sqrt{\sigma_{22} - \frac{\sigma_{12}}{\sigma_{11}}}, \cdots
 +
 +
</math>
 +
 +
Рекуррентная формула для расчета любого элемента матрицы преобразования <math>B</math> имеет вид:
 +
 +
:<math>
 +
 +
b_{ij} = \frac {\sigma_{ij} - \sum_{k=1}^{j-1} b_{ik} b_{jk}} {\sqrt{\sigma_{jj} - \sum_{k=1}^{j-1} b_{jk}^2}}, \quad 1 \leqslant j \leqslant i \leqslant n
 +
 +
</math>
 +
 +
(суммы с верхним нулевым пределом считаются равными нулю)
 +
 +
Таким образом, матрица <math>B</math> получается с помощью [[Разложение_Холецкого_(метод_квадратного_корня)|разложения Холецкого]] матрицы <math>\Sigma</math>.
 +
 +
==== Генерация случайного вектора Y ====
 +
 +
Пусть имеется генератор псевдослучайных чисел (ГПСЧ, PRNG), с помощью которого можно получить реализацию случайной величины <math>u \sim U(0,1) </math>. Описанный выше случайный вектор <math>Y=(y_1, \dots, y_n)</math> с независимыми компонентами составим из <math>n</math> реализаций случайной величины <math>\eta \sim N(0,1)</math>. Каждую такую реализацию <math>y_i</math>, в свою очередь, получим с помощью приближения по ЦПТ <math>\kappa</math> случайными величинами, распределенными равномерно на отрезке [0,1]:
 +
 +
:<math>
 +
 +
y_i = \sqrt{\frac{12}{\kappa}} \left(\sum_{j=1}^{\kappa}{u_j^i} - \frac{\kappa}{2}\right), \quad i = \overline{1, n} \\
 +
</math>
 +
где <math>u_j^i - </math> реализации случайной величины <math>u</math>, а <math>\kappa - </math> параметр, обеспечивающий качество приближения.
 +
 +
Как правило, <math>\kappa</math> берут равным <math>12</math> и считают, что для подавляющего числа практических задач обеспечивается должная точность вычислений<ref name = "VVV">Балдин К.В., Уткин В.Б. Информационные системы в экономике. — М.:Дашков и Кo, 2008. — 395 с.</ref>.
  
 
=== Вычислительное ядро алгоритма ===
 
=== Вычислительное ядро алгоритма ===
 +
 +
Вычислительное ядро алгоритма можно составить из множественных (всего  <math>\frac{n(n−1)}{2}</math>) вычислений следующих выражений, включающих скалярные произведения строк матрицы:
 +
 +
<math>\sigma_{ij} - \sum_{k=1}^{j-1} b_{ik} b_{jk}, \quad 1 \leqslant j \leqslant i \leqslant n</math>
  
 
=== Макроструктура алгоритма ===
 
=== Макроструктура алгоритма ===
 +
 +
* Заполнение матрицы <math>B</math>
 +
* Генерация вектора <math>Y</math>
 +
* Вычисление вектора <math>X = BY + M</math>
  
 
=== Схема реализации последовательного алгоритма ===
 
=== Схема реализации последовательного алгоритма ===
 +
 +
1) Заполнение матрицы <math>B</math>.
 +
 +
: (a) <math>b_{11}= \sqrt{\sigma_{11}}</math>
 +
 +
: (b) <math>b_{k1}= \frac{\sigma_{k1}}{l_{11}}</math>, <math>\quad k = \overline{2,n}</math>
 +
 +
: Следующие два пункта выполняются циклически, друг за другом для <math>i = \overline{2,n}</math>.
 +
 +
: (c) <math>b_{ii} = \sqrt{\sigma_{ii} - \sum_{p = 1}^{i - 1} b_{ip}^2}</math>
 +
 +
: (d) <math>b_{ji} = \frac{\sigma_{ji} - \sum_{p = 1}^{i - 1} b_{ip} b_{jp}} {l_{ii}}</math>, <math>\quad i \neq n, j = \overline{i+1,n}</math> <br>
 +
 +
2) Генерация вектора <math>Y</math>.
 +
 +
: Следующие два пункта выполняются циклически, друг за другом для <math>i = \overline{1,n}</math>.
 +
 +
: (a) <math>u_j^i \leftarrow PRNG,</math> <math>\quad j = \overline{1,k}</math>
 +
 +
: (b) <math>y_i = \sqrt{\frac{12}{\kappa}} \left(\sum_{j=1}^{\kappa}{u_j^i} - \frac{\kappa}{2}\right)</math>
 +
 +
3) Вычисление вектора <math>X</math>.
 +
 +
: <math>x_i = m_{x_i} + \sum_{k = 1}^{i} b_{ik} y_{k}, \quad i = \overline{1,n}</math>.
  
 
=== Последовательная сложность алгоритма ===
 
=== Последовательная сложность алгоритма ===
 +
 +
Для генерации гауссовского случайного вектора порядка <math>n</math> в последовательном варианте требуется:
 +
* <math> n(n+1) + \frac {n(n-1)(n-2)}{6} </math> (умножений)
 +
* <math> \frac {n(n+1)}{2} + \frac {n(n-1)(n-2)}{6} </math> (вычитаний)
 +
* <math> \frac {n(n-1+2\kappa)}{2} </math> (сложений)
 +
* <math> \frac {n(n+3)}{2} </math> (делений)
 +
* <math> 2n </math> (вычислений квадратного корня)
 +
 +
<math> \frac{n}{6}\left(2n^2+9n+31+6\kappa\right) </math> - итого операций.
 +
 +
Таким образом, последовательная сложность алгоритма равна <math>O(n^3)</math>.
  
 
=== Информационный граф ===
 
=== Информационный граф ===
 +
 +
==== Заполнение матрицы B====
 +
 +
Граф этапа 1 макроструктуры алгоритма ([[Разложение_Холецкого_(метод_квадратного_корня)|''разложения Холецкого'']]) состоит из трёх групп вершин, расположенных в целочисленных узлах трёх областей разной размерности.
 +
 +
'''Первая''' группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию SQRT.
 +
Единственная координата каждой из вершин <math>i</math> меняется в диапазоне от <math>1</math> до <math>n</math>, принимая все целочисленные значения.
 +
 +
Аргумент этой функции
 +
 +
* при <math>i = 1</math> — элемент ''входных данных'' <math>\sigma_{11}</math>;
 +
* при <math>i > 1</math> — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами <math>i - 1</math>, <math>i</math>, <math>i - 1</math>.
 +
Результат срабатывания операции является ''выходным данным'' <math>b_{ii}</math>.
 +
 +
'''Вторая''' группа вершин расположена в двумерной области, соответствующая ей операция <math>a / b</math>.
 +
Естественно введённые координаты области таковы:
 +
* <math>i</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
 +
* <math>j</math> — меняется в диапазоне от <math>i+1</math> до <math>n</math>, принимая все целочисленные значения.
 +
 +
Аргументы операции следующие:
 +
*<math>a</math>:
 +
** при <math>i = 1</math> — элемент ''входных данных'' <math>\sigma_{j1}</math>;
 +
** при <math>i > 1</math> — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами <math>i - 1, j, i - 1</math>;
 +
* <math>b</math> — результат срабатывания операции, соответствующей вершине из первой группы, с координатой <math>i</math>.
 +
 +
Результат срабатывания операции является ''выходным данным'' <math>b_{ji}</math>.
 +
 +
'''Третья''' группа вершин расположена в трёхмерной области, соответствующая ей операция  <math>a - b * c</math>.
 +
Естественно введённые координаты области таковы:
 +
* <math>i</math> — меняется в диапазоне от <math>2</math> до <math>n</math>, принимая все целочисленные значения;
 +
* <math>j</math> — меняется в диапазоне от <math>i</math> до <math>n</math>, принимая все целочисленные значения;
 +
* <math>p</math> — меняется в диапазоне  от <math>1</math> до <math>i - 1</math>, принимая все целочисленные значения.
 +
 +
Аргументы операции следующие:
 +
*<math>a</math>:
 +
** при <math>p = 1</math> элемент входных данных <math>\sigma_{ji}</math>;
 +
** при <math>p > 1</math> — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами <math>i, j, p - 1</math>;
 +
*<math>b</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>p, i</math>;
 +
*<math>c</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>p, j</math>;
 +
 +
Результат срабатывания операции является ''промежуточным данным'' этапа.
 +
 +
Описанный граф можно посмотреть на рис.1, выполненном для случая <math>n = 4</math>.
 +
 +
[[file:Cholesky full_in_out.png|thumb|center|1400px|Рисунок 1. Граф этапа 1 макроструктуры алгоритма. SQ — вычисление квадратного корня, F — операция a-bc, Div — деление, In — входные данные, Out — выходные.]]
 +
 +
====Генерация вектора Y====
 +
 +
Граф этапа 2 макроструктуры алгоритма состоит из двух групп вершин, расположенных в целочисленных узлах двумерной области.
 +
 +
Естественно введённые координаты области таковы:
 +
* <math>i</math> — меняется в диапазоне от <math>1</math> до <math>n</math>, принимая все целочисленные значения;
 +
* <math>j</math> — меняется в диапазоне от <math>1</math> до <math>\kappa</math>, принимая все целочисленные значения.
 +
 +
'''Первая''' группа вершин. Соответствующая ей операция вычисляет функцию <math>S = a+b</math>.
 +
 +
Аргументы операции следующие:
 +
*<math>a:</math>
 +
** при <math>j=1</math> — константа <math>0</math>;
 +
** при <math>j>1</math> — результат срабатывания операции, соответствующей вершине из первой группы, с координатами <math>i,j-1</math>;
 +
*<math>b</math> — элемент ''входных данных'' — значение, поступающее с <math>PRNG</math>.
 +
 +
Результат срабатывания операции является ''промежуточным данным'' этапа.
 +
 +
'''Вторая''' группа вершин. Соответствующая ей операция вычисляет функцию <math>F = \left(a-\frac{\kappa}{2}\right)\sqrt{\frac{12}{\kappa}}</math>.
 +
 +
Аргументы операции следующие:
 +
*<math>a</math> —  результат срабатывания операции, соответствующей вершине из первой группы, с координатами <math>i,\kappa</math>.
 +
 +
Результат срабатывания операции является ''выходным данным'' этапа <math>y_i</math>.
 +
 +
Описанный граф можно посмотреть на рис.2, выполненном для случая <math>n = 2, \kappa = 2</math>.
 +
 +
[[Файл:LonGraf2.png|550px|thumb|center|Рисунок 2. Граф этапа 2 макроструктуры алгоритма]]
 +
 +
====Вычисление вектора X====
 +
 +
Граф этапа 3 макроструктуры алгоритма (''умножение нижней треугольной матрицы на вектор и сложение результата с другим вектором'') состоит из одной группы вершин, расположенной в целочисленных узлах двумерной области, соответствующая ей операция <math>a+bc</math>.
 +
 +
Естественно введённые координаты области таковы:
 +
* <math>i</math> — меняется в диапазоне от <math>1</math> до <math>n</math>, принимая все целочисленные значения;
 +
* <math>j</math> — меняется в диапазоне от <math>1</math> до <math>n</math>, принимая все целочисленные значения.
 +
 +
Аргументы операции следующие:
 +
*<math>a</math>:
 +
** при <math>j = 1</math> — компонента <math>m_{x_i}</math> вектора <math>M</math>;
 +
** при <math>j > 1</math> — результат срабатывания операции, соответствующей вершине с координатами <math>i, j-1</math>;
 +
*<math>b</math> — элемент ''входных данных'' <math>b_{ij}</math>;
 +
*<math>c</math> — элемент ''входных данных'' <math>y_{j}</math>;
 +
 +
Результат срабатывания операции является:
 +
* при <math>j < i</math> — ''промежуточным данным'' этапа;
 +
* при <math>j = i</math> — ''выходным данным'' алгоритма <math>x_i</math>.
 +
 +
Описанный граф можно посмотреть на рис.3, выполненном для случая <math> n = 4</math>.
 +
 +
[[Файл:LonGraf3.png|650px|thumb|center|Рисунок 3. Граф этапа 3 макроструктуры алгоритма,  F — операция a+bc]]
  
 
=== Ресурс параллелизма алгоритма ===
 
=== Ресурс параллелизма алгоритма ===
 +
 +
Для разложения Холецкого <math>(n \times n)</math> матрицы <math>\Sigma</math> в параллельном варианте требуется последовательно выполнить следующие ярусы:
 +
* <math>n</math> ярусов с вычислением квадратного корня (единичные вычисления в каждом из ярусов),
 +
* <math>n - 1</math> ярус делений (в каждом из ярусов — линейное количество операций, в зависимости от яруса — от <math>1</math> до <math>n - 1</math>),
 +
* по <math>n - 1</math> ярусов умножений и сложений/вычитаний (в каждом из ярусов — квадратичное количество операций, от <math>1</math> до <math>\frac{n^2 - n}{2}</math>).
 +
 +
Для генерации <math>(n)</math> вектора <math>Y</math> в параллельном варианте требуется последовательно выполнить следующие ярусы:
 +
* <math>\kappa</math> ярусов сложений (<math>n</math> операций в каждом из ярусов),
 +
* <math>1</math> ярус вычитаний, делений, умножений, вычислений квадратного корня (<math>n, 2n, n, n</math> операций соответственно),
 +
 +
Для вычисления <math>(n)</math> вектора <math>X</math> в параллельном варианте требуется последовательно выполнить следующие ярусы:
 +
* <math>n</math> ярусов умножений и сложений (в каждом из ярусов — линейное количество операций, в зависимости от яруса — от <math>1</math> до <math>n</math>)
 +
 +
Учитывая <math>n \gg \kappa</math>, можно заключить: при классификации по высоте ЯПФ алгоритм имеет сложность <math>O(n)</math>, при классификации по ширине — <math>O(n^2)</math>.
  
 
=== Входные и выходные данные алгоритма ===
 
=== Входные и выходные данные алгоритма ===
 +
 +
* Входные данные: вещественная всюду плотная положительно определенная симметрическая  <math>(n \times n)</math> матрица <math>\Sigma</math> и вещественный <math>(n)</math> вектор <math>M</math>;
 +
* Выходные данные: вещественный <math>(n)</math> вектор <math>X</math>.
  
 
=== Свойства алгоритма ===
 
=== Свойства алгоритма ===
Строка 36: Строка 284:
  
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
 +
 +
<syntaxhighlight lang = c>
 +
 +
// compute matrix B (core)
 +
for (int j = 0; j < n; j++){
 +
    double sum_down = matrix_cov[j][j];
 +
    for (int k = 0; k < j; k++){
 +
        double cur = matrix_cov[j][k];
 +
        sum_down -= cur * cur;
 +
    }
 +
    matrix_cov[j][j] = sqrt(sum_down);
 +
    for (int i = j + 1; i < n; i++){
 +
        double sum_up = matrix_cov[i][j];
 +
        for (int k = 0; k < j; k++){
 +
            sum_up -= matrix_cov[i][k] * matrix_cov[j][k];
 +
        }
 +
        matrix_cov[i][j] = sum_up / matrix_cov[j][j];
 +
    }
 +
}
 +
// generate vector Y
 +
for (int i = 0; i < n; i++){
 +
    double sum = generate_uniform(n);
 +
    for (int j = 0; j < CLT_NUM - 1; j++){
 +
        sum += generate_uniform(n);
 +
    }
 +
    vec_Y[i] = (sum - CLT_NUM / 2) * sqrt(12.0 / CLT_NUM);
 +
}
 +
// linear transformation
 +
for (int i = 0; i < n; i++){ 
 +
    vec_X[i] = vec_exp[i];
 +
    for (int k = 0; k <= i; k++){
 +
        vec_X[i] += matrix_cov[i][k] * vec_Y[k];
 +
    }
 +
}
 +
 +
</syntaxhighlight>
  
 
=== Локальность данных и вычислений ===
 
=== Локальность данных и вычислений ===
  
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
 +
В данной реализации алгоритма весь ресурс параллелизма использован в следующих этапах:
 +
* вычисление элементов матрицы <math>B</math> в рамках столбца (под диагональю),
 +
* генерация вектора <math>Y</math>,
 +
* вычисление вектора <math>X</math> (умножение матрицы на вектор, сложение векторов).
  
 
=== Масштабируемость алгоритма и его реализации ===
 
=== Масштабируемость алгоритма и его реализации ===
 +
 +
==== Масштабируемость алгоритма ====
 +
 +
==== Масштабируемость реализации алгоритма ====
 +
Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета]. Алгоритм был реализован на языке C++ с использованием Intel MPI Library и запущен на кластере regular4.
 +
 +
Параметры компиляции и запуска:
 +
module add slurm
 +
module add impi/5.0.1
 +
mpicxx -std=c++0x task.c -o _scratch/task
 +
sbatch -n256 impi ./task <matrix_size>
 +
 +
Параметры запуска:
 +
* размер матрицы [1024 : 4864] с шагом 256,
 +
* число ядер [2 : 256] с шагом по степеням 2.
 +
<br>
 +
 +
[[file:Ops_2.png|thumb|center|700px|Рисунок 4. Изменение производительности в зависимости от числа ядер и размера матрицы.]]
 +
 +
[[file:Eff_2.png|thumb|center|700px|Рисунок 5. Изменение эффективности реализации в зависимости от числа ядер и размера матрицы.]]
 +
 +
[[file:Eff_2_up.png|thumb|center|700px|Рисунок 6. Изменение эффективности реализации, вид сверху.]]
 +
 +
Оценка масштабируемости:
 +
 +
* '''По размеру задачи''' — при увеличении размерности матрицы эффективность реализации постепенно увеличивается, при этом тем быстрее, чем большее число ядер задействовано.
 +
* '''По числу ядер''' — при увеличении числа mpi-процессов эффективность реализации снижается. Это объясняется ростом накладных расходов на организацию параллельного выполнения. При этом с ростом размерности матрицы эффективность снижается с такой же скоростью, но уже при больших значениях числа mpi-процессов.
 +
* '''По двум направлениям''' — эффективность реализации возрастает (на рассматриваемой области значений параметров), но крайне медленно.
 +
 +
[https://github.com/mensigo/skipod/blob/master/mpi_gauss_vec_gen.c Параллельная реализация на C++]
  
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
=== Динамические характеристики и эффективность реализации алгоритма ===
Строка 48: Строка 367:
  
 
=== Существующие реализации алгоритма ===
 
=== Существующие реализации алгоритма ===
 +
 +
MS Excel и распространённые статистические пакеты (напр., SPSS, Statistica) позволяют моделировать только одномерные статистические распределения. Имеется возможность "собрать" многомерное распределение из нескольких одномерных при условии, что компоненты независимы. Однако если необходимо исследовать данные с зависящими друг от друга переменными, то приходится писать собственную программу. При этом удобно воспользоваться готовой реализацией разложения Холецкого (точечного метода), представленной в таких прикладных пакетах, как LINPACK, LAPACK, SCALAPACK.
  
 
== Литература ==
 
== Литература ==

Текущая версия на 09:54, 8 декабря 2018

Автор описания: Меньших И. М.

Содержание

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

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

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

В статье приведен алгоритм генерации n-мерного гауссовского случайного вектора с помощью метода линейных преобразований[1]. Известно, что случае нормально распределенного случайного вектора, ковариационная матрица вместе с математическим ожиданием этого вектора полностью определяют его распределение. Поэтому для полного статистического соответствия моделируемого и теоретического распределения гауссовского вектора достаточно обеспечить требуемые значения указанных параметров[2].

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

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

1.2.1 Метод линейных преобразований

Даны ковариационная матрица [math]\Sigma[/math] и вектор математических ожиданий [math]M[/math]:

[math] \Sigma = \|\sigma_{ij}\| = \| \mathbb{E}[(X_{i} - m_{x_{i}})(X_{j} - m_{x_{j}})]\|, \\ M = (m_{x_{1}}, m_{x_{2}}, ..., m_{x_{n}})^T. [/math]

Требуется найти такую матрицу [math]B[/math], которая позволяла бы получить искомый вектор [math]X[/math] с требуемыми характеристиками в результате линейного преобразования [math]X = BY + M[/math], где [math]Y[/math] — n-мерный случайный вектор с независимыми нормально распределенными компонентами со стандартными параметрами.

Будем искать матрицу [math]B[/math] в виде нижней треугольной матрицы. Перейдем от матричной записи к системе алгебраических уравнений:

[math] \begin{pmatrix} X_{1} \\ X_{2} \\ \vdots \\ X_{n} \end{pmatrix} = \begin{pmatrix} b_{11} & 0 & \cdots & 0 \\b_{21} & b_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ b_{n1} & b_{n2} & \cdots & b_{nn} \end{pmatrix} \times \begin{pmatrix} Y_{1} \\ Y_{2} \\ \vdots \\ Y_{n} \end{pmatrix} + \begin{pmatrix} m_{x_{1}} \\ m_{x_{2}} \\ \vdots \\ m_{x_{n}} \end{pmatrix} \Rightarrow [/math]


[math] \begin{cases}X_{1} - m_{x_{1}} = b_{11}Y_{1} \\X_{2} - m_{x_{2}} = b_{21}Y_{1} + b_{22}Y_{2} \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \\X_{n} - m_{x_{n}} = b_{n1}Y_{1} + b_{n2}Y_{2} + \cdots + b_{nn}Y_{n} \end{cases} [/math]


Поскольку компоненты вектора [math]Y[/math] независимы и имеют стандартные параметры, справедливо выражение:

[math] \mathbb{E}[Y_{i}Y_{j}] = \left\{\begin{matrix} 1, &i = j, \\ 0, &i \not= j. \end{matrix}\right. [/math]

Почленно перемножив сами на себя и между собой соответственно левые и правые части уравнений системы и взяв от результатов перемножения математическое ожидание, получаем систему уравнений вида:

[math] \begin{cases} \mathbb{E}[(X_{1} - m_{x_{1}})(X_{1} - m_{x_{1}})] = \mathbb{E}[b_{11}Y_{1}b_{11}Y_{1}], \\ \mathbb{E}[(X_{1} - m_{x_{1}})(X_{2} - m_{x_{2}})] = \mathbb{E}[(b_{21}Y_{1} + b_{22}Y_{2})b_{11}Y_{1}], \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \end{cases} [/math]

Таким образом, в левых частях полученной системы уравнений располагаются элементы заданной ковариационной матрицы [math]\Sigma[/math], а в правых — элементы искомой матрицы [math]B[/math]. Последовательно решая эту систему, получаем формулы для расчета элементов [math]\sigma_{ij}[/math]:

[math] b_{11}=\sqrt{\sigma_{11}}; b_{21}=\frac{\sigma_{12}}{\sqrt{\sigma_{11}}}; b_{22} = \sqrt{\sigma_{22} - \frac{\sigma_{12}}{\sigma_{11}}}, \cdots [/math]

Рекуррентная формула для расчета любого элемента матрицы преобразования [math]B[/math] имеет вид:

[math] b_{ij} = \frac {\sigma_{ij} - \sum_{k=1}^{j-1} b_{ik} b_{jk}} {\sqrt{\sigma_{jj} - \sum_{k=1}^{j-1} b_{jk}^2}}, \quad 1 \leqslant j \leqslant i \leqslant n [/math]

(суммы с верхним нулевым пределом считаются равными нулю)

Таким образом, матрица [math]B[/math] получается с помощью разложения Холецкого матрицы [math]\Sigma[/math].

1.2.2 Генерация случайного вектора Y

Пусть имеется генератор псевдослучайных чисел (ГПСЧ, PRNG), с помощью которого можно получить реализацию случайной величины [math]u \sim U(0,1) [/math]. Описанный выше случайный вектор [math]Y=(y_1, \dots, y_n)[/math] с независимыми компонентами составим из [math]n[/math] реализаций случайной величины [math]\eta \sim N(0,1)[/math]. Каждую такую реализацию [math]y_i[/math], в свою очередь, получим с помощью приближения по ЦПТ [math]\kappa[/math] случайными величинами, распределенными равномерно на отрезке [0,1]:

[math] y_i = \sqrt{\frac{12}{\kappa}} \left(\sum_{j=1}^{\kappa}{u_j^i} - \frac{\kappa}{2}\right), \quad i = \overline{1, n} \\ [/math]

где [math]u_j^i - [/math] реализации случайной величины [math]u[/math], а [math]\kappa - [/math] параметр, обеспечивающий качество приближения.

Как правило, [math]\kappa[/math] берут равным [math]12[/math] и считают, что для подавляющего числа практических задач обеспечивается должная точность вычислений[3].

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

Вычислительное ядро алгоритма можно составить из множественных (всего [math]\frac{n(n−1)}{2}[/math]) вычислений следующих выражений, включающих скалярные произведения строк матрицы:

[math]\sigma_{ij} - \sum_{k=1}^{j-1} b_{ik} b_{jk}, \quad 1 \leqslant j \leqslant i \leqslant n[/math]

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

  • Заполнение матрицы [math]B[/math]
  • Генерация вектора [math]Y[/math]
  • Вычисление вектора [math]X = BY + M[/math]

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

1) Заполнение матрицы [math]B[/math].

(a) [math]b_{11}= \sqrt{\sigma_{11}}[/math]
(b) [math]b_{k1}= \frac{\sigma_{k1}}{l_{11}}[/math], [math]\quad k = \overline{2,n}[/math]
Следующие два пункта выполняются циклически, друг за другом для [math]i = \overline{2,n}[/math].
(c) [math]b_{ii} = \sqrt{\sigma_{ii} - \sum_{p = 1}^{i - 1} b_{ip}^2}[/math]
(d) [math]b_{ji} = \frac{\sigma_{ji} - \sum_{p = 1}^{i - 1} b_{ip} b_{jp}} {l_{ii}}[/math], [math]\quad i \neq n, j = \overline{i+1,n}[/math]

2) Генерация вектора [math]Y[/math].

Следующие два пункта выполняются циклически, друг за другом для [math]i = \overline{1,n}[/math].
(a) [math]u_j^i \leftarrow PRNG,[/math] [math]\quad j = \overline{1,k}[/math]
(b) [math]y_i = \sqrt{\frac{12}{\kappa}} \left(\sum_{j=1}^{\kappa}{u_j^i} - \frac{\kappa}{2}\right)[/math]

3) Вычисление вектора [math]X[/math].

[math]x_i = m_{x_i} + \sum_{k = 1}^{i} b_{ik} y_{k}, \quad i = \overline{1,n}[/math].

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

Для генерации гауссовского случайного вектора порядка [math]n[/math] в последовательном варианте требуется:

  • [math] n(n+1) + \frac {n(n-1)(n-2)}{6} [/math] (умножений)
  • [math] \frac {n(n+1)}{2} + \frac {n(n-1)(n-2)}{6} [/math] (вычитаний)
  • [math] \frac {n(n-1+2\kappa)}{2} [/math] (сложений)
  • [math] \frac {n(n+3)}{2} [/math] (делений)
  • [math] 2n [/math] (вычислений квадратного корня)

[math] \frac{n}{6}\left(2n^2+9n+31+6\kappa\right) [/math] - итого операций.

Таким образом, последовательная сложность алгоритма равна [math]O(n^3)[/math].

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

1.7.1 Заполнение матрицы B

Граф этапа 1 макроструктуры алгоритма (разложения Холецкого) состоит из трёх групп вершин, расположенных в целочисленных узлах трёх областей разной размерности.

Первая группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию SQRT. Единственная координата каждой из вершин [math]i[/math] меняется в диапазоне от [math]1[/math] до [math]n[/math], принимая все целочисленные значения.

Аргумент этой функции

  • при [math]i = 1[/math] — элемент входных данных [math]\sigma_{11}[/math];
  • при [math]i \gt 1[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]i - 1[/math], [math]i[/math], [math]i - 1[/math].

Результат срабатывания операции является выходным данным [math]b_{ii}[/math].

Вторая группа вершин расположена в двумерной области, соответствующая ей операция [math]a / b[/math]. Естественно введённые координаты области таковы:

  • [math]i[/math] — меняется в диапазоне от [math]1[/math] до [math]n-1[/math], принимая все целочисленные значения;
  • [math]j[/math] — меняется в диапазоне от [math]i+1[/math] до [math]n[/math], принимая все целочисленные значения.

Аргументы операции следующие:

  • [math]a[/math]:
    • при [math]i = 1[/math] — элемент входных данных [math]\sigma_{j1}[/math];
    • при [math]i \gt 1[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]i - 1, j, i - 1[/math];
  • [math]b[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатой [math]i[/math].

Результат срабатывания операции является выходным данным [math]b_{ji}[/math].

Третья группа вершин расположена в трёхмерной области, соответствующая ей операция [math]a - b * c[/math]. Естественно введённые координаты области таковы:

  • [math]i[/math] — меняется в диапазоне от [math]2[/math] до [math]n[/math], принимая все целочисленные значения;
  • [math]j[/math] — меняется в диапазоне от [math]i[/math] до [math]n[/math], принимая все целочисленные значения;
  • [math]p[/math] — меняется в диапазоне от [math]1[/math] до [math]i - 1[/math], принимая все целочисленные значения.

Аргументы операции следующие:

  • [math]a[/math]:
    • при [math]p = 1[/math] элемент входных данных [math]\sigma_{ji}[/math];
    • при [math]p \gt 1[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]i, j, p - 1[/math];
  • [math]b[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]p, i[/math];
  • [math]c[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]p, j[/math];

Результат срабатывания операции является промежуточным данным этапа.

Описанный граф можно посмотреть на рис.1, выполненном для случая [math]n = 4[/math].

Рисунок 1. Граф этапа 1 макроструктуры алгоритма. SQ — вычисление квадратного корня, F — операция a-bc, Div — деление, In — входные данные, Out — выходные.

1.7.2 Генерация вектора Y

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

Естественно введённые координаты области таковы:

  • [math]i[/math] — меняется в диапазоне от [math]1[/math] до [math]n[/math], принимая все целочисленные значения;
  • [math]j[/math] — меняется в диапазоне от [math]1[/math] до [math]\kappa[/math], принимая все целочисленные значения.

Первая группа вершин. Соответствующая ей операция вычисляет функцию [math]S = a+b[/math].

Аргументы операции следующие:

  • [math]a:[/math]
    • при [math]j=1[/math] — константа [math]0[/math];
    • при [math]j\gt 1[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатами [math]i,j-1[/math];
  • [math]b[/math] — элемент входных данных — значение, поступающее с [math]PRNG[/math].

Результат срабатывания операции является промежуточным данным этапа.

Вторая группа вершин. Соответствующая ей операция вычисляет функцию [math]F = \left(a-\frac{\kappa}{2}\right)\sqrt{\frac{12}{\kappa}}[/math].

Аргументы операции следующие:

  • [math]a[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатами [math]i,\kappa[/math].

Результат срабатывания операции является выходным данным этапа [math]y_i[/math].

Описанный граф можно посмотреть на рис.2, выполненном для случая [math]n = 2, \kappa = 2[/math].

Рисунок 2. Граф этапа 2 макроструктуры алгоритма

1.7.3 Вычисление вектора X

Граф этапа 3 макроструктуры алгоритма (умножение нижней треугольной матрицы на вектор и сложение результата с другим вектором) состоит из одной группы вершин, расположенной в целочисленных узлах двумерной области, соответствующая ей операция [math]a+bc[/math].

Естественно введённые координаты области таковы:

  • [math]i[/math] — меняется в диапазоне от [math]1[/math] до [math]n[/math], принимая все целочисленные значения;
  • [math]j[/math] — меняется в диапазоне от [math]1[/math] до [math]n[/math], принимая все целочисленные значения.

Аргументы операции следующие:

  • [math]a[/math]:
    • при [math]j = 1[/math] — компонента [math]m_{x_i}[/math] вектора [math]M[/math];
    • при [math]j \gt 1[/math] — результат срабатывания операции, соответствующей вершине с координатами [math]i, j-1[/math];
  • [math]b[/math] — элемент входных данных [math]b_{ij}[/math];
  • [math]c[/math] — элемент входных данных [math]y_{j}[/math];

Результат срабатывания операции является:

  • при [math]j \lt i[/math] — промежуточным данным этапа;
  • при [math]j = i[/math] — выходным данным алгоритма [math]x_i[/math].

Описанный граф можно посмотреть на рис.3, выполненном для случая [math] n = 4[/math].

Рисунок 3. Граф этапа 3 макроструктуры алгоритма, F — операция a+bc

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

Для разложения Холецкого [math](n \times n)[/math] матрицы [math]\Sigma[/math] в параллельном варианте требуется последовательно выполнить следующие ярусы:

  • [math]n[/math] ярусов с вычислением квадратного корня (единичные вычисления в каждом из ярусов),
  • [math]n - 1[/math] ярус делений (в каждом из ярусов — линейное количество операций, в зависимости от яруса — от [math]1[/math] до [math]n - 1[/math]),
  • по [math]n - 1[/math] ярусов умножений и сложений/вычитаний (в каждом из ярусов — квадратичное количество операций, от [math]1[/math] до [math]\frac{n^2 - n}{2}[/math]).

Для генерации [math](n)[/math] вектора [math]Y[/math] в параллельном варианте требуется последовательно выполнить следующие ярусы:

  • [math]\kappa[/math] ярусов сложений ([math]n[/math] операций в каждом из ярусов),
  • [math]1[/math] ярус вычитаний, делений, умножений, вычислений квадратного корня ([math]n, 2n, n, n[/math] операций соответственно),

Для вычисления [math](n)[/math] вектора [math]X[/math] в параллельном варианте требуется последовательно выполнить следующие ярусы:

  • [math]n[/math] ярусов умножений и сложений (в каждом из ярусов — линейное количество операций, в зависимости от яруса — от [math]1[/math] до [math]n[/math])

Учитывая [math]n \gg \kappa[/math], можно заключить: при классификации по высоте ЯПФ алгоритм имеет сложность [math]O(n)[/math], при классификации по ширине — [math]O(n^2)[/math].

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

  • Входные данные: вещественная всюду плотная положительно определенная симметрическая [math](n \times n)[/math] матрица [math]\Sigma[/math] и вещественный [math](n)[/math] вектор [math]M[/math];
  • Выходные данные: вещественный [math](n)[/math] вектор [math]X[/math].

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

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

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

// compute matrix B (core)
for (int j = 0; j < n; j++){
    double sum_down = matrix_cov[j][j];
    for (int k = 0; k < j; k++){
        double cur = matrix_cov[j][k];
        sum_down -= cur * cur;
    }
    matrix_cov[j][j] = sqrt(sum_down);
    for (int i = j + 1; i < n; i++){
        double sum_up = matrix_cov[i][j];
        for (int k = 0; k < j; k++){
            sum_up -= matrix_cov[i][k] * matrix_cov[j][k];
        }
        matrix_cov[i][j] = sum_up / matrix_cov[j][j];
    }
}
// generate vector Y
for (int i = 0; i < n; i++){
    double sum = generate_uniform(n);
    for (int j = 0; j < CLT_NUM - 1; j++){
        sum += generate_uniform(n);
    }
    vec_Y[i] = (sum - CLT_NUM / 2) * sqrt(12.0 / CLT_NUM);
}
// linear transformation
for (int i = 0; i < n; i++){  
    vec_X[i] = vec_exp[i];
    for (int k = 0; k <= i; k++){
        vec_X[i] += matrix_cov[i][k] * vec_Y[k];
    }
}

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

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

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

  • вычисление элементов матрицы [math]B[/math] в рамках столбца (под диагональю),
  • генерация вектора [math]Y[/math],
  • вычисление вектора [math]X[/math] (умножение матрицы на вектор, сложение векторов).

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

2.4.1 Масштабируемость алгоритма

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

Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета. Алгоритм был реализован на языке C++ с использованием Intel MPI Library и запущен на кластере regular4.

Параметры компиляции и запуска: 

module add slurm
module add impi/5.0.1
mpicxx -std=c++0x task.c -o _scratch/task
sbatch -n256 impi ./task <matrix_size>

Параметры запуска:

  • размер матрицы [1024 : 4864] с шагом 256,
  • число ядер [2 : 256] с шагом по степеням 2.


Рисунок 4. Изменение производительности в зависимости от числа ядер и размера матрицы.
Рисунок 5. Изменение эффективности реализации в зависимости от числа ядер и размера матрицы.
Рисунок 6. Изменение эффективности реализации, вид сверху.

Оценка масштабируемости:

  • По размеру задачи — при увеличении размерности матрицы эффективность реализации постепенно увеличивается, при этом тем быстрее, чем большее число ядер задействовано.
  • По числу ядер — при увеличении числа mpi-процессов эффективность реализации снижается. Это объясняется ростом накладных расходов на организацию параллельного выполнения. При этом с ростом размерности матрицы эффективность снижается с такой же скоростью, но уже при больших значениях числа mpi-процессов.
  • По двум направлениям — эффективность реализации возрастает (на рассматриваемой области значений параметров), но крайне медленно.

Параллельная реализация на C++

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

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

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

MS Excel и распространённые статистические пакеты (напр., SPSS, Statistica) позволяют моделировать только одномерные статистические распределения. Имеется возможность "собрать" многомерное распределение из нескольких одномерных при условии, что компоненты независимы. Однако если необходимо исследовать данные с зависящими друг от друга переменными, то приходится писать собственную программу. При этом удобно воспользоваться готовой реализацией разложения Холецкого (точечного метода), представленной в таких прикладных пакетах, как LINPACK, LAPACK, SCALAPACK.

3 Литература

  1. ↑ Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование. Методы Монте-Карло — М.: Академия, 2006. — 368 с.
  2. ↑ https://ru.wikipedia.org/wiki/Многомерное_нормальное_распределение
  3. ↑ Балдин К.В., Уткин В.Б. Информационные системы в экономике. — М.:Дашков и Кo, 2008. — 395 с.
Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Lonalone/Генерация_гауссовского_вектора_методом_линейных_преобразований&oldid=26783