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

Алгоритм Ланцоша с полной переортогонализацией
Последовательный алгоритм
Последовательная сложность	[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 Общее описание алгоритма

Алгоритм Ланцоша - итерационный метод, используемый для вычисления части собственных значений и соответствующих им собственных векторов матрицы [math]A[/math] размера [math]n*n[/math], изначально разработанный Корнелием Ланцошем. Преимуществами использования метода является относительно небольшое потребление памяти и вычислительных ресурсов, а также наличие параметра [math]k[/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]. Ее собственные значения приближают собственные значения исходной матрицы. Пусть T_k=V[math]\Lambda[/math]V^T - спектральное разложение матрицы T_k, тогда столбцы матрицы Q_kV рассматриваются как приближения к соответствующим собственным векторам матрицы A и называются векторами Ритца. Числа и векторы Ритца являются оптимальными приближениями к собственным значениям и собственным векторам матрицы A.

Поиск собственных значений матрицы T намного легче, чем для исходной матрицы, так как предполагается, что [math]k \lt \lt n[/math], и матрица T - трехдиагональная.

Полная переортогонализация необходима для того, чтобы гарантировать, что каждый полученный вектор q_j+1 ортогонален уже имеющимся векторам q_1..j. Без этого процесса будут накапливаться существенные вычислительный ошибки.

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

[math] \begin{align} q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\ for \; & j = 1 \; to \; k \\ & z = Aq_j\\ & \alpha_j = q^T_j z\\ & z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\ & 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\\ & if \; \beta_j =0\; break\\ & q_{j+1} = \frac{z}{\beta_j}\\ end \; & for \end{align} [/math]

2 Литература

<references \>