Участник:Nikita270a
Алгоритм Ланцоша с полной переортогонализацией | |
Последовательный алгоритм | |
Последовательная сложность | O(n^2)+O(nk^2) |
Объём входных данных | \frac{n(n + 1)}{2} |
Объём выходных данных | k(n + 1) |
Авторы: Алексейчук Н.Н 616, Ястребов К.С. 609.
Содержание
- 1 Общее описание алгоритма.
- 2 Математическое описание алгоритма
- 3 Вычислительное ядро алгоритма
- 4 Макроструктура алгоритма
- 5 Схема
- 6 Вычислительная сложность
- 7 Информационный граф
- 8 Входные и выходные данные
- 9 Свойства алгоритма
- 10 Реусурс параллелизма алгоритма
- 11 Библиотеки реализующие алгоритм
- 12 Литература
1 Общее описание алгоритма.
Алгоритм Ланцоша ищет собвственные значения и собственные векторы для симетричной матрицы A вещественных чисел. Является итерацонным алгоритмом. Алгоритм Ланцоша использует степенной метод ( b_{k+1} = \frac{Ab_k}{||Ab_k||} ) для поиска наибольших собственных значений и векторов матриц.
В отличие от прямых алгоритмов требует мешьше памяти и мощности, что является несомненным плюсом для больших матриц.
Несмотря на свою скорость работы и экономию памяти, сначала не был популярным алгоритмом из – за недостаточной вычислительной устойчивости. В 1970 году Ojalvo и Newman [1] показали способ сделать алгоритм достаточно устойчивым. В этой же статье алгоритм был применен к расчету инженерной конструкции с большим количеством узлов, которые подвергались динамической нагрузке.
2 Математическое описание алгоритма
Памятка: Степенной метод нахождения наибольшего собственного числа матрицы можно сформулировать в предельном виде: если b_0 – случайный вектор, и b_n+1 = Ab_n , тогда для больших чисел n предел x_n/||x_n|| стремится к нормированному наибольшему собственному вектору.
Алгоритм Ланцоша комбинирует метод Ланцоша для нахождения крыловского подпространства и метод Релэя – Ритца.
Подпространство Крылова для степенного метода: K_m(v,A) = span[x_1, Ax_1, A^2x_1, ..., A^{k-1}x_1]
В качестве входных данных для алгоритма Ланцоша подаются квадратная матрица размерности nXn: A=A^T; а так же вектор начального приближения b.
Метод осуществляет поиск трехдиагональной симметричной матрицы T_k=Q_k^TAQ_k. Причем собственные значения T_k таковы, что приближают собственные значения исходной матрицы A. То есть на каждом k-м шаге из ортонормированных векторов Ланцоша строится матрица Q_k = [q_1,q_2,...,q_k] и в качестве приближенных собственных значений матрицы A принимаются числа Ритца.
Из-за ошибок округления вектора A^{k-1}x_1 , формирующие подпространство Крылова, становятся неортогональными. Чтобы решить данную прорблему, проводят переортогонализацию методом Грамма-Шмидта.
3 Вычислительное ядро алгоритма
Вычислительным ядром данного алгоритма являются следующие шаги:
- 1. 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).
- 2. z=z-\sum\nolimits_{i=1}^{k}(z^Tq_i)q_i.
В данном случае первая операция выполняет умножение симметричной матрицы A размерности nXn на n-мерный вектор q, вследствие чего вычислительная сложность выполнения заключается в n^2 умножений и n^2-n сложений. Второе действие является процессом ортогонализации Грама-Шмидта. В последнем действии вычисляются k^2n+k(n+2) умножений и k^2n + k(n + 1) + 2 операций сложения.
4 Макроструктура алгоритма
- 1. Скалярное произведение: (x,y) .
- 2. Суммирование: x+\alpha y .
5 Схема
Input: A, b (random vector with unit norm)
- \begin{align} q_1 = b/||b||_2, \beta_0 = 0, q_0 = 0 \\ j = 1 ,...,k \\ q_1&=b/||b||,\beta_0=0,q_o=0. \\ z&=Aq_j, \\ \alpha_j&=q_j^Tz, \\ z&=z-\alpha_jq_j-\beta_{j-1}q_{j-1}, \\ \beta&=||z||,\\ q_{j+1}&=z/\beta_j, \quad j \in [1, k]. \end{align}
6 Вычислительная сложность
Количество операций складывается из количества операций для классического метода Ланцоша[2] и и количества операций для переортогонализации методом грамма-Шмидта[3].
Итоговая сложность составляет O(n^2)+O(nk^2) .
7 Информационный граф
8 Входные и выходные данные
На вход принимается сама матрица A \in R^{n \times n}, случайный вектор b (впрочем возможен вариант, при котором этот вектор генерируется случайным образом). Ввиду симметричности входной матрицы достаточно передать лишь ее верхнюю треугольную матрицу, что дает объем входных данных \frac{n(n + 1)}{2}.
Выходными данными для jой итерации являются вектор Aq_j и собственное число \lambda . Общий объем выходных данных k(n+1).
9 Свойства алгоритма
Если A эрмитова матрица, то алгоритм Ланцоша и Bi-Lanczos сходятся к одинаковым трехдиагональным матрицам ритца[4].
При реализации классического алгоритма Ланцоша возникает большая погрешность при округлении. Вариант с полной переортогонализацией позволяет избегать больших погрешностей, однако является более ресурсоемким. Существует промежуточный вариант с частичной переортогонализацией.
Алгоритм может завершить свою работу досрочно, когда найденные собственные значения будут достаточно близки к целевым.
10 Реусурс параллелизма алгоритма
Хотя алгоритм является итерационным, возможно распараллелить внутри каждой итерации умножение матрицы на вектор и процесс переортогонализации Грамма-Шмидта.
Процесс умножения матриц matvec можно распараллелить несколькими способами[5].
11 Библиотеки реализующие алгоритм
The IETL Project http://www.comp-phys.org/software/ietl/ C++
NAG Library http://www.nag.com/content/nag-library C, C++, Fortran, C#, MATLAB, R
ARPACK https://people.sc.fsu.edu/~jburkardt/m_src/arpack/arpack.html MATLAB
GrapLab https://turi.com/products/create/open_source.html C++
С частичной переортаганализацией
LANSO/PLANSO http://web.cs.ucdavis.edu/~bai/ET/lanczos_methods/overview_PLANSO.html Fortran (уже распараллелена)
12 Литература
Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).