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

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показано 35 промежуточных версий этого же участника)
Строка 6: Строка 6:
 
Из многомерных распределений особый интерес представляет нормальное. Этому закону подчиняются все результаты воздействия большого числа случайных факторов, среди которых нет превалирующих.
 
Из многомерных распределений особый интерес представляет нормальное. Этому закону подчиняются все результаты воздействия большого числа случайных факторов, среди которых нет превалирующих.
  
В статье приведен алгоритм генерации n-мерного гауссовского случайного вектора с помощью ''метода линейных преобразований''<ref name="V">(п. 1.10.4) Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование. Методы Монте-Карло — М.: Академия, 2006. — 368 с.</ref>. Известно, что случае нормально распределенного случайного вектора, ковариационная матрица вместе с математическим ожиданием этого вектора полностью определяют его распределение. Поэтому для полного статистического соответствия моделируемого и теоретического распределения гауссовского вектора достаточно обеспечить требуемые значения указанных параметров<ref name="VV"> https://ru.wikipedia.org/wiki/Многомерное_нормальное_распределение</ref>.
+
В статье приведен алгоритм генерации n-мерного гауссовского случайного вектора с помощью ''метода линейных преобразований''<ref name="V">Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование. Методы Монте-Карло — М.: Академия, 2006. — 368 с.</ref>. Известно, что случае нормально распределенного случайного вектора, ковариационная матрица вместе с математическим ожиданием этого вектора полностью определяют его распределение. Поэтому для полного статистического соответствия моделируемого и теоретического распределения гауссовского вектора достаточно обеспечить требуемые значения указанных параметров<ref name="VV"> https://ru.wikipedia.org/wiki/Многомерное_нормальное_распределение</ref>.
  
 
Идея алгоритма заключается в линейном преобразовании n-мерного случайного вектора <math>Y</math>, компоненты которого независимы и одинаково распределены по нормальному закону со стандартными параметрами, в случайный вектор <math>X</math> с требуемыми ковариационной матрицей и вектором математических ожиданий.
 
Идея алгоритма заключается в линейном преобразовании n-мерного случайного вектора <math>Y</math>, компоненты которого независимы и одинаково распределены по нормальному закону со стандартными параметрами, в случайный вектор <math>X</math> с требуемыми ковариационной матрицей и вектором математических ожиданий.
Строка 88: Строка 88:
 
:<math>
 
:<math>
  
b_{ij} = \frac {\sigma_{ij} - \sum_{k=1}^{j-1} b_{ik} b_{jk}} {\sqrt{\sigma_{ij} - \sum_{k=1}^{j-1} b_{jk}^2}}, \quad 1 \leqslant j \leqslant i \leqslant n
+
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>
Строка 98: Строка 98:
 
==== Генерация случайного вектора Y ====
 
==== Генерация случайного вектора Y ====
  
Пусть имеется генератор псевдослучайных чисел (ГПСЧ, PRNG), с помощью которого можно получить реализацию случайной величины <math>u \sim U(0,1) </math>. Описанный выше случайный вектор <math>Y=(y_1, \cdots, y_n)</math> с независимыми компонентами составим из <math>n</math> реализаций случайной величины <math>\eta \sim N(0,1)</math>. Каждую такую реализацию <math>y_i</math>, в свою очередь, получим с помощью приближения по ЦПТ <math>\kappa</math> случайными величинами, распределенными равномерно на отрезке [0,1]:
+
Пусть имеется генератор псевдослучайных чисел (ГПСЧ, 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>
 
:<math>
Строка 106: Строка 106:
 
где <math>u_j^i - </math> реализации случайной величины <math>u</math>, а <math>\kappa - </math> параметр, обеспечивающий качество приближения.
 
где <math>u_j^i - </math> реализации случайной величины <math>u</math>, а <math>\kappa - </math> параметр, обеспечивающий качество приближения.
  
Как правило, <math>\kappa</math> берут равным 12 и считают, что для подавляющего числа практических задач обеспечивается должная точность вычислений<ref name = "VVV">Балдин К.В., Уткин В.Б. Информационные системы в экономике. — М.:Дашков и Кo, 2008. — 395 с.</ref>.
+
Как правило, <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>\sigma_{ij} - \sum_{k=1}^{j-1} b_{ik} b_{jk}, \quad 1 \leqslant j \leqslant i \leqslant n</math>
Строка 146: Строка 148:
 
=== Последовательная сложность алгоритма ===
 
=== Последовательная сложность алгоритма ===
  
* <math> \quad \frac {n(n-1)}{2} + \sum_{j=1}^{n-2}j(n-j-1) + n + \frac {n(n+1)}{2} =
+
Для генерации гауссовского случайного вектора порядка <math>n</math> в последовательном варианте требуется:
  n(n+1) + \frac {n(n-1)(n-2)}{6}</math> (умножений)
+
* <math> n(n+1) + \frac {n(n-1)(n-2)}{6} </math> (умножений)
* <math> \quad \frac {n(n-1)}{2} + \sum_{j=1}^{n-2}j(n-j-1) + n =
+
* <math> \frac {n(n+1)}{2} + \frac {n(n-1)(n-2)}{6} </math> (вычитаний)
  \frac {n(n+1)}{2} + \frac {n(n-1)(n-2)}{6}</math> (вычитаний)
+
* <math> \frac {n(n-1+2\kappa)}{2} </math> (сложений)
* <math> \quad (\kappa-1)n + \frac {n(n+1)}{2} = \frac {n(n-1+2\kappa)}{2}</math> (сложений)
+
* <math> \frac {n(n+3)}{2} </math> (делений)
* <math> \quad \frac {n(n-1)}{2} + 2n = \frac {n(n+3)}{2}</math> (делений)
+
* <math> 2n </math> (вычислений квадратного корня)
* <math> \quad n + n  = 2n </math> (взятия квадратного корня)
 
  
 
<math> \frac{n}{6}\left(2n^2+9n+31+6\kappa\right) </math> - итого операций.
 
<math> \frac{n}{6}\left(2n^2+9n+31+6\kappa\right) </math> - итого операций.
 +
 +
Таким образом, последовательная сложность алгоритма равна <math>O(n^3)</math>.
  
 
=== Информационный граф ===
 
=== Информационный граф ===
Строка 160: Строка 163:
 
==== Заполнение матрицы B====
 
==== Заполнение матрицы B====
  
Граф этапа 1 макроструктуры алгоритма (''разложение Холецкого'') состоит из трёх групп вершин, расположенных в целочисленных узлах трёх областей разной размерности.
+
Граф этапа 1 макроструктуры алгоритма ([[Разложение_Холецкого_(метод_квадратного_корня)|''разложения Холецкого'']]) состоит из трёх групп вершин, расположенных в целочисленных узлах трёх областей разной размерности.
  
 
'''Первая''' группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию SQRT.  
 
'''Первая''' группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию SQRT.  
Строка 201: Строка 204:
 
Описанный граф можно посмотреть на рис.1, выполненном для случая <math>n = 4</math>.  
 
Описанный граф можно посмотреть на рис.1, выполненном для случая <math>n = 4</math>.  
  
[[file:Cholesky full_in_out.png|thumb|center|1400px|Рисунок 1. Граф этапа 1 макроструктуры алгоритма. SQ - вычисление квадратного корня, F - операция a-bc, Div - деление, In - входные данные, Out - выходные.]]
+
[[file:Cholesky full_in_out.png|thumb|center|1400px|Рисунок 1. Граф этапа 1 макроструктуры алгоритма. SQ вычисление квадратного корня, F операция a-bc, Div деление, In входные данные, Out выходные.]]
  
 
====Генерация вектора Y====
 
====Генерация вектора 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====
 
====Вычисление вектора X====
  
Граф этапа 3 макроструктуры алгоритма (умножение нижней треугольной матрицы на вектор и сложение результата с другим вектором) состоит из одной группы вершин, расположенной в целочисленных узлах двумерной области, соответствующая ей операция <math>a+bc</math>.  
+
Граф этапа 3 макроструктуры алгоритма (''умножение нижней треугольной матрицы на вектор и сложение результата с другим вектором'') состоит из одной группы вершин, расположенной в целочисленных узлах двумерной области, соответствующая ей операция <math>a+bc</math>.  
  
 
Естественно введённые координаты области таковы:  
 
Естественно введённые координаты области таковы:  
Строка 226: Строка 256:
 
Описанный граф можно посмотреть на рис.3, выполненном для случая <math> n = 4</math>.
 
Описанный граф можно посмотреть на рис.3, выполненном для случая <math> n = 4</math>.
  
[[Файл:LonGraf3.png|750px|thumb|center|Рисунок 3. Граф этапа 3 макроструктуры алгоритма,  F - операция a+bc]]
+
[[Файл:LonGraf3.png|650px|thumb|center|Рисунок 3. Граф этапа 3 макроструктуры алгоритма,  F операция a+bc]]
  
 
=== Ресурс параллелизма алгоритма ===
 
=== Ресурс параллелизма алгоритма ===
  
Для распараллеливания доступны следующие опции:
+
Для разложения Холецкого <math>(n \times n)</math> матрицы <math>\Sigma</math> в параллельном варианте требуется последовательно выполнить следующие ярусы:
* вычисление элементов матрицы <math>B</math> в рамках столбца (под диагональю),
+
* <math>n</math> ярусов с вычислением квадратного корня (единичные вычисления в каждом из ярусов),
* генерация вектора <math>Y</math>,
+
* <math>n - 1</math> ярус делений (в каждом из ярусов — линейное количество операций, в зависимости от яруса — от <math>1</math> до <math>n - 1</math>),
* вычисление вектора <math>X</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>.
  
 
=== Входные и выходные данные алгоритма ===
 
=== Входные и выходные данные алгоритма ===
Строка 285: Строка 324:
  
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
 +
В данной реализации алгоритма весь ресурс параллелизма использован в следующих этапах:
 +
* вычисление элементов матрицы <math>B</math> в рамках столбца (под диагональю),
 +
* генерация вектора <math>Y</math>,
 +
* вычисление вектора <math>X</math> (умножение матрицы на вектор, сложение векторов).
  
 
=== Масштабируемость алгоритма и его реализации ===
 
=== Масштабируемость алгоритма и его реализации ===
  
 +
==== Масштабируемость алгоритма ====
 +
 +
==== Масштабируемость реализации алгоритма ====
 
Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета]. Алгоритм был реализован на языке C++ с использованием Intel MPI Library и запущен на кластере regular4.
 
Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета]. Алгоритм был реализован на языке C++ с использованием Intel MPI Library и запущен на кластере regular4.
  
Строка 294: Строка 341:
 
  module add impi/5.0.1
 
  module add impi/5.0.1
 
  mpicxx -std=c++0x task.c -o _scratch/task
 
  mpicxx -std=c++0x task.c -o _scratch/task
cd _scratch
 
 
  sbatch -n256 impi ./task <matrix_size>
 
  sbatch -n256 impi ./task <matrix_size>
  
 
Параметры запуска:
 
Параметры запуска:
 
* размер матрицы [1024 : 4864] с шагом 256,
 
* размер матрицы [1024 : 4864] с шагом 256,
* число ядер [2 : 256] с шагом-умножением в 2 раза.
+
* число ядер [2 : 256] с шагом по степеням 2.
 
<br>
 
<br>
  
[[file:Ops_2.png|thumb|center|700px|Рисунок 1. Изменение производительности в зависимости от числа процессоров и размера матрицы.]]
+
[[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++]
  
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
=== Динамические характеристики и эффективность реализации алгоритма ===
Строка 310: Строка 368:
 
=== Существующие реализации алгоритма ===
 
=== Существующие реализации алгоритма ===
  
В отличие от разложения Холецкого, которое реализовано (точечный метод) во многих прикладных пакетах (напр., LINPACK, LAPACK, SCALAPACK), моделирование гауссовского случайного вектора не так "популярно". MS Excel и распространённые статистические пакеты (напр., SPSS, Statistica) позволяют моделировать только одномерные статистические распределения.
+
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 с.