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

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

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

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

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

//TODO

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

//TODO

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

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

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

// TODO about heisenberg

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

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

1.1. Метод Гивенса (вращений) QR-разложения квадратной матрицы или 1.2. Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант

2. Перемножение плотных неособенных матриц

// TODO about heisenberg

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

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

На итерации [math]\displaystyle{ n }[/math]:

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

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

// TODO about heisenberg form and algorithm with shift

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

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

// TODO about heisenberg complexity

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

//TODO

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

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

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

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

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

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

// TODO for Heisenberg

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

int 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

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

// не нужно???