Участник:AleksLevin/Алгоритм Ланцоша вычисления собственных значений симметричной матрицы для точной арифметики (без переортогонализации)
Основные авторы описания: Левин А.Д. (студент, кафедра вычислительных методов, 604 группа)
Свойства и структура алгоритмов
Общее описание алгоритма
Данный алгоритм появился в 1950 г. и носит имя венгерского физика и математика Корнелия Ланцоша (венг. Lánczos Kornél). Алгоритм Ланцоша относится к итерационным методам вычисления собственных значений для матриц столь больших порядков [math]\displaystyle{ n }[/math], что к ним нельзя применить прямые методы из-за ограничений по времени и памяти.
Сам Ланцош указывал, что его метод предназначен для отыскания нескольких собственных векторов симметричных матриц, хотя к методу сразу было обращено внимание, как к способу приведения всей матрицы к трёхдиагональному виду. Двадцатью годами позже канадский математик Крис Пэж показал, что, несмотря на чувствительность к округлениям, алгоритм Ланцоша - эффективное средство вычисления некоторых [math]\displaystyle{ k }[/math] собственных чисел и соответствующих им собственных векторов [1, c.276].
Алгоритм Ланцоша для вычисления собственных значений симметричной матрицы [math]\displaystyle{ A }[/math] соединяет в себе метод Ланцоша для построения последовательности подпространств Крылова и ортонормированных векторов Ланцоша и процедуру Рэлея-Ритца получения оптимальных приближений собственных значений и соответствующих собственных векторов исходной матрицы [math]\displaystyle{ A }[/math] [2, с.381].
В данной статье рассматривается вариант алгоритма Ланцоша, в котором опущено влияние ошибок округления на вычислительный процесс, хотя на практике этому посвящается отдельное внимание и существуют различные методы решения данной проблемы, такие как частичная или полная переортогонализация. Этим модифицированным методам посвящены отдельные статьи на данном ресурсе.
Математическое описание алгоритма
Для лучшего понимания описания, данного в этом пункте статьи, рекомендуется ознакомиться с параграфом 6.6 Методы Крыловского подпространства [2, с.313]. Здесь же дано краткое описание всех переменных, математических операций и необходимый теоретический минимум.
Алгоритм Ланцоша для вычисления [math]\displaystyle{ k }[/math] собственных значений и собственных векторов вещественной симметричной матрицы [math]\displaystyle{ A=A^T }[/math] в точной арифметике [2, с.381]:
[math]\displaystyle{ \begin{align} q_1 = & 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} = 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]\displaystyle{ b }[/math] - заданный вещественный вектор. Также полагается известным алгоритм вычисления произведения матрицы [math]\displaystyle{ A }[/math] на вектор [math]\displaystyle{ x }[/math].
Введём матрицу Крылова, определяемую следующим соотношением: [math]\displaystyle{ K_j = [b,Ab,A^2b,...,A^{j-1}b] }[/math].
Далее, на практике, матрица [math]\displaystyle{ K }[/math] заменяется матрицей [math]\displaystyle{ Q }[/math], такой, что при любом числе [math]\displaystyle{ k }[/math] линейные оболочки первых [math]\displaystyle{ k }[/math] столбцов в [math]\displaystyle{ K }[/math] и [math]\displaystyle{ Q }[/math] являются одним и тем же подпространством [2, c.315]. Тогда матрица [math]\displaystyle{ Q }[/math], в отличие от матрицы [math]\displaystyle{ K }[/math], хорошо обусловлена и легко обратима. В результате получаем матрицу [math]\displaystyle{ Q_j = [q_1, q_2, \dots, q_j] }[/math] размерности [math]\displaystyle{ n\times j }[/math], столбцы которой ортогональны, являются базисом подпространства Крылова и называются векторами Ланцоша.
В алгоритме Ланцоша вычислению подлежит столько первых столбцов в матрице [math]\displaystyle{ Q_j }[/math], сколько необходимо для получения требуемого приближения к решению [math]\displaystyle{ A\,x\,=b\,\,(A\,x=\lambda \, x) }[/math].
Получение нулевого [math]\displaystyle{ \beta_j }[/math] в итерационном процессе - событие, желательное в том смысле, что оно свидетельствует о том, что найдено некоторое подпространство, являющееся в точности инвариантным. На практике же нулевое и близкие к нулю значения получаются редко [3, стр. 429].
Затем на каждом шаге цикла формируется симметричная трёхдиагональная матрица [math]\displaystyle{ T_j = Q^T_j A Q }[/math], к которой применяется процесс Рэлея-Ритца для поиска её собственных значений (на практике для поиска этих собственных значений можно использовать любой из специальных методов, изложенных в [3, [math]\displaystyle{ \S }[/math] 8.4]). Эти собственные значения, они же числа Ритца, и полагаются приближёнными собственными значениями исходной матрицы [math]\displaystyle{ A }[/math].
Вычислительное ядро алгоритма
Вычислительное ядро рассматриваемого алгоритма (наиболее ресурсно-затратная часть алгоритма) - формирование на каждом шаге цикла промежуточного вектора [math]\displaystyle{ z }[/math], вычисляемого по формуле: [math]\displaystyle{ z = A\,q_j }[/math]
Макроструктура алгоритма
Как уже было упомянуто в данной статье в описании алгоритма, рассматриваемый алгоритм состоит из двух последовательно выполненных алгоритмов: метод Ланцоша для построения последовательности подпространств Крылова и ортонормированных векторов Ланцоша и процедура Рэлея-Ритца получения оптимальных приближений собственных значений и соответствующих собственных векторов исходной матрицы.
В первой части алгоритма, где итерационно строятся Крыловское подпространство и трёхдиагональная матрица [math]\displaystyle{ T_j }[/math], можно выделить следующие макрооперации:
- умножение матрицы на вектор (в свою очередь представляет собой умножение вектора на число и сложение векторов);
- вычисление нормы (скалярное произведение векторов + вычисление квадратного корня);
- скалярное произведение векторов;
- сложение векторов и их умножение на вещественные числа
Умножение матрицы на вектор - весьма дешевая и нетривиальная операция, если предполагать, что исходная матрица [math]\displaystyle{ A }[/math] разреженная. Именно она подлежит оптимизации с использованием техник распараллеливания, поскольку она является вычислительным ядром алгоритма, а весь упомянутый процесс построения матрицы [math]\displaystyle{ T_j }[/math] выполняется последовательно.
Вторая часть алгоритма - поиск собственных значений матрицы [math]\displaystyle{ T_j }[/math], в нашем случае с помощью процедуры Рэлея-Ритца, может рассматриваться как отдельная макрооперация.
Схема реализации последовательного алгоритма
С целью полностью не дублировать пункт [math]\displaystyle{ \mathbf{1.2} }[/math] данной статьи, где представлен очень наглядный и доступный для понимания псевдокод алгоритма, постараемся лишь коротко просуммировать и обобщить весь перечень последовательных действий в рассматриваемом алгоритме.
В качестве входных данных имеем: вещественная симметричная матрица [math]\displaystyle{ A }[/math] размера [math]\displaystyle{ n\times n }[/math], известный вещественный вектор [math]\displaystyle{ b }[/math] длины [math]\displaystyle{ n }[/math] и число [math]\displaystyle{ k\lt n }[/math].
Далее:
- Инициализируем вещественную константу [math]\displaystyle{ \beta_0 }[/math] и векторы [math]\displaystyle{ q_0 }[/math] и [math]\displaystyle{ q_1 }[/math]
- Итерационно для всех [math]\displaystyle{ j=\overline {1,k} }[/math]:
- Вычисляем вспомогательный вектор [math]\displaystyle{ z }[/math], используя умножение матрицы на вектор
- Вычисляем, используя скалярное произведение векторов, значение [math]\displaystyle{ \alpha_j }[/math] - элемент на главной диагонали трёхдиагональной матрицы [math]\displaystyle{ T_j }[/math]
- Вычисляем по формуле новое значение и перезаписываем вектор [math]\displaystyle{ z }[/math]
- Вычисляем значение [math]\displaystyle{ \beta_j }[/math] - элементы непосредственно над и под главной диагональю той же матрицы [math]\displaystyle{ T_j }[/math]
- Вычисляем с вектора [math]\displaystyle{ z }[/math] очередной столбец матрицы [math]\displaystyle{ Q }[/math], он же вектор Ланцоша : [math]\displaystyle{ \,\,q_{j+1} }[/math]
- Теперь, когда матрица [math]\displaystyle{ T_j }[/math] сформирована на очередном шаге, находим её собственные значения и вектора, например с помощью процедуры Рэлея-Ритца (о используемых на практике методах нахождения подробнее было сказано ранее в математическом описании алгоритма)
- Переходим на следующую итерацию: [math]\displaystyle{ j }[/math]++
Последовательная сложность алгоритма
Оценим количество операций и сложность последовательного алгоритма, изложенного в предыдущем пункте.
- Для вычисления вектора [math]\displaystyle{ z }[/math] умножается матрица размера [math]\displaystyle{ n\times n }[/math] на вектор длины [math]\displaystyle{ n }[/math], и потребуется [math]\displaystyle{ n^2 }[/math] операций умножения и [math]\displaystyle{ n^2 }[/math] операций сложения.
- Умножение векторов [math]\displaystyle{ q^T_j }[/math] и [math]\displaystyle{ z }[/math] длины [math]\displaystyle{ n }[/math] потребует [math]\displaystyle{ n }[/math] операций умножения и [math]\displaystyle{ n }[/math] операций сложения.
- При вычислении нового значения вектора [math]\displaystyle{ z }[/math] поэлементное сложение двух векторов длины [math]\displaystyle{ n }[/math] требует [math]\displaystyle{ n }[/math] операций сложения.
- Умножение вектора длины [math]\displaystyle{ n }[/math] на вещественное число требует [math]\displaystyle{ n }[/math] операций умножения.
- Нахождение квадратичной нормы вектора [math]\displaystyle{ z }[/math] длины [math]\displaystyle{ n }[/math] требует [math]\displaystyle{ n }[/math] умножений и [math]\displaystyle{ n }[/math] сложений и ещё 1 операцию извлечения квадратного корня.
- Заключительное нахождение собственных значений и собственных векторов трёхдиагональной матрицы [math]\displaystyle{ T_j }[/math] размера [math]\displaystyle{ j\times j }[/math] с помощью QR-алгоритма потребует в худшем случае [math]\displaystyle{ O(j^3) }[/math] операций[1].
Посчитаем теперь суммарное количество простейших операций на одной итерации алгоритма (без учёта записей в память):
- [math]\displaystyle{ n^2+5\,n }[/math] операций умножения/деления
- [math]\displaystyle{ n^2+4\,n }[/math] операций сложения/вычитания
- 1 операция извлечения квадратного корня
И отдельно после всех итераций при [math]\displaystyle{ j=k }[/math] :
- [math]\displaystyle{ O(k^3) }[/math] операций для поиска собственных значений и собственных векторов трёхдиагональной матрицы [math]\displaystyle{ T_k }[/math]
Итого можно утверждать следующее:
Последовательная сложность одной итерации рассматриваемого алгоритма: [math]\displaystyle{ O(n ^ 2) }[/math]
Суммарная сложность алгоритма ([math]\displaystyle{ k }[/math] итераций): [math]\displaystyle{ O(k\,n ^ 2) + O(k ^ 3) }[/math] [2]
На практике же обычно рассматривается число [math]\displaystyle{ k\lt \lt n }[/math], поэтому второе слагаемое можно опустить.
Информационный граф
Ресурс параллелизма алгоритма
Входные и выходные данные алгоритма
Свойства алгоритма
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
Локальность данных и вычислений
Возможные способы и особенности параллельной реализации алгоритма
Масштабируемость алгоритма и его реализации
Динамические характеристики и эффективность реализации алгоритма
Выводы для классов архитектур
Существующие реализации алгоритма
Литература
[1] Парлетт Б. - Симметричная проблема собственных значений. Численные методы //М.: Мир. - 1983. - С. 276-294
[2] James W. Demmel - Вычислительная линейная алгебра. Теория и приложения //Мир. - 2001. С. 381-391
[3] Голуб Дж., Ван Лоун Ч. - Матричные вычисления//М.: Мир. - 1999. - С. 426-436
[5] http://num-anal.srcc.msu.ru/lib_na/int_ae/int_ae6.htm
- ↑ http://num-anal.srcc.msu.ru/lib_na/int_ae/int_ae6.htm (пункт 6.2)
- ↑ Здесь мы полагаем, что мы не вычисляем собственные значения и векторы матрицы [math]\displaystyle{ T_j }[/math] на каждом промежуточном шаге, а лишь единожды затрачиваем [math]\displaystyle{ O(k ^ 3) }[/math] действий в конце для нахождения [math]\displaystyle{ k }[/math] собственных значений и векторов матрицы [math]\displaystyle{ T_k }[/math]