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

Файл:Symbol confirmed.svg

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

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

Основные авторы описания:

• Смирнова А.С. - описание теоретической части
• Киямова А. - описание программистской части

Содержательно все пункты обсуждались авторами вместе, и их вклад был равным.

---

Задача нахождения собственных значений и собственных векторов для матрицы [math]\displaystyle{ A }[/math] заключается в поиске таких чисел [math]\displaystyle{ \lambda }[/math], которые удовлетворяют уравнению:

[math]\displaystyle{ Ax=\lambda x }[/math], при этом, числа [math]\displaystyle{ \lambda }[/math] называются собственными значениями, а вектора [math]\displaystyle{ x }[/math] - собственными векторами.

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

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

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

---

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

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

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

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

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

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

---

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

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

[math]\displaystyle{ A=QR }[/math], где [math]\displaystyle{ Q }[/math] - ортогональная матрица ([math]\displaystyle{ Q^T=Q^{-1} }[/math]), [math]\displaystyle{ R }[/math] - верхняя треугольная матрица.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Заметим, что [math]\displaystyle{ 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 }[/math].

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

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

Предположим, что для [math]\displaystyle{ \forall m }[/math] выполнены следующие условия:

1. [math]\displaystyle{ A=X\Lambda X^{-1} }[/math], где [math]\displaystyle{ \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} }[/math].

2. [math]\displaystyle{ \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 }[/math], где [math]\displaystyle{ \{\lambda_1,...,\lambda_m\} = \lambda(\Lambda_1), \{\lambda_{m+1},...,\lambda_{m+r}\} = \lambda(\Lambda_2) }[/math].

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

Тогда при [math]\displaystyle{ k \rightarrow \infty }[/math] последовательность матриц [math]\displaystyle{ A_k }[/math] сходится к матрице с верхней треугольной формой.

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

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

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

[math]\displaystyle{ 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} }[/math].

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

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

Основная статья: Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме

Матрицей, имеющей форму Хессенберга, называется такая матрица, у которой все элементы, находящиеся ниже первой поддиагонали, равны нулю ([math]\displaystyle{ a_{ij}=0 }[/math] при [math]\displaystyle{ i\lt j+1 }[/math]). Пример такой матрицы приведен ниже:

[math]\displaystyle{ \begin{bmatrix} 1& 2& 3& 4\\ 2& 5& 6& 7\\ 0& 3& 8& 9\\ 0& 0& 4& 1 \end{bmatrix} }[/math]

Любую матрицу [math]\displaystyle{ A }[/math] можно привести к подобной ей матрице [math]\displaystyle{ A_H }[/math], имеющей форму Хессенберга (в силу подобия данные матрицы будут иметь одинаковые собственные значения). Наличие нулевых элементов в данной матрице позволяет ускорить процесс QR-разложения, причем данное ускорение будет иметь место на каждой итерации алгоритма, т.к. матрица с формой Хессенберга инвариантна относительно QR-итерации. Ускорение можно получить за счет использования модифицированного метода Гивенса(вращений) QR-разложения, который из-за изначального наличия нулевых элементов будет состоять не из [math]\displaystyle{ \frac{n(n-1)}{2} }[/math] шагов, а из [math]\displaystyle{ n-1 }[/math] шагов (будет происходить домножение только на те матрицы вращения [math]\displaystyle{ T_{ij} }[/math], у которых [math]\displaystyle{ i=j+1 }[/math]).

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

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

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

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

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

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

• На каждом шаге подбирается число [math]\displaystyle{ \nu_k }[/math] и ищется следующее QR-разложение: [math]\displaystyle{ A_k-\nu_kE=Q_kR_k }[/math].
• Вычисляется матрица [math]\displaystyle{ A_{k+1} }[/math]: [math]\displaystyle{ A_{k+1} = R_kQ_k+\nu_kE }[/math].

При этом сохраняется свойство подобия матриц [math]\displaystyle{ A_k }[/math] и [math]\displaystyle{ A_{k+1} }[/math]:

[math]\displaystyle{ A_{k+1} = R_kQ_k+\nu_kE=Q_{k}^{T}Q_kR_kQ_k+\nu_kE=Q_{k}^{T}(A_k-\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-Q_{k}^{T}(\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-\nu_kE+\nu_kE=Q_{k}^{T}A_kQ_k }[/math].

В качестве параметра [math]\displaystyle{ \nu_k }[/math] может быть взят последний диагональный элемент матрицы [math]\displaystyle{ A_k }[/math].

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

---

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

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

1. Вычисление QR-разложения матрицы: [math]\displaystyle{ A_k=Q_kR_k }[/math]. Существует несколько методов решения данной задачи:
   • Метод Хаусхолдера (отражений): вычислительное ядро данного алгоритма состоит из операций скалярного произведения, необходимых для вычисления матрицы отражения, и из операций скалярного произведения, необходимых для пересчета матрицы на каждом шаге.
   • Метод Гивенса (вращений): вычислительное ядро данного алгоритма состоит из численных арифметических операций, необходимых для вычисления параметров матрицы вращения, и из операций (эквивалентных перемножению комплексных чисел), необходимых для пересчета матрицы на каждом шаге.
2. Перемножение двух плотных матриц: [math]\displaystyle{ A_{k+1}=R_kQ_k }[/math].

Модифицированный алгоритм

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

1. Приведение изначальной матрицы к форме Хессенберга: вычислительное ядро данного алгоритма состоит из операций скалярного произведения, необходимых для вычисления матрицы отражения, и из операций скалярного произведения, необходимых для пересчета матрицы на каждом шаге.
2. Вычисление QR-разложения матрицы: [math]\displaystyle{ A_k-\nu_kE=Q_kR_k }[/math] при помощи модифицированного метода Гивенса (вращений). Описание вычислительного ядра не отличается от приведенного в описании базового алгоритма, единственная разница заключается в значительно меньшем количестве матриц вращения.
3. Перемножение двух плотных матриц: [math]\displaystyle{ A_{k+1}-\nu_kE=R_kQ_k }[/math].

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

---

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

1. Вычисление QR-разложения матрицы: [math]\displaystyle{ A_k=Q_kR_k }[/math].
2. Перемножение двух плотных матриц: [math]\displaystyle{ A_{k+1}=R_kQ_k }[/math].

В случае модифицированного QR-алгоритма появляется еще макрооперация приведения матрицы к форме Хессенберга.

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

---

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

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

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
}

Модифицированный алгоритм

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

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

curA = A
nextA = A
while ( not (triangular(nextA) & difference(curA,nextA,eps)) )
{
    curA = nextA
    nu = computeNu(curA)
    shiftA = curA - nu*E
    findQRdecomposition(shiftA,Q,R)
    nextA = R*Q + nu*E
}

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

---

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

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

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

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

Модифицированный алгоритм

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

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

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

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

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

---

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

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

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

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

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

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

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

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

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

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

Модифицированный алгоритм

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

• Shift - операция сдвига
• H - операция приведения матрицы к форме Хессенберга

Рисунок 5. Информационный граф модифицированного QR-алгоритма

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

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

---

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

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

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

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

Модифицированный алгоритм

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

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

Таким образом, итоговые характеристики всего модифицированного алгоритма будут следующие: высота ярусно-параллельной формы [math]\displaystyle{ O(n^2)+N*O(n) }[/math] и ширина ярусно-параллельной формы [math]\displaystyle{ O(n^2) }[/math].

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

---

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

квадратная вещественная плотная матрица [math]\displaystyle{ A }[/math]: [math]\displaystyle{ A \in \R^{n \times n} }[/math].

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

[math]\displaystyle{ n^2 }[/math] (необходимо хранить все элементы матрицы).

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

собственные значения матрицы [math]\displaystyle{ A }[/math].

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

[math]\displaystyle{ n }[/math] (квадратная матрица размера [math]\displaystyle{ n \times n }[/math] имеет ровно [math]\displaystyle{ n }[/math] собственных значений при этом некоторые из них могут быть комплексными).

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

---

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

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

---

---

---

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

---

Масштабируемость алгоритма (рис.7) исследовалась на основе реализации алгоритма в библиотеке ScaLAPACK v2.0.0: функция pdhseqr. Время работы было замерено для различных размеров матрицы, а также различного числа процессоров.

Рисунок 7. Масштабируемость QR-алгоритма

При увеличении числа процессоров время работы программы увеличивается при малом размере матрицы, распараллеливание алгоритма неэффективно. Однако при больших размерах матрицы (более [math]\displaystyle{ 10^3\times10^3 }[/math]) время работы параллельной реализации значительно меньше последовательной.

Код программы

Все вычисления были проведены на суперкомпьютере Regatta. Библиотека ScaLAPACK v2.0.0 была установлена вручную. Компиляция проводилась с использованием компилятора mpicxx с опциями -lscalapack и -llapack с указанием полного пути до используемых библиотек.

---

---

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

---

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

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.
   • Перемножение матриц: операторы * и *=.

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

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.