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 для решения проблемы собственных значений у матриц общего вида в настоящее время использует тот подход к выбору сдвигов, который не подтверждён теоретически, но имел лучшие практические показатели на момент создания библиотеки LAPACK. Это выполнение процедуры "двойного сдвига" с величинами, равными собственным значениям последней диагональной клетки размера 2х2 у всех диагональных блоков, больших размером. Поскольку существуют матрицы, для которых этот выбор сдвигов не даёт сходимости, то после 10й итерации при отсутствии уменьшения поддиагональных элементов выполняется некий "особый" сдвиг.

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

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