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

Основные авторы описания: В.С.Шумовская

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

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

Пусть [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]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Эти параметры меняются в коде

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

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

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

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

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

Все тесты проводились на суперкомпьютере Ломоносов. Характеристики вычислительного кластера: link.

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

#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(); //конец параллеливания
}

3 Литература

• Boostrapping (wiki)
• Developer Reference for Intel® Math Kernel Library 2017 - C(pdf)