QR-алгоритм для хессенберговой матрицы, используемый в SCALAPACK: различия между версиями
[досмотренная версия] | [досмотренная версия] |
Frolov (обсуждение | вклад) |
ASA (обсуждение | вклад) |
||
(не показана 1 промежуточная версия этого же участника) | |||
Строка 3: | Строка 3: | ||
'''QR-алгоритм для хессенберговой матрицы, используемый в SCALAPACK''' - вторая часть [[QR-алгоритм, используемый в SCALAPACK|QR-алгоритма для квадратных несимметричных матриц, используемого в SCALAPACK]]<ref name="SCALAEig">http://www.netlib.org/scalapack/slug/node63.html#SECTION04335100000000000000</ref>. | '''QR-алгоритм для хессенберговой матрицы, используемый в SCALAPACK''' - вторая часть [[QR-алгоритм, используемый в SCALAPACK|QR-алгоритма для квадратных несимметричных матриц, используемого в SCALAPACK]]<ref name="SCALAEig">http://www.netlib.org/scalapack/slug/node63.html#SECTION04335100000000000000</ref>. | ||
− | = Свойства и структура алгоритма = | + | == Свойства и структура алгоритма == |
− | == Общее описание алгоритма == | + | === Общее описание алгоритма === |
Суть [[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> и поиска собственных значений для нее. | Суть [[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> и поиска собственных значений для нее. | ||
Строка 11: | Строка 11: | ||
В описываемом алгоритме исходная матрица ''хессенбергова''. Благодаря инвариантности хессенберговой формы к QR-алгоритму видно, что так как <math>A_0 = A</math> имеет эту форму, то имеют её и все матрицы <math>A_k</math>. | В описываемом алгоритме исходная матрица ''хессенбергова''. Благодаря инвариантности хессенберговой формы к QR-алгоритму видно, что так как <math>A_0 = A</math> имеет эту форму, то имеют её и все матрицы <math>A_k</math>. | ||
− | == Математическое описание == | + | === Математическое описание === |
Известно<ref name="VOLA">Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.</ref>, что произвольная квадратная матрица может быть представлена в виде произведения унитарной (''в вещественном случае ортогональной'') и верхней треугольной матриц. Такое разложение называется [[QR-разложения плотных неособенных матриц|QR-разложением]]. | Известно<ref name="VOLA">Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.</ref>, что произвольная квадратная матрица может быть представлена в виде произведения унитарной (''в вещественном случае ортогональной'') и верхней треугольной матриц. Такое разложение называется [[QR-разложения плотных неособенных матриц|QR-разложением]]. | ||
Строка 24: | Строка 24: | ||
Существуют доказательства<ref name="VOLA" /> сходимости получающихся матриц к клеточной правой треугольной матрице с диагональными клетками, размеры которых зависят от размеров канонических ящиков Жордана исходной матрицы. Таким образом, проблема собственных значений матрицы сводится к проблемам собственных значений матриц меньшего размера. | Существуют доказательства<ref name="VOLA" /> сходимости получающихся матриц к клеточной правой треугольной матрице с диагональными клетками, размеры которых зависят от размеров канонических ящиков Жордана исходной матрицы. Таким образом, проблема собственных значений матрицы сводится к проблемам собственных значений матриц меньшего размера. | ||
− | === Особенности вещественного варианта QR-алгоритма === | + | ==== Особенности вещественного варианта QR-алгоритма ==== |
Если вещественная матрица <math>A</math> имеет различные вещественные собственные значения, то QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения, и тогда алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка<ref name = "Dem"> Деммель Д. Вычислительная линейная алгебра. – 2001. - С.261-264. </ref>. Блоки 1-го порядка будут содержать различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений<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</math> имеет различные вещественные собственные значения, то QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения, и тогда алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка<ref name = "Dem"> Деммель Д. Вычислительная линейная алгебра. – 2001. - С.261-264. </ref>. Блоки 1-го порядка будут содержать различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений<ref>R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems".</ref>. | ||
Строка 41: | Строка 41: | ||
\end{bmatrix}</math>. | \end{bmatrix}</math>. | ||
− | === Приёмы ускорения сходимости, используемые в алгоритме === | + | ==== Приёмы ускорения сходимости, используемые в алгоритме ==== |
− | ==== Одновременная обработка диагональных блоков и дефляция ==== | + | ===== Одновременная обработка диагональных блоков и дефляция ===== |
Известно, что в случае, если хессенбергова матрица содержит на нижней кодиагонали только ненулевые элементы, то все кратные собственные значения будут соответствовать одному ящику Жордана<ref name="VOLA" />. Это означает, что каждое собственное значение с геометрической кратностью выше 1 даст после точного приведения матрицы к хессенбергову виду нули на нижней кодиагонали, что сразу даёт возможность параллельно работать с несколькими диагональными блоками. Эта возможность используется в описываемом алгоритме. | Известно, что в случае, если хессенбергова матрица содержит на нижней кодиагонали только ненулевые элементы, то все кратные собственные значения будут соответствовать одному ящику Жордана<ref name="VOLA" />. Это означает, что каждое собственное значение с геометрической кратностью выше 1 даст после точного приведения матрицы к хессенбергову виду нули на нижней кодиагонали, что сразу даёт возможность параллельно работать с несколькими диагональными блоками. Эта возможность используется в описываемом алгоритме. | ||
Строка 49: | Строка 49: | ||
После того, как в процессе итераций под диагональю какой-либо элемент становится менее порогового значения, его можно считать равным нулю, в чём и заключается суть приёма, называемого '''дефляцией'''. Это даёт возможность применить дальнейшее разделение диагональных блоков матрицы, соседствующих с ним, и раздельное вычисление их собственных чисел. | После того, как в процессе итераций под диагональю какой-либо элемент становится менее порогового значения, его можно считать равным нулю, в чём и заключается суть приёма, называемого '''дефляцией'''. Это даёт возможность применить дальнейшее разделение диагональных блоков матрицы, соседствующих с ним, и раздельное вычисление их собственных чисел. | ||
− | ==== Сдвиги QR-алгоритма ==== | + | ===== Сдвиги QR-алгоритма ===== |
QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости<ref name="bahvalov">Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.</ref>. Пусть у нас есть матрица <math>A_k</math>, тогда процесс перехода к матрице <math>A_{k+1}</math> выглядит следующим образом: | QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости<ref name="bahvalov">Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.</ref>. Пусть у нас есть матрица <math>A_k</math>, тогда процесс перехода к матрице <math>A_{k+1}</math> выглядит следующим образом: | ||
Строка 66: | Строка 66: | ||
[[File:Выступn-1.png|thumb|left|250px|Рисунок 3. Область действия последнего шага (вытеснение <math>n-2</math>го выступа) итерации, состоящей из левого преобразования Хаусхолдера, меняющего 2 строки, и правого преобразования Хаусхолдера, меняющего 2 столбца матрицы.]] | [[File:Выступn-1.png|thumb|left|250px|Рисунок 3. Область действия последнего шага (вытеснение <math>n-2</math>го выступа) итерации, состоящей из левого преобразования Хаусхолдера, меняющего 2 строки, и правого преобразования Хаусхолдера, меняющего 2 столбца матрицы.]] | ||
− | ==== Неявное выполнение сдвигов ==== | + | ===== Неявное выполнение сдвигов ===== |
Чтобы не выполнять итерации в комплексной арифметике, выполняют т.н. неявный двойной сдвиг, давно описанный в литературе<ref name="VOLA" /><ref name = "Dem" />. При этом используются двусторонние подобные операции Хаусхолдера (отражения) 3го порядка. На первом шаге итерации формируют у хессенберговой формы выступ размера 2х2, а на последующих шагах с их же помощью выступ вытесняют вправо вниз. На последнем шаге выступ размера 1х1. | Чтобы не выполнять итерации в комплексной арифметике, выполняют т.н. неявный двойной сдвиг, давно описанный в литературе<ref name="VOLA" /><ref name = "Dem" />. При этом используются двусторонние подобные операции Хаусхолдера (отражения) 3го порядка. На первом шаге итерации формируют у хессенберговой формы выступ размера 2х2, а на последующих шагах с их же помощью выступ вытесняют вправо вниз. На последнем шаге выступ размера 1х1. | ||
− | == Вычислительное ядро алгоритма == | + | === Вычислительное ядро алгоритма === |
Основная часть алгоритма состоит из выполняемых над диагональными блоками, оставшимися с размерами более, чем 2, очередной итерации. При этом очередная итерация состоит из вычисления и применения к блоку двусторонних подобных преобразований Хаусхолдера (отражения) размерности 3 (последнее - размерности 2), выполняемыми для стартового создания выступа 2х2 и затем для его вытеснения. | Основная часть алгоритма состоит из выполняемых над диагональными блоками, оставшимися с размерами более, чем 2, очередной итерации. При этом очередная итерация состоит из вычисления и применения к блоку двусторонних подобных преобразований Хаусхолдера (отражения) размерности 3 (последнее - размерности 2), выполняемыми для стартового создания выступа 2х2 и затем для его вытеснения. | ||
− | == Макроструктура алгоритма == | + | === Макроструктура алгоритма === |
На макроуровне алгоритм состоит из последовательности отдельных итераций. На каждой итерации одна и та же пошаговая процедура формирования и затем вытеснения выступа проводится независимо над разными диагональными блоками, которые разделены в "текущей" матрице нулями на кодиагонали. При этом непосредственно при вытеснении выступа в каждом блоке очередные новые кодиагональные элементы над вытесненным выступом могут быть сразу проверены на их малость, и при выполнении критерия малости к ним применяется процедура дефляции. | На макроуровне алгоритм состоит из последовательности отдельных итераций. На каждой итерации одна и та же пошаговая процедура формирования и затем вытеснения выступа проводится независимо над разными диагональными блоками, которые разделены в "текущей" матрице нулями на кодиагонали. При этом непосредственно при вытеснении выступа в каждом блоке очередные новые кодиагональные элементы над вытесненным выступом могут быть сразу проверены на их малость, и при выполнении критерия малости к ним применяется процедура дефляции. | ||
− | == Схема реализации последовательного алгоритма == | + | === Схема реализации последовательного алгоритма === |
− | == Последовательная сложность алгоритма == | + | === Последовательная сложность алгоритма === |
− | == Информационный граф == | + | === Информационный граф === |
− | == Ресурс параллелизма алгоритма == | + | === Ресурс параллелизма алгоритма === |
− | == Входные и выходные данные алгоритма == | + | === Входные и выходные данные алгоритма === |
− | == Свойства алгоритма == | + | === Свойства алгоритма === |
− | = Программная реализация алгоритма = | + | == Программная реализация алгоритма == |
+ | === Особенности реализации последовательного алгоритма === | ||
+ | === Возможные способы и особенности параллельной реализации алгоритма === | ||
+ | === Результаты прогонов === | ||
+ | === Выводы для классов архитектур === | ||
+ | == Литература == | ||
− | |||
<references /> | <references /> | ||
+ | |||
+ | [[Категория:Статьи в работе]] | ||
+ | |||
+ | [[en:Hessenberg QR algorithm as implemented in SCALAPACK]] |
Текущая версия на 09:16, 9 июля 2022
QR-алгоритм для хессенберговой матрицы, используемый в SCALAPACK - вторая часть QR-алгоритма для квадратных несимметричных матриц, используемого в SCALAPACK[1].
Содержание
- 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 Общее описание алгоритма
Суть 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] и поиска собственных значений для нее.
В описываемом алгоритме исходная матрица хессенбергова. Благодаря инвариантности хессенберговой формы к QR-алгоритму видно, что так как [math]A_0 = A[/math] имеет эту форму, то имеют её и все матрицы [math]A_k[/math].
1.2 Математическое описание
Известно[2], что произвольная квадратная матрица может быть представлена в виде произведения унитарной (в вещественном случае ортогональной) и верхней треугольной матриц. Такое разложение называется 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] и имеют те же собственные значения.
Существуют доказательства[2] сходимости получающихся матриц к клеточной правой треугольной матрице с диагональными клетками, размеры которых зависят от размеров канонических ящиков Жордана исходной матрицы. Таким образом, проблема собственных значений матрицы сводится к проблемам собственных значений матриц меньшего размера.
1.2.1 Особенности вещественного варианта QR-алгоритма
Если вещественная матрица [math]A[/math] имеет различные вещественные собственные значения, то QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения, и тогда алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка[3]. Блоки 1-го порядка будут содержать различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений[4].
[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 Одновременная обработка диагональных блоков и дефляция
Известно, что в случае, если хессенбергова матрица содержит на нижней кодиагонали только ненулевые элементы, то все кратные собственные значения будут соответствовать одному ящику Жордана[2]. Это означает, что каждое собственное значение с геометрической кратностью выше 1 даст после точного приведения матрицы к хессенбергову виду нули на нижней кодиагонали, что сразу даёт возможность параллельно работать с несколькими диагональными блоками. Эта возможность используется в описываемом алгоритме.
После того, как в процессе итераций под диагональю какой-либо элемент становится менее порогового значения, его можно считать равным нулю, в чём и заключается суть приёма, называемого дефляцией. Это даёт возможность применить дальнейшее разделение диагональных блоков матрицы, соседствующих с ним, и раздельное вычисление их собственных чисел.
1.2.2.2 Сдвиги QR-алгоритма
QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости[5]. Пусть у нас есть матрица [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 =[/math] [math]= 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]. При алгоритме подбора параметров, известном, как правило Фрэнсиса[3], в качестве сдвигов выбираются для каждого из диагональных блоков собственные значения их нижних правых клеток размером 2х2, и среднее количество итераций для большинства матриц получается порядка 2 на одно собственное значение. Существуют, однако, матрицы, для которых правило Фрэнсиса не даёт сходимости, поэтому если по истечении 10 итераций с обычными сдвигами не получилось сходимости, для такого диагонального блока применяют так называемые "особые сдвиги".
1.2.2.3 Неявное выполнение сдвигов
Чтобы не выполнять итерации в комплексной арифметике, выполняют т.н. неявный двойной сдвиг, давно описанный в литературе[2][3]. При этом используются двусторонние подобные операции Хаусхолдера (отражения) 3го порядка. На первом шаге итерации формируют у хессенберговой формы выступ размера 2х2, а на последующих шагах с их же помощью выступ вытесняют вправо вниз. На последнем шаге выступ размера 1х1.
1.3 Вычислительное ядро алгоритма
Основная часть алгоритма состоит из выполняемых над диагональными блоками, оставшимися с размерами более, чем 2, очередной итерации. При этом очередная итерация состоит из вычисления и применения к блоку двусторонних подобных преобразований Хаусхолдера (отражения) размерности 3 (последнее - размерности 2), выполняемыми для стартового создания выступа 2х2 и затем для его вытеснения.
1.4 Макроструктура алгоритма
На макроуровне алгоритм состоит из последовательности отдельных итераций. На каждой итерации одна и та же пошаговая процедура формирования и затем вытеснения выступа проводится независимо над разными диагональными блоками, которые разделены в "текущей" матрице нулями на кодиагонали. При этом непосредственно при вытеснении выступа в каждом блоке очередные новые кодиагональные элементы над вытесненным выступом могут быть сразу проверены на их малость, и при выполнении критерия малости к ним применяется процедура дефляции.
1.5 Схема реализации последовательного алгоритма
1.6 Последовательная сложность алгоритма
1.7 Информационный граф
1.8 Ресурс параллелизма алгоритма
1.9 Входные и выходные данные алгоритма
1.10 Свойства алгоритма
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Возможные способы и особенности параллельной реализации алгоритма
2.3 Результаты прогонов
2.4 Выводы для классов архитектур
3 Литература
- ↑ http://www.netlib.org/scalapack/slug/node63.html#SECTION04335100000000000000
- ↑ 2,0 2,1 2,2 2,3 Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.
- ↑ 3,0 3,1 3,2 Деммель Д. Вычислительная линейная алгебра. – 2001. - С.261-264.
- ↑ R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems".
- ↑ Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.