Эта работа успешно выполнена
Преподавателю: в основное пространство, в подстраницу
Данное задание было проверено и зачтено.
Проверено Frolov и VadimVV.

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

Основные авторы описания: О.Д.Баев (пп. 1.1 - 1.5 базовый, 1.9, 2.4, 2.7), А.С.Шевелев (пп. 1.2-1.5 хессенберг, 1.6 - 1.8, 1.10)

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

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

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

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

Базовый QR-алгоритм

Пусть [math]\displaystyle{ A }[/math] — вещественная квадратная матрица, для которой мы хотим найти собственные числа и векторы. Положим [math]\displaystyle{ A_0 = A }[/math]. На k-ом шаге (начиная с k = 0) вычислим QR-разложение [math]\displaystyle{ A_k = Q_k R_k }[/math], где [math]\displaystyle{ Q_k }[/math] — ортогональная матрица (то есть [math]\displaystyle{ Q_k^{T} = Q_k^{-1} }[/math]), а [math]\displaystyle{ R_k }[/math] — верхняя треугольная матрица. Затем мы определяем [math]\displaystyle{ A_{k+1} = R_k Q_k }[/math].

Заметим, что

[math]\displaystyle{ A_{k+1} = R_k Q_k = Q_k^{-1} Q_k R_k Q_k = Q_k^{-1} A_k Q_k = Q_k^{T} A_k Q_k, }[/math]

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

Пусть все диагональные миноры матрицы [math]\displaystyle{ A }[/math] не вырождены. Тогда последовательность матриц [math]\displaystyle{ A_k }[/math] при [math]\displaystyle{ k \rightarrow \infty }[/math] сходится по форме к клеточному правому треугольному виду, соответствующему клеткам с одинаковыми по модулю собственными значениями. ^[1]

Для того, чтобы получить собственные векторы матрицы, нужно перемножить все матрицы [math]\displaystyle{ Q_k }[/math].

Приведение матрицы к форме Хессенберга

Матрица [math]\displaystyle{ A }[/math] имеет форму Хессенберга^[2], если все элементы [math]\displaystyle{ a_{ij} }[/math] под первой поддиагональю равны нулю ([math]\displaystyle{ a_{ij} = 0, i \gt j + 1 }[/math]).

Любая матрица [math]\displaystyle{ A }[/math] может быть приведена к подобной матрице, имеющей форму Хессенберга. Подобная манипуляция позволяет значительно ускорить дальнейший процесс QR-разложения. Одним из способов приведения матрицы к форме Хессенберга является метод Хаусхолдера.

Идея метода состоит в последовательном приведении исходной матрицы [math]\displaystyle{ A }[/math] к форме Хессенберга за счет умножения на матрицы отражений, каждая из которых может быть определена одним вектором и имеет вид [math]\displaystyle{ U = E - 2ww^* }[/math], где [math]\displaystyle{ w }[/math] - вектор, удовлетворяющий равенству [math]\displaystyle{ w^*w = 1 }[/math]. Цель домножения исходной матрицы на матрицу отражения на [math]\displaystyle{ i }[/math]-ом шаге заключается в обнулении поддиагональных элементов [math]\displaystyle{ i }[/math]-го столбца. Таким образом за [math]\displaystyle{ n-1 }[/math] итерацию исходная матрица будет приведена к форме Хессенберга.

При этом матрица отражений на каждом из шагов имеет не стандартный вид, описанный ранее, а [math]\displaystyle{ U=E-\frac{1}{\gamma}vv^* }[/math], где [math]\displaystyle{ v }[/math] находится для левых матриц отражения через координаты текущего [math]\displaystyle{ i }[/math]-го столбца при помощи операции скалярного произведения.

QR-алгоритм со сдвигами

Скорость сходимости QR-алгоритма зависит от того, насколько близки друг к другу собственные числа матрицы [math]\displaystyle{ A }[/math] - чем ближе они друг к другу, тем медленнее алгоритм сходится. Поэтому для ускорения сходимости используют так называемые сдвиги.

Суть метода заключается в следующем: Пусть на [math]\displaystyle{ k }[/math]-й итерации есть матрица [math]\displaystyle{ A_k }[/math], тогда переход к матрице [math]\displaystyle{ A_{k+1} }[/math] происходит следующим образом:

Выбирается число [math]\displaystyle{ v_k }[/math] и выполняется QR-разложение: [math]\displaystyle{ A_k - v_kE = Q_kR_k }[/math].
[math]\displaystyle{ A_{k+1} = R_kQ_k + v_kE }[/math]

При этом матрицы [math]\displaystyle{ A_k }[/math] и [math]\displaystyle{ A_{k+1} }[/math] подобны, поэтому набор собственных чисел не изменяется.

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

Базовый QR-алгоритм

Основными вычислительными ядрами алгоритма являются:

QR-разложение матрицы [math]\displaystyle{ A_k = Q_k R_k }[/math]
Умножение плотных матриц [math]\displaystyle{ A_{k+1} = R_k Q_k }[/math]

Существуют различные алгоритмы вычисления QR-разложения. Можно выделить такие, как метод Гивенса и метод Хаусхолдера.

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

QR-алгоритм с приведением к форме Хессенберга и сдвигами

Основными вычислительными ядрами данного варианта алгоритма являются:

Приведение исходной матрицы [math]\displaystyle{ A }[/math] к матрице в форме Хессенберга. Данная операция производится методом Хаусхолдера.
Выполнение QR-разложения матрицы [math]\displaystyle{ A_k - v_kE = Q_kR_k }[/math] модифицированным методом Гивенса.
Перемножение двух плотных матриц [math]\displaystyle{ R_kQ_k }[/math].

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

Базовый QR-алгоритм

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

QR-разложение
- Метод Гивенса (вращений)
- Метод Хаусхолдера (отражений)
Умножение матриц

QR-алгоритм с приведением к форме Хессенберга и сдвигами

В данном варианте макроструктура алгоритма выглядит следующим образом:

Приведение исходной матрицы к матрице в форме Хессенберга.
Итеративное выполнение QR-разложения методом Гивенса^[3] и перемножения матриц.

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

Базовый QR-алгоритм

QR-алгоритм является итерационным.

На каждой итерации [math]\displaystyle{ k }[/math]:

Строится QR-разложение матрицы [math]\displaystyle{ A_k = Q_k R_k }[/math]
Получается матрица [math]\displaystyle{ A_{k+1} = R_k Q_k }[/math], используемая на следующей итерации

Алгоритм выполняется до сходимости матрицы [math]\displaystyle{ A_k }[/math] к треугольному виду.

Псевдокод алгоритма:

algorithm QR is
  input: Matrix A
  output: Eigenvalues of the matrix A
    
  repeat
    Q, R ← QR_factorization(A)
    A ← R×Q
  until convergence
    
  return diag(A)

QR-алгоритм с приведением к форме Хессенберга и сдвигами

Как уже было описано ранее, данный алгоритм состоит из следуюших этапов:

Приведение исходной матрицы [math]\displaystyle{ A }[/math] к матрице в форме Хессенберга [math]\displaystyle{ H }[/math].
Выполнение [math]\displaystyle{ k }[/math] итераций ([math]\displaystyle{ A_0 = H }[/math]):
- Выбор коэффициента [math]\displaystyle{ v_k }[/math].
- QR-разложение [math]\displaystyle{ A_k - v_kE = Q_kR_k }[/math].
- [math]\displaystyle{ A_{k+1} = R_kQ_k + v_kE }[/math].

Псевдокод алгоритма:

algorithm QR is
  input: Matrix A
  output: Eigenvalues of the matrix A
  
  H ← To_Hessenberg(A)
  A ← H
  repeat
    v_k ← Get_coeff(A)
    Q, R ← QR_factorization(A - v_k*E)
    A ← R×Q + v_k*E 
  until convergence
    
  return diag(A)

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

Базовый QR-алгоритм

Пусть дана матрица [math]\displaystyle{ A \in \mathbb{R}^{n \times n} }[/math]. До момента приведения матрицы к треугольной форме произведено [math]\displaystyle{ N }[/math] итераций алгоритма. На каждой итерации алгоритма производится три действия: QR-разложение, перемножение матриц и проверка, является ли матрица треугольной.

Распишем последовательную сложность каждого действия:

QR-разложение:
- Методом Гивенса (вращений) - [math]\displaystyle{ 2*n^3 }[/math]^[3].
- Методом Хаусхолдера (отражений) - [math]\displaystyle{ \frac{4}{3}n^3 }[/math]^[4].
Перемножение матриц - [math]\displaystyle{ n^3 }[/math] операций^[5].
Проверка, имеет ли матрица треугольную форму - [math]\displaystyle{ \frac{n^2 - n}{2} }[/math] операций.

Таким образом, на каждой итерации последовательный алгоритм выполняет [math]\displaystyle{ O(n^3) }[/math] операций. Для всего алгоритма потребуется выполнить [math]\displaystyle{ N * O(n^3) }[/math] операций.

QR-алгоритм с приведением к форме Хессенберга и сдвигами

Пусть как и для базового алгоритма дана матрица [math]\displaystyle{ A \in \mathbb{R}^{n \times n} }[/math]. Как уже было описано выше, в данном варианте алгоритма выполняются три основные операции: приведение исходной матрицы к матрице в форме Хессенберга, QR-разложение методом Гивенса и перемножение матриц.

Распишем последовательную сложность каждой операции:

Приведение к матрице в форме Хессенберга методом Хаусхолдера потребует [math]\displaystyle{ O(n^3) }[/math] операций.
QR-разложение методом Гивенса матрицы в форме Хессенберга - [math]\displaystyle{ 6n^2 + O(n) }[/math] операций ^[6].
Перемножение матриц потребует столько же операций, как и для базового алгоритма - [math]\displaystyle{ n^3 }[/math] операций.
Вычисление сдвигов не является затратной операцией по сравнению с другими.

Таким образом, полная последовательная сложность QR-алгоритма с матрицей в форме Хессенберга и сдвигами - [math]\displaystyle{ O(n^3) + N*(6n^2+O(n)) }[/math] или [math]\displaystyle{ O(n^3) + N*O(n^2) }[/math]

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

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

В случае QR-алгоритма, в котором применяется приведение исходной матрицы к форме Хессенберга и сдвиги, представленные графы будут немного изменены: на Рисунке 1 перед первой итерацией появится элемент с приведением матрицы, а на Рисунке 2 добавятся элементы, отображающие сдвиги.

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

Базовый QR-алгоритм

Базовый алгоритм является итерационным и выполняется последовательно. Возможность распараллеливания алгоритма предоставляется при реализации операций таких, как QR-разложение, перемножение матриц и проверка, имеет ли матрица треугольную форму.

Рассчитаем параллельную сложность для каждой из указанных операций:

QR-разложение:
- Методом Гивенса (вращений) - [math]\displaystyle{ 11n - 16 }[/math]^[3].
- Методом Хаусхолдера (отражений) - [math]\displaystyle{ O(n^2) }[/math]^[4].
Перемножение матриц - [math]\displaystyle{ n }[/math]^[5].
Проверка, имеет ли матрица треугольную форму - [math]\displaystyle{ 1 }[/math].

Результирующая параллельная сложность одной итерации [math]\displaystyle{ 12n-15 }[/math] (для метода вращений). А всего алгоритма из [math]\displaystyle{ N }[/math] итераций - [math]\displaystyle{ N*(12n-15) }[/math].

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

Приведением исходной матрицы к форме Хессенберга - [math]\displaystyle{ O(n^2) }[/math].
Отдельно каждая итерация - [math]\displaystyle{ O(n) }[/math].

Итоговая параллельная сложность алгоритма - [math]\displaystyle{ O(n^2) + N*O(n) }[/math].

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

Входные данные: плотная квадратная матрица [math]\displaystyle{ A }[/math] размера [math]\displaystyle{ n }[/math].

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

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

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

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

Описанный алгоритм обладает следующими свойствами:

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

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

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

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

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

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

Используется процедура PDHSEQR из библиотеки ScaLAPACK.
Данный метод работает с матрицами в форме Хессенберга, поэтому сначала происходит приведение матрицы к данному виду с помощью процедуры PDGEHRD

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

Время работы алгоритма (сек)
Процессы\Матрица	600	1200	1800	2400
1	4.30549	25.2566	78.0507	177.054
4	3.63899	17.6559	51.2146	127.386
9	3.81662	9.19687	23.9181	53.7105
16	3.11883	7.26239	16.0958	32.3561

Эффективность распараллеливания [math]\displaystyle{ E = S/p }[/math], где [math]\displaystyle{ S = T_1/T_p }[/math] - полученное ускорение работы программы, а [math]\displaystyle{ p }[/math] - число используемых процессов.

Эффективность распараллеливания (%)
Процессы\Матрица	600	1200	1800	2400
4	29.57	35.76	38.09	34.74
9	12.53	30.51	36.25	36.62
16	8.62	21.73	30.30	34.20

Оценка масштабируемости

Минимальная эффективность: 8.62%
Максимальная эффективность: 38.09%
По числу процессов: при увеличении числа процессов эффективность в целом уменьшается
По размеру задачи: при увеличении размера задачи эффективность в целом увеличивается

Уменьшение или увеличение эффективности происходит не интенсивно, существенные изменения заметны только при небольшом размере задачи.

Библиотека ScaLAPACK 2.0.0 была установлена самостоятельно. Использовался компилятор mpicxx с опциями -lscalapack -llapack и указанием пути до установленной библиотеки.

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

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

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

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

Полный алгоритм:

LAPACK: функция DGEEV (последовательная реализация)
ScaLAPACK: функция PDHSEQR (параллельная реализация)
ALGLIB: реализация алгоритма для различных языков программирования
Python NumPy: функции linalg.eig, linalg.eigvals
C++ Eigen: Eigenvalues модуль

QR-разложение:

LAPACK: функция DGEQRF (последовательная реализация)
ScaLAPACK: функция PDGEQRF (параллельная реализация)
MATLAB: функция qr
Mathematica: функция QRDecomposition
ALGLIB: реализация QR-разложения для различных языков программирования
Python NumPy: функция linalg.qr
C++ Eigen: QR модуль

Умножение матриц:

LAPACK: функция DGEMM (последовательная реализация)
ScaLAPACK: функция PDGEMM (параллельная реализация)
MATLAB: функции mtimes, *
Mathematica: функция Dot
Python NumPy: функция numpy.dot

Литература

↑ Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы" — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.
↑ | Матрица Хессенберга
↑ ^3,0 ^3,1 ^3,2 Метод Гивенса (вращений) QR-разложения квадратной матрицы
↑ ^4,0 ^4,1 Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы
↑ ^5,0 ^5,1 Перемножение плотных неособенных матриц
↑ | QR algorithm

[1] Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы" — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.

[2] | Матрица Хессенберга

[QRGivens-3] 3,0 ^3,1 ^3,2 Метод Гивенса (вращений) QR-разложения квадратной матрицы

[QRHouseholder-4] 4,0 ^4,1 Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы

[MatrixMultiplication-5] 5,0 ^5,1 Перемножение плотных неособенных матриц

[6] | QR algorithm

[1]

[2]

[3]

[4]

[5]

[6]

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

Содержание

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

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

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

Базовый QR-алгоритм

Приведение матрицы к форме Хессенберга

QR-алгоритм со сдвигами

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

Базовый QR-алгоритм

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

Базовый QR-алгоритм

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

Базовый QR-алгоритм

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

Базовый QR-алгоритм

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

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

Базовый QR-алгоритм

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

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

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

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

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

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

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

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

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

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

Литература

Навигация

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

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

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

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

Базовый QR-алгоритм

Приведение матрицы к форме Хессенберга

QR-алгоритм со сдвигами

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

Базовый QR-алгоритм

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

Базовый QR-алгоритм

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

Базовый QR-алгоритм

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

Базовый QR-алгоритм

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

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

Базовый QR-алгоритм

QR-алгоритм с приведением к форме Хессенберга и сдвигами

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

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

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

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

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

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

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

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

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

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

Литература

Навигация

Поиск