Участник:F-morozov/Нахождение собственных чисел квадратной матрицы методом QR разложения (4)

Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Frolov и VadimVV.

Нахождение собственных чисел квадратной матрицы методом QR-разложения
Последовательный алгоритм
Последовательная сложность	$O(n^3 + N \times n^2)$
Объём входных данных	$n^2$
Объём выходных данных	$n$
Параллельный алгоритм
Высота ярусно-параллельной формы	$N \times O(n)$
Ширина ярусно-параллельной формы	$O(n^2)$

Основные авторы описания: Ф. В. Морозов, Н. Ф. Пащенко

Описание составлено совместно, точно выделить вклад отдельных авторов невозможно.

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

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

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

Нахождение собственных чисел матрицы $A$ методом QR-разложения (QR-алгоритм) заключается в построении последовательности матриц, сходящейся по форме к клеточному правому треугольному виду. Матрицы $A_k$ данной последовательности строятся с использованием QR-разложения таким образом, что они подобны между собой и подобны исходной матрице $A$, поэтому их собственные значения равны. Характеристический многочлен клеточной правой треугольной матрицы равен произведению характеристических многочленов ее диагональных клеток. Его корни — собственные значения матрицы, которые согласно вышесказанному и являются искомыми.

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

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

1.2.1 QR-алгоритм

Пусть $A_0 = A$ — исходная матрица. Для $k = 1, 2, \ldots$:

• $A_{k-1} = Q_kR_k$, где $Q_k$ — унитарная (ортогональная) матрица, $R_k$ — верхняя треугольная матрица;
• $A_k = R_kQ_k$.

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

1.2.2 Сходимость QR-алгоритма

Сходимость QR-алгоритма трактуется как сходимость к нулю поддиагонального блока матрицы $A_k$.

Предположим, что для матрицы $A \in \mathbb{C}^{n \times n}$ выполнены следующие условия:

1. $A = X \Lambda X^{-1}, \quad \Lambda = \begin{bmatrix} \Lambda_1 & 0 \\ 0 & \Lambda_2 \end{bmatrix}, \quad \Lambda_1 \in \mathbb{C}^{m \times m}, \quad \Lambda_2 \in \mathbb{C}^{r \times r}$;
2. $\left | \lambda_1 \right | \ge \ldots \ge \left | \lambda_m \right | \gt \left | \lambda_{m+1} \right | \ge \ldots \ge \left | \lambda_{m+r} \right | \gt 0, \quad \{ \lambda_1, \ldots, \lambda_m \} = \lambda(\Lambda_1), \quad \{ \lambda_{m+1}, \ldots, \lambda_{m+r} \} = \lambda(\Lambda_2)$;
3. ведущая подматрица порядка $m$ в $X^{-1}$ невырожденная.

И пусть QR-алгоритм порождает последовательность матриц

$A_k = \begin{bmatrix}A_{11}^{(k)} & A_{12}^{(k)} \\ A_{21}^{(k)} & A_{22}^{(k)}\end{bmatrix}$.

Тогда $A_{21}^{(k)} \rightarrow O$ при $k \rightarrow \infty$.

Если блок $A_{21}^{(k)}$ достаточно мал, то он заменяется нулевым, и собственные значения ищутся для диагональных блоков $A_{11}^{(k)}$ и $A_{22}^{(k)}$.

Если условия 1-3 выполнены для всех $1 \le m \le n - 1$, то к нулю сходятся все поддиагональные блоки. А это значит, что все поддиагональные элементы матриц $A_k$ стремятся к нулю, и диагональные элементы матриц $A_k$ сходятся к искомым собственным значениям матрицы $A$.

1.2.3 Сокращение затрат на одну итерацию

Одна QR-итерация для матрицы общего вида требует $O(n^3)$ арифметических операций. Это очень большие затраты, даже если итераций не очень много (обычно их требуется не больше 5 на каждое собственное значение). Для решения данной проблемы используют тот факт, что с помощью отражений или вращений матрицу $A$ можно привести к унитарно подобной верхней хессенберговой (почти треугольной) матрице:

• $H = (h_{ij}), \quad h_{ij} = 0$ при $i \gt j + 1$;
• $H = PAP^*$, где $P$ — произведение конечного числа отражений (или вращений).

Затем QR-алгоритм применяется к матрице $A_0 = H$.

Приведение к хессенберговой форме требует $O(n^3)$ операций, но после него любая QR-итерация будет выполняться за $O(n^2)$ операций. Сокращение затрат на одну итерацию объясняется инвариантностью хессенберговой формы по отношению к QR-итерациям.

1.2.4 Уменьшение числа итераций

Скорость сходимости алгоритма зависит от отношения модулей собственных значений: чем больше отношение $|\lambda_{m+1}| / |\lambda_m|$, тем медленнее убывает поддиагональный $(n-m) \times m$ блок. Для увеличения скорости сходимости алгоритма можно перейти от матрицы $A$ к матрице $(A - sI)$, где $s \in \mathbb{C}$, а $I$ — единичная матрица, в надежде, что $|\lambda_{m+1} - s| / |\lambda_m - s| \ll |\lambda_{m+1}| / |\lambda_m|$.

Такой подход называется QR-алгоритмом со сдвигами и имеет следуюший вид:

$A_0 = A$ — исходная матрица. Для $k = 1, 2, \ldots$:

• $A_{k-1} - s_kI= Q_kR_k$ (QR-разложение);
• $A_k = R_kQ_k + s_kI$.

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

• $f_k(A_{k-1}) = Q_kR_k$;
• $A_k = Q_k^{-1}A_{k-1}Q_k$.

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

Основное время работы алгоритма приходится на вычисление QR-разложения на итерациях алгоритма.

Для получения QR-разложения могут использоваться метод Хаусхолдера, метод Гивенса или процесс ортогонализации Грама-Шмидта.

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

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

Основные элементы алгоритма:

• (опционально) приведение к хессенберговой форме;
• QR-разложение;
• умножение плотных матриц.

1.5 Схема реализации последовательного алгоритма

Пример реализации на языке Python:

def get_eigenvalues(A, use_hessenberg_form):
    """Вычисление собственных значений матрицы методом QR-разложения

    Параметры:
    A -- квадратная матрица
    use_hessenberg_form -- использовать ли приведение к хессенберговой форме
    """
    if use_hessenberg_form:
        A = to_hessenberg_form(A) #Матрица приводится к хессенберговой форме
    while not is_upper_triangular(A): #Является ли матрица верхней треугольной
        Q, R = QR_decomposition(A) #Производится QR-разложение матрицы
        A = R * Q
    #Собственными значениями являются диагональные элементы полученной матрицы
    return diagonal_items(A)

1.6 Последовательная сложность алгоритма

Рассмотрим сложность отдельных элементов алгоритма:

• QR-разложение имеет сложность $O(n^3)$ арифметических операций;
• сложность перемножения квадратных матриц равна $n^3 + O(n^2)$;
• для проверки, является ли матрица верхней треугольной, необходимо $O(n^2)$ операций.

Таким образом, каждая итерация имеет кубическую сложность и, если алгоритм сходится через $N$ итераций, общая сложность составит $N \times O(n^3)$ арифметических операций.

Сложность построения хессенберговой формы составляет $O(n^3)$, одна итерация алгоритма для хессенберговой матрицы имеет сложность $O(n^2)$. Общая сложность QR-алгоритма с предварительным приведением матрицы к хессенберговой форме равна $O(n^3 + N \times n^2)$ арифметических операций, где $N$ — число итераций алгоритма.

1.7 Информационный граф

Рисунок 1. Информационный граф.

Узлы графа:

• $QR$ — QR-разложение;
• $\times$ — перемножение матриц;
• $isUT$ — проверка, является ли матрица верхней треугольной.

1.8 Ресурс параллелизма алгоритма

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

• параллельная сложность QR-разложения составляет $O(n)$ при использовании метода Гивенса или $O(n^2)$ при использовании метода Хаусхолдера;
• параллельная сложность умножения матриц равна $O(n)$;
• для проверки, является ли матрица верхней треугольной, достаточно одного яруса.

Таким образом, при классификации по высоте ярусно-параллельной формы сложность равна $N \times O(n)$ или $N \times O(n^2)$ в зависимости от используемого метода QR-разложения.

При классификации по ширине ярусно-параллельной формы, сложность QR-алгоритма составляет $O(n^2)$.

1.9 Входные и выходные данные алгоритма

Входные данные: плотная квадратная матрица $A$ порядка $n$.

Объем входных данных: $n^2$.

Выходные данные: собственные числа матрицы $A$.

Объем выходных данных: $n$.

1.10 Свойства алгоритма

• Соотношение последовательной и параллельной сложности является квадратичным или линейным в зависимости от метода QR-разложения.
• Вычислительная мощность алгоритма равна $N \times O(n)$.
• Число итераций $N$ заранее неизвестно, то есть алгоритм недетерминирован.
• Если рассмотренные в разделе «Сходимость QR-алгоритма» условия 1-3 не выполнены, то, возможно, ни один из поддиагональных блоков не сходится к нулю. Это означает, что в общем случае QR-алгоритм не обязан сходиться. Однако внеся в элементы матрицы сколь угодно малые возмущения, всегда можно удовлетворить указанным условиям.
• Скорость сходимости алгоритма зависит от отношения собственных значений, чем ближе по модулю собственные значения, тем медленее скорость сходимости.
• Если матрица $A$ эрмитова, то хессенбергова матрица $H$ будет трехдиагональной, и сложность одной QR-итерации составит $O(n)$ операций.

2 Программная реализация алгоритма

2.4 Масштабируемость алгоритма и его реализации

Исследование проводилось на узлах раздела test суперкомпьютера "Ломоносов" Суперкомпьютерного комплекса Московского университета.

Исследовалась реализация QR-алгоритма из библиотеки Intel MKL версии 11.2.0 (функции LAPACKE_sgehrd и LAPACKE_shseqr). Была написана программа на языке C, использовался компилятор Intel C++. Пример компиляции программы и запуска для матрицы размера $8000 \times 8000$:

 $ mpicc qr.c -o _scratch/qr -L${MKLROOT}/lib/intel64 -lmkl_rt -lpthread -lm -ldl
 $ sbatch -n 1 -p test run _scratch/qr 8000

Оценивалось время вычисления собственных чисел произвольной квадратной матрицы для различных значений порядка матрицы (Matrix size) и числа потоков (Threads). Ввиду того, что количество операций зависит от числа итераций алгоритма, оценить производительность не удается.

Набор и границы значений изменяемых параметров запуска реализации алгоритма:

• число потоков [1 : 8];
• размер матрицы [2000 : 8000] с шагом 500.

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

Рисунок 2. Производительность QR-алгоритма.

Из рисунка видно, что алгоритм хорошо масштабируется. Время выполнения в зависимости от размера задачи растет значительно быстрее при использовании одного потока по сравнению с использованиемм нескольких потоков. Увеличение числа потоков при фиксированном размере задачи приводит к уменьшению времени выполнения. Для матрицы размера $8000 \times 8000$ время работы однопоточной реализации более чем в два раза превосходит время работы параллельной.

Исходный код программы на языке C и скрипт для обработки результатов доступны на github.

2.7 Существующие реализации алгоритма

1. LAPACK (Fortran, C). Имеется реализация в Intel MKL:
   • ?geev - нахождение собственных значений матрицы;
   • ?gehrd - приведение к хессенберговой форме;
   • ?hseqr - нахождение собственных значений матрицы в хессенберговой форме.
2. ScaLAPACK (Fortran, C). Имеется реализация в Intel MKL:
   • p?gehrd- приведение к хессенберговой форме;
   • p?hseqr - нахождение собственных значений матрицы в хессенберговой форме.
3. numpy.linalg.eigvals (Python) - использует ?geev из LAPACK.
4. Eigen (C++).
5. IBM ESSL (Fortran, C, C++).