Участник:Kozlov Vladimir/Алгоритм Ланцоша для арифметики с плавающей точкой с полной переортогонализацией: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
Строка 15: Строка 15:
 
\end{array} \right].</math>
 
\end{array} \right].</math>
  
Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math> — приближённое спектральное разложение матрицы <math>A</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>. Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно.
+
Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.  
  
В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама-Шмидта.
+
Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
 +
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
 +
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.
 +
 
 +
В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама Шмидта.
  
 
=== Математическое описание алгоритма ===
 
=== Математическое описание алгоритма ===

Версия 18:30, 15 октября 2016

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

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

Алгоритм Ланцоша — это итерационный алгоритм поиска k приближённых собственных значений симметричной вещественной матрицы A размера n \times n. Алгоритм применяется, когда матрица A слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм метод Рэлея — Ритца поиска приближённых собственных значений и метод Ланцоша построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска k приближённых собственных значений симметричной вещественной матрицы A размера n \times n. Если Q = [Q_k, Q_u] — ортонормированная матрица размера n \times n, Q_k имеет размер n \times k, Q_u имеет размер n \times n - k, то можно записать равенство

T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] = \left[ \begin{array}{cc} Q_k^T A Q_k & Q_k^T A Q_u\\ Q_u^T A Q_k & Q_u^T A Q_u \end{array} \right] = \left[ \begin{array}{cc} T_{k} & T_{ku}^T\\ T_{ku} & T_{u} \end{array} \right].

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы T_k = Q_k^T A Q_k объявляются приближёнными собственными значениями матрицы A. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если T_k = V \Lambda V^{-1} — спектральное разложение T_k, то пара (Q_k V, \Lambda) минимизирует функционал L(P_k, D) = \Vert A P_k - P_k D \Vert_2, причём L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2, то есть A \approx (Q_k V) \Lambda (Q_k V)^{-1}. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы A. Более того, можно показать, что собственные значения T отличаются от некоторых собственных значений A не более чем на \Vert T_{ku} \Vert_2.

Метод Ланцоша — это метод построения матрицы Q, при использовании которого, во-первых, матрица T оказывается симметричной трёхдиагональной, во-вторых, столбцы Q и T вычисляются последовательно. Трёхдиагональность T приводит к следующим явлениям:

  1. матрица T_k является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
  2. матрица T_{ku} имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца q_{j + 1} матрицы Q достаточно знать только q_{j - 1} и q_{j} в силу трёхдиагональности матрицы T. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов q_{1}, \dots, q_{k} перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

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

Исходные данные: симметрическая матрица A — матрица, для которой будут вычисляться собственные значения, вектор b — начальное приближение для метода Ланцоша.

Формулы метода:

\begin{array}{l} q_1 = b / \Vert b \Vert_2, \beta_0 = 0, q_0 = 0\\ j = \overline{1, k}:\\ \quad z_j = A q_j \\ \quad \alpha_j = q_j^T A q_j = q_j^T z_j \\ \quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\ \quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\ \quad \beta_j = \Vert z_j'' \Vert_2\\ \quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j \end{array}

Вычисление z_j', z_j'' — это полная переортогонализация z_j методом Грама — Шмидта. Двойной запуск практически гарантирует, что z_j'' будет ортогонален q_1, \dots, q_j.

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

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

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

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

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

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

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

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

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

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

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

3 Литература