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

Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Konshin и ASA.

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

Алгоритм Ланцоша для точной арифметики
Последовательный алгоритм
Последовательная сложность	[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 Общее описание алгоритма

Алгоритм Ланцоша является итерационным алгоритмом поиска собственных значений симметричной матрицы, разработанным венгерским физиком и математиком 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]q=b/||b||_2[/math], [math]\beta_0=0[/math], [math]q_0=0[/math]

[math] ~ {\rm for } j=1 ~ {\rm to} k[/math]

[math]z = Aq_j[/math]

[math]\alpha_j = q_j^T z[/math]

[math] z = z - \alpha_jq_j -\beta_{j-1}q_{j-1}[/math]

[math] \beta_j = ||z||_2[/math]

[math] ~ {\rm if} \beta_j = 0 ~ {\rm quit}[/math]

[math] q_{j+1} = z/\beta_j[/math]

Вычисление собственных значений и векторов матрицы [math]T_j[/math]

[math] ~ {\rm end} ~ {\rm for} [/math]

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

• Прежде чем войти в цикл, происходит инициализация констант [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 Информационный граф

В данном разделе представлен информационный граф алгоритма, дающий представление о переходе от одной итерации к другой.

Схема состоит из трех частей, каждая из которых символизирует свой итерационный слой. Первый и второй слои являются стартовыми - на них задаются начальные приближенные значения параметров, с которыми алгоритм может быть запущен. Далее работа с уточнением значений этих параметров идет в рамках итерационного цикла на третьем слое. Такие входные данные, как матрица оператора [math]A[/math] используются на каждом итерационном слое, поэтому обозначена как входные данные на втором и третьем слоях на схеме. Те данные, которые передаются на следущий итерационный слой с предыдущего - обозначены стрелками между вторым и третьим слоями.

[math]||x||[/math] - операция вычисления нормы; [math]/[/math] - операция деления; • - операция скалярного умножения векторов; Ax - операция применения оператора A к вектору x; F - операция вычисления линейной комбинации. Желтый ромб символизирует проверку логического условия выхода из итерационного цикла (в случае если посчитанная норма равна нулю). Input,Out - входные и выходные данные.

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

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

• умножение матрицы на вектор,
• процесс поиска собственных значений. Этот процесс является самостоятельной задачей. Для поиска собственных значений могут использоваться различные методы, например такие, как QR-алгоритм или метод "Разделяй и властвуй". Соответственно, в зависимости от выбранного алгоритма поиска с./знач. могут иметь различные свойства.

Умножение квадратной матрицы размерности [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]\frac{k(2n^2+8n-1)+3n}{n^2+2n+1}[/math]. Следовательно порядок мощности алгоритма примерно равен [math]2k[/math], при [math]k\lt \lt n[/math]

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

Алгоритм Ланцоша может быть модифицирован так, чтобы быть пригодным для применения к несимметричным системам.

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

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

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

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

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

Для компиляции программы использовались компиляторы intel/15.0.090 и Openmpi/1.8.4-icc.

Для проведения расчетов и получения полноценной картины поведения алгоритма в зависимости от входных данных и числа процессоров, программа была запущена на следующих параметрах:

• в качестве размера входной матрицы подавались значения в диапазоне [2000:30000] c шагом 2000.
• Число процессоров варьировалось от 1 до 256.
• В качестве числа, отвечающего за количество выполняемых методом итераций бралось значение 100.

Ниже на рисунке изображен трехмерный график, показывающий зависимость времени выполнения программы от входных данных/

График зависимости времени выполнения программы от входных данных

Также была исследована эффективность распараллеливания:

• Средняя эффективность выполнения алгоритма составила порядка 12%
• Наивысший показатель равен 43%
• Минимальное значение составило <1%

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

[5] "Parallel Scientific Computing in C++ and MPI", A seamless approach to parallel algorithms and their implementation, George Em Karniadakis and Robert M. Kirby II, Cambridge University Press