QR-алгоритм

Авторы описания: Смирнова А.С. (Раздел 2), Фролов А.В. (общая редактура и правка)

Задача нахождения собственных значений и собственных векторов для матрицы [math]A[/math] заключается в поиске таких соответствующих друг другу чисел [math]\lambda[/math] и ненулевых векторов [math]x[/math], которые удовлетворяют уравнению [math]Ax=\lambda x[/math], при этом числа [math]\lambda[/math] называются собственными значениями, а вектора [math]x[/math] - собственными векторами.

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

1 Общее описание метода

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

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

2 Математическое описание базового метода

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

Пусть [math]A_0 = A[/math] — исходная матрица. Для [math]k = 0, 1, 2, \ldots[/math] выполняется QR-разложение:

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

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

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

2.1 Особенности вещественного несимметричного варианта QR-алгоритма

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

[math]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].

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

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

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

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

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

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

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

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

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

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

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

[math]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]\nu_k[/math] осуществляется в зависимости от элементов матрицы [math]A_k[/math]. При алгоритме подбора параметров, известном, как правило Фрэнсиса, среднее количество итераций для большинства матриц получается порядка 2 на одно собственное значение. Существуют, однако, матрицы, для которых правило Фрэнсиса не даёт сходимости, поэтому если по истечении 10 итераций с обычным сдвигом не получилось сходимости, применяют так называемые "особые сдвиги". Для симметрического же случая используется стратегия Уилкинсона.

Указанные стратегии выбора сдвигов используется в программах библиотек LAPACK и SCALAPACK. Но они не являются наилучшими, поскольку для выбора новых сдвигов нужно провести предшествующие итерации до конца. Поэтому в публикациях XXI века уже появились описания т.н. многосдвиговых итераций. В этой стратегии выбор сдвигов также выполняется с помощью собственных чисел подматриц из нижнего правого угла обрабатываемых блоков, но их размерность уже берётся больше, чем 2.

3.3 Неявное проведение итераций

Q-теорема помогает не вычислять в явном виде требуемые матрицы с их последующим перемножением, а выполнять эти операции неявно. При этом в соответствии с особенностями вещественного случая несимметричных матриц удаётся выполнять т.н. неявный двойной сдвиг с использованием двух комплексно сопряжённых сдвигов (например, из правила Фрэнсиса), оставаясь в рамках вещественной арифметрики, с помощью выполнения двусторонними преобразованиями Хаусхолдера т.н. вытеснения выступа. В вещественном случае аналогично (но с преимуществами симметрии) выполняется неявный одинарный сдвиг, вытеснением выступа двусторонними преобразованиями вращения (Гивенса).

3.4 Обнуление поддиагональных элементов (дефляция)

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

4 QR-алгоритм в различных вариантах

В классическом (без приёмов ускорения) виде QR-алгоритм не применяется из-за медленности. Варианты метода, в основном, используются для матриц разного вида:

4.1 Метод для матриц общего вида

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

4.2 Метод для симметричных и эрмитовых матриц

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