Участник:Alexbashirov/Алгоритм Ланцоша для точной арифметики

Авторы: Баширов А.Д.

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

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

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

1.1.1 История

Алгоритм Ланцоша является прямым алгоритмом поиска собственных значений симметричной матрицы, разработанным венгерским физиком и математиком XX столетия Корнелием Ланцошем в 1950 году. Алгоритм совмещает в себе два других алгоритма: алгоритм построения подпространства Крылова (который также был создан Ланцошем) и процедуру Рэлея-Ритца приближения собственных значений матрицы. Алгоритм является вариацией степенного метода, позволяющего найти собственные значения и собственные вектора матрицы [math]n[/math]-го порядка. Однако в алгоритме число итераций ограничено некоторым [math]k[/math], где [math]k\lt \lt n[/math]. Несмотря на заявленную вычислительную эффективность, алгоритм не получил широкого применения в своём первоначальном виде, ввиду его неустойчивости. 20 лет спустя, в 1970 году ученые Ojalvo и Newman показали доработанный алгоритм Ланцоша, который теперь стал устойчивым а также предложили способ выбора вектора начального приближения [math]v_0[/math]. Полученный последователями Ланцоша алгоритм получил название "Алгоритм Ланцоша с полной переортогонализацией".

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

Алгоритм Ланцоша работает с вещественной симметричной матрицей [math]A=A^T[/math], таким образом в памяти достаточно хранить лишь немногим более половины элементов матрицы [math]A[/math].

Прежде чем переходить к описанию алгоритма, следует упомянуть следующие понятия: подпространство Крылова и процесс Рэлея-Ритца, о которых уже было сказано в исторической справке выше. Подпространством Крылова называется набор векторов[math]K_n =[b,Ab,A^2b,...,A^{n-1}b][/math]. Начав свою работу с произвольного вектора [math]b[/math], степенной метод вычисляет вектора [math]Ab, A^2b,...,A^{n-1}b[/math] итерационно сохраняя полученный результат в вектор [math]b[/math]. Полученная последовательность векторов сходится к собственному вектору, отвечающему максимальному собственному значению [math]\lambda_{max}[/math]. Однако, использование лишь последнего полученного вектора и потеря всех предшествующих посчитанных величин не всегда является наилучшим выбором. Поэтому сохранив все вектора и объединив их в одну матрицу мы получим так называемую матрицу Крылова:

[math]K_j = [b,Ab,A^2b,...,A^{j-1}b][/math].

Столбцы получившейся матрицы не ортогональны, однако их можно ортогонализовать, применив процесс Грамма-Шмидта. Получившиеся вектора будут составлять базис подпространства Крылова [math]Q_j[/math]. Можно ожидать, что вектора полученного базиса будут достаточно хорошо приближать [math]j[/math] собственных векторов, отвечающих [math]j[/math] наибольшим собственным значения.

Затем с использованием полученной матрицы из ортогональных векторов [math]Q_j[/math] мы получим трехдиагональную матрицу [math]T_j[/math] на основе исходной матрицы [math]A[/math]: [math]T_j=Q_j^TAQ_j[/math]. Наконец, к полученной матрице [math]T_j[/math] будет применяться процесс Рэлея-Ритца, который найдет собственные значения трехдиагональной матрицы. В последствии эти значения и будут приближениями искомых собственных значений исходной матрицы [math]A[/math].

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

Вычислительным ядром алгоритма Ланцоша, то есть наиболее ресурсно-затратной частью алгоритма, является шаг на котором на каждой итерации происходит перемножение матрицы [math]A[/math] на вектор [math]q_j[/math], полученный на предыдущей итерации. Полученный в результате умножения вектор обозначается [math]z[/math]:

[math] z=Aq_j[/math].

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

Как уже говорилось в математическом описании, алгоритм объединяет в себе два алгоритма: метод построения подпространства Крылова и процедуру поиска собственных значений трехдиагональной матрицы (процес Рэлея-Ритца).

Первую составную часть (алгоритм построения подпространства Крылова) представляет итерационный процесс построения матрицы [math]Q_j=[q_1,q_2,...q_j][/math] из [math]j[/math] ортонормированных вектров Ланцоша и последующее формирование трехдиагональной матрицы [math]T_j=Q_j^TAQ_j[/math]. Число таких векторов и, следовательно, размер матрицы растет с каждой итерацией. Итерационный процесс завершается при [math]j=k[/math].

На долю второй части алгоритма приходится поиск собственных значений и соответствующих им собственных векторов матрицы [math]T_j=Q_j^TAQ_j[/math], полученной в первой части.

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

На рисунке выше представлен псевдокод, отражающий последовательность действий в алгоритме метода. Перечислим основные ступени в данном методе:

• Прежде чем войти в цикл, происходит инициализация констант [math]\Beta_0,q_0[/math], вычисление нормы вектора начального приближения [math]b[/math] и нормализация этого вектора.
• Затем, в цикле, происходит: поиск вектора [math]z[/math], являющегося результатом применения матричного оператора [math]A[/math] к вектору [math]q_j[/math]; скалярное произведение двух векторов [math]q_j[/math] и [math]z[/math]; вычисляется линейная комбинация векторов; вычисляется норма вектора [math]z[/math]; проверяется условие выхода из цикла (равенство нулю нормы вектора [math]z[/math]); в случае если норма не нулевая, то производим нормировку вектора; затем происходит процедура поиска собственных значений и векторов матрицы [math]T_j[/math] и, наконец, переход на следующий шаг цикла.

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

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

• [math]n^2[/math] действий для умножения матрицы на вектор
• [math]n[/math] действия для перемножения двух векторов
• [math]n[/math] операций для поэлементного сложения или вычитания двух векторов
• для отыскания нормы вектора помимо [math]n[/math] сложений и [math]n[/math] умножений требуется 1 операция извлечения корня
• для поиска собственных значений и векторов матрицы [math]T_j[/math] требуется [math]O(k^3)[/math] операций.

Суммарно для одной итерации алгоритма требуется совершить [math]O(n^2)+O(k^3)[/math] операций.

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

В данном разделе приведен информационный граф одной итерации алгоритма Ланцоша. Для простоты на схеме подразумевается матрица размерности 5x5, однако это не нарушает общности данной схемы.

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

Поскольку алгоритм является итерационным, то выполнять параллельные операции возможно только внутри каждой итерации алгоритма. Сами итерации выполняются в строгой последовательности и не могут быть параллелизованны. Однако внутри каждой итерации происходят два процесса, которые могут быть параллелизованны:

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

Умножение квадратной матрицы размерности [math]n\times n[/math] на вектор размерности [math]n[/math] требует [math]n[/math] ярусов умножений и сложений. В свою очередь ресурс параллелизма операции поиска собственных значений зависит от выбранного способа и не рассматривается в данной статье.

Таким образом, с точки зрения ресурса параллелизма, алгоритм Ланцоша обладает линейной сложностью по ширине ЯПФ и квадратичной сложностью по высоте ЯПФ. Основное время работы алгоритма будет уделено операции умножения матрицы на вектор, имеющей ширину яруса [math]n^2[/math].

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

Алгоритм призван найти фиксированное число [math]k[/math] первых собственных значений матрицы [math]A[/math], следовательно на вход данный алгоритм во всяком случае получает 1 число и 1 матрицу размерности [math]n\times n[/math]. Помимо этого для начала работы алгоритму нужен некоторый вектор [math]b_0[/math] начального приближения размерности [math]n[/math], который тоже является входными данными.

После окончания работы алгоритма мы получаем набор собственных значений матрицы [math]A[/math], состоящий из [math]k[/math] элементов. Если помимо собственных значений выводятся еще и полученные собственные векторы, то суммарный объём всех таких векторов составляет [math]k\times n[/math] элементов.

Таким образом, совокупный объём входных данных составляет: [math]n^2+n+1[/math] элемент, а на выходе получаем [math]k(n+1)[/math] элементов.

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

Преимуществом алгоритма является то, что он начинает поиск собственных значений матрицы [math]T_j[/math] (а значит и матрицы [math]A[/math]) начиная с максимального в абсолютном смысле значения. Это выгодно отличает данный алгоритм с точки зрения поиска собственных значений матрицы, находящихся по краям её спектра.

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

В случае использования метода вычисления собственных значений с линейной сложностью, соотношение последовательной и параллельной сложностей будет равно [math](k ^ 2 + n ^ 2) / (k + n)[/math].

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

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

TODO

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

В настоящее время алгоритм Ланцоша поиска собственных значений квадратной симметричной матрицы включён в несколько библиотек и реализован на различных языках, среди которых такие наиболее распространенные как C, C++, FORTRAN77/90, MathLab. Так, например, в языке MathLab во встроенном пакете ARPACK найти собственные значения методом Ланцоша можно при помощи вызова функции eigs().

На языке С существует библиотека PRIMME, название которой расшифровывается как PReconditioned Iterative MultiMethod Eigensolver (итерационные методы поиска собственных значений с предусловием).

Библиотека NAG Numerical Library также содержит обширный набор функций, позволяющих проводить математический анализ, среди которых есть и реализация алгоритма Ланцоша. Библиотека доступна в трех пакетах: NAG C Library, NAG Fortran Library и NAG Library for .NET, совместна со многими языками и с основными операционными системами (Windows, Linux and OS X, а также Solaris, AIX и HP-UX).

3 Литература

[1] Дж. Деммель «Вычислительная линейная алгебра», с 381.

[2] Wikipedia https://en.wikipedia.org/wiki/Lanczos_algorithm

[3] "Investigation and Implementation of an Eigensolver Method", Jorge Moreira, The University of Edinburgh, August 2014

[4] "Lanczos Vector Procedures", Travis Hoppe