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

Основные авторы описания: Левин А.Д. (студент, кафедра вычислительных методов, 604 группа) - все описанные разделы самостоятельно

Алгоритм Ланцоша для точной арифметики (без переортогонализации)
Последовательный алгоритм
Последовательная сложность	[math]O(k\,n^2)[/math]
Объём входных данных	[math]\frac{n^2+3n}{2} + 1[/math]
Объём выходных данных	[math]k\;(n + 1)[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(\log_2 n)[/math] для каждой итерации алгоритма. Максимальное число таких итераций = [math]k[/math], если не будет осуществлён ранний выход из цикла.
Ширина ярусно-параллельной формы	[math]O(n^2)[/math]

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

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

Данный алгоритм появился в 1950 г. и носит имя венгерского физика и математика Корнелия Ланцоша (венг. Lánczos Kornél). Алгоритм Ланцоша относится к итерационным методам вычисления собственных значений для матриц столь больших порядков [math] n[/math], что к ним нельзя применить прямые методы из-за ограничений по времени и памяти.

Сам Ланцош указывал, что его метод предназначен для отыскания нескольких собственных векторов симметричных матриц, хотя к методу сразу было обращено внимание, как к способу приведения всей матрицы к трёхдиагональному виду. Двадцатью годами позже канадский математик Крис Пэж показал, что, несмотря на чувствительность к округлениям, алгоритм Ланцоша - эффективное средство вычисления некоторых [math]k[/math] собственных чисел и соответствующих им собственных векторов [1, c.276]. Обычно число [math]k[/math] берут в два раза больше, чем количество собственных значений матрицы, которые поставлена цель найти, и оно имеет порядок не более [math]10^2[/math]

Алгоритм Ланцоша для вычисления собственных значений симметричной матрицы [math]A[/math] соединяет в себе метод Ланцоша для построения последовательности подпространств Крылова и ортонормированных векторов Ланцоша и процедуру Рэлея-Ритца получения оптимальных приближений собственных значений и соответствующих собственных векторов исходной матрицы [math]A[/math] [2, с.381].

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

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

Для лучшего понимания описания, данного в этом пункте статьи, рекомендуется ознакомиться с параграфом 6.6 Методы Крыловского подпространства [2, с.313]. Здесь же дано краткое описание всех переменных, математических операций и необходимый теоретический минимум.

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

[math] \begin{align} q_1 = & \,\frac{b}{\|b\|_2},\; \beta_0 = 0,\; q_0 = 0\\ for \; & j = 1 \; to \; k \\ & z = A\,q_j\\ & \alpha_j = q^T_j z\\ & z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\ & \beta_j = \|z\|_2\\ & If \; (\beta_j == 0) \; then\; stop\; the\; algorithm \\ & \; q_{j+1} = \frac{z}{\beta_j} \\ & compute\; eigenvalues\; and \;eigenvectors\;of \;matrix \;\;T_j= Q^T_j A Q\;\;and\;estimate \;the\; errors\\ end \; & for \end{align} [/math]

В продемонстрированном выше алгоритме [math]b[/math] - заданный вещественный вектор. Также полагается известным алгоритм вычисления произведения матрицы [math]A[/math] на вектор [math]x[/math]. Введём матрицу Крылова, определяемую следующим соотношением: [math]K_j = [b,Ab,A^2b,...,A^{j-1}b][/math].

Далее, на практике, матрица [math]K[/math] заменяется матрицей [math]Q[/math], такой, что при любом числе [math]k[/math] линейные оболочки первых [math]k[/math] столбцов в [math]K[/math] и [math]Q[/math] являются одним и тем же подпространством [2, c.315]. Тогда матрица [math]Q[/math], в отличие от матрицы [math]K[/math], хорошо обусловлена и легко обратима. В результате получаем матрицу [math]Q_j = [q_1, q_2, \dots, q_j][/math] размерности [math]n\times j[/math], столбцы которой ортогональны, являются базисом подпространства Крылова и называются векторами Ланцоша.

В алгоритме Ланцоша вычислению подлежит столько первых столбцов в матрице [math]Q_j[/math], сколько необходимо для получения требуемого приближения к решению [math]A\,x\,=b\,\,(A\,x=\lambda \, x)[/math].

Получение нулевого [math]\beta_j[/math] в итерационном процессе - событие, желательное в том смысле, что оно свидетельствует о том, что найдено некоторое подпространство, являющееся в точности инвариантным, а все собственные значения матрицы [math]T_j[/math] являются собственными значениями матрицы [math]A[/math]. На практике же нулевое и близкие к нулю значения получаются редко [3, стр. 429].

Затем на каждом шаге цикла формируется симметричная трёхдиагональная матрица [math]T_j = Q^T_j A Q[/math], к которой применяется процесс Рэлея-Ритца для поиска её собственных значений (на практике для поиска этих собственных значений можно использовать любой из специальных методов, изложенных в [3, [math]\S[/math] 8.4]). Эти собственные значения, они же числа Ритца, и полагаются приближёнными собственными значениями исходной матрицы [math]A[/math].

Если требуется вычислить собственные векторы, то необходимо хранить вычисляемые векторы Ланцоша. Если работа при этом производится только с оперативной памятью, то при большом порядке исходной матрицы (а именно для таких матриц и предназначен метод Ланцоша) число шагов [math]j[/math], которые могут быть выполнены до того момента, как оперативная память будет исчерпана, будет очень незначительно, а следовательно, искомые собственные пары могут не успеть сойтись с требуемой точностью. В этом случае метод Ланцоша можно использовать итерационным образом: после выполнения заданного числа шагов процесс прерывается и, если еще не все требуемые пары Ритца сошлись с требуемой точностью, то формируется новый начальный вектор [math]q_1[/math], являющийся линейной комбинацией векторов Ритца, еще не сошедшихся с заданной точностью, и процесс начинается заново.

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

Вычислительное ядро рассматриваемого алгоритма (наиболее ресурсно-затратная часть алгоритма) - формирование на каждом шаге цикла промежуточного вектора [math]z[/math], вычисляемого по формуле: [math]z = A\,q_j[/math]

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

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

В первой части алгоритма, где итерационно строятся Крыловское подпространство и трёхдиагональная матрица [math]T_j[/math], можно выделить следующие макрооперации:

• умножение матрицы на вектор (в свою очередь представляет собой умножение вектора на число и сложение векторов);
• вычисление нормы (скалярное произведение векторов + вычисление квадратного корня);
• скалярное произведение векторов;
• сложение векторов и их умножение на вещественные числа

Умножение матрицы на вектор - весьма дешевая и нетривиальная операция, если предполагать, что исходная матрица [math]A[/math] разреженная. Именно она подлежит оптимизации с использованием техник распараллеливания, поскольку она является вычислительным ядром алгоритма, а весь упомянутый процесс построения матрицы [math]T_j[/math] выполняется последовательно.

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

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

С целью полностью не дублировать пункт [math]\mathbf{1.2}[/math] данной статьи, где представлен очень наглядный и доступный для понимания псевдокод алгоритма, постараемся лишь коротко просуммировать и обобщить весь перечень последовательных действий в рассматриваемом алгоритме.

В качестве входных данных имеем: вещественная симметричная матрица [math]A[/math] размера [math]n\times n[/math], известный вещественный вектор [math]b[/math] длины [math]n[/math] и число [math]k\lt n[/math].

Далее:

• Инициализируем вещественную константу [math]\beta_0[/math] и векторы [math]q_0[/math] и [math]q_1[/math]
• Итерационно для всех [math]j=\overline {1,k}[/math]:
  • Вычисляем вспомогательный вектор [math]z[/math], используя умножение матрицы на вектор
  • Вычисляем, используя скалярное произведение векторов, значение [math]\alpha_j[/math] - элемент на главной диагонали трёхдиагональной матрицы [math]T_j[/math]
  • Вычисляем по формуле новое значение и перезаписываем вектор [math]z[/math]
  • Вычисляем значение [math]\beta_j[/math] - элементы непосредственно над и под главной диагональю той же матрицы [math]T_j[/math]
  • Вычисляем с вектора [math]z[/math] очередной столбец матрицы [math]Q[/math], он же вектор Ланцоша : [math]\,\,q_{j+1}[/math]
  • Теперь, когда матрица [math]T_j[/math] сформирована на очередном шаге, находим её собственные значения и вектора, например с помощью процедуры Рэлея-Ритца (о используемых на практике методах нахождения подробнее было сказано ранее в математическом описании алгоритма)
  • Переходим на следующую итерацию: [math]j[/math]++

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

Оценим количество операций и сложность последовательного алгоритма, изложенного в предыдущем пункте.

• Для вычисления вектора [math]z[/math] умножается матрица размера [math]n\times n[/math] на вектор длины [math]n[/math], и потребуется [math]n^2[/math] операций умножения и [math]n^2[/math] операций сложения.
• Умножение векторов [math]q^T_j[/math] и [math]z[/math] длины [math]n[/math] потребует [math]n[/math] операций умножения и [math]n[/math] операций сложения.
• При вычислении нового значения вектора [math]z[/math] поэлементное сложение двух векторов длины [math]n[/math] требует [math]n[/math] операций сложения.
• Умножение вектора длины [math]n[/math] на вещественное число требует [math]n[/math] операций умножения.
• Нахождение квадратичной нормы вектора [math]z[/math] длины [math]n[/math] требует [math]n[/math] умножений и [math]n[/math] сложений и ещё 1 операцию извлечения квадратного корня.
• Заключительное нахождение собственных значений и собственных векторов трёхдиагональной матрицы [math]T_j[/math] размера [math]j\times j[/math] с помощью QR-алгоритма потребует в худшем случае [math]O(j^3)[/math] операций.

Посчитаем теперь суммарное количество простейших операций на одной итерации алгоритма (без учёта записей в память):

• [math]n^2+5\,n[/math] операций умножения/деления
• [math]n^2+4\,n[/math] операций сложения/вычитания
• 1 операция извлечения квадратного корня
  

Отдельно после всех итераций при [math]j=k[/math] :

• [math]O(k^3)[/math] операций для поиска собственных значений и собственных векторов трёхдиагональной матрицы [math]T_k[/math]

Итого можно утверждать следующее:

Последовательная сложность одной итерации рассматриваемого алгоритма: [math]O(n ^ 2) [/math]
Суммарная сложность алгоритма ([math]k[/math] итераций): [math]O(k\,n ^ 2) + O(k ^ 3)[/math]

На практике же обычно рассматривается число [math]k\lt \lt n[/math], поэтому второе слагаемое можно опустить.

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

На рисунке продемонстрирован граф одной итерации (при каком-либо произвольном значении итератора [math]j[/math]) описанного алгоритма без отображения входных и выходных данных. Наглядно продемонстрированы связи между всеми действиями, и с помощью деления на [math]n[/math] потоков показано, какие операции возможно выполнять параллельно на каждом шаге с целью оптимизации и экономии машинного времени. В результате, после параллельно выполненных операций Division на заключительном ярусе, на выходе имеем значение следующего вектора Ланцоша [math]q_{j+1}[/math], после чего следует переход на следующую итерацию, и процесс повторяется.

Simp (от англ. simple) - операции сложения и умножения вещественных чисел при вычислении произведения матрицы на вектор;

Vec mult (от англ. vectors multiplication) - операция скалярного умножения двух векторов;

Vec subtr (от англ. vectors subtraction, из-за использованных в операции вычитаний векторов) - операция умножения вектора на число (выполняется дважды) и формирование нового значения [math]z_{new}[/math] по формуле: [math]z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}[/math]

[math]||\cdot||[/math] - операция вычисления нормы вектора;

Division - операция деления вектора на вещественное число [math]\beta_j[/math]

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

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

Параллельные операции возможны только внутри итераций, так как рассматриваемый алгоритм является итерационным, и сами итерации строго последовательны.
Рассмотрим в этом пункте подробнее, какой выигрыш мы можем получить при распараллеливании операций внутри одной итерации согласно графу из предыдущего пункта.

Умножение вещественной симметричной матрицы [math]A[/math] размера [math]n\times n[/math] на вектор [math]b[/math] длины [math]n[/math] требует последовательного выполнения [math]n[/math] ярусов умножений и сложений. Сложение элементов вектора длины [math]n[/math] на каждом шаге перемножения можно оптимизировать и выполнить за [math]\log_2 n[/math] действий, если применять метод сдваивания.

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

Ресурс параллелизма алгоритма вычисления собственных значений построенной трёхдиагональной матрицы [math]T_k[/math] размера [math]k\times k[/math] зависит от выбранного алгоритма (об этих алгоритмах подробнее было сказано в пункте 1.2). Будет очень затратно с точки зрения объёма материала пытаться рассмотреть пусть даже часть этих методов в этой статье, а так же тонкости их распараллеливания. Предлагаем читателям при желании самим ознакомиться с этими алгоритмами на данном ресурсе или прибегнув к иным источникам.

При классификации по высоте ЯПФ, таким образом, алгоритм Ланцоша без переортогонализации относится к алгоритмам с логарифмической сложностью, и его сложность равна [math]O(\log_2 n)[/math] с важным замечанием, что это сложность только одной итерации, а максимум в алгоритме может быть выполнено последовательно [math]k[/math] таких однотипных итераций.

Несмотря на то, что на графе в предыдущем пункте для лучшей визуализации и понимания было изображено [math]n[/math] потоков (никаких допущений мы не делали и нарушений с точки зрения логики и математики нет, можно изобразить большее число потоков, просто часть времени некоторые из них могут бездействовать и не использоваться), на самом деле правильным и разумным будет реализовать начальное параллельное выполнение [math]n^2[/math] операций умножения всех вещественных элементов матрицы [math]A[/math] на соответствующие элементы вектора [math]q_j[/math], после чего применять тот же самый упомянутый метод сдваивания.

Поэтому при классификации по ширине ЯПФ его сложность будет квадратичной и равна [math]O(n^2)[/math].

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

Входные данные алгоритма: вещественная симметричная [math]A[/math] размера [math]n\times n[/math], вещественный вектор [math]b[/math] длины [math]n[/math], и одно целое число [math]k[/math], означающее количество итераций.

Объём входных данных: [math]n^2 +n + 1 [/math]. Если учесть, что матрица симметричная и достаточно для её однозначного определения задать элементы на главной диагонали и элементы над или под ней, то можно сократить объём входных данных с целью экономии памяти до: [math]\ \frac{n^2-n}{2}+n+n+1 =\frac{n^2+3n}{2} + 1[/math]

Выходные данные алгоритма: вектор собственных значений [math]\Theta^k[/math] длины [math]k[/math], матрица [math]S^k[/math] из соответствующих собственных векторов размера [math]n\times k[/math].

Объём выходных данных: [math]k\;(n + 1)[/math].

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

В данном пункте перечислены без какой-либо хронологии или связи друг с другом наиболее важные и разнообразные свойства алгоритма, на которые стоит и полезно обратить внимание, а также которые выгодно отличают его от схожих алгоритмов, решающих ту же задачу.
1) Отношение последовательной сложности алгоритма к параллельной,с учётом [math]k[/math] итераций, равно [math]\frac{kn^2}{k \log_2{n}}=\frac{n^2}{\log_2{n}}[/math];
2) Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, равна [math]\frac{k\,n^2}{0.5\,n^2+1.5\,n+k\,n+k+1}[/math] и [math]\approx 2\,k[/math], так как на практике верно отношение [math]k\lt \lt n[/math] ;

3) Для собственных значений формируемых матриц [math]T_j[/math] верно отношение (следствие теоремы Коши о перемежаемости): [math]\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})[/math] ; [2, c.383]

4) Крайние (наибольшие и наименьшие) собственные значения матриц [math]T_j[/math] сходятся к крайним собственным значениям матрицы [math]A[/math] первыми (это часто происходит уже на ранних шагах при [math]j\approx 2\,\sqrt{n}[/math]), а собственные значения в середине спектра - последними;

5) Из-за ограниченной точности постепенно теряется ортогональность вектора Ланцоша [math]q_j[/math] к предыдущим векторам Ланцоша, поэтому необходимо примерно каждые 10-15 итераций проводить дорогостоящую операцию реортогонализации, но это уже будет иной алгоритм с очень малым различием в одной операции и в коде, реализующем его, соответственно;

6) Эффективность и быстрота выполнения алгоритма очень сильно зависит от скорости памяти (диска). Согласно тесту только 6% времени работы программы приходится на вычисления, остальное время затрачивается на ожидание отклика диска;

7) Метод Ланцоша удобен тем, что не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому метод Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.

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

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

Несмотря на то, что алгоритм является далеко не самым популярным, алгоритмом, про который нельзя сказать, что он всегда «на слуху», на сегодняшний день существует множество реализаций на различных популярных языках.
Перечислим далее некоторые из них, указывая язык программирования, лежащий в их основе, и некоторые важные характеристики, которые обсуждались в статье ранее:

• Code on netlib.org by Jane Cullum and Ralph Willoughby (Fortran 77; без реортогонализации векторов; без реализации распараллеливания)
• Lanczos algorithm on freesourcecode.net (Matlab; присутствует реортогонализация векторов; без реализации распараллеливания)
• Lanczos for eigenvalues in Numerical Linear Algebra, Spring 2007 (Matlab; частичная или полная реортогонализация - доступны оба варианта; без реализации распараллеливания)
• Planso (название происходит от parallel Lanso.f) code by Kesheng Wu and Horst D. Simon (Fortran 77; частичная реортогонализация векторов; использование MPI для коммуникаций)
• Lanczos.cpp in the Vienna Computing Library (C++, частичная реортогонализация, без реализации распараллеливания)
• Parallel Eigen Values code using C++ AMP on microsoft msdn blog (C++; ?(нужно проанализировать код на наличие ортогонализации)? ; распараллеливание есть и реализовано через C++ AMP)

3 Литература

[1] Парлетт Б. - Симметричная проблема собственных значений. Численные методы //М.: Мир. - 1983. - С. 276-294

[2] James W. Demmel - Вычислительная линейная алгебра. Теория и приложения //Мир. - 2001. С. 381-391

[3] Голуб Дж., Ван Лоун Ч. - Матричные вычисления//М.: Мир. - 1999. - С. 426-436

[4] https://ru.wikipedia.org/wiki/Алгоритм_вычисления_собственных_значений#.D0.9C.D0.B0.D1.82.D1.80.D0.B8.D1.86.D1.8B_.D0.A5.D0.B5.D1.81.D1.81.D0.B5.D0.BD.D0.B1.D0.B5.D1.80.D0.B3.D0.B0_.D0.B8_.D1.82.D1.80.D1.91.D1.85.D0.B4.D0.B8.D0.B3.D0.BE.D0.BD.D0.B0.D0.BB.D1.8C.D0.BD.D1.8B.D0.B5_.D0.BC.D0.B0.D1.82.D1.80.D0.B8.D1.86.D1.8B

[5] http://num-anal.srcc.msu.ru/lib_na/int_ae/int_ae6.htm