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

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

Материал из Алговики
Перейти к навигации Перейти к поиску
[досмотренная версия][досмотренная версия]
(Новая страница: «{{level-a}}»)
 
 
(не показано 13 промежуточных версий 3 участников)
Строка 1: Строка 1:
 
{{level-a}}
 
{{level-a}}
 +
 +
'''QR-алгоритм, используемый в SCALAPACK''' - алгоритм, использующий все проверенные на настоящее время приёмы ускорения [[QR-алгоритм|QR-алгоритма]] для несимметричных ''вещественных'' матриц<ref>Алгоритм для неэрмитовых комплексных матриц в SCALAPACK отсутствует, есть только их приведение к хессенберговой форме</ref>. Включён своими частями в разные подпрограммы пакета SCALAPACK<ref name="SCALAEig">http://www.netlib.org/scalapack/slug/node63.html#SECTION04335100000000000000</ref>. Состоит из двух основных частей: [[Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме|ортогонально подобного приведения матрицы к хессенберговому виду]] и [[QR-алгоритм для хессенберговой матрицы, используемый в SCALAPACK|QR-итераций со сдвигами для хессенберговой матрицы]].
 +
 +
== Свойства и структура алгоритма ==
 +
 +
=== Общее описание алгоритма ===
 +
 +
Суть базового QR-алгоритма заключается в итерационном приведении матрицы <math>A</math> к некоторой унитарно подобной ей матрице <math>A_N</math> при помощи QR-разложения. Матрица <math>A_N</math> является правой верхней треугольной (или блочно треугольной) матрицей, а значит ее диагональ содержит собственные значения (в блочном случае собственные значения матрицы являются собственными значениями диагональных блоков). В силу подобия матриц <math>A</math> и <math>A_N</math> их наборы собственных значений совпадают. Таким образом, задача поиска собственных значений матрицы <math>A</math> сводится к задаче выведения матрицы <math>A_N</math> и поиска собственных значений для нее.
 +
 +
=== Математическое описание ===
 +
 +
Известно<ref name="VOLA">Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.</ref>, что произвольная квадратная матрица может быть представлена в виде произведения унитарной (''в вещественном случае ортогональной'') и верхней треугольной матриц. Такое разложение называется [[QR-разложения плотных неособенных матриц|QR-разложением]].
 +
 +
Пусть <math>A_0 = A</math> — исходная матрица.
 +
Для <math>k = 0, 1, 2, \ldots</math> выполняется QR-разложение:
 +
* <math>A_{k} = Q_kR_k</math>, где <math>Q_k</math> — унитарная (''ортогональная'') матрица, <math>R_k</math> — верхняя треугольная матрица, и затем найденные матрицы перемножаются в обратном порядке:
 +
* <math>A_{k+1} = R_kQ_k</math>.
 +
 +
Поскольку <math>A_{k+1} = R_kQ_k = Q_k^{*}A_{k}Q_k</math>, то матрицы <math>A_{k+1}</math> и <math>A_k</math> унитарно подобны для любого <math>k</math>. Поэтому матрицы <math>A_1, A_2, \ldots</math> унитарно подобны исходной матрице <math>A</math> и имеют те же собственные значения.
 +
 +
Существуют доказательства<ref name="VOLA" /> сходимости получающихся матриц к клеточной правой треугольной матрице с диагональными клетками, размеры которых зависят от размеров канонических ящиков Жордана исходной матрицы. Таким образом, проблема собственных значений матрицы сводится к проблемам собственных значений матриц меньшего размера.
 +
 +
==== Особенности вещественного варианта QR-алгоритма ====
 +
 +
Если вещественная матрица <math>A</math> имеет различные вещественные собственные значения, то QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения, и тогда алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка. Блоки 1-го порядка будут содержать различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений<ref name = "Dem"> Деммель Д. Вычислительная линейная алгебра. – 2001. - С.261-264. </ref><ref>R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems".</ref>.
 +
 +
<math>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>.
 +
 +
 +
==== Приёмы ускорения сходимости, используемые в алгоритме ====
 +
 +
Применение приёмов ускорения является существенной частью рассматриваемого здесь алгоритма. Приёмы делятся на две группы - приёмы ускорения одной итерации алгоритма и приёмы ускорения сходимости итерационного процесса, то есть уменьшения числа итераций.
 +
 +
===== Использование хессенберговой формы матрицы =====
 +
 +
Хессенбергова (почти треугольная) форма матрицы интересна в приложении к QR-алгоритму тем, что в случае, если <math>A_0 = A</math> имеет эту форму, то имеют её и все матрицы <math>A_k</math>. Если теперь посмотреть на методы [[QR-разложения плотных неособенных матриц]], то видно, что выполнение разложения <math>A_k = Q_kR_k</math> в случае плотной неособенной <math>A_k</math>требует (в последовательном варианте) <math>O(n^3)</math> операций, как и последующий процесс вычисления <math>A_{k+1} = R_kQ_k</math>. В случае же, если матрицы хессенберговы, то как [[Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме|QR-разложения плотных хессенберговых матриц]], так и вычисления <math>A_{k+1} = R_kQ_k</math> потребуют уже по <math>O(n^2)</math> операций. Это позволяет экономить довольно большие вычислительные ресурсы.
 +
 +
Кроме этого, известно, что в случае, если хессенбергова матрица содержит на нижней кодиагонали только ненулевые элементы, то все кратные собственные значения будут соответствовать одному ящику Жордана<ref name="VOLA" />. Это означает, что каждое собственное значение с геометрической кратностью выше 1 даст после точного приведения матрицы к хессенбергову виду нули на нижней кодиагонали, что сразу даёт возможность параллельно работать с несколькими диагональными блоками.
 +
 +
Поэтому начальное [[Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме|ортогонально подобное приведение матрицы к хессенбергову виду]] является органически необходимой стартовой частью алгоритма<ref name="SCALAEig" />, несмотря на то, что оно выделено в пакете SCALAPACK в отдельную подпрограмму.
 +
 +
===== Сдвиги QR-алгоритма =====
 +
 +
QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости<ref name="bahvalov">Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.</ref>. Пусть у нас есть матрица <math>A_k</math>, тогда процесс перехода к матрице <math>A_{k+1}</math> выглядит следующим образом:
 +
 +
* На каждом шаге подбирается число <math>\nu_k</math> и ищется следующее QR-разложение: <math>A_k-\nu_kE=Q_kR_k</math>.
 +
* Вычисляется матрица <math>A_{k+1}</math>: <math>A_{k+1} = R_kQ_k+\nu_kE</math>.
 +
 +
При этом сохраняется свойство подобия матриц <math>A_k</math> и <math>A_{k+1}</math>:
 +
 +
<math>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>\nu_k</math> в пакете SCALAPACK осуществляется в зависимости от элементов матрицы <math>A_k</math>. При алгоритме подбора параметров, известном, как ''правило Фрэнсиса''<ref name = "Dem" />, в качестве сдвигов выбираются для каждого из диагональных блоков собственные значения их нижних правых клеток размером 2х2, и среднее количество итераций для большинства матриц получается порядка 2 на одно собственное значение. Существуют, однако, матрицы, для которых правило Фрэнсиса не даёт сходимости, поэтому если по истечении 10 итераций с обычными сдвигами не получилось сходимости, применяют так называемые "особые сдвиги".
 +
 +
===== Неявное выполнение сдвигов =====
 +
 +
Чтобы не выполнять итерации в комплексной арифметике, выполняют т.н. неявный двойной сдвиг, давно описанный в литературе<ref name="VOLA" /><ref name = "Dem" />. При этом используются двусторонние подобные операции Хаусхолдера (отражения) 3го порядка. На первом шаге итерации формируют у хессенберговой формы выступ размера 2х2, а на последующих шагах с их же помощью выступ вытесняют вправо вниз.
 +
 +
===== Обнуление поддиагональных элементов (дефляция) и распараллеливание по блокам =====
 +
 +
После того, как под диагональю какой-либо элемент становится менее порогового значения, его можно считать равным нулю. Это даёт возможность применить разделение диагональных блоков матрицы, соседствующих с ним, и раздельное вычисление их собственных чисел. Этот приём называется '''дефляцией'''.
 +
Благодаря обнулению поддиагональных элементов далее нахождение собственных чисел разделённых этими нулями хессенберговых блоков, расположенных на диагонали блочно-треугольной формы, можно производить раздельно. Если в библиотеке LAPACK, предшествующей библиотеке SCALAPACK, на каждой из итераций обрабатывался текущий диагональный блок максимальных размеров, то в алгоритме, используемом в SCALAPACK, все блоки размера больше, чем 3, обрабатываются на очередной итерации параллельно, что позволило ускорить процесс.
 +
 +
=== Вычислительное ядро алгоритма ===
 +
 +
=== Макроструктура алгоритма ===
 +
 +
=== Схема реализации последовательного алгоритма ===
 +
 +
=== Последовательная сложность алгоритма ===
 +
 +
=== Информационный граф ===
 +
 +
=== Ресурс параллелизма алгоритма ===
 +
 +
=== Входные и выходные данные алгоритма ===
 +
 +
=== Свойства алгоритма ===
 +
 +
== Программная реализация алгоритма ==
 +
=== Особенности реализации последовательного алгоритма ===
 +
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
=== Результаты прогонов ===
 +
=== Выводы для классов архитектур ===
 +
== Литература ==
 +
 +
<references />
 +
 +
[[Категория:Статьи в работе]]
 +
[[Категория:Алгоритмы для матриц]]
 +
[[Категория:ScaLaPACK]]
 +
 +
[[en:QR algorithm as implemented in SCALAPACK]]

Текущая версия на 09:13, 9 июля 2022


QR-алгоритм, используемый в SCALAPACK - алгоритм, использующий все проверенные на настоящее время приёмы ускорения QR-алгоритма для несимметричных вещественных матриц[1]. Включён своими частями в разные подпрограммы пакета SCALAPACK[2]. Состоит из двух основных частей: ортогонально подобного приведения матрицы к хессенберговому виду и QR-итераций со сдвигами для хессенберговой матрицы.

Содержание

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

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

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

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

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

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

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

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

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

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

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

[math]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].


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

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

1.2.2.1 Использование хессенберговой формы матрицы

Хессенбергова (почти треугольная) форма матрицы интересна в приложении к QR-алгоритму тем, что в случае, если [math]A_0 = A[/math] имеет эту форму, то имеют её и все матрицы [math]A_k[/math]. Если теперь посмотреть на методы QR-разложения плотных неособенных матриц, то видно, что выполнение разложения [math]A_k = Q_kR_k[/math] в случае плотной неособенной [math]A_k[/math]требует (в последовательном варианте) [math]O(n^3)[/math] операций, как и последующий процесс вычисления [math]A_{k+1} = R_kQ_k[/math]. В случае же, если матрицы хессенберговы, то как QR-разложения плотных хессенберговых матриц, так и вычисления [math]A_{k+1} = R_kQ_k[/math] потребуют уже по [math]O(n^2)[/math] операций. Это позволяет экономить довольно большие вычислительные ресурсы.

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

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

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

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

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

При этом сохраняется свойство подобия матриц [math]A_k[/math] и [math]A_{k+1}[/math]:

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

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

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

1.2.2.4 Обнуление поддиагональных элементов (дефляция) и распараллеливание по блокам

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

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

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

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

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

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

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

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

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

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

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

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

2.3 Результаты прогонов

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

3 Литература

  1. ↑ Алгоритм для неэрмитовых комплексных матриц в SCALAPACK отсутствует, есть только их приведение к хессенберговой форме
  2. ↑ 2,0 2,1 http://www.netlib.org/scalapack/slug/node63.html#SECTION04335100000000000000
  3. ↑ 3,0 3,1 3,2 3,3 Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.
  4. ↑ 4,0 4,1 4,2 Деммель Д. Вычислительная линейная алгебра. – 2001. - С.261-264.
  5. ↑ R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems".
  6. ↑ Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.
Источник — https://algowiki-project.org/w/ru/index.php?title=QR-алгоритм,_используемый_в_SCALAPACK&oldid=29098