Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией: различия между версиями
Sfedotov (обсуждение | вклад) (Новая страница: «= Свойства и структура алгоритмов = == Общее описание алгоритма == '''Алгоритм Ланцоша с по…») |
Sfedotov (обсуждение | вклад) |
||
| Строка 1: | Строка 1: | ||
| + | {{algorithm | ||
| + | | name = Алгоритм Ланцоша с полной переортогонализацией | ||
| + | | serial_complexity = <math>O(n^2)</math> | ||
| + | | pf_height = <math>O(m)</math> | ||
| + | | pf_width = <math>O(n^2)</math> | ||
| + | | input_data = <math>\frac{n(n + 1)}{2}</math> | ||
| + | | output_data = <math>m(n + 1)</math> | ||
| + | }} | ||
| + | |||
| + | Автор описания: [[Участник:Sfedotov|Федотов С.А.]] | ||
| + | |||
| + | ---- | ||
| + | |||
= Свойства и структура алгоритмов = | = Свойства и структура алгоритмов = | ||
== Общее описание алгоритма == | == Общее описание алгоритма == | ||
| − | '''Алгоритм Ланцоша | + | '''Алгоритм Ланцоша''' --- итерационный метод вычисления собственных значений и собственных векторов веществнной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности -[[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов. |
| − | + | ||
== Математическое описание алгоритма == | == Математическое описание алгоритма == | ||
| + | Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=A<sup>T<sup>''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''T<sub>k</sub>=Q<sub>k</sub><sup>T</sup>AQ<sub>k</sub>'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Q<sub>k</sub> = [q<sub>1</sub>,q<sub>2</sub>,...,q<sub>k</sub>,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''T<sub>k</sub>=V<math>\Lambda</math>V<sup>T</sup>'' есть спектральное разложение матрицы ''T<sub>k</sub>'', столбцы матрицы ''Q<sub>k</sub>V'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''><ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>: | ||
| + | <math> | ||
| + | \begin{align} | ||
| + | q_1 = & \frac{b}{\|b\|_2}\\ | ||
| + | for \; & j = 1 \; to \; k \\ | ||
| + | & z = Aq_j\\ | ||
| + | & \alpha_j = q^T_j z\\ | ||
| + | & z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\ | ||
| + | & \beta_j = \|z\|_2\\ | ||
| + | & q_{j+1} = \frac{z}{\beta_j}\\ | ||
| + | end \; & for | ||
| + | \end{align} | ||
| + | </math> | ||
| + | |||
| + | Диагональные элементы обозначены как <math>\alpha=t_{jj}</math>, а элементы побочной диагонали - <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица | ||
| + | <math>T_{jj} = \begin{pmatrix} | ||
| + | \alpha_1 & \beta_2 & & & & 0 \\ | ||
| + | \beta_2 & \alpha_2 & \beta_3 & & & \\ | ||
| + | & \beta_3 & \alpha_3 & \ddots & & \\ | ||
| + | & & \ddots & \ddots & \beta_{j-1} & \\ | ||
| + | & & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\ | ||
| + | 0 & & & & \beta_j & \alpha_j \\ | ||
| + | \end{pmatrix}</math> | ||
| + | |||
| + | Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q<sub>1</sub>,q<sub>2</sub>,...,q<sub>j-1</sub>''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''T<sub>k<sub>'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов такой матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''<sup>2</sup>) флопов. | ||
== Вычислительное ядро алгоритма == | == Вычислительное ядро алгоритма == | ||
| + | За один шаг алгоритма можно выделить следующие вычислительные ядра: | ||
| + | *умножение исходной матрицы на вектор <math>z=Aq_j</math> | ||
| + | *процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math> | ||
== Макроструктура алгоритма == | == Макроструктура алгоритма == | ||
| − | + | Основные состовляющие алгоритма это умножение матрицы на вектор, процесс переортогонализации и поиск всех собственных значений и всех собственных векторов симметричной трехдиагональной матрицы | |
== Схема реализации последовательного алгоритма == | == Схема реализации последовательного алгоритма == | ||
| + | [[file:Lanczos1.jpg|thumb|center|450px|Блок-схема алгоритма]] | ||
== Последовательная сложность алгоритма == | == Последовательная сложность алгоритма == | ||
| − | + | Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализации требуется выполнить: | |
| + | *''n+2k'' умножений вектора на вектор | ||
| + | *''k+1'' умножение вектора на константу | ||
| + | *''2'' сложения векторов | ||
| + | *''O''(''k''<sup>2</sup>) операций на разложение матрицы ''T<sub>k<sub>'' | ||
| + | Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''<sup>2</sup>''n'') флопов и ''O''(''mn'') слов памяти | ||
== Информационный граф == | == Информационный граф == | ||
| − | + | [[file:Lanczos2.jpg|thumb|center|300px|Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]] | |
== Ресурс параллелизма алгоритма == | == Ресурс параллелизма алгоритма == | ||
| + | Алгоритм Ланцоша с полной переортогонализацией в целом последовательный. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q<sub>j+1</sub></math>, которые служат входными данными для ''j+1''-го шага. Нетрудно заметить, что процесс переортогонализации ''j''-го приближения может быть вычислен параллельно | ||
| + | == Входные и выходные данные алгоритма == | ||
| + | '''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>. | ||
| + | Дополнительные ограничения: | ||
| + | * <math>A</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>. | ||
| − | + | '''Объём входных данных''': <math>\frac{n (n + 1)}{2}</math> | |
| + | '''Выходные данные''': | ||
| + | *вектор приближенных собственных значений <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>). | ||
| + | *приближения к собственным векторам ''<math>q_i,i=1, m </math> | ||
| + | '''Объём выходных данных''': <math> m(n+1) </math>. | ||
== Свойства алгоритма == | == Свойства алгоритма == | ||
Версия 20:11, 15 октября 2016
| Алгоритм Ланцоша с полной переортогонализацией | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]O(n^2)[/math] |
| Объём входных данных | [math]\frac{n(n + 1)}{2}[/math] |
| Объём выходных данных | [math]m(n + 1)[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]O(m)[/math] |
| Ширина ярусно-параллельной формы | [math]O(n^2)[/math] |
Автор описания: Федотов С.А.
Содержание
- 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 Общее описание алгоритма
Алгоритм Ланцоша --- итерационный метод вычисления собственных значений и собственных векторов веществнной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный Корнелием Ланцошом, адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности n с ограниченным числом итераций, m, где m много меньше n. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman[1] показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций m. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности -Алгоритм Ланцоша с выборочной ортогонализацией. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.
1.2 Математическое описание алгоритма
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица A=ATи вектор начального приближения b. Мы пытаемся найти трехдиагональную симметричную матрицу Tk=QkTAQk, собственные значения которой приближают собственные значения матрицы A. Иными словами на k-м шаге из ортонормированных векторов Ланцоша строится матрица Qk = [q1,q2,...,qk,] и в качестве приближенных собственных значений матрицы A принимаются числа Ритца. Пусть Tk=V[math]\Lambda[/math]VT есть спектральное разложение матрицы Tk, столбцы матрицы QkV рассматриваются как приближения к соответсвующим собственным векторам матрицы A>[2]: [math] \begin{align} q_1 = & \frac{b}{\|b\|_2}\\ for \; & j = 1 \; to \; k \\ & z = Aq_j\\ & \alpha_j = q^T_j z\\ & z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\ & \beta_j = \|z\|_2\\ & q_{j+1} = \frac{z}{\beta_j}\\ end \; & for \end{align} [/math]
Диагональные элементы обозначены как [math]\alpha=t_{jj}[/math], а элементы побочной диагонали - [math]\beta_j=t_{j-1,j}=t_{j,j-1}[/math]. После каждой итерации мы вычисляем [math]\alpha_j, \beta_j[/math], из которых строится матрица [math]T_{jj} = \begin{pmatrix} \alpha_1 & \beta_2 & & & & 0 \\ \beta_2 & \alpha_2 & \beta_3 & & & \\ & \beta_3 & \alpha_3 & \ddots & & \\ & & \ddots & \ddots & \beta_{j-1} & \\ & & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\ 0 & & & & \beta_j & \alpha_j \\ \end{pmatrix}[/math]
Полная переортогонализация соответствует повторному проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что [math]z[/math] будет ортогонален векторам q1,q2,...,qj-1. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем повторные копии сошедшихся чисел Ритца. То есть для больших k матрица Tk может иметь не одно собственное значение, очень близкое к [math]\lambda_i(A)[/math], но много таких собственных значений. Нахождение собственных значений и собственных векторов такой матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи MRRR алгоритма они могут быть вычислены за O(k2) флопов.
1.3 Вычислительное ядро алгоритма
За один шаг алгоритма можно выделить следующие вычислительные ядра:
- умножение исходной матрицы на вектор [math]z=Aq_j[/math]
- процесс ортогонализации Грама-Шмидта [math]z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i[/math]
1.4 Макроструктура алгоритма
Основные состовляющие алгоритма это умножение матрицы на вектор, процесс переортогонализации и поиск всех собственных значений и всех собственных векторов симметричной трехдиагональной матрицы
1.5 Схема реализации последовательного алгоритма
1.6 Последовательная сложность алгоритма
Для вычисления k-й итерации алгоритма Ланцоша с полной переортогонализации требуется выполнить:
- n+2k умножений вектора на вектор
- k+1 умножение вектора на константу
- 2 сложения векторов
- O(k2) операций на разложение матрицы Tk
Легко проверить, что m шагов этого алгоритма потребуют O(m2n) флопов и O(mn) слов памяти
1.7 Информационный граф
1.8 Ресурс параллелизма алгоритма
Алгоритм Ланцоша с полной переортогонализацией в целом последовательный. На j-м шаге вычисляются числа [math]\alpha_j, \beta_j[/math] и ортонормированный вектор [math]q\lt sub\gt j+1\lt /sub\gt [/math], которые служат входными данными для j+1-го шага. Нетрудно заметить, что процесс переортогонализации j-го приближения может быть вычислен параллельно
1.9 Входные и выходные данные алгоритма
Входные данные: матрица [math]A\in\mathbb{R}^{n\times n}[/math]. Дополнительные ограничения:
- [math]A[/math] – симметрическая матрица, т. е. [math]a_{ij}= a_{ji}, i, j = 1, \ldots, n[/math].
Объём входных данных: [math]\frac{n (n + 1)}{2}[/math]
Выходные данные:
- вектор приближенных собственных значений [math]\Lambda[/math] (элементы [math]\lambda_{ii}[/math]).
- приближения к собственным векторам [math]q_i,i=1, m [/math]
Объём выходных данных: [math] m(n+1) [/math].
