Алговики - Вклад [ru]

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T12:32:19Z

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

Основные авторы описания: [[Участник:v|В.С.Шумовская]]

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

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

Пусть <math>X_{1},\dots, X_{n}</math> -- независимые и нормально распределенные случайные вектора в <math>\R^p</math> с нулевым средним и матрицей ковариацией <math>\Sigma</math>, она лежит в <math>\R^{pxp}</math> и такая, что ее собственные значения быстро убывают, т.е. 3-5 больших, а остальные, к примеру, в диапазоне от 1 до 3.

К этой выборке применим бутстрэп и вычислим в мире бутстрэпа <math>M</math> матриц ковариаций сигма <math>\Sigma^{o}_{j}, j = 1,\dots,M.</math>

Далее фиксируем <math>r</math>, обозначим за <math>P_{r}, P^{o}_{j}, j = 1,\dots,M</math> - проекторы на r-ое подпространство и вычислим следующие статистики:

<math>S^{o}_{j} = ||P_{r} - P^{o}_{r}||^{2}_{2}</math>

Задача -- вычислить большое количество этих статистик для визуализации их распределения.

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

#Генерация матрицы ковариации.
##Возьмем нужный нам набор собственных значений и поместим их на диагонали матрицы <math>\Lambda.</math>
##Далее генерируем базис <math>p</math>-мерного пространства и поместим его в матрицу <math>U.</math>
##Тогда <math>\Sigma = U\Lambda U^{T}.</math>
#Проделываем следующий цикл N раз:
##Генерация векторов <math>X</math> ~ <math>N(\theta, \Sigma).</math>
##Проделываем следующий цикл M раз:
###Генерация весов <math>w</math> ~ <math>N(1, 1).</math>
###Вычисление матриц ковариаций в бутстрэпе <math>\Sigma^{o} = \frac{1}{n}\sum_{k = 1}^{n}w_{k}X_{k}X^{T}_{k}. </math>
###Вычисление собственных векторов <math>u_{1},\dots,u_{p}</math> (в порядке убывания).
###Вычисление проекторов (считаем, что одномерные) по формуле <math>u = u_{r}u^{T}_{r}.</math>
###Вычисление статистик <math>S^{o} = ||P_{r} - P^{o}_{r}||^{2}_{2}.</math>

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

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

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

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

[[Файл:послед.png|center|500px]]

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

Сложность последовательного алгоритма порядка <math>O(n^{3})</math>

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

[[Файл:Парал.png|center|700px]]

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

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

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

* Размерности матриц
* Количество повторений циклов
* Собственные значения
Эти параметры меняются в коде

=== Выходные данные ===

От каждого вычислительного узла будет получена матрица значений статистик

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

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

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

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

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

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

Для исследования масштабируемости были проведены серии тестов с различным числом узлов (в степенях 2). Чем больше узлов, тем меньше загруженность каждого, если хотим получить одинаковый совокупный объем данных, а, следовательно, программа выполняется быстрее.
График для фиксированных n = 100, p = 10.
[[Файл:Graphhh.jpg|center]]

Все тесты проводились на [https://ru.wikipedia.org/wiki/Ломоносов_(суперкомпьютер) суперкомпьютере Ломоносов].
Характеристики вычислительного кластера: [http://users.parallel.ru/wiki/pages/22-config link].

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

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

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

<source lang="c++">
#include "mkl.h"
#include <ctime>
#include <cstdlib>
#include <iostream>
#include <string.h>
#include <sstream>
#include "mpi.h"
#include <fstream>
#include <algorithm>
#include <iterator>
#include <string>
#include <vector>

using namespace std;

void diagMatrix(int p, int k, float *m);
void ortBasis(int p, float *m);
void out(int p, float *m);
void choleskyFactorization(int p, float *m);
void copy(int p, float *from, float *to);
void samplesGenerate(int n, int p, float *mean, float *chol, float *res);
void wGenerate(int n, float* w);
void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result);
void sampleCovariance(int n, int p, float* Samples, float* Result);
void eigenDecomposition(int p, float* m, float *values);
void rProjector(int p, int r, float *vectors, float* projector);
float S(int n, int p, float *m1, float *m2);

void diagMatrix(int p, int k, float *m)
{
float numbers[5] = {36, 30, 25, 19, 15};
int ind = 0;
for (int i = 0; i < p; ++i)
{
for (int j = 0; j < p; ++j)
m[i*p + j] = 0.0;
if (k > 0)
{
m[i*p + i] = numbers[ind];
ind++;
k--;
}
else
m[i*p + i] = 1.0 + 5.0*rand() / float(RAND_MAX);
}
}

void ortBasis(int p, float *m)
{
for (int i = 0; i < p; i++)
for(int j = 0; j < p; j++)
m[i*p + j] = rand();
float* tau = (float*) calloc(p, sizeof(float));
LAPACKE_sgeqrf(LAPACK_ROW_MAJOR, p, p, m, p, tau);
LAPACKE_sorgqr(LAPACK_ROW_MAJOR, p, p, 1, m, p, tau);
free(tau);
}

void out(int r, int c, float *m)
{
for (int i = 0; i < r; i++)
{
cout << "{ ";
for (int j = 0; j < c; j++)
cout << m[i*c + j] << " ";
cout << "}" << endl;
}
}

void choleskyFactorization(int p, float *m)
{
lapack_int q = LAPACKE_spotrf(LAPACK_ROW_MAJOR, 'L', p, m, p);
}

void copy(int p, float *from, float *to)
{
mkl_somatcopy('R', 'N', p, p, 1.0, from, p, to, p);
}

void samplesGenerate(int n, int p, float *mean, float *chol, float *res)
{
VSLStreamStatePtr stream;
int i = (int) rand();
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, n, res, p, VSL_MATRIX_STORAGE_FULL, mean, chol);
}

void wGenerate(int n, float* w)
{
float *mean = (float *) calloc(n, sizeof(float));
float *cov = (float*) calloc(n * n, sizeof(float));
for (int i = 1; i < n; i++)
{
mean[i] = 1;
cov[i*n + i] = 1;
}
int i = (int) rand();
VSLStreamStatePtr stream;
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, 1, w, n, VSL_MATRIX_STORAGE_FULL, mean, cov);
free(mean);
free(cov);
}

void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *W = (float *) calloc(n, sizeof(float));
float *X = (float *) calloc(p, sizeof(float));
wGenerate(n, W);
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, W[k], X, p, X, p, 1.0, Result, p);
}
free(W);
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}

free(X);
}

void sampleCovariance(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *X = (float *) calloc(p, sizeof(float));
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, Result, p);
}
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}
free(X);
}

void eigenDecomposition(int p, float* m, float *values) //decreasing orderб столбыц!
{
LAPACKE_ssyevd(LAPACK_ROW_MAJOR, 'V', 'U', p, m, p, values);
}

void rProjector(int p, int r, float *vectors, float* projector)
{
int num = p - r;
float *X = (float *) calloc(p, sizeof(float));
for (int i = 0; i < p; i++)
{
X[i] = vectors[p*i + num];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, projector, p);
free(X);
}

float S(int n, int p, float *m1, float *m2)
{
float *m = (float *) calloc(p * p, sizeof(float));
for (int i = 0; i < p; i++)
for (int j = 0; j < p; j++)
m[i*p + j] = m1[i*p + j] - m2[i*p + j];
float norm = LAPACKE_slange(LAPACK_ROW_MAJOR, 'F', p, p, m, p);
free(m);
return n*norm*norm;
}

int main(int argc, char **argv)
{
//задаем размерности
int paral = 16;
int p = 20;
int k = 3;
int M = 3000;
int r = 1; //нумерация с 1
int N = 3; //45
int n = 100;
int rank, size, prev, next;

//инициализация MPI
MPI_Request reqs[4];
MPI_Status stats[4];
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//
double MPI_Wtime(void);
double t1, t2, dt;
//
t1 = MPI_Wtime();
prev = rank - 1;
next = rank + 1;
if(rank==0) prev = size - 1;
if(rank==size - 1) next = 0;

//чтобы у всех сгенерировалась одинаковая матрица ковариации srand сработает одновременно
srand(time(NULL));
//файлы, в который будем писать результаты, у каждого параллельного процесса свой (номер процесса + .txt)
ofstream myfile;
char intStr[2];
if (rank < 10)
{
intStr[0] = '0';
intStr[1] = '0' + rank;
}
else
{
intStr[0] = '1';
intStr[1] = '0' + rank - 10;
}
char *txt = ".txt";
char *str = strcat(intStr,txt);

myfile.open(str);

//файлы, в которые будем отправлять размерности и тд, за это будет отвечать нулевой процесс
ofstream dim;
ofstream val;

if (rank == 0){
dim.open("dimensions.txt");
val.open("values.txt");
}

if (rank == 0)
{
dim << "p = " << p << endl;
dim << "k = " << k << endl;
dim << "M = " << M << endl;
dim << "r = " << r << endl;
dim << "N = " << N*(paral-1) << endl << endl;
dim << "n = " << n;
}

//генерация матрицы ковариации
float *COV = (float *) calloc(p * p, sizeof(float));
float *MEAN = (float *) calloc(p, sizeof(float));
float *COV_Chol = (float *) calloc(p * p, sizeof(float));
float *Basis = (float *) calloc(p * p, sizeof(float));
float *COV_L = (float *) calloc(p * p, sizeof(float));

ortBasis(p, Basis);

diagMatrix(p, k, COV_L);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, Basis, p, COV_L, p, 0.0, COV, p);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, COV, p, Basis, p, 0.0, COV_L, p);
copy(p, COV_L, COV);

copy(p, COV, COV_Chol);
choleskyFactorization(p, COV_Chol);
free(Basis);
free(COV_L);

float *values_pr = (float *) calloc(p, sizeof(float));
float *vectors_pr = (float *) calloc(p * p, sizeof(float));

copy(p, COV, vectors_pr);
eigenDecomposition(p, vectors_pr, values_pr);

//для исследования нам понадобятся собственные значения, поэтому их тоже выводит нулевой процесс
if(rank ==0){
for (int i = 0; i < p; i++)
{
val << values_pr[i] << " ";
}
val << endl << endl;
}

free(values_pr);
free(vectors_pr);

//вычислим проектор для реальной матрицы ковариации
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorR = (float *) calloc(p * p, sizeof(float));

copy(p, COV, vectors);
eigenDecomposition(p, vectors, values);
rProjector(p, r, vectors, ProjectorR);
//ProjectorR - нужный проектор
free(values);
free(vectors);

//для каждого процесса запустим свой рандом (зависящий от ранка), иначе все значения у каждого процесса будут одинаковые
srand(time(NULL) + rank);

float* NormsBootstrapWorld = (float *) calloc(N * M, sizeof(float));
for (int i = 0; i < N; i++)
{
float *Samples = (float *) calloc(n * p, sizeof(float)); //Каждый Xj - строка в матрице
for (int j = 0; j < M; j++)
{
float *COV_SampleBootstrap = (float *) calloc(p * p, sizeof(float));
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorSB = (float *) calloc(p * p, sizeof(float));

sampleCovarianceBootstrap(n, p, Samples, COV_SampleBootstrap);
copy(p, COV_SampleBootstrap, vectors);
eigenDecomposition(p, vectors, values); //вектора в столбцах
rProjector(p, r, vectors, ProjectorSB);
//ProjectorSB - нужный проектор

float s = S(n, p, ProjectorR, ProjectorSB);
NormsBootstrapWorld[i*M + j] = s;

free(COV_SampleBootstrap);
free(values);
free(vectors);
free(ProjectorSB);
}
free(Samples);
}

free(ProjectorR);

for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
myfile << NormsBootstrapWorld[j*M + i] << " ";
}
myfile << endl;
}
//в myfile N столбцов

myfile << endl;
free(NormsBootstrapWorld);

free(COV);
free(MEAN);
free(COV_Chol);
t2 = MPI_Wtime();
dt = t2 - t1;
if (rank == 0)
cout << dt << endl;
MPI_Finalize(); //конец параллеливания
}
</source>

= Литература =

* [https://en.wikipedia.org/wiki/Bootstrapping_(statistics) Boostrapping (wiki)]
* [https://software.intel.com/sites/default/files/managed/ff/c8/mkl-2017-developer-reference-c_0.pdf Developer Reference for Intel® Math Kernel Library 2017 - C(pdf)]

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T11:45:03Z

Основные авторы описания: [[Участник:v|В.С.Шумовская]]

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

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

Пусть <math>X_{1},\dots, X_{n}</math> -- независимые и нормально распределенные случайные вектора в <math>\R^p</math> с нулевым средним и матрицей ковариацией <math>\Sigma</math>, она лежит в <math>\R^{pxp}</math> и такая, что ее собственные значения быстро убывают, т.е. 3-5 больших, а остальные, к примеру, в диапазоне от 1 до 3.

К этой выборке применим бутстрэп и вычислим в мире бутстрэпа <math>M</math> матриц ковариаций сигма <math>\Sigma^{o}_{j}, j = 1,\dots,M.</math>

Далее фиксируем <math>r</math>, обозначим за <math>P_{r}, P^{o}_{j}, j = 1,\dots,M</math> - проекторы на r-ое подпространство и вычислим следующие статистики:

<math>S^{o}_{j} = ||P_{r} - P^{o}_{r}||^{2}_{2}</math>

Задача -- вычислить большое количество этих статистик для визуализации их распределения.

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

#Генерация матрицы ковариации.
##Возьмем нужный нам набор собственных значений и поместим их на диагонали матрицы <math>\Lambda.</math>
##Далее генерируем базис <math>p</math>-мерного пространства и поместим его в матрицу <math>U.</math>
##Тогда <math>\Sigma = U\Lambda U^{T}.</math>
#Проделываем следующий цикл N раз:
##Генерация векторов <math>X</math> ~ <math>N(\theta, \Sigma).</math>
##Проделываем следующий цикл M раз:
###Генерация весов <math>w</math> ~ <math>N(1, 1).</math>
###Вычисление матриц ковариаций в бутстрэпе <math>\Sigma^{o} = \frac{1}{n}\sum_{k = 1}^{n}w_{k}X_{k}X^{T}_{k}. </math>
###Вычисление собственных векторов <math>u_{1},\dots,u_{p}</math> (в порядке убывания).
###Вычисление проекторов (считаем, что одномерные) по формуле <math>u = u_{r}u^{T}_{r}.</math>
###Вычисление статистик <math>S^{o} = ||P_{r} - P^{o}_{r}||^{2}_{2}.</math>

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

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

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

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

[[Файл:послед.png|center|500px]]

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

Сложность последовательного алгоритма порядка <math>O(n^{3})</math>

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

[[Файл:Парал.png|center|700px]]

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

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

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

* Размерности матриц
* Количество повторений циклов
* Собственные значения
Эти параметры меняются в коде

=== Выходные данные ===

От каждого вычислительного узла будет получена матрица значений статистик

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

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

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

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

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

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

Для исследования масштабируемости были проведены серии тестов с различным числом узлов (в степенях 2). Чем больше узлов, тем меньше загруженность каждого, если хотим получить одинаковый совокупный объем данных, а, следовательно, программа выполняется быстрее.
График для фиксированных n = 100, p = 10.
[[Файл:Graphhh.jpg|center]]

Все тесты проводились на [https://ru.wikipedia.org/wiki/Ломоносов_(суперкомпьютер) суперкомпьютере Ломоносов].
Характеристики вычислительного кластера: [http://users.parallel.ru/wiki/pages/22-config link].

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

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

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

<source lang="c++">
#include "mkl.h"
#include <ctime>
#include <cstdlib>
#include <iostream>
#include <string.h>
#include <sstream>
#include "mpi.h"
#include <fstream>
#include <algorithm>
#include <iterator>
#include <string>
#include <vector>

using namespace std;

void diagMatrix(int p, int k, float *m);
void ortBasis(int p, float *m);
void out(int p, float *m);
void choleskyFactorization(int p, float *m);
void copy(int p, float *from, float *to);
void samplesGenerate(int n, int p, float *mean, float *chol, float *res);
void wGenerate(int n, float* w);
void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result);
void sampleCovariance(int n, int p, float* Samples, float* Result);
void eigenDecomposition(int p, float* m, float *values);
void rProjector(int p, int r, float *vectors, float* projector);
float S(int n, int p, float *m1, float *m2);
float **alloc2d(int n, int m);

float **alloc2d(int n, int m)
{
float *data = (float*) malloc(n*m*sizeof(float));
float **array = (float **)malloc(n*sizeof(float *));
for (int i=0; i<n; i++)
array[i] = &(data[i*m]);
return array;
}

struct aboutMatrix
{
float* matrix;
float* l;
float* u;
float* cholesky;
};

void diagMatrix(int p, int k, float *m)
{
float numbers[5] = {36, 30, 25, 19, 15};
int ind = 0;
for (int i = 0; i < p; ++i)
{
for (int j = 0; j < p; ++j)
m[i*p + j] = 0.0;
if (k > 0)
{
m[i*p + i] = numbers[ind];
ind++;
k--;
}
else
m[i*p + i] = 1.0 + 5.0*rand() / float(RAND_MAX);
}
}

void ortBasis(int p, float *m)
{
for (int i = 0; i < p; i++)
for(int j = 0; j < p; j++)
m[i*p + j] = rand();
float* tau = (float*) calloc(p, sizeof(float));
LAPACKE_sgeqrf(LAPACK_ROW_MAJOR, p, p, m, p, tau);
LAPACKE_sorgqr(LAPACK_ROW_MAJOR, p, p, 1, m, p, tau);
free(tau);
}

void out(int r, int c, float *m)
{
for (int i = 0; i < r; i++)
{
cout << "{ ";
for (int j = 0; j < c; j++)
cout << m[i*c + j] << " ";
cout << "}" << endl;
}
}

void choleskyFactorization(int p, float *m)
{
lapack_int q = LAPACKE_spotrf(LAPACK_ROW_MAJOR, 'L', p, m, p);
}

void copy(int p, float *from, float *to)
{
mkl_somatcopy('R', 'N', p, p, 1.0, from, p, to, p);
}

void samplesGenerate(int n, int p, float *mean, float *chol, float *res)
{
VSLStreamStatePtr stream;
int i = (int) rand();
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, n, res, p, VSL_MATRIX_STORAGE_FULL, mean, chol);
}

void wGenerate(int n, float* w)
{
float *mean = (float *) calloc(n, sizeof(float));
float *cov = (float*) calloc(n * n, sizeof(float));
for (int i = 1; i < n; i++)
{
mean[i] = 1;
cov[i*n + i] = 1;
}
int i = (int) rand();
VSLStreamStatePtr stream;
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, 1, w, n, VSL_MATRIX_STORAGE_FULL, mean, cov);
free(mean);
free(cov);
}

void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *W = (float *) calloc(n, sizeof(float));
float *X = (float *) calloc(p, sizeof(float));
wGenerate(n, W);
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, W[k], X, p, X, p, 1.0, Result, p);
}
free(W);
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}

free(X);
}

void sampleCovariance(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *X = (float *) calloc(p, sizeof(float));
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, Result, p);
}
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}
free(X);
}

void eigenDecomposition(int p, float* m, float *values) //decreasing orderб столбыц!
{
LAPACKE_ssyevd(LAPACK_ROW_MAJOR, 'V', 'U', p, m, p, values);
}

void rProjector(int p, int r, float *vectors, float* projector)
{
int num = p - r;
float *X = (float *) calloc(p, sizeof(float));
for (int i = 0; i < p; i++)
{
X[i] = vectors[p*i + num];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, projector, p);
free(X);
}

float S(int n, int p, float *m1, float *m2)
{
float *m = (float *) calloc(p * p, sizeof(float));
for (int i = 0; i < p; i++)
for (int j = 0; j < p; j++)
m[i*p + j] = m1[i*p + j] - m2[i*p + j];
float norm = LAPACKE_slange(LAPACK_ROW_MAJOR, 'F', p, p, m, p);
free(m);
return n*norm*norm;
}

int main(int argc, char **argv)
{
//задаем размерности
int paral = 16;
int p = 20;
int k = 3;
int M = 3000;
int r = 1; //нумерация с 1
int N = 3; //45
int n = 100;
int rank, size, prev, next;

//инициализация MPI
MPI_Request reqs[4];
MPI_Status stats[4];
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//
double MPI_Wtime(void);
double t1, t2, dt;
//
t1 = MPI_Wtime();
prev = rank - 1;
next = rank + 1;
if(rank==0) prev = size - 1;
if(rank==size - 1) next = 0;

//чтобы у всех сгенерировалась одинаковая матрица ковариации srand сработает одновременно
srand(time(NULL));
//файлы, в который будем писать результаты, у каждого параллельного процесса свой (номер процесса + .txt)
ofstream myfile;
char intStr[2];
if (rank < 10)
{
intStr[0] = '0';
intStr[1] = '0' + rank;
}
else
{
intStr[0] = '1';
intStr[1] = '0' + rank - 10;
}
char *txt = ".txt";
char *str = strcat(intStr,txt);

myfile.open(str);

//файлы, в которые будем отправлять размерности и тд, за это будет отвечать нулевой процесс
ofstream dim;
ofstream val;

if (rank == 0){
dim.open("dimensions.txt");
val.open("values.txt");
}

if (rank == 0)
{
dim << "p = " << p << endl;
dim << "k = " << k << endl;
dim << "M = " << M << endl;
dim << "r = " << r << endl;
dim << "N = " << N*(paral-1) << endl << endl;
dim << "n = " << n;
}

//генерация матрицы ковариации
float *COV = (float *) calloc(p * p, sizeof(float));
float *MEAN = (float *) calloc(p, sizeof(float));
float *COV_Chol = (float *) calloc(p * p, sizeof(float));
float *Basis = (float *) calloc(p * p, sizeof(float));
float *COV_L = (float *) calloc(p * p, sizeof(float));

ortBasis(p, Basis);

diagMatrix(p, k, COV_L);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, Basis, p, COV_L, p, 0.0, COV, p);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, COV, p, Basis, p, 0.0, COV_L, p);
copy(p, COV_L, COV);

copy(p, COV, COV_Chol);
choleskyFactorization(p, COV_Chol);
free(Basis);
free(COV_L);

float *values_pr = (float *) calloc(p, sizeof(float));
float *vectors_pr = (float *) calloc(p * p, sizeof(float));

copy(p, COV, vectors_pr);
eigenDecomposition(p, vectors_pr, values_pr);

//для исследования нам понадобятся собственные значения, поэтому их тоже выводит нулевой процесс
if(rank ==0){
for (int i = 0; i < p; i++)
{
val << values_pr[i] << " ";
}
val << endl << endl;
}

free(values_pr);
free(vectors_pr);

//вычислим проектор для реальной матрицы ковариации
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorR = (float *) calloc(p * p, sizeof(float));

copy(p, COV, vectors);
eigenDecomposition(p, vectors, values);
rProjector(p, r, vectors, ProjectorR);
//ProjectorR - нужный проектор
free(values);
free(vectors);

//для каждого процесса запустим свой рандом (зависящий от ранка), иначе все значения у каждого процесса будут одинаковые
srand(time(NULL) + rank);

float* NormsBootstrapWorld = (float *) calloc(N * M, sizeof(float));
for (int i = 0; i < N; i++)
{
float *Samples = (float *) calloc(n * p, sizeof(float)); //Каждый Xj - строка в матрице
for (int j = 0; j < M; j++)
{
float *COV_SampleBootstrap = (float *) calloc(p * p, sizeof(float));
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorSB = (float *) calloc(p * p, sizeof(float));

sampleCovarianceBootstrap(n, p, Samples, COV_SampleBootstrap);
copy(p, COV_SampleBootstrap, vectors);
eigenDecomposition(p, vectors, values); //вектора в столбцах
rProjector(p, r, vectors, ProjectorSB);
//ProjectorSB - нужный проектор

float s = S(n, p, ProjectorR, ProjectorSB);
NormsBootstrapWorld[i*M + j] = s;

free(COV_SampleBootstrap);
free(values);
free(vectors);
free(ProjectorSB);
}
free(Samples);
}

free(ProjectorR);

for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
myfile << NormsBootstrapWorld[j*M + i] << " ";
}
myfile << endl;
}
//в myfile N столбцов

myfile << endl;
free(NormsBootstrapWorld);

free(COV);
free(MEAN);
free(COV_Chol);
t2 = MPI_Wtime();
dt = t2 - t1;
if (rank == 0)
cout << dt << endl;
MPI_Finalize(); //конец параллеливания
}
</source>

= Литература =

* [https://en.wikipedia.org/wiki/Bootstrapping_(statistics) Boostrapping (wiki)]
* [https://software.intel.com/sites/default/files/managed/ff/c8/mkl-2017-developer-reference-c_0.pdf Developer Reference for Intel® Math Kernel Library 2017 - C(pdf)]

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T11:34:05Z

Основные авторы описания: [[Участник:v|В.С.Шумовская]]

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

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

Пусть <math>X_{1},\dots, X_{n}</math> -- независимые и нормально распределенные случайные вектора в <math>\R^p</math> с нулевым средним и матрицей ковариацией <math>\Sigma</math>, она лежит в <math>\R^{pxp}</math> и такая, что ее собственные значения быстро убывают, т.е. 3-5 больших, а остальные, к примеру, в диапазоне от 1 до 3.

К этой выборке применим бутстрэп и вычислим в мире бутстрэпа <math>M</math> матриц ковариаций сигма <math>\Sigma^{o}_{j}, j = 1,\dots,M.</math>

Далее фиксируем <math>r</math>, обозначим за <math>P_{r}, P^{o}_{j}, j = 1,\dots,M</math> - проекторы на r-ое подпространство и вычислим следующие статистики:

<math>S^{o}_{j} = ||P_{r} - P^{o}_{r}||^{2}_{2}</math>

Задача -- вычислить большое количество этих статистик для визуализации их распределения.

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

#Генерация матрицы ковариации.
##Возьмем нужный нам набор собственных значений и поместим их на диагонали матрицы <math>\Lambda.</math>
##Далее генерируем базис <math>p</math>-мерного пространства и поместим его в матрицу <math>U.</math>
##Тогда <math>\Sigma = U\Lambda U^{T}.</math>
#Проделываем следующий цикл N раз:
##Генерация векторов <math>X</math> ~ <math>N(\theta, \Sigma).</math>
##Проделываем следующий цикл M раз:
###Генерация весов <math>w</math> ~ <math>N(1, 1).</math>
###Вычисление матриц ковариаций в бутстрэпе <math>\Sigma^{o} = \frac{1}{n}\sum_{k = 1}^{n}w_{k}X_{k}X^{T}_{k}. </math>
###Вычисление собственных векторов <math>u_{1},\dots,u_{p}</math> (в порядке убывания).
###Вычисление проекторов (считаем, что одномерные) по формуле <math>u = u_{r}u^{T}_{r}.</math>
###Вычисление статистик <math>S^{o} = ||P_{r} - P^{o}_{r}||^{2}_{2}.</math>

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

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

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

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

[[Файл:послед.png|center|500px]]

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

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

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

[[Файл:Парал.png|center|700px]]

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

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

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

* Размерности матриц
* Количество повторений циклов
* Собственные значения
Эти параметры меняются в коде

=== Выходные данные ===

От каждого вычислительного узла будет получена матрица значений статистик

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

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

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

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

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

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

Для исследования масштабируемости были проведены серии тестов с различным числом узлов (в степенях 2). Чем больше узлов, тем меньше загруженность каждого, если хотим получить одинаковый совокупный объем данных, а, следовательно, программа выполняется быстрее.
График для фиксированных n = 100, p = 10.
[[Файл:Graphhh.jpg|center]]

Все тесты проводились на [https://ru.wikipedia.org/wiki/Ломоносов_(суперкомпьютер) суперкомпьютере Ломоносов].
Характеристики вычислительного кластера: [http://users.parallel.ru/wiki/pages/22-config link].

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

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

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

<source lang="c++">
#include "mkl.h"
#include <ctime>
#include <cstdlib>
#include <iostream>
#include <string.h>
#include <sstream>
#include "mpi.h"
#include <fstream>
#include <algorithm>
#include <iterator>
#include <string>
#include <vector>

using namespace std;

void diagMatrix(int p, int k, float *m);
void ortBasis(int p, float *m);
void out(int p, float *m);
void choleskyFactorization(int p, float *m);
void copy(int p, float *from, float *to);
void samplesGenerate(int n, int p, float *mean, float *chol, float *res);
void wGenerate(int n, float* w);
void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result);
void sampleCovariance(int n, int p, float* Samples, float* Result);
void eigenDecomposition(int p, float* m, float *values);
void rProjector(int p, int r, float *vectors, float* projector);
float S(int n, int p, float *m1, float *m2);
float **alloc2d(int n, int m);

float **alloc2d(int n, int m)
{
float *data = (float*) malloc(n*m*sizeof(float));
float **array = (float **)malloc(n*sizeof(float *));
for (int i=0; i<n; i++)
array[i] = &(data[i*m]);
return array;
}

struct aboutMatrix
{
float* matrix;
float* l;
float* u;
float* cholesky;
};

void diagMatrix(int p, int k, float *m)
{
float numbers[5] = {36, 30, 25, 19, 15};
int ind = 0;
for (int i = 0; i < p; ++i)
{
for (int j = 0; j < p; ++j)
m[i*p + j] = 0.0;
if (k > 0)
{
m[i*p + i] = numbers[ind];
ind++;
k--;
}
else
m[i*p + i] = 1.0 + 5.0*rand() / float(RAND_MAX);
}
}

void ortBasis(int p, float *m)
{
for (int i = 0; i < p; i++)
for(int j = 0; j < p; j++)
m[i*p + j] = rand();
float* tau = (float*) calloc(p, sizeof(float));
LAPACKE_sgeqrf(LAPACK_ROW_MAJOR, p, p, m, p, tau);
LAPACKE_sorgqr(LAPACK_ROW_MAJOR, p, p, 1, m, p, tau);
free(tau);
}

void out(int r, int c, float *m)
{
for (int i = 0; i < r; i++)
{
cout << "{ ";
for (int j = 0; j < c; j++)
cout << m[i*c + j] << " ";
cout << "}" << endl;
}
}

void choleskyFactorization(int p, float *m)
{
lapack_int q = LAPACKE_spotrf(LAPACK_ROW_MAJOR, 'L', p, m, p);
}

void copy(int p, float *from, float *to)
{
mkl_somatcopy('R', 'N', p, p, 1.0, from, p, to, p);
}

void samplesGenerate(int n, int p, float *mean, float *chol, float *res)
{
VSLStreamStatePtr stream;
int i = (int) rand();
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, n, res, p, VSL_MATRIX_STORAGE_FULL, mean, chol);
}

void wGenerate(int n, float* w)
{
float *mean = (float *) calloc(n, sizeof(float));
float *cov = (float*) calloc(n * n, sizeof(float));
for (int i = 1; i < n; i++)
{
mean[i] = 1;
cov[i*n + i] = 1;
}
int i = (int) rand();
VSLStreamStatePtr stream;
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, 1, w, n, VSL_MATRIX_STORAGE_FULL, mean, cov);
free(mean);
free(cov);
}

void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *W = (float *) calloc(n, sizeof(float));
float *X = (float *) calloc(p, sizeof(float));
wGenerate(n, W);
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, W[k], X, p, X, p, 1.0, Result, p);
}
free(W);
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}

free(X);
}

void sampleCovariance(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *X = (float *) calloc(p, sizeof(float));
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, Result, p);
}
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}
free(X);
}

void eigenDecomposition(int p, float* m, float *values) //decreasing orderб столбыц!
{
LAPACKE_ssyevd(LAPACK_ROW_MAJOR, 'V', 'U', p, m, p, values);
}

void rProjector(int p, int r, float *vectors, float* projector)
{
int num = p - r;
float *X = (float *) calloc(p, sizeof(float));
for (int i = 0; i < p; i++)
{
X[i] = vectors[p*i + num];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, projector, p);
free(X);
}

float S(int n, int p, float *m1, float *m2)
{
float *m = (float *) calloc(p * p, sizeof(float));
for (int i = 0; i < p; i++)
for (int j = 0; j < p; j++)
m[i*p + j] = m1[i*p + j] - m2[i*p + j];
float norm = LAPACKE_slange(LAPACK_ROW_MAJOR, 'F', p, p, m, p);
free(m);
return n*norm*norm;
}

int main(int argc, char **argv)
{
//задаем размерности
int paral = 16;
int p = 20;
int k = 3;
int M = 3000;
int r = 1; //нумерация с 1
int N = 3; //45
int n = 100;
int rank, size, prev, next;

//инициализация MPI
MPI_Request reqs[4];
MPI_Status stats[4];
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//
double MPI_Wtime(void);
double t1, t2, dt;
//
t1 = MPI_Wtime();
prev = rank - 1;
next = rank + 1;
if(rank==0) prev = size - 1;
if(rank==size - 1) next = 0;

//чтобы у всех сгенерировалась одинаковая матрица ковариации srand сработает одновременно
srand(time(NULL));
//файлы, в который будем писать результаты, у каждого параллельного процесса свой (номер процесса + .txt)
ofstream myfile;
char intStr[2];
if (rank < 10)
{
intStr[0] = '0';
intStr[1] = '0' + rank;
}
else
{
intStr[0] = '1';
intStr[1] = '0' + rank - 10;
}
char *txt = ".txt";
char *str = strcat(intStr,txt);

myfile.open(str);

//файлы, в которые будем отправлять размерности и тд, за это будет отвечать нулевой процесс
ofstream dim;
ofstream val;

if (rank == 0){
dim.open("dimensions.txt");
val.open("values.txt");
}

if (rank == 0)
{
dim << "p = " << p << endl;
dim << "k = " << k << endl;
dim << "M = " << M << endl;
dim << "r = " << r << endl;
dim << "N = " << N*(paral-1) << endl << endl;
dim << "n = " << n;
}

//генерация матрицы ковариации
float *COV = (float *) calloc(p * p, sizeof(float));
float *MEAN = (float *) calloc(p, sizeof(float));
float *COV_Chol = (float *) calloc(p * p, sizeof(float));
float *Basis = (float *) calloc(p * p, sizeof(float));
float *COV_L = (float *) calloc(p * p, sizeof(float));

ortBasis(p, Basis);

diagMatrix(p, k, COV_L);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, Basis, p, COV_L, p, 0.0, COV, p);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, COV, p, Basis, p, 0.0, COV_L, p);
copy(p, COV_L, COV);

copy(p, COV, COV_Chol);
choleskyFactorization(p, COV_Chol);
free(Basis);
free(COV_L);

float *values_pr = (float *) calloc(p, sizeof(float));
float *vectors_pr = (float *) calloc(p * p, sizeof(float));

copy(p, COV, vectors_pr);
eigenDecomposition(p, vectors_pr, values_pr);

//для исследования нам понадобятся собственные значения, поэтому их тоже выводит нулевой процесс
if(rank ==0){
for (int i = 0; i < p; i++)
{
val << values_pr[i] << " ";
}
val << endl << endl;
}

free(values_pr);
free(vectors_pr);

//вычислим проектор для реальной матрицы ковариации
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorR = (float *) calloc(p * p, sizeof(float));

copy(p, COV, vectors);
eigenDecomposition(p, vectors, values);
rProjector(p, r, vectors, ProjectorR);
//ProjectorR - нужный проектор
free(values);
free(vectors);

//для каждого процесса запустим свой рандом (зависящий от ранка), иначе все значения у каждого процесса будут одинаковые
srand(time(NULL) + rank);

float* NormsBootstrapWorld = (float *) calloc(N * M, sizeof(float));
for (int i = 0; i < N; i++)
{
float *Samples = (float *) calloc(n * p, sizeof(float)); //Каждый Xj - строка в матрице
for (int j = 0; j < M; j++)
{
float *COV_SampleBootstrap = (float *) calloc(p * p, sizeof(float));
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorSB = (float *) calloc(p * p, sizeof(float));

sampleCovarianceBootstrap(n, p, Samples, COV_SampleBootstrap);
copy(p, COV_SampleBootstrap, vectors);
eigenDecomposition(p, vectors, values); //вектора в столбцах
rProjector(p, r, vectors, ProjectorSB);
//ProjectorSB - нужный проектор

float s = S(n, p, ProjectorR, ProjectorSB);
NormsBootstrapWorld[i*M + j] = s;

free(COV_SampleBootstrap);
free(values);
free(vectors);
free(ProjectorSB);
}
free(Samples);
}

free(ProjectorR);

for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
myfile << NormsBootstrapWorld[j*M + i] << " ";
}
myfile << endl;
}
//в myfile N столбцов

myfile << endl;
free(NormsBootstrapWorld);

free(COV);
free(MEAN);
free(COV_Chol);
t2 = MPI_Wtime();
dt = t2 - t1;
if (rank == 0)
cout << dt << endl;
MPI_Finalize(); //конец параллеливания
}
</source>

= Литература =

* [https://en.wikipedia.org/wiki/Bootstrapping_(statistics) Boostrapping (wiki)]
* [https://software.intel.com/sites/default/files/managed/ff/c8/mkl-2017-developer-reference-c_0.pdf Developer Reference for Intel® Math Kernel Library 2017 - C(pdf)]

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T11:32:48Z

Основные авторы описания: [[Участник:v|В.С.Шумовская]]

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

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

Пусть <math>X_{1},\dots, X_{n}</math> -- независимые и нормально распределенные случайные вектора в <math>\R^p</math> с нулевым средним и матрицей ковариацией <math>\Sigma</math>, она лежит в <math>\R^{pxp}</math> и такая, что ее собственные значения быстро убывают, т.е. 3-5 больших, а остальные, к примеру, в диапазоне от 1 до 3.

К этой выборке применим бутстрэп и вычислим в мире бутстрэпа <math>M</math> матриц ковариаций сигма <math>\Sigma^{o}_{j}, j = 1,\dots,M.</math>

Далее фиксируем <math>r</math>, обозначим за <math>P_{r}, P^{o}_{j}, j = 1,\dots,M</math> - проекторы на r-ое подпространство и вычислим следующие статистики:

<math>S^{o}_{j} = ||P_{r} - P^{o}_{r}||^{2}_{2}</math>

Задача -- вычислить большое количество этих статистик для визуализации их распределения.

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

#Генерация матрицы ковариации.
##Возьмем нужный нам набор собственных значений и поместим их на диагонали матрицы <math>\Lambda.</math>
##Далее генерируем базис <math>p</math>-мерного пространства и поместим его в матрицу <math>U.</math>
##Тогда <math>\Sigma = U\Lambda U^{T}.</math>
#Проделываем следующий цикл N раз:
##Генерация векторов <math>X</math> ~ <math>N(\theta, \Sigma).</math>
##Проделываем следующий цикл M раз:
###Генерация весов <math>w</math> ~ <math>N(1, 1).</math>
###Вычисление матриц ковариаций в бутстрэпе <math>\Sigma^{o} = \frac{1}{n}\sum_{k = 1}^{n}w_{k}X_{k}X^{T}_{k}. </math>
###Вычисление собственных векторов <math>u_{1},\dots,u_{p}</math> (в порядке убывания).
###Вычисление проекторов (считаем, что одномерные) по формуле <math>u = u_{r}u^{T}_{r}.</math>
###Вычисление статистик <math>S^{o} = ||P_{r} - P^{o}_{r}||^{2}_{2}.</math>

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

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

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

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

[[Файл:послед.png|center|500px]]

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

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

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

[[Файл:Парал.png|center|700px]]

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

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

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

* Размерности матриц
* Количество повторений циклов
* Собственные значения
Эти параметры меняются в коде

=== Выходные данные ===

От каждого вычислительного узла будет получена матрица значений статистик

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

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

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

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

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

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

Для исследования масштабируемости были проведены серии тестов с различным числом узлов (в степенях 2). Чем больше узлов, тем меньше загруженность каждого, если хотим получить одинаковый совокупный объем данных, а, следовательно, программа выполняется быстрее.
График для фиксированных n = 100, p = 10.
[[Файл:Graphhh.jpg|center]]

Все тесты проводились на [https://ru.wikipedia.org/wiki/Ломоносов_(суперкомпьютер) суперкомпьютере Ломоносов].
Характеристики вычислительного кластера: [http://users.parallel.ru/wiki/pages/22-config link].

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

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

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

<source lang="c++">
#include "mkl.h"
#include <ctime>
#include <cstdlib>
#include <iostream>
#include <string.h>
#include <sstream>
#include "mpi.h"
#include <fstream>
#include <algorithm>
#include <iterator>
#include <string>
#include <vector>

using namespace std;

void diagMatrix(int p, int k, float *m);
void ortBasis(int p, float *m);
void out(int p, float *m);
void choleskyFactorization(int p, float *m);
void copy(int p, float *from, float *to);
void samplesGenerate(int n, int p, float *mean, float *chol, float *res);
void wGenerate(int n, float* w);
void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result);
void sampleCovariance(int n, int p, float* Samples, float* Result);
void eigenDecomposition(int p, float* m, float *values);
void rProjector(int p, int r, float *vectors, float* projector);
float S(int n, int p, float *m1, float *m2);
float **alloc2d(int n, int m);

float **alloc2d(int n, int m)
{
float *data = (float*) malloc(n*m*sizeof(float));
float **array = (float **)malloc(n*sizeof(float *));
for (int i=0; i<n; i++)
array[i] = &(data[i*m]);
return array;
}

struct aboutMatrix
{
float* matrix;
float* l;
float* u;
float* cholesky;
};

void diagMatrix(int p, int k, float *m)
{
float numbers[5] = {36, 30, 25, 19, 15};
int ind = 0;
for (int i = 0; i < p; ++i)
{
for (int j = 0; j < p; ++j)
m[i*p + j] = 0.0;
if (k > 0)
{
m[i*p + i] = numbers[ind];
ind++;
k--;
}
else
m[i*p + i] = 1.0 + 5.0*rand() / float(RAND_MAX);
}
}

void ortBasis(int p, float *m)
{
for (int i = 0; i < p; i++)
for(int j = 0; j < p; j++)
m[i*p + j] = rand();
float* tau = (float*) calloc(p, sizeof(float));
LAPACKE_sgeqrf(LAPACK_ROW_MAJOR, p, p, m, p, tau);
LAPACKE_sorgqr(LAPACK_ROW_MAJOR, p, p, 1, m, p, tau);
free(tau);
}

void out(int r, int c, float *m)
{
for (int i = 0; i < r; i++)
{
cout << "{ ";
for (int j = 0; j < c; j++)
cout << m[i*c + j] << " ";
cout << "}" << endl;
}
}

void choleskyFactorization(int p, float *m)
{
lapack_int q = LAPACKE_spotrf(LAPACK_ROW_MAJOR, 'L', p, m, p);
}

void copy(int p, float *from, float *to)
{
mkl_somatcopy('R', 'N', p, p, 1.0, from, p, to, p);
}

void samplesGenerate(int n, int p, float *mean, float *chol, float *res)
{
VSLStreamStatePtr stream;
int i = (int) rand();
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, n, res, p, VSL_MATRIX_STORAGE_FULL, mean, chol);
}

void wGenerate(int n, float* w)
{
float *mean = (float *) calloc(n, sizeof(float));
float *cov = (float*) calloc(n * n, sizeof(float));
for (int i = 1; i < n; i++)
{
mean[i] = 1;
cov[i*n + i] = 1;
}
int i = (int) rand();
VSLStreamStatePtr stream;
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, 1, w, n, VSL_MATRIX_STORAGE_FULL, mean, cov);
free(mean);
free(cov);
}

void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *W = (float *) calloc(n, sizeof(float));
float *X = (float *) calloc(p, sizeof(float));
wGenerate(n, W);
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, W[k], X, p, X, p, 1.0, Result, p);
}
free(W);
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}

free(X);
}

void sampleCovariance(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *X = (float *) calloc(p, sizeof(float));
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, Result, p);
}
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}
free(X);
}

void eigenDecomposition(int p, float* m, float *values) //decreasing orderб столбыц!
{
LAPACKE_ssyevd(LAPACK_ROW_MAJOR, 'V', 'U', p, m, p, values);
}

void rProjector(int p, int r, float *vectors, float* projector)
{
int num = p - r;
float *X = (float *) calloc(p, sizeof(float));
for (int i = 0; i < p; i++)
{
X[i] = vectors[p*i + num];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, projector, p);
free(X);
}

float S(int n, int p, float *m1, float *m2)
{
float *m = (float *) calloc(p * p, sizeof(float));
for (int i = 0; i < p; i++)
for (int j = 0; j < p; j++)
m[i*p + j] = m1[i*p + j] - m2[i*p + j];
float norm = LAPACKE_slange(LAPACK_ROW_MAJOR, 'F', p, p, m, p);
free(m);
return n*norm*norm;
}

int main(int argc, char **argv)
{
//задаем размерности
int paral = 16;
int p = 20;
int k = 3;
int M = 3000;
int r = 1; //нумерация с 1
int N = 3; //45
int n = 100;
int rank, size, prev, next;

//инициализация MPI
MPI_Request reqs[4];
MPI_Status stats[4];
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//
double MPI_Wtime(void);
double t1, t2, dt;
//
t1 = MPI_Wtime();
prev = rank - 1;
next = rank + 1;
if(rank==0) prev = size - 1;
if(rank==size - 1) next = 0;

//чтобы у всех сгенерировалась одинаковая матрица ковариации srand сработает одновременно
srand(time(NULL));
//файлы, в который будем писать результаты, у каждого параллельного процесса свой (номер процесса + .txt)
ofstream myfile;
char intStr[2];
if (rank < 10)
{
intStr[0] = '0';
intStr[1] = '0' + rank;
}
else
{
intStr[0] = '1';
intStr[1] = '0' + rank - 10;
}
char *txt = ".txt";
char *str = strcat(intStr,txt);

myfile.open(str);

//файлы, в которые будем отправлять размерности и тд, за это будет отвечать нулевой процесс
ofstream dim;
ofstream val;

if (rank == 0){
dim.open("dimensions.txt");
val.open("values.txt");
}

if (rank == 0)
{
dim << "p = " << p << endl;
dim << "k = " << k << endl;
dim << "M = " << M << endl;
dim << "r = " << r << endl;
dim << "N = " << N*(paral-1) << endl << endl;
dim << "n = " << n;
}

//генерация матрицы ковариации
float *COV = (float *) calloc(p * p, sizeof(float));
float *MEAN = (float *) calloc(p, sizeof(float));
float *COV_Chol = (float *) calloc(p * p, sizeof(float));
float *Basis = (float *) calloc(p * p, sizeof(float));
float *COV_L = (float *) calloc(p * p, sizeof(float));

ortBasis(p, Basis);

diagMatrix(p, k, COV_L);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, Basis, p, COV_L, p, 0.0, COV, p);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, COV, p, Basis, p, 0.0, COV_L, p);
copy(p, COV_L, COV);

copy(p, COV, COV_Chol);
choleskyFactorization(p, COV_Chol);
free(Basis);
free(COV_L);

float *values_pr = (float *) calloc(p, sizeof(float));
float *vectors_pr = (float *) calloc(p * p, sizeof(float));

copy(p, COV, vectors_pr);
eigenDecomposition(p, vectors_pr, values_pr);

//для исследования нам понадобятся собственные значения, поэтому их тоже выводит нулевой процесс
if(rank ==0){
for (int i = 0; i < p; i++)
{
val << values_pr[i] << " ";
}
val << endl << endl;
}

free(values_pr);
free(vectors_pr);

//вычислим проектор для реальной матрицы ковариации
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorR = (float *) calloc(p * p, sizeof(float));

copy(p, COV, vectors);
eigenDecomposition(p, vectors, values);
rProjector(p, r, vectors, ProjectorR);
//ProjectorR - нужный проектор
free(values);
free(vectors);

//для каждого процесса запустим свой рандом (зависящий от ранка), иначе все значения у каждого процесса будут одинаковые
srand(time(NULL) + rank);

float* NormsBootstrapWorld = (float *) calloc(N * M, sizeof(float));
for (int i = 0; i < N; i++)
{
float *Samples = (float *) calloc(n * p, sizeof(float)); //Каждый Xj - строка в матрице
for (int j = 0; j < M; j++)
{
float *COV_SampleBootstrap = (float *) calloc(p * p, sizeof(float));
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorSB = (float *) calloc(p * p, sizeof(float));

sampleCovarianceBootstrap(n, p, Samples, COV_SampleBootstrap);
copy(p, COV_SampleBootstrap, vectors);
eigenDecomposition(p, vectors, values); //вектора в столбцах
rProjector(p, r, vectors, ProjectorSB);
//ProjectorSB - нужный проектор

float s = S(n, p, ProjectorR, ProjectorSB);
NormsBootstrapWorld[i*M + j] = s;

free(COV_SampleBootstrap);
free(values);
free(vectors);
free(ProjectorSB);
}
free(Samples);
}

free(ProjectorR);

for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
myfile << NormsBootstrapWorld[j*M + i] << " ";
}
myfile << endl;
}
//в myfile N столбцов

myfile << endl;
free(NormsBootstrapWorld);

free(COV);
free(MEAN);
free(COV_Chol);
t2 = MPI_Wtime();
dt = t2 - t1;
if (rank == 0)
cout << dt << endl;
MPI_Finalize(); //конец параллеливания
}
</source>

= Литература =

* [https://en.wikipedia.org/wiki/Bootstrapping_(statistics)]
* [https://software.intel.com/sites/default/files/managed/ff/c8/mkl-2017-developer-reference-c_0.pdf]

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T11:27:54Z

Основные авторы описания: [[Участник:v|В.С.Шумовская]]

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

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

Пусть <math>X_{1},\dots, X_{n}</math> -- независимые и нормально распределенные случайные вектора в <math>\R^p</math> с нулевым средним и матрицей ковариацией <math>\Sigma</math>, она лежит в <math>\R^{pxp}</math> и такая, что ее собственные значения быстро убывают, т.е. 3-5 больших, а остальные, к примеру, в диапазоне от 1 до 3.

К этой выборке применим бутстрэп и вычислим в мире бутстрэпа <math>M</math> матриц ковариаций сигма <math>\Sigma^{o}_{j}, j = 1,\dots,M.</math>

Далее фиксируем <math>r</math>, обозначим за <math>P_{r}, P^{o}_{j}, j = 1,\dots,M</math> - проекторы на r-ое подпространство и вычислим следующие статистики:

<math>S^{o}_{j} = ||P_{r} - P^{o}_{r}||^{2}_{2}</math>

Задача -- вычислить большое количество этих статистик для визуализации их распределения.

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

#Генерация матрицы ковариации.
##Возьмем нужный нам набор собственных значений и поместим их на диагонали матрицы <math>\Lambda.</math>
##Далее генерируем базис <math>p</math>-мерного пространства и поместим его в матрицу <math>U.</math>
##Тогда <math>\Sigma = U\Lambda U^{T}.</math>
#Проделываем следующий цикл N раз:
##Генерация векторов <math>X</math> ~ <math>N(\theta, \Sigma).</math>
##Проделываем следующий цикл M раз:
###Генерация весов <math>w</math> ~ <math>N(1, 1).</math>
###Вычисление матриц ковариаций в бутстрэпе <math>\Sigma^{o} = \frac{1}{n}\sum_{k = 1}^{n}w_{k}X_{k}X^{T}_{k}. </math>
###Вычисление собственных векторов <math>u_{1},\dots,u_{p}</math> (в порядке убывания).
###Вычисление проекторов (считаем, что одномерные) по формуле <math>u = u_{r}u^{T}_{r}.</math>
###Вычисление статистик <math>S^{o} = ||P_{r} - P^{o}_{r}||^{2}_{2}.</math>

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

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

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

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

[[Файл:послед.png|center|500px]]

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

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

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

[[Файл:Парал.png|center|700px]]

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

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

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

* Размерности матриц
* Количество повторений циклов
* Собственные значения
Эти параметры меняются в коде

=== Выходные данные ===

От каждого вычислительного узла будет получена матрица значений статистик

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

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

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

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

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

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

Для исследования масштабируемости были проведены серии тестов с различным числом узлов (в степенях 2). Чем больше узлов, тем меньше загруженность каждого, если хотим получить одинаковый совокупный объем данных, а, следовательно, программа выполняется быстрее.
График для фиксированных n = 100, p = 10.
[[Файл:Graphhh.jpg|center]]

Все тесты проводились на [https://ru.wikipedia.org/wiki/Ломоносов_(суперкомпьютер) суперкомпьютере Ломоносов].
Характеристики вычислительного кластера: [http://users.parallel.ru/wiki/pages/22-config link].

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

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

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

<source lang="c++">
#include "mkl.h"
#include <ctime>
#include <cstdlib>
#include <iostream>
#include <string.h>
#include <sstream>
#include "mpi.h"
#include <fstream>
#include <algorithm>
#include <iterator>
#include <string>
#include <vector>

using namespace std;

void diagMatrix(int p, int k, float *m);
void ortBasis(int p, float *m);
void out(int p, float *m);
void choleskyFactorization(int p, float *m);
void copy(int p, float *from, float *to);
void samplesGenerate(int n, int p, float *mean, float *chol, float *res);
void wGenerate(int n, float* w);
void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result);
void sampleCovariance(int n, int p, float* Samples, float* Result);
void eigenDecomposition(int p, float* m, float *values);
void rProjector(int p, int r, float *vectors, float* projector);
float S(int n, int p, float *m1, float *m2);
float **alloc2d(int n, int m);

float **alloc2d(int n, int m)
{
float *data = (float*) malloc(n*m*sizeof(float));
float **array = (float **)malloc(n*sizeof(float *));
for (int i=0; i<n; i++)
array[i] = &(data[i*m]);
return array;
}

struct aboutMatrix
{
float* matrix;
float* l;
float* u;
float* cholesky;
};

void diagMatrix(int p, int k, float *m)
{
float numbers[5] = {36, 30, 25, 19, 15};
int ind = 0;
for (int i = 0; i < p; ++i)
{
for (int j = 0; j < p; ++j)
m[i*p + j] = 0.0;
if (k > 0)
{
m[i*p + i] = numbers[ind];
ind++;
k--;
}
else
m[i*p + i] = 1.0 + 5.0*rand() / float(RAND_MAX);
}
}

void ortBasis(int p, float *m)
{
for (int i = 0; i < p; i++)
for(int j = 0; j < p; j++)
m[i*p + j] = rand();
float* tau = (float*) calloc(p, sizeof(float));
LAPACKE_sgeqrf(LAPACK_ROW_MAJOR, p, p, m, p, tau);
LAPACKE_sorgqr(LAPACK_ROW_MAJOR, p, p, 1, m, p, tau);
free(tau);
}

void out(int r, int c, float *m)
{
for (int i = 0; i < r; i++)
{
cout << "{ ";
for (int j = 0; j < c; j++)
cout << m[i*c + j] << " ";
cout << "}" << endl;
}
}

void choleskyFactorization(int p, float *m)
{
lapack_int q = LAPACKE_spotrf(LAPACK_ROW_MAJOR, 'L', p, m, p);
}

void copy(int p, float *from, float *to)
{
mkl_somatcopy('R', 'N', p, p, 1.0, from, p, to, p);
}

void samplesGenerate(int n, int p, float *mean, float *chol, float *res)
{
VSLStreamStatePtr stream;
int i = (int) rand();
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, n, res, p, VSL_MATRIX_STORAGE_FULL, mean, chol);
}

void wGenerate(int n, float* w)
{
float *mean = (float *) calloc(n, sizeof(float));
float *cov = (float*) calloc(n * n, sizeof(float));
for (int i = 1; i < n; i++)
{
mean[i] = 1;
cov[i*n + i] = 1;
}
int i = (int) rand();
VSLStreamStatePtr stream;
vslNewStream(&stream,VSL_BRNG_MCG31, i);
vsRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, stream, 1, w, n, VSL_MATRIX_STORAGE_FULL, mean, cov);
free(mean);
free(cov);
}

void sampleCovarianceBootstrap(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *W = (float *) calloc(n, sizeof(float));
float *X = (float *) calloc(p, sizeof(float));
wGenerate(n, W);
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, W[k], X, p, X, p, 1.0, Result, p);
}
free(W);
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}

free(X);
}

void sampleCovariance(int n, int p, float* Samples, float* Result)
{
//Result p*p
float *X = (float *) calloc(p, sizeof(float));
for(int k = 0; k < n; k++)
{
for (int i = 0; i < p; i++)
{
X[i] = Samples[k*p + i];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, Result, p);
}
for(int k = 0; k < p; k++)
{
for (int i = 0; i < p; i++)
{
Result[k*p + i] /= n;
}
}
free(X);
}

void eigenDecomposition(int p, float* m, float *values) //decreasing orderб столбыц!
{
LAPACKE_ssyevd(LAPACK_ROW_MAJOR, 'V', 'U', p, m, p, values);
}

void rProjector(int p, int r, float *vectors, float* projector)
{
int num = p - r;
float *X = (float *) calloc(p, sizeof(float));
for (int i = 0; i < p; i++)
{
X[i] = vectors[p*i + num];
}
cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, p, p, 1, 1.0, X, p, X, p, 1.0, projector, p);
free(X);
}

float S(int n, int p, float *m1, float *m2)
{
float *m = (float *) calloc(p * p, sizeof(float));
for (int i = 0; i < p; i++)
for (int j = 0; j < p; j++)
m[i*p + j] = m1[i*p + j] - m2[i*p + j];
float norm = LAPACKE_slange(LAPACK_ROW_MAJOR, 'F', p, p, m, p);
free(m);
return n*norm*norm;
}

int main(int argc, char **argv)
{
//задаем размерности
int paral = 16;
int p = 20;
int k = 3;
int M = 3000;
int r = 1; //нумерация с 1
int N = 3; //45
int n = 100;
int rank, size, prev, next;

//инициализация MPI
MPI_Request reqs[4];
MPI_Status stats[4];
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//
double MPI_Wtime(void);
double t1, t2, dt;
//
t1 = MPI_Wtime();
prev = rank - 1;
next = rank + 1;
if(rank==0) prev = size - 1;
if(rank==size - 1) next = 0;

//чтобы у всех сгенерировалась одинаковая матрица ковариации srand сработает одновременно
srand(time(NULL));
//файлы, в который будем писать результаты, у каждого параллельного процесса свой (номер процесса + .txt)
ofstream myfile;
char intStr[2];
if (rank < 10)
{
intStr[0] = '0';
intStr[1] = '0' + rank;
}
else
{
intStr[0] = '1';
intStr[1] = '0' + rank - 10;
}
char *txt = ".txt";
char *str = strcat(intStr,txt);

myfile.open(str);

//файлы, в которые будем отправлять размерности и тд, за это будет отвечать нулевой процесс
ofstream dim;
ofstream val;

if (rank == 0){
dim.open("dimensions.txt");
val.open("values.txt");
}

if (rank == 0)
{
dim << "p = " << p << endl;
dim << "k = " << k << endl;
dim << "M = " << M << endl;
dim << "r = " << r << endl;
dim << "N = " << N*(paral-1) << endl << endl;
dim << "n = " << n;
}

//генерация матрицы ковариации
float *COV = (float *) calloc(p * p, sizeof(float));
float *MEAN = (float *) calloc(p, sizeof(float));
float *COV_Chol = (float *) calloc(p * p, sizeof(float));
float *Basis = (float *) calloc(p * p, sizeof(float));
float *COV_L = (float *) calloc(p * p, sizeof(float));

ortBasis(p, Basis);

diagMatrix(p, k, COV_L);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, Basis, p, COV_L, p, 0.0, COV, p);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, p, p, p, 1.0, COV, p, Basis, p, 0.0, COV_L, p);
copy(p, COV_L, COV);

copy(p, COV, COV_Chol);
choleskyFactorization(p, COV_Chol);
free(Basis);
free(COV_L);

float *values_pr = (float *) calloc(p, sizeof(float));
float *vectors_pr = (float *) calloc(p * p, sizeof(float));

copy(p, COV, vectors_pr);
eigenDecomposition(p, vectors_pr, values_pr);

//для исследования нам понадобятся собственные значения, поэтому их тоже выводит нулевой процесс
if(rank ==0){
for (int i = 0; i < p; i++)
{
val << values_pr[i] << " ";
}
val << endl << endl;
}

free(values_pr);
free(vectors_pr);

//вычислим проектор для реальной матрицы ковариации
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorR = (float *) calloc(p * p, sizeof(float));

copy(p, COV, vectors);
eigenDecomposition(p, vectors, values);
rProjector(p, r, vectors, ProjectorR);
//ProjectorR - нужный проектор
free(values);
free(vectors);

//для каждого процесса запустим свой рандом (зависящий от ранка), иначе все значения у каждого процесса будут одинаковые
srand(time(NULL) + rank);

float* NormsBootstrapWorld = (float *) calloc(N * M, sizeof(float));
for (int i = 0; i < N; i++)
{
float *Samples = (float *) calloc(n * p, sizeof(float)); //Каждый Xj - строка в матрице
for (int j = 0; j < M; j++)
{
float *COV_SampleBootstrap = (float *) calloc(p * p, sizeof(float));
float *values = (float *) calloc(p, sizeof(float));
float *vectors = (float *) calloc(p * p, sizeof(float));
float *ProjectorSB = (float *) calloc(p * p, sizeof(float));

sampleCovarianceBootstrap(n, p, Samples, COV_SampleBootstrap);
copy(p, COV_SampleBootstrap, vectors);
eigenDecomposition(p, vectors, values); //вектора в столбцах
rProjector(p, r, vectors, ProjectorSB);
//ProjectorSB - нужный проектор

float s = S(n, p, ProjectorR, ProjectorSB);
NormsBootstrapWorld[i*M + j] = s;

free(COV_SampleBootstrap);
free(values);
free(vectors);
free(ProjectorSB);
}
free(Samples);
}

free(ProjectorR);

for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
myfile << NormsBootstrapWorld[j*M + i] << " ";
}
myfile << endl;
}
//в myfile N столбцов

myfile << endl;
free(NormsBootstrapWorld);

free(COV);
free(MEAN);
free(COV_Chol);
t2 = MPI_Wtime();
dt = t2 - t1;
if (rank == 0)
cout << dt << endl;
MPI_Finalize(); //конец параллеливания
}
</source>

= Литература =

ываыва

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T11:25:09Z

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T11:24:43Z

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T11:24:25Z

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T11:24:11Z

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T11:23:38Z

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-30T08:54:37Z

Участница:V/Вычисление статистик квадрата норм разностей спектральных проекторов случайных матриц

2016-11-29T22:30:26Z