Участник:Odbaev/Нахождение собственных чисел квадратной матрицы методом QR разложения (2): различия между версиями
Frolov (обсуждение | вклад) |
|||
(не показано 89 промежуточных версий 3 участников) | |||
Строка 1: | Строка 1: | ||
+ | {{Assignment|VadimVV|Frolov}} | ||
+ | |||
{{algorithm | {{algorithm | ||
| name = Нахождение собственных чисел квадратной матрицы методом QR разложения | | name = Нахождение собственных чисел квадратной матрицы методом QR разложения | ||
− | | serial_complexity = <math>N * O(n^ | + | | serial_complexity = <math>O(n^3) + N * O(n^2)</math> |
− | | pf_height = <math>N*( | + | | pf_height = <math>O(n^2)+N*O(n)</math> |
| pf_width = <math>O(n^2)</math> | | pf_width = <math>O(n^2)</math> | ||
| input_data = <math>n^2</math> | | input_data = <math>n^2</math> | ||
Строка 8: | Строка 10: | ||
}} | }} | ||
− | Основные авторы описания: [[Участник:Odbaev|О.Д.Баев]], [[Участник:AlexShevelev|А.С.Шевелев]] | + | Основные авторы описания: [[Участник:Odbaev|О.Д.Баев]] (пп. 1.1 - 1.5 базовый, 1.9, 2.4, 2.7), [[Участник:AlexShevelev|А.С.Шевелев]] (пп. 1.2-1.5 хессенберг, 1.6 - 1.8, 1.10) |
= Свойства и структура алгоритма = | = Свойства и структура алгоритма = | ||
Строка 16: | Строка 18: | ||
== Математическое описание алгоритма == | == Математическое описание алгоритма == | ||
+ | |||
+ | === Базовый QR-алгоритм === | ||
Пусть <math>A</math> — вещественная квадратная матрица, для которой мы хотим найти собственные числа и векторы. Положим <math>A_0 = A</math>. На ''k''-ом шаге (начиная с ''k'' = 0) вычислим QR-разложение <math>A_k = Q_k R_k</math>, где <math>Q_k</math> — ортогональная матрица (то есть <math>Q_k^{T} = Q_k^{-1}</math>), а <math>R_k</math> — верхняя треугольная матрица. Затем мы определяем <math>A_{k+1} = R_k Q_k</math>. | Пусть <math>A</math> — вещественная квадратная матрица, для которой мы хотим найти собственные числа и векторы. Положим <math>A_0 = A</math>. На ''k''-ом шаге (начиная с ''k'' = 0) вычислим QR-разложение <math>A_k = Q_k R_k</math>, где <math>Q_k</math> — ортогональная матрица (то есть <math>Q_k^{T} = Q_k^{-1}</math>), а <math>R_k</math> — верхняя треугольная матрица. Затем мы определяем <math>A_{k+1} = R_k Q_k</math>. | ||
Строка 22: | Строка 26: | ||
то есть все матрицы <math>A_k</math> являются подобными и их собственные значения равны. | то есть все матрицы <math>A_k</math> являются подобными и их собственные значения равны. | ||
− | Пусть все диагональные миноры матрицы <math>A</math> не вырождены. Тогда последовательность матриц <math>A_k</math> при <math>k \rightarrow \infty</math> сходится по форме к клеточному правому треугольному виду, соответствующему клеткам с одинаковыми по модулю собственными значениями. <ref> Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с. </ref> | + | Пусть все диагональные миноры матрицы <math>A</math> не вырождены. Тогда последовательность матриц <math>A_k</math> при <math>k \rightarrow \infty</math> сходится по форме к клеточному правому треугольному виду, соответствующему клеткам с одинаковыми по модулю собственными значениями. <ref> Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы" — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с. </ref> |
Для того, чтобы получить собственные векторы матрицы, нужно перемножить все матрицы <math>Q_k</math>. | Для того, чтобы получить собственные векторы матрицы, нужно перемножить все матрицы <math>Q_k</math>. | ||
+ | |||
+ | === Приведение матрицы к форме Хессенберга === | ||
+ | Матрица <math>A</math> имеет форму Хессенберга<ref>[https://en.wikipedia.org/wiki/Hessenberg_matrix | Матрица Хессенберга]</ref>, если все элементы <math>a_{ij}</math> под первой поддиагональю равны нулю (<math>a_{ij} = 0, i > j + 1</math>). | ||
+ | |||
+ | Любая матрица <math>A</math> может быть приведена к подобной матрице, имеющей форму Хессенберга. Подобная манипуляция позволяет значительно ускорить дальнейший процесс QR-разложения. Одним из способов приведения матрицы к форме Хессенберга является метод Хаусхолдера. | ||
+ | |||
+ | Идея метода состоит в последовательном приведении исходной матрицы <math>A</math> к форме Хессенберга за счет умножения на матрицы отражений, каждая из которых может быть определена одним вектором и имеет вид <math>U = E - 2ww^*</math>, где <math>w</math> - вектор, удовлетворяющий равенству <math>w^*w = 1</math>. | ||
+ | Цель домножения исходной матрицы на матрицу отражения на <math>i</math>-ом шаге заключается в обнулении поддиагональных элементов <math>i</math>-го столбца. Таким образом за <math>n-1</math> итерацию исходная матрица будет приведена к форме Хессенберга. | ||
+ | |||
+ | При этом матрица отражений на каждом из шагов имеет не стандартный вид, описанный ранее, а <math>U=E-\frac{1}{\gamma}vv^*</math>, где <math>v</math> находится для левых матриц отражения через координаты текущего <math>i</math>-го столбца при помощи операции скалярного произведения. | ||
+ | |||
+ | === QR-алгоритм со сдвигами === | ||
+ | Скорость сходимости QR-алгоритма зависит от того, насколько близки друг к другу собственные числа матрицы <math>A</math> - чем ближе они друг к другу, тем медленнее алгоритм сходится. Поэтому для ускорения сходимости используют так называемые сдвиги. | ||
+ | |||
+ | Суть метода заключается в следующем: | ||
+ | Пусть на <math>k</math>-й итерации есть матрица <math>A_k</math>, тогда переход к матрице <math>A_{k+1}</math> происходит следующим образом: | ||
+ | * Выбирается число <math>v_k</math> и выполняется QR-разложение: <math>A_k - v_kE = Q_kR_k</math>. | ||
+ | * <math>A_{k+1} = R_kQ_k + v_kE</math> | ||
+ | При этом матрицы <math>A_k</math> и <math>A_{k+1}</math> подобны, поэтому набор собственных чисел не изменяется. | ||
== Вычислительное ядро алгоритма == | == Вычислительное ядро алгоритма == | ||
+ | === Базовый QR-алгоритм === | ||
Основными вычислительными ядрами алгоритма являются: | Основными вычислительными ядрами алгоритма являются: | ||
− | * QR-разложение | + | * QR-разложение матрицы <math>A_k = Q_k R_k</math> |
− | * Умножение матриц | + | * Умножение плотных матриц <math>A_{k+1} = R_k Q_k</math> |
+ | |||
+ | Существуют различные алгоритмы вычисления QR-разложения. Можно выделить такие, как метод Гивенса и метод Хаусхолдера. | ||
+ | |||
+ | Подробное описание вычислительных ядер содержится в соответствующих статьях, указанных в пункте 1.4 (Макроструктура алгоритма). | ||
+ | |||
+ | === QR-алгоритм с приведением к форме Хессенберга и сдвигами === | ||
+ | Основными вычислительными ядрами данного варианта алгоритма являются: | ||
+ | * Приведение исходной матрицы <math>A</math> к матрице в форме Хессенберга. Данная операция производится методом Хаусхолдера. | ||
+ | * Выполнение QR-разложения матрицы <math>A_k - v_kE = Q_kR_k</math> модифицированным методом Гивенса. | ||
+ | * Перемножение двух плотных матриц <math>R_kQ_k</math>. | ||
== Макроструктура алгоритма == | == Макроструктура алгоритма == | ||
+ | === Базовый QR-алгоритм === | ||
QR-алгоритм на каждой итерации использует следующие макрооперации: | QR-алгоритм на каждой итерации использует следующие макрооперации: | ||
# QR-разложение | # QR-разложение | ||
Строка 37: | Строка 72: | ||
#* [[Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант|Метод Хаусхолдера (отражений)]] | #* [[Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант|Метод Хаусхолдера (отражений)]] | ||
# [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Умножение матриц]] | # [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Умножение матриц]] | ||
+ | |||
+ | === QR-алгоритм с приведением к форме Хессенберга и сдвигами === | ||
+ | В данном варианте макроструктура алгоритма выглядит следующим образом: | ||
+ | # Приведение исходной матрицы к матрице в форме Хессенберга. | ||
+ | # Итеративное выполнение QR-разложения методом Гивенса<ref name="QRGivens">[[Метод Гивенса (вращений) QR-разложения квадратной матрицы (вещественный вариант)|Метод Гивенса (вращений) QR-разложения квадратной матрицы]]</ref> и перемножения матриц. | ||
== Схема реализации последовательного алгоритма == | == Схема реализации последовательного алгоритма == | ||
− | + | === Базовый QR-алгоритм === | |
QR-алгоритм является итерационным. | QR-алгоритм является итерационным. | ||
На каждой итерации <math style="vertical-align:0%"> k </math>: | На каждой итерации <math style="vertical-align:0%"> k </math>: | ||
− | # Строится QR-разложение матрицы <math | + | # Строится QR-разложение матрицы <math> A_k = Q_k R_k </math> |
− | # Получается матрица <math | + | # Получается матрица <math> A_{k+1} = R_k Q_k </math>, используемая на следующей итерации |
− | Алгоритм выполняется до сходимости матрицы <math | + | Алгоритм выполняется до сходимости матрицы <math> A_k</math> к треугольному виду. |
Строка 58: | Строка 98: | ||
''Q'', ''R'' ← QR_factorization(''A'') | ''Q'', ''R'' ← QR_factorization(''A'') | ||
''A'' ← ''R''×''Q'' | ''A'' ← ''R''×''Q'' | ||
+ | '''until''' convergence | ||
+ | |||
+ | '''return''' diag(A) | ||
+ | |||
+ | === QR-алгоритм с приведением к форме Хессенберга и сдвигами === | ||
+ | Как уже было описано ранее, данный алгоритм состоит из следуюших этапов: | ||
+ | * Приведение исходной матрицы <math>A</math> к матрице в форме Хессенберга <math>H</math>. | ||
+ | * Выполнение <math>k</math> итераций (<math>A_0 = H</math>): | ||
+ | ** Выбор коэффициента <math>v_k</math>. | ||
+ | ** QR-разложение <math>A_k - v_kE = Q_kR_k</math>. | ||
+ | ** <math>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 | '''until''' convergence | ||
Строка 63: | Строка 127: | ||
== Последовательная сложность алгоритма == | == Последовательная сложность алгоритма == | ||
− | Пусть дана матрица <math | + | === Базовый QR-алгоритм === |
+ | Пусть дана матрица <math> A \in \mathbb{R}^{n \times n}</math>. До момента приведения матрицы к треугольной форме произведено <math>N</math> итераций алгоритма. | ||
На каждой итерации алгоритма производится три действия: QR-разложение, перемножение матриц и проверка, является ли матрица треугольной. | На каждой итерации алгоритма производится три действия: QR-разложение, перемножение матриц и проверка, является ли матрица треугольной. | ||
Распишем последовательную сложность каждого действия: | Распишем последовательную сложность каждого действия: | ||
* QR-разложение: | * QR-разложение: | ||
− | ** Методом Гивенса (вращений) - <math | + | ** Методом Гивенса (вращений) - <math>2*n^3</math><ref name="QRGivens">[[Метод Гивенса (вращений) QR-разложения квадратной матрицы (вещественный вариант)|Метод Гивенса (вращений) QR-разложения квадратной матрицы]]</ref>. |
− | ** Методом Хаусхолдера (отражений) - <math | + | ** Методом Хаусхолдера (отражений) - <math>\frac{4}{3}n^3</math><ref name="QRHouseholder">[[Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы, вещественный вариант|Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы]]</ref>. |
− | * Перемножение матриц - <math | + | * Перемножение матриц - <math>n^3</math> операций<ref name="MatrixMultiplication">[[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение плотных неособенных матриц]]</ref>. |
− | * Проверка, имеет ли матрица треугольную форму - <math | + | * Проверка, имеет ли матрица треугольную форму - <math>\frac{n^2 - n}{2}</math> операций. |
− | Таким образом, на каждой итерации последовательный алгоритм выполняет <math | + | Таким образом, на каждой итерации последовательный алгоритм выполняет <math>O(n^3)</math> операций. Для всего алгоритма потребуется выполнить <math>N * O(n^3)</math> операций. |
+ | |||
+ | === QR-алгоритм с приведением к форме Хессенберга и сдвигами === | ||
+ | Пусть как и для базового алгоритма дана матрица <math> A \in \mathbb{R}^{n \times n}</math>. Как уже было описано выше, в данном варианте алгоритма выполняются три основные операции: приведение исходной матрицы к матрице в форме Хессенберга, QR-разложение методом Гивенса и перемножение матриц. | ||
+ | |||
+ | Распишем последовательную сложность каждой операции: | ||
+ | * Приведение к матрице в форме Хессенберга методом Хаусхолдера потребует <math>O(n^3)</math> операций. | ||
+ | * QR-разложение методом Гивенса матрицы в форме Хессенберга - <math>6n^2 + O(n)</math> операций <ref>[https://en.wikipedia.org/wiki/QR_algorithm | QR algorithm]</ref>. | ||
+ | * Перемножение матриц потребует столько же операций, как и для базового алгоритма - <math>n^3</math> операций. | ||
+ | * Вычисление сдвигов не является затратной операцией по сравнению с другими. | ||
+ | |||
+ | Таким образом, полная последовательная сложность QR-алгоритма с матрицей в форме Хессенберга и сдвигами - <math>O(n^3) + N*(6n^2+O(n))</math> или <math>O(n^3) + N*O(n^2)</math> | ||
== Информационный граф == | == Информационный граф == | ||
Строка 79: | Строка 155: | ||
[[Файл:macro_shev_graph.png|thumb|center|700px|Рисунок 1. Макрограф базового QR-алгоритма.]] | [[Файл:macro_shev_graph.png|thumb|center|700px|Рисунок 1. Макрограф базового QR-алгоритма.]] | ||
[[Файл:Screenshot - 15.10.2016 - 19-02-55.png|thumb|center|800px|Рисунок 2. Информационный граф i-й итерации QR-алгоритма.]] | [[Файл:Screenshot - 15.10.2016 - 19-02-55.png|thumb|center|800px|Рисунок 2. Информационный граф i-й итерации QR-алгоритма.]] | ||
+ | |||
+ | В случае QR-алгоритма, в котором применяется приведение исходной матрицы к форме Хессенберга и сдвиги, представленные графы будут немного изменены: на Рисунке 1 перед первой итерацией появится элемент с приведением матрицы, а на Рисунке 2 добавятся элементы, отображающие сдвиги. | ||
== Ресурс параллелизма алгоритма == | == Ресурс параллелизма алгоритма == | ||
+ | === Базовый QR-алгоритм === | ||
Базовый алгоритм является итерационным и выполняется последовательно. | Базовый алгоритм является итерационным и выполняется последовательно. | ||
Возможность распараллеливания алгоритма предоставляется при реализации операций таких, как QR-разложение, перемножение матриц и проверка, имеет ли матрица треугольную форму. | Возможность распараллеливания алгоритма предоставляется при реализации операций таких, как QR-разложение, перемножение матриц и проверка, имеет ли матрица треугольную форму. | ||
Строка 86: | Строка 165: | ||
Рассчитаем параллельную сложность для каждой из указанных операций: | Рассчитаем параллельную сложность для каждой из указанных операций: | ||
* QR-разложение: | * QR-разложение: | ||
− | ** Методом Гивенса (вращений) - <math | + | ** Методом Гивенса (вращений) - <math>11n - 16</math><ref name="QRGivens">[[Метод Гивенса (вращений) QR-разложения квадратной матрицы (вещественный вариант)|Метод Гивенса (вращений) QR-разложения квадратной матрицы]]</ref>. |
− | ** Методом Хаусхолдера (отражений) - <math | + | ** Методом Хаусхолдера (отражений) - <math>O(n^2)</math><ref name="QRHouseholder">[[Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы, вещественный вариант|Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы]]</ref>. |
− | * Перемножение матриц - <math | + | * Перемножение матриц - <math>n</math><ref name="MatrixMultiplication">[[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение плотных неособенных матриц]]</ref>. |
− | * Проверка, имеет ли матрица треугольную форму - <math | + | * Проверка, имеет ли матрица треугольную форму - <math>1</math>. |
− | Результирующая параллельная сложность одной итерации <math | + | Результирующая параллельная сложность одной итерации <math>12n-15</math> (для метода вращений). А всего алгоритма из <math>N</math> итераций - <math>N*(12n-15)</math>. |
+ | |||
+ | === QR-алгоритм с приведением к форме Хессенберга и сдвигами === | ||
+ | Как и для базового алгоритма подсчитаем параллельную сложность для каждой операции, а в конце для всего алгоритма в целом: | ||
+ | * Приведением исходной матрицы к форме Хессенберга - <math>O(n^2)</math>. | ||
+ | * Отдельно каждая итерация - <math>O(n)</math>. | ||
+ | |||
+ | Итоговая параллельная сложность алгоритма - <math>O(n^2) + N*O(n)</math>. | ||
== Входные и выходные данные алгоритма == | == Входные и выходные данные алгоритма == | ||
− | '''Входные данные''': плотная квадратная матрица <math | + | '''Входные данные''': плотная квадратная матрица <math>A</math> размера <math>n</math>. |
− | '''Объём входных данных''': <math | + | '''Объём входных данных''': <math>n^2</math>. |
− | '''Выходные данные''': <math | + | '''Выходные данные''': <math>n</math> собственных чисел матрицы <math>A</math>. |
− | '''Объём выходных данных''': <math | + | '''Объём выходных данных''': <math>n</math>. |
== Свойства алгоритма == | == Свойства алгоритма == | ||
Строка 108: | Строка 194: | ||
* Скорость сходимости зависит от собственных чисел. Чем ближе собственные числа друг к другу, тем ниже скорость сходимости. | * Скорость сходимости зависит от собственных чисел. Чем ближе собственные числа друг к другу, тем ниже скорость сходимости. | ||
* Может быть хорошо распараллелен, так как соотношение последовательной и параллельной сложности алгоритма является квадратичным. | * Может быть хорошо распараллелен, так как соотношение последовательной и параллельной сложности алгоритма является квадратичным. | ||
− | * Перемещение данных необходимых для выполнения алгоритма не затратно, так как отношение числа операций к суммарному объему входных и выходных данных | + | * Перемещение данных необходимых для выполнения алгоритма не затратно, так как вычислительная мощность <math> = \frac{N * O(n^3)}{n^2 + n}</math> (отношение числа операций к суммарному объему входных и выходных данных) линейна на каждой итерации. |
= Программная реализация алгоритма = | = Программная реализация алгоритма = | ||
+ | == Особенности реализации последовательного алгоритма == | ||
+ | == Локальность данных и вычислений == | ||
+ | == Возможные способы и особенности параллельной реализации алгоритма == | ||
== Масштабируемость алгоритма и его реализации == | == Масштабируемость алгоритма и его реализации == | ||
+ | |||
+ | Используется процедура [http://www.netlib.org/scalapack/explore-html/d3/d3a/pdhseqr_8f.html PDHSEQR] из библиотеки [http://www.netlib.org/scalapack/ ScaLAPACK]. <br> | ||
+ | Данный метод работает с матрицами в форме Хессенберга, поэтому сначала происходит приведение матрицы к данному виду с помощью процедуры [http://www.netlib.org/scalapack/explore-html/d2/d03/pdgehrd_8f.html PDGEHRD] | ||
+ | |||
+ | Вычисления производились на суперкомпьютере [http://hpc.cmc.msu.ru/regatta Regatta]. | ||
+ | |||
+ | {| class="wikitable" | ||
+ | |+ Время работы алгоритма (сек) | ||
+ | |- | ||
+ | |Процессы\Матрица||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 | ||
+ | |} | ||
+ | |||
+ | [[Файл:baev_qr_execution.png|thumb|center|700px|Рисунок 3. Время выполнения QR-алгоритма.]] | ||
+ | |||
+ | |||
+ | Эффективность распараллеливания <math>E = S/p</math>, где <math>S = T_1/T_p</math> - полученное ускорение работы программы, а <math>p</math> - число используемых процессов. | ||
+ | |||
+ | {| class="wikitable" | ||
+ | |+ Эффективность распараллеливания (%) | ||
+ | |- | ||
+ | |Процессы\Матрица||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 | ||
+ | |} | ||
+ | |||
+ | [[Файл:baev_qr_efficiency.png|thumb|center|700px|Рисунок 4. Эффективность QR-алгоритма.]] | ||
+ | |||
+ | Оценка масштабируемости | ||
+ | * Минимальная эффективность: 8.62% | ||
+ | * Максимальная эффективность: 38.09% | ||
+ | * По числу процессов: при увеличении числа процессов эффективность в целом уменьшается | ||
+ | * По размеру задачи: при увеличении размера задачи эффективность в целом увеличивается | ||
+ | Уменьшение или увеличение эффективности происходит не интенсивно, существенные изменения заметны только при небольшом размере задачи. | ||
+ | |||
+ | Библиотека [http://www.netlib.org/scalapack/ ScaLAPACK 2.0.0] была установлена самостоятельно. | ||
+ | Использовался компилятор mpicxx с опциями -lscalapack -llapack и указанием пути до установленной библиотеки. | ||
+ | |||
+ | [https://algowiki-project.org/ru/Файл:Odbaev_eigenvalues.docx Реализация алгоритма] | ||
+ | |||
+ | == Динамические характеристики и эффективность реализации алгоритма == | ||
+ | == Выводы для классов архитектур == | ||
== Существующие реализации алгоритма == | == Существующие реализации алгоритма == | ||
Полный алгоритм: | Полный алгоритм: |
Текущая версия на 15:53, 14 декабря 2016
Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Frolov и VadimVV. |
Нахождение собственных чисел квадратной матрицы методом QR разложения | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(n^3) + N * O(n^2)[/math] |
Объём входных данных | [math]n^2[/math] |
Объём выходных данных | [math]n[/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math]O(n^2)+N*O(n)[/math] |
Ширина ярусно-параллельной формы | [math]O(n^2)[/math] |
Основные авторы описания: О.Д.Баев (пп. 1.1 - 1.5 базовый, 1.9, 2.4, 2.7), А.С.Шевелев (пп. 1.2-1.5 хессенберг, 1.6 - 1.8, 1.10)
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
QR-алгоритм — это численный метод в линейной алгебре, предназначенный для решения полной проблемы собственных значений, то есть отыскания всех собственных чисел и собственных векторов матрицы. Был разработан в конце 1950-х годов независимо В.Н. Кублановской и Дж. Фрэнсисом.
1.2 Математическое описание алгоритма
1.2.1 Базовый QR-алгоритм
Пусть [math]A[/math] — вещественная квадратная матрица, для которой мы хотим найти собственные числа и векторы. Положим [math]A_0 = A[/math]. На k-ом шаге (начиная с k = 0) вычислим QR-разложение [math]A_k = Q_k R_k[/math], где [math]Q_k[/math] — ортогональная матрица (то есть [math]Q_k^{T} = Q_k^{-1}[/math]), а [math]R_k[/math] — верхняя треугольная матрица. Затем мы определяем [math]A_{k+1} = R_k Q_k[/math].
Заметим, что
- [math] 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]A_k[/math] являются подобными и их собственные значения равны.
Пусть все диагональные миноры матрицы [math]A[/math] не вырождены. Тогда последовательность матриц [math]A_k[/math] при [math]k \rightarrow \infty[/math] сходится по форме к клеточному правому треугольному виду, соответствующему клеткам с одинаковыми по модулю собственными значениями. [1]
Для того, чтобы получить собственные векторы матрицы, нужно перемножить все матрицы [math]Q_k[/math].
1.2.2 Приведение матрицы к форме Хессенберга
Матрица [math]A[/math] имеет форму Хессенберга[2], если все элементы [math]a_{ij}[/math] под первой поддиагональю равны нулю ([math]a_{ij} = 0, i \gt j + 1[/math]).
Любая матрица [math]A[/math] может быть приведена к подобной матрице, имеющей форму Хессенберга. Подобная манипуляция позволяет значительно ускорить дальнейший процесс QR-разложения. Одним из способов приведения матрицы к форме Хессенберга является метод Хаусхолдера.
Идея метода состоит в последовательном приведении исходной матрицы [math]A[/math] к форме Хессенберга за счет умножения на матрицы отражений, каждая из которых может быть определена одним вектором и имеет вид [math]U = E - 2ww^*[/math], где [math]w[/math] - вектор, удовлетворяющий равенству [math]w^*w = 1[/math]. Цель домножения исходной матрицы на матрицу отражения на [math]i[/math]-ом шаге заключается в обнулении поддиагональных элементов [math]i[/math]-го столбца. Таким образом за [math]n-1[/math] итерацию исходная матрица будет приведена к форме Хессенберга.
При этом матрица отражений на каждом из шагов имеет не стандартный вид, описанный ранее, а [math]U=E-\frac{1}{\gamma}vv^*[/math], где [math]v[/math] находится для левых матриц отражения через координаты текущего [math]i[/math]-го столбца при помощи операции скалярного произведения.
1.2.3 QR-алгоритм со сдвигами
Скорость сходимости QR-алгоритма зависит от того, насколько близки друг к другу собственные числа матрицы [math]A[/math] - чем ближе они друг к другу, тем медленнее алгоритм сходится. Поэтому для ускорения сходимости используют так называемые сдвиги.
Суть метода заключается в следующем: Пусть на [math]k[/math]-й итерации есть матрица [math]A_k[/math], тогда переход к матрице [math]A_{k+1}[/math] происходит следующим образом:
- Выбирается число [math]v_k[/math] и выполняется QR-разложение: [math]A_k - v_kE = Q_kR_k[/math].
- [math]A_{k+1} = R_kQ_k + v_kE[/math]
При этом матрицы [math]A_k[/math] и [math]A_{k+1}[/math] подобны, поэтому набор собственных чисел не изменяется.
1.3 Вычислительное ядро алгоритма
1.3.1 Базовый QR-алгоритм
Основными вычислительными ядрами алгоритма являются:
- QR-разложение матрицы [math]A_k = Q_k R_k[/math]
- Умножение плотных матриц [math]A_{k+1} = R_k Q_k[/math]
Существуют различные алгоритмы вычисления QR-разложения. Можно выделить такие, как метод Гивенса и метод Хаусхолдера.
Подробное описание вычислительных ядер содержится в соответствующих статьях, указанных в пункте 1.4 (Макроструктура алгоритма).
1.3.2 QR-алгоритм с приведением к форме Хессенберга и сдвигами
Основными вычислительными ядрами данного варианта алгоритма являются:
- Приведение исходной матрицы [math]A[/math] к матрице в форме Хессенберга. Данная операция производится методом Хаусхолдера.
- Выполнение QR-разложения матрицы [math]A_k - v_kE = Q_kR_k[/math] модифицированным методом Гивенса.
- Перемножение двух плотных матриц [math]R_kQ_k[/math].
1.4 Макроструктура алгоритма
1.4.1 Базовый QR-алгоритм
QR-алгоритм на каждой итерации использует следующие макрооперации:
1.4.2 QR-алгоритм с приведением к форме Хессенберга и сдвигами
В данном варианте макроструктура алгоритма выглядит следующим образом:
- Приведение исходной матрицы к матрице в форме Хессенберга.
- Итеративное выполнение QR-разложения методом Гивенса[3] и перемножения матриц.
1.5 Схема реализации последовательного алгоритма
1.5.1 Базовый QR-алгоритм
QR-алгоритм является итерационным.
На каждой итерации [math] k [/math]:
- Строится QR-разложение матрицы [math] A_k = Q_k R_k [/math]
- Получается матрица [math] A_{k+1} = R_k Q_k [/math], используемая на следующей итерации
Алгоритм выполняется до сходимости матрицы [math] 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)
1.5.2 QR-алгоритм с приведением к форме Хессенберга и сдвигами
Как уже было описано ранее, данный алгоритм состоит из следуюших этапов:
- Приведение исходной матрицы [math]A[/math] к матрице в форме Хессенберга [math]H[/math].
- Выполнение [math]k[/math] итераций ([math]A_0 = H[/math]):
- Выбор коэффициента [math]v_k[/math].
- QR-разложение [math]A_k - v_kE = Q_kR_k[/math].
- [math]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)
1.6 Последовательная сложность алгоритма
1.6.1 Базовый QR-алгоритм
Пусть дана матрица [math] A \in \mathbb{R}^{n \times n}[/math]. До момента приведения матрицы к треугольной форме произведено [math]N[/math] итераций алгоритма. На каждой итерации алгоритма производится три действия: QR-разложение, перемножение матриц и проверка, является ли матрица треугольной.
Распишем последовательную сложность каждого действия:
- QR-разложение:
- Перемножение матриц - [math]n^3[/math] операций[5].
- Проверка, имеет ли матрица треугольную форму - [math]\frac{n^2 - n}{2}[/math] операций.
Таким образом, на каждой итерации последовательный алгоритм выполняет [math]O(n^3)[/math] операций. Для всего алгоритма потребуется выполнить [math]N * O(n^3)[/math] операций.
1.6.2 QR-алгоритм с приведением к форме Хессенберга и сдвигами
Пусть как и для базового алгоритма дана матрица [math] A \in \mathbb{R}^{n \times n}[/math]. Как уже было описано выше, в данном варианте алгоритма выполняются три основные операции: приведение исходной матрицы к матрице в форме Хессенберга, QR-разложение методом Гивенса и перемножение матриц.
Распишем последовательную сложность каждой операции:
- Приведение к матрице в форме Хессенберга методом Хаусхолдера потребует [math]O(n^3)[/math] операций.
- QR-разложение методом Гивенса матрицы в форме Хессенберга - [math]6n^2 + O(n)[/math] операций [6].
- Перемножение матриц потребует столько же операций, как и для базового алгоритма - [math]n^3[/math] операций.
- Вычисление сдвигов не является затратной операцией по сравнению с другими.
Таким образом, полная последовательная сложность QR-алгоритма с матрицей в форме Хессенберга и сдвигами - [math]O(n^3) + N*(6n^2+O(n))[/math] или [math]O(n^3) + N*O(n^2)[/math]
1.7 Информационный граф
На Рисунке 1 изображен макрограф описываемого QR-алгоритма, представляющего собой последовательность итераций алгоритма. На рисунке 2 представлен граф, на котором более подробно отображено содержимое каждой итерации. Информационные графы QR-разложения и матричного умножения представлены в соответствующих статьях.
В случае QR-алгоритма, в котором применяется приведение исходной матрицы к форме Хессенберга и сдвиги, представленные графы будут немного изменены: на Рисунке 1 перед первой итерацией появится элемент с приведением матрицы, а на Рисунке 2 добавятся элементы, отображающие сдвиги.
1.8 Ресурс параллелизма алгоритма
1.8.1 Базовый QR-алгоритм
Базовый алгоритм является итерационным и выполняется последовательно. Возможность распараллеливания алгоритма предоставляется при реализации операций таких, как QR-разложение, перемножение матриц и проверка, имеет ли матрица треугольную форму.
Рассчитаем параллельную сложность для каждой из указанных операций:
- QR-разложение:
- Перемножение матриц - [math]n[/math][5].
- Проверка, имеет ли матрица треугольную форму - [math]1[/math].
Результирующая параллельная сложность одной итерации [math]12n-15[/math] (для метода вращений). А всего алгоритма из [math]N[/math] итераций - [math]N*(12n-15)[/math].
1.8.2 QR-алгоритм с приведением к форме Хессенберга и сдвигами
Как и для базового алгоритма подсчитаем параллельную сложность для каждой операции, а в конце для всего алгоритма в целом:
- Приведением исходной матрицы к форме Хессенберга - [math]O(n^2)[/math].
- Отдельно каждая итерация - [math]O(n)[/math].
Итоговая параллельная сложность алгоритма - [math]O(n^2) + N*O(n)[/math].
1.9 Входные и выходные данные алгоритма
Входные данные: плотная квадратная матрица [math]A[/math] размера [math]n[/math].
Объём входных данных: [math]n^2[/math].
Выходные данные: [math]n[/math] собственных чисел матрицы [math]A[/math].
Объём выходных данных: [math]n[/math].
1.10 Свойства алгоритма
Описанный алгоритм обладает следующими свойствами:
- Недетерминирован - не известно заранее число итераций алгоритма, которое необходимо произвести до момента схождения матрицы к полностью треугольной форме или с некоторым значением точности.
- Скорость сходимости зависит от собственных чисел. Чем ближе собственные числа друг к другу, тем ниже скорость сходимости.
- Может быть хорошо распараллелен, так как соотношение последовательной и параллельной сложности алгоритма является квадратичным.
- Перемещение данных необходимых для выполнения алгоритма не затратно, так как вычислительная мощность [math] = \frac{N * O(n^3)}{n^2 + n}[/math] (отношение числа операций к суммарному объему входных и выходных данных) линейна на каждой итерации.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Используется процедура 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]E = S/p[/math], где [math]S = T_1/T_p[/math] - полученное ускорение работы программы, а [math]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 и указанием пути до установленной библиотеки.
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
Полный алгоритм:
- 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
3 Литература
- ↑ Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы" — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.
- ↑ | Матрица Хессенберга
- ↑ 3,0 3,1 3,2 Метод Гивенса (вращений) QR-разложения квадратной матрицы
- ↑ 4,0 4,1 Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы
- ↑ 5,0 5,1 Перемножение плотных неособенных матриц
- ↑ | QR algorithm