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

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

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

Основные авторы описания: О.Д.Баев (пп. 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-алгоритм

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

Заметим, что

$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,$

то есть все матрицы $A_k$ являются подобными и их собственные значения равны.

Пусть все диагональные миноры матрицы $A$ не вырождены. Тогда последовательность матриц $A_k$ при $k \rightarrow \infty$ сходится по форме к клеточному правому треугольному виду, соответствующему клеткам с одинаковыми по модулю собственными значениями.

Для того, чтобы получить собственные векторы матрицы, нужно перемножить все матрицы $Q_k$.

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

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

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

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

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

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

Скорость сходимости QR-алгоритма зависит от того, насколько близки друг к другу собственные числа матрицы $A$ - чем ближе они друг к другу, тем медленнее алгоритм сходится. Поэтому для ускорения сходимости используют так называемые сдвиги.

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

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

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

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

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

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

• QR-разложение матрицы $A_k = Q_k R_k$
• Умножение плотных матриц $A_{k+1} = R_k Q_k$

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

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

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

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

• Приведение исходной матрицы $A$ к матрице в форме Хессенберга. Данная операция производится методом Хаусхолдера.
• Выполнение QR-разложения матрицы $A_k - v_kE = Q_kR_k$ модифицированным методом Гивенса.
• Перемножение двух плотных матриц $R_kQ_k$.

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

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

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

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

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

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

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

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

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

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

На каждой итерации $k$:

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

Алгоритм выполняется до сходимости матрицы $A_k$ к треугольному виду.

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

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-алгоритм с приведением к форме Хессенберга и сдвигами

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

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

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

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-алгоритм

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

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

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

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

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

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

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

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

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

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

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

Рисунок 1. Макрограф базового QR-алгоритма.

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

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

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

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

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

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

• QR-разложение:
  • Методом Гивенса (вращений) - $11n - 16$.
  • Методом Хаусхолдера (отражений) - $O(n^2)$.
• Перемножение матриц - $n$.
• Проверка, имеет ли матрица треугольную форму - $1$.

Результирующая параллельная сложность одной итерации $12n-15$ (для метода вращений). А всего алгоритма из $N$ итераций - $N*(12n-15)$.

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

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

• Приведением исходной матрицы к форме Хессенберга - $O(n^2)$.
• Отдельно каждая итерация - $O(n)$.

Итоговая параллельная сложность алгоритма - $O(n^2) + N*O(n)$.

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

Входные данные: плотная квадратная матрица $A$ размера $n$.

Объём входных данных: $n^2$.

Выходные данные: $n$ собственных чисел матрицы $A$.

Объём выходных данных: $n$.

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

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

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

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

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-алгоритма.

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

Эффективность распараллеливания (%)

Процессы\Матрица	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.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