Участник: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.4 Масштабируемость алгоритма и его реализации

To be done soon!

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