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

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

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

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

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

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

Вычислительное ядро алгоритма Ланцоща состоит в получении на каждой итерации очередного промежуточного вектора [math]z[/math], получаемого путём умножения исходной матрицы [math]A[/math] на вектор Ланцоша [math]q_{j}[/math], полученный на предыдущей итерации.

Также при значениях [math]k[/math], сопоставимых с [math]n[/math], процедура вычисления собственных значений и собственных векторов симметричной трёхдиагональной матрицы по некоторому выбранному алгоритму подразумевает значительный объём расчётов и может считаться вычислительным ядром. Однако на практике метод используется прежде всего при малых [math]k[/math].

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

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

• Умножение квадратной матрицы порядка [math]n[/math] на вектор длины [math]n[/math] и наоборот требует по [math]n^2[/math] умножений и сложений.
• Перемножение векторов длины [math]n[/math] требует по [math]n[/math] умножений и сложений.
  
• Поэлементное сложение векторов длины [math]n[/math] требует [math]n[/math] сложений.
• Умножение вектора длины [math]n[/math] на число требует [math]n[/math] умножений.
• Нахождение квадратичной нормы вектора длины [math]n[/math] требует по [math]n[/math] умножений и сложений, а также одну операцию извлечения квадратного корня.
• Вычисление собственных значений матрицы порядка [math]k[/math] QR-алгоритмом требует [math]O(k^2)[/math] операций, вычисление также и собственных векторов требует примерно [math]6 * k ^ 3[/math], то есть [math]O(k^3)[/math] операций; при использовании метода «Разделяй-и-властвуй» для вычисления собственных значений и векторов аналогичной матрицы в среднем затрачивается [math]O(k^{2.3})[/math] операций.

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

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

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

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

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

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

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

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

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

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

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

Существует несколько как последовательных, так и параллельных реализаций алгоритма Ланцоша, включенных в программные библиотеки для поиска собственных значений матриц. Свойства этих реализаций можно увидеть в таблицах ниже. Первая таблица — относительно недавние реализации, вторая - «музейные экспонаты». Следует уделить особенное внимание столбцам «тип метода» и «параллелизация». В первом из них значение только 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