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

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

Основные авторы описания: Смирнова А.С., Киямова А.

---

Задача нахождения собственных значений и собственных векторов для матрицы $A$ заключается в поиске таких чисел $\lambda$, которые удовлетворяют уравнению:

$Ax=\lambda x$, при этом, числа $\lambda$ называются собственными значениями, а вектора $x$ - собственными векторами

Данная задача является одной из важнейших задач линейной алгебры. Собственные вектора и собственные значения применяются в различных областях науки: в аналитической геометрии, при решении систем интегральных уравнений, в математической физике. Однако не существует простых алгоритмов прямого вычисления собственных значений для матриц в общем случае, поэтому данная задача на практике решается численными методами. Одним из таких методов является QR-алгоритм.

1 ЧАСТЬ. Свойства и структура алгоритма

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

---

QR-алгоритм — это численный метод в линейной алгебре, предназначенный для решения полной проблемы собственных значений, то есть отыскания всех собственных чисел матрицы. При этом алгоритм позволяет найти и собственные вектора матрицы. Он был разработан в конце 1950-х годов независимо В. Н. Кублановской(Россия) и Дж. Фрэнсисом(Англия). Открытию QR-алгоритма предшествовал LR-алгоритм, который использовал LU-разложение вместо QR-разложения. В настоящее время LR-алгоритм используется очень редко ввиду своей меньшей эффективности, однако он был важным шагом на пути к открытию QR-алгоритма.

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

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

• Перед итерациями привести матрицу $A$ к подобной ей матрице $A_H$, которая будет иметь форму Хессенберга. Данный шаг позволит ускорить процесс QR-разложения
• Использовать QR-алгоритм со сдвигами

В данной статье будет рассмотрен только базовый QR-алгоритм.

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

---

1.2.1 QR-разложение матрицы

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

$A=QR$, где $Q$ - ортогональная матрица ($Q^T=Q^{-1}$), $R$ - верхняя треугольная матрица

Данное разложение называется QR-разложением.

Есть несколько алгоритмов вычисления QR-разложения матрицы . Кратко опишем их в данной статье.

1.2.1.1 Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы

Основная статья: Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы

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

На $i$-ом шаге задача $i$-ой матрицы отражения заключается в обнулении всех поддиагональных элементов $i$-го столбца матрицы $A$ (при этом столбцы левее $i$-го не изменяются). Таким образом, алгоритм состоит из $n$ шагов, на каждом из которых вычисляется очередная матрица отражения, после чего найденное отражение применяется к матрице, которая является результатом предыдущего шага.

Матрица отражения имеет вид $E-\frac{1}{\gamma }vv^*$, где вектор $v$ вычисляется из текущего $i$-го столбца матрицы при помощи использования операции скалярного произведения векторов. Данное представление матриц отражения позволяет хранить их в виде одного вектора и сводит домножение матрицы отражения на текущую матрицу к арифметическим операциям над векторами (скалярное произведение и сложение векторов).

1.2.1.2 Метод Гивенса (вращений) QR-разложения квадратной матрицы

Основная статья: Метод Гивенса (вращений) QR-разложения квадратной матрицы

Суть метода Гивенса заключается в последовательном приведении матрицы $A$ к верхней треугольной форме при помощи домножения ее на так называемые матрицы вращения. Получившаяся треугольная матрица будет искомой матрицей $R$, а матрица $Q$ будет равна произведению сопряженных матриц вращения.

На каждом шаге задачей матрицы вращения является обнуление одного поддиагонального элемента. Вначале обнуляются все поддиагональные элементы 1-го столбца, затем 2-го и так далее до $(n-1)$-го. Таким образом, алгоритм состоит из $\frac{n(n-1)}{2}$ шагов на каждом из которых вычисляется очередная матрица вращения, после чего она применяется к матрице, которая является результатом предыдущего шага.

Матрицы вращения $T_{ij}$ по ствоей структуре очень близки к единичным матрицам, за исключением четырех элементов: элементы с номерами $ii$ и $jj$ содержат некоторое число-параметр $c$, элементы с номерами $ij$ и $ji$ содержат числа-параметры $-s$ и $s$ соответственно. Вычисление параметров $c$ и $s$ происходит на каждом шаге в зависимости от текущей матрицы и состоит из простых численных арифметических операций. Умножение матрицы вращения на текущую матрицу может быть представлено как последовательное изменение элементов с номерами $ik$ и $kk$ для всех столбцов $k$ находящихся правее столбца $i$. Каждое такое изменение по своей структуре эквивалентно операции перемножения двух комплексных чисел.

1.2.2 QR-алгоритм нахождения собственных чисел

Пусть матрица $A$ - матрица, для которой мы хотим найти собственные значения. Тогда итерационный процесс строится следующим образом:

1. $A_0=A$
2. Пусть имеется матрица $A_k$, тогда матрица $A_{k+1}$ строится следующим образом:

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

Заметим, что $A_{k+1}=R_kQ_k={Q_{k}^{-1}}Q_kR_kQ_k={Q_{k}^{-1}}A_kQ_k={Q_{k}^{T}}A_kQ_k$

Таким образом матрицы $A_{k+1}$ и $A_{k}$ подобны для $\forall k$, а значит, в силу транзитивности подобия, все матрицы $A_{k}$ подобны матрице $A$ и имеют одинаковый набор собственных значений.

1.2.3 Сходимость алгоритма

Предположим, что для $\forall m$ выполнены следующие условия:

1. $A=X\Lambda X^{-1}$, где $\Lambda =\begin{bmatrix} \Lambda_1 & 0\\ 0 & \Lambda_2 \end{bmatrix}, \Lambda_1\in\mathbb{C}^{m\times m},\Lambda_2\in\mathbb{C}^{r\times r}$

2. $\left | \lambda_1 \right | \geq ...\geq \left | \lambda_m \right | \gt \left | \lambda_{m+1} \right | \geq ...\geq \left | \lambda_{m+r} \right | \gt 0$, где $\{\lambda_1,...,\lambda_m\} = \lambda(\Lambda_1), \{\lambda_{m+1},...,\lambda_{m+r}\} = \lambda(\Lambda_2)$

3. Ведущая подматрица порядка $m$ в $X^{-1}$ невырождена

Тогда при $k \rightarrow \infty$ последовательность матриц $A_k$ сходится к матрице с верхней треугольной формой.

Таким образом, на практике необходимо выполнять итерации до тех пор пока матрица $A_k$ не станет треугольной (также можно продолжать выполнять их пока искомая матрица не будет найдена с некоторой заранее заданной точностью $\varepsilon$). Если итерации закончились на шаге $N$, то числа на диагонали матрицы $A_N$ будут считаться собственными значениями матрицы $A$

1.2.4 Вещественный вариант QR-алгоритма

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

$A_N= \begin{bmatrix} \blacksquare& \bullet& \bullet& \cdots& \cdots& \cdots& \cdots& \cdots& \bullet\\ 0& \blacksquare& \blacksquare& \bullet& \ddots& \ddots& \ddots& \ddots& \vdots\\ 0& \blacksquare& \blacksquare& \bullet& \bullet& \ddots& \ddots& \ddots& \vdots\\ \vdots& 0& 0& \blacksquare& \bullet& \bullet& \ddots& \ddots& \vdots\\ \vdots& \ddots& 0& 0& \blacksquare& \bullet& \bullet& \ddots& \vdots\\ \vdots& \ddots& \ddots& 0& 0& \blacksquare& \blacksquare& \ddots& \vdots\\ \vdots& \ddots& \ddots& \ddots& 0& \blacksquare& \blacksquare& \ddots& \bullet\\ \vdots& \ddots& \ddots& \ddots& \ddots& \ddots& \ddots& \ddots& \bullet\\ 0& \cdots& \cdots& \cdots& \cdots& \cdots& 0& 0& \blacksquare \end{bmatrix}$

В дальнейшем, в данной статье, матрицы, имеющие вышеописанную форму, будут называться квазитреугольными.

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

---

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

1. Вычисление QR-разложения матрицы: $A_k=Q_kR_k$. Существует несколько методов решения данной задачи:

• Метод Хаусхолдера (отражений)

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

• Метод Гивенса (вращений)

Вычислительное ядро данного алгоритма состоит из численных арифметических операций, необходимых для вычисления параметров матрицы вращения, и из операций (эквивалентных перемножению комплексных чисел), необходимых для пересчета матрицы на каждом шаге.

2. Перемножение двух плотных матриц: $A_{k+1}=R_kQ_k$

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

---

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

1.Вычисление QR-разложения матрицы: $A_k=Q_kR_k$.

2.Перемножение двух плотных матриц: $A_{k+1}=R_kQ_k$

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

---

Опишем необходимый для реализации цикл при помощи псевдокода:

A - исходная матрица
curA - текущая матрица, на основе которой будет вычисляться QR-разложение на каждом шаге
nextA - новая матрица, полученная после перемножения матриц R и Q 
triangular - функция, проверяющая, имеет ли матрица квазитреугольную форму.
difference - функция, проверяющая, что элементы двух матриц, стоящие на одинаковых местах, различаются не более чем на eps (данная функция нужна чтобы проверять не только сходимость матрицы к квазитреугольной форме, но и сходимость ее ненулевых элементов).

curA = A
nextA = A
while ( not (triangular(nextA) & difference(curA,nextA,eps)) )
{
    curA = nextA
    findQRdecomposition(curA,Q,R)
    nextA = R*Q
}

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

---

Подсчитаем сложность одной итерации QR-алгоритма (расчет сложностей для QR-разложения и перемножения матриц представлен в статьях, посвященных данным алгоритмам)

1. QR-разложение матрицы
   • Метод Хаусхолдера (отражений) имеет сложность $\frac{4}{3}n^3$
   • Метод Гивенса (вращений) имеет сложность $2n^3$
2. Перемножение двух плотных матриц имеет сложность $n^3$
3. Проверка матрицы на квазитреугольную форму состоит из набора сравнений элементов с номерами $ij$ ($i\gt j+1$) с нулем (таких элементов $\frac{(n-2)(n-1)}{2}$). Поддиагональные элементы с номерами $ij$ ($i=j+1$) должны быть проверены на соответствие необходимой квазитреугольной форме. Для этого достаточно для каждого такого элемента проверить следующее условие: $a_{i+1,i}==0 \vee a_{i,i-1}==0 \wedge a_{i+2,i+1}==0$ (либо поддиагональный элемент равен 0, либо, в противном случае, соседние поддиагональные элементы равны 0, чтобы текущий элемент соответствовал блоку 2-го порядка). Таких проверок поддиагональных элементов будет $n-1$. После данных проверок следует набор логических операций "И" между результатами всех сравнений (таких операций будет $\frac{n(n-1)}{2}-1$)
4. Сравнение новой матрицы с предыдущей состоит из операций вычитания и сравнения для каждой пары ненулевых соответствующих элементов (такие пары имеют номера элементов $ij$ ($i \geq j+1$), количество таких пар равно $\frac{n(n+1)}{2}+(n-1)$), а также из набора логических операций "И" между результатами сравнений (таких операций будет $\frac{n(n+1)}{2}+(n-1)-1$)

Итого, в сумме получаем $O(n^3)$ - сложность алгоритма на каждой итерации. Если алгоритм остановился на итерации с номером $N$, то общая сложность алгоритма будет равна $N*O(n^3)$

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

---

На рисунке ниже изображен информационный граф QR-алгоритма. Вершины данного графа обозначают следующее:

• QR - операция QR-разложения матрицы
• *_M - операция перемножения матриц
• Triang - операция проверки матрицы на квазитреугольную форму
• Dif - операция проверки различия элементов двух матриц не более чем на некоторое заданное число

Информационный граф QR-алгоритма

Подробные графы операций QR-разложения (Метод Хаусхолдера (отражений), Метод Гивенса (вращений)) и перемножения матриц можно найти в статьях, посвященных этим алгоритмам. Далее рассмотрим информационные графы операций Triang и Dif на примере матрицы размера $5 \times 5$. Графы для других матриц выглядят аналогичным образом.

Информационный граф операции Triang, которая проверяет соответствие квазитреугольной форме

Вершины V соответствуют операции сравнения с 0. Вершины F соответствуют проверке поддиагональных элементов на соответствие квазитреугольной форме, которая была описана в предыдущем разделе. Вершины & соответствуют логической операции "И".

Информационный граф операции Dif, которая проверяет сходимость ненулевых элементов

Вершины -V соответствуют операции вычитания и сравнения с 0. Вершины & соответствуют логической операции "И".

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

---

На макроуровне (который можно увидеть на информационном графе QR-алгоритма) алгоритм не обладает ресурсами параллелизма. Все макрооперации внутри итерации, а также сами итерации выполняются последовательно (за исключением операций Triag и Dif, которые могут выполняться параллельно). Основной ресурс параллелизма заложен отдельно в каждой из макроопераций. На каждой итерации алгоритм имеет следующие параллельные характеристики:

1. QR-разложение матрицы (описание ресурсов параллелизма для алгоритмов QR-разложения можно найти в статьях, посвященных этим алгоритмам)
   • Метод Хаусхолдера (отражений) имеет высоту ярусно-параллельной формы $O(n^2)$ и ширину ярусно-параллельной формы $O(n)$
   • Метод Гивенса (вращений) имеет высоту ярусно-параллельной формы $11n-16$ и ширину ярусно-параллельной формы $O(n^2)$
2. Перемножение двух плотных матриц имеет высоту ярусно-параллельной формы $n$ и ширину ярусно-параллельной формы $n^2$
3. Проверка матрицы на выходе из итерации
   • Проверка матрицы на квазитреугольную форму состоит из одного яруса сравнений для каждого элемента и последующих ярусов, вычисляющих итоговый результат при помощи логической операции "И". Для вычисления итогового результата можно использовать метод сдваивания, который дает высоту ярусно-параллельной формы порядка логарифма от количества элементов, над которыми совершается операция. Таким образов высота ярусно-параллельной формы будет равна $O(log_2n)$. Ширина ярусно-параллельной формы достигается на ярусе ссравнений для каждого элемента и равна $O(n^2)$.
   • Сравнение новой матрицы с предыдущей состоит из одного яруса сравнений для каждого элемента и последующих ярусов, вычисляющих итоговый результат при помощи логической операции "И". Высота ярусно-параллельной формы будет равна $O(log_2n)$. Ширина ярусно-параллельной формы достигается на ярусе сравнений для каждого элемента и равна $O(n^2)$.

Таким образом, основной вклад в высоту ярусно-параллельной формы одной итерации вносит операция QR-разложения матрицы и она будет равна $O(n)$, если использовать метод Гивенса. Ширина ярусно-параллельной формы будет равна $O(n^2)$. Если алгоритм остановился на итерации с номером $N$, то параллельные характеристики для всего алгоритма будут равны $N*O(n)$ для высоты и $O(n^2)$ для ширины ярусно-параллельной формы.

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

---

Входные данные:

квадратная вещественная плотная матрица $A$: $A \in \R^{n \times n}$

Объем входных данных:

$n^2$ (необходимо хранить все элементы матрицы)

Выходные данные:

собственные значения матрицы $A$

Объем выходных данных:

$n$ (квадратная матрица размера $n \times n$ имеет ровно $n$ собственных значений при этом некоторые из них могут быть комплексными)

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

---

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

2 ЧАСТЬ. Программная реализация алгоритма

---

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

2.2.1 Последовательные реализации

1. LAPACK (Linear Algebra PACKage) — библиотека с открытым исходным кодом, содержащая методы для решения основных задач линейной алгебры. Написана на языке Fortran с использованием библиотеки BLAS и является развитием пакета LINPACK. Находится в открытом доступе в соответствии с модифицированной лицензией BSD, в том числе и для коммерческого использования.
   • Полный алгоритм нахождения собственных значений: функция dhseqr
   • QR-разложение: функция dgeqrf
   • Перемножение матриц: функция dgemm
2. ALGLIB - это кросс-платформенная библиотека численного анализа, поддерживающая несколько языков программирования (C++, C#, Pascal, VBA) и несколько операционных систем (Windows, Linux, Solaris). ALGLIB является свободным программным обеспечением, которое использует двойное лицензирование: оно может быть использовано в соответствии с лицензией GPL (версии 2+), а для использования в коммерческих целях необходимо купить отдельную лицензию.
   • Полный алгоритм нахождения собственных значений: подпрограмма RMatrixEVD
   • QR-разложение: подпрограмма rmatrixqr реализует QR-разложение для вещественных матриц, cmatrixqr – для комплексных матриц
   • Перемножение матриц: для перемножения матриц ALGLIB использует соответствующую реализацию библиотеки BLAS
3. Eigen – библиотека шаблонов для линейной алгебры, написанная на языке C++. Eigen – свободно распространяемое программное обеспечение. Начиная с версии 3.1.1, оно лицензируется MPL2, на ранние версии распространяется действие лицензии LGPL3+.
   • Полный алгоритм нахождения собственных значений: модуль Eigenvalues module
   • QR-разложение: модуль QR module
   • Перемножение матриц: операторы * и *=

2.2.2 Параллельные реализации

1. ScaLAPACK (Scalable Linear Algebra PACKage) — библиотека с открытым исходным кодом, включающая в себя подмножество процедур LAPACK, переработанных для использования на MPP-компьютерах. ScaLAPACK разработана с использованием PBLAS и BLACS, и предназначена для вычислений на любом компьютере или кластере поддерживающим MPI или PVM. Библиотека в настоящее время написана на языке Fortran. Находится в открытом доступе в соответствии с модифицированной лицензией BSD, в том числе и для коммерческого использования.
   • Полный алгоритм нахождения собственных значений: функция pdhseqr
   • QR-разложение: функция pdgeqrf
   • Перемножение матриц: функция pdgemm
2. PLAPACK (Parallel Linear Algebra Package) — пакет функций LAPACK для параллельного решения задач линейной алгебры. Пакет функций PLAPACK является альтернативой библиотеке ScaLAPACK. Для осуществления межпроцессорных коммуникаций в PLAPACK использован интерфейс передачи сообщений MPI. При передаче сообщений в PLAPACK в основном используются коллективные операции, такие, как обобщенная передача данных от одного процесса всем процессам (MPI_Scatter), обобщенная передача данных от всех процессов одному процессу (MPI_Gather), широковещательная рассылка (MPI_Bcast) и другие. PLAPACK включает интерфейсы для языков Fortran и С.
   • QR-разложение: функция PLA_QR
   • Перемножение матриц: функция PLA_Gemm