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

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

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

Метод Ланцоша — процедура, применяемая для нахождения [math]\displaystyle{ k }[/math] собственных значений симметричной матрицы произвольного размера путём вычисления их приближений. Суть метода заключается в итеративном построении специальной матрицы [math]\displaystyle{ T_{j} }[/math], а затем вычислении её собственных значений и собственных векторов.

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

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

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

• Умножение квадратной матрицы порядка [math]\displaystyle{ n }[/math] на вектор длины [math]\displaystyle{ n }[/math] и наоборот требует по [math]\displaystyle{ n^2 }[/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 }[/math] QR-алгоритмом требует [math]\displaystyle{ O(k^2) }[/math] операций, вычисление также и собственных векторов требует примерно [math]\displaystyle{ 6 * k ^ 3 }[/math], то есть [math]\displaystyle{ O(k^3) }[/math] операций; при использовании метода «Разделяй-и-властвуй» для вычисления собственных значений и векторов аналогичной матрицы в среднем затрачивается [math]\displaystyle{ O(k^{2.3}) }[/math] операций.

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

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

Входные данные: квадратная матрица [math]\displaystyle{ A }[/math] ( элементы [math]\displaystyle{ a_{ij} }[/math] ), вектор [math]\displaystyle{ b }[/math], число [math]\displaystyle{ k }[/math]
Объём входных данных: [math]\displaystyle{ n^2 + n + 1 }[/math]
Выходные данные: вектор [math]\displaystyle{ L }[/math], матрица [math]\displaystyle{ E }[/math] ( элементы [math]\displaystyle{ e_{pq} }[/math] )
Объём выходных данных: [math]\displaystyle{ k + k * n }[/math]

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

В случае неограниченности ресурсов соотношение последовательной и параллельной сложности алгоритма Ланцоща по сумме всех итераций представляет собой [math]\displaystyle{ k ^ 3 + k * n ^ 2 }[/math] к [math]\displaystyle{ X + k * n }[/math], где [math]\displaystyle{ X }[/math] - параллельная сложность алгоритма, используемого для вычисления собственных значений. В случае использования метода вычисления собственных значений с линейной сложностью при классисифкации по высоте ЯПФ это соотношение будет равно [math]\displaystyle{ (k ^ 2 + n ^ 2) / (k + n) }[/math].

Алгоритм Ланцоша не является полностью детерминированным в том смысле, что возможно выполнение числа итераций алгоритма, меньшего, чем заданное ( из-за того, что вычислены все ненулевые собственные значения матрицы [math]\displaystyle{ A }[/math] ).

Важное свойство метода Ланцоша состоит в том, что первыми в матрице [math]\displaystyle{ T_{j} }[/math] появляются собственные значения с максимальной величиной по модулю. Таким образом, метод особенно хорошо подходит для вычисления собственных значений матрицы [math]\displaystyle{ A }[/math], находящихся на краях её спектра.

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

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

Для проверки масштабируемости алгоритма из реализации была исключена непосредственно задача поиска собственных значений матрицы [math]\displaystyle{ T_{j} }[/math] на каждой итерации, так как это является предметом отдельных статей ( QR-алгоритм, метод «Разделяй-и-властвуй», итерации Арнольди и.т.д. ). Таким образом, проверяемая на масштабируемость задача сводится к итерационному вычислению коэффициентов матрицы [math]\displaystyle{ T_{j} }[/math]. Матрица [math]\displaystyle{ A }[/math] в рассматриваемой реализации хранится построчно и построчно же разделяется между процессорами.

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

Существует несколько как последовательных, так и параллельных реализаций алгоритма Ланцоша, включенных в программные библиотеки для поиска собственных значений матриц. Свойства этих реализаций можно увидеть в таблицах ниже. Первая таблица — относительно недавние реализации, вторая - «музейные экспонаты». Следует уделить особенное внимание столбцам «тип метода» и «параллелизация». В первом из них значение только N соответствует описанному варианту без реортогонализации, F – полной реортогонализации, P – частичной реортогонализации, S – выборочной реортогонализации. Во втором значние «none» соответствует отсутствию параллельной реализации, M – параллельной реализации посредством MPI, O – параллельной реализации посредством OpenMP. Их приведение целесообразно, так как на практике алгоритм Ланцоша без переортогонализации неустойчив. Также указываются библиотеки, в которых реализованы более глубокие модификации метода Ланцоша, с указанием изменений в графе «тип метода».

Таблица 1 - «современные» реализации.

Название библиотеки	Язык программирования	Дата появления, версия библиотки	Тип метода	Параллелизация
BLKLAN	C/Matlab	2003	P	none
BLZPACK	F77	2000, 04/00	P + S	M
IETL	C++	2006, 2.2	N	none
SLEPc	C/F77	2009, 3.0.0	All	M
TRLAN	F90	2006	Dynamic thick-restart	M
PROPACK	F77/Matlab	2005, 2.1 / 1.1	P, finds SVD	O
IRBLEIGS	Matlab	2000, 1.0	Indefinitie symmetric	none

Таблица 2 - «исторические» реализации.

Название библиотеки	Язык программирования	Дата появления, версия библиотки	Тип метода	Параллелизация
ARPACK	F77	1995, 2.1	Arnoldi iterations, impliicit restart	M
ARPACK++	C++	1998, 1.1	Arnoldi iterations, impliicit restart	none
LANCZOS	F77	1992	N	none
LANZ	F77	1991, 1.0	P	none
LASO	F77	1983, 2	S	none
NAPACK	F77	1987	N	none
QMRPACK	F77	1996	Designed for nonsymmetric matrices ( lookahead )	none
SVDPACK	C/F77	1992	P, finds SVD	none
Underwood	F77	1975	F, block version	none