QR-алгоритм для хессенберговой матрицы, используемый в SCALAPACK

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

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

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

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

В описываемом алгоритме исходная матрица хессенбергова. Благодаря инвариантности хессенберговой формы к QR-алгоритму видно, что так как $A_0 = A$ имеет эту форму, то имеют её и все матрицы $A_k$.

1.2 Математическое описание

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

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

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

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

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

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

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

$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}$.

1.2.2 Приёмы ускорения сходимости, используемые в алгоритме

1.2.2.1 Одновременная обработка диагональных блоков и дефляция

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

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

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

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

Рисунок 1. Область действия стартового шага (формирование $1$го выступа) итерации, состоящей из левого преобразования Хаусхолдера, меняющего 3 строки, и правого преобразования Хаусхолдера, меняющего 3 столбца матрицы.

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

При этом сохраняется свойство подобия матриц $A_k$ и $A_{k+1}$: $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$.

Рисунок 2. Область действия шага (вытеснение $i$го выступа) итерации, состоящей из левого преобразования Хаусхолдера, меняющего 3 строки, и правого преобразования Хаусхолдера, меняющего 3 столбца матрицы.

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

Рисунок 3. Область действия последнего шага (вытеснение $n-2$го выступа) итерации, состоящей из левого преобразования Хаусхолдера, меняющего 2 строки, и правого преобразования Хаусхолдера, меняющего 2 столбца матрицы.

1.2.2.3 Неявное выполнение сдвигов

Чтобы не выполнять итерации в комплексной арифметике, выполняют т.н. неявный двойной сдвиг, давно описанный в литературе. При этом используются двусторонние подобные операции Хаусхолдера (отражения) 3го порядка. На первом шаге итерации формируют у хессенберговой формы выступ размера 2х2, а на последующих шагах с их же помощью выступ вытесняют вправо вниз. На последнем шаге выступ размера 1х1.

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

Основная часть алгоритма состоит из выполняемых над диагональными блоками, оставшимися с размерами более, чем 2, очередной итерации. При этом очередная итерация состоит из вычисления и применения к блоку двусторонних подобных преобразований Хаусхолдера (отражения) размерности 3 (последнее - размерности 2), выполняемыми для стартового создания выступа 2х2 и затем для его вытеснения.

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

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