Уровень алгоритма

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

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показаны 144 промежуточные версии 5 участников)
Строка 1: Строка 1:
 +
{{Assignment|VadimVV|Frolov}}
 +
 
{{algorithm
 
{{algorithm
 
| name              = Нахождение собственных чисел квадратной матрицы методом QR разложения
 
| name              = Нахождение собственных чисел квадратной матрицы методом QR разложения
−
| serial_complexity =  
+
| serial_complexity = <math>O(n^3) + N * O(n^2)</math>
−
| pf_height        =  
+
| pf_height        = <math>O(n^2)+N*O(n)</math>
−
| pf_width          =  
+
| pf_width          = <math>O(n^2)</math>
 
| input_data        = <math>n^2</math>
 
| input_data        = <math>n^2</math>
 
| output_data      = <math>n</math>
 
| output_data      = <math>n</math>
 
}}
 
}}
  
−
Основные авторы описания: [[Участник: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:
  
 
== Математическое описание алгоритма ==
 
== Математическое описание алгоритма ==
−
Пусть ''A'' — вещественная матрица, для которой мы хотим найти собственные числа и векторы. Положим ''A''<sub>0</sub>=''A''. На ''k''-ом шаге (начиная с ''k'' = 0) вычислим QR-разложение ''A''<sub>''k''</sub>=''Q''<sub>''k''</sub>''R''<sub>''k''</sub>, где ''Q''<sub>''k''</sub> — ортогональная матрица (то есть ''Q''<sub>''k''</sub><sup>''T''</sup> = ''Q''<sub>''k''</sub><sup>−1</sup>), а ''R''<sub>''k''</sub> — верхняя треугольная матрица. Затем мы определяем ''A''<sub>''k''+1</sub> = ''R''<sub>''k''</sub>''Q''<sub>''k''</sub>.
+
 
 +
=== Базовый 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+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>
−
то есть все матрицы ''A''<sub>''k''</sub> являются подобными и их собственные значения равны.
+
то есть все матрицы <math>A_k</math> являются подобными и их собственные значения равны.
 +
 
 +
Пусть все диагональные миноры матрицы <math>A</math> не вырождены. Тогда последовательность матриц <math>A_k</math> при <math>k \rightarrow \infty</math> сходится по форме к клеточному правому треугольному виду, соответствующему клеткам с одинаковыми по модулю собственными значениями. <ref> Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы" — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с. </ref>
 +
 
 +
Для того, чтобы получить собственные векторы матрицы, нужно перемножить все матрицы <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> итерацию исходная матрица будет приведена к форме Хессенберга.
  
−
Пусть все диагональные миноры матрицы ''A'' не вырождены. Тогда последовательность матриц ''A''<sub>''k''</sub> при ''k''→∞ сходится по форме к клеточному правому треугольному виду, соответствующему клеткам с одинаковыми по модулю собственными значениями. <ref> Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с. </ref>
+
При этом матрица отражений на каждом из шагов имеет не стандартный вид, описанный ранее, а <math>U=E-\frac{1}{\gamma}vv^*</math>, где <math>v</math> находится для левых матриц отражения через координаты текущего <math>i</math>-го столбца при помощи операции скалярного произведения.
  
−
Для того, чтобы получить собственные векторы матрицы, нужно перемножить все матрицы ''Q''<sub>''k''</sub>.
+
=== 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-алгоритм является итерационным.
 +
 +
На каждой итерации <math style="vertical-align:0%"> 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'' &larr; QR_factorization(''A'')
 +
    ''A'' &larr; ''R''&times;''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'' &larr; To_Hessenberg(''A'')
 +
  ''A &larr; H''
 +
  '''repeat'''
 +
    ''v_k'' &larr; Get_coeff(''A'')
 +
    ''Q'', ''R'' &larr; QR_factorization(''A - v_k*E'')
 +
    ''A'' &larr; ''R''&times;''Q + v_k*E''
 +
  '''until''' convergence
 +
   
 +
  '''return''' diag(A)
 +
 
== Последовательная сложность алгоритма ==
 
== Последовательная сложность алгоритма ==
 +
=== Базовый QR-алгоритм ===
 +
Пусть дана матрица <math> A \in \mathbb{R}^{n \times n}</math>. До момента приведения матрицы к треугольной форме произведено <math>N</math> итераций алгоритма.
 +
На каждой итерации алгоритма производится три действия: QR-разложение, перемножение матриц и проверка, является ли матрица треугольной.
 +
 +
Распишем последовательную сложность каждого действия:
 +
* QR-разложение:
 +
** Методом Гивенса (вращений) - <math>2*n^3</math><ref name="QRGivens">[[Метод Гивенса (вращений) QR-разложения квадратной матрицы (вещественный вариант)|Метод Гивенса (вращений) QR-разложения квадратной матрицы]]</ref>.
 +
** Методом Хаусхолдера (отражений) - <math>\frac{4}{3}n^3</math><ref name="QRHouseholder">[[Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы, вещественный вариант|Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы]]</ref>.
 +
* Перемножение матриц - <math>n^3</math> операций<ref name="MatrixMultiplication">[[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение плотных неособенных матриц]]</ref>.
 +
* Проверка, имеет ли матрица треугольную форму - <math>\frac{n^2 - n}{2}</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>
 +
 
== Информационный граф ==
 
== Информационный граф ==
 +
На Рисунке 1 изображен макрограф описываемого QR-алгоритма, представляющего собой последовательность итераций алгоритма. На рисунке 2 представлен граф, на котором более подробно отображено содержимое каждой итерации.
 +
Информационные графы [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)#.D0.98.D0.BD.D1.84.D0.BE.D1.80.D0.BC.D0.B0.D1.86.D0.B8.D0.BE.D0.BD.D0.BD.D1.8B.D0.B9_.D0.B3.D1.80.D0.B0.D1.84|QR-разложения]] и [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)#.D0.98.D0.BD.D1.84.D0.BE.D1.80.D0.BC.D0.B0.D1.86.D0.B8.D0.BE.D0.BD.D0.BD.D1.8B.D0.B9_.D0.B3.D1.80.D0.B0.D1.84|матричного умножения]] представлены в соответствующих статьях.
 +
[[Файл:macro_shev_graph.png|thumb|center|700px|Рисунок 1. Макрограф базового QR-алгоритма.]]
 +
[[Файл:Screenshot - 15.10.2016 - 19-02-55.png|thumb|center|800px|Рисунок 2. Информационный граф i-й итерации QR-алгоритма.]]
 +
 +
В случае QR-алгоритма, в котором применяется приведение исходной матрицы к форме Хессенберга и сдвиги, представленные графы будут немного изменены: на Рисунке 1 перед первой итерацией появится элемент с приведением матрицы, а на Рисунке 2 добавятся элементы, отображающие сдвиги.
 +
 
== Ресурс параллелизма алгоритма ==
 
== Ресурс параллелизма алгоритма ==
 +
=== Базовый QR-алгоритм ===
 +
Базовый алгоритм является итерационным и выполняется последовательно.
 +
Возможность распараллеливания алгоритма предоставляется при реализации операций таких, как QR-разложение, перемножение матриц и проверка, имеет ли матрица треугольную форму.
 +
 +
Рассчитаем параллельную сложность для каждой из указанных операций:
 +
* QR-разложение:
 +
** Методом Гивенса (вращений) - <math>11n - 16</math><ref name="QRGivens">[[Метод Гивенса (вращений) QR-разложения квадратной матрицы (вещественный вариант)|Метод Гивенса (вращений) QR-разложения квадратной матрицы]]</ref>.
 +
** Методом Хаусхолдера (отражений) - <math>O(n^2)</math><ref name="QRHouseholder">[[Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы, вещественный вариант|Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы]]</ref>.
 +
* Перемножение матриц - <math>n</math><ref name="MatrixMultiplication">[[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение плотных неособенных матриц]]</ref>.
 +
* Проверка, имеет ли матрица треугольную форму - <math>1</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>A</math> (элементы <math>a_{ij}</math>).
+
'''Входные данные''': плотная квадратная матрица <math>A</math> размера <math>n</math>.
  
 
'''Объём входных данных''': <math>n^2</math>.
 
'''Объём входных данных''': <math>n^2</math>.
  
−
'''Выходные данные''': собственные числа матрицы <math>A</math>.
+
'''Выходные данные''': <math>n</math> собственных чисел матрицы <math>A</math>.
  
 
'''Объём выходных данных''': <math>n</math>.
 
'''Объём выходных данных''': <math>n</math>.
  
 
== Свойства алгоритма ==
 
== Свойства алгоритма ==
 +
Описанный алгоритм обладает следующими свойствами:
 +
* Недетерминирован - не известно заранее число итераций алгоритма, которое необходимо произвести до момента схождения матрицы к полностью треугольной форме или с некоторым значением точности.
 +
* Скорость сходимости зависит от собственных чисел. Чем ближе собственные числа друг к другу, тем ниже скорость сходимости.
 +
* Может быть хорошо распараллелен, так как соотношение последовательной и параллельной сложности алгоритма является квадратичным.
 +
* Перемещение данных необходимых для выполнения алгоритма не затратно, так как вычислительная мощность <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

Symbol confirmed.svgЭта работа успешно выполнена
Преподавателю: в основное пространство, в подстраницу
Данное задание было проверено и зачтено.
Проверено 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 Общее описание алгоритма

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. QR-разложение
  2. Умножение матриц

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

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

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

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

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

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

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

  1. Строится QR-разложение матрицы [math] A_k = Q_k R_k [/math]
  2. Получается матрица [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]2*n^3[/math][3].
    • Методом Хаусхолдера (отражений) - [math]\frac{4}{3}n^3[/math][4].
  • Перемножение матриц - [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-разложения и матричного умножения представлены в соответствующих статьях.

Рисунок 1. Макрограф базового QR-алгоритма.
Рисунок 2. Информационный граф i-й итерации QR-алгоритма.

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

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

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

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

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

  • QR-разложение:
    • Методом Гивенса (вращений) - [math]11n - 16[/math][3].
    • Методом Хаусхолдера (отражений) - [math]O(n^2)[/math][4].
  • Перемножение матриц - [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
Рисунок 3. Время выполнения QR-алгоритма.


Эффективность распараллеливания [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
Рисунок 4. Эффективность QR-алгоритма.

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

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

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

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

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

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

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

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

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

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

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

3 Литература

  1. ↑ Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы" — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.
  2. ↑ | Матрица Хессенберга
  3. ↑ 3,0 3,1 3,2 Метод Гивенса (вращений) QR-разложения квадратной матрицы
  4. ↑ 4,0 4,1 Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы
  5. ↑ 5,0 5,1 Перемножение плотных неособенных матриц
  6. ↑ | QR algorithm
Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Odbaev/Нахождение_собственных_чисел_квадратной_матрицы_методом_QR_разложения_(2)&oldid=20916