Учacтник:Malikovmt/Алгоритм Ланцоша для арифметики с плавающей точкой с полной переортогонализацией: различия между версиями
Malikovmt (обсуждение | вклад) |
Malikovmt (обсуждение | вклад) |
||
Строка 36: | Строка 36: | ||
\quad z_j = A q_j \\ | \quad z_j = A q_j \\ | ||
\quad \alpha_j = q_j^T A q_j = q_j^T z_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 | + | \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 \beta_j = \Vert z_j \Vert_2\\ | ||
− | \quad q_{j+1} = z_j / \Vert z_j | + | \quad q_{j+1} = z_j / \Vert z_j \Vert_2 = z_j/\beta_j |
\end{array} | \end{array} | ||
</math> | </math> |
Версия 20:31, 1 февраля 2017
Алгоритм Ланцоша с полной переортогонализацией | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(n^2k)[/math] |
Объём входных данных | [math]\frac{n(n + 1)}{2}[/math] |
Объём выходных данных | [math]k(n + 1)[/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math]O(k)[/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 Свойства алгоритма
- 1.11 Локальность данных и вычислений
- 1.12 Возможные способы и особенности параллельной реализации алгоритма
- 1.13 Масштабируемость алгоритма и его реализации
- 1.14 Динамические характеристики и эффективность реализации алгоритма
- 1.15 Выводы для классов архитектур
- 1.16 Существующие реализации алгоритма
- 2 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Алгоритм Ланцоша - итерационный метод, используемый для вычисления части собственных значений и соответствующих им собственных векторов матрицы [math]A[/math] размера [math]n*n[/math], изначально разработанный Корнелием Ланцошем. Преимуществами использования метода является относительно небольшое потребление памяти и вычислительных ресурсов, а также наличие параметра [math]k \lt \lt n[/math], контролирующего количество итераций. Несмотря на то, что алгоритм является вычислительно эффективным, первоначально сформулированный метод был плохо применим из-за численной неустойчивости - метод хорошо работал на целочисленных значениях, однако в арифметике с плавающей точкой ошибки округления давали большую погрешность. В 1970 году Ojalvo и Newman показали, как сделать метод численно стабильным и применили его для расчета крупных инженерных сооружений, подверженных динамическим нагрузкам. Кроме того, они показали способ выбора начального приближения (с использованием ГПСЧ), а также эмпирический способ для выбора числа [math]k[/math] (примерно в полтора раза больше искомого числа собственных векторов). В данный момент существует две основных модификации метода (с полной и выборочной переортогонализацией), а также большое количество модификаций, использующихся в различных технических областях. Алгоритм используется для больших [math]n[/math].
1.2 Математическое описание алгоритма
Первый этап алгоритма - использование метода Ланцоша для построения крыловского подпространства: [math] K_k(A,x) = span[x_1, Ax_1, A^2x_1, ..., A^{k-1}x_1] [/math]. Входные данные алгоритма: квадратная симметричная матрица [math]A[/math] размерности [math]n*n[/math], вектор начального приближения [math]b[/math], а так же число итераций [math]k[/math]. Метод осуществляет поиск трехдиагональной симметричной матрицы [math]T_k=Q_k^TAQ_k[/math].
[math]T_k=\begin{bmatrix} \alpha_1 & \beta_2 \\ \beta_2 & \alpha_2 & \beta_3 &\\ &. & . & .\\ &&\beta_{k-1} & \alpha_{k-1} & \beta_k\\ &&&\beta_k & \alpha_k \end{bmatrix}[/math]
Описание метода:
[math] \begin{array}{l} q_1 = b / \Vert b \Vert_2\\ 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 \beta_j = \Vert z_j \Vert_2\\ \quad q_{j+1} = z_j / \Vert z_j \Vert_2 = z_j/\beta_j \end{array} [/math]
Следующий шаг алгоритма - процедура Рэлея-Ритца. Она зкалючается в интерпретации собственных значений матрицы [math] T_k=Q_k^TAQ_k[/math]. Ее собственные значения приближают собственные значения исходной матрицы. Пусть Tk=V[math]\Lambda[/math]VT - спектральное разложение матрицы Tk, тогда столбцы матрицы QkV рассматриваются как приближения к соответствующим собственным векторам матрицы A и называются векторами Ритца. Числа и векторы Ритца являются оптимальными приближениями к собственным значениям и собственным векторам матрицы A.
Поиск собственных значений матрицы T намного легче, чем для исходной матрицы, так как предполагается, что [math]k \lt \lt n[/math], и матрица T - трехдиагональная.
Полная переортогонализация необходима для того, чтобы гарантировать, что каждый полученный вектор qj+1 ортогонален уже имеющимся векторам q1..j. Без этого процесса будут накапливаться существенные вычислительный ошибки.
1.3 Вычислительное ядро алгоритма
Вычислительным ядро алгоритма состоит ииз двух основных частей:
- [math]Aq=( \sum\nolimits_{i=^n}a_{1i}q_i, \sum\nolimits_{i=2}^na_{2i}q_i, ..., \sum\nolimits_{i=1}^na_{ni}q_i)[/math] - умножение симметричной матрицы [math]A[/math] размерности [math]n*n[/math] на вектор q размерности n.
- [math]z=z-\sum\nolimits_{i=1}^{k}(z^Tq_i)q_i.[/math] - процесс ортогонализации Грама-Шмидта.
1.4 Макроструктура алгоритма
В указанном алгоритме можно выделить следующие основные операции:
- умножение матрицы на вектор;
- вычисление нормы вектора;
- скалярное произведение векторов;
- вычисление собственных значений и соответствующих им векторов симметричной трехдиагональной матрицы.
1.5 Схема реализации последовательного алгоритма
В параграфе 1.2 приводится полная схема последовательного алгоритма.
Заполняем начальные значения алгоритма (b - начальное преближение).
[math] \begin{align} & q_1=b/||b||,\\ & \beta_1=0,\\ &q_0=0, \\ \end{align} [/math]
Для всех [math]j=1..k[/math]:
- 1. Вычисляется j-й диагональный элемент матрицы [math]T_k[/math]: [math]z=Aq_j; \alpha_j=q_j^Tz;[/math]
- 2. Проводится полная переортогонализация Грамма-Шмидта: [math]z =z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i;[/math]
- 3. Вычисляются значения [math]\beta_{j+1}[/math] матрицы [math]T_k[/math]: [math]\beta_{j+1}=||z||;[/math]
- 4. Если [math]\beta_{j+1}=0[/math], то алгоритм завершается;
- 5. Сохраняем значения для следующей итерации [math]q_{j+1}=z/\beta_{j+1}.[/math]
1.6 Последовательная сложность алгоритма
- 1. умножение симметричной матрицы [math]A[/math] размерности [math]n*n[/math] на вектор q размерности n - вычислительная сложность: [math]n^2[/math] умножений и [math]n^2-n[/math] сложений
- 2. процесс ортогонализации Грама-Шмидта - вычислительная сложность: [math]k^2n+k(n+2)[/math] умножений и [math]k^2n + k(n + 1) + 2[/math] сложений.
1.7 Информационный граф
1.8 Ресурс параллелизма алгоритма
1.9 Входные и выходные данные алгоритма
1.10 Свойства алгоритма
1.11 Локальность данных и вычислений
1.11.1 Локальность реализации алгоритма
1.11.1.1 Структура обращений в память и качественная оценка локальности
1.11.1.2 Количественная оценка локальности
1.12 Возможные способы и особенности параллельной реализации алгоритма
1.13 Масштабируемость алгоритма и его реализации
1.13.1 Масштабируемость алгоритма
1.13.2 Масштабируемость реализации алгоритма
1.14 Динамические характеристики и эффективность реализации алгоритма
1.15 Выводы для классов архитектур
1.16 Существующие реализации алгоритма
2 Литература
<references \>