QR-алгоритм, используемый в SCALAPACK

QR-алгоритм, используемый в SCALAPACK - алгоритм, использующий все проверенные на настоящее время приёмы ускорения QR-алгоритма для несимметричных вещественных матриц. Включён своими частями в разные подпрограммы пакета SCALAPACK. Состоит из двух основных частей: ортогонально подобного приведения матрицы к хессенберговому виду и 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-алгоритму тем, что в случае, если [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] операций. Это позволяет экономить довольно большие вычислительные ресурсы.

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

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

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

Особенности вещественного варианта 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].

Поэтому выбор сдвигов в алгоритме учитывает, что конечные матрицы будут представленного вида. Чтобы не выполнять итерации в комплексной арифметике, выполняют т.н. неявный двойной сдвиг, давно описанный в литературе.

Обнуление поддиагональных элементов

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