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

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

Основные авторы описания: Д.Р.Слюсарь

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

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

Алгоритм Ланцоша поиска собственных значений был опубликован Корнелием Ланцошем в 1950 году [1]. Этот итерационный алгоритм применим только к эрмитовым матрицам [math]\displaystyle{ A }[/math]. Метод позволяет за [math]\displaystyle{ k }[/math] итераций вычислять [math]\displaystyle{ k }[/math]-ое приближение собственных значений и собственных векторов исходной матрицы [math]\displaystyle{ A }[/math].

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

Данный алгоритм является неустойчивым, вследствие чего на практике применяется модифицированный алгоритм Ланцоша с полной переортогонализацией.

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

На вход алгоритма подается эрмитова матрица [math]\displaystyle{ A = A^\dagger }[/math] (в вещественном случае - симметрическая) ,

[math]\displaystyle{ A = \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1\ n-1} & a_{1\ n} \\ a_{12} & a_{22} & a_{23} & \cdots & a_{2\ n-1} & a_{2\ n} \\ a_{13} & a_{23} & a_{33} & \cdots & a_{3\ n-1} & a_{3\ n} \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ a_{1\ n-1} & \cdots & \cdots & a_{n-2\ n-1} & a_{n-1\ n-1} & a_{n-1\ n} \\ a_{1\ n} & \cdots & \cdots & a_{n-2\ n} & a_{n-1\ n} & a_{n\ n} \\ \end{pmatrix} }[/math]

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца. Иными словами, из оргонормированных векторов Ланцоша [??] на каждой итерации строится матрица [math]\displaystyle{ Q_k = [q_1, q_2, \dots, q_k] }[/math] размерности [math]\displaystyle{ n \times k }[/math]. В качестве приближенных собственных значений матрицы [math]\displaystyle{ A }[/math] берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы [math]\displaystyle{ T_k = Q^T_k A Q }[/math]:

[math]\displaystyle{ T_k = \begin{pmatrix} \alpha_1 & \beta_1 & 0 & \dots & 0 \\ \beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\ 0 & \beta_2 & \ddots & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\ 0 & \dots & \dots & \beta_{k-1} & \alpha_k \end{pmatrix} }[/math]

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы [math]\displaystyle{ A }[/math] на вектор [math]\displaystyle{ q_i }[/math] с предыдущей итерации

[math]\displaystyle{ z = Aq_i }[/math]

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

Макрооперациями в алгоритме являются:

• Процедура иттеративного построения трехдиагональной симметричной матрицы, включающая:
  • умножение матрицы на вектор;
  • скалярное произведение векторов;
  • деление (умножение) вектора на вещественное число;
• Вычисление собственного значение и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

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

Исходные данные: симметричная матрица [math]\displaystyle{ A }[/math], случайный вектор [math]\displaystyle{ b }[/math].

Вычисляемые данные: собственные вектора матрицы [math]\displaystyle{ T_k }[/math] являющиеся столбцами матрицы [math]\displaystyle{ Q_k V }[/math], и матрица собственных значений [math]\displaystyle{ \Lambda }[/math], где [math]\displaystyle{ V, \Lambda }[/math] из спектрального разложения [math]\displaystyle{ T_k = V\Lambda V^T }[/math].

Алгоритм на псевдокоде:

[math]\displaystyle{ \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} }[/math]

После этого вычисляются собственные значения и собственные вектора симметричной трехдиагональной матрицы [math]\displaystyle{ T_k }[/math] наиболее удобным образом.

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:

• Умножение квадратной матрицы [math]\displaystyle{ n * n }[/math] на вектор длины [math]\displaystyle{ n }[/math]. Требует [math]\displaystyle{ n * n }[/math] умножений и сложений;
• Скалярное произведение векторов длины [math]\displaystyle{ n }[/math]. Требует [math]\displaystyle{ n }[/math] умножений и сложений;
• Сложение векторов длины [math]\displaystyle{ n }[/math]. Требует [math]\displaystyle{ n }[/math] сложений;
• Умножение вектора длины [math]\displaystyle{ n }[/math] на скаляр. Требует [math]\displaystyle{ n }[/math] умножений;
• Нахождение квадратичной нормы вектора длины [math]\displaystyle{ n }[/math]. Требует [math]\displaystyle{ n }[/math] умножений и сложений + извлечения квадратного корня;
• Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера [math]\displaystyle{ k \times k }[/math]. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует [math]\displaystyle{ ~O(k^{2.3}) }[/math] операций.

Таким образом, в худшем случае, алгоритм Ланшоца имеет сложность [math]\displaystyle{ O(k * n^2) }[/math] и односится к квадратичным алгоритмам.

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

Информационный граф алгоритма можно разбить на две части:

• Общий граф алгоритма, от входа Init до функции вычислений собственных значений полученной трехдиагональной матрицы (Fin);
• Подграф на каждой итерации внутри цикла;

Общий граф алгоритма. Init - вход, Iter - одна итерация алгоритма, Fin - процедура вычисления собственных значений

Подграф алгоритма, описывающий одну итерацию. IN - вход в очедедную итерацию, NEXT - переход к следующей итерации или выход цикла, между ними возможен параллелизм при операциях умножения/сложения, вычитания из q(i) q(i-1)-ого, а также деления нового вектора на норму (полученную в блоке normalising)

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

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

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

Входные данные: симметричная вещественная матрица [math]\displaystyle{ A }[/math], случайный вектор [math]\displaystyle{ b }[/math], число итераций [math]\displaystyle{ k }[/math]

Объём входных данных: [math]\displaystyle{ n * (n + 1) + 1 }[/math]

Выходные данные: вектор собственных значений [math]\displaystyle{ \Lambda }[/math], матрица собственных векторов [math]\displaystyle{ E }[/math]

Объём выходных данных: [math]\displaystyle{ k * (n + 1) }[/math]

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

TDB

Программная реализация алгоритма

Масштабируемость алгоритма и его реализации

TDB next week

Существующие реализации алгоритма

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (ARPACK), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project [2]

2. MatLab [3]

3. ARPACK [4]

4. Julia Math [5]

...

Литература

1. Алгоритм Ланцоша (Википедия) [6]

2. Деммель Д. Вычислительная линейная алгебра/Пер. с англ. ХД Икрамова. – 2001.