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

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

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

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

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

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

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) \\ 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.5 Схема реализации последовательного алгоритма

1. Находим с помощью кругов Гершгорина интервал [math](a, b)[/math] на котором содержатся все собственные значения матрицы. Во избежание проблем с точным попаданием в собственные значения к [math]a[/math] и [math]b[/math] можно добавить случайные возмущения, расширяющие интервал.
2. Инициализируем очередь интервалов интервалом [math](a, b, 1, n)[/math], где [math]a[/math] и [math]b[/math] обозначают границы интервалов, а [math]1[/math] и [math]n[/math] номера собственных значений, которые на нем находятся.
3. Извлекаем интервал [math](l, r, p, q)[/math] из очереди
4. Если [math]p \ne q[/math], переходим к шагу 999
5. Вычисляем пары [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], то добавляем интервал в список локализованных интервалов и переходим к шагу 999
6. Вычисляем [math]\lambda = (l + r) / 2[/math] и [math]neg_n(\lambda)[/math]
7. Если [math]neg_n(\lambda) \le p - 1[/math], то добавляем интервал [math](\lambda, r, p, q)[/math] в очередь и переходим к шагу 999
8. Если [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), q)[/math].
9. Если очередь не пуста, переходим к шагу 999