Участник:Elijah/Нахождение собственных чисел квадратной матрицы методом QR разложения

Нахождение собственных чисел квадратной матрицы методом QR разложения
Последовательный алгоритм
Последовательная сложность	[math]N * O(n^3)[/math]
Объём входных данных	[math]n^2[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math] N * (12n - 15) [/math]
Ширина ярусно-параллельной формы	[math] O(n^2) [/math]

Основные авторы описания: И.В.Афанасьев, В.А.Шишватов

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

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

//TODO

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

//TODO

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

У базового QR алгоритма есть два вычислительных ядра:

1. операция поиска QR разложения с использованием метода Гивенса
   
2. матричное умножения полученных матриц

// TODO about heisenberg

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

Как уже описано в описании ядра алгоритма, базовая версия QR-алгоритма на каждой итерации использует следующие алгоритмы:

1. Метод Гивенса (вращений) QR-разложения квадратной матрицы или Метод Хаусхолдера(отражений) QR-разложения_квадратной_матрицы
2. Перемножение плотных неособенных матриц

// TODO about heisenberg

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

Последовательная реализация простейшего алгоритма сострит из некоторого числа итераций.

На итерации [math] n [/math]:

1. Для матрицы [math] A_{n} [/math] строится QR разложение любым доступным последовательным алгоритмом на матрицы [math] Q_{n} [/math] и [math] R_{n} [/math]
   
2. Матрицы [math] R_{n} [/math] и [math] Q_{n} [/math] перемножаются, таким образом получается матрица [math] A_{n+1} = R_{n} * Q_{n} [/math] для следующей итерации [math] n + 1 [/math]

По окончании каждой итерации проверятся, приведена ли матрица к диагональной форме, и если да, то производится выход из цикла.

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

На итерации [math] n [/math]:

1. Для матрицы [math] A_{n} [/math] строится QR разложение с помощью метода вращения Гивенса. При этом занулять нужно не все элементы, а только [math] a_{i + 1, i} [/math], где [math]i \in {1, n - 1}[/math]. При этом полученная матрица будет опять являться матрицей Хессенберга.

По окончании каждой итерации проверятся, приведена ли матрица к диагональной форме, и если да, то производится выход из цикла.

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

Сложность базового алгоритма:
Рассчитаем последовательную сложность базового алгоритма. Пусть [math] A \in \mathbb{R}^{n \times n}[/math], и для приведения матрицы к диагональной форме необходимо произвести N итераций алгоритма.
На каждой итерации алгоритма производится QR разложение (сложностью [math] 2 * n^3 [/math]) и матричное умножение (сложностью [math] n^3 [/math]). Проверка того, имеет ли матрица диагональную форму, может быть проведена за [math] n^2 [/math] операций.
Таким образом итоговая сложность одной итерации составляет: [math] 3*n^3 + n^2 = O(n^3)[/math]
Общая сложность алгоритма при [math] N [/math] итерациях составляет [math] N * O(n^3) [/math]

Сложность алгоритма с приведением матрицы к форме Гейзенберга:
// TODO about heisenberg complexity

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

Макрограф алгоритма изображён на рисунке 1. Он представляет из себя верхний уровень алгоритма, на котором в цикле итеративно вызываются следующие операции:

1. Инициализация матриц (init matrices)
2. QR разложение (QR decomposition)
3. Матричное умножение (matrix multiplication)
4. Проверка, приведена ли матрица к верхней треугольной форме (check if matrix is upper triangular)
5. Получение собственных значений из диагонали матрицы (get eigenvalues from matrix diagonal)

Информационные графы QR разложения и матричного умножения могут быть найдены по приведенным ссылкам. Информационный граф проверки того, приведена ли матрица к верхнетреугольной форме, будет приведен далее на рисунке 2. Информационные графы инициализации и получения собственных значений представляют собой работу с вводом и выводом входных/выходных данных, и, поэтому, рассматриваться не будут.

Рисунок 1. Информационный граф базовой версии QR-алгоритма без отображения входных и выходных данных. N итераций.

Рисунок 2. Информационный граф проверки приведенности матрицы к верхнетреугольной форме.

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

Все итерации базового QR-алгоритма производятся последовательно, поэтому на верхнем уровне алгоритм чисто последователен.

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

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

1. Параллельная сложность QR разложения методом вращений составляет [math] 11n - 16 [/math]
   
2. Параллельная сложность матричного умножения составляет [math] n [/math].
   
3. Параллельная сложность проверки, является ли матрица диагональной, равна 1.
   

Таким образом параллельная сложность каждой итерации составляет [math] 12n - 15 [/math]. При N итерациях, которые необходимо производить последовательно, итоговая сложность базового алгоритма составляет [math] N * (12n - 15) [/math]

Параллельная сложность алгоритма с приведением матрицы к форме Гейзенберга:
// TODO for Heisenberg

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

• Входные данные: плотная квадратная матрица [math]A[/math] (элементы [math]a_{ij}[/math]).

• Объём входных данных: [math]n^2[/math].

• Выходные данные: [math]n[/math] вещественных собственных чисел [math] | l_{i} | [/math] матрицы [math]A[/math]

• Объём выходных данных: [math]n[/math].

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

Соотношение последовательной и параллельной сложности как для базового, так и для оптимизированного алгоритма, является квадратичным, что даёт хороший стимул для распараллеливания.

Вычислительная мощность, равная отношению числа операций [math] N * O(n^3) [/math] к суммарному объему входных и выходных данных [math] n^2 + n [/math], для каждой итерации линейна, а для всего алгоритма в целом равна [math] N*n [/math], что позволяет сделать вывод о том, что перемещение данных для обработки не играет важной роли в данном алгоритме.

Вычислительная погрешность растет линейно, из-за использования метода вращений для QR разложения.

Сбалансированность (??????)

Данный алгоритм недетерминирован, так как число итераций зависит от значений элементов матрицы и выход из основного цикла происходит по достижению некоторой точности. В зависимости от этого, вычислительный поток имеет различную длину. На каждой отдельной итерации в свою очередь алгоритм полностью детерминирован, так как детерминированы его обе базовые макрооперации - QR разложение и матричное умножение.

Степнь исхода вершины информационного графа (??????)

"Длинные" дуги в информационном графе (??????)

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

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

Основным типом данных для QR алгоритма является матрица. В данной статье представлена реализация на языке C++, что накладывает определенные ограничения на формат хранени матриц: все матрицы хранятся в одномерном массиве в построчном (ROW_MAJOR) формате.

В простейшем случае QR-алгоритм выглядит следующим образом:

int sequential_QR_algorithm(float *A, int size)
{
    float *Q = new double[size * size];
    float *R = new double[size * size];
    while(!check_if_upper_triangular(A, size)))
    {
        sgeqrf (ROW_MAJOR, size, size, A, size, tau);
        set_Q_and_R(A, Q, R, size);
        sorgqr(ROW_MAJOR, size, size, size, Q, size, tau);
        sgemm (RowMajor, NoTrans, NoTrans, size, size, size, 1, R, size, Q, size, 0, A, size);
    }
}

В последовательном случае функции dgeqrf, dorgqr и dgemm представляют собой последовательные функции стандарта BLAS, вычисляющие QR разложение(geqrf), матрицу Q после разложения (dorgqr) и матричное умножение (dgemm). Так же здесь используются простейшие функции set_Q_and_R и check_if_upper_triangular. Функция set_Q_and_R получает после QR разложения из матрицы A данные о матрицах Q и R. Функция check_if_upper_triangular проверяет, приведена ли матрица A к треугольном виду, и если да, происходит выход из цикла.

Для работы программы в памяти необходимо хранить 3 матрицы, таким образом общий объем составляет 3 * n^2. Так же необходим вспомогательный массив tau для QR разложения, однако его размер составляет всего n. Так же важно заметить, что различные реализации функций dgeqrf, dorgqr и dgemm обычно используют внутренние вспомогательные массивы для блочной обработки, однако, их размер не превышает n * C, где C - некоторая архитектурно зависимая константа.

Все вычисления рекомендуется производить в одинарной точности.

Как уже говорилось ранее, внешний цикл представленной программы сугубо последователен (итеративен), и основной ресурс параллелизма расположен в вызываемых функциях (geqrf, gemm). Их реализации бывают крайне различны, поэтому данный алгоритм можно легко преобразовать для любой архитектуры, используя набор стандартных библиотечных функций. К примеру, для реализации алгоритма на GPU можно использовать библиотеку MAGMA, на многоядерных CPU - Intel MKL.

// TODO code with heisinberg

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

Как уже говорилось в разделе 2.1., при использовании различных библиотек можно реализовать данный алгоритм практически на всех архитектурах. Далее более подробно будут рассмотрены примеры для некоторых из них. Особенности параллельных реализаций для базовых макроопераций (QR разложения и матричного умножения), уже описаны в соответствующих статьях (см. ссылки в разделе 1.4 Макроструктура алгоритма)

Реализация на многоядерных CPU (одном узле)

Для распарарллеливания на многоядерных CPU используются функции LAPACKE_sgeqrf, LAPACKE_sorgqr и cblas_sgemm. Распараллеливание оставшихся функций может быть произведено с использованием openMP, так как внутренние циклы в них независимы по данным.

Исходный код предложенной реализации представлен далее:

int CPU_QR_algorithm(float *A, int size)
{
    float *Q = new double[size * size];
    float *R = new double[size * size];
    while(!check_if_upper_triangular(A, size)))
    {
        LAPACKE_sgeqrf (LAPACK_ROW_MAJOR, size, size, A, size, tau);
        set_Q_and_R(A, Q, R, size);
        LAPACKE_sorgqr(LAPACK_ROW_MAJOR, size, size, size, Q, size, tau);
        cblas_sgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans, size, size, size, 1, R, size, Q, size, 0, A, size);
    }
}

Реализация на GPU (и multi-GPU)

Для реализации алгоритма на GPU и multi-GPU можно использовать аналогичные функции dgemm и geqrf из библиотеки magma.