Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

Алгоритм Ланцоша для точной арифметики (без ортогонализации)
Последовательный алгоритм
Последовательная сложность	$O(kn^2)$
Объём входных данных	$\frac{n(n+1)}{2} + 1$
Объём выходных данных	$kn+k$
Параллельный алгоритм
Высота ярусно-параллельной формы	$O(k*log(n))$
Ширина ярусно-параллельной формы	$O(n^2)$

Авторы:

• А.В. Шостик.

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

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

Алгоритм Ланцоша - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача $Ax = \lambda x$). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

Алгоритм был предложен в 1950 году венгерским физиком и математиком Корнелием Ланцошем. Основан на понятии построения крыловского подпространства и процедуре Рэлея-Ритца: строит последовательность трехдиагональных матриц с тем свойством, что экстремальные собственные значения для каждой последующей трехдиагональной матрицы дают все более точные оценки экстремальных собственных значений для А.

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией).

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

Исходные данные:

• симметрическая матрица $A=A^T$ размерности $n$, элементы матрицы $a_{ij}=a_{ji}$
• начальный вектор $v = {v_1, v_2, ..., v_n}, v \neq \theta$

Так как матрица является симметрической, достаточно хранить $\frac{n(n+1)}{2}$ ее элементов.

Для исходной матрицы $A$ и вектора $v$ строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов $[v, Av, A^2v, ..., A^{k-1}v]$, называемых векторами Ланцоша. Каждый из векторов Ланцоша ортонормируется: $q_k=\frac{A^kv}{\|A^kv\|}$, и на шаге $k$ алгоритма имеется Крыловская матрица $Q = {q_0, ..., q_{k-1}}$ размерности $n \times k$.

Соответствующий текущей итерации $k$ вектор Ланцоша считается $k$-ым приближением собственного вектора исходной матрицы $A$. В качестве приближенных собственных значений матрицы $A$ берутся собственные значения симметричной трехдиагональной матрицы $T_k = Q^T_k A Q$ (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует $O(k^3)$ операций..

Выходные данные:

• совокупность собственных векторов матрицы $T_k$ (т.е. $k$-ых приближений собственных векторов матрицы $A$)
• совокупность собственных значений матрицы $T_k$ (т.е. $k$-ых приближений собственных значений матрицы $A$)

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

Вычислительное ядро. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора $z=Aq_i$

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

Итерация алгоритма:

• вычисление нормы вектора,
• деление вектора на число,
• умножение матрицы на вектор,
• линейная комбинация векторов (умножение на число и сложение).

Финальный расчет:

• вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

Псевдокод алгоритма:

Вход  : размерность матрицы [math]n[/math], элементы симметричной матрицы [math]A[/math], количество итераций [math]k[/math]
Выход : собственные вектора матрицы [math]T_k[/math], матрица собственных значений 

$\begin{align} q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\ for \; & i = 1 \; to \; k \\ & z = Aq_i\\ & \alpha_i = q^T_i z\\ & z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\ & \beta_i = \|z\|_2\\ & If \; \beta_i == 0 \; then \\ & \; \; \; \; exit\\ & else \\ & \; \; \; \; q_{i+1} = z / \beta_i \\ end \; & for \end{align};$

$procedure(T_k)$;

$procedure()$ - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

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

Глобально алгоритм состоит из инициализации, $k$ итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:

• вычисление нормы вектора - $(n-1)$ сложений, $n$ умножений, операция извлечения корня
• деление вектора на число - $n$ операций
• умножение матрицы на вектор - $n^2$ операций
• скалярное произведение векторов - $n$ операций
• линейная комбинация векторов - $2n$ умножений и вычитаний
• вычисление собственных векторов и собственных значений матрицы $T_k$ методом "Разделяй и властвуй" - $O(k^3)$ операций.

Инициализация:

• вычисление нормы вектора,
• деление вектора на число.

Одна итерация:

• умножение матрицы на вектор,
• скалярное произведение векторов,
• линейная комбинация векторов,
• вычисление нормы вектора,
• деление вектора на число.

Финальный расчет:

• вычисление собственных векторов и собственных значений матрицы

Итого последовательная сложность алгоритма составляет $O(3n+k(n^2+n+2n+2n+3n))+O(k^3)$.

Так как на практике порядок матрицы $n$ намного превышает количество итераций $k$, сложность алгоритма Ланцоша можно рассматривать $O(kn^2)$.

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

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

Рисунок 1. Информационный граф алгоритма Ланцоша.
$\|.\|$ — операция вычисления нормы,
$/$ — операция деления,
$A$ — операция умножения матрицы на вектор,
$*$ — операция скалярного произведения,
$L$ — операция вычисления линейной комбинации векторов,
▽ - проверка условия $\mathsf{(\beta_i = 0)}$ и выход из цикла в случае выполнения условия

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

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации $n$ умножений и $n-1$ сложение. В таком виде сложение $n$ элементов можно выполнить за $\log_2 n$ действий.

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

Входные данные: симметрическая матрица $A=A^T$ размерности $n$, количество итераций $k$.

Объем входных данных: $\frac{n(n+1)}{2} + 1$.

Выходные данные: $k$ собственных векторов, собственных значений матрицы на $k$-ом приближении.

Объем выходных данных: $kn+k$.

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

• Сложность последовательной реализации алгоритма $O(kn^2)$.
• Сложность параллельного алгоритма Ланцоша по высоте ЯПФ $O(k\log n)$, по ширине ЯПФ $O(n^2)$.
  • В умножении матрицы на вектор сложение n элементов выполняется в $\log_2 n$ ярусов шириной $\frac{n}{2}, \frac{n}{4}, ... , 1$, остальные операции внутри итерации выполняются последовательно
• Отношение последовательной сложности к параллельной $\frac{kn^2}{k \log{n}}$.
• Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма $\frac{k(2n^2+8n-1)+3n}{n^2+2k}$. Учитывая $k$ много меньше $n$, вычислительная мощность ≈ $2k$.
• В алгоритме Ланцоша возможно выполнение числа итераций менее $k$ если все собственные значения матрицы вычисляются раньше.
• Для найденных собственных значений справедливо: $\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})$, т.е. в первую очередь находятся максимальные по модулю собственные значения.
• Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.