Уровень алгоритма

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

Материал из Алговики
Перейти к навигации Перейти к поиску
Файл:Symbol confirmed.svgЭта работа успешно выполнена
Преподавателю: в основное пространство, в подстраницу
Данное задание было проверено и зачтено.
Проверено Frolov и VadimVV.



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


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

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

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

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

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

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

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

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

QR-алгоритм

Пусть [math]\displaystyle{ A_0 = A }[/math] — исходная матрица. Для [math]\displaystyle{ k = 1, 2, \ldots }[/math]:

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

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

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

Сходимость QR-алгоритма трактуется как сходимость к нулю поддиагонального блока матрицы [math]\displaystyle{ A_k }[/math].

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

  1. [math]\displaystyle{ 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];
  2. [math]\displaystyle{ \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];
  3. ведущая подматрица порядка [math]\displaystyle{ m }[/math] в [math]\displaystyle{ X^{-1} }[/math] невырожденная.

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

[math]\displaystyle{ A_k = \begin{bmatrix}A_{11}^{(k)} & A_{12}^{(k)} \\ A_{21}^{(k)} & A_{22}^{(k)}\end{bmatrix} }[/math].

Тогда [math]\displaystyle{ A_{21}^{(k)} \rightarrow O }[/math] при [math]\displaystyle{ k \rightarrow \infty }[/math].

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

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

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

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

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

Затем QR-алгоритм применяется к матрице [math]\displaystyle{ A_0 = H }[/math].

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

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

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

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

[math]\displaystyle{ A_0 = A }[/math] — исходная матрица. Для [math]\displaystyle{ k = 1, 2, \ldots }[/math]:

  • [math]\displaystyle{ A_{k-1} - s_kI= Q_kR_k }[/math] (QR-разложение);
  • [math]\displaystyle{ A_k = R_kQ_k + s_kI }[/math].

QR-алгоритм со сдвигами является частным случаем обобщенного QR-алгоритма (QR-алгоритма с мультисдвигами), в котором на каждой итерации используется полином [math]\displaystyle{ f_k }[/math]:

  • [math]\displaystyle{ f_k(A_{k-1}) = Q_kR_k }[/math];
  • [math]\displaystyle{ A_k = Q_k^{-1}A_{k-1}Q_k }[/math].

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

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

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

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

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

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

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

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

Пример реализации на языке 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)

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

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

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

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

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

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

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

Узлы графа:

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

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

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

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

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

При классификации по ширине ярусно-параллельной формы, сложность QR-алгоритма составляет [math]\displaystyle{ O(n^2) }[/math].

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

Входные данные: плотная квадратная матрица [math]\displaystyle{ A }[/math] порядка [math]\displaystyle{ n }[/math].

Объем входных данных: [math]\displaystyle{ n^2 }[/math].

Выходные данные: собственные числа матрицы [math]\displaystyle{ A }[/math].

Объем выходных данных: [math]\displaystyle{ n }[/math].

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

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

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

Особенности реализации последовательного алгоритма

Локальность данных и вычислений

Возможные способы и особенности параллельной реализации алгоритма

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

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

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

 $ 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-алгоритма.

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

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

Динамические характеристики и эффективность реализации алгоритма

Выводы для классов архитектур

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

  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++).

Литература

  1. ↑ 1,0 1,1 Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. — 6-е изд. — М.: БИНОМ. Лаборатория знаний, 2008. — 636 с.
  2. ↑ 2,0 2,1 2,2 Тыртышников Е. Е. Методы численного анализа. — М.: Академия, 2007. — 320 c.
  3. ↑ Воеводин Вл., Жуматий С., Соболев С., Антонов А., Брызгалов П., Никитенко Д., Стефанов К., Воеводин Вад. Практика суперкомпьютера «Ломоносов» // Открытые системы, 2012, N 7, С. 36-39.
Источник — https://xn--80aecvaldz.xn--80ahdrnmvk.xn--p1ai/w/ru/index.php?title=Участник:F-morozov/Нахождение_собственных_чисел_квадратной_матрицы_методом_QR_разложения_(4)&oldid=21128