QR-алгоритм

Задача нахождения собственных значений и собственных векторов для матрицы [math]\displaystyle{ A }[/math] заключается в поиске таких соответствующих друг другу чисел [math]\displaystyle{ \lambda }[/math] и ненулевых векторов [math]\displaystyle{ x }[/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_0 = A }[/math] — исходная матрица. Для [math]\displaystyle{ k = 0, 1, 2, \ldots }[/math] выполняется QR-разложение:

• [math]\displaystyle{ A_{k} = Q_kR_k }[/math], где [math]\displaystyle{ Q_k }[/math] — унитарная (ортогональная) матрица, [math]\displaystyle{ R_k }[/math] — верхняя треугольная матрица, и затем найденные матрицы перемножаются в обратном порядке:
• [math]\displaystyle{ A_{k+1} = R_kQ_k }[/math].

Поскольку [math]\displaystyle{ A_{k+1} = R_kQ_k = Q_k^{*}A_{k}Q_k }[/math], то матрицы [math]\displaystyle{ A_{k+1} }[/math] и [math]\displaystyle{ A_k }[/math] унитарно подобны для любого [math]\displaystyle{ k }[/math]. Поэтому матрицы [math]\displaystyle{ A_1, A_2, \ldots }[/math] унитарно подобны исходной матрице [math]\displaystyle{ A }[/math] и имеют те же собственные значения.

В специальной литературе по численным методам приводится доказательство сходимости получающихся матриц к клеточной правой треугольной матрице с диагональными клетками, размеры которых зависят от размеров канонических ящиков Жордана исходной матрицы. Таким образом, проблема собственных значений матрицы сводится к проблемам собственных значений матриц меньшего размера. QR-алгоритм в узком смысле и есть эта процедура. Однако обычно к нему относят не только её, но и весь комплекс приёмов ускорения этого итерационного метода.

Ускорение сходимости QR-алгоритма

Применение приёмов ускорения даёт QR-алгоритму существенные преимущества перед альтернативными методами решения проблемы собственных значений. Приёмы делятся по своей сути на две группы - приёмы ускорения одной итерации алгоритма и приёмы ускорения сходимости итерационного процесса, то есть уменьшения числа итераций.

Использование хессенберговой формы матрицы

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

Поэтому естественным приёмом ускорения итераций QR-алгоритма является начальное унитарно подобное приведение матрицы к хессенбергову виду. После этого все процедуры 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]. При удачном алгоритме подбора параметров среднее количество итераций получается порядка 5 на одно собственное значение.