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

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

Материал из Алговики
Перейти к навигации Перейти к поиску
[досмотренная версия][досмотренная версия]
м (Frolov переименовал страницу Наиболее употребимый QR-алгоритм в QR-алгоритм, используемый в SCALAPACK: уточнение по смыслу)
Строка 1: Строка 1:
 
{{level-a}}
 
{{level-a}}
 +
 +
'''QR-алгоритм, используемый в SCALAPACK''' - алгоритм, использующий все проверенные на настоящее время приёмы ускорения [[QR-алгоритм|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> и поиска собственных значений для нее.
 +
 +
== Математическое описание ==
 +
 +
Известно<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_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-разложения плотных хессенберговых матриц|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"> Деммель Д. Вычислительная линейная алгебра. – 2001. - С.261-264. </ref>, в качестве сдвигов выбираются для каждого из диагональных блоков собственные значения их нижних правых клеток размером 2х2, и среднее количество итераций для большинства матриц получается порядка 2 на одно собственное значение. Существуют, однако, матрицы, для которых правило Фрэнсиса не даёт сходимости, поэтому если по истечении 10 итераций с обычными сдвигами не получилось сходимости, применяют так называемые "особые сдвиги".
 +
 +
===== Особенности вещественного варианта QR-алгоритма =====
 +
 +
Если вещественная матрица <math>A</math> имеет различные вещественные собственные значения, то QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения, и тогда алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка. Блоки 1-го порядка будут содержать различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений<ref name = "Dem" /><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>.
 +
 +
Поэтому выбор сдвигов в алгоритме учитывает, что конечные матрицы будут представленного вида. Чтобы не выполнять итерации в комплексной арифметике, выполняют т.н. неявный двойной сдвиг, давно описанный в литературе<ref name="VOLA" /><ref name = "Dem" />.
 +
 +
==== Обнуление поддиагональных элементов ====
 +
 +
После того, как под диагональю какой-либо элемент становится менее порогового значения, его можно считать равным нулю. Это даёт возможность применить разделение диагональных блоков матрицы, соседствующих с ним, и раздельное вычисление их собственных чисел.

Версия 14:45, 4 августа 2017


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

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] и поиска собственных значений для нее.

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] сходимости получающихся матриц к клеточной правой треугольной матрице с диагональными клетками, размеры которых зависят от размеров канонических ящиков Жордана исходной матрицы. Таким образом, проблема собственных значений матрицы сводится к проблемам собственных значений матриц меньшего размера.

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

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

2.1.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] операций. Это позволяет экономить довольно большие вычислительные ресурсы.

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

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

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

QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости[3]. Пусть у нас есть матрица [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 итераций с обычными сдвигами не получилось сходимости, применяют так называемые "особые сдвиги".

2.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].

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

2.1.3 Обнуление поддиагональных элементов

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

  1. 1,0 1,1 http://www.netlib.org/scalapack/slug/node63.html#SECTION04335100000000000000
  2. 2,0 2,1 2,2 2,3 Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.
  3. Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.
  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".