Участник:Stanis-morozov/Алгоритм Тренча вычисления спектра симметрической теплицевой матрицы

Алгоритм Тренча вычисления спектра симметрической теплицевой матрицы. Автор описания: С.В. Морозов.

Задача нахождения собственных значений и собственных векторов для матрицы [math]A[/math] заключается в поиске таких соответствующих друг другу чисел [math]\lambda[/math] и ненулевых векторов [math]x[/math], которые удовлетворяют уравнению [math]Ax=\lambda x[/math]. Числа [math]\lambda[/math] называются собственными значениями, а вектора [math]x[/math] - соответствующими им собственными векторами.^[1]

Теплицевой матрицей называется матрица, в которой на всех диагоналях, параллельных главной, стоят равные элементы.^[2]

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

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

Задача эффективного вычисления собственных значений матриц является одной из основных задач линейной алгебры. Известно, что собственные значения произвольной матрицы невозможно вычислить за конечное число операций, однако существует множество итерационных методов, позволяющих решать эту задачу, например, QR-алгоритм. Однако, для некоторых специальных классов матриц можно предложить более эффективные методы, чем QR-алгоритм. Так, в 1989 была опубликована работа В. Тренча ^[3], в которой был предложен метод вычисления собственных значений симметрической теплицевой матрицы, учитывающий особенности структуры таких матриц, и допускающий эффективную параллельную реализацию ^[4]^[5]^[6].

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

Пусть имеется матрица [math]T_{n}[/math] порядка [math]n\times n[/math] следующего вида:

[math] T_{n} = \begin{bmatrix} t_{0} & t_{1} & t_{2} & \cdots & t_{n-2} & t_{n-1} \\ t_{1} & t_{0} & t_{1}& \cdots & t_{n-3} & t_{n-2} \\ t_{2} & t_{1} & t_{0} & \cdots & t_{n-4} & t_{n-3} \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ t_{n-2} & \cdots & \cdots & t_{1} & t_{0} & t_{1} \\ t_{n-1} & \cdots & \cdots & t_{2} & t_{1} & t_{0} \\ \end{bmatrix} [/math]

Требуется найти ее собственные значения. Будем также предполагать, что ни одно из собственных значений [math]T_{n}[/math] не является собственным значением для [math]T_{1}, \dots, T_{n-1}[/math]. Однако, практические эксперименты, проведенные во многих работах, показывают, что данное ограничение существенно только в теории и не дает о себе знать в практических вычислениях.

Введем следующие обозначения. Обозначим

[math]p_{m}(\lambda) = det[T_{m} - \lambda I_{m}][/math]

[math]q_{m}(\lambda) = p_{m}(\lambda) / p_{m-1}(\lambda)[/math]

Метод Тренча состоит в нахождении корней рациональной функции [math]q_{n}(\lambda)[/math]. При сделанных предположениях относительно матрицы [math]T_{n}[/math], эта задача равносильна задаче нахождения собственных значений матрицы [math]T_{n}[/math].

Обозначим также [math] v_{m} = \begin{bmatrix} t_{1} \\ t_{2} \\ \vdots \\ t_{m} \\ \end{bmatrix} [/math]

Можно доказать, что если [math]\lambda[/math] не является собственным значением [math]T_{n}[/math], а [math]z_{m}(\lambda)[/math] является решением

[math](T_{m} - \lambda I_{m})z_{m}(\lambda) = v_{m}[/math],

то

[math]q_{m}(\lambda) = t_{0} - \lambda - v_{m-1}^{T} z_{m-1}(\lambda)[/math].

Для решения системы можно использовать следующий модифицированный алгоритм Левинсона-Дурбина, эквивалентный в данном случае построению [math]LDL^{T}[/math] разложения для матрицы [math]T_{m} - \lambda I_{m}[/math]. Вычисления алгоритма проходят по следующим формулам:

[math] q_{1}(\lambda) = t_{0} - \lambda \\ z_{1,1}(\lambda) = t_{1} / q_{1}(\lambda) \\ q_{m}(\lambda) = (1 - z_{m-1, m-1}^2(\lambda)) q_{m-1}(\lambda) \\ z_{m,m}(\lambda) = q_{m}^{-1}(\lambda)\big(t_{m} - \sum\limits_{j=1}^{m-1} t_{m-j}z_{j,m-1}(\lambda)\big) \\ z_{j,m}(\lambda) = z_{j,m-1}(\lambda) - z_{m,m}(\lambda)z_{m-j,m-1}(\lambda), 1 \le j \le m-1 [/math]

С другой стороны, заметив, что на диагонали матрицы [math]D[/math] в полученном [math]LDL^{T}[/math] разложении стоят значения [math]q_{j}(\lambda)[/math], и применяя закон инерции квадратичных форм получаем, что количество элементов [math]q_{j}(\lambda)[/math] меньших нуля, в точности равно количеству собственных значений матрицы [math]T_n[/math], меньших [math]\lambda[/math].

Кроме того, приведем важный для эффективного вычисления результат:

Теорема. Предположим, что [math]\alpha[/math] и [math]\beta[/math] не являются собственными значениями ни для какой из матриц [math]T_{1}, \dots, T_{n}[/math], и что [math](\alpha, \beta)[/math] содержит в точности одно собственное значение [math]T_{n}[/math] кратности единица. Тогда [math](\alpha, \beta)[/math] не содержит собственных значений матрицы [math]T_{n-1}[/math] тогда и только тогда, когда [math]q_{n}(\alpha) \gt 0[/math] и [math]q_{n}(\beta) \lt 0[/math].

Отсюда вытекает следующий алгоритм нахождения [math]j[/math]-го по величине собственного значения матрицы [math]T_{n}[/math].

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

2. Локализуем [math]j[/math]-ое собственное значение. Для этого используем метод деления отрезка пополам. Обозначим [math]neg_{n}(\lambda)[/math] количество собственных значений [math]T_{n}[/math] меньших [math]\lambda[/math]. Имея процедуру вычисления [math]neg_{n}(\lambda)[/math] и [math]q_{n}(\lambda)[/math] легко построить алгоритм, сходящийся к [math]\alpha[/math] и [math]\beta[/math], таким что

[math]neg_{n}(\alpha) = j - 1[/math] и [math]neg_{n}(\beta) = j[/math]

[math]q_{n}(\alpha) \gt 0[/math] и [math]q_{n}(\beta) \lt 0[/math]

3. Тогда на этом интервале функция [math]q_{n}(\lambda)[/math] будет непрерывна (и более того, бесконечно дифференцируема) и можно применять другие, более эффективные методы нахождения корня уравнения.

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

В алгоритме два основных вычислительных ядра. Первое состоит в локализации всех собственных значений, т.е. нахождении для каждого собственного значения матрицы [math]\lambda_j[/math] такого отрезка [math](a, b)[/math], что он содержит в себе [math]\lambda_j[/math] и не содержит других собственных значений матрицы, а также собственных значений матрицы меньшей размерности.

Второе состоит в приближении локализованного на интервале собственного значения с заданной точностью.

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

Введем базовые макрооперации, необходимые для описания алгоритма. Одна из базовых операций в алгоритме - это нахождение интервала, содержащего все собственные значения матрицы с помощью кругов Гершгорина. На псевдокоде эта задача может быть реализована следующим образом (за t_j обозначен элемент теплицевой матрицы t_kl для которого k - l = j):

   mx = 0
   for i = 1, n:
       s₁ = s₂ = 0
       for j = 2, i:
           s1 += |t_j|
       for j = 2, n - i:
           s2 += |t_j|
       mx = max(mx, s1 + s2)
   a = t₀ - mx
   b = t₀ + mx

Еще одна базовая макрооперация - это вычисление величин [math]neg_n(\lambda)[/math] и [math]q_n(\lambda)[/math] происходит одновременно по следующим формулам:

   q₁(λ) = t₀ - λ
   z_1,1(λ) = t₁ / q₁(λ)
   q_m(λ) = (1 - z²_{m-1, m-1}(λ)) q_m-1(λ)
   tmp = 0
   for j = 1, m - 1:
       tmp += t_m-jz_j,m-1(λ)
   z_m,m(λ) = q^-1_m(λ)(t_m - tmp)
   for j = 1, m:
       z_j,m(λ) = z_j,m-1(λ) - z_m,m(λ)z_m-j,m-1(λ)
   neg_n(λ) = 0
   for j = 1, n:
       if q_j(λ) < 0:
           neg_n(λ) += 1

Теперь, опишем операцию локализации собственного значения. Пусть имеется интервал [math](a, b)[/math]. Будем дробить его пополам, и с помощью величины [math]neg_n((a+b) / 2)[/math] находить количество собственных значений на каждом из получившихся интервалов. Как утверждает теорема, сформулированная выше, рано или поздно такое дробление приведет нас к интервалам, на которых выполняются следующие неравенства: [math]neg_n(a) = j - 1[/math], [math]neg_n(b) = j[/math], [math]q_n(a) \gt 0[/math], [math]q_n(b) \lt 0[/math]. На этом локализация заканчивается. Более детальное описание алгоритма локализации приведено в следующем разделе.

Кроме этого, в алгоритме присутствует операция уточнения собственного значения. Она состоит в последовательном дроблении локализованного интервала, с последующей проверкой с помощью [math]neg_n(\lambda)[/math] в какой из частей лежит собственное значения. Интервал сжимается до достижения требуемой точности.

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

С помощью кругов Гершгорина находим интервал, на котором содержатся все собственные значения
Локализуем все собственные значения
Уточняем все локализованные собственные значения

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

Находим с помощью кругов Гершгорина интервал [math](a, b)[/math] на котором содержатся все собственные значения матрицы. Во избежание проблем с точным попаданием в собственные значения к [math]a[/math] и [math]b[/math] можно добавить случайные возмущения, расширяющие интервал.
Инициализируем очередь интервалов интервалом [math](a, b, 1, n)[/math], где [math]a[/math] и [math]b[/math] обозначают границы интервалов, а [math]1[/math] и [math]n[/math] номера собственных значений, которые на нем находятся.
Извлекаем интервал [math](l, r, p, q)[/math] из очереди
Если [math]p \ne q[/math], переходим к шагу 6
Вычисляем пары [math](neg_n(l), q_n(l))[/math] и [math](neg_n(r), q_n(r))[/math]. Если [math]neg_n(l) = p - 1[/math], [math]neg_n(r) = p[/math], [math]q_n(l) \gt 0[/math] и [math]q_n(r) \lt 0[/math], то добавляем интервал [math](l, r, p)[/math] в список локализованных интервалов и переходим к шагу 3
Вычисляем [math]\lambda = (l + r) / 2[/math] и [math]neg_n(\lambda)[/math]
Если [math]neg_n(\lambda) \le p - 1[/math], то добавляем интервал [math](\lambda, r, p, q)[/math] в очередь и переходим к шагу 3
Если [math]neg_n(\lambda) \ge q[/math], то добавляем интервал [math](l, \lambda, p, q)[/math] в очередь. Иначе добавляем в очередь [math](l, \lambda, p, neg_n(\lambda))[/math] и [math](\lambda, r, neg_n(\lambda) + 1, q)[/math].
Если очередь не пуста, переходим к шагу 3
Извлекаем интервал [math](l, r, j)[/math] из списка локализованных интервалов
Если [math]r - l \le \epsilon[/math] переходим к шагу 14
Вычисляем [math]\lambda = (l + r) / 2[/math] и [math]neg_n(\lambda)[/math]
Если [math]neg_n(\lambda) \le p - 1[/math], то [math]l = \lambda[/math], иначе [math]r = \lambda[/math]. Переходим к шагу 11
Запоминаем [math]j[/math]-ое собственное значение как [math](l + r) / 2[/math]. Если список локализованных интервалов не пуст, переходим к шагу 10.

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

Число операций, необходимых для вычисления спектра матрицы зависит от требуемой точности вычислений. Легко видеть, что число операций, необходимых для вычисления [math]neg_n(\lambda)[/math] и [math]q_n(\lambda)[/math] равно [math]O(N^2)[/math]. Вычисление начально интервала, содержащего все собственные значения также требует [math]O(N^2)[/math] операций. Также нетрудно понять, что число операций, необходимых для приближения локализованного собственного значения с заданной точностью равно [math]O(N^2\log \frac{M}{\epsilon})[/math], где [math]\epsilon[/math] - требуемая точность, а [math]M[/math] - величина исходного интервала ([math]O(N^2)[/math] на каждую итерацию и [math]\log \frac{M}{\epsilon}[/math] итераций). Теоретическая оценка числа итераций, необходимых для локализации одного собственного значения неизвестна, однако при достаточно высокой требуемой точности (такой, чтобы после локализации было необходимо дальнейшее приближение собственного значения), это число итераций, очевидно, меньше [math]\log \frac{M}{\epsilon}[/math]. Далее, так как необходимо вычислить каждое собственное значение, итоговое число операций равно [math]O(N^3\log \frac{M}{\epsilon})[/math], где [math]N[/math] - размерность матрицы, [math]\epsilon[/math] - требуемая точность приближения, а [math]M = \|T_n\|_{\infty}[/math]. Стоит заметить также, что алгоритм позволяет вычислять не весь спектр сразу, а лишь необходимые [math]k[/math] собственных значений. Тогда для этого потребуется [math]O(k N^2\log \frac{M}{\epsilon})[/math] операций.

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

Рисунок 1. Информационный граф вычисления собственных значений симметрической теплицевой матрицы

Опишем значения вершин информационного графа.

Вершина [math]T_n[/math] обозначает входные данные алгоритма - матрицу [math]T_n[/math]
Вершина [math]G[/math] обозначает вычисление интервала, содержащего все собственные значения с помощью кругов Гершгорина
Вершина [math]C[/math] обозначает условный переход при завершении локализации некоторого собственного значения, к его уточнению. Передача данных вниз в графе происходит только при срабатывании условия. Передача данных вправо происходит всегда
Вершина [math]P[/math] обозначает процесс дробления следующего интервала из очереди при локализации собственных значений
Вершина [math]F[/math] обозначает процесс уточнения локализованного собственного значения
Вершина [math]\lambda_j[/math] обозначает выходные данные алгоритма - вычисленное собственное значение

Стоит заметить, что число вершин типов [math]C[/math] и [math]P[/math] не равно [math]n[/math], а равно количеству необходимых итераций дробления для достижения локализации собственных значений.

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

Глядя на последовательную версию алгоритма, становится ясно, что локализовывать собственные значения можно не в одной последовательной очереди, а разбив номера от [math]1[/math] до [math]n[/math] на [math]d[/math] равных частей (где [math]d[/math] - число доступных узлов) и инициализировать очередь лишь подотрезком на каждом узле. Также, приближения локализованных собственных значений можно производить независимо на каждом узле. Это означает полную распараллеливаемость алгоритма. Тогда число операций, необходимое для вычисления спектра равно [math]O\Big(\frac{N^3\log \frac{M}{\epsilon}}{d}\Big)[/math] при [math]d \lt N[/math]. При наличии бесконечного числа узлов, асимптотическая сложность будет равна [math]O(N^2\log \frac{M}{\epsilon})[/math].

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

На вход алгоритму подается размер матрицы [math]n[/math] и симметрическая теплицева матрица [math]T_n[/math] размера [math]n \times n[/math]. Также предполагается, что ни одно из собственных значений [math]T_{n}[/math] не является собственным значением для [math]T_{1}, \dots, T_{n-1}[/math], однако утверждается, что это ограничение не затрудняет практического применения. Также на вход подается число [math]\epsilon \gt 0[/math] - требуемая абсолютная погрешность вычисления собственных значений.

Результатом работы алгоритма является массив из [math]n[/math] вещественных чисел, [math]j[/math]-ое из которых с точностью до [math]\epsilon[/math] равно [math]j[/math]-ому собственному значению матрицы [math]T_n[/math].

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

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

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

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

Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета. Исходный код был написан на языке C++ с использованием библиотеки Intel MPI (Intel(R) MPI Library for Linux* OS, Version 5.0 Update 1). Компиляция производилась с опциями -O3 -march=native -std=c++0x.

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

число процессоров [20 : 240] с шагом 20;

размер матрицы [200 : 2400] c шагом 200.

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

Рисунок 2. Параллельная реализация алгоритма вычисления собственных значений симметрической теплицевой матрицы

Исходный код тестируемой реализации на C++

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

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

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

3 Литература

↑ Е.Е. Тыртышников Основы алгебры. М.: Физматлит, 2017.
↑ Е.Е. Тыртышников Методы численного анализа. М.: Издательский центр "Академия", 2007, 320 с.
↑ Trench, W. F. (1989). Numerical solution of the eigenvalue problem for Hermitian Toeplitz matrices. SIAM J. Matrix Anal. Appl., 10(2): 135-146
↑ F. Noor and S.D. Morgera, Recursive and iterative algorithms for computing eigenvalues of Hermitian Toeplitz matrices, IEEE Trans. Signal Process. 41 (1993) 1272–1280.
↑ J. M. Badía and A. M. Vidal, Parallel Algorithms to Compute the Eigenvalues and Eigenvectors of Symmetric Toeplitz Matrices, Parallel Algorithms and Applications, Vol. 13, pp. 75-93. (2000).
↑ Noor F., Misbahuddin S. (2010) Using MPI on PC Cluster to Compute Eigenvalues of Hermitian Toeplitz Matrices. In: Hsu CH., Yang L.T., Park J.H., Yeo SS. (eds) Algorithms and Architectures for Parallel Processing. ICA3PP 2010. Lecture Notes in Computer Science, vol 6081. Springer, Berlin, Heidelberg

[1] Е.Е. Тыртышников Основы алгебры. М.: Физматлит, 2017.

[2] Е.Е. Тыртышников Методы численного анализа. М.: Издательский центр "Академия", 2007, 320 с.

[3] Trench, W. F. (1989). Numerical solution of the eigenvalue problem for Hermitian Toeplitz matrices. SIAM J. Matrix Anal. Appl., 10(2): 135-146

[4] F. Noor and S.D. Morgera, Recursive and iterative algorithms for computing eigenvalues of Hermitian Toeplitz matrices, IEEE Trans. Signal Process. 41 (1993) 1272–1280.

[5] J. M. Badía and A. M. Vidal, Parallel Algorithms to Compute the Eigenvalues and Eigenvectors of Symmetric Toeplitz Matrices, Parallel Algorithms and Applications, Vol. 13, pp. 75-93. (2000).

[6] Noor F., Misbahuddin S. (2010) Using MPI on PC Cluster to Compute Eigenvalues of Hermitian Toeplitz Matrices. In: Hsu CH., Yang L.T., Park J.H., Yeo SS. (eds) Algorithms and Architectures for Parallel Processing. ICA3PP 2010. Lecture Notes in Computer Science, vol 6081. Springer, Berlin, Heidelberg

[1]

[2]

[3]

[4]

[5]

[6]