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

Материал из Алговики
Перейти к навигации Перейти к поиску

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ \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]\displaystyle{ B }[/math], которая позволяла бы получить искомый вектор [math]\displaystyle{ X }[/math] с требуемыми характеристиками в результате линейного преобразования [math]\displaystyle{ X = BY + M }[/math], где [math]\displaystyle{ Y }[/math] — n-мерный случайный вектор с независимыми нормально распределенными компонентами со стандартными параметрами.

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

[math]\displaystyle{ \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]\displaystyle{ \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]\displaystyle{ Y }[/math] независимы и имеют стандартные параметры, справедливо выражение:

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

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

[math]\displaystyle{ \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]\displaystyle{ \Sigma }[/math], а в правых — элементы искомой матрицы [math]\displaystyle{ B }[/math]. Последовательно решая эту систему, получаем формулы для расчета элементов [math]\displaystyle{ \sigma_{ij} }[/math]:

[math]\displaystyle{ 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]\displaystyle{ B }[/math] имеет вид:

[math]\displaystyle{ 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 }[/math]

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

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

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

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

[math]\displaystyle{ 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]\displaystyle{ u_j^i - }[/math] реализации случайной величины [math]\displaystyle{ u }[/math], а [math]\displaystyle{ \kappa - }[/math] параметр, обеспечивающий качество приближения.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Для распараллеливания доступны следующие опции:

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

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

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

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

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

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

// 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];
    }
}

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

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

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

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

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

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

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

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


Рисунок 1. Изменение производительности в зависимости от числа процессоров и размера матрицы.

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

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

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

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

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

Литература

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