Участник:F-morozov/Нахождение собственных чисел квадратной матрицы методом QR разложения (4)
Нахождение собственных чисел квадратной матрицы методом QR-разложения | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(N \times n^3)[/math] |
Объём входных данных | [math]n^2[/math] |
Объём выходных данных | [math]n[/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math][/math] |
Ширина ярусно-параллельной формы | [math][/math] |
Основные авторы описания: Ф. В. Морозов, Н. Ф. Пащенко
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Для решения ряда задач механики, физики, химии требуется получение всех собственных значений (собственных чисел), а иногда и всех собственных векторов некоторых матриц. Эту задачу называют полной проблемой собственных значений[1]. Если рассматриваются матрицы общего вида, порядок которых не больше тысячи (нескольких тысяч), то для вычисления всех собственных значений (и собственных векторов) можно рекомендовать QR-алгоритм. Он был разработан в начале 1960-х годов независимо В. Н. Кублановской (Россия) и Дж. Фрэнсисом (Великобритания)[2].
Нахождение собственных чисел матрицы [math]A[/math] методом QR-разложения (QR-алгоритм) заключается в построении последовательности матриц, сходящейся по форме к клеточному правому треугольному виду. Матрицы [math]A_k[/math] данной последовательности строятся с использованием QR-разложения таким образом, что они подобны между собой и подобны исходной матрице [math]A[/math], поэтому их собственные значения равны. Характеристический многочлен клеточной правой треугольной матрицы равен произведению характеристических многочленов ее диагональных клеток[1]. Его корни — собственные значения матрицы, которые согласно вышесказанному и являются искомыми.
1.2 Математическое описание алгоритма
Известно, что произвольная квадратная матрица может быть представлена в виде произведения унитарной (в вещественном случае ортогональной) и верхней треугольной матриц. Такое разложение называется QR-разложением.
1.2.1 QR-алгоритм
Пусть [math]A_0 = A[/math] — исходная матрица. Для [math]k = 1, 2, \ldots[/math]:
- [math]A_{k-1} = Q_kR_k[/math], где [math]Q_k[/math] — унитарная (ортогональная) матрица, [math]R_k[/math] — верхняя треугольная матрица;
- [math]A_k = R_kQ_k[/math].
Заметим, что [math]A_k = R_kQ_k = Q_k^{-1}A_{k-1}Q_k[/math], то есть матрицы [math]A_k[/math] и [math]A_{k-1}[/math] подобны для любого [math]k[/math]. Таким образом, матрицы [math]A_1, A_2, \ldots[/math] подобны исходной матрице [math]A[/math] и имеют те же собственные значения. При выполнении некоторых условий последовательность [math]\{A_k\}[/math] сходится к верхней треугольной матрице [math]\hat{A}[/math], на диагонали которой расположены искомые собственные значения.
1.2.2 Сходимость QR-алгоритма
Сходимость QR-алгоритма трактуется как сходимость к нулю поддиагонального блока матрицы [math]A_k[/math].
Предположим, что для матрицы [math]A \in \mathbb{C}^{n \times n}[/math] выполнены следующие условия:
- [math]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}[/math];
- [math]\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)[/math];
- ведущая подматрица порядка [math]m[/math] в [math]X^{-1}[/math] невырожденная.
И пусть QR-алгоритм порождает последовательность матриц
- [math]A_k = \begin{bmatrix}A_{11}^{(k)} & A_{12}^{(k)} \\ A_{21}^{(k)} & A_{22}^{(k)}\end{bmatrix}[/math].
Тогда [math]A_{21}^{(k)} \rightarrow O[/math] при [math] k \rightarrow \infty[/math].
Если блок [math]A_{21}^{(k)}[/math] достаточно мал, то он заменяется нулевым, и собственные значения ищутся для диагональных блоков [math]A_{11}^{(k)}[/math] и [math]A_{22}^{(k)}[/math].
Если условия 1-3 выполнены для всех [math]1 \le m \le n - 1[/math], то к нулю сходятся все поддиагональные блоки. А это значит, что все поддиагональные элементы матриц [math]A_k[/math] стремятся к нулю, и диагональные элементы матриц [math]A_k[/math] сходятся к искомым собственным значениям матрицы [math]A[/math][2].
1.2.3 Сокращение затрат на одну итерацию
В общем случае QR-алгоритм сходится медленно - одна QR-итерация для матрицы общего вида требует [math]O(n^3)[/math] арифметических операций. Это очень большие затраты, даже если итераций не очень много (обычно их требуется не больше 5 на каждое собственное значение). Для решения данной проблемы используют тот факт, что с помощью отражений или вращений матрицу [math]A[/math] можно привести к унитарно подобной верхней хессенберговой матрице [math]H = (h_{ij})[/math] (то есть [math]h_{ij} = 0[/math] при [math]i \gt j + 1[/math], [math]H = PAP^*[/math], где [math]P[/math] — произведение конечного числа отражений или вращений). Затем QR-алгоритм применяется к матрице [math]A_0 = H[/math]. Приведение к хессенберговой форме требует [math]O(n^3)[/math] операций, но после него любая QR-итерация будет выполняться за [math]O(n^2)[/math] операций. Сокращение затрат на одну итерацию объясняется инвариантностью хессенберговой формы по отношению к QR-итерациям[2].
1.3 Вычислительное ядро алгоритма
Основное время работы алгоритма приходится на вычисление QR-разложения на итерациях алгоритма.
Для получения QR-разложения могут использоваться метод Хаусхолдера, метод Гивенса или процесс ортогонализации Грамма-Шмидта.
При использовании хессенберговой формы, приведение к ней имеет сложность, сравнимую с последующими итерациями QR-алгоритма.
1.4 Макроструктура алгоритма
Основные элементы алгоритма:
- (Опционально) Приведение к хессенберговой форме
- QR-разложение
- Умножение плотных матриц
- Проверка, имеет ли матрица верхнюю треугольную форму
1.5 Схема реализации последовательного алгоритма
if use_hessenberg_form: A = to_hessenberg_form(A) while not is_upper_triangular(A): Q, R = QR_decomposition(A) A = R * Q return diagonal_items(A)
1.6 Последовательная сложность алгоритма
QR-разложение имеет сложность [math]O(n^3)[/math] арифметических операций, сложность перемножения квадратных матриц равна [math]n^3 + O(n^2)[/math]. Таким образом, каждая итерация имеет кубическую сложность и, если алгоритм сходится через [math]N[/math] итераций, общая сложность составит [math]O(N \times n^3)[/math] арифметических операций.
Сложность построения хессенберговой формы составляет [math]O(n^3)[/math], QR-разложение хессенберговой матрицы имеет сложность [math]O(n^2)[/math]. Общая сложность QR-алгоритма с предварительным приведением матрицы к хессенберговой форме равна [math]O(n^3 + N \times n^2)[/math] арифметических операций, где [math]N[/math] — число итераций алгоритма.
1.7 Информационный граф
1.8 Ресурс параллелизма алгоритма
1.9 Входные и выходные данные алгоритма
Входные данные: плотная квадратная матрица [math]A[/math] порядка [math]n[/math].
Объем входных данных: [math]n^2[/math].
Выходные данные: собственные числа матрицы [math]A[/math].
Объем выходных данных: [math]n[/math].
1.10 Свойства алгоритма
2 Программная реализация алгоритма
2.1 Масштабируемость алгоритма и его реализации
2.2 Существующие реализации алгоритма
- LAPACK (последовательная). Имеется реализация в Intel MKL.
- ?geev - нахождение собственных значений матрицы.
- p?gehrd - приведение к хессенберговой форме.
- p?lahqr - нахождение собственных значений матрицы в хессенберговой форме.
- ScaLAPACK (параллельная). Имеется реализация в Intel MKL.
- ?gehrd- приведение к хессенберговой форме.
- ?hseqr - нахождение собственных значений матрицы в хессенберговой форме.
- numpy.linalg.eigvals (Python) - использует ?geev из LAPACK.
- Eigen (C++)