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

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

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

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

История

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

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

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

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

[math]\displaystyle{ K_n = [b,Ab,A^2b,...,A^{n-1}b] }[/math].

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

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

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

[math]\displaystyle{ z=Aq_j }[/math].

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

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

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

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

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

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

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

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

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

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

TODO

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

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

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

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

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

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